Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Molecular Signatures in Arabidopsis thaliana in Response to Insect Attack and Bacterial Infection

  • Pankaj Barah,

    Affiliation Department of Biology, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

  • Per Winge,

    Affiliation Department of Biology, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

  • Anna Kusnierczyk,

    Affiliation Department of Biology, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

  • Diem Hong Tran,

    Affiliation Department of Biology, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

  • Atle M. Bones

    atle.bones@bio.ntnu.no

    Affiliation Department of Biology, Norwegian University of Science and Technology (NTNU), Trondheim, Norway

Abstract

Background

Under the threat of global climatic change and food shortages, it is essential to take the initiative to obtain a comprehensive understanding of common and specific defence mechanisms existing in plant systems for protection against different types of biotic invaders. We have implemented an integrated approach to analyse the overall transcriptomic reprogramming and systems-level defence responses in the model plant species Arabidopsis thaliana (A. thaliana henceforth) during insect Brevicoryne brassicae (B. brassicae henceforth) and bacterial Pseudomonas syringae pv. tomato strain DC3000 (P. syringae henceforth) attacks. The main aim of this study was to identify the attacker-specific and general defence response signatures in A. thaliana when attacked by phloem-feeding aphids or pathogenic bacteria.

Results

The obtained annotated networks of differentially expressed transcripts indicated that members of transcription factor families, such as WRKY, MYB, ERF, BHLH and bZIP, could be crucial for stress-specific defence regulation in Arabidopsis during aphid and P. syringae attack. The defence response pathways, signalling pathways and metabolic processes associated with aphid attack and P. syringae infection partially overlapped. Components of several important biosynthesis and signalling pathways, such as salicylic acid (SA), jasmonic acid (JA), ethylene (ET) and glucosinolates, were differentially affected during the two the treatments. Several stress-regulated transcription factors were known to be associated with stress-inducible microRNAs. The differentially regulated gene sets included many signature transcription factors, and our co-expression analysis showed that they were also strongly co-expressed during 69 other biotic stress experiments.

Conclusions

Defence responses and functional networks that were unique and specific to aphid or P. syringae stresses were identified. Furthermore, our analysis revealed a probable link between biotic stress and microRNAs in Arabidopsis and, thus gives indicates a new direction for conducting large-scale targeted experiments to explore the detailed regulatory links between them. The presented results provide a comparative understanding of ArabidopsisB. brassicae and ArabidopsisP. syringae interactions at the transcriptomic level.

Introduction

Plants are sessile organisms that are unable to escape biotic and abiotic stresses. As a result, they have evolved flexibility in their responses to changing environmental conditions, such as light, drought, temperature, the available nutritional supply and biotic invasion. Different types of biotic invasions, such as insect, bacterial, fungal and viral invasions, represent a severe threat to agricultural production worldwide [1]. Some responses of host plants to different stress conditions are very general and provide protection from a variety of invading organisms, whereas others are more specific and target particular types of attackers. Highly complex and often connected signalling pathways, regulating numerous metabolic networks, coordinate plant responses to different stress conditions. Over the last decade or so, clear advances have been made in understanding how defence responses are orchestrated in higher plants. The development of microarray technology has allowed monitoring of expressional changes in thousands of genes simultaneously, and this technology has now become a major tool for examining plant stress biology. Most of these studies have adopted A. thaliana as a model plant organism because of the vast amount of genomic information made available for this species with the completion of the A. thaliana genome sequence and advanced annotation of A. thaliana genes [2]. Analysing the regulation of gene expression under various stress conditions has revealed that the early defence responses of a plant to different stress factors often overlap and engage the same sets of genes [3]. It has also become evident that different types of plant invaders may induce substantially different changes in the host plant transcriptome. Furthermore, studies on plants subjected to various treatments indicate that the induced defences can be both general – being commonly manifested regardless of the type of applied treatment; and specific – providing protection from a certain type of stress [4], [5]. In many cases, however, the multidimensional level of network crosstalk makes it challenging to recognise which of the observed responses are general and which are more stress specific [6], [7].

Aphids are one of the world’s major insect pests, causing serious economic damage to a range of temperate and tropical crops [8]. Aphids use their mouthparts, formed into a stylet-like structure, to pierce plant tissue in the search for sieve elements (SEs) containing their primary food source: phloem sap [9], [10]. Feeding by an aphid causes minimal wounding, as its stylet proceeds mostly intercellulary and is inserted only into selected cells on its way to the phloem tissue [11]. However, the disruption of cell walls and membranes of the pierced cells is likely to be the first factor triggering a plant response. In addition, the salivary secretions lubricating the stylet throughout its pathway through plants tissues and injected into SEs during feeding contain molecular signatures that activate plant defences. Therefore, despite their stealthy feeding, aphids are strong inducers of plant defences against them. Recently Kuśnierczyk et al. reported the timing and dynamics of early Arabidopsis defence responses [12] to an aphid attack.

P. syringae is a bacterial leaf pathogen that causes extensive chlorosis and necrotic spots [13]. Many strains of P. syringae are pathogenic in the model plant A. thaliana, and P. syringae is therefore widely used to study plant – pathogen interactions under laboratory conditions. P. syringae enters host tissues through wounds or natural openings such as stomata, and in susceptible plants, it multiplies to high concentrations in intercellular spaces [14]. The ability of P. syringae to multiply endophytically is dependent on its type III secretion pathway enabling the secretion of proteins into the apoplast. These proteins interact with the cell wall and plasma membrane and are directly translocated into the cytoplasm of host cells [15]. Several strains of P. syringae produce coronatine, a molecule that mimics endogenous plant jasmonyl-L-isoleucine and an activator of the jasmonic acid signalling pathway [16]. By doing so, the bacteria manipulate host responses, suppressing salicylic acid defences through the activation of jasmonic acid signalling [17], [18].

A great number of experiments conducted to assess plant responses to different stresses have made substantial contributions to our understanding of the induced defences of plants. However, the comparison of independent experiments and extraction of meaningful information from such comparisons is complicated and difficult in most cases, mainly due to the lack of common standards regarding how to grow plants, conduct expression profile experiments, and finally, how to evaluate the resulting gene expression data [19]. In recent years, integrated approaches, such as systems biology methods, have been evolving, providing promising tools for studying plant stress responses [20], [21]. Scientists intend to go beyond simple functional enrichment analyses to understand the molecular basis of genome-scale microarray experiments. Methods inspired by systems biology utilise lists of differentially expressed genes ranked by biological criteria to search for the distribution of blocks of functionally related genes without imposing any artificial threshold. Such ranked lists of genes can be arranged into functional classes, pathways and biological processes. Co-expression or co-regulation of particular genes can indicate their involvement in similar biological processes, meaning that individual modules of genes can be attributed to specific biological processes. Using this basic concept, modular network topology-based analysis has been proven to be useful in identifying functional modules of genes [22]. In a recent co-expression study, Weston and co-workers showed how a co-expression network-based analysis can be used for understanding population-level adaptive physiological responses of plants to abiotic stress [23].

MicroRNAs (microRNAs) are small, non-coding RNAs that play critical roles in post-transcriptional gene regulation and stress-inducible transcriptional regulation in Arabidopsis [24]. In plants, mature microRNAs pair with complementary sites on mRNAs, subsequently leading to the cleavage and degradation of the mRNAs. Many microRNAs target mRNAs that encode transcription factors and, thus, influence the expression of many genes whose regulation is controlled by these transcription factors [25]. The identification, detection, regulation and functional analysis of microRNAs associated with biotic stress remains a great challenge. In contrast, information about plant stress-responsive genes and their transcription factor binding sites is available to some extent in several databases [26], [27], [28], [29], [30]. Integration of such publicly available knowledge bases with experimental approaches would provide useful insights in understanding the plant defence responses to different biotic stresses.

In this manuscript, we present such an integrated approach to explore the common (general) and attacker-specific defence responses of A. thaliana subjected to two different types of biotic invaders: phloem-feeding aphids (B. brassicae) and pathogenic bacteria (P. syringae). To allow comparison between the obtained gene expression profiles and the observed regulation of gene pathways involved in defence against the aphid and the bacterium, the same growth and experimental conditions were used in the two simultaneous experimental setups. Transcriptional changes resulting either from infestation with B. brassicae or infection with P. syringae were assessed with the use of full-genome Arabidopsis microarrays (the data have been deposited in GEO with accession numbers GSE39245 and GSE39246).

