Gene 555 (2015) 127–139

Contents lists available at ScienceDirect

Gene journal homepage: www.elsevier.com/locate/gene

Synergistic regulatory networks mediated by microRNAs and transcription factors under drought, heat and salt stresses in Oryza Sativa spp. Deepti Nigam, Sanjeev Kumar ⁎, D.C. Mishra, Anil Rai, Shuchi Smita, Arijit Saha Centre for Agricultural Bio-Informatics (CAB-in), Indian Agricultural Statistics Research Institute (IASRI), Library Avenue, PUSA Campus, New Delhi, India

a r t i c l e

i n f o

Article history: Received 20 March 2014 Received in revised form 12 September 2014 Accepted 26 October 2014 Available online 31 October 2014 Keywords: miRNA Transcription factors Gene regulatory network Algorithm for the reconstruction of accurate cellular networks

a b s t r a c t Background: Transcription factors (TFs) and microRNAs (miRNAs) are primary gene regulators within the cell. Regulatory mechanisms of these two main regulators are of great interest to biologists and may provide insights into the abiotic and biotic stresses. However, the interaction between miRNAs and TFs in a gene regulatory network (GRN) still remains uncovered. Previous research has been mostly directed at inferring either miRNA or TF regulatory networks from data. However, networks involving a single type of regulator may not fully reveal the complex gene regulatory mechanisms, therefore study of interplay among these two regulators in gene regulation is important towards explaining the mechanism of different abiotic stresses. Result: Oligonucleotide microarrays containing 51,279 transcripts were used to identify total 133 salt responsive target genes regulated by 11 TFs that are also differentially regulated by miRNA under salinity, heat and drought stresses in Oryza sativa. TF's-target interactions which are most enriched in their downstream regulation were also identified. Many genes whose encoded proteins are implicated in response to light and radiation stimulus, hormone stimuli, oxidative stress, copper ion binding and electron transport were found to be enriched. However the majority were novel for the combined abiotic stress, which indicates that there are a great number of genes induced after the exposure these abiotic stresses and regulated by miRNA. Conclusion: Analysis of the expression profile data of Oryza provides clues regarding some putative cellular and molecular processes that are undertaken in response to these stresses. The study also identified a large number of candidate functional genes that appear to be constitutively involved in salt, drought and heat stresses tolerance. Further examination of these genes may enable the molecular basis of abiotic stress tolerance in Oryza, to be elucidated. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Abiotic stress causes substantial loss to agricultural production globally. Rice (Oryza sativa) is a staple food crop that is cultivated worldwide and constitutes a primary source of food especially in Asia. Abiotic stresses basically drought, extreme temperatures and salinity are the major limiting factor in rice cultivation. Reflecting the diminishing environmental conditions, more often than not crops are exposed simultaneously to multiple stresses resulting in vast changes in the molecular scene within the cell. Therefore, comprehensive understanding of the regulatory networks that govern the dynamic adaptive changes in a

Abbreviations: GRN, Gene regulatory network; ARACNe, Algorithm for the reconstruction of accurate cellular networks; CLR, Context likelihood ratio; RN, Relevance network; TFs, Transcription factors; miRNA, microRNA; CAN, Comparative network analysis; TFBSs, Transcription factor binding sites; CRMs, Cis-regulatory modules. ⁎ Corresponding author. E-mail address: [email protected] (S. Kumar).

http://dx.doi.org/10.1016/j.gene.2014.10.054 0378-1119/© 2014 Elsevier B.V. All rights reserved.

plant responding to stress is critical to meet future energy needs. Besides high agricultural importance, rice is model plant organisms representing monocots and has extensive biological resources including complete genome sequence, highest number of microarray and transcriptome studies. Therefore it is important to study the molecular mechanism of different kinds of abiotic stress within rice as it may reveal a number of pivotal attributes that may behave similarly for other crops of the same family. MicroRNAs (miRNAs) are approximately 21-nucleotide (nt) noncoding RNAs that play critical roles in gene expression regulation at the post-transcriptional level. In plants, cleavage of the target mRNA appears to be the prevalent method of posttranscriptional regulation. Brodersen et al. (2008) provided evidence that plant miRNA-guided silencing has a widespread translational inhibitory component. Plant miRNA-guided gene regulation has been shown to be involved in multiple developmental processes including organ polarity (Bowman, 2004), leaf growth (Chuck et al., 2007), and male or female sterility (Millar and Gubler, 2005). Also it has been observed that some plant's

128

D. Nigam et al. / Gene 555 (2015) 127–139

miRNA respond to abiotic stress conditions and some miRNA targets are stress-related genes suggests that miRNAs may play important roles in the plant stress response (Phillips et al., 2007; Schommer et al., 2012) Several stress-specific miRNAs have also been identified in model plants under various biotic and abiotic stress conditions, including high salinity (Sunkar et al., 2008), drought (Zhao et al., 2007), cold (Lv et al., 2010), oxidative stress (Sunkar et al., 2006), UV-B radiation (Zhou et al., 2007) and mechanical stress (Lu et al., 2005). Advancements in high throughput technologies have resulted in generation of huge omics data addressing different aspects of temporal and spatial response in variety of stresses in plants. Among them microarray technology has revolutionized the identification of global transcriptomic changes. Multiple transcriptomic studies have been conducted on same or related stress conditions. Thus meta-analysis of related microarray studies is increasingly becoming popular to validate conclusions for hypothesis under consideration (Tseng et al., 2012). However, limited meta-analysis studies are available in plant systems (Ghanekar et al., 2008; Meier et al., 2008; Shaik and Ramakrishna, 2013). Meta-analysis of microarray data of Arabidopsis related to infection by eight different viruses revealed hub genes that are highly connected, organized in modules and are central to plant defence response (Rodrigo et al., 2012). It is reported that in plants which are subjected to multiple stresses, there exists extensive cross-talk between different stress responses via hormonal signalling pathways (Kreps et al., 2002; Seki et al., 2002). Thus, it is imperative to correlate and analyse different kinds of stress responses to find the genes, proteins and metabolites that are common as well as specific to different kinds of abiotic and biotic stress conditions (Atkinson et al., 2013). Meta-analysis of microarray studies involving samples from a wide range of tissues, developmental stages and different levels of stresses but specific to one stress condition would unravel interesting features related to the stress response (Altmann et al., 2004; Schmid et al., 2005). Also comparative analysis of such interesting molecular profiles from different tissue and stresses would allow the identification of unique and shared features. Thus, in this study, large scale comparative analysis was performed via meta-analysis of microarray data on correlated abiotic stresses mainly salt, heat and drought within rice. In order to elucidate the crosstalk between these three stresses conditions, knowledge of the expression status of miRNAs involved in these stress responses is critical. Therefore, this analysis is focused on miRNA genes and their target TFs which are common in all three stresses. Synergetic regulatory network has been constructed taking above abiotic stress induced miRNAs, their target TFs and genes. The constructed network has been validated using various in-silico approaches. Extensive analysis of various gene sets based on Gene ontologies (GO) and metabolic pathway was done to unravel the underlying biological mechanisms related to different stresses under consideration. The performance of co-expression network of miRNA regulated TF's-targets genes which divides the stress responsive genes as co-expressed modules was studied to reveal organization of “synergistic stress transcriptome” for rice.

