Metatranscriptomic Analysis of Multiple Environmental Stresses Identifies RAP2.4 Gene Associated with Arabidopsis Immunity to Botrytis cinerea

In this study, we aimed to identify common genetic components during stress response responsible for crosstalk among stresses, and to determine the role of differentially expressed genes in Arabidopsis-Botrytis cinerea interaction. Of 1,554 B. cinerea up-regulated genes, 24%, 1.4% and 14% were induced by biotic, abiotic and hormonal treatments, respectively. About 18%, 2.5% and 22% of B. cinerea down-regulated genes were also repressed by the same stress groups. Our transcriptomic analysis indicates that plant responses to all tested stresses can be mediated by commonly regulated genes; and protein-protein interaction network confirms the cross-interaction between proteins regulated by these genes. Upon challenges to individual or multiple stress(es), accumulation of signaling molecules (e.g. hormones) plays a major role in the activation of downstream defense responses. In silico gene analyses enabled us to assess the involvement of RAP2.4 (related to AP2.4) in plant immunity. Arabidopsis RAP2.4 was repressed by B. cinerea, and its mutants enhanced resistance to the same pathogen. To the best of our knowledge, this is the first report demonstrating the role of RAP2.4 in plant defense against B. cinerea. This research can provide a basis for breeding programs to increase tolerance and improve yield performance in crops.