Two sets of differentially expressed genes, corresponding to the plant responses to either aphid or bacterial treatment, were created as the outcome of the microarray data analysis. In an attempt to integrate the resulting data with publicly available knowledge extracted from several different databases as well as from published results of other experiments, these two differentially regulated gene sets were subsequently analysed through a set of computational approaches. The following analyses were incorporated into the presented work: an analysis of enriched functional categories or processes; exploration of potential connections between microRNAs and biotic stress-inducible transcriptional regulation during insect and bacterial attack; cross-validation of the aphid- and Pseudomonas-regulated genes using a co-expression network constructed from a compendium of 69 other biotic stress microarray datasets complied in the CORNET tool [31] (https://cornet.psb.ugent.be/).

Results and Discussion

Overall Changes in the Arabidopsis Transcriptome in Response to Insect and Bacterial Attack

To explore the complexity of the transcriptional changes induced by the different examined A. thaliana attackers, we compared the overlap between the obtained gene sets. From the results, it is evident that the transcriptional responses of A. thaliana to these very different attackers are massive. Aphid infestation and P. syringae infection resulted in significant differential regulation of 4,979 (2,803 up-regulated, 2,176 down-regulated) and 3,199 (1,634 up, 1,565 down) genes, respectively (Table 1 and Tables S4, S5). Although aphids and bacteria exhibit very different modes of action and trigger a highly dissimilar signal signature, a large number of Arabidopsis genes were expressed in response to both attackers. There were 1,597 common genes affected after both aphid infestation and P. syringae treatment. A total of 3,382 genes (1,963 up, 1,419 down) showed aphid-specific expression, while 1602 genes (842 up, 760 down) showed P. syringae-specific expression (Table S6). In the common set of genes, there were a total of 186 genes that showed opposite expression patterns in the two experiments. Of these genes, 117 were up-regulated under aphid and down-regulated under P. syringae attack, while 69 genes were down-regulated under aphid and up-regulated under P. syringae attack. Out of the 117 genes that were up-regulated in the aphid and down-regulated in the P. syringae experiment, 17 have been reported to be transcription factors. Six of these transcription factors are members of the ERF/AP2 transcription factor family. Among them ERF104, which is regulated by MPK6, is a key controller of innate immunity and dehydration stress [32].

thumbnail
Table 1. Overall summary of the differentially regulated genes in A. thaliana during Brevicoryne brassicae (aphid) attack or P. syringae (bacteria) infection.

https://doi.org/10.1371/journal.pone.0058987.t001

In total, 303 transcription factors were found to be affected by the aphid treatment, while 191 transcription factors showed altered expression under P. syringae infection. The common category (differentially expressed during both aphid infestation and P. syringae treatment) included 87 known Arabidopsis transcription factors. The analysis also identified 216 transcription factors that were differentially regulated only during the aphid treatment and 104 transcription factors that were differentially regulated only during P. syringae infection. The annotated network of these transcripts showed that some of the differentially expressed transcription factors could be crucial for stress-specific defence responses in A. thaliana plants.

Analysis of overrepresented gene ontologies (GO) in A. thaliana indicates rigorous reprogramming of several biological processes. As seen from the Table 1, a large number of genes were differentially regulated in A. thaliana during both the aphid and P. syringae experiments, which indicated that intense transcriptional reprogramming took place. A network-based analysis of the corresponding GO terms under the Biological Process classification using ClueGO (correction method = Bonferroni, kappa score ≥0.3) in the common aphid-specific and P. syringae-specific transcript dataset was performed.

When this analysis was applied to the list of 1,597 common genes whose expression was affected during both of the experiments, 17 significantly overrepresented categories were identified (some of these categories are shown in Figure 1.) Most of the cellular and metabolic processes were clustered in distinctly separate modules, and there were few highly interconnecting overrepresented processes. More than half of the genes from the common list were involved in central metabolic and cellular processes, such as electron transport and energy pathways located in the plastid. Some of the most significant categories were indole-containing compound metabolic processes, host localised cell death, cellular responses to starvation, downregulation of photosynthesis, responses to jasmonic acid, sulphur compound biosynthetic processes, and negative regulation of cellular processes. Analysis of the modules showed that the majority of the jasmonic acid responsive genes were up-regulated by both treatments, but the number of genes and their degree of induction were markedly higher in the P. syringae-treated plants, which may be due to the effects of coronatine (file S11). It has been previously reported that P. syringae uses the virulence factor coronatine (COR) as a mimic of jasmonyl-l-isoleucine (JA-Ile) [16], [33]. The coronatin-regulated A. thaliana genes reported in Thilmony R. et al., 2006 [34] show strong overlap with our P. syringae data. More than 450 genes reported as coronatin-regulated by Thilmony R. et al. show highly similar expression patterns in the two datasets (data not shown).

thumbnail
Figure 1. Over-represented GO-categories in the common gene list.

Network representations of enriched GO categories among the genes that were differentially regulated during both experiments. Figure generated from the functionally grouped networks of enriched GO categories among genes whose expression is induced by both the aphid and pathogenic bacterium treatments. GO terms are represented as nodes based on their kappa score (≥0.3); only networks with at least three nodes are represented. The node size indicates the significance of the term’s enrichment. The edges are related to the relationships between the selected terms, which are defined based on the genes that are shared in a similar way. The label of the most significant term is used as the leading group term. Visualisation was conducted using Cytoscape 2.7.0.

https://doi.org/10.1371/journal.pone.0058987.g001

Tryptophan-derived indolic compounds, such as indolic glucosinolates (iGS) and indolic-derived phytoalexins, are an important component elicitor-induced responses in Arabidopsis plants [35], [36]. The biosynthesis of tryptophan-derived indolic compounds was up-regulated under both treatments but was stronger induced by aphid infestation (file S4 and S5). Two of the affected modules, cellular responses to starvation and sugar-mediated signalling pathways, further indicated that both treatments resulted in cells experiencing a nutrient deficiency. Although we did not analyse cellular nutrient deficiency in the plants during our experiments, the profiles observed here are in agreement with existing information in annotation databases such as TAIR (release 10) and Gene Ontology, which are derived from the published literature.

Localised host programmed cell death is a crucial mechanism through which plants respond to pathogen and insect attack. This phenomenon regulates multiple physiological processes, including terminal differentiation, senescence, and disease resistance [37]. Several of the genes involved in the localised host programmed cell death categories were up-regulated during both treatments. These genes are also known to be induced by senescence and salicylic acid treatment, including the PR genes (PATHOGENESIS-RELATED GENE) PR1, PR2, PR4 and PR5.

Visualisation of the networks of GO terms based on the aphid-specific responses (Figure 2) and P. syringae-specific responses (Figure 3) demonstrated the massive transcriptional responses evoked in A. thaliana. Most of the significant processes were related to responses to stimuli, biosynthesis of secondary metabolites, and transcriptional and posttranscriptional regulation. Superposition of the two GO term networks generated from the aphid-specific gene list and P. syringae-specific gene-list showed significant differences in the overrepresented GO terms. The superimposed network diagram has not been included in this manuscript, but all three networks (.cys file) have been provided as additional files (files S1, S2, and S3). The interested reader can locally open these files in Cytoscape and conduct interactive exploration. (For local visualisation, download cytoscape software from http://www.cytoscape.org/, and load the.cys files on the software. Please note that the view of the annotated network presented in this manuscript has been manually simplified for representation purposes.).

thumbnail
Figure 2. Over-represented GO-categories in the aphid-specific gene list.

Network representations of enriched GO categories among genes that were differentially regulated only during the aphid experiment. Figure generated by ClueGO showing functionally grouped networks of enriched GO categories among genes whose expression was induced only in the aphid experiment. GO terms are represented as nodes based on their kappa score (≥0.3); only networks with at least three nodes are represented. The node size represents the significance of the term’s enrichment. The edges are related to the relationships between the selected terms, which are defined based on the genes that are shared in a similar way. The label for the most significant term is used as the leading group term. Visualisation was conducted using Cytoscape 2.7.0.

https://doi.org/10.1371/journal.pone.0058987.g002

thumbnail
Figure 3. Over-represented GO-categories in the P. syringae-specific gene list.

Network representations of enriched GO categories among genes that were differentially regulated only during the P. syringae experiment. Figure generated by ClueGO showing functionally grouped networks of enriched GO categories among genes whose expression was induced only in the Pseudomonas experiment. GO terms are represented as nodes based on their kappa score (≥0.3); only networks with at least three nodes are represented. The node size represents the significance of the term’s enrichment. The edges are related to the relationships between the selected terms, which are defined based on the genes that are shared in a similar way. The label for the most significant term is used as the leading group term. Visualisation was conducted using Cytoscape 2.7.0.

https://doi.org/10.1371/journal.pone.0058987.g003

Mapping the Insect- and Bacterial-specific Responses on Pathways and Processes

To structure the genes present on the A. thaliana whole-genome microarray, they were assigned to functional categories using the pathway analysis program MapMan (http://gabi.rzpd.de/projects/MapMan, version 3.5.0). MapMan is a user-driven tool that displays large datasets such as gene expression data from Arabidopsis microarrays in diagrams of metabolic pathways or other processes. After the normalisation of expression values, differential fold-change values were calculated with statistical tests, as described in the Materials and Methods section. The ratios in the 4 biological replicates were averaged and converted to a log 2 scale, then imported into MapMan as ‘.xls’ files (files S4, S5). MapMan converts the values to a false colour scale and displays them in diagrams. Transcripts that increase, decrease or change less than a given threshold are shown in blue, red and white, respectively. Some of the important categories (or functional BINs as per MapMan definition) identified via MapMan analysis are explained below.

Metabolism Overview Map

An overview of the transcriptional responses affecting genes coupled to metabolic processes showed that many genes connected to photosynthesis and energy metabolism were down-regulated after P. syringae and aphid attack (Figure 4). P. syringae infection resulted in leaf senescence and leaf yellowing, which had a major effect on chloroplast function and processes connected to the chloroplast, such as fatty acid biosynthesis, carotenoid production, chlorophyll biosynthesis, carbon fixation and others. Genes related to these processes showed clear down-regulation following P. syringae treatment. Secondary metabolism was strongly affected during both treatments, particularly regarding the phenylpropanoid and glucosinolate pathways. The results of P. syringae treatment also showed that genes connected to the terpenoid and alkaloid pathways were up-regulated, including DXPS1, TPS10, GES/TPS04, SS2, SQE6 and LAS1. In general, the stress associated with the activation and continuation of defence responses is metabolically expensive, and the plant must reallocate a significant amount of the resources that would normally be used in plant growth and reproduction to the production of defence-related compounds [38], [39]. However, in a recent work, Foyer et al. [40] explained that the decreases in growth and photosynthesis in response to stress are more likely the result of programmed down-regulation. Our experimental results showed that exposure to two biotic stresses resulted in the down-regulation of genes linked to auxin, gibberelin and cytokinin responses as well as genes coupled to cell wall modifications and cell division. The infected plants might also compensate for the depletion of sugars and amino acids, resulting in increased carbon assimilation and mobilisation of carbon, mannitol and nitrogen reserves. The plants may have degraded proteins/amino acids to generate energy (glycolysis) and re-assimilate nitrogen, through the glutamate dehydrogenase GDH2 or lysine-ketoglutarate reductase (At4g33150). There were also genes connected to starch degradation/sugar responses induced, indicating that the plants might be degrading starches, e.g., BAM5 and GPT2, used in glycolysis. Starch biosynthesis genes were generally down-regulated. The degradation of starch and maltose may also generate an osmotic force that balances water losses. During an aphid infestation, plants suffer from osmotic stress as the insect sucks large quantities of liquids from them. To counteract this situation, the transcription of genes involved in the regulation of water balance was observed to be induced, such as the WRKY40, CYP707A3 (ABA-biosynthesis), ZAT10 and ZAT7.

thumbnail
Figure 4. Metabolic overview map.

Metabolic pathways associated with the transcriptional changes affecting A. thaliana during aphid and P. syringae attack. Overview of the expression changes related to metabolic pathways observed in A. thaliana plants during the (A) aphid and (B) P. syringae treatments using MapMan software. The represented spots are only for genes showing a significant (P = 0.01) change in expression between the treatment and the untreated control that were attributed to the respective bins by MapMan. Genes whose expression levels were increased are indicated with an increasingly blue colour, while decreasing expression is indicated in red. The graduation can be seen on the scale presented in the top right corner of each subfigure. A change in expression of log2 = 2.0 scale was selected as giving full saturation.

https://doi.org/10.1371/journal.pone.0058987.g004

Comparative Overview of the Response to Biotic Stress during the Aphid and P. syringae Treatments

A plant’s reaction to biotic stress involves several steps: after the initial signal input from the pathogen, which is recognised by the corresponding receptors (putative R genes), transcription of the cascade associated with the plant defence mechanism is triggered, including changes related to oxidative stress. Inside the cell, signals are transmitted and lead to the production of defence molecules (PR proteins, heat shock proteins and secondary metabolites). A large number of signalling genes were activated during both the aphid and P. syringae treatments (Figure 5). Most of these genes encode receptor kinases, leucine-rich receptor kinases, MAP kinases, calcium-binding proteins and proteins regulating oxidative stress, such as peroxidases (details in file S6). The number of signalling proteins that were differentially expressed during the aphid experiment was more than four times higher compared to the P. syringae treatment. There were 278 aphid-specific signalling genes, but only 62 P. syringae-specific signalling genes. Thirty-one heat shock proteins were differentially expressed only during the aphid treatment (file S7), the majority of which were of the DnaJ/Hsp40 type chaperones and were induced. Large numbers of proteolytic enzymes were differentially expressed during both the aphid (220) and P. syringae (89) treatments. The majority of these enzymes were ubiquitin proteases, F-box proteins, cysteine proteases, serine proteases, C3HC4-type RING fingers, and metalloproteases (file S8). Several of the down-regulated proteolytic enzymes were chloroplast localised or predicted to be located in the plastid/chloroplast, while most of the C3HC4-type RING finger proteins were induced. Secondary metabolites play a crucial role during plant defences. Sixty-three genes related to secondary metabolic processes were differentially regulated during the aphid and P. syringae treatments. Some of these secondary processes include the biosynthesis of isoprenoids, phenylpropanoids, glucosinolates and flavonoids. A detailed analysis of the differentially regulated secondary metabolic processes can be found in file S9 and in a later section of this article. There were 76 differentially regulated genes connected to cell wall-related processes identified during the aphid treatment, but only 36 in the P. syringae experiment. These genes included components involved in cell wall precursor synthesis, cellulose synthases, cell wall structural proteins such as AGPs (arabinogalactan protein), LRR (leucine-rich repeat) extensin-like proteins, and HRGPs (hydroxyproline-rich glycoproteins) (details in file S10). In general, aphid attack appeared to affect cell wall-related processes to a greater extent than P. syringae infection. In particular, a large number of APGs and xyloglucan:xyloglucosyl transferases were observed to be down-regulated. Xyloglucan endotransglucosylases are known to play an important role during cell elongation and cell wall modifications during shade avoidance [41], [42]. The effects observed on genes encoding pathogenesis-related proteins (PR proteins) showed a clear bias between the treatments: while only 11 P. syringae-specific PR proteins were affected, 56 genes encoding aphid-specific PR proteins showed differential expression. The PR proteins include a wide variety of protein types, such as ß-1,3-glucαnases, chitinases, thaumatin-like protein, proteinase inhibitors, plant defensins and others. The PR1 protein, which is often used as a marker for salicylic acid responses, was more than ten-fold higher induced by the aphid attack than by P. syringae infection. Another class of proteins that was induced and significantly overrepresented after aphid attack corresponded to a large number of disease resistance proteins belonging to the TIR-NBS-LRR (Toll/Interleukin1 receptor–nucleotide binding site–leucine-rich repeat) proteins.

thumbnail
Figure 5. Biotic stress response overview map.

This figure shows the changes in the expression of biotic stress-responsive genes in A. thaliana plants during the response to the aphid and P. syringae treatments. Genes that have been experimentally indicated to be involved in biotic stress are collected in the main panel (coloured with dark grey), while genes and pathways that are putatively involved in biotic stress pathways are shown on the left and right sides (coloured in light grey). (A) Aphid infestation. (B) P. syringae infection. In both cases, the signal after infection is expressed as a ratio relative to the signal in uninfected controls, which was converted to a log2 scale and displayed. The scale is shown in the figures. Only the genes showing a significant (P = 0.01) change in expression between the treatment and the untreated control that were attributed to the respective bins by MapMan are shown. Genes whose expression was increased are indicated with increasingly intense blue and red colours. The gradation can be seen in the scale presented in the top right corner of each subfigure.

https://doi.org/10.1371/journal.pone.0058987.g005

Among the biotic stress-related transcription factors, some WRKY and bZIP proteins were expressed differentially only during the aphid experiment, while some MYB proteins were expressed differentially only during P. syringae infection. Other differentially regulated classes of transcription factors included ERF/AP2, NAC, bHLH and DOF. Details regarding the differentially regulated transcription factors are provided in a separate section of this article. The plant defence responses associated with P. syringae and aphid attack induced and repressed various hormonal signalling pathways. The most affected of these pathways during our experiments were the JA, SA, ABA, ethylene and auxin pathways. Among the hormonal signalling pathways, some components of the ethylene, JA, SA, ABA, auxin and brassinosteroid pathways appeared to specifically be regulated during the aphid and P. syringae treatments. There were relatively few ethylene responses observed in general, but such effects were clearly stronger after the aphid than the Pseudomonas treatment. Examples of ethylene responses included ACS6, ERF11, which may modulate ABA-regulated ethylene biosynthesis, ORA59, which integrates JA and ethylene signals during plant defence, EFE (ethylene forming enzyme) and ATARD3 (methionine recycling during ethylene synthesis). Some proteins involved in the biosynthesis of ethylene were also affected. JA was more strongly induced by P. syringae, but the SA response was stronger following aphid attack. The details of the differentially regulated genes involved in hormone-mediated signalling pathways are provided in a separate section of this article.

Regulatory Overview Map

The categories that included most of the induced regulatory genes were TFs, receptor kinases, protein degradation and protein modification. In addition, several genes involved in overrepresented induced biological processes, such as the auxin signalling pathway and autophagy, were included in the regulatory categories (Figure 6). A comparative list of the differentially expressed genes (both aphid-specific and P. syringae-specific genes) involved in hormonal pathways and their corresponding log2-transformed expression values are provided in file S11. The ethylene pathway was up-regulated after the aphid treatment, and genes such as ACS6, ATARD3, ERS1 and a number of ethylene responsive element binding factors (ERFs) were induced. Genes belonging to the ERF/AP2 family are induced by many biotic and abiotic factors, among which ethylene ERFs not only control a subset of ET-mediated responses but might also integrate ET with other signalling pathways. Increased ethylene production is a common defence response after herbivore attack and has been reported in several plant species [43]. Both ABA and JA responses were up-regulated by the P. syringae treatment, but few known SA-responsive genes were induced. Two categories, receptor kinases and calcium regulation (in Figure 6), appeared to be quite highly represented according to the MapMan annotation during the aphid experiment. Nevertheless, two other categories, light signalling and redox control, included fewer transcripts and gene families, respectively. These were some key differences between the aphid and P. syringae treatment. Genes encoding receptor kinases and proteins coupled to calcium signalling were overrepresented following the aphid treatment. These genes include a large number of cysteine-rich receptor-like protein kinases, such as, CRK7, CRK37, CRK36, CRK23, CRK14, CRK11, CRK15, CRK6, CRK28 and calmodulin-like proteins (CML40, CML47, CML11, TCH3, TCH2/CML24, CML44, CML45, CML30, CML37, CML38 and others) as well as several calmodulin-binding IQ-domain proteins. Together with the MAP kinases (MPK11, MKK9, ATMPK3, MEKK3, MEK1, MEKK1, MKK2, MKK4, MPK4 and others), they constitute a large network that activates various plant defence responses, resulting in the activation of key transcription factors.

thumbnail
Figure 6. Regulatory overview map.

MapMan regulatory overview map showing differences in transcript levels between aphid-specific and P. syringae-specific genes. Aphid-specific and P. syringae-specific bins are marked as ‘A’, and P. syringae-specific bins are marked as ‘P’. In the colour scale, blue represents higher gene expression, and red represents lower gene expression. IAA, Indole-3-acetic acid; ABA, abscisic acid; BA, brassinosteroid; SA, salicylic acid; MAP, mitogen-activated protein.

https://doi.org/10.1371/journal.pone.0058987.g006

Differences Observed in the Jasmonic Acid (JA) Biosynthesis Pathway during the Aphid and P. syringae Treatments

The jasmonic acid signalling pathway is a highly conserved, powerful regulator of plant defence signalling that is activated during infection by various pathogenic microorganisms as well as upon insect attack [44]. Kuśnierczyk et al. reported that more than 200 genes are dependent on the plant’s jasmonate status, irrespective of external stimuli, and that the aphid-induced response of more than 800 transcripts is regulated by jasmonate signalling [45]. The release of linolenic acid from membrane lipids initiates a series of enzymatic reactions known as the octadecanoid pathway, leading to accumulation of JA and related compounds. Additionally, 12-oxophytodieonic acid (OPDA) is a biosynthetic precursor of JA signalling molecules, which activate the expression of related-related genes. The selection of transcripts induced by JA and OPDA varies to some extent. This difference can be attributed to the electrophilic activities of the cyclopentanone ring of JA [46]. A number of enzymes coupled to oxylipin/JA biosynthesis, such as AOC3, OPR3, OPCL1, LOX2 and LOX3, were up-regulated by both treatments, while AOS, AOC1, AOC2, AOC4, ACX1, ACX5 and LOX1 were mainly induced by Pseudomonas. Two OPR-related genes, At1g18020 and At1g17990, as well as the lipoxygenases LOX4, LOX5 and LOX6 were only induced by aphid attack. Almost none of the genes encoding proteins potentially linked to oxylipin biosynthesis were down-regulated, with the exception of OPR1, which was down-regulated by P. syringae infection.

SA Regulates the Expression of Aphid-specific Defence Proteins, and Methyl Salicylate Activates P. syringae-specific Defence Proteins

Salicylic acid is another stimulator of plant defence responses and is an important trigger of systemic acquired resistance (SAR), resulting in increased defence against a variety of pathogens. Methyl salicylate (MeSA) has been identified as one of the mobile signals required for SAR. MeSA is translocated from the site of infection through the vascular system to distal (systemic) tissues, where it activates specific defence responses. The SAR response results in a complex chain of events and is regulated by various transcription factors. In higher plants, SA can be synthesised from phenylalanine via cinamic acid or from isochorismate. During pathogen attack, SA signalling leads to accumulation of various pathogenesis-related proteins (PR proteins), which can possess antimicrobial and anti-insect activities. Interestingly, MeSA released by the attacked plants can be detected by insects and changes their plant preferences [47]. In our analysis, expression of a methyltransferase gene (At3g21950) related to salicylate O-methyltransferases was down-regulated during aphid treatment. At3g21950 encodes a S-adenosyl-L-methionine:salicylic acid carboxyl methyltransferase related to BSMT1 that may convert SA to MeSA. In contrast, other methyltransferase genes BSMT1 (At3g11480), which converts SA to MeSA, and UGT74E2 (At1g05680) were up-regulated during P. syringae treatment (file S11). UGT74E2 is hydrogen peroxide responsive and may be involved in water stress responses. There were relatively few known genes coupled to the biosynthesis of SA found in both datasets. However, BSMT1 might be a key enzyme.

Although relatively few genes connected to the biosynthesis of SA and MeSA were found in the obtained datasets, several genes induced by SA were identified. Additionally, isochorismate synthase 1 (ICS1) and one of its transcriptional regulators, WRKY46, were induced by aphid infestation. Another gene induced after aphid treatment that may be under the regulation of WRKY46 is PBS3. PBS3 most likely encodes an enzyme producing SA-glucoside, a putative storage form of SA, and pbs3 mutant plants exhibit impaired activation of defence genes such as PR1. The PR1 gene, a common marker for SA-induced genes, was strongly up-regulated by the aphid treatment (log2 = 5.5) and slightly less induced by P. syringae treatment (log2 = 1.7). The WRKY53 gene, which is known to be up-regulated by SA [48], was only induced in the aphid treatment. A number of genes, such as ALD1 and BAP1, coupled to systemic defence responses were uniquely induced by aphids.

Overview of Differences in Secondary Metabolism

Plants have evolved many secondary metabolites involved in plant defence, which are collectively known as antiherbivory compounds and can be classified into three sub-groups: nitrogen compounds (including alkaloids, cyanogenic glycosides and glucosinolates), terpenoids, and phenolics [49]. In addition to the three larger groups of substances mentioned above, fatty acid derivatives, amino acids and even peptides are used in defence. The terpene synthase genes GES (geranyllinalool synthase, At1g61120), LAS1 (Lanosterol synthase, At3g45130), and TPS10 (Terpene Synthase 10, At2g24210) were highly up-regulated during P. syringae treatment. The elicitor-activated gene CAD-B2 (At4g37990), belonging to the phenylpropanoid metabolism category, was strongly up-regulated in a P. syringae-specific manner. In the alkaloid-like compound biosynthesis category, strictosidine synthase genes (At1g74010, At1g74020) were highly up-regulated in the P. syringae experiment. In the flavonoids category, two genes SRG1(Senescence-Related Gene 1; At1g17020) and 2-oxoacid-dependent oxidas (At3g50210), were also up-regulated in a P. syringae-specific manner. Significant differences were observed in genes coupled to glucosinolate biosynthesis and hydrolysis. Glucosinolates (GS) are secondary metabolites typical of the order Brassicales [50]. Most of the aphid-specifically expressed aliphatic GS genes were repressed, whereas most of the Pseudomonas-specifically expressed genes were positively regulated. The lists of these genes are provided in Tables 2 and 3. Following P. syringae treatment, two myrosinase-associated proteins (At1g52040, At1g54020) and a nitrile-specific protein AtNSP5 (At5g48180) were highly up-regulated. It was reported by Kissen et al., that the nitrile specifier proteins involved in glucosinolate hydrolysis in Arabidopsis thaliana and products generated after hydrolysis, such as isothiocyanates, play multiple roles in growth regulation and defence [51].

thumbnail
Table 2. Genes involved in glucosinolate metabolism affected by aphid infestation.

https://doi.org/10.1371/journal.pone.0058987.t002

thumbnail
Table 3. Genes involved in glucosinolate metabolism affected by P. syringae infection, with log2 fold-change values.

https://doi.org/10.1371/journal.pone.0058987.t003

A Large Number of Transcription Factors are Differentially Regulated, Many of which are Unique to Insect or Bacterial Stress

Transcription factors are the key regulators of gene expression changes and, thus, represent important part of a complex regulatory network allowing plants to adjust to changes in their environment [52]. Members of several Arabidopsis transcription factor families have been linked to plant stress responses, and a significant overlap in the expression profiles of many of these genes corresponding to a range of stress conditions has been reported. TFs are often induced by signalling phytohormones such as JA, SA or ET. The TFs that were differentially expressed during the aphid and Pseudomonas treatments are reported in Table 1, and their names are given in file S12. Additionally, pictorial representations of the aphid-specific and P. syringae-specific TFs produced using MapMan software are shown in Figure 7. There were 16 WRKY TFs that were up-regulated in an aphid-specific manner (WRKY20, WRKY22, WRKY39, WRKY21, WRKY40, WRKY26, WRKY50, WRKY25, WRKY38, WRKY51, WRKY53, WRKY47, WRKY46, WRKY69, WRKY33, WRKY16). WRKY TFs can act as both positive and negative regulators of plant defence pathways. The mechanisms activating WRKY TFs can involve the MAP kinase cascade and calcium signalling. It has been demonstrated that a subgroup of WRKY TFs can act as calcium concentration sensors, being activated by the increase in the Ca2+ concentration that occurs under inducer attack [53]. Mechanical penetration of cells by aphid stylets changes the plasma membrane potential and increases in intracellular Ca2+ concentrations. Fluctuations in the cytosolic Ca2+ concentration resulting from the opening of membrane-bound calcium channels are further decoded by several Ca2+-binding proteins, including the WRKY TFs. The up-regulated aphid-specific TFs also include C2H2 zinc finger proteins.

thumbnail
Figure 7. Transcription overview map.

(A) Aphid specific; (B) P. syringae specific.

https://doi.org/10.1371/journal.pone.0058987.g007

The MYB family, which is another large family of TFs characterised by a conserved MYB DNA-binding domain, bind to a variety of different DNA sequences. Among the P. syringae-specific TFs, there are 9 MYBs (MYB95, MYB112, MYB90, MYB102, MYB32, MYB114, MYB59, MYB60, MYB20), of which the first 7 are significantly up-regulated. MYB90 is also known as PAP2 (PRODUCTION OF ANTHOCYANIN PIGMENT 2), suggesting that some of these MYBs are likely to be involved in anthocyanin biosynthesis. A few members of this family specifically activate genes related to tryptophan aliphatic glucosinolate and indoyl glucosinolate synthesis [54].

Integrated Information from Available Public Domains Reveals a Pattern of Potential MicroRNA-mediated Post-transcriptional Regulation during Insect and Bacterial Attack

We constructed a genetic network of the differentially regulated gene lists using the Gene network tool in VirtualPlant. First, individual genes belonging to the common category were grouped into a “super node” based on shared functional properties, such as GO terms, KEGG pathways, Gene families and even similar annotations. The functional annotations were categorised in a hierarchical manner, where the functional terms and pathways were themselves grouped into higher, more generic categories (details are given in the Materials and Methods). During this analysis, we used post-transcriptional regulation, protein-protein interactions, and transcriptional regulation information from both experimental and predicted databases. The ‘Regulated Edges’ are predicted interactions based on the presence of known transcription factor cis-acting binding sites located in the 3 kbp upstream region of annotated transcripts. Interestingly, some of the key stress-regulated transcription factors are reported in publications, or have been computationally predicted to be regulated by different microRNAs. Thus, we were able to hypothesise that the activation of microRNA genes under biotic stresses would lead to the repression of many downstream protein-coding genes and affect physiological responses. This analysis indicates a new direction for conducting large-scale experiments and subsequent bioinformatics analyses to explore the regulatory links between biotic stress and microRNAs in A. thaliana.

Connection of microRNAs to Genes from the Common Category

Supernode analysis of the differentially expressed common genes using the VirtualPlant tool revealed a supernode, or cluster, of 66 genes known to show connections with 27 microRNAs (Figure 8A, marked as a blue-coloured cluster). Further analysis on this cluster of 66 genes identified 9 genes (Table 4) with experimentally validated microRNA-binding sites (Figure 8B). Six of these genes encode known Arabidopsis transcription factors. Manually retrieved related literature references for each of this microRNA are provided in Table 5. We then studied all of the microRNA genes curated in the microRNA Registry database (microrna.sanger.ac.uk/sequences/). Out of the 66 genes in this cluster (Figure 8A, supernode annotated as ‘none’), At1g20510 (OPCL1) and At4g05160 (putative 4-coumarate-CoA ligase/4-coumaroyl-CoA synthase) are known to be involved in jasmonic acid biosynthetic processes. Two genes, At1g50670 and At5g53160 (SPL4), showed the maximum number of connections to microRNAs. RD26 (RESPONSIVE TO DESSICATION 26, At4g27410), BZIP25 (BASIC LEUCINE ZIPPER 25, At3g54620), JIN1 (JASMONATE INSENSITIVE 1, At1G32640) and BGLU11 (BETA GLUCOSIDASE 11, At1g02850) contain putative microRNA binding sites, though they have not yet been verified experimentally. Out of these 13 genes, 6 are known to be TFs. Details are given in Tables 4 and 5.

thumbnail
Figure 8. Retrieved micro-RNA connections of the common genes.

A) Super node analysis using the Gene networks tool in VirtualPlant, visualised with Cytoscape 2.7.0. Individual genes in the common category were grouped into a supernode (red-coloured nodes) based on shared functional properties, such as GO terms, KEGG pathways, gene families and even similar annotations. Each supernode size corresponds to the number of genes present in that category. The edges represent connections among different functionally grouped supernodes. The top 6 most highly populated supernodes are filled with green colour. A supernode consisting of 66 genes known to show connections with 27 microRNAs (cluster of blue-coloured nodes). microRNA binding sites have been reported in existing literature for 9 of these genes, and 6 of them are known transcription factors. B) Details of the 9 genes mentioned above, which are known to be regulated by 27 microRNAs. MicroRNAs are shown as green-coloured circles, whereas target genes are depicted as red-coloured triangles. Edges represent the interactions between microRNAs and their target genes. Please also refer to Table 4 and Table 5 for detailed information and related evidence in the literature.

https://doi.org/10.1371/journal.pone.0058987.g008

thumbnail
Table 4. The 9 genes in the common set known to be regulated by biotic stresses and their association with stress-inducible microRNAs (Refer to Figure 8 B).

https://doi.org/10.1371/journal.pone.0058987.t004

thumbnail
Table 5. Functional targets of the microRNA families in the common set of genes (retrieved from literature searches).

https://doi.org/10.1371/journal.pone.0058987.t005

Connection of microRNAs to Genes Showing Aphid-specific Responses

Among the 3,382 transcripts showing an aphid-specific response, the GO enrichment category ‘response to stimuli (biotic and abiotic stress)’ included 242 stress-responsive genes. Among these genes, 42 are known to exhibit transcription factor activity (Figure 9A). Additionally, out of these 242 stress-regulated genes, 21 genes have been reported to be associated with microRNAs based on literature and database searches, as described in the Materials and Methods section (Table 6). The reported target gene families for these microRNAs were retrieved through a manual literature search and are listed in Table 7. Many of the genes that were differentially regulated by aphid attack belong to these reported gene families.

thumbnail
Figure 9. Retrieved micro-RNA connections of aphid-specific and P. syringae-specific genes. A red triangle represents a target gene, and a green circle represents a microRNA.

A) Among the transcripts showing aphid-specific responses, 42 genes are known to contain microRNA binding sites. Please also refer to Table 6 and Table 7 for detailed information and related evidence from the literature. B) Among the transcripts showing P. syringae-specific responses, 9 genes are known to contain validated microRNA binding sites. We were able to find related references in the literature for the reported 23 microRNAs. Please also refer to Table 8 and Table 9 for detailed information and related evidence from the literature.

https://doi.org/10.1371/journal.pone.0058987.g009

thumbnail
Table 6. The 21 genes in the aphid-specific gene set known to be regulated by biotic stresses and their association with stress-inducible microRNAs.

https://doi.org/10.1371/journal.pone.0058987.t006

thumbnail
Table 7. Functional targets of the microRNA families in the aphid-specific set of genes (retrieved from the existing literature).

https://doi.org/10.1371/journal.pone.0058987.t007

Connection of microRNAs to Genes Showing P. syringae-specific Responses

Among the 1602 transcripts showing P. syringae-specific responses, the GO enrichment category ‘response to stimuli (biotic and abiotic stress)’ included 146 genes. Out of these 146 stress-responsive genes, 24 are known to exhibit transcription factor activity (Figure 9B), and 6 have been reported to be associated with microRNAs based on literature and database searches (Table 8). The reported target gene families for these microRNAs were retrieved through a manual literature search and are listed in Table 9. Many of the genes that were differentially regulated by Pseudomonas attack belong to these reported gene families.

thumbnail
Table 8. The 6 genes in the Pseudomonas-specific gene set known to be regulated by biotic stresses and their association with stress-inducible microRNAs.

https://doi.org/10.1371/journal.pone.0058987.t008

thumbnail
Table 9. Functional targets of the microRNA families in the P. syringae-specific set of genes (retrieved from the existing literature).

https://doi.org/10.1371/journal.pone.0058987.t009

Cross-validation of Differentially Regulated Aphid and Pseudomonas-specific Transcription Factors via Co-expression Analysis of the Multiple Biotic Stress Dataset

The differentially regulated gene sets included many signature transcription factors known for their involvement in stress responses. A co-expression analysis based on a compendium of 69 ATH1 biotic stress experiments, generated using the CORNET tool, showed that many of these TFs have been found to be strongly co-expressed during various biotic stress experiments. From the 66-gene supernode cluster in the common group, the co-expression analysis produced a network of 26 nodes with 25 edges (Figure 10A). One module consisted of 9 genes CPK6, TCH3, BZIP25, AOX1D, RD26, ERD2, MPK1, GDH2 and HSF4 that were strongly co-expressed. The extended module contained 16 genes, 5 of which are involved in calcium-mediated signalling: CPK6, TCH2, TCH3 and two EF-hand proteins. Functional annotation revealed that these genes are known to be involved in several different biotic and abiotic stress responsive processes. The calcium-dependent protein kinase (CPK6) is a positive regulator of methyl jasmonate signalling in guard cells and represents an important gene involved in methyl jasmonate signalling and signal crosstalk between methyl jasmonate and abscisic acid in guard cells [55]. TCH3 is a calmodulin-like protein that is up-regulated in response to various environmental stimuli, including mechanical stimuli [56]. Responsive to desiccation 26 (RD26) encodes an NAC transcription factor that may be coupled to an ABA-dependent stress-signalling pathway [57], while the heat shock protein-70 cognate protein Early-responsive to dehydration (ERD2) which is induced by heat and dehydration is a key element in defence response signalling pathways [58]. The MAP-kinase gene MPK1 participates in pathogen signalling, and its kinase activity increases in response to mechanical injury [59]. Glutamate dehydrogenase 2 (GDH2), the alpha-subunit of glutamate dehydrogenase, is a mitochondrial protein that has been reported to be responsive to diverse environmental stresses [60]. Arabidopsis heat shock factor (HSF4) regulates the expression of heat shock proteins [61]. The genes in the aphid-specific and pseudomonas-specific co-expression module have been discussed in previous sections.

thumbnail
Figure 10. Co-expression network.

Co-expression networks generated by CORNET using AtGenExpress biotic stress compendia based on a Pearson’s correlation coefficient threshold ≥0.7. The networks were visualised using Cytoscape2.7.0. Pink-coloured edges represent a strong correlation of ≥0.9, and cyan-coloured edges represent a correlation of ≥0.7 to 0.9. A) Co-expression network analysis among the 66-supernode cluster in the common group resulted in a network of 26 nodes 25 edges. B) Co-expression network analysis of the aphid-specific TFs resulted in a network of 24 tightly co-expressed TF modules. C) Co-expression network analysis among 104 Pseudomonas-specific TFs resulted in a tightly co-expressed modular network consisting of 55 nodes and 94 edges.

https://doi.org/10.1371/journal.pone.0058987.g010

Conclusions

We generated and analysed data from two different biotic stress experiments conducted in Arabidopsis thaliana in which the plants were challenged with the aphid Brevicoryne brassicae and the bacterium P. syringae syringae. Our data showed that the transcriptional response of Arabidopsis to these very different attackers resulted in the differential regulation of a diverse range of biological processes. Transcriptional responses and networks unique to insect or bacterial stress conditions were identified, as were sets of genes showing similar a response under both stresses. By examining the responding genes and the functional network characteristics of each stress response, we found that a significant number of the transcripts encode transcription factors. Most of these transcription factors have shown to be involved in stress responses and regulatory processes. Some WRKY and bZIP genes were expressed differentially only during the aphid experiment, whereas some MYB genes were expressed differentially only during P. syringae infection. A Gene Ontology-based overrepresentation analysis revealed that half of the genes from the common list were involved in central metabolic and cellular processes, such as electron transport and energy pathways localised to the plastid. Secondary metabolism was strongly affected during both treatments, particularly the phenylpropanoid and glucosinolate pathways. Processes connected to the chloroplast, such as fatty acid biosynthesis, carotenoid production, chlorophyll biosynthesis, carbon fixation and others were down-regulated following P. syringae treatment. Starch biosynthesis genes were generally down-regulated, and an indication was found that the plants were degrading starch, which could help the plants to maintain the osmotic balance. Components of the ethylene, JA, SA, ABA, auxin and brassinosteroid pathways appeared to be specifically regulated during the aphid and P. syringae treatments. Ethylene responses were clearly induced during aphid feeding, while JA was more strongly induced by Pseudomonas. The number of signalling proteins that were differentially expressed during the aphid experiment was more than four times higher compared to the P. syringae treatment. By integrating secondary information from most available public sources, we further explored the regulatory links between biotic stress and microRNAs associated with aphid-and P. syringae -specific differentially regulated processes in A. thaliana, and the corresponding genes are briefly summarised in Table 10.

thumbnail
Table 10. Summary of aphid-specific and P. syringae-specific genes associated with differentially regulated processes during both of the treatments.

https://doi.org/10.1371/journal.pone.0058987.t010

This study therefore demonstrates that the integration of heterogeneous publicly available information from multiple databases with experimental results can help plant biologists develop a better understanding of stress-associated processes in plants. Due to logistics and costs we examined only a single time point during the A. thaliana (Col-0) - P. syringae treatment. We are fully aware that comparing single time point restricts some analyses and is a potential limiting factor as demonstrated by Bricchi et al (2012) [62]. Although several datasets reporting temporal responses of A.thaliana to P. syringae infection were available from previously published independent studies [33], [63], [64], [65], we decided not to combine them in the current analysis while making comparisons with our own B. brassicae data [66] to maintain the homogeneity of the comparisons. The analysis presented here will therefore not explain the comparative temporal dynamics of A. thalianaB. brassicae and A. thaliana – P. syringae interactions.

Materials and Methods

To overcome the problem of the incompatibility of independent microarray experiments, a genome-wide expression analysis involving 2 different biotic stresses was conducted, in which Arabidopsis thaliana plants were infested with aphids (Brevicoryne brassicae) [45] or infected with P. syringae bacteria (4 biological replicates, and an untreated control were used for each comparison). The microarray data from the aphid experiment was part of a larger plant-insect study [45]. The Pseudomonas data were generated for the present study using the same technology platform to reduce experimental variation. All data have been deposited in GEO (GSE39245 and GSE39246). A systems biology approach was followed to understand common and specific responses in terms of different pathways and processes in Arabidopsis during insect and bacterial attack. A simplified flow chart diagram of the applied methodology is provided in Figure 11.

Plant Material and Cultivation

The Arabidopsis thaliana Columbia-0 ecotype (Col-0) plants used in the experiment were derived from seeds produced by Lehle Seeds (Round Rock, USA; Catalog No. WT-2-8, Seed Lot No. GH195-1). The seeds were sterilised according to standard procedures and grown on agar medium containing an MS basal salt mixture (Sigma), 3% (v/w) sucrose, and 0.7% (v/w) agar (pH 5.7) to assure uniform germination. After 15 days, the seedlings were transferred to 6 cm diameter pots (3 seedlings per pot) filled with a sterile soil mix (1.0 part soil and 0.5 parts horticultural perlite). The plants were kept in Vötsch VB 1514 growth chambers (Vötsch Industrietechnik GmbH, Germany) under the following conditions: 8 h/16 h (light/dark) photoperiod, 22°C/18°C, 40%/70% relative humidity, and 70/0 µmol m−2s−1 light intensity. A short day length was applied to prevent the plants from bolting.

Infestation Experiments

At 32 days of age (17 days after being transferred to soil), the plants had 8 fully developed leaves. Each plant was infested with 32 wingless aphids (4 per leaf), which were transferred to the leaves with a fine paintbrush. Infested plants and aphid-free controls were maintained in Plexi-glass cylinders, as described previously [66]. The plants were harvested 72 h after infestation between the 6th and 8th hours of the light photoperiod. Four biological replicates were produced from the control and infested plants, with each being sampled from 15 individual plants. Whole rosettes were cut at the hypocotyl, and aphids were removed by washing with Milli-Q-filtered water. The harvested material was immediately frozen in liquid nitrogen.

A P. syringae culture was grown overnight in 10 ml of Kings B solution (King et al., 1954) supplemented with the antibiotics rifampicin (50 µg ml−1) and kanamycin (25 µg ml−1). The overnight culture was washed once in 10 mm MgCl2, and the final cell densities were adjusted to an OD of approximately 0.20 at 600 nm (approximately 1.5×108 cfu ml−1) in 10 mm MgCl2. Plants were grown as described in the Plant material and cultivation section. Then, 30-dayold plants were mock-challenged with 10 mm MgCl2 or inoculated with the DC3000 strain of P. syringae by infiltrating 3–4 leaves on the abaxial surface with a needleless 1 ml syringe. Four biological replicates of infested leaves and leaves obtained from control plants grown under identical conditions were harvested after 3 days (between the 6th and 8th hours of the light photoperiod). The leaf material was immediately frozen in liquid nitrogen. Leaves from 15 plants were included in each replicate.

RNA Isolation, cDNA Synthesis, Labelling and Hybridisation

Total RNA was isolated from cauline leaf tissue from plants from each experiment. Each experiment consisted of four infested samples and four control samples. Total RNA was extracted from 100 mg of cauline leaf material using the RNeasy Plant Minikit (Qiagen, Hilden, Germany) and eluted in 2×50 µl of RNAse-free water. Any residual DNA in the RNA samples was removed by on-column treatment with RNAse-free DNase. The eluted RNA was concentrated to 10–20 µl using a 30 kDa cut-off Microcon spin filter unit (Amicon, Bedford, USA). To protect the RNA from degradation, the RNasin Plus RNase inhibitor (Promega, Madison, USA) was added to a final concentration of 1 unit µl-1. The purity and quantity of the obtained RNA was determined using a Nanodrop ND 1000 instrument (Nanodrop Technologies, Wilmington, DE, USA). RNA integrity was analysed via formaldehyde agarose gel electrophoresis. First-strand cDNA was generated from total RNA (15 µg) using the Superscript III Reverse Transcriptase (Invitrogen, Carlsbad, USA) and oligo dT primers with a 3DNA capture sequence from the 3DNA Array 350TM kit (Genisphere, Hatfield, USA). RNA samples were labelled with either the Cy5-capture primer or Cy3-capture primers (sample dye-swapping). The cDNAs were hybridised to the microarray slides at 58°C using a Slide Booster Hybridisation Station (Advalytix, Brunnthal, Germany) together with Cy3- and Cy5-labelled dendrimers from Genisphere. The slides were washed according to the manufacturers’ instructions (Genisphere and Advalytix).

Microarrays

The microarray slides contained 31811 unique 70-mer oligos with a C6-amino linker, corresponding to a total of 33696 spots, covering 26624 genes. Of these oligos, 29110 were from the Qiagen-Operon Arabidopsis Genome Array Ready Oligo Set (AROS), Version 3.0, while the others were custom made and produced by Operon (Alameda, CA, USA) or MWG (Ebersberg, Germany). The sequences of all of the custom-made probes on the chip have been deposited in GEO and are available under accession GPL15699. The oligonucleotides were dissolved in MQ grade water and 50% DMSO (20 pmol/µl) and spotted on aminosilane-coated UltraGaps slides (Corning, NY, USA) using a BioRobotics MicroGrid II robot (Genomic Solutions, MI, USA). Printing of the microarray slides was performed at the Norwegian Microarray Consortium (Trondheim, Norway). Hybridisations were conducted using a Slide Booster Hybridization Station (Advalytix, Brunnthal, Germany), and the slides were washed according to the manufacturer’s instructions (Genisphere and Advalytix). The slides were scanned at a 10 µm resolution on a G2505B Agilent DNA microarray scanner (Agilent Technologies). The resulting images were processed using GenePix 5.1 software (Axon Instruments, Union City, USA).

Statistical Analysis of the Microarray Data

Each dataset obtained from the aphid and Pseudomonas treatments corresponded to 4 microarray slides, where the controls and treated samples were alternately labelled with Cy5 and Cy3. The GenePix-processed data were filtered to remove spots that had been flagged as ‘Absent’, ‘Not Found’ or ‘Bad’, or exhibited median foreground intensity below the local median background intensity. The R statistical program (version 2.10.1) was used for all statistical analyses [67]. No background subtraction was performed. The data from each array were log-transformed and normalised using the printtip-loess approach (Yang et al. 2001). Within-array replicated measurements for the same gene were merged by taking the average over the replicates. The data were then scaled so that all array datasets presented the same median absolute deviation. Genes showing dye-biased responses due to Cy5 and Cy3 labelling were identified and excluded. During data processing, we focused on genes for which at least 3 out of 4 biological replicates for the examined time points passed the quality control criteria suggested by Jørstad et al. [68], [69]. To make statistical inferences about differentially regulated genes, the Limma package [70] was used. The Limma approach is based on fitting a linear model to the expression data from each probe on a microarray. Genes showing an adjusted p-value of less than 0.05 were considered to be significantly differentially expressed. All of the genes discussed in this paper were found to be significantly differentially expressed in one of the two treatments (aphid or Pseudomonas).

GO Enrichment Analysis of Common Genes

We employed a simple set theory-based operation in R to find common and specific transcriptional responses that occurred in both experiments. To conduct automated GO [71], TAIR [72] annotations, we simultaneously used three programs: ClueGO [73], BiNGO [74] and VirtualPlant [75]. Only the ClueGO results were included in this manuscript. Transcription factors were classified according to the ‘The Database of Arabidopsis Transcription Factors’ [76]. In ClueGO, to calculate enrichment values for terms and groups, we used two-sided (enrichment/depletion) tests based on the hypergeometric distribution to calculate doubling for two-sided tests to address discreetness and conservatism effects, as suggested by Rivals et al. [77]. To correct the P-values for multiple testing, the Bonferroni method was used to control the type I error (false positive) rate [78]. ClueGO employs a new kappa statistic. To link the terms in the network, ClueGO first creates a binary gene-term matrix with the selected terms and their associated genes. Based on this matrix, a term–term similarity matrix is calculated using chance-corrected kappa statistics to determine the strength of the associations between the terms. Because the term–term matrix is of categorical origin, using a kappa statistic was found to be the most suitable method. Finally, the created network represents the terms as nodes, which are linked based on a predefined kappa score level. The kappa score threshold can initially be adjusted on a positive scale from 0 to 1 to restrict the network connectivity in a customised way. In our analysis, we used a kappaScore threshold of 0.3. The size of the nodes reflects the enrichment significance of the terms. The functional groups represented by their most significant (leading) term are visualised in the network, providing an insightful view of their interrelationships. Furthermore, other ways of selecting the group-leading term, e.g., based on the number or percentage of genes per term, are also provided.

VirtualPlant [75] integrates genome-wide data regarding the known and predicted relationships among genes, proteins, and molecules as well as genome-scale experimental measurements. This warehouse includes descriptions of molecular entities (e.g., gene annotations and functional classifications), molecular interactions (metabolic associations, regulatory interactions, and other interaction data from public databases), and publicly available microarray data (including more than 1,800 gene chip hybridisations from the ATH1 Affymetrix platform obtained from the European Arabidopsis Stock Center [NASC] using the Affywatch subscription service). VirtualPlant also provides visualisation techniques that render multivariate information in visual formats that facilitate the extraction of biological concepts.

Co-expression Analysis of Common Genes using CORNET

The construction of co-expression networks for multiple input genes was conducted using the CORNET tool [31]. The co-expression tool calculates the correlation between gene expression profiles using one or more precompiled expression datasets and, as such, identifies possible functional associations between genes. Out of all of the available expression data, we selected the subgroup consisting of 69 ATH1 AtGenExpress biotic stress compendium expression data. All the expression data were processed using RMA from the R BioConductor package and making use of the CDF described in Casneuf et al. [79]. Pearson’s correlation coefficients were calculated between the given genes. Correlation coefficients higher and lower than a certain value are reported. Pearson’s correlation coefficient (PCC) was used at a cut off ≥0.7. Networks and associated evidence were visualised in Cytoscape 2.7.0.

Gene Networks, microRNAs and Connections to Post-transcriptional Gene Regulation

The Gene networks tool in VirtualPlant groups individual genes into a supernode based on shared functional properties, such as GO terms, KEGG pathways, gene families and even similar annotations. Edges were drawn between two supernodes when at least one gene or gene product in each supernode showed a molecular interaction. To improve the regulatory interaction predictions, we filtered the transcription factor:target gene predictions to include only the transcription factor and target pairs whose expression values were correlated in the microarray experiment [80]. The selected statistic for the calculation of correlations in this analysis was the Pearson’s correlation coefficient, with a cut-off value of less than or equal to 0.7. The results were then cross-compared with all of the microRNA genes curated in the microRNA Registry (microrna.sanger.ac.uk/sequences) [81] and in the Arabidopsis Small RNA Project (ASRP) Database [82]. In certain cases, we also compared the results with microRNAs and precursor candidates predicted for the A. thaliana genome by the algorithm findMicroRNA [83]. We followed specific criteria required for the annotation of plant microRNAs, including experimental and computational data as well as refinements of standard nomenclature, as described in [84] [85].

Supporting Information

File S1.

GO-annotation cytoscape network file (.cys) for common genes.

https://doi.org/10.1371/journal.pone.0058987.s001

(CYS)

File S2.

GO-annotation cytoscape network file (.cys) for aphid-specific genes.

https://doi.org/10.1371/journal.pone.0058987.s002

(CYS)

File S3.

GO-annotation cytoscape network file (.cys) for P. syringae-specific genes.

https://doi.org/10.1371/journal.pone.0058987.s003

(CYS)

File S4.

MapMan input file for all aphid-responsive genes and corresponding log2-transformed expression values.

https://doi.org/10.1371/journal.pone.0058987.s004

(XLS)

File S5.

MapMan input file for all P. syringae-responsive genes and corresponding log2-transformed expression values.

https://doi.org/10.1371/journal.pone.0058987.s005

(XLS)

File S6.

Lists of aphid-specific and P. syringae-specific genes and their log2-transformed expression values related to biotic stress signalling processes.

https://doi.org/10.1371/journal.pone.0058987.s006

(XLS)

File S7.

List of the 31 heat shock protein (HSP) genes that were differentially expressed only during aphid treatment and their log2-transformed expression values.

https://doi.org/10.1371/journal.pone.0058987.s007

(XLS)

File S8.

List of the proteolyitc enzymes differentially expressed during the aphid and P. syringae treatments.

https://doi.org/10.1371/journal.pone.0058987.s008

(XLS)

File S9.

List of genes related to secondary metabolic processes that were differentially regulated during the aphid and P. syringae treatments.

https://doi.org/10.1371/journal.pone.0058987.s009

(XLS)

File S10.

List of genes involved in cell wall precursor synthesis that were differentially regulated during the aphid and P. syringae treatments.

https://doi.org/10.1371/journal.pone.0058987.s010

(XLS)

File S11.

A comparative list of the differentially expressed genes (aphid-specific and P. syringae-specific genes) involved in hormonal pathways with their corresponding log2-transformed expression values.

https://doi.org/10.1371/journal.pone.0058987.s011

(XLS)

File S12.

Details of all of the TFs that were differentially expressed during the aphid and Pseudmonas treatments.

https://doi.org/10.1371/journal.pone.0058987.s012

(XLS)

Acknowledgments

The authors thank Torfinn Sparstad for excellent technical assistance.

Author Contributions

Coordinated the study: AMB. Conceived and designed the experiments: PB AMB. Performed the experiments: PB AK DHT. Analyzed the data: PB PW. Contributed reagents/materials/analysis tools: PW AMB. Wrote the paper: PB PW AK AMB.

References

  1. 1. Ahuja I, de Vos RC, Bones AM, Hall RD (2010) Plant molecular stress responses face climate change. Trends Plant Sci 15: 664–674.
  2. 2. Koornneef M, Meinke D (2010) The development of Arabidopsis as a model plant. Plant Journal 61: 909–921.
  3. 3. Pitzschke A, Hirt H (2010) Bioinformatic and systems biology tools to generate testable models of signaling pathways and their targets. Plant Physiol 152: 460–469.
  4. 4. Reymond P, Bodenhausen N, Van Poecke RM, Krishnamurthy V, Dicke M, et al. (2004) A conserved transcript pattern in response to a specialist and a generalist herbivore. Plant Cell 16: 3132–3147.
  5. 5. Kreps JA, Wu Y, Chang HS, Zhu T, Wang X, et al. (2002) Transcriptome changes for Arabidopsis in response to salt, osmotic, and cold stress. Plant Physiol 130: 2129–2141.
  6. 6. Dicke M, van Loon JJ, Soler R (2009) Chemical complexity of volatiles from plants induced by multiple attack. Nat Chem Biol 5: 317–324.
  7. 7. Knight H, Knight MR (2001) Abiotic stress signalling pathways: specificity and cross-talk. Trends Plant Sci 6: 262–267.
  8. 8. Van Emden HF, Harrington R (2007) Aphids as crop pests/edited by Helmut F. van Emden and Richard Harrington. Wallingford, UK; Cambridge, MA: CABI. xxviii, 717 p., 716 p. of plates p.
  9. 9. Goggin FL (2007) Plant-aphid interactions: molecular and ecological perspectives. Curr Opin Plant Biol 10: 399–408.
  10. 10. Broekgaarden C, Poelman EH, Steenhuis G, Voorrips RE, Dicke M, et al. (2008) Responses of Brassica oleracea cultivars to infestation by the aphid Brevicoryne brassicae: an ecological and molecular approach. Plant Cell Environ 31: 1592–1605.
  11. 11. Spiller NJ, Kimmins FM, Llewellyn M (1985) Fine-Structure of Aphid Stylet Pathways and Its Use in Host Plant-Resistance Studies. Entomologia Experimentalis Et Applicata 38: 293–295.
  12. 12. Kusnierczyk A, Winge P, Midelfart H, Armbruster WS, Rossiter JT, et al. (2007) Transcriptional responses of Arabidopsis thaliana ecotypes with different glucosinolate profiles after attack by polyphagous Myzus persicae and oligophagous Brevicoryne brassicae. J Exp Bot 58: 2537–2552.
  13. 13. Whalen MC, Innes RW, Bent AF, Staskawicz BJ (1991) Identification of Pseudomonas syringae pathogens of Arabidopsis and a bacterial locus determining avirulence on both Arabidopsis and soybean. Plant Cell 3: 49–59.
  14. 14. Fumiaki Katagiri RT, Sheng Yang He (2002) The Arabidopsis Thaliana-Pseudomonas Syringae Interaction. The Arabidopsis Book: American Society of Plant Biologists. pp. e012.
  15. 15. Preston GM (2000) Pseudomonas syringae pv. tomato: the right pathogen, of the right plant, at the right time. Mol Plant Pathol 1: 263–275.
  16. 16. Weiler EW, Kutchan TM, Gorba T, Brodschelm W, Niesel U, et al. (1994) The Pseudomonas Phytotoxin Coronatine Mimics Octadecanoid Signaling Molecules of Higher-Plants (Vol 345, Pg 9, 1994). Febs Letters 349: 317–317.
  17. 17. Brooks DM, Bender CL, Kunkel BN (2005) The Pseudomonas syringae phytotoxin coronatine promotes virulence by overcoming salicylic acid-dependent defences in Arabidopsis thaliana. Mol Plant Pathol 6: 629–639.
  18. 18. Melotto M, Underwood W, Koczan J, Nomura K, He SY (2006) Plant stomata function in innate immunity against bacterial invasion. Cell 126: 969–980.
  19. 19. Moreau Y, Aerts S, De Moor B, De Strooper B, Dabrowski M (2003) Comparison and meta-analysis of microarray data: from the bench to the computer desk. Trends Genet 19: 570–577.
  20. 20. Cramer GR, Urano K, Delrot S, Pezzotti M, Shinozaki K (2011) Effects of abiotic stress on plants: a systems biology perspective. BMC Plant Biol 11: 163.
  21. 21. Konika Chawla PB, Kuiper M, Bones AM (2010) Systems Biology: a promising tool to study abiotic stress responses. Omics and Plant Abiotic Stress Tolerance. USA Bentham Publishers.
  22. 22. Williams EJ, Bowles DJ (2004) Coexpression of neighboring genes in the genome of Arabidopsis thaliana. Genome Res 14: 1060–1067.
  23. 23. Weston DJ, Gunter LE, Rogers A, Wullschleger SD (2008) Connecting genes, coexpression modules, and molecular signatures to environmental stress phenotypes in plants. BMC Syst Biol 2: 16.
  24. 24. Bartel DP (2004) MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116: 281–297.
  25. 25. Zhou X, Wang G, Sutoh K, Zhu JK, Zhang W (2008) Identification of cold-inducible microRNAs in plants by transcriptome analysis. Biochim Biophys Acta 1779: 780–788.
  26. 26. Davuluri RV, Sun H, Palaniswamy SK, Matthews N, Molina C, et al. (2003) AGRIS: Arabidopsis gene regulatory information server, an information resource of Arabidopsis cis-regulatory elements and transcription factors. BMC Bioinformatics 4: 25.
  27. 27. Bulow L, Bolivar JC, Ruhe J, Brill Y, Hehl R (2012) ‘MicroRNA Targets’, a new AthaMap web-tool for genome-wide identification of miRNA targets in Arabidopsis thaliana. BioData Min 5: 7.
  28. 28. Bulow L, Engelmann S, Schindler M, Hehl R (2009) AthaMap, integrating transcriptional and post-transcriptional data. Nucleic Acids Res 37: D983–986.
  29. 29. Higo K, Ugawa Y, Iwamoto M, Korenaga T (1999) Plant cis-acting regulatory DNA elements (PLACE) database: 1999. Nucleic Acids Res 27: 297–300.
  30. 30. Katari MS, Nowicki SD, Aceituno FF, Nero D, Kelfer J, et al. (2010) VirtualPlant: a software platform to support systems biology research. Plant Physiol 152: 500–515.
  31. 31. De Bodt S, Carvajal D, Hollunder J, Van den Cruyce J, Movahedi S, et al. (2010) CORNET: a user-friendly tool for data mining and integration. Plant Physiol 152: 1167–1179.
  32. 32. Bethke G, Scheel D, Lee J (2009) Sometimes new results raise new questions: the question marks between mitogen-activated protein kinase and ethylene signaling. Plant Signal Behav 4: 672–674.
  33. 33. Lin W, Katagiri F, Glazebrook J (2009) Arabidopsis defense response against Pseudomonas syringae - Effects of major regulatory genes and the impact of coronatine; 2009 17–21 May. 1–4.
  34. 34. Thilmony R, Underwood W, He SY (2006) Genome-wide transcriptional analysis of the Arabidopsis thaliana interaction with the plant pathogen Pseudomonas syringae pv. tomato DC3000 and the human pathogen Escherichia coli O157 : H7. Plant Journal 46: 34–53.
  35. 35. Ahuja I, Kissen R, Bones AM (2012) Phytoalexins in defense against pathogens. Trends Plant Sci 17: 73–90.
  36. 36. Burow M, Zhang ZY, Ober JA, Lambrix VM, Wittstock U, et al. (2008) ESP and ESM1 mediate indol-3-acetonitrile production from indol-3-ylmethyl glucosinolate in Arabidopsis. Phytochemistry 69: 663–671.
  37. 37. Pajerowska-Mukhtar K, Dong X (2009) A kiss of death–proteasome-mediated membrane fusion and programmed cell death in plant defense against bacterial infection. Genes Dev 23: 2449–2454.
  38. 38. Theis N, Lerdau M (2003) The evolution of function in plant secondary metabolites. International Journal of Plant Sciences 164: S93–S102.
  39. 39. Strauss SY, Rudgers JA, Lau JA, Irwin RE (2002) Direct and ecological costs of resistance to herbivory. Trends in Ecology & Evolution 17: 278–285.
  40. 40. Foyer CH, Noctor G, van Emden HF (2007) An evaluation of the costs of making specific secondary metabolites: Does the yield penalty incurred by host plant resistance to insects result from competition for resources? International Journal of Pest Management 53: 175–182.
  41. 41. Sasidharan R, Chinnappa CC, Staal M, Elzenga JT, Yokoyama R, et al. (2010) Light quality-mediated petiole elongation in Arabidopsis during shade avoidance involves cell wall modification by xyloglucan endotransglucosylase/hydrolases. Plant Physiol 154: 978–990.
  42. 42. Rose JK, Braam J, Fry SC, Nishitani K (2002) The XTH family of enzymes involved in xyloglucan endotransglucosylation and endohydrolysis: current perspectives and a new unifying nomenclature. Plant Cell Physiol 43: 1421–1435.
  43. 43. De Vos M, Van Oosten VR, Van Poecke RM, Van Pelt JA, Pozo MJ, et al. (2005) Signal signature and transcriptome changes of Arabidopsis during pathogen and insect attack. Mol Plant Microbe Interact 18: 923–937.
  44. 44. Browse J, Howe GA (2008) New weapons and a rapid response against insect attack. Plant Physiol 146: 832–838.
  45. 45. Kusnierczyk A, Tran DH, Winge P, Jorstad TS, Reese JC, et al. (2011) Testing the importance of jasmonate signalling in induction of plant defences upon cabbage aphid (Brevicoryne brassicae) attack. BMC Genomics 12: 423.
  46. 46. Farmer EE, Ryan CA (1992) Octadecanoid Precursors of Jasmonic Acid Activate the Synthesis of Wound-Inducible Proteinase Inhibitors. Plant Cell 4: 129–134.
  47. 47. Ulland S, Ian E, Mozuraitis R, Borg-Karlson AK, Meadow R, et al. (2008) Methyl salicylate, identified as primary odorant of a specific receptor neuron type, inhibits oviposition by the moth Mamestra brassicae L. (Lepidoptera, Noctuidae). Chemical Senses 33: 35–46.
  48. 48. Luna E, Bruce TJ, Roberts MR, Flors V, Ton J (2011) Next Generation Systemic Acquired Resistance. Plant Physiol.
  49. 49. Grant M, Truman B, Zabala MD, Bennett M, Turnbull C (2006) Towards a systems approach to plant defense responses. Comparative Biochemistry and Physiology a-Molecular & Integrative Physiology 143: S137–S137.
  50. 50. Kissen R, Rossiter JT, Bones AM (2009) The ‘mustard oil bomb’: not so easy to assemble?! Localization, expression and distribution of the components of the myrosinase enzyme system. Phytochemistry Reviews 8: 69–86.
  51. 51. Kissen R, Bones AM (2009) Nitrile-specifier proteins involved in glucosinolate hydrolysis in Arabidopsis thaliana. J Biol Chem 284: 12057–12070.
  52. 52. Eulgem T (2005) Regulation of the Arabidopsis defense transcriptome. Trends Plant Sci 10: 71–78.
  53. 53. Reddy AS, Ali GS, Celesnik H, Day IS (2011) Coping with stresses: roles of calcium- and calcium/calmodulin-regulated gene expression. Plant Cell 23: 2010–2032.
  54. 54. Gigolashvili T, Berger B, Mock HP, Muller C, Weisshaar B, et al. (2007) The transcription factor HIG1/MYB51 regulates indolic glucosinolate biosynthesis in Arabidopsis thaliana. Plant J 50: 886–901.
  55. 55. Munemasa S, Hossain MA, Nakamura Y, Mori IC, Murata Y (2011) The Arabidopsis calcium-dependent protein kinase, CPK6, functions as a positive regulator of methyl jasmonate signaling in guard cells. Plant Physiol 155: 553–561.
  56. 56. Wright AJ, Knight H, Knight MR (2002) Mechanically stimulated TCH3 gene expression in Arabidopsis involves protein phosphorylation and EIN6 downstream of calcium. Plant Physiol 128: 1402–1409.
  57. 57. Fujita M, Fujita Y, Maruyama K, Seki M, Hiratsu K, et al. (2004) A dehydration-induced NAC protein, RD26, is involved in a novel ABA-dependent stress-signaling pathway. Plant J 39: 863–876.
  58. 58. Tamura K, Fukao Y, Iwamoto M, Haraguchi T, Hara-Nishimura I (2010) Identification and characterization of nuclear pore complex components in Arabidopsis thaliana. Plant Cell 22: 4084–4097.
  59. 59. Doczi R, Brader G, Pettko-Szandtner A, Rajh I, Djamei A, et al. (2007) The Arabidopsis mitogen-activated protein kinase kinase MKK3 is upstream of group C mitogen-activated protein kinases and participates in pathogen signaling. Plant Cell 19: 3266–3279.
  60. 60. Tarasenko VI, Garnik EY, Shmakov VN, Konstantinov YM (2009) Induction of Arabidopsis gdh2 gene expression during changes in redox state of the mitochondrial respiratory chain. Biochemistry (Mosc) 74: 47–53.
  61. 61. Prandl R, Hinderhofer K, Eggers-Schumacher G, Schoffl F (1998) HSF3, a new heat shock factor from Arabidopsis thaliana, derepresses the heat shock response and confers thermotolerance when overexpressed in transgenic plants. Mol Gen Genet 258: 269–278.
  62. 62. Bricchi I, Bertea CM, Occhipinti A, Paponov IA, Maffei ME (2012) Dynamics of membrane potential variation and gene expression induced by Spodoptera littoralis, Myzus persicae, and Pseudomonas syringae in Arabidopsis. Plos One 7: e46673.
  63. 63. Sato M, Tsuda K, Wang L, Coller J, Watanabe Y, et al. (2010) Network modeling reveals prevalent negative regulatory relationships between signaling sectors in Arabidopsis immune signaling. PLoS Pathog 6: e1001011.
  64. 64. Bhardwaj V, Meier S, Petersen LN, Ingle RA, Roden LC (2011) Defence responses of Arabidopsis thaliana to infection by Pseudomonas syringae are regulated by the circadian clock. PLoS One 6: e26968.
  65. 65. Thilmony R, Underwood W, He SY (2006) Genome-wide transcriptional analysis of the Arabidopsis thaliana interaction with the plant pathogen Pseudomonas syringae pv. tomato DC3000 and the human pathogen Escherichia coli O157: H7. Plant J 46: 34–53.
  66. 66. Kusnierczyk A, Winge P, Jorstad TS, Troczynska J, Rossiter JT, et al. (2008) Towards global understanding of plant defence against aphids–timing and dynamics of early Arabidopsis defence responses to cabbage aphid (Brevicoryne brassicae) attack. Plant Cell Environ 31: 1097–1115.
  67. 67. Hornik K (2011) The R FAQ.
  68. 68. Jorstad TS, Midelfart H, Bones AM (2008) A mixture model approach to sample size estimation in two-sample comparative microarray experiments. BMC Bioinformatics 9: 117.
  69. 69. Jorstad TS, Langaas M, Bones AM (2007) Understanding sample size: what determines the required number of microarrays for an experiment? Trends Plant Sci 12: 46–50.
  70. 70. Smyth GK (2005) Limma: linear models for microarray data. In: Gentleman VC, Dudoit S, Irizarry R, Huber W, editors. Bioinformatics and Computational Biology Solutions using R and Bioconductor. New York: Springer.
  71. 71. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, et al. (2000) Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25: 25–29.
  72. 72. Lamesch P, Berardini TZ, Li D, Swarbreck D, Wilks C, et al. (2012) The Arabidopsis Information Resource (TAIR): improved gene annotation and new tools. Nucleic Acids Res 40: D1202–1210.
  73. 73. Bindea G, Mlecnik B, Hackl H, Charoentong P, Tosolini M, et al. (2009) ClueGO: a Cytoscape plug-in to decipher functionally grouped gene ontology and pathway annotation networks. Bioinformatics 25: 1091–1093.
  74. 74. Maere S, Heymans K, Kuiper M (2005) BiNGO: a Cytoscape plugin to assess overrepresentation of gene ontology categories in biological networks. Bioinformatics 21: 3448–3449.
  75. 75. Tardieu F (2003) Virtual plants: modelling as a tool for the genomics of tolerance to water deficit. Trends Plant Sci 8: 9–14.
  76. 76. Guo A, He K, Liu D, Bai S, Gu X, et al. (2005) DATF: a database of Arabidopsis transcription factors. Bioinformatics 21: 2568–2569.
  77. 77. Rivals I, Personnaz L, Taing L, Potier MC (2007) Enrichment or depletion of a GO category within a class of genes: which test? Bioinformatics 23: 401–407.
  78. 78. Ge Y, Dudoit S., Speed T.P. (2003) Resampling-based multiple testing for microarray data analysis. UC Berkeley Dept. of Statistics. 633 p.
  79. 79. Casneuf T, Van de Peer Y, Huber W (2007) In situ analysis of cross-hybridisation on microarrays and the inference of expression correlation. BMC Bioinformatics 8: 461.
  80. 80. Gutierrez RA, Lejay LV, Dean A, Chiaromonte F, Shasha DE, et al. (2007) Qualitative network models and genome-wide expression data define carbon/nitrogen-responsive molecular machines in Arabidopsis. Genome Biol 8: R7.
  81. 81. Griffiths-Jones S (2004) The microRNA Registry. Nucleic Acids Res 32: D109–111.
  82. 82. Backman TW, Sullivan CM, Cumbie JS, Miller ZA, Chapman EJ, et al. (2008) Update of ASRP: the Arabidopsis Small RNA Project database. Nucleic Acids Res 36: D982–985.
  83. 83. Adai A, Johnson C, Mlotshwa S, Archer-Evans S, Manocha V, et al. (2005) Computational prediction of miRNAs in Arabidopsis thaliana. Genome Res 15: 78–91.
  84. 84. Ambros V, Bartel B, Bartel DP, Burge CB, Carrington JC, et al. (2003) A uniform system for microRNA annotation. RNA 9: 277–279.
  85. 85. Meyers BC, Axtell MJ, Bartel B, Bartel DP, Baulcombe D, et al. (2008) Criteria for annotation of plant MicroRNAs. Plant Cell 20: 3186–3190.
  86. 86. Hsieh LC, Lin SI, Shih AC, Chen JW, Lin WY, et al. (2009) Uncovering small RNA-mediated responses to phosphate deficiency in Arabidopsis by deep sequencing. Plant Physiol 151: 2120–2132.
  87. 87. Allen E, Xie Z, Gustafson AM, Carrington JC (2005) microRNA-directed phasing during trans-acting siRNA biogenesis in plants. Cell 121: 207–221.
  88. 88. Jones-Rhoades MW, Bartel DP (2004) Computational identification of plant microRNAs and their targets, including a stress-induced miRNA. Mol Cell 14: 787–799.
  89. 89. Kurihara Y, Watanabe Y (2004) Arabidopsis micro-RNA biogenesis through Dicer-like 1 protein functions. Proc Natl Acad Sci U S A 101: 12753–12758.
  90. 90. Sieber P, Wellmer F, Gheyselinck J, Riechmann JL, Meyerowitz EM (2007) Redundancy and specialization among plant microRNAs: role of the MIR164 family in developmental robustness. Development 134: 1051–1060.
  91. 91. Nikovics K, Blein T, Peaucelle A, Ishida T, Morin H, et al. (2006) The balance between the MIR164A and CUC2 genes controls leaf margin serration in Arabidopsis. Plant Cell 18: 2929–2945.
  92. 92. Merchan F, Boualem A, Crespi M, Frugier F (2009) Plant polycistronic precursors containing non-homologous microRNAs target transcripts encoding functionally related proteins. Genome Biol 10: R136.
  93. 93. Rhoades MW, Reinhart BJ, Lim LP, Burge CB, Bartel B, et al. (2002) Prediction of plant microRNA targets. Cell 110: 513–520.
  94. 94. Park W, Li J, Song R, Messing J, Chen X (2002) CARPEL FACTORY, a Dicer homolog, and HEN1, a novel protein, act in microRNA metabolism in Arabidopsis thaliana. Curr Biol 12: 1484–1495.
  95. 95. Palatnik JF, Allen E, Wu X, Schommer C, Schwab R, et al. (2003) Control of leaf morphogenesis by microRNAs. Nature 425: 257–263.
  96. 96. Kurihara Y, Watanabe Y (2004) Arabidopsis micro-RNA biogenesis through Dicer-like 1 protein functions. Proceedings of the National Academy of Sciences of the United States of America 101: 12753–12758.
  97. 97. Jones-Rhoades MW, Bartel DP (2004) Computational identification of plant micro-RNAs and their targets, including a stress-induced miRNA. Mol Cell 14: 787–799.
  98. 98. Nag A, King S, Jack T (2009) miR319a targeting of TCP4 is critical for petal growth and development in Arabidopsis. Proc Natl Acad Sci U S A 106: 22534–22539.
  99. 99. Sunkar R, Zhu JK (2004) Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell 16: 2001–2019.
  100. 100. Cuperus JT, Carbonell A, Fahlgren N, Garcia-Ruiz H, Burke RT, et al. (2010) Unique functionality of 22-nt miRNAs in triggering RDR6-dependent siRNA biogenesis from target transcripts in Arabidopsis. Nat Struct Mol Biol 17: 997–1003.
  101. 101. Jagadeeswaran G, Saini A, Sunkar R (2009) Biotic and abiotic stress down-regulate miR398 expression in Arabidopsis. Planta 229: 1009–1014.