microarray, usually provided in .cel file format, were log transformed and median cantered per array. Probe sets with expression mean μ b 50 and S.D. i.e. σ b 0.3 μ were considered non-informative and therefore excluded, leaving 35,254 probe sets for the analysis. 2.2. Identification of stress responsive miRNA and their target TFs In the present study miRNAs and their differentially expressed target genes were filtered from the salt, heat and drought stresses from PMTED database (http://pmted.agrinome.org/). The fold change ≥2.0 was selected during analysis. Constitutively expressed (i.e. commonly expressed) miRNA and target genes were identified by Venn analysis for these three stresses. We have used RiceSRTFDB (http://www.nipgr. res.in/RiceSRTFDB.html) for the identification of stress responsive transcription factors among miRNA target genes. 2.3. Network inference using ARACNE Network generation for TFs and their regulated genes done through ARACNe (Algorithm for the Reconstruction of Accurate Cellular Networks) (http://wiki.c2b2.columbia.edu/califanolab/index.php/ Software/ARACNE). This tool was developed and successfully used to infer a GRN in human B cells (Basso et al., 2005; Margolin et al., 2006). Recently ARACNe algorithm has been used for interactome analysis for different stages of fibre development and drought stress condition in Gossypium hirsutum especially for AUX/IAA and brassinosteroid TFs (Nigam and Sawant, 2013; Nigam et al., 2014). ARACNe calculates expression correlations between genes based on mutual information and identify statistically significant correlate genes. ARACNe provides direct interactions (1st neighbours) and 2nd neighbours for selected hub genes. We have used two other methods for our network validation like CLR (Context Likelihood Ratio) and RN (Relevance Network). 2.4. Validation of synergetic regulatory network

In this study first we collected and processed gene expression profile data. The stress responsive miRNA and their target TFs were identified. The network for TFs and their regulated genes were generated. The TFs of the network were further annotated. Also, their meta-analysis was performed for the summarization of expression levels. Also their motif analysis was performed. Finally, different techniques were used for the validation of constructed network. Detail methodology is given in following sub-sections.

The constructed network using ARACNe algorithm was further authenticated by two different methods based on mutual information i.e. CLR (Context Likelihood Ratio) (Roy et al., 2014) and RN (Relevance Network) (Nazri and Lio, 2012). The performance of these three mutual information based approaches was compared on the basis of number of nodes and edges obtained. For further validation of constructed network, a TF's knock-down approach was used that mimics like in-vitro knockdown experiment (Olsen et al., 2014). For this validation, 10 TFs (Nuclear transcription factor Y subunit A-10, ADO3/flavin-binding, AGL16, OsFBDUF57, OsNAC5, bZIP, OsMYB, ARF2, ARF8, SPL2) out of 11 TFs were used. Each of these TFs was knocked down one by one systematically at the time of modelling network by removing the expression values of the concerned TF. The construction of the perturbed network was performed using reduced adjacency matrix resulted in reduced network for each of the knocked TF with less number of targets. Given a perturbed constructed network and a list of target genes significantly affected by a specific knock-down of TF concerned, the descendants of the knocked down TF were classified in the perturbed network. In this case, each of these descendants was named as true positive (TP), as false positive (FP) if it was not affected by the perturbation, and finally a gene was affected by the knock-down experiment but it was not inferred as a descendant in the network, it was named as false negative (FN). This classification then allowed one to compute a quality measure such as the F-score for the perturbed network (Olsen et al., 2014). (Tables 2a and 2b, Supplementary file 3)

2.1. Gene expression profile data collection and processing

F ¼ 2  TP=ð2  TP þ FP þ FNÞ;

Total 216 microarray data sets were downloaded from GEO database (http://www.ncbi.nlm.nih.gov/geo/). Data sets of oligonucleotide

Where, F = 1 corresponds to the perfect classification of the affected genes and F = 0 to no correctly identified affected genes.

2. Material and methods

F∈½0; 1 :

ð1Þ

D. Nigam et al. / Gene 555 (2015) 127–139

129

Table 1 The quality measurement of the constructed network for 11 TFs used in the study. Transcription factors

Total number of the targets under knock-down condition

Total number of the targets (control condition)

TP

FP

FN

F-score

ADO3/flavin-binding AGL16|AGAMOUS-like 16 OsFBDUF57 OsNAC5 bZIP transcription factor family OsMYB ARF2 ARF8 SPL2 WRKY family transcription factor

129 104 120 115 102 107 121 110 113 128

133 133 133 133 133 133 133 133 133 133

129 104 120 115 102 107 121 110 113 128

100 36 49 95 87 59 73 64 44 61

29 60 32 38 46 74 58 32 58 58

0.66 0.68 0.75 0.63 0.61 0.62 0.65 0.70 0.69 0.76

2.5. Annotation analyses of rice gene chip Constitutive TF's target genes were analysed using r functional categorization based on three GO categories at p-values ≥ 0.05. The agriGO tool (http://bioinfo.cau.edu.cn/agriGO/) was used to perform the enrichment analysis using SEA (Singular Enrichment Analysis) coupled with available background data of rice probes. Gene percentage analysis was calculated for each agriGO annotation in the GO category. Rice gene chip annotation was based on the top hits using the PLEXdb tool. 2.6. Meta-profiling GENEVESTIGATOR was used for meta-profiling (i.e. summarizing expression levels according to the biological contexts of the samples) of target genes. Here the purpose of meta-profiling is to address the hypothesis that a selected set of differential expression signatures i.e. miRNA regulated TFs shares a significant intersection with their target genes (a meta-signature) within rice abiotic stress expression datasets, thus inferring a biological relatedness between miRNA-TFs and TF'stargets. For this purpose the identified targets were cross-checked for their occurrence in miRNA regulated differentially expressed target genes deposited in PMTED database these genes were further used for meta-profiling. 2.7. Identification of over-represented motifs in promoters of target genes Transcription factor's target genes were selected for motif analysis. One kb upstream regions from these genes were extracted from rice genome database (TIGR). MEME (http://meme.sdsc.edu/meme/website/ meme.html), a publicly available motif discovery program (Bailey and Elkan, 1994) was used for motif analysis. MEME was invoked to identify motifs 6 to 50 bp long that look highly conserved among promoter sequences of the selected genes. 3. Results 3.1. Pipeline of the analytic procedure There may be some biological relationship between the unique set of miRNA and targeted transcripts which has been explored in this study.

The schematic workflow of the experiment is given in Fig. 1. A model regulatory network was designed for the miRNA regulated TFs reported by different studies (Gutierrez et al. 2009; Wang et al., 2011; Atkinson and Urwin, 2012). Firstly miRNA and their targets were extracted from PMTED database. Annotation of target genes was carried out using PLEXdb annotation tool. RiceSRTFDB was used for the identification of TFs. Ideally the identification of putative target of these TFs requires expression profile dataset which should be ≥ 100 in number. ARACNe algorithm was applied at this step for the generation of network model. Using this network model, we extracted downstream target genes. Further we have in-silico validated our data with the GENEVESTIGATOR and thereafter we have tried to figure out the enriched cis-regulatory element using MEME and TOMTOM. Therefore, our study is more suitable to evaluate whether a regulatory model can successfully identify key downstream regulators purely based on expression profiles without depending on external knowledge. 3.2. Construction of synergistic gene regulatory networks and identification of downstream regulators for salt, heat and drought stresses In order to identify the putative target genes regulated by TF in harmony with miRNA, 35,254 probe sets were filtered out from 216 expression profiles taken from GEO database for salt, heat and drought stresses. The filtered miRNAs and their differentially expressed target genes under salt, heat and drought stresses, obtained from PMTED database, were examined. It was found that only 4 miRNAs were unique in salt stress whereas 10 miRNAs were unique for heat and drought stresses separately. It was also observed that 26 miRNAs were commonly expressed in all three stresses whereas 4, 5 and 14 miRNAs were common in heat and salt, salt and drought and drought and heat respectively (Fig. 2a). Likewise 54 transcripts were unique in heat stress, 16 in salt stresses and 47 were unique in drought stress whereas 24 transcripts (targets) are common among all these stresses (Fig. 2b). It was hypothesized here that there may be some relationship in between these unique and common sections of miRNA and their targets (Fig. 2c). Therefore, the relationship between common miRNAs and their targets were further explored here. RiceSRTFDB annotation revealed that there were 11 transcription factors (Os.22707.1.S1_at, Os.52157.1.S1_x_at, Os.49762.1.S1_at, Os.53273.1.S1_at, Os.9492.1.S1_a_at, OsAffx.20062.1.S1_at, OsAffx. 21972.1.S1_at, Os.11696.2.S1_x_at, Os.7177.2.S1_a_at, OsAffx.2442.1.S1_

Table 2a Comparative Network Analysis (CAN) using ARACNe, CLR and RN methods. Methods

No. of nodes

No. of edges

ARACNe (AR) CLR RN Common

133 110 180 No. of nodes = AR (105) CLR, AR (112) RN

1200 1605 16,287 No. of edges = AR(1013) CLR, AR (832) RN

Table 2b Percentage similarity between the ARACNe and other two methods. Methods Comparison

Percentage

AR–CLR AR–RN

(105/133) 78.94 (112/133) 84.21

130

D. Nigam et al. / Gene 555 (2015) 127–139

Fig. 1. The procedure for screening and analysing the target genes of TFs regulated by miRNA.

x_at, Os.10472.1.S1_a_at) among 24 transcripts (targets). These 11 TFs were further checked in the miRNA-target data set and were found to be regulated by total 5 miRNAs namely osa-miR528, osa-miR408, osamiR156, osa-miR164 and osa-miRi159. These 5 miRNAs were also found to be represented in the commonly expressing 26 miRNAs (Fig. 2a) in all three stresses thus conferring their relationship with each other (Fig. 2c). Therefore, these 11 transcription factors were further used as a “Hub genes” and mutual information based network learning algorithm ARACNe was used for gene regulatory network inference. A combined stringent cut-off of an error tolerance ε = 0.2 and a p-value threshold of mutual information (MI) at 0.05 were used in inference of global gene networks. It was observed that 10 TFs out of 11 TF genes resulted in high confidence 133 downstream targets (Fig. 3, Supplementary file 1), whereas the TF (Os.10472.1.S1_a_at) has not shown any connectivity with other target genes in the developed network. During analysis the bootstrapping was applied to generate 100 bootstrapped networks and consensus network was obtained. 3.3. Validation of the constructed network The performance of three mutual information based approaches i.e. ARACNe algorithm, CLR and RN was compared on the basis of number of nodes and edges obtained. Network obtained from ARACNe algorithm consists of 133 nodes and 1200 edges (Supplementary file 2), whereas, networks form CLR and RN based method are having 110 and 180 nodes and 1605 and 1628 edges respectively. Further, it was found that ARACNe (AR) with CLR shares 105 overlapping nodes and 1031 edges indicating 78.94% similarity whereas, with RN it was 84.21% as they share 112 nodes and 732 edges (Table 2a and 2b). In the case of systematic validation based on a knock-down approach, it was found that the F score for different networks resulted as a consequence of the knock down of each TF ranges between 0.61

and 0.76 and on an average each perturbed network fit with score of 0.70 (Supplementary file 3). This indicates that the constructed network correctly identifies the effected target genes governed by each of the 10 TFs. 3.4. Singular Enrichment Analysis (SEA) for the identification of enriched GO terms in miRNA-TF regulated target genes Gene ontology analysis was carried out using the Singular Enrichment Analysis (SEA) tool offered by agriGO (Du et al., 2010) at default settings of Fisher t-test (p b 0.05), False Discovery Rate (FDR) correction by Hochberg method and five minimum numbers of mapping entries against species specific pre-computed background reference (Fig. 4). SEA method revealed a number of common metabolism in Oryza spp. under the salt, heat, and drought conditions (Supplementary file 4; Fig. 4). We observed a cohort of significant GO processes like response to light stimulus, response to radiation, response to hormone stimulus, oxidoreductase activity, and metal ion binding along with cellular nitrogen metabolic process that are well reported for their role in individual as well as combined abiotic stresses (Xue et al., 2004; Takashi & Kazuo, 2010; Ahmad & Prasad, 2012; Kumar et al., 2013). 3.5. Meta-profiling of target genes The target genes which overlap with the PMTED database were filtered with their fold expression value and further used for signature analysis using GENEVESTIGATOR. For meta-profiling analysis of only top 50 target genes (based on their fold change expression, FC ≥ 2.0) were used (Figs. 5, 6 and 7). Meta-profiling was performed separately for all three stress conditions (i.e. salt, heat and drought stresses). Interestingly for all three cases similar type of profiling was observed in related type of conditions i.e. strong differential expression of target genes was observed in conditions related to stress like drought, salt and

D. Nigam et al. / Gene 555 (2015) 127–139

a)

131

b)