The least number of DEGs was observed in plants under the two treatments that belong to oxidative stress and wounding (Fig. 1A,B). This could be attributed to the natural adaptation of Arabidopsis to the abiotic stress group in comparison to other stress conditions 10,11 . In Arabidopsis, the gene expression levels upon individual treatments of SA, MeJA, ACC and ABA were altered for 3,051, 2,918, 2,590 and 3,231 transcripts, respectively, from which 1,340 (43.9%), 1,397 (47.9%), 1,328 (51.3%) and 1,564 (48.4%) genes were stress-induced genes. On average, most FC of the DEGs in all stresses ranged between two-fold or three-fold. Interestingly, we noticed that some genes were induced >10 fold or repressed <10 fold. Although each dataset may have its unique DEGs, this does not rule out the possibility that common genes can be identified across the categories of stresses.
To visualize the gene expression data of B. cinerea and other examined stresses, a heatmap of the top 50 BUGs and BDGs was generated (Fig. 1C). The heatmap was combined with clustering methods, displaying a group of genes/samples together based on the similarity of their gene expression pattern. Together, this can be the first step for identifying genes that are commonly regulated or biological signatures that are associated with multiple conditions.
Highly conserved expression of common DEGs to multiple stress responses. A scatter plot was constructed to compare the transcript level of the DEGs of each dataset with that altered by B. cinerea infection (Fig. 2). Clearly, our results demonstrated similar patterns of gene expression levels between Arabidopsis plants infected with B. cinerea and any individual stress at the time of treatment. This suggests that the mode of action of some DEGs common to multiple stresses may contribute in functions or processes that are common among responses to stress acclimation. We also determined the number (and percentage) of genes that were co-regulated in response to B. cinerea and single treatments of each stress type. In response to both the virulent and avirulent strains of P. syrinage, more than 2/3 and 1/2 of the genes were also induced and repressed, respectively, to B. www.nature.com/scientificreports www.nature.com/scientificreports/ cinerea (Table 1). Upon inoculation with B. cinerea, 40-50% of the genes were expressed by either the fungal pathogen A. brassicicola or the herbivory insect P. rapae. We also noticed that between 443 (28.5% in ACC) and 562 (36.2% in ABA) of the induced, and 429 (35.6% in ACC) and 532 (44.1% in SA) of the repressed genes in response to single hormonal stresses were co-regulated with infections by B. cinerea. To lesser extent, 193 (12.4%) and 69  (4.4%) of BUGs and 65 (5.4%) and 47 (3.9%) of BDGs were also up-and down-regulated by paraquat (oxidative stress) and wounding treatments, respectively. Only one gene was induced, and three genes were repressed in response to all individual stresses ( Table 1). The list of identified DEGs that were co-regulated by B. cinerea and the rest of single stresses can be found in Datasets S5-S7 (Supporting Information).  Table 2). The DEGs common to all stresses were further investigated. For example, up-regulated gene, coronatine induced 3/JA responsive 2 (CORI3/JR2; At4g23600), was common in response to all stresses (Table 2; Dataset S8, Supporting Information). The three altered expressed genes by the 11 tested stresses, RAP2.4 (At1g22190), At1g72060 and At2g20670, were found to be commonly down-regulated. Overall, our data emphasize on the complex nature of multiple stress responses and support the importance of studying plant stresses in combination.
The GO annotation was established on commonalities of Arabidopsis upon inoculation with B. cinerea and other pathogens, exposure to abiotic and hormone challenges. Based on the functional similarities of their encoded proteins, BUGs or BDGs were grouped with those up-or down-regulated belonging to other stress classes. According to AGI locus identifiers, 45 functional categories were classified into three major categories upon their biological processes, molecular functions and cellular components (Fig. S4, Supporting Information).  www.nature.com/scientificreports www.nature.com/scientificreports/ Our analysis revealed that the general "response to stresses" was the category of most up-regulated clustered genes when plants were inoculated with B. cinerea and abiotic stress-affected group. The dominant subcategory 'signal transduction' was highly associated with plant defense against pathogens and abiotic cues. The ABA insensitive 1 (ABI1) 31 was up-regulated by B. cinerea and bacterial pathogens (Dataset S5, Supporting Information) as well as the SA and ACC hormones (Dataset S7, Supporting Information). This suggests that plant hormones are tightly associated with defense against B. cinerea and other environmental stresses. Enzymatic activities including kinases were also among the dominant subcategories in BUGs and other groups (Fig. S4, Supporting Information). In addition, "cell wall" term within the cellular component was highly up-regulated in the tested groups; and the wall-associated kinase 1 (WAK1) 11 was also induced by multiple stresses (Datasets S5 and S7, Supporting Information).
Different GO terms of BDGs encoding for structural activity, receptor binding/activity or enzymatic activity proteins were highly down-regulated in response to biotic, abiotic or hormonal stresses, respectively (Fig. S4, Supporting Information). The "electron transport/energy pathways" was the category that most down-regulated genes belonged in response to B. cinerea/biotic stress group and B. cinerea/hormonal stress group; whereas, "plastids" and "ribosomes" were the dominant subcategories in the cellular component. Consistent with previous findings, our data suggest rapid metabolic repression of photosynthetic proteins when plants are under stress including infection with B. cinerea (Fig. S4, Supporting Information) 10 . Many reports have shown that B. cinerea induces/represses several genes particularly those encoding developmental, structural and regulatory proteins in Arabidopsis 5,10,32 . Our data suggest connections between gene expression alternation and mechanisms underlying stress resistance/tolerance during multiple stress exposure.
Our knowledge of biological mechanisms related to B. cinerea infection is still meager 17,22 . Therefore, we constructed interactome networks of commonly regulated genes in silico. Arabidopsis PPI network of DEGs after exposure to groups that belong to multiple environmental challenges and hormonal treatments were annotated by calculating their interactive degrees obtained from the STRING database 33 . The interactive networks displayed clique in the sub-network, suggesting that the hub proteins may form a super complex to play important roles in stress response. As a result, a total of 258 nodes were demonstrated to be involved in network construction and 250 edges were established in the network. Key DEGs, such as those encoding RAP2.4 protein possessed degrees of 10 were markedly more compared with those of other proteins ( Fig. 3; Dataset S9, Supporting Information); and functional analysis of RAP2.4 was further studied. Together, this indicates that the dynamic interactive networks between commonly regulated genes may help decide whether certain network topologies can explain experimental observations. qRt-pcR validation for the microarray results. To confirm the changes in gene expression revealed by the microarray analyses, 13 DEGs were selected based on the microarray data for verification tests (Fig. 4). www.nature.com/scientificreports www.nature.com/scientificreports/ We performed qRT-PCR on Arabidopsis leaves infected with B. cinerea at 0 and 18 hpi. The transcript levels of the DEGs (9 up-regulated; 4 down-regulated) with altered expression in response to various stress treatments (Table 2), were quantified and compared with the obtained microarray analyses. Although the FC values observed in the two expression methods differed somewhat, the tested DEGs displayed comparable patterns in transcript accumulation in the analyses of the two approaches (Table 3; Fig. 4). In response to B. cinerea. For example, the transcript levels of the nine BUGs (CCR2, CYP71A13, α-DOX1, PDF1.2, At1g754300, At3g51660, NATA1, SRG1 and ELI3-2) were increased (Fig. 4), similar to that of the same set of genes identified in the microarray analyses. On the other hand, RAP2.4, DIR-Like, At1g5490 and HAD, which were considered as BDGs, were also repressed by B. cinerea at 18 hpi. In general, a similar trend of expression was found in both the microarrays and the qRT-PCR.

Mutations in RAP2.4 enhanced resistance to B. cinerea.
From the microarray data, RAP2.4 was selected on the basis of its expression profiles for further analysis. Two mutant alleles in RAP2.4 gene displayed very low basal and repressed RAP2.4 expression compared to wild-type plants (Fig. 5A). The rap2.4-1 showed the lowest transcript levels of RAP2.4 at 18 hpi with B. cinerea; and therefore, this T-DNA insertion line was used in further experiments. We assumed that RAP2.4 more likely plays a role in defense. Two homozygous T-DNA insertion lines were inoculated with the fungal pathogen B. cinerea, evaluated for disease resistance and compared with inoculated Arabidopsis wild-type (Col-0) plants -a relatively resistant ecotype.
After inoculations, disease lesions remained restricted in both T-DNA insertion mutant lines of RAP2.4 and exhibited a more resistant phenotype than the wild-type disease phenotype at all time points of infection (Fig. 5B). In wild-type plants, lesions expanded until 4 days post inoculation (dpi), with chlorosis surrounding them. Obviously, when we measured the diameter of the lesions, the rap2.4-1 and rap2.4-2 mutants demonstrated smaller lesions than the wild-type plants (Fig. 5C), confirming that in the observed disease phenotype. To measure disease development more precisely, fungal biomass accumulation of B. cinerea ActinA (Bc ActinA) per leaf   www.nature.com/scientificreports www.nature.com/scientificreports/ was measured at 3 and 6 dpi for each of the mutant lines. Compared to growth in wild-type plants, there was a clear decrease in the growth of fungal biomass in the rap2.4 mutants (Fig. 5D). This suggests that RAP2.4 gene contributes to plant immunity toward B. cinerea.

Down-regulation of RAP2.4 alters the expression of defense-regulated genes in response to
B. cinerea. In order to link RAP2.4 function in defense to specific pathway(s), we assessed the response of B. cinerea-infected tissues to molecular markers of different signaling pathways. The transcript levels of the SA-mediated defense-associated genes, PR1, β-1,3-glucanase (BGL2/PR2) and phytoalexin deficient 4 (PAD4) 34 were determined in rap2.4-1 mutant. The basal expression levels of all SA-associated genes in uninfected plants revealed significant reduction in rap2.4 compared to wild type ( Fig. 6A-C). Although the transcript of the SA marker gene, PR-1, increased at 18 hpi with B. cinerea in wild-type plants, the expression was significantly higher in rap2.4 (Fig. 6A). This suggests that rap2.4 plants respond faster in terms of PR-1 expression; thus, the increased PR-1 levels may be an indirect consequence of the decreased rate of fungal growth in rap2.4 plants. After B. cinerea inoculation, PR2 and PAD4 showed a similar pattern of expression between rap2.4 and wild type plants infected with B. cinerea although the repression was much lower in the mutant plants (Fig. 6B,C). This suggests that the basal expression level of SA-related defense genes might be dependent on RAP2.4.
At 18 hpi with B. cinerea, the transcript levels of JA/ET pathway defense marker genes, PDF1.2, PR3 (β-chitinase) and PR4 (hevein-like) were significantly induced in both wild-type and rap2.4 mutant plants ( Fig. 6D-F). The increase in the expression of PR3 and PR4, but not PDF1.2, was comparable in wild-type and rap2.4 plants. Hence, these B. cinerea-induced genes which are activated by the oxylipins, JA and/or OPDA, are COI1-dependent 35,36 and are negatively regulated by the TF, MYC2 37 . In wild-type plants, there was repression in VSP2 in response to B. cinerea; thus, the transcript levels of VSP2 were markedly increased in rap2.4 after inoculation (Fig. 6G). Similar to VSP2, healthy rap2.4 plants expressed lower basal levels of MYC2, but the transcript accumulated to high levels within 18 hpi (Fig. 6H). This is in contrast to the suppression of MYC2 in B. cinerea-inoculated wild-type plants. This suggests that the down-regulation of RAP2.4 may require oxylipins during pathogen infection.
There was an increase in transcript levels of 12-oxophytodienoate reductase 1 (OPR1) when wild-type and rap2.4 plants were challenged with the necrotrophic fungal pathogen B. cinerea (Fig. 6I). The induction was significantly increased in rap2.4 mutant. On the other hand, the expression of OPR3 was altered in rap2.4 mutant plants after infection (Fig. 6J). The expression of genes that were responsive to the cyclopentenones, phytoprostanes www.nature.com/scientificreports www.nature.com/scientificreports/ (PPs) or 12-oxo-phytodienoic acid (OPDA) 38,39 was also analyzed. The Arabidopsis detoxification-related genes, GST6 and GSTU19, were up-regulated in in wild-type-infected plants (Fig. 6K,L). No induction of the same genes was observed in rap2.4 mutant after infection, indicating that there is a common regulation between electrophilic oxylipins and B. cinerea, and that RAP2.4 plays a major role in this regulation. Overall, these results suggest that the down-regulation of RAP2.4 gene alters the expression of JA-and/or cyclopentenone-mediated response genes, which may result in the increased resistance of rap2.4 mutants to the necrotrophic pathogen B. cinerea.

Discussion
Plants are often affected by diseases mainly caused by pathogens 10,40 with a consequent serious reduction in plant growth, development and productivity. Abiotic stress conditions have a direct impact on plant-microbe interactions, and alternation of plant physiology and defense responses. Plants retain a range of defense mechanisms to combat stress conditions 41 , which involve a variety of metabolic reprogramming events and cellular pathways. Once sensing a single or multiple stress(es), plants show an immediate and evoking response to initiate a complex stress-specific signaling by synthesizing hormones and accumulating phenolic compounds 14 . In this work, we performed a meta-analysis of transcriptomic data related to environmental challenges in Arabidopsis publicly available at present to explore the complexity of the transcriptional changes of Arabidopsis. For that reason, we analyzed the transcriptome of Arabidopsis at one time point after pathogen infection, insect infestation or abiotic stress using Affymetrix ATH1 whole-genome GeneChips. We aimed to determine which genes, pathways, gene set categories and predicted PPI networks may play a key role in specific responses to environmental stresses (pathogen infections and abiotic cues). Our focus in the current study was not only on plant responses to the necrotrophic pathogen B. cinerea, but also to the virulent and avirulent bacterial strains of P. syringae pv. tomato, the other fungal pathogen A. brassicicola and the herbivore insect P. rapae. We also extended our analysis to include oxidative stress and wounding on Arabidopsis. Because alterations in the level of phytohormones have some prominent responses against environmental stresses in planta, we studied the effect of SA, MeJA, ACC and ABA treatments.
There is an increasing demand for more meta-analysis of transcriptomic studies 7,41 , which can be featured for many reasons. First, the transcript levels are highly affected by changing environmental conditions. Second, the inconsistency in the generated results, coming from field studies, can be affected by other disturbing factors. www.nature.com/scientificreports www.nature.com/scientificreports/ Third, the integration of high-throughput technologies (transcriptomics, proteomics and metabolomics) is of a concern nowadays due to the high costs and qualified experts in "omic" analyses. Moreover, the commonalities of independent studies will identify genes strongly associated with the studied stresses in order to focus on the functional analysis of these common DEGs 7,42 . Here, the well-defined analysis pipelines, standardized approaches for data analysis in addition to the availability of datasets to public were among the reasons that microarrays were selected over any other transcriptomic approaches.
In Arabidopsis, we found 6.8% and 5.3% of the transcriptome was considered up-and down-regulated in response to B. cinerea, respectively, which is in agreement with previous data in a number of studies 5,10,26 . The time point for global expression profiling was set to 18 hpi with B. cinerea, due to the fact that most changes in gene expression occur between 18 and 30 hpi 5,10,11,43 . In addition, the findings of our microarray meta-analyses on Arabidopsis plants individually treated with a tested single biotic stress, abiotic stress or hormone, are comparable to those in the literature 5,10,18,44,45 . It must be noted that in all Arabidopsis-stress combinations, many genes showed a more than twofold change in expression at the time point tested. Although several studies have reported meta-transcriptomic analysis combining different environmental stresses, this report linked B. cinerea-infected plants with other biotic and abiotic stresses at the transcriptional level. More than 50% of BUGs were also induced by any of the P. syringae pv. tomato strains. In addition, 40-50% of all consistent changes elicited by A. brassicicola and P. rapae were consistently triggered by B. cinerea. This suggests that these genes are commonly activated or repressed during these Arabidopsis-attacker interactions.
On the other hand, the number of co-regulated genes in response to B. cinerea and the tested Arabidopsis-abiotic challenge combinations was much lower. We also investigated the role of hormones in the regulation of the overlapping gene sets of Arabidopsis-B. cinerea interaction. We identified probesets representing individual hormone-responsive genes among the selected B. cinerea-responsive genes. Comparison of the hormone-responsive genes with that of B. cinerea-responsive probesets revealed that 34%, 33%, 29% and 36% of the BUGs are responsive to SA-, JA-, ACC-and ABA, respectively ( Table 1). The percentages of all hormonal responsive genes with the B. cinerea-repressed changes were even higher, indicating that hormones play a dominant role in the transcriptional reprogramming of Arabidopsis in response to B. cinerea infection. This is confirmed by another study of which the expression of BUGs is also affected by ethylene-insensitive 2 (ein2), coronatine insensitive 1 (coi1) and SA-deficient (nahG) mutations 5 , suggesting a regulatory role of hormones in mediating gene expression by B. cinerea which may have effects on disease responses. It has been reported that a number of P. syringae-, A. brassicicola-, and P. rapae-induced genes are also considered as JA-responsive genes 44 . Although pathogen/insect attackers with very different modes of action (e.g., B. cinerea, A. brassicicola, P. syringae and P. rapae) may induce similar sets of responsive genes, the majority of these genes are affected by a specific attacker. This suggests that hormones and other factors may enable fine-tuning the transcriptional machinery of defense response.
Interestingly, the eleven stress treatments induced only one gene and repressed three genes that could be considered as common regulators of the overlapping gene sets. The comparative microarray data analysis demonstrated that CORI3/JR2 was commonly up-regulated and RAP2.4, At1g72060 and At2g20670 were found to be commonly down-regulated (Fig. S3, Supporting Information). The CORI3/JR2 gene encoding cystine lyase is an enzyme that generates an ET precursor, and has been previously reported to be induced by B. cinerea, cold, drought and oxidative stresses 10 , P. syringae 40 , wounding 46 and MeJA 47 . Arabidopsis RAP2.4d, a common repressed gene to all stresses in test, is a member of RAP2.4 (AP2/DREB-type TF) family, was down-regulated in response to cold, light and ET but up-regulated in response to oxidative stress, increased salt and drought 48,49 , suggesting a role of RAP2.4 to coordinately regulate multiple abiotic stress responses. In addition to B. cinerea, the transcript levels of At1g72060 was reduced by cold, drought and oxidative stress 10,32 ; and At2g20670 was repressed by heat, salinity and osmotic stress 11 .
To get an insight into the function of commonalities of B. cinerea and other stress groups, we categorized the biological function according to the GO tool. Some of these functional categories covered a relatively large proportion of the Arabidopsis genome e.g. "electron transfer/energy pathways", representing most of the annotated genes commonly regulated by B. cinerea-biotic stresses and B. cinerea-hormonal stresses, but not B. cinerea-abiotic stress group. Thus, most of B. cinerea-abiotic stress group belonged to the functional category "biological processes" and "response to stress/stimulus". The number of up-regulated genes predicted to be involved in "response to biotic and abiotic stress" in all Arabidopsis biotic stress combinations tested was similar to observations done by the four biotic stresses of P. syringae, A. brassicicola, P. rapae and the feeding thrips (Frankliniella occidentalis) 44 . This distribution of the identified DEG sets over the various functional classes within the genome makes it possible to better understand the importance of the functional category in plant response. Analysis of GO data revealed that several B. cinerea-and abiotic stress-induced genes were in favor in "receptor binding" and "TF activity". Another large functional group of DEGs was enriched in the regulation of cellular components. These categories were consistent with the response of the plant to environmental stress known to derive to strong metabolic readjustment, defense mechanism and turning off the photosynthesis machinery. Our data are in accordance with previous findings in response to biotic and abiotic stresses in Arabidopsis 50,51 .
Modern 'omic' technologies can generate countless lists of biological identifiers/genes. Visualization of the identified overlapping genes/proteins can provide new insights into phenotypes and better understanding of the biological significance 52 . The co-expression data and PPI networks were incorporated and visualized to identify the regulation of the up-and down-regulated genes/proteins (Table 2). Our analysis of the interactive networks of proteins matches with what has been previously reported that biological networks may elucidate the roles of hub proteins in distinct pathway(s) involved in response to multiple environmental changes 53 . STRING is a distinctive database in a way that it not only allows the visualization of gene/protein lists of an interactome network, but also links it directly to a comprehensive annotation of the gene/protein lists; thereby providing biological implications of these overlapping as well as non-overlapping gene/protein sets 33  www.nature.com/scientificreports www.nature.com/scientificreports/ the interactome network can provide further biological information associated with the function of these genes in response to multiple environmental challenges.
The genes involved in co-expression analysis have been underlying in molecular network formation 54 . The co-expressed genes might be validated by their regulation is such having similar cis-regulatory elements for a TF. Thus, co-regulation studies validate the relationship of correlated genes. Among these is RAP2.4, which is involved in ET-activated signaling pathway 49 and was repressed by B. cinerea (Fig. 5). As most RAP2.4 TFs, RAP2.4d protein mediates the fine-control and adjusts the availability of three chloroplast peroxidases 48 . In addition, RAP2.4 can bind to both the dehydration-responsive element (DRE) and the ET-responsive GCC-box 49 . The T-DNA insertion line of RAP2.4 showed increased resistance to B. cinerea. The B. cinerea-induced RAP2.2 which also relies on ET signaling pathway, showed, on the other hand, decreased resistance to this fungal pathogen 55 . This indicates that RAP2.4 appears to be important in response to abiotic and biotic stresses, particularly in the tolerance to drought 49 and pathogenesis of B. cinerea. Although resistance mediated by SA has been reported to be mainly against biotrophic pathogens, other reports have shown the involvement of SA-mediated genes in response to necrotrophic pathogens including B. cinerea 14,22 . Similar to our results, PAD4 gene was repressed but PR1 gene was induced after B. cinerea infection 22 , regardless of the inoculated genotype (Fig. 6). In comparison to wild-type, rap2.4 mutant plants exhibited significantly decreased expression levels of PR1, PR2 and PAD4, indicating that the basal expression level of SA-responsive genes may depend on RAP2.4. In addition, JA/ET-mediated response contributes to plant defense against necrotrophic pathogens 14,22,56 . The JA/ET responsive genes, PDF1.2, PR3 and PR4, were induced by B. cinerea in wild-type and rap2.4 plants. The VSP2 gene that is regulated by the MYC2 can be activated by JA and wound but suppressed by ET 37 . Our results showed that the expression of VSP2 and MYC2 was enhanced in rap2.4 but repressed in wild type plants after B. cinerea infection. These findings are in agreement with a previous study that infections with B. cinerea induce Arabidopsis VSP2 and MYC2 at 14 hpi, followed by a significant decrease at 24 hpi 28 . This might be attributed the increased levels of the stress hormone ET leading to down-regulation of VSP2 and MYC2 in the wild-type-infected plants 57 . Interestingly, our qRT-PCR analysis showed that the reduced transcript levels of VSP2 (Fig. 6) seems to contradict the results obtained from the microarray data analysis ( Table 2). This discordance could be due to the different plant growth conditions and pathogen infection procedures in this study and the datasets obtained from publicly available microarrays.
Previous studies have reported that the cyclopentenone OPDA regulates gene expression in concert with JA to fine-tune the expression of defense genes 16,35 . Arabidopsis EXLA2 and OPR3 genes can modulate gene expression through COI1 or the electrophile effect of the cyclopentenones under biotic stress conditions 16,58,59 . In addition, Arabidopsis and tomato mutants of OPR3 enhanced resistance to B. cinerea 60,61 . Consistent with that, it was found that the transcript levels of OPR3 were dramatically reduced in rap2.4-infected plants, indicating a common regulation of gene expression in response to electrophilic oxylipins and B. cinerea. The detoxification-related genes, GST6 and GSTU19, that are known to be highly up-regulated by OPDA and PPs 38,39 were also induced in the wild-type, but not in the rap2.4 mutant, after infection. This confirms that the regulation of these genes is affected by the down-regulation of RAP2.4 toward B. cinerea. We speculate that rap2.4 mutant accumulates JA as well as cyclopentenone oxylipins upon infection with B. cinerea. To test this hypothesis, quantification of JA, OPDA and phytoprostane A 1 (PPA 1 ) in rap2.4 mutant before and after infection is among our priorities. Future research to dissect the importance of RAP2.4 in the regulation of JA-and/or cyclopentenone-mediated gene responses and to map all implicated elements in stress signal transduction pathways should be a major focus.
Our data provided a detailed picture of the highly dynamic interactive network involving several signaling proteins, in which inoculation with B. cinerea served as major control parameter to multiple stress responses in Arabidopsis. A recent genome-wide association study has revealed the combinatorial effect of P. rapae and drought on the regulation of genes in response to B. cinerea; thus by differentially regulating the level of resistance to B. cinerea 1 . The meta-analysis of the current study revealed that upon individual or multiple stress(es), several common genes play a significant role in the defense mechanism of the plant. We conclude that RAP2.4 has the potential to serve in plant defense through the regulation of endogenous signal molecules and/or pathogenderived effectors. Future investigations into the identification of pathogen-suppressed RAP2.4 gene expression and the relationship with membrane-associated microbe pattern (MAMP)-triggered defense will help explain the functions of RAP2.4 in innate immunity against B. cinerea. The long-term goal is to identify target genes, which can be helpful in breeding programs and agricultural biotechnology, and to generate crop varieties resistant to pathogen challenges, climatic changes and hormonal imbalances.
Methods plant growth and treatment conditions. We analyzed data from publicly available microarray datasets on Arabidopsis plants (ecotype Col-0) either infected or treated with the necrotrophic fungus pathogen B. cinerea, other biotic (P. syringae pv tomato DC3000 and avrRpm1, A. brassicicola and P. rapae), abiotic (oxidative stress and wounding) or hormonal (SA, MeJA, ACC or ABA) stresses. The microarray datasets were downloaded from NASCArrays at the BAR for Arabidopsis functional genomics database (http://bbc.botany.utoronto.ca/) 30 . The reference numbers for the corresponding stresses can be found in Table 3. Briefly, five-week-old Arabidopsis plants were inoculated by placing four 5-μl drops of 5 × 10 5 spore mL −1 solution of B. cinerea on detached leaves. Responses to B. cinerea infection were assayed at 0 and 18 hpi of adult leaves. For bacterial infections, P. syringae pv. tomato DC3000 or avrRpm1, were infiltrated with 1 × 10 8 CFU mL −1 . Bioassays using drop-inoculation of 3-µl drops of 1 × 10 6 spores mL −1 of A. brassicicola were carried out on detached leaves of 5-week-old plants. Infestations with P. rapae were performed by transferring 5 chewing larvae per plant 42 . Responses to these biotic stresses were tested at 0 and 24 hpi.
For the oxidative stress or wounding treatment, 18-day old seedlings were exposed to 10 µM methyl viologen (paraquat) or punctured with pins, respectively, at 0 and 24 hours post treatment (hpt) as previously described 8 . Arabidopsis seedlings grown in MS media were treated with 10 µM SA 8 , MeJA, ACC and ABA 62

. Expressions
Microarray procedure and normalization method. Affymetrix Expression Console (EC) software was used to treat the publicly available raw CEL files 63 . This included probeset signal integration, background correction and quantile normalization. Transcriptome Analysis Console (TAC) software was used to analyze DEGs.
Affymetrix GeneChip ATH1 genome arrays representing 22,810 Arabidopsis genes were considered in each dataset from which, those with the absence (A) and medium (M) signal detection calls were removed and those which showed present (P) were only considered for further analysis. Depending on the treatment, the number of tested samples (n) for each replicate varies (Table 3). Three technical replicates were used for the analysis for each gene and their average was taken as the signal intensity. The FC for each gene was calculated by dividing the signal intensity (SI) of that particular gene in the treated data by the SI of that gene in the non-treated 10,11 . The FC gene expression of ≥2 (up-regulated) or ≤0.5 (down-regulated) was set as default filter criteria for significant DEGs at P ≤ 0.05. Log 2 -transformed expression level data were used to generate scatter plots to detect the effect of each stress at 0 and the corresponding 3 hpt for hormones, 18 hpi for B. cinerea or 24 hpi/hpt for other biotic and abiotic stresses on plant gene expression. The gene identities and GO across microarray datasets were established using The Arabidopsis Information Resources (TAIR; www.arabidopsis.com). Although we compared gene sets across experiments and identified overlapping patterns of DEGs, the raw data was re-normalized under identical platforms regardless of spatiotemporal variability in the data.
The heatmap was generated using Morpheus (https://software.broadinstitute.org/morpheus/) to visualize large data matrices. To analyze the generated gene expression datasets, the top 50 BUGs and BDGs were plotted along with their corresponding DEGs of the other tested datasets. Log 4 values of FC were scaled as −4 to +4; and denoted by color intensity from red (down-regulated genes) to green (up-regulated genes). Black was denoted for the absence of gene expression in the specific dataset. To identify the common DEGs across the datasets, Venn diagrams were generated using Venny 2.1 software (http://bioinfogp.cnb.csic.es/tools/venny/index.html).
Arabidopsis ppi network. The PPI dataset was downloaded from A. thaliana Protein Interaction Network (AtPIN) 64 . The interactome data includes the experimentally identified PPIs and computationally predicted interactions. The PPI network was created and visualized using Cytoscape Version-3.7.0 65 . Node-edge attributes for the genes/proteins of interest were downloaded from STRING database (http://string-db.org/) 33 which contains the known and predicted PPIs of the query. The network was modified based on these genes/proteins using a Cytoscape plugin Style (formerly known as Vizmapper). In the network, nodes represent the proteins/edges; whereas lines represent the interactions/connections. The attributes of the nodes were altered for the visualization of the network in respect to the genes/proteins of interest. Another plugin, BINGO v3.0.3, was used to determine the GO categories, which were statistically overrepresented in the selected set of proteins 66 . fungal culture, plant inoculation and t-DnA insertion lines. For disease assays, B. cinerea strain BO5-10 was grown on V8 agar media (36% V8 juice, 0.2% CaCO 3 and 2% Agar). The fungal cultures were sub-cultured by transferring a piece of agar containing the mycelium to a fresh plate of V8 agar and incubated at 25 °C. In order to study the response of B. cinerea, we followed the procedure of disease assays coming from the microarray reference obtained from the BAR. Four rosette leaves from five-week-old soil grown Arabidopsis plants were drop-inoculated by placing a 3 µl drop of 2.5 × 10 5 conidiospores mL −1 solution of B. cinerea to each leaf 5 . Plants were further kept under a sealed transparent covered plate to maintain high humidity in a growth chamber with fluorescent lighting (186 μE m −2 sec −1 ) at 21 °C/18 ± 2 °C day/night temperature and a 16/8 hr light/dark cycle. Responses to B. cinerea infection were assayed at 0 and 18 hpi on detached leaves.
Homozygous T-DNA insertion lines of rap2.4-1 (SALK_139727; N677156) and rap2.4-1 (SALK_110897; N667030) were obtained from the NASC (Nottingham, UK). Mutant lines were in Col-0 ecotype background. T-DNA lines were confirmed via qRT-PCR using primers provided in Table S1 (Supporting Information). Fungal biomass was assessed by accumulation of Bc ActinA relative to Arabidopsis Actin2 (AtActin2). Lesion size was determined by measuring the diameter of the necrotic area.

RNA extraction and expression analysis.
Samples of B. cinerea-inoculated leaves were collected at 0 and 18 hpt for the qRT-PCR analysis. RNA extraction and qRT-PCR expression analyses were performed as described previously 16 . The qRT-PCR amplification was performed using Bio-Rad CFX96 Real-Time PCR System C1000 Thermal Cycler (Bio-Rad, Hercules, California, USA) in triplicates. Reaction mixture (20 μL) contains 100 ng total RNA as a template, 10 μL 2X GoTaq qPCR Master Mix, 0.4 μL 50X GoScript RT Mix for 1-Step qRT-PCR (Promega, Madison, Wisconsin, USA), 0.3 μM each specific left and right primer (Table S1, Supporting Information). The reaction condition was as follows: 40 °C for 15 min for the reverse transcription (RT) followed by 95 °C for 10 min of RT inactivation, 40 cycles of 95 °C for 10 s, 60 °C for 30 s, and 72 °C for 30 s. AtActin2 was used as an endogenous reference for normalization. Expression levels were calculated by the comparative cycle threshold method, and normalization to the control was performed.
Statistical analysis. For qRT-PCR assay, three technical replicates of each sample were used with a minimum of three biological replicates. Results were expressed as means ± standard deviation (SD) of the number of experiments. A Student's t-test for the values was performed at P < 0.05.
Data of B. cinerea growth and lesion size in inoculated plants represent the mean ± SD (n = 20). Analysis of variance (ANOVA) and Duncan's multiple range test were performed to determine the statistical significance 67 . Mean values followed by a different letter are significantly different from the corresponding control (P < 0.05). All experiments were carried out independently in triplicate with similar results.