Meta-analysis of transcriptomic responses to biotic and abiotic stress in tomato

Selection microarray studies

In this study, a meta-analysis of stress response was performed in tomato using publicly available microarray gene expression data implemented on a single platform (Affymetrix Tomato Genome Array, Santa Clara, CA, USA). Expression data of tomato plants exposed to BS and AS were combined. The raw expression data of various experiments were obtained from the Array Express from EBI (https://www.ebi.ac.uk/arrayexpress/).

Identification of differentially expressed genes

Microarray expression data from each source study was pre-processed separately as individual datasets. The raw data were normalized using Robust Multi-array Average background correction and quantile normalization. The data were processed in R version 3.2.2 (R Core Team, 2015). Meta-analysis of AS and BS was performed separately to identify DEGs involved in stress conditions. Genes that showed significant P-value (false discovery rate (FDR) of 5%), meaning that an estimated five percent of the DEGs are false positive, were considered as DEGs (Campain & Yang, 2010; Chang et al., 2013). Two methods of the combination of P-value have used: Fisher’s method and maxP (Wilkinson, 1951; Fisher, 1992). The Fisher’s combined probability test sums the logarithm-transformed P-values obtained from individual studies, under the null hypothesis, and follows a chi-squared distribution with 2k degrees of freedom, where k is the number of studies being combined (Rhodes et al., 2002). A very small P-value in just one study can be adequate to make statistical significance, even if the same genes are not significant in any other study. The maxP method takes maximum P-values across studies. It follows a beta distribution with parameters K (number of transcriptomic studies which are combined for a meta-analysis where each study contains G genes) and 1. This method targets on DEGs that have small P-values in all studies. (Chang et al., 2013; Song & Tseng, 2014).

Gene analysis

The genes selected by meta-analysis were further analyzed and characterized. In the first instance, the results were compared with the plant genes in the Plant Transcription factor & Protein Kinase Identifier and Classifier, iTAK (http://bioinfo.bti.cornell.edu/cgi-bin/itak/index.cgi) and the Plant Transcription Factor, PlantTFDB (http://plntfdb.bio.uni-potsdam.de/v3.0/) databases. The list includes 1,845 genes of tomato which encoding the TFs that directly or indirectly involved in signaling pathway and response to AS and BS. Venn’s diagrams (http://bioinfogp.cnb.csic.es/tools/venny/) were used to represent co-occurrence of DEGs. GO enrichment analysis of the DEGs (such as identifying biological processes) was conducted using the AgriGO platform (Du et al., 2010).

Network analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins, STRING 10.5 (http://www.mybiosoftware.com/string-9-0-search-tool-retrieval-interacting-genesproteins.html) (Szklarczyk et al., 2014) database was used to figure out all functional relations between the genes and their occurrence patterns across various genomes. The database makes relations based on several lines of evidences: the empirical evidence from protein–protein interaction assays, co-expression data from the NCBI Gene Expression Omnibus database, the extraction of information from other databases, coexistence of the genes in the same organisms, conserved gene neighborhood in known genomes, gene fusion events, pathway annotation from other resources such as the Kyoto Encyclopedia of Genes and Genomes or GO databases, and automated text-mining tools (Szklarczyk et al., 2014). STRING calculates a confidence value for those interactions according to the pieces of evidences from above 0.4 to 0.9, as the medium to highest score, respectively.

Identification of genes involved in biotic and abiotic stresses

After searching in the database, microarray data from 23 different experiments were selected and analyzed. To obtain a global analysis, 391 microarray samples including 213 AS from different categories: drought, heat, salinity, light, hormone, mineral and heavy metals deficiencies/toxicities and 178 BS that caused by viruses, nematode, fungi, bacteria and pests in tomato (wild type and transgenic) were chosen (Table S1). We found a total of 1,862 genes differentially regulated in response to BS, was far greater than that induced by AS, 835 genes, with a FDR ≤ 0.05 by Fisher’s statistical method (Tables S1–S3). We also figured out 15.4% or 361 DEGs with the conserved expression between AS and BS by Fisher’s statistical method, suggesting that these genes and their associated cell-signaling pathways are regulated in a similar way in a wide range of stresses (Fig. 1A; Table S5). Using maxP method, a total of 14 genes differentially regulated in response to AS and a total of 220 genes differentially regulated in response to BS were found. We figured out 0.9% or 2 DEGs (LES.1399.1 and LES.3969.1) with the conserved expression between AS and BS by maxP method (Fig. 1A; Tables S2–S4). These methods used in other studies to identify important genes in different physiological processes in other plants. Shaar-Moshe, Hübner & Peleg (2015) identified 225 DEGs across-species by using 17 microarray experiments of drought stress at the reproductive stage of Arabidopsis, rice, wheat and barley. Fisher’s and maxP methods are two popular microarray meta-analysis methods used in the literature that have different strength for detecting different types of DEGs (Tseng, Ghosh & Feingold, 2012; Wang et al., 2012). Our results also showed that the two used methods detected different sets of DEGs, suggesting different assumptions behind the methods. Among these methods, Fisher’s method is more popular in the microarray meta-analysis studies (Wang et al., 2012); hence, we were more interested in detecting DEGs across all studies through Fisher’s statistical method.

Gene ontology

We focused on the 361 DEGs common to both AS and BS and subjected them to GO and network analysis to explore other possible functions of them. GO enrichment analysis of the DEGs showed that shared DEGs represented genes involved in main biological and cellular processes. Among them, the greatest percentage of genes involved in the cellular process (>76% genes), metabolic process (>76% genes) and response to stimulus (>50%). One of largest functional group of DEGs was associated with metabolic processes (e.g., carbohydrate metabolic process, lipid metabolic process, protein metabolic process, nucleobase, nucleoside, nucleotide and nucleic acid metabolic process, cellular amino acid and derivative metabolic process and regulation of primary metabolic process).

In addition, another large functional group of DEGs was associated with cellular processes (e.g., cellular metabolic process, regulation of cellular process and cellular response to stimulus), the show a significant rearrangement in plant metabolism as part of stress adaptation (Shaar-Moshe, Hübner & Peleg, 2015).

On the other hand, the lowest percentage of genes involved in positive regulation of biological process (5%) (Fig. 2; Table S6). The most significant GO terms for cellular process were cellular metabolic process (FDR: 1.8E-37), organic acid metabolic process (FDR: 9.2E-29), cellular ketone metabolic process (FDR: 1.2E-28) in metabolic process were primary metabolic process (FDR: 7.8E-39), cellular metabolic process (FDR: 1.8E-37) and organic acid metabolic process (FDR: 9.2E-29). Moreover, in response to stimulus were response to abiotic stimulus (FDR: 3.4E-40), response to biotic stimulus (FDR: 9.4E-26), and response to stress (FDR: 9.4E-26). Some of these genes possess important roles in DEGs shared by AS and BS networks (Table S7). A genome-wide transcriptional analysis was performed just for Fe deficiency and was identified 97 DEGs by comparing Fe-deficient and Fe-sufficient in roots of tomato. These transcripts are related to the physiological responses to the nutrient stress resulting in an improved iron uptake, including regulatory aspects, translocation, root morphological modification, and adaptation in primary metabolic pathways, such as glycolysis and TCA cycle (Zamboni et al., 2012). Another research project showed that the largest functional groups of DEGs was associated respectively to metabolic processes (e.g., amino acid and carbohydrate metabolism), regulator function (e.g., protein degradation and transcription) and response to the stimulus that was identified by the research of Shaar-Moshe, Hübner & Peleg (2015).

Network analysis

Network analysis using STRING 10.5 identified the connections among DEGs common to AS and BS (Fig. 3; Table S7). In addition, some genes that possess an important role in the three biological pathways; cellular process, metabolic process, response to stimulus were detected. For example, Solyc02g090890.2.1 (ZE: plays an important role in resistance to stresses, seed development and dormancy), Solyc04g049350.2.1 (CS1: biosynthesis of aromatic amino acids), Solyc07g052480.2.1 (LOC544276: involved in storage lipid mobilization during the growth of higher plant seedling), Solyc07g049530.2.1 (ACO1: involved in ethylene biosynthesis via S-adenosyl-L-methionine), Solyc01g099630.2.1 (XTH1: an essential component of the primary cell wall, and thereby participates in cell wall construction of growing tissues), Solyc05g055230.1.1 (RPS17: ribosomal small subunit assembly), Solyc12g098350.1.1 (infA: no specific function has so far been attributed to this initiation factor; however, it seems to stimulate more or less all the activities of the other two initiation factors, IF-2 and IF-3), Solyc02g082920.2.1 (CHI3: defense against chitin containing fungal pathogens), Solyc06g082600.2.1 (UBC3: mediates the selective degradation of abnormal and short-lived proteins), Solyc01g096480.2.1 (LOC543657: probable thiol-disulfide oxidoreductase that may play a role in proper chloroplast development and negatively regulates Cf-9/Avr9-mediated cell death and defense responses), Solyc07g047850.2.1 (CAB4: it receives and transfers excitation energy to photosystems).

The results of network analysis also showed some of the TFs have connections to other molecules. Based on pieces of literatures reviewing, detected TFs in this study, Jasmonate Ethylene Response Factor 1 (JERF1), EIL3, EIL2, SlGRAS6, Protein LATE ELONGATED HYPOCOTYL (LHY) and MYB 48 affected by a wide range of stresses and possibly play essential and significant roles in multiple stress responses which is in agreement with their functional roles in BS and AS responses. Expression of JERF1 was induced by salt, ethylene, MeJA and abscisic acid (ABA) treatments in tomato that may be pointed to a central role for JERF1 in different signal transduction pathways (Zhang et al., 2004). It has been also reported that JERF1 activates ABA biosynthesis related genes such as NtSDR (short-chain dehydrogenase/reductase) and enhanced cold and salt stress tolerance in transgenic tobacco (Wu et al., 2007; Ouyang et al., 2007). It showed that the expression of EIL2 and EIL3 was induced by ethylene and salinity treatments (Ying et al., 2004) and SlGRAS6 is one of the most important members of the tomato GRAS family in disease resistance and mechanical stress (Mayrose et al., 2006). Different studies have shown the large family of MYB has essential roles in plant growth, development, primary and secondary metabolism and response to BS and AS (Oh, Park & Han, 2003; Van der Ent et al., 2008; Dou et al., 2016). It was also reported that LHY controls the expression of a large number of ABA-responsive genes by binding directly to their promoters (Huang et al., 2007). These results support that the detected TFs in our study might play an important role in BS and AS responses.

Identification of transcription factors involved in biotic and abiotic stresses

In this research, because of the importance of TFs as a powerful tool for the manipulation of complex metabolic pathways, we have tried to identify genes encoding them. According to the PlantTFDB, 1,845 genes from tomato encoding TFs were identified and classified into 58 families that directly or indirectly involved in signaling and response to AS and BS (Table S8). Based on our results, from the 1,862 genes response to BS selected by Fisher’s statistical method, 93 genes (5%) encode TFs (Fig. 1C; Table S9). Ethylene-responsive transcription factors (ERF) family has the greatest number of genes (18 genes), and C3H, FAR1, GATA, LBD, MIKC-MADS, PHD, SLP protein, TALE, TIFE and YABBY families have the lowest number of genes (each family has only one gene) (Table 1). From the 835 genes response to AS selected by Fisher’s statistical method, 31 of those (3.7%) encode TFs (Fig. 1B; Table S10). ERF and NAC families have the greatest number of genes (five genes), and ARF, BBR-BPC, GRAS, NF-NY and VOZ families have the lowest number of genes (each family has only one gene) (Table 1).

Since TFs are the key regulators of plant growth, development and metabolisms via up- and down-regulation of various genes, identification of the genes encoding these proteins and transferring them into plants to improve their stress tolerance would be interesting. These TFs (JERF1, Pti6, SlERF1, SlERF4, SlGRAS6, MYB48, protein LHY, EIL2, EIL3, WRKY 11, WRKY 26, inducer of CBF expression 1 (ICE1)-like and bZIP TF) can be used as candidate genes in molecular breeding programs because these TFs are common to AS and BS; also, our results indicate they are being affected by a wide range of stresses. In addition, they possibly play essential and significant roles in multiple stress responses, which is in agreement with their functional roles in BS and AS responses (Gu et al., 2002; Ying et al., 2004; Zhang et al., 2004, 2010; Borrás-Hidalgo et al., 2006; Mayrose et al., 2006; Wu et al., 2007; Yáñez et al., 2009; Sharma et al., 2010; Chen et al., 2012; Kim, Stork & Mudgett, 2013; Feng et al., 2013; Li, Luan & Liu, 2015). The results of network analysis also showed some of these proteins have central roles in BS and AS responses (Fig. 3). The two-way Venn diagram indicated that the number of common TFs between BS and AS was 13 while 18 and 80 TFs were uniquely identified in AS and BS, respectively (Fig. 4).

ERF family

The ERF has been shown to play a key role during plant development and stress tolerance mechanisms. Overexpression of ERF family genes in different plant species such as Arabidopsis, tobacco, rice and tomato has been shown to increase tolerance to a wide range of BS and AS (Sharma et al., 2010). The tomato has 112 AP2/ERF superfamily genes (Sharma et al., 2010, 2013). Our results showed that four of 13 TFs shared between AS and BS belonged to ERF family, including JERF1, SlERF4, SlERF1 and pti6 (Table 2).

Jasmonate Ethylene Response Factor 1, has critical roles in plant abiotic stress responses, such as salinity, low temperature, dehydration and various BS in different plant species (Zhang et al., 2004, 2010; Wu et al., 2007). It has been shown that the expression of JERF1 in tomato was induced by salt, ethylene, MeJA and ABA treatments that may be pointed to a central role for JERF1 in different signal transduction pathways (Zhang et al., 2004). Overexpression of JERF1 enhanced the tolerance of transgenic tobacco and rice plants to high salinity, low temperature, osmotic stress and pathogen attack by regulating the expression of downstream stress-responsive genes and a number of pathogenesis-related (PR) genes due to binding to DRE/CRT and GCC-box cis-elements (Wu et al., 2007; Zhang et al., 2010). Expressing JERF1 in tobacco induced expression of GCC box-containing genes such as CHN50, Prb-1b, GLA, NtSDR and osmotin that subsequently resulted in enhanced tolerance to salt stress and low temperature (Zhang et al., 2004; Wu et al., 2007). JERF1 up-regulated the expression of genes related to osmolyte synthesis such as OsP5CS and OsSPDS2, which, has increased the production of osmolytes. It also induced the expression of OsABA2 and Os03g0810800 genes (containing DRE element in their promoter regions), which encode two key enzymes involved in ABA biosynthesis and increased the synthesis of ABA in transgenic rice that resulted in elevated tolerance to drought (Zhang et al., 2010). JERF1 may also be located downstream of EIN3 in the ethylene signaling pathway (Zhang et al., 2004). It has been previously reported that SlERF4, an ERF in tomato, is a XopD, a type III secretion effector from Xanthomonas campestris, destabilize SlERF4 and this task has been proposed to promote bacterial growth (Kim, Stork & Mudgett, 2013). SlERF4 reduces ethylene production and, in turn, limits disease symptom development (Kim, Stork & Mudgett, 2013) (Table 2). Pti6 belonged to the ERF family in the plant that binds specifically to the GCC-box cis-element present in the promoter regions of many PR genes and induced the expression of a wide range of these genes that have key roles in the plant defense. It was shown that the expression of jasmonic acid- and ethylene-regulated genes, such as PR3, PR4, PDF1.2 and Thi2.1, was being differently affected by Pti4 or other similar genes, Pti5, or Pti6 (Gu et al., 2002; Chakravarthy et al., 2003) (Table 2). The gene encoding SlERF1 was another detected gene that probably has an important role in response to various stresses. It has been shown that the overexpression of SlERF1 in Arabidopsis can confer resistance to necrotrophic fungi like Plectosphaerella cucumerina and Botrytis cinerea (Mayrose et al., 2006) (Table 2).

GRAS family

SlGRAS6 was another identified TF that is one of the most important members of tomato GRAS family in disease resistance and mechanical stress (Mayrose et al., 2006). SlGRAS6 is required for tomato disease resistance to the bacterial pathogen Pseudomonas syringae pv. tomato. It has been shown that the suppression of SlGRAS6 gene expression through virus-induced gene silencing reduced resistance of tomato to Pseudomonas syringae pv. tomato (Mayrose et al., 2006) (Table 2).

MYB family and MYB-related family

The MYB is a large TF family in plants that have essential roles in plant growth, development, primary and secondary metabolism and response to BS and AS (Oh, Park & Han, 2003; Van der Ent et al., 2008; Dou et al., 2016). Most of MYB genes involved in response to various AS, belong to the R2R3-type group. With the comparison of OsMYB48-1 expression with other MYB proteins, a special role for this gene in response to the AS in rice can be suggested. It has been shown that its expression was strongly induced by ABA, H₂O₂, PEG, and dehydration stresses, while slightly induced by salt and cold stresses (Xiong et al., 2014). It was reported that the overexpression of OsMYB48-1 regulated the expression level of some stress responsive-genes, such as RAB21, RAB16D, RAB16C and LEA3, under drought stress conditions (Xiong et al., 2014) (Table 2). The LHY TF plays a critical role in the regulation of low temperature stress response via controlling expression of CBF 1, 2 and 3, COR27 and COL1 genes and contribute to plant defense through the circadian regulation of stomatal aperture (Grundy, Stoker & Carré, 2015; Seo & Mas, 2015). It was also shown that LHY controls the expression of a large number of ABA-responsive genes by binding directly to their promoters (Huang et al., 2007).

EIL family

It has been shown that the expression of EIL2 and EIL3, which encode TFs in ethylene signal transduction, was induced by ethylene and salinity treatments (Ying et al., 2004). It also reported that its ortholog in tobacco is required for resistance against Peronospora hyoscyami (Borrás-Hidalgo et al., 2006) (Table 2).

WRKY family

Different studies have shown that the members of WRKY family contribute multiple biological processes such as plant growth and development. They also have important roles in a wide range of BS and AS such as drought, salt, invasion of pathogens, implying that the members of this family may be putative regulators in response to various BS and AS (Huang et al., 2012; Li, Luan & Jin, 2012). Two of 13 TFs, WRKY11 and WRKY26, which are involved in both AS and BS, belong to the WRKY family (Table 2). In other plant species, it was also shown that the members of this family have a key role in AS and BS responses (Li et al., 2011; Chen et al., 2012). It was reported that OsWRKY11 involved in drought and heat tolerance in rice (Wu et al., 2009). It was shown that, in Arabidopsis thaliana, WRKY26 positively regulates the cooperation between heat shock proteins and ethylene-activated signaling pathways that mediate responses to heat stress (Li et al., 2011). It has also been shown that the promoters of some heat tolerance-related genes including Hsp and Hsf genes such as HsfA2, HsfB1, Hsp101, and MBF1c contain W-box sequences that are recognized by WRKY proteins (Li et al., 2011).

ICE1-like bHLH family

Inducer of CBF expression 1-like encodes a MYC-type basic helix-loop-helix (bHLH) TF (Chinnusamy et al., 2003). ICE1 is an important inducer of CBF3/DREB1A which regulates cold stress tolerance (Chinnusamy et al., 2003; Huang et al., 2015; Feng et al., 2013) (Table 2). In Arabidopsis, ice1 mutants showed decreased cold and chilling tolerance, whereas ICE1 overexpression in Arabidopsis increased the cold tolerance of transgenic plants (Chinnusamy et al., 2003). Furthermore, osmotic and salt tolerance were also much markedly increased in transgenic tobacco carrying SlICE1 from tomato by increasing the expression level of DREB1/CBF and their target genes (Feng et al., 2013) (Table 2).

bZIP family

Transcription factors of the bZIP family regulate various biological processes in plant growth and development, including morphogenesis and seed formation, as well as in AS and BS responses by binding to promoters of specific target genes to up- or down-regulate their expression (Li et al., 2015). WRKY TFs, bZIP domain and MYB TFs involved in plant defense against pathogens (Alves et al., 2013) (Table 2). TF bZIP11 is a member of the S1 class of the bZIP family, and it was shown that bZIP11 mRNA translation is repressed in response to sucrose (Hummel et al., 2009). It has been shown that sucrose repression of bZIP11 mRNA translation depends on the presence of the 5′-leader sequence upstream of bZIP11 gene (Hummel et al., 2009). It has also been shown that sucrose signaling has an important role in the translational control of bZIP11 in response to stress conditions (Hummel et al., 2009).

These reports implied that these TFs are important regulators of biological processes and may be useful in the breeding of plants to improve their tolerance against adverse stresses. The presence or absence of the TFs across various plant species and other organisms including bacteria, archaea and other eukaryotes is shown in Fig. 5. A high homology was observed for all TFs except protein LHY between S. lycopersicum and S. tuberosum and the highest percentage of similarity was found between Vitis vinifera, Populus trichocarpa and Glycine max. The results showed that GRAS6 is the only TF present in bacteria with a low similarity to that of tomato GRAS6. Just three of the identified TFs, ICE1-like, protein LHY and MYB48, are present in Opisthokonta that the homology between MYB48 and those of Opisthokonta was the highest. MYB48 is the only protein, which is present in all organisms except bacteria and archaea (Fig. 5). It has been shown that MYB is one of the largest and most diverse TF families, which are the key factors in regulatory networks controlling metabolism, development and responses to BS and AS (Roy, 2016). The members of MYB family found in all eukaryotes and molecular studies in Arabidopsis has been shown the role of this family in the epigenetic regulation of stress responses in plants (Roy, 2016). The results showed MYB48 as one of the members of MYB family possesses an important role in the network of abiotic and biotic responses (Fig. 3). In addition to MYB48, three other TFs, WRKY11, WRKY26 and protein LHY, were observed in Amoebozoa with low similarity to those of tomato (Fig. 5). As well as Amoebozoa, the WRKYs were observed in Giardia lamblia, also known as Giardia intestinalis (Fig. 5). It was previously thought that WRKYs are specific to plants and algae (Eulgem & Somssich, 2007) but recent studies have shown their presence in human protozoan parasite G. lamblia and Dictyostelium discoideum (slime mold) (Wu et al., 2005). It has been reported that WRKYs also have important roles in other physiological processes. It was shown that the expression of genes encoding proteins of cell wall increased, when a WRKY-like gene overexpressed in G. lamblia (Pan et al., 2009).

The results showed that protein LHY present in all eukaryotes except G. lamblia, Trichomonas vaginalis and Trypanosomatidea (Fig. 5) and also suggested a central role for this protein in the network of abiotic and biotic responses (Fig. 3). The results of other studies also support an important role for this protein and suggested that this protein involved in circadian clock control (Lu et al., 2009).