c)

Fig. 2. a) Venn analysis for identification of unique and commonly regulated miRNAs during heat, drought and salt stresses. b) Venn analysis for identification unique and commonly regulated miRNA targets during heat, drought and salt stresses. c) Relationship between unique and common miRNAs and their targets during heat, salt and drought stresses.

cold. On comparison of separate meta-profiling of target genes under each of three conditions reveals that green patches are predominant under heat stress as compared to other two stresses, where a balanced

ratio of down-regulated (green patches) to up-regulated (red patches) genes was observed. The bias towards the down regulation of the genes indicates that miRNA may have higher impact on gene regulation

Transcription Factors TF’s Targets

Fig. 3. Transcriptional regulatory network miRNA-TF's target genes under combined stress (i.e. salinity, drought and heat stresses). A total number of 133 transcripts participated in the regulatory program, clustered in 9 non-overlapping modules (main frame). Commonly expressed miRNA regulated TFs and their targets are colour-coded as indicated in the left-down inset. The directionality of the regulation is indicated through arrowed edges heading from the TF to the target gene.

132

D. Nigam et al. / Gene 555 (2015) 127–139

Fig. 4. Singular Enrichment Analysis of miRNA target genes for the identification of significant GO (gene ontology) terms i.e. BP (biological process), MF (molecular function), CC (cellular component). Enriched GO terms revealed their involvement in process related to abiotic stress.

under heat stress as compared to other two stresses (Fig. 7). The relative similarity value for each experiment was found to be N 0.5 which is nothing but the correlation value of one experiment its nearest neighbour experiment. 3.6. Identification of stress responsive genes under combined stress Numbers of unreported target genes were identified in Oryza under combined abiotic stress but their functionally annotated homologue genes exist in Arabidopsis (Supplementary file 2). A number of interesting genes were found such as GTP-binding family protein, HSP90.5, MADSbox transcription factor 20, alcohol dehydrogenase family protein, sorbitol dehydrogenase, O-methyltransferase 1, NADH-ubiquinone oxidoreductase 20 kDa subunit, NAC domain containing protein 36, SAUR-like auxin-responsive protein family, IAA14, RCI3A (peroxidase superfamily protein), RAB GTPASE HOMOLOG B18, and ATM3 (ABC transporter of the mitochondrion 3). Recently it was reported that the overexpression of stress-inducible small GTP-binding protein AhRab7 (AhRabG3f) in peanut (Arachis hypogaea L.) enhances abiotic stress tolerance (Lin et al., 2012) whereas the constitutive overexpression of a stress-inducible small GTP-binding protein PgRab7 from Pennisetum glaucum enhances abiotic stress tolerance in transgenic tobacco (Agarwal et al., 2008). Similarly, MADS box family members are also reported to play an important role in different abiotic stress responses. As in the case of rice, OsMADS26, an AGL12 class gene, was also reported to be involved in stress tolerances. Recently, flowering related MADS-box genes, such as rice SHORT VEGETATIVE PHASE (SVP)-like genes OsMADS22 and OsMADS55 and Arabidopsis SUPPRESSOR OF CONSTANS 1 (SOC1) are found to be involved in stress tolerances, revealing novel roles for MADS-box genes. Similarly, Imai and Yahara 2000, have reported that in yeast HSP90 provides the salt stress tolerance via stabilization and regulation of calcineurin. Moreover, IAA14 is a member of Auxin family transcription factors that are also known to be involved in abiotic stress resistance (Song et al., 2012). 3.7. Enrichment analysis of combined stress inducible cis-regulatory elements in putative target genes Fundamentally, the regulation of gene expression is mediated by the binding of TFs to specific cis-regulatory elements present in the

promoter sequence of their target genes. One or more TF(s) can regulate the expression of other TF(s) by binding to their cis-regulatory elements, which result in complex regulatory networks for highly controlled and specific pattern of gene expression. A few such transcriptional regulatory networks operative in response to abiotic stresses have been defined in Arabidopsis and rice (Nakashima et al., 2009). Since, our interest lies here in the mechanism of TF's interplay in combined stress therefore, the identification of commonly elicited cis-regulatory elements in the targets becomes imperative. Thus to reveal the presence of various cis-regulatory elements in the promoter, 1 kb upstream region of each target genes was analysed by MEME suite (http://meme.nbcr.net/meme/) search. Based on the available knowledge in the literature, the identified cis-regulatory elements were categorized as stress-responsive and other motifs using TOMTOM. Interestingly it was observed that these 5 motifs were clustered (co-occurrence) together giving a clue of their similar as well as synchronous role in response to stress taken under consideration. The significance of co-occurrence of these five motifs on each target can be observed by the combined p-value (Fig. 8). Gene regulation under synergistic effect of combined abiotic stress can result from the combinatorial interplay of transcription factors (TFs) with cis-regulatory elements (i.e. cis-regulatory module), which are guided by clusters of TF binding sites (TFBSs) (Arnone and Davidson, 1997; Ghazi and VijayRaghavan, 2000). However, the presence of TFBSs along with their structural arrangements within cis-regulatory modules (CRMs) of 100–10,000 bp is critical for comprehensive gene expression regulation. In present study a number of stress-responsive elements were identified in the promoters of miRNA regulated TF's target genes. Their probability of occurrence at that site is represented in terms of significant E-value (Fig. 9a, b and c). Motif 1 was found to be enriched with E-value of 1.3e− 063 and hit with three well know TFs i.e. abi4, bZIP911 and Gamyb. Similarly Motif 2 shows the E-value of 4.9e − 031 and it hits with total eight TF's binding sites namely abi4, bZIP911, Gamyb, HMG-1, bZIP910, id1, EmBP-1 and Dof3. Likewise Motif 3 is overrepresented with E-value of 1.2e− 031 with a total of six hits to abi4, bZIP911, Gamyb, HMG-1, bZIP910 and EmBP. Similarly, Motif 4 gives hits with only four TF's binding sites i.e. abi4, bZIP911, bZIP910 and id1 with E-value of 1.8e−009. Lastly, in the case of Motif 5 we observed the highest number of hits i.e. nine hits with abi4, bZIP911, bZIP910, Gamyb, HMG-I/Y, id1, Gamyb, PEND and Dof. Among these, abi4, bZIP911 were common

D. Nigam et al. / Gene 555 (2015) 127–139 Fig. 5. Signature analysis for the identification of expression pattern of miRNA regulated TF's target genes during drought stress. Red dotted circles represent signature gene expression in PMTED and blue dotted circle represents their contrast expression in GENEVESTIGATOR database. It shows the contrast gene expression pattern for drought stress induced target genes regulated by miRNA-TFs in different experiments. The relative similarity between different experiments is presented in the secondary axis.

133

134 D. Nigam et al. / Gene 555 (2015) 127–139 Fig. 6. Signature analysis for identification of expression pattern of miRNA regulated TF's target genes during heat stress. Signature genes expressions in PMTED are majorly up regulated whereas these genes are down regulated in GENEVESTIGATOR database. It shows the contrast gene expression pattern for heat stress induced target genes regulated by miRNA-TFs in different experiments. The relative similarity between different experiments is presented in the secondary axis.

D. Nigam et al. / Gene 555 (2015) 127–139 Fig. 7. Signature analysis for identification of expression pattern of miRNA regulated TF's target genes during salt stress. Red dotted circles represent signature gene expression in PMTED and blue dotted circle represents their contrast expression in GENEVESTIGATOR database. It shows the contrast gene expression pattern for salt stress induced target genes regulated by miRNA-TFs in different experiments. The relative similarity between different experiments is presented in the secondary axis.

135

136

D. Nigam et al. / Gene 555 (2015) 127–139

Fig. 8. Co-occurrence of the identified top five motifs in miRNA regulated TF's target genes that commonly (constitutive) expressed in all three stresses (salt, heat and drought stresses). Nearby occurrence can be due to the instant requirement for “Hub TF's” that are directly guided by miRNA in response to any of the three stresses.

among all five enriched motifs. Interestingly the majority of the identified motifs are also found in other studies where the coordination of different TFs is required (Leister et al., 2011). ABI4 is a member of the AP2/ERF family and members of this family bind specifically to the ABRE elements in the promoters of abiotic stress-responsive genes and regulate their expression (Wind et al., 2013). Further in Arabidopsis, group A bZIP-type proteins such as AREB/ABF proteins are a major factor in ABA-responsive gene expression under osmotic stress conditions (Yamaguchi-Shinozaki and Shinozaki, 2006). Several recent studies have indicated that the miR159–Gamyb pathway might also be involved in the abiotic stress response (Xin et al., 2010; Zhang et al., 2013). The integration of expression data and cis-regulatory elements along with experimental data will provide insights into the regulatory interactions among miRNA, TFs and other proteins regulating stress responses in rice. 4. Discussion Microarray technology has enabled us to decipher the genes regulated either by miRNA or directly by TFs. There is very less information known about the mechanism of how miRNA regulated TF's target genes employ to effect co-ordinated regulation within the multiple stress effected cell. We propose an in-silico analysis pipeline to find transcriptional modules for regulatory gene expression signatures like miRNA and TFs. Our method only requires miRNAtarget information and expression profiles in the appropriate context such as tissue type or stress condition. This procedure was applied to identify key downstream targets for the 11 prognostic signature genes (TFs) for abiotic stress. According to the previous studies these transcription factors are miRNA regulated and found to be crucial for tolerance against different abiotic stresses (Chinnusamy et al., 2004; Jain et al., 2007; Fang et al., 2008; Ding et al., 2009; Hossain et al., 2010; Takasaki et al., 2010; Leyva-González et al., 2012; Lima et al., 2012; Todaka et al., 2012). In this regard, a global transcriptional network was constructed using the ARACNe algorithm that resulted in the identification of 133

targets. Validation of these targets from two other network inference methods (CLR and RN) has shown higher percentage of similarity between their identified targets, suggesting accuracy among the used methods as well as of our modelled network. Similarly the TF's knockdown approach of validation has shown network accuracy of ~ 70%. These downstream targets include known regulators for abiotic stress inducible effects such as osmolyte accumulation, senescence and stomata development. The utility of our method may be expandable to other types of signatures such as the combination of correlated biotic stress. Signature analysis of top 50 target genes shows contrast expression pattern when compared with GENEVESTIGATOR tool. It was observed that in all three cases i.e. drought, heat and salt, similar type of experimental conditions was observed and their relative similarity value was found to be ≥ 0.5 (Figs. 5–7). In case of drought as well as in salt stress, affected target genes regulated by miRNA-TF were grouped into four and five clusters respectively. It indicates that these target genes behave in similar fashion across experimental conditions. However, in the case of heat stress no visible cluster was observed. Moreover the core binding sites that seem to be associated with transcription factor families, there may be a degree of subtlety in the cis-regulatory elements important for transcription factors to facilitate their unique regulatory effects in response to diverse adverse conditions. Here in the present study, different combinations of TFs and their cis-regulatory elements were identified, which may function through the recruitment of their target in different combinations during combined stress (i.e. salt, drought and heat stresses). Recently it was found that ABI4 binds the specific element (CCACGT) at the promoter of their target genes to prevent their transcription by competing with other competitive TFs that leads to long-lasting signalling. Therefore mutant of ABI4, showed weaker drought-tolerance after a herbicide norflurazon treatment. Plant bZIP transcription factors have been shown to homodimerize or heterodimerize. Wiese et al. (2004) have shown that in snapdragon (Antirrhinum majus) bZIP910 and bZIP911 bind to hybrid c-box/g-box ACGT elements as homodimers, whereas, heterodimers of these transcription factors show low affinity for these

D. Nigam et al. / Gene 555 (2015) 127–139

137

(a)

(b)

(c) Fig. 9. (a–c) Motifs identified in miRNA-TF's target genes using MEME and TOMTOM tools. Five consensus motifs are displayed as motif logos, and information about each motif is shown in left hand side. These motifs were found to be correlated with different abiotic stress responses.

boxes (Martinez-Garcia et al., 1998). Thus, such unique system is involved in controlling the activity of bZIP transcription factors, including transcription, translation, post-translational modifications, and homo-

dimerization or hetero-dimerization with other transcription factors. This type of arrangement provides a flexible, multi-responsive regulation of bZIP target genes.

138

D. Nigam et al. / Gene 555 (2015) 127–139

Fig. 10. Hypothetical model showing miRNAs, TFs and their target relationship in response to drought, heat and salt stresses governing resistance in a broad spectrum manner.

In conclusion during the coexistence of multiple stresses some combinatorial effectors work together and may provoke the major regulators such as miRNAs and TFs. These regulators in turn regulate their targets with some refined binding sites and stimulate the downstream processing and thus may provide broad spectrum abiotic resistance (Fig. 10). We have to explore the relationship between the target genes and their regulators (TF and miRNA) in vivo so that actual mechanism by which they are regulated in response to the combined abiotic stress can be elucidated. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.gene.2014.10.054.

Accession codes The gene expression profiles used in this study were series.

Funding sources We financial support from, the National Fund for Basic and Strategic Research in Agricultural Sciences, Indian Council for Agricultural Research (ICAR), PUSA, New Delhi (NFBSFARA/ASP-4009/2013-14).

Conflict of interest The authors declare that there is no conflict of interest.

Acknowledgement We acknowledge Indian Agricultural Statistics Research Institute, PUSA, and New Delhi for providing all facility to conduct this study.

References Agarwal, P.K., Agarwal, P., Jain, P., Jha, B., Reddy, M.K., Sopory, S.K., 2008. Constitutive overexpression of a stress-inducible small GTP-binding protein PgRab7 from Pennisetum glaucum enhances abiotic stress tolerance in transgenic tobacco. Plant Cell Rep. 27, 105–115. Ahmad, P., Prasad, M.N.V. (Eds.), 2012. Abiotic stress responses in plants: metabolism, productivity and sustainability. Springer Science + Business Media, LLC http://dx. doi.org/10.1007/978-1-4614-0634-1_8. Altmann, T., Weigel, D., Nover, L., 2004. AtGenExpress—Ein multinational koordiniertes Programm zur Erforschung des Arabidopsis Transkriptoms. GenomXpress 3, 13–14. Arnone, M.I., Davidson, E.H., 1997. The hardwiring of development: organization and function of genomic regulatory systems. Development 124, 1851–1864. Atkinson, N.J., Urwin, P.E., 2012. The interaction of plant biotic and abiotic stresses: from genes to the field. J. Exp. Bot. 63 (10), 3523–3543. Atkinson, N.J., Lilley, C.J., Urwin, P.E., 2013. Identification of genes involved in the response of Arabidopsis to simultaneous biotic and abiotic stresses. Plant Physiol. 162, 2028–2041. Bailey, T.L., Elkan, C., 1994. Fitting a mixture model by expectation maximization to discover motifs in bipolymers. Department of Computer Science and Engineering, University of California, San Diego. Basso, K., Margolin, A.A., Stolovitzky, G., Klein, U., Dalla-Favera, R., Califano, A., 2005. Reverse enginering of regulatory networks in human B cels. Nat. Genet. 37 (4), 382–390. Bowman, J.L., 2004. Class III HD‐Zip gene regulation, the golden fleece of ARGONAUTE activity? Bioessays 26, 938–942. Brodersen, P., Sakvarelidze-Achard, L., Bruun-Rasmussen, M., Dunoyer, P., Yamamoto, Y.Y., Sieburth, L., Voinnet, O., 2008. Widespread translational inhibition by plant miRNAs and siRNAs. Science 320 (5880), 1185–1190. Chinnusamy, V., Schumaker, K., Zhu, J.K., 2004. Molecular genetic perspectives on cross‐talk and specificity in abiotic stress signalling in plants. J. Exp. Bot. 55, 225–236. Chuck, G., Meeley, R., Irish, E., Sakai, H., Hake, S., 2007. The maize tasselseed4 microRNA controls sex determination and meristem cell fate by targeting Tasselseed6/indeterminate spikelet1. Nat. Genet. 39, 1517–1521. Ding, D., Zhang, L., Wang, H., Liu, Z., Zhang, Z., Zheng, Y., 2009. Differential expression of miRNAs in response to salt stress in maize roots. Ann. Bot. 103, 29–38. Du, Z., Zhou, X., Ling, Y., Zhang, Z., Su, Z., 2010. agriGO: a GO analysis toolkit for the agricultural community. Nucleic Acids Res. 38, W64–W70. Fang, Y., You, J., Xie, K., Xie, W., Xiong, L., 2008. Systematic sequence analysis and identification of tissue-specific or stress-responsive genes of NAC transcription factor family in rice. Mol. Gen. Genomics 280, 547–563.

D. Nigam et al. / Gene 555 (2015) 127–139 Ghanekar, R., Srinivasasainagendra, V., Page, G.P., 2008. Cross-chip probe matching tool: a web-based tool for linking microarray probes within and across plant species. Int. J. Plant Genomics 2008, 7 (Article ID 451327). Ghazi, A., VijayRaghavan, K., 2000. Developmental biology: control by combinatorial codes. Nature 408, 419–420. Gutierrez, L., Bussell, J.D., Păcurar, D.I., Schwambach, J., Păcurar, M., Bellini, C., 2009. Phenotypic plasticity of adventitious rooting in arabidopsis Is controlled by complex regulation of AUXIN RESPONSE FACTOR transcripts and microRNA abundance. Plant Cell. 21 (10), 3119–3132. Hossain, M.A., Lee, Y., Cho, J.-I., Ahn, C.-H., Lee, S.-K., Jeon, J.-S., Kang, H., Lee, C.-H., An, G., Park, P.B., 2010. The bZIP transcription factor OsABF1 is an ABA responsive element binding factor that enhances abiotic stress signaling in rice. Plant Mol. Biol. 72, 557–566. Imai, J., Yahara, I., 2000. Role of HSP90 in salt stress tolerance via stabilization and regulation of calcineurin. Mol. Cell Biol. 20 (24), 9262–9670. Jain, M., Nijhawan, A., Arora, R., Agarwal, P., Ray, S., Sharma, P., Kapoor, S., Tyagi, A.K., Khurana, J.P., 2007. F-box proteins in rice. Genome-wide analysis, classification, temporal and spatial gene expression during panicle and seed development, and regulation by light and abiotic stress. Plant Physiol. 143, 1467–1483. Kreps, J.A., Wu, Y., Chang, H.-S., Zhu, T., Wang, X., Harper, J.F., 2002. Transcriptome changes for Arabidopsis in response to salt, osmotic, and cold stress. Plant Physiol. 130, 2129–2141. Kumar, K., Kumar, M., Kim, S.-R., Ryu, H., Cho, Y.-G., 2013. Insights into genomics of salt stress response in rice. Rice 6 (27). Leister, D., Wang, X., Haberer, G., Mayer, K.F., Kleine, T., 2011. Intracompartmental and intercompartmental transcriptional networks coordinate the expression of genes for organellar functions. Plant Physiol. 157, 386–404. Leyva-González, M.A., Ibarra-Laclette, E., Cruz-Ramírez, A., Herrera-Estrella, L., 2012. Functional and transcriptome analysis reveals an acclimatization strategy for abiotic stress tolerance mediated by Arabidopsis NF-YA family members. PLoS One 7, e48138. Lima, J.C.d., Loss-Morais, G., Margis, R., 2012. MicroRNAs play critical roles during plant development and in response to abiotic stresses. Genet. Mol. Biol. 35, 1069–1077. Lin, S., Li, R., Xiang, X., Wang, J., Qiao, L., Song, X., Pei, Y., Sui, J., 2012. Overexpression of stress-inducible small GTP-binding protein AhRab7 (AhRabG3f) in peanut (Arachis hypogaea L.) enhances abiotic stress tolerance. J. Food Agric. Environ 10 (3&4), 888–894. Lu, S., Sun, Y.-H., Shi, R., Clark, C., Li, L., Chiang, V.L., 2005. Novel and mechanical stress– responsive microRNAs in Populus trichocarpa that are absent from Arabidopsis. Plant Cell Online 17, 2186–2203. Lv, D.-K., Bai, X., Li, Y., Ding, X.-D., Ge, Y., Cai, H., Ji, W., Wu, N., Zhu, Y.-M., 2010. Profiling of cold-stress-responsive miRNAs in rice by microarrays. Gene 459, 39–47. Margolin, A.A., Nemenman, I., Basso, K., Wiggins, C., Stolovitzky, G., Califano, A., 2006. ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian celular context. BMC Bioinf. 7 (Supl1), S7. Martinez-Garcia, J.F., Moyano, E., Alcocer, M.J., Martin, C., 1998. Two bZIP proteins from Antirrhinum flowers preferentially bind a hybrid C-box/G-box motif and help to define a new sub-family of bZIP transcription factors. Plant J. 13, 489–505. Meier, S., Bastian, R., Donaldson, L., Murray, S., Bajic, V., Gehring, C., 2008. Co-expression and promoter content analyses assign a role in biotic and abiotic stress responses to plant natriuretic peptides. BMC Plant Biol. 8, 24. Millar, A.A., Gubler, F., 2005. The Arabidopsis GAMYB-like genes, MYB33 and MYB65, are microRNA-regulated genes that redundantly facilitate anther development. Plant Cell Online 17, 705–721. Nakashima, K., Ito, Y., Yamaguchi-Shinozaki, K., 2009. Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 149, 88–95. Nazri, A., Lio, P., 2012. Investigating meta-approaches for reconstructing gene networks in a mammalian cellular context. PLoS One 7, e28713. Nigam, D., Sawant, S.V., 2013. Identification and analyses of AUX-IAA target genes controlling multiple pathways in developing fiber cells of Gossypium hirsutum L. Bioinformation 9, 996. Nigam, D., Kavita, P., Tripathi, R.K., Ranjan, A., Goel, R., Asif, M., Shukla, A., Singh, G., Rana, D., Sawant, S.V., 2014. Transcriptome dynamics during fibre development in contrasting genotypes of Gossypium hirsutum. Plant Biotechnol J. 12 (2), 204–218. Olsen, C., Fleming, K., Prendergast, N., Rubio, R., Emmert-Streib, F., Bontempi, G., HaibeKains, B., Quackenbush, J., 2014. Inference and validation of predictive gene networks from biomedical literature and gene expression data. Genomics 103, 329–336.

139

Phillips, J.R., Dalmay, T., Bartels, D., 2007. The role of small RNAs in abiotic stress. FEBS Lett. 581, 3592–3597. Rodrigo, G., Carrera, J., Ruiz-Ferrer, V., del Toro, F.J., Llave, C., Voinnet, O., Elena, S.F., 2012. A meta-analysis reveals the commonalities and differences in Arabidopsis thaliana response to different viral pathogens. PLoS One 7, e40526. Roy, S., Bhattacharyya, D.K., Kalita, J.K., 2014. Reconstruction of gene co-expression network from microarray data using local expression patterns. BMC Bioinformatics 15, 1–14. Schmid, M., Davison, T.S., Henz, S.R., Pape, U.J., Demar, M., Vingron, M., Schölkopf, B., Weigel, D., Lohmann, J.U., 2005. A gene expression map of Arabidopsis thaliana development. Nat. Genet. 37, 501–506. Schommer, C., Bresso, E.G., Spinelli, S.V., Palatnik, J.F., 2012. Role of microRNA miR319 in plant development. MicroRNAs in Plant Development and Stress Responses. Springer, pp. 29–47. Seki, M., Narusaka, M., Ishida, J., Nanjo, T., Fujita, M., Oono, Y., Kamiya, A., Nakajima, M., Enju, A., Sakurai, T., 2002. Monitoring the expression profiles of 7000 Arabidopsis genes under drought, cold and high‐salinity stresses using a full‐length cDNA microarray. Plant J. 31, 279–292. Shaik, R., Ramakrishna, W., 2013. Genes and co-expression modules common to drought and bacterial stress responses in Arabidopsis and rice. PLoS One 8, e77261. Song, Y., Wang, Z., Bo, W., Ren, Y., Zhang, Z., Zhang, D., 2012. Transcriptional profiling by cDNA-AFLP analysis showed differential transcript abundance in response to water stress in Populus hopeiensis. BMC Genomics 13, 286. Sunkar, R., Kapoor, A., Zhu, J.-K., 2006. Posttranscriptional induction of two Cu/Zn superoxide dismutase genes in Arabidopsis is mediated by downregulation of miR398 and important for oxidative stress tolerance. Plant Cell Online 18, 2051–2065. Sunkar, R., Zhou, X., Zheng, Y., Zhang, W., Zhu, J.-K., 2008. Identification of novel and candidate miRNAs in rice by high throughput sequencing. BMC Plant Biol. 8, 25. Takashi, H., Kazuo, S., 2010. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 61, 1041–1052. Takasaki, H., Maruyama, K., Kidokoro, S., Ito, Y., Fujita, Y., Shinozaki, K., YamaguchiShinozaki, K., Nakashima, K., 2010. The abiotic stress-responsive NAC-type transcription factor OsNAC5 regulates stress-inducible genes and stress tolerance in rice. Mol. Gen. Genomics. 284, 173–183. Todaka, D., Nakashima, K., Shinozaki, K., Yamaguchi-Shinozaki, K., 2012. Toward understanding transcriptional regulatory networks in abiotic stress responses and tolerance in rice. Rice 5, 1–9. Tseng, G.C., Ghosh, D., Feingold, E., 2012. Comprehensive literature review and statistical considerations for microarray meta-analysis. Nucleic Acids Res. 40, 3785–3799. Wang, Z., Lu, Y., Yang, B., 2011. MicroRNAs and atrial fibrillation: new fundamentals. Cardiovasc. Res. 89, 710–721. Wiese, A., Elzinga, N., Wobbes, B., Smeekens, S., 2004. A conserved upstream open reading frame mediates sucrose-induced repression of translation. Plant Cell 16 (7), 1717–1729. Wind, J.J., Peviani, A., Snel, B., Hanson, J., Smeekens, S.C., 2013. ABI4: versatile activator and repressor. Trends Plant Sci. 18 (3), 125–132. Xin, M., Wang, Y., Yao, Y., Xie, C., Peng, H., Ni, Z., Sun, Q., 2010. Diverse set of microRNAs are responsive to powdery mildew infection and heat stress in wheat (Triticum aestivum L.). BMC Plant Biol. 10, 123. Xue, Z.-Y., Zhi, D.-Y., Xue, G.-P., Zhang, H., Zhao, Y.-X., Xia, G.-M., 2004. Enhanced salt tolerance of transgenic wheat (Tritivum aestivum L.) expressing a vacuolar Na+/H+ antiporter gene with improved grain yields in saline soils in the field and a reduced level of leaf Na+. Plant. Sci 167 (4), 849–859. Yamaguchi-Shinozaki, K., Shinozaki, K., 2006. Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu. Rev. Plant Biol. 57, 781–803. Zhang, Y., Zhao, L., Li, H., Gao, Y., Li, Y., Wu, X., Teng, W., Han, Y., Zhao, X., Li, W., 2013. GmGBP1, a homolog of human ski interacting protein in soybean, regulates flowering and stress tolerance in Arabidopsis. BMC Plant Biol. 13, 21. Zhao, B., Liang, R., Ge, L., Li, W., Xiao, H., Lin, H., Ruan, K., Jin, Y., 2007. Identification of drought-induced microRNAs in rice. Biochem. Biophys. Res. Commun. 354, 585–590. Zhou, X., Wang, G., Zhang, W., 2007. UV-B responsive microRNA genes in Arabidopsis thaliana. Mol. Syst. Biol. 3.