Transcriptional responses of Arabidopsis thaliana to chewing and sucking insect herbivores

We tested the hypothesis that Arabidopsis can recognize and respond differentially to insect species at the transcriptional level using a genome wide microarray. Transcriptional reprogramming was characterized using co-expression analysis in damaged and undamaged leaves at two times in response to mechanical wounding and four insect species. In all, 2778 (10.6%) of annotated genes on the array were differentially expressed in at least one treatment. Responses differed mainly between aphid and caterpillar and sampling times. Responses to aphids and caterpillars shared only 10% of up-regulated and 8% of down-regulated genes. Responses to two caterpillars shared 21 and 12% of up- and down-regulated genes, whereas responses to the two aphids shared only 7 and 4% of up-regulated and down-regulated genes. Overlap in genes expressed between 6 and 24 h was 3–15%, and depended on the insect species. Responses in attacked and unattacked leaves differed at 6 h but converged by 24 h. Genes responding to the insects are also responsive to many stressors and included primary metabolism. Aphids down-regulated amino acid catabolism; caterpillars stimulated production of amino acids involved in glucosinolate synthesis. Co-expression analysis revealed 17 response networks. Transcription factors were a major portion of differentially expressed genes throughout and responsive genes shared most of the known or postulated binding sites. However, cis-element composition of genes down regulated by the aphid M. persicae was unique, as were those of genes down-regulated by caterpillars. As many as 20 cis-elements were over-represented in one or more treatments, including some from well-characterized classes and others as yet uncharacterized. We suggest that transcriptional changes elicited by wounding and insects are heavily influenced by transcription factors and involve both enrichment of a common set of cis-elements and a unique enrichment of a few cis-elements in responding genes.


INTRODUCTION
Approximately a quarter of all described eukaryotic species are insects that feed on plants, and as a group they are thought to exert strong selection on plants to detect and repel them (Ehrlich and Raven, 1964;Futuyma and Agrawal, 2009). Any individual plant may be attacked by many species of herbivores feeding in a variety of ways. Some species feed by chewing tissues, which involves wounding and water loss, while others feed by inserting needle-like mouthparts (stylets) between and into plant cells, causing less overall damage. Some vector plant diseases whereas others inadvertently introduce microbial pathogens from leaf surfaces or gut contents. Although field studies of plant evolution in response to herbivory are rare, they indicate substantial selection for plant resistance traits Zust et al., 2012).
Many potential plant resistance traits are modified or enhanced in response to insect attack and can comprise a significant barrier to insect feeding and concomitant pathogen introduction (Karban and Baldwin, 1997;Dicke and Hilker, 2003). Since insects wound as they feed, there is always a question of whether these changes are general responses to wounding or are specific to the insect species. Although there are currently few examples of gene-for-gene recognition systems in plant-herbivore interactions analogous to pathogen effector triggered immunity (Rossi et al., 1998;Aggarwal et al., 2014), specificity in plant phenotypic responses to different herbivores is commonplace (Ali and Agrawal, 2012;Barrett and Heil, 2012). Not surprisingly, responses to sucking insects usually differ from those elicited by chewing insects. More surprising is that biochemical responses can be herbivore species-specific, such that some insects suppress or fail to elicit defense responses (Alba et al., 2012). Although the basis of species-specific responses is largely unknown, the different responses elicited by insects and mechanical wounding support the view that elicitors in saliva or regurgitant cue them (Bonaventure et al., 2011;Maffei et al., 2012). Based on the pathogen recognition system of plants, a model of insect recognition by plants has been proposed with herbivore-associated molecular patterns (HAMPs) and damage-associated molecular patterns (DAMPs) cueing downstream HAMPs-triggered immunity (HTI) (Heil, 2009).
Since transcriptional reprogramming underlies defense responses, HAMPs should reflect underlying transcriptional changes. There are now many published studies of transcriptional change in plants responding to insect herbivores (Heidel-Fischer et al., 2014). Drawing conclusions among them is confounded by the many differences in experimental approaches, including herbivore species and treatments, plant species and tissues, sampling times, and gene expression platforms. Replicated, whole-genome expression profiling is important to comparative studies because it can reveal the full range of transcriptional responses, unlike partial genome arrays or targeted quantitative reverse-transcription polymerase chain reaction (qRT-PCR). The age of plant tissue treated and sampled should be consistent to address the large effect that age can have on gene expression profiles. Use of multiple herbivores and sampling times is also important for comparative studies because plant responses to different herbivores can differ qualitatively and quantitatively over time.
To identify HAMPs that address these experimental issues, we report here a fully replicated study of genome-wide transcriptional responses by Arabidopsis thaliana (L.) to wounding, to chewing insects (caterpillars) and to stylet feeders (aphids). We used a genome-wide Arabidopsis microarray to examine all transcriptional changes at two time intervals after mechanical wounding or being attacked by four insect species. The insect species were two aphids (phloem feeders) and two caterpillars (leaf chewers). Both feeding types included one with an extremely broad diet (generalist) and another with a narrow diet (specialist) focused on Brassicaceae, the plant family to which A. thaliana belongs. This design enabled testing the hypothesis that Arabidopsis can recognize and respond individually to these insects at the transcriptional level. Common patterns in fundamental reprogramming of plant metabolism were identified using gene enrichment tests. Candidate genes that may determine species-specific responses were identified with coexpression network analysis. Frequencies determined for coregulated gene promoters linked to genes differentially expressed suggest that reprogramming of expression patterns elicited by wounding and different insect species involve both a common set and a unique, small set of cis-elements.

INSECTS AND PLANTS
We assessed the transcriptional responses in rosette leaves of 4-week-old A. thaliana ecotype Columbia Wild-Type (Col WT) to attack by larvae of two leaf chewing caterpillar species and by adults and nymphs of two phloem feeding aphid species ( Table 1). The caterpillars were Spodoptera exigua (Hübner) (Noctuidae), which feeds on as many as 20 plant families including Brassicaceae (Greenberg et al., 2001), and Pieris rapae (L.) (Pieridae), which feeds exclusively on plants in the family Brassicaceae, to which Arabidopsis belongs (Renwick and Lopez, 1999). The aphids (both in Aphididae) were Myzus persicae (Sulzer), a broad generalist feeding on species in many plant families including Brassicaceae, and Brevicoryne brassicae (L.), whose feeding is limited to the Brassicaceae (Blackman and Eastop, 1994). We henceforth refer to S. exigua and M. persicae as "dietary generalists" and P. rapae and B. brassicae as "dietary specialists." Both aphids were maintained as plant virus free clones on pakchoi plants (Brassica campestris L. ssp. chinensis cv. Black Behi). Eggs of S. exigua were obtained from Benzon Research and larvae were reared on artificial diet (Bio-Serv). Pieris rapae was maintained as a culture in our lab on pak-choi and originated from the Carolina Biological Supply Company. Both caterpillar species were transferred to Col WT plants one day before the experiments for acclimation to the new host. Plant seeds were vernalized in 2% agar and sown into 6 × 5 cm pots containing sterile Metromix 200 soil (Sun Gro Horticulture). Plants were chamber grown at 22 ± 1 • C, 65 ± 5% relative humidity, and 200 μmol

EXPERIMENTAL PLANT TREATMENTS
The treatments are summarized in Table 1. The caterpillar treatment was designed to capture early gene expression events and minimize variation due to leaf age and amount of insect damage. All leaves selected for treatment and harvest were fully expanded mature leaves. Second and third instar (N = 6-10) caterpillars of both species were allowed to feed for 2-4 h to generate 6 leaves of similar age per plant with ∼20% leaf area removed. Caterpillars were wrangled as needed with a size 0 camel's hair brush to concentrate their feeding on the 6 leaves sampled. Caterpillars were removed when sufficient damage was achieved and the plants were returned to the growth chamber. The mechanical wounding treatment was designed to approximate insect damage to tissues by running a sterile pattern wheel across both sides of the midrib of 6 leaves of similar age on each plant, once at the beginning of the caterpillar treatment and again half way through the caterpillar wrangling period. Control plants were jiggled with the brush to simulate leaf movement caused by caterpillar wrangling or mechanical wounding. Damaged leaves ("local" to the attack) were harvested for gene expression 6 and 24 h after the start of caterpillar damage or wounding. Unwounded leaves ("systemic") were harvested separately from size-matched damaged or wounded leaves. Leaves from 3 to 4 plants per treatment (caterpillar, wounding, and controls) were pooled for each of the four biological replicates. Treatment of plants with aphids was different from that for caterpillars because aphids have effects on plants that are much weaker and slower to develop than those of caterpillars (Mewis et al., 2005(Mewis et al., , 2006 and aphids cannot be readily contained on individual leaves. Sub-adult (final instar) and adult aphids (N = 20) were placed on plants whose rosettes were caged at the soil line by transparent mylar cylinders (5 cm diameter, 9 cm high) with tops of fine mesh gauze (<0.01 mm mesh wide) to maintain air exchange. Controls were caged plants without aphids and all plants were returned to the growth chamber. After 1 week, all cages and aphids were removed and control plants were jiggled with a camel's hair brush to simulate leaf movement caused by aphid removal. Plants were returned to the growth chamber and whole plants were harvested per treatment (aphid and control) for gene expression 6 and 24 h after aphid removal.

RNA ISOLATION
Total RNA was isolated from leaves using a modified TRIZOL extraction method as follows. Approximately 0.5 g of plant material was ground to a powder in liquid nitrogen using mortar and pestle, resuspended without thawing in 6 ml TRIZOL reagent (Invitrogen) by vortexing, and incubated at 65 • C for 5 min with regular mixing. Cell debris was pelleted by centrifugation (30 min, 12,000 g, 4 • C) and the supernatant was extracted twice with 3 ml chloroform with the aqueous phase recovered each time after centrifugation (20 min, 12,000 × g, 4 • C). RNA was precipitated from this phase at room temperature for 5 min with 0.5 volumes each of 0.8 M sodium citrate and isopropanol. The RNA pellet obtained after centrifugation (30 min, 12,000 × g, 4 • C) was washed with 70% ethanol, recovered again by centrifugation, air dried for 5 min and resuspended in 200 μl nuclease free water. The RNA was further purified by standard ethanolic sodium acetate precipitation at −20 • C overnight. Following a wash with 70% ethanol the pelleted RNA (30 min, 12,000 g, 4 • C) was air dried and resuspended in nuclease free water to an approximate concentration of 5 μg/μl whose actual concentration was determined spectrophotometrically. Quality of randomly selected samples of RNA was determined using a 2100 Bioanalyzer (Agilent Technologies).

MICROARRAYS AND THE PREPARATION AND HYBRIDIZATION OF cDNA TO ARRAYS
Design and production of microarrays with 26,090 Arabidopsis oligonucleotide targets (Qiagen Operon), 12 housekeeping gene oligos (Qiagen Operon) as positive controls and 16 oligos with no similarity to any Arabidopsis gene (4 synthesized human genes and 12 others) as negative controls and internal, spike controls was previously described (Ehlting et al., 2005). Location of all oligos and orientation marker was also previously described (Ehlting et al., 2005) and provided in the platform file deposited to the NCBI GEO database (http://www.ncbi.nlm.nih.gov/geo/; series accession GSE62287).
Total RNA isolated as previously described (Ehlting et al., 2005) was used for a direct labeling procedure that generated the cDNA hybridized microarray slides. All procedures for labeling, microarray hybridization, image scanning, and the identification and quantification of spots were performed as previously described (Ehlting et al., 2005(Ehlting et al., , 2008. Microarray experiments involved hybridizing three or four replicate arrays per treatment. The RNA from control plants harvested with each treatment was pooled within time points to obtain sufficient control RNA. Labeled cDNA derived from this control RNA was co-hybridized with labeled cDNA derived from RNA isolated from independent biological replicates receiving insects or mechanical wounding, resulting in four biological replicates per treatment. Dye bias was accounted for by swapping dye labeling between treatment and control samples among bioreplicates within each treatment group. Hybridization order was randomized to avoid biases due to the hybridization time.

MICROARRAY DATA ANALYSIS
Analyses of gene specific elements used customized scripts for R and Bioconductor (Team, 2014). Background correction performed as previously described (Ehlting et al., 2008) excluded on average 19% (12-27%) of all spots from further analyses as non-detectable. Background corrected signal intensities were used for subsequent Loess normalization and statistical testing (Student's t-test, analysis of variance, false discovery rate of expression ratios) between treatment and corresponding control, as previously described (Ehlting et al., 2008). The expression data was initially filtered to obtain genes with a t-test p < 0.05 and fold-change >3 between treatment and control in at least one sample. Normalized expression ratios for these genes from all treatments were then used to perform an analysis of variance (ANOVA) and an estimate of the false discovery rate (FDR) based on the distribution of parametric p-values (Supplemental Figure 1). Means of the normalized expression ratios were subjected to a hierarchical clustering analysis with average linkage using Genesis v1.2 (Institute for Biomedical Engineering).
A high estimated FDR arose from our experimental design intended to capture the biological variation inherent in plant herbivore interactions (biological replicates treated with a small number of herbivores that behaved individually). Lower p-values from t-tests were associated with a lower false discovery rate (Supplemental Table 1). However, a relatively large number of array probes were associated with high p-values which still contain a substantial number of truly differentially expressed genes as estimated from the higher frequency of genes in these p-value bins compared to the frequency expected if no genes were differentially expressed (indicated by the horizontal line in Supplemental Figure 1). Although a p-value cut-off of 0.01 would reduce the number of falsely discovered genes, a substantial number of truly differentially expressed genes would also be missed. Therefore, we assumed that high fold change difference is associated with a lower likelihood of being a false positive (Pylatuik and Fobert, 2005) to initially obtain 3123 genes as "differentially expressed" (i.e., genes with treatment-induced change in transcript abundance) as those genes for each time point that were associated with a t-test p < 0.05 (accepting a false discovery rate of up to 0.3) and also displayed a more than two-fold change between treatment and control (Supplemental Table 1). After removing genes identified merely as chromosome loci, BAC clones, and other annotations not clearly known to produce functional proteins, 2778 genes remained in the differentially expressed category (Supplemental Table 2). All further analyses were done on this set of genes.

TARGETED EVALUATION OF MICROARRAY EXPRESSION DATA
To compare statistically significant patterns of gene expression from the array with those from qPCR, we measured qPCR expression of AP2-ERF transcription factors (N = 17) in the caterpillar treatments at 6 and 24 h in both local and systemic tissue because of their species-specific expression pattern, (Rehrig et al., 2011) for a total of 136 measurements. There were many cases (36%) in which the more sensitive qPCR detected statistically significant changes in expression not detected by the array. However, a majority of those identified by the array as significant were confirmed by qPCR as significant (20 of 26). Four of the 6 array false positives had qPCR values for relative expression in the same direction as those of the array even though they failed the test of significance, a more stringent criterion than most authors apply (Rehrig et al., 2011).

CLUSTERING AND OVERLAP ANALYSES
The overall pattern of similarity and difference in gene expression among the treatments was identified with a simple clustering algorithm (VARCLUS, SAS) that used the centroid method and the maximum number of possible clusters set to 16, or the number of treatments. The percentage overlap in gene expression between treatments was calculated as the (# genes in common/(sum # genes elicited by both treatments-# genes in common).
To identify functional patterns in the large number of genes differentially expressed in response to insects, we used the DAVID gene enrichment analysis tool (Huang et al., 2008). The gene lists for each treatment consisted of only those that were statistically up or down regulated in that treatment compared to controls (Supplemental Table 3) and the classification stringency was set to medium. Functional clusters were identified based on the modified Fisher exact p-value (EASE score) and clusters with a Group Enrichment Score > 1 were examined further. A Benjamini-Hochberg p < 0.05 correcting for multiple comparisons was considered to be statistically enriched; values significant at the uncorrected p-value were also indicated in some comparisons. Similar analyses were run using MAPMAN and PAGEMAN but no significant enrichments were detected with a Benjamin-Hochberg correction for multiple comparisons. Thus, we concluded that the large number of genes of unknown function in this study made DAVID a more useful tool because it incorporates a larger number of databases.
The Hormonometer Tool (Volodarsky et al., 2009) was used to evaluate the similarity in expression profiles elicited by insects with the published, indexed list of those elicited by exogenous application of plant hormones. Although exogenously applied hormones may not represent in vivo levels of hormones, this kind of comparison is frequently used to identify the importance of hormone signaling pathways in plant responses. Arabidopsis gene identities (AGI) were converted to Affymetrix GeneChip identities using the "at to AGI converter" tool (The Bio-Analytic Resource for Plant Biology, http://bar.utoronto.ca/) (Supplemental Table 4). In a few cases there was no correspondence and the AGI data were omitted, while in a few other cases there were two GeneChip IDs for one AGI and the lines were duplicated and retained.
Groups of genes that are coregulated in response to treatments were identified with a weighted gene coexpression network analysis (WGCNA) as described previously (Zhang and Horvath, 2005;Weston et al., 2008Weston et al., , 2011. Background corrected signal intensities underwent variance stabilization and normalization followed by log 2 transformation before entering network construction (Huber et al., 2002). The WGCNA analysis consisted of 4 steps: (1) creation of a pair-wise Pearson correlation matrix for all genes across all treatments; (2) transformation of correlations to connection strengths (connectivity) using a signed power adjacency function (Zhu-Salzman et al., 2004); (3) identification of modules, or groups of highly correlated gene expression patterns by coupling hierarchical clustering with topological overlap matrix; and (4) relating external gene or treatment information to network properties. Correlations were corroborated using a random seed permutation t-test with 106 iterations.
The coexpression landscape of genes in functionally enriched categories in DAVID analysis was depicted with comparison to the ATTED-II database whose software generates networks by coexpressed ranks calculated from the 1388 GeneChip data (Obayashi and Kinoshita, 2010). AGI codes for genes of interest were pasted into the Network Drawer and the default settings were used to depict the network.

IDENTIFICATION OF PROMOTER REGIONS
Several on-line databases and bioinformatics tools were used to conduct three separate analyses to identify potential cis-elements involved in the signaling pathways after insect attack. Gene sequences up to 1000 bp upstream of the AGI transcription start site were downloaded from The Arabidopsis Information Network (TAIR; www.arabidopsis.org) Sequence Database for all genes whose expression was significantly affected by insect feeding or mechanical wounding. First, the database of motifs found in plant cis-acting regulatory DNA elements (PLACE; Higo et al., 1999) was used to search the downloaded sequences for all known transcription factor binding sites and cis-elements in known gene promoters. A customized Perl script was used to tally all elements found in each differentially regulated gene. Principal component analysis (PCA) and cluster analysis was performed with the tallied data to identify specific cis-element fingerprints possibly unique to insect treatments. Treatments were clustered on the basis of similarities of their cis-element distributions using the SAS VARCLUS procedure (SAS Institute).
Second, the ATHENA search tool (O'Connor et al., 2005), which uses a library of 105 known Arabidopsis transcription factor binding sites and 30,067 predicted promoters, was used to search the downloaded sequences for enriched, known transcription factor binding sites and cis-elements in gene promoters. The ATHENA algorithm conducts a student's t-test to determine whether motifs in a given sample set are significantly different from a random distribution in the genome. Motif occurrences with a p < 0.00001 were designated as "enriched." Third, MotifSampler (Thijs et al., 2001(Thijs et al., , 2002 was used to identify all potential (known and putative) transcription factor (TF) binding sites and enriched sequence consensus within the promoter regions in the downloaded, upstream sequences of differentially regulated genes. This tool successfully identified potential TF binding sites in a number of different species (Singh, 1998;Chen et al., 2002). Due to high nucleotide substitution rates, we wrote a short PHP script to place all combinations of A, C, G, or T's in place of degenerative nucleotide indicators (e.g., K, S, N, V). Resolved sequences were then compared against the ATHENA results.

TRANSCRIPTIONAL PROFILES WERE DISTINCT, EVEN WITHIN FEEDING TYPES
We assessed the transcriptional responses in rosette leaves of 4week-old A. thaliana Col plants to artificial wounding, attack by larvae of two species of leaf chewing caterpillar species and by adults and nymphs of two phloem feeding aphid species (Table 1). In all, 2778 genes on the full genome wide array were differentially expressed in one or more of our treatments, representing approximately 10.6% of the genes on the array having AGI annotation (Supplemental Table 2). Hierarchical clustering based on similarities (correlations) among the resulting transcriptional profiles (Figure 1) revealed three major clusters, one containing all aphid treatments, a second comprising wounding plus caterpillar treatments at 24 h, and the third comprising wounding plus caterpillar treatments at 6 h. The aphid cluster contained two subclusters based on sampling time. The wounding plus FIGURE 1 | Hierarchical cluster tree of A. thaliana genes differentially expressed in response to insect feeding and mechanical wounding. Abbreviations for treatment names are the same as those used in Table 1. caterpillar cluster divided into two subclusters, one for wounding and another for the two caterpillar species; the latter contained sub subclusters for each species. The third cluster, responses to wounding and caterpillars at 6 h, divided into two subclusters, one for wounding and another for the insects, with the insect sub subcluster further divided by species. Hence feeding type (sucking vs. chewing) and sampling time (6 and 24 h) contributed most (about 40-60% of variance explained) to differentiating the transcriptional response profiles, whereas species contributed some (about 20% of variance explained) to differentiation and tissue treatment (wounded or not) contributed least (<10% of the variance) (Figure 1).
Arabidopsis responses to treatments were highly dynamic (Figure 2). Since sampling was done 6 and 24 h after the removal of insects, the differences in expression between insects and control plants presumably represents a combination of induction and relaxation of gene expression responses. Across all treatments and tissues, overlap in genes differentially expressed at 6 and 24 h ranged from 3 to 15% and 1 to 8% for upregulated and downregulated genes, respectively (Figure 2A). The aphid M. persicae elicited far more differentially expressed genes at 6 h (n = 815, 45% up and 54% down) than at 24 h (n = 222, 74% up and 38% down) (Figure 3). In contrast, the aphid B. brassicae elicited fewer differentially expressed genes at 6 h (n = 188, 44% up and 56% down) than at 24 h (n = 249, 76% up and 23% down) Transcriptional responses to the two aphids across both times only shared a 4-8% overlap in differentially expressed genes. Both caterpillars elicited similar numbers of differentially expressed genes at 6 h and 24 h (Figures 2A, 3). However, the degree of overlap in these genes between the caterpillars differed with time, 24% up vs. 9% down at 6 h and 8% up vs. 4% down at 24 h for S. exigua and P. rapae, respectively. Feeding by S. exigua elicited the longestlasting changes in gene expression and hence the most overlap in genes differentially expressed at both time points.
Response profiles differed between feeding type treatments. Aphids clustered with aphids and caterpillars with caterpillars plus wounding (Figure 1). Chewing caterpillars and sucking aphids elicited in common only 10% of upregulated genes and 8% of downregulated genes (Figure 3). Although profiles elicited by a particular caterpillar species clustered together (Figure 1), responses of specific genes elicited were species-specific to a large extent (Figure 3) with only 21 and 12% of up-and downregulated genes shared, respectively. Expression profiles elicited by the two aphids shared only 7 and 4% of all upregulated and downregulated genes, respectively. Caterpillars elicited roughly 3 times more upregulated genes and twice as many downregulated genes as the mechanical wounding. The amount of overlap between responses to caterpillars and mechanical wounding was similar to that observed between caterpillars and aphids (Figure 3).
The fewest transcriptional differences in response to caterpillars occurred between leaf types (attacked or "local" and unattacked or "systemic") ( Figure 2B). More genes were differentially expressed in leaves attacked by S. exigua compared with unattacked leaves on the same plant, whereas P. rapae differentially expressed roughly equal numbers of genes in both leaf types. Overlap in genes differentially expressed in the leaf type by the caterpillars also differed, with 19 and 9% overlap in local leaves vs. 16 and 13% overlap in systemic leaves for upregulated and downregulated genes, respectively. Data from both caterpillar species was combined by leaf type and compared to leaves wounded locally and systemic unwounded leaves. Few of the same genes elicited by caterpillars in the attached and unattacked leaves were elicited in both wounded and unwounded leaves at the same time points; 11-16% of upregulated and 2-15% of downregulated genes were common depending on caterpillar species leaf type and sampling time.
Transcriptional responses of Arabidopsis to herbivores were related to the insects' dietary breadth. The two dietary generalists upregulated more genes than did their respective dietary specialists (p < 0.0219, Supplemental Table 5). There was no difference (p > 0.2715) between them for downregulated genes.

TRANSCRIPTIONAL PROFILES WERE ENRICHED IN STRESS RESPONSES AND SHIFTS IN PRIMARY METABOLISM
The DAVID Functional Clustering Tool was used to identify biological functions significantly enriched among differentially expressed genes. About two thirds of the genes differentially expressed in response to our treatments have some functional or structural information in TAIR. DAVID used TAIR and added information from other species for genes with little annotation in TAIR. Although most of these ontological data were collected in contexts other than responses to insects, we can use them to infer some functions of proteins encoded by genes influenced by our treatments.
Transcriptional responses to herbivory overlapped broadly with responses to 14 other biotic and abiotic stressors and stimuli, including responses to wounding, pathogens, cold, starvation, nutrient supply, metal ions, osmotic stress, reactive oxygen species, radiation, and light quantity and quality ( Table 2). Caterpillars elicited the greatest number of different stress/stimulus responses, including many that were not elicited by wounding or aphids at either time point. Pieris rapae elicited transcriptional changes in 8 and S. exigua in 10 stimulus pathways linked to our treatments by DAVID at least once; each species elicited only two changes by the 6 h sampling. The two aphid species elicited fewer stress-related responses than did the caterpillars. Brevicoryne brassicae elicited changes in 4 stress-related pathways while M. persicae altered expression in 6 of the 14 relevant pathways. The transcriptional gene set most often altered by all the insects was response to hexose stimulus (6 of 16 possible interactions), but no single response gene set in this analysis characterized responses to insect elicitation as different from wounding or other stressors.
Herbivory caused widespread changes in expression of genes involved in primary metabolism ( Table 3). Both caterpillars and aphids elicited changes in expression of genes associated with metabolism of organic acids, fatty acids, lipids, and amino acids (Tables 3A-C). Expression of genes in metabolic pathways related to 12 amino acids was altered at least once by our treatments. P. rapae elicited changes in 11 of the 12; no other insect treatment elicited changes in more than five. Both caterpillars and aphids elicited changes in expression of genes associated with cell wall metabolism, although their impacts were different (Table 3D; Supplemental Figure 4). Aphids upregulated five times as many cell wall-related genes as did caterpillars, including a large number encoding extensins and extensin-like proteins, several peroxidases, and many cell wall-associated heat shock proteins and cognates. In contrast, caterpillars downregulated twice as many cell wall genes as did aphids, including those encoding galactosidases, xylosidases, and pectinases.
Caterpillars and aphids caused very different changes in the expression of genes associated with amino acid metabolism. Caterpillars, but not aphids, altered the expression of genes involved in biosynthesis of indole and aromatic amino acids (Table 3C). Caterpillars upregulated three genes directly involved in TRP biosynthesis, three involved in indole glucosinolate (GS) synthesis, and one in auxin biosynthesis (Tables 5A,B). One gene involved in indole GS synthesis was downregulated by both aphids at 24 h (Table 5B).
Consistent with upregulation of genes involved in biosynthesis of amino acids from which GS are derived, caterpillars also uniquely influenced expression of genes specifically associated with GS biosynthesis and activity (Table 3D). These included a wide variety of genes in nitrogen and sulfur assimilation and modulation of GS activity (Supplemental Table 6). Responses to caterpillars were functionally enriched in genes www.frontiersin.org November 2014 | Volume 5 | Article 565 | 9  involved in glucosinolate metabolism; responses to aphids were not (Tables 3G, 5). Genes involved in glucosinolate biosynthesis were upregulated by S. exigua in damaged leaves at both time points and by P. rapae at 24 h in both leaf types (Tables 5B-D). Genes involved in glucosinolate catabolism, e.g., nitrile biosynthesis genes that promote simple nitrile, but not epithionitrile or thiocyanate formation, were upregulated by both caterpillar species in damaged tissue at 6 h ( Table 5E). Aphids failed to alter the expression of enough glucosinolate-related genes to generate a statistically significant enrichment pattern.

TRANSCRIPTIONAL PROFILES WERE ENRICHED IN GENES ASSOCIATED WITH HORMONE SIGNALING
To determine whether particular hormone signaling pathways were involved in responses to insect attack, the Hormonometer Tool was used to evaluate the similarity in expression profiles elicited by insects to those elicited by exogenous application of plant hormones. The early, 6 h wound treatments and all of the caterpillar treatments and times strongly upregulated genes associated with responses to methyl jasmonate (MeJA) in all leaf types (Figure 4). There was some evidence that the aphid B. brassicae downregulated MeJA-related responses. Caterpillar treatments (but not wounding or aphids) elicited expression of genes associated with responses to abscisic acid (ABA) in all leaf types at both times (Figure 4). All of the insect treatments tended to downregulate genes associated with early responses to ethylene and cytokinins and upregulated genes associated with salicylic acid (SA; Figure 4). These patterns were generally supported by the DAVID analysis (Table 3G) although that analysis found fewer significant associations.

TRANSCRIPTIONAL PROFILES INDUCED BY CATERPILLARS DIFFERED IN LOCAL AND SYSTEMIC LEAVES
We used the DAVID tool to examine functional differences between transcriptional responses of two leaf types; leaves damaged directly ("local") by caterpillars and age-matched undamaged ("systemic") leaves on the same plants. DAVID identified 8 functional categories that were significantly enriched in one or more combinations of leaf type, caterpillar and time ( Table 6, Supplemental Table 6). Two of these were found in a single combination; response of genes involved in water deprivation were upregulated only at 6 h in leaves damaged by S. exigua and responses of genes associated with lectins were upregulated only at 24 h in leaves damaged by P. rapae. Ethylene signaling pathway genes were upregulated in both local and systemic leaves only at 24 h by S. exigua. Responses of genes involved in sulfur metabolism were upregulated at both time points in leaves damaged by S. exigua. Both species elicited responses in genes involved in auxin signaling at 6 h, but in different leaf types, systemic for P. rapae and local for S. exigua. Most of the auxin signaling-related FIGURE 4 | Hormonometer analysis of differential gene expression by A. thaliana in response to insect feeding and mechanical wounding. The response in gene expression of Arabidopsis in our treatments is compared with that of Arabidopsis at 30, 60, and 180 min, or 3, 6, and 9 h, after hormone application. The magnitude of correlation in gene expression is indicated by the color scale at bottom and correlation values of 0.4 and above are considered statistically significant. Treatment abbreviations are the same as in Table 1. MJ, methyl jasmonate; ACC, 1-aminocyclopropane-1-caroxylic acid; ABA, abscisic acid; IAA, indole-3-acetic acid; GA3, gibberellic acid 3; zeatin, cytokinin; brassino, brassinosteroid.   genes in the enriched functional category were downregulated by insects (Supplemental Table 6). Both caterpillars elicited expression of genes associated with responses to wounding and jasmonates in all leaf types at 6 h ( (Supplemental Table 6). Genes associated with responses to both wound and JA continued to be elicited in all leaves at 24 h in S. exigua treatments. Expression of a subset of jasmonateand wound-responsive genes was consistently upregulated by S. exigua, but not P. rapae, at both time points, including ERF4 (AT3G15210), ORA47 (AT1G74930), WRKY40 (AT1G80840), AP2C1 (AT2G28900), CCR-4-NOT (AT3G44260), SYP122 (AT3G52400), and STZ (AT1G27730) (Figures 5A,B). This duration of ethylene-responsive ERF4 and ORA47 expression implicates ethylene signaling. Indeed, both leaf types at 24 h for S. exigua elicited expression of genes associated with the ethylene signaling pathway at higher levels than did P. rapae (Figure 5C).

COEXPRESSION NETWORK ANALYSIS IDENTIFIED MODULES OF COEXPRESSION AND PUTATIVE HUB GENES
The weighted gene coexpression analysis identified a nonrandom network architecture of 18 modules of coexpressed genes (Supplemental Figure 4a). All but one of the modules was statistically associated with treatments, and there was no overlap in modules responding to caterpillars, mechanical wounding, and aphids ( Table 7). For caterpillars, the "cyan," "blue," "light yellow" and "tan" modules were associated with responses to P. rapae, whereas the "light cyan," "green yellow," "red," "gray60" and "dark red" modules were associated with responses to S. exigua. Only one of the modules (turquoise) was associated with mechanical wounding (local and systemic responses at 6 h). For aphids, the "yellow," "salmon" and "brown" modules were associated with responses to B. brassicae, whereas the "black," "dark turquoise," "midnight blue," "brown" and "dark green" modules were associated with responses to M. persicae.
We then looked to see if any species-specific modules were also characteristic of caterpillars or aphids in general. The "red," "tan" and "gray60" modules were associated with 6 h responses to caterpillars, and the "light cyan" and "light yellow" modules were associated with 24 h responses to caterpillars. The "black," "dark turquoise," "midnight blue," "yellow" and "salmon modules were associated with 6 h responses to aphids," and the "brown" and "dark green" modules were associated with the 24 h responses to aphids. These patterns are reflected in the three main arms of the network architecture, comprising the 6 h responses to wounding ("turquoise"), caterpillars ("red"), and aphids ("midnightblue" ) (Supplemental Figure 4b).
Patterns of gene function within modules were determined by examining the ontological enrichment of their member genes. Only a few of the 17 modules associated with treatments were enriched in specific gene functions, which may reflect the absence of functional annotations. To identify regulatory hubs that may direct treatment-specific responses, only module membership values ≥0.8 were used to recognize highly connected module members. These results are summarized for all modules in Supplemental Table 7. Although all modules contained one or more highly connected genes, only a few of those genes are characterized, and are discussed below. A majority of the potential hub genes are either of unknown function or have only an hypothesized function based on sequence similarity.
The "green yellow" module was statistically associated with responses to S. exigua at 6 h in unattacked, systemic leaves  Table 1.
( Table 7). Ontological analysis indicated that this module was enriched in genes associated with responses to JA, bacteria, and chitin (Supplemental Table 8). Highly connected genes in this module were examined for their potential to serve as Only modules with eigenvalues > 0.35 were considered significantly associated with a treatment. Asterisks denote the p-value for the given eigenvalue: *p < 0.05, **p < 0.01, ***p < 0.005, ****p < 0.001, *****p < 0.0005, ******p < 0.0001. Abbreviations for plant treatment names are the same as those used in Table 1. www.frontiersin.org November 2014 | Volume 5 | Article 565 | 13 hub genes directing treatment-specific plant responses, and four (At1g80840, At5g54145, At5g35475, At1g57990) had module memberships ≥ 0.8 suggesting that they have potential regulatory importance. WRKY40 (At1g80840) is a transcription factor that binds W-box sequences and forms protein complexes with itself and with WRKY18 and WRKY60 to influence Arabidopsis susceptibility to biotic and abiotic stresses (Xu et al., 2006;Chen et al., 2012). The gene At5g54145 is an expressed protein of unknown function, while At5g35475 is similar to a hAT dimerization domain-containing protein/transposase-related. The At1g57990 gene is related to the purine transporter PUP1, which may be involved in the transport of purines and purine derivatives like cytokinins across the plasma membrane.
To infer the functions of these potential hub genes, we used ATTED to look for gene subnetworks they anchor. However, only WRKY40 was present on the Affymmetrix arrays used to construct the ATTED network (Obayashi and Kinoshita, 2010). In this network, WRKY40 interacts with three genes identified with plant responses to wounding (CCR4-NOT, AT3G44260; BCB, At5g20230; RHL41/ZAT1, AT5G59820) and two genes identified with plant responses to ethylene (ERF11, At1g28370; ACS6, At4g11280) (Supplemental Figure 5). WRKY40 links many such genes elicited by S. exigua but not the other insects in our study, and is thus a candidate hub gene.
The network coexpression analysis also identified several other potential hub genes for early Arabidopsis responses to M. persicae. Two modules ("dark turquoise" and "black") were statistically associated with responses to M. persicae at 6 h ( Table 7). The dark turquoise module was enriched in genes associated with extensin-like activity, including many of the cell wall associated genes upregulated by M. persicae at 6 h (Supplemental Figure 4, Supplemental Table 8). There were two highly connected, potential hub genes in the dark turquoise module, a calcium-binding EF-hand family protein (At3g47480) about which little is known, and VBF1 (At1g47056), a VIER F-box protein known to be a positive regulator of auxin response and cell wall metabolism (Supplemental Table 8; Schwager et al., 2007). Although there was no functional enrichment of all of the genes in the black module, the potential hub genes were enriched in an annotation of "extracellular region." These genes included a putative peroxidase (At4g37520), two disease resistance proteins (SOBIR1, At2g31880; CC-NBS-LRR, At4g33300) and three genes upregulated by M. persicae at 6 h and associated with cell wall metabolism (Supplemental Figure 4; EXP4, At2g39700; HSP81.4, At5g56000; GUS2, At5g07830).

Cis-ELEMENTS OF INSECT REGULATED GENES EXPLAIN ONLY SOME OF THE TRANSCRIPTIONAL DIFFERENCES
We conducted three separate analyses to identify potential ciselements that are involved in TF signaling pathways responding to our treatments. First, we found that most genes whose expression was altered by insect feeding shared the great majority of known or postulated TF binding sites (Figure 6). A cluster analysis found that the great majority of these genes contained all of the known cis-elements. However, 2 clades appeared distinct. The cis-element composition of genes down-regulated by M. persicae formed a single clade, and cis-regulatory sequences of genes down-regulated by the two caterpillars, S. exigua and P. rapae, comprised another. Cis-elements of genes involved in responses to wounding did not form a distinct clade.
ATHENA enabled identification of eight cis-elements that were over-represented (p < 0.00001) in genes co-expressed in response to individual treatments (Table 8). These elements may serve as important regulatory components of signaling pathways responding to insect attack. Several of the enriched cis-elements belong to well characterized classes, including ABREs, G-Boxes, TATA boxes, WRKY boxes, and I-boxes. TATA boxes were overrepresented in many of the genes upregulated by caterpillars, while I-boxes, found in many light-regulated genes (Giuliano et al., 1988), were enriched in genes downregulated by caterpillars. Interestingly, most treatments repressed genes that were enriched in "Evening Element" motifs. This cis-element is characteristic of many genes whose expression peaks at dusk (Harmer et al., 2000).
We also conducted a cis-element analysis using the MotifSampler tool, which unlike ATHENA, can detect enrichment of both unknown and known motifs (Supplemental Table 9) (Thijs et al., 2001(Thijs et al., , 2002. However, there was a high nucleotide substitution rate due to sequence consensuses, so we wrote a script that resolved all degenerate nucleotide sequences, predicted all possible 6-, 7-, and 8-mers identified by MotifSampler, and then matched them against the Athena cis-element database. This process enabled us to confirm 6 of the 20 putative elements identified by MotifSampler and relate them to our ATHENA results ( Table 8). The reason why there is not 100% confirmation is due to differences in methodology between MotifSampler and ATHENA. Motif Sampler tries to find over-represented motifs in the upstream region of a set of co-regulated genes. This type of motif finding-algorithm uses Gibbs sampling to find the position probability matrix that represents the motif and uses background models to improve the robustness of the motif finding. Therefore, if no significant difference exists between background and entered list, the motifs will not be predicted due to lack of higher scores. Nonetheless, these motifs would have to be experimentally verified to confirm their role in insect responses.

SOURCES OF SPECIFICITY IN PLANT RESPONSES TO HERBIVORES AND WOUNDING
Transcriptional responses to insects are distinct from responses to wounding in both number and identities of responsive genes.
Here we focused our analysis on gene enrichment to identify changes in gene expression of groups of functionally related genes, to address the limitations of microarray technology (false positives and false negatives) which require that expression changes of individual genes be confirmed. The larger number of genes upregulated by the two dietary generalists is consistent with the commonly held hypothesis that dietary specialists have evolved a stealthier way to exploit their host plants than dietary generalists (Ali and Agrawal, 2012). But, as others have noted, it is unwise to draw conclusions from studies such as ours that lack many examples of insects of each dietary type (Bidart-Bouzat and Kliebenstein, 2011).
Insect herbivores provide potential signals to plants that reveal their identity by feeding at different times of the day, on different kinds of plant tissue, on tissue of different ages, causing different damage patterns and amounts, and/or any chemistry they introduce with oral secretions (Bonaventure et al., 2011;Maffei et al., 2012). These signals generate HAMPs that plants can use to identify their attackers and respond appropriately. Evidently HAMPs set insects apart even though insect feeding may involve extensive wounding. Therefore, given their very different feeding styles, it is not surprising that aphids and caterpillars used in our study elicited changes in gene expression with only about 10% of genes in common in this study.
Aphids were given free range of the plant in our study and usually chose to feed on phloem of younger leaves. Although there was little difference between the two aphid species in feeding sites, about 90% of the genes whose expression they elicited were different from each other. Thus, the specific searching behavior of their stylets within the plants and the oral secretions they release are likely sources of differentiating signals . Aphid stylets secrete a variety of proteins as they make their way through and around other cells in the leaf to reach the phloem. The list of proteins is extensive and includes lipases, peroxidases, pectinases, and glucosidases as well as many proteins of unidentified function (Miles, 1999;Elzinga and Jander, 2013). Several have been identified from M. persicae, but the paucity of functional studies makes it premature to suggest which may be involved in the different responses we observed.
In the caterpillar experiments, we kept tissue and age constant by restricting feeding to specific, similarly aged leaves and had caterpillars feed until they removed similar amounts of leaf area. Despite this, only a little more than 20 and 10% of upregulated and downregulated gene expression, respectively, were elicited in common by both species. The caterpillars exhibited some differences in feeding pattern, with P. rapae eating continuously in a single location, usually on a leaf margin, whereas S. exigua usually fed at several locations on the leaf and rarely from the margin. As a result, the observed differences in the effect of these two species on gene expression could have been caused by differences in their patterns of damage, their oral secretions deposited at the feeding site or both. Oral secretions of caterpillars are known to modulate plant responses to mechanical wounding. Fatty acidamino acid conjugates (FACs) and glucose oxidase (GOX) are in www.frontiersin.org November 2014 | Volume 5 | Article 565 | 15 FIGURE 6 | Heat map of Cis-element distributions in genes differentially expressed by insect herbivory and mechanical wounding. Transcription factor binding sites up to 1000 bp upstream of all differentially expressed genes were located using the PLACE database and clustered by treatment. Copies of a motif found is indicated by the color scale; e.g., red, many; black, none.
Considerably more work will be required to identify the many HAMPs likely present in diverse insect species.

STRESS RESPONSES AND RECONFIGURATION OF METABOLISM
Expression of genes typical of responses to other stressors is a common feature of plant responses to herbivory. Most of these responses involve a reconfiguration of primary metabolism (Schwachtje and Baldwin, 2008). The changes we observed in the Arabidopsis transcriptome generally reflect downregulation of key metabolic components following insect attack. The detection of stress caused by insect feeding is thought to be mediated by SNF1/AMPK/SnRK1 protein kinases, which as energy sensors serve as hubs for generalized stress signaling (Polge et al., 2008;Crozet et al., 2014). In Nicotiana attenuata, a SnRKkinase regulates reallocation of photoassimilates in response to chewing herbivores (Schwachtje et al., 2006). SnRK1 also likely mediated downregulation of four TPS genes by M. persicae in our study. Although the SnRK1 gene itself was not differentially expressed in our experiments, M. persicae feeding at 6 h downregulated the expression of 9 of 24 known targets of SnRK1 (Supplemental Figure 3; Baena-Gonzalez and Sheen, 2008). Downregulation of photosynthesis is a common response to biotic stress and has been previously reported to also occur in response to herbivory (Zangerl et al., 2002;Giri et al., 2006;Tang et al., 2006;Bilgin et al., 2010). Although we did not measure photosynthesis in our study, we found caterpillars and aphids to have different impacts on gene expression associated with photosynthesis, with M. persicae upregulating Calvin Cycle genes and the caterpillars having little effect.
Some of the plant transcriptional responses to insects involved hexoses. Previous work also has shown that chewing herbivores can stimulate localized sink strength in the tissues they attack (Arnold et al., 2004;Appel et al., 2012). Hexoses provide both signals and substrates for these and other defense responses (Schultz et al., 2013;Tauzin and Giardina, 2014). While we did not measure sink strength in these studies, the results support the view that hexose signaling and modified source-sink relationships are part of the plant-insect interaction.
Downregulation of protein synthesis is a common plant response to stress and is accomplished by phosphorylation of the alpha subunit of EUKARYOTIC INITIATION FACTOR 2 (eIF2alpha) by the kinase GENERAL CONTROL NON-REPRESSIBLE 2 (GCN2; Immanuel et al., 2012). Upregulation of GCN2 by M. persicae may represent a core regulatory response, but the downstream effects of eIF2alpha are not well enough understood in Arabidopsis to confirm its role in plant responses to insects (Li et al., 2013;Luna et al., 2014). Downregulation of protein translation is triggered by stress on the endoplasmic reticulum and likely involves the unfolded protein response (Duwi Fanata et al., 2013). This scenario may explain the increase in transcription of heat shock proteins observed in many transcriptome studies, including our own (Heidel-Fischer et al., 2014). However, downregulation of all protein synthesis following tissue damage is unlikely to be adaptive for the plant. Indeed, we observed in response to caterpillar feeding the upregulation of defense genes encoding proteins involved in biosynthesis and activity of glucosinolate biosynthesis, including their amino acid precursors. Similarly, the upregulation by aphids and downregulation by caterpillars of genes associated with cell wall metabolism are likely to reflect their differential impact on the TOR pathway controlling plant cell growth (Leiber et al., 2010;Henriques et al., 2014).
Fine tuning of metabolic programming to modulate investment in growth vs. defense has been reported previously (Schwachtje and Baldwin, 2008;Heidel-Fischer et al., 2014) and references therein. In N. attenuata, deactivation of ribulose bisphosphate carboxylase (RuBPCase) by RUBPCAse activase mediates the change in primary metabolism which in turn attenuates JA-induced defense responses following chewing insect attack (Mitra and Baldwin, 2014). In Arabidopsis, DELLA proteins are also involved in attenuation of the JA response (Lan et al., 2014) but a role for RuBCase has not yet been examined.

TRANSCRIPTIONAL SIGNATURES OF HORMONE SIGNALING
It has become almost a truism to attribute plant responses to sucking insects as arising primarily from SA signaling and plant responses to chewing insects as arising primarily from signaling by JA and ethylene. However, the reality is much less clear-cut (Verhage et al., 2011;Erb et al., 2012) and signaling differences may in part be a matter of relative timing of herbivory along with importance of the herbivore Rehrig et al., 2014) and a matter of ABA signaling. In this study, both chewing and sucking insects elicited the expression of genes associated with SA signaling, consistent with some reports of positive crosstalk between the two pathways (Schweiger et al., 2014) and the susceptibility of sucking insects to changes in JA signaling (Ellis et al., 2002;Mewis et al., 2005). However, the pattern was not reciprocal for JA; only caterpillars and mechanical wounding elicited expression of genes associated with JA signaling. Caterpillars also elicited the expression of genes associated with ethylene signaling, with evidence of differential induction of ethylene by the two species used in this study (Rehrig et al., 2014 this volume). Caterpillars uniquely elicited expression of genes associated with ABA signaling, which in N. attenuata is due in part to the inhibitory effect of an elicitor in caterpillar oral secretions on ABA catabolism (Dinh et al., 2013). Since ABA interacts with the SnRK1 sensors (Rodrigues et al., 2013), it is another likely source of hormonal crosstalk in plant responses to herbivores.

Cis-ELEMENTS OF INSECT REGULATED GENES EXPLAIN ONLY SOME OF THE TRANSCRIPTIONAL DIFFERENCES
Given evidence of insect-specific plant responses, we hypothesized that the cis-binding element composition of treatmentaffected genes might display "fingerprints" characteristic of the treatment. No known regulatory elements have been found to be insect-specific in plants and little evidence is available about gene regulatory networks in plant-insect interactions (Zou et al., 2011). Some have proposed that a substantial part of the patterns in any gene expression data could be explained by "ciselement profiles" (Beer and Tavazoie, 2004;Segal et al., 2005). Using computational predictions with yeast stress data, some have argued that motif profiles in conjunction with microarray analysis can help identify important regulatory networks including those involved in the pathology of cancer development (Segal et al., 2005). There are only a few studies where this has been applied to plants. For example, using computational analysis and publically available Arabidopsis microarray databases, over 53 putative motifs involved in phytohormone signaling in plants were discovered (Yamamoto et al., 2011). The authors suggest using computational models as starting points for hypothesis formation before beginning wet-bench experiments. With this in mind, we characterized all of the cis-elements in promoter regions 1000 bp upstream of all genes affected by treatments using three different bioinformatic tools. Our results showed that the genes differentially expressed in any treatment differed in only a few known cis-elements in their promoter regions. However, a few elements were associated specifically with particular treatments and differed from randomly-selected TF gene profiles. First, cis-elements of TF genes down-regulated by M. persicae were quite different from any other set, consistent with the distinctive patterns of gene expression it elicited. Second, ABRE and ABRE-like elements were over-represented in most insect response treatments. This association is consistent with the role of ABREs in regulating the expression of genes that are ABA and drought-responsive (Fujita et al., 2005) and with the enrichment of ABA-responsive genes we observed in our caterpillar treatments. Although we did not detect a significant change in gene expression associated with ABA signaling in the aphid treatments of our study, others have reported them (Morkunas et al., 2011).
Third, TATA boxes were enriched in genes upregulated by both caterpillars; these elements may also be important regulatory components in biotic stress. Genes of Saccharomyces sp. with TATA boxes were found to be overrepresented in stress responses to abnormal osmolarity, pH balance, or nutrient availability (Basehoar et al., 2004). Genes without TATA boxes performed more constitutive housekeeping functions and may not need as much transcriptional regulation. The presence of TATAbinding motifs may identify Arabidopsis genes that are especially responsive to environmental stressors, especially chewing insects.
Taken together, these data suggest that several cis-elements distinguish the Arabidopsis gene regulatory networks involved in responses to insect attack, but that the promoter regions of responsive genes contain mostly the same elements. It is also likely that regulation of transcriptional responses to insects via TFs could be combinatorial, requiring multiple TFs to initiate particular expression patterns (Singh, 1998;Lindlöf et al., 2009;Zou et al., 2011).
In summary, our results provide a comprehensive overview of Arabidopsis transcriptional responses to chewing and sucking insects. We identified expression patterns in genes and cis-elements over-represented in co-expressed genes after insect attack that provide insight into the complex networks and regulatory pathways insects elicit in plants during herbivory. The results provide definitive evidence that plant responses to insects are distinct from and considerably more complex than are responses to wounding alone, even though insects wound as they feed. Thus, HAMPs are more important than DAMPS in determining plant responses to insects. Analyses of coexpression patterns and known interaction networks suggest that characterizing "master regulators" of plant responses to stress, including NPR1, Mediator Complex (Med25), and SnRK1 and TOR kinases is likely to provide additional insight (Robaglia et al., 2012;Balderas-Hernández et al., 2013;Sheen, 2014).

AUTHOR CONTRIBUTIONS
Heidi M. Appel-Helped design array and insect experiments, conducted insect experiments, prepared manuscript, analyzed and interpreted data. Howard Fescemyer-Helped design microarray experiments, prepared RNA for microarray analysis, edited manuscript. Juergen Ehlting-Helped design and conduct microarray experiments. David Weston-Conducted network cluster analysis. Erin Rehrig-Conducted qRT-PCR validation of microarray and cis-element analysis, edited manuscript. Trupti Joshi-Provided advice on cis-element analysis, edited manuscript. Dong Xu -Provided advice on cis-element analysis, edited manuscript. Joerg Bohlmann-Helped design microarray experiments, edited manuscript. Jack Schultz-Helped design microarray and insect experiments, assisted with the analysis and interpretation of data, conducted hierarchical clustering, edited manuscript, authored funding for project.