Transcriptional regulatory networks in Arabidopsis thaliana during single and combined stresses

Abstract

Differentially evolved responses to various stress conditions in plants are controlled by complex regulatory circuits of transcriptional activators, and repressors, such as transcription factors (TFs). To understand the general and condition-specific activities of the TFs and their regulatory relationships with the target genes (TGs), we have used a homogeneous stress gene expression dataset generated on ten natural ecotypes of the model plant Arabidopsis thaliana, during five single and six combined stress conditions. Knowledge-based profiles of binding sites for 25 stress-responsive TF families (187 TFs) were generated and tested for their enrichment in the regulatory regions of the associated TGs. Condition-dependent regulatory sub-networks have shed light on the differential utilization of the underlying network topology, by stress-specific regulators and multifunctional regulators. The multifunctional regulators maintain the core stress response processes while the transient regulators confer the specificity to certain conditions. Clustering patterns of transcription factor binding sites (TFBS) have reflected the combinatorial nature of transcriptional regulation, and suggested the putative role of the homotypic clusters of TFBS towards maintaining transcriptional robustness against cis-regulatory mutations to facilitate the preservation of stress response processes. The Gene Ontology enrichment analysis of the TGs reflected sequential regulation of stress response mechanisms in plants.

Figure 1.

Methodology flow chart. (A) Top 500 stress-regulated genes were chosen from each of the eleven stress experiments performed on 10 Arabidopsis thaliana natural ecotypes. The united lists of stress-regulated genes contained 3429 genes, out of which 294 were TFs. A gene co-expression network was constructed using Pearson correlation coefficient (PCC ≥ 0.70), which contained 10 534 bipartite relations between 240 TFs and 1638 target genes (TGs). This network was later fed into the Network Component Analysis (NCA) algorithm, which then calculated the activity profiles of 194 stress-regulated TFs. Regulatory signal strength based filtering criteria (│r│≥0.75) was applied, which retained 5380 interactions between 182 TFs and 1199 TGs. Eleven regulatory sub-networks were extracted. Predicted connections were compared with independent interactions obtained from STIF analysis, AGRIS, AthaMap, AraNet and independent literature. (B) Flowchart of STIF algorithm, which was used for predicting transcription factor (TF) binding sites using HMM-based profile search. Total 187 stress-regulated TFs were classified into 39 TF Families and corresponding HMM profiles were generated. STIF algorithm predicted 10 474 interactions between 25 TF families and 1021 TGs at the Z-score threshold of 1.5. (adapted and modified with permission from (39))

Figure 2.

Overlapping statistics of stress-regulated TFs. Number of stress specific differentially regulated TFs are shown inside each circle along the periphery of the hendecagon. The small circles in the diagram represent overlapping TFs among different stress conditions. The arrow marks on the inner circles point towards the stress conditions in which they are significantly regulated. The size of the inner circles represents the number of TFs. The filled inner circle represents CONTSTANS-LIKE 5, which was significantly regulated in seven conditions. (Also see Supplementary Table S2).

Figure 3.

Topology of the predicted bipartite regulatory network visualized using Cytoscape software version 3.2.1 (169). (A) A topological overview of the predicted bipartite regulatory network between each pair of TF and TGs using PCC (│r│> 0.7). TFs are represented as triangles and TGs are represented as ellipses. This network is composed of 10 534 interactions (8450 positive interactions and 2084 negative interactions) between 240 TFs and 1638 TGs. The network topology has been simplified for presentation purpose. Detailed TF–TG interactions are provided in Supplementary Table S4. (B) The degree distribution of the nodes in the predicted regulatory network showed a power-law. (C) The degree distribution of the control network showed a Poisson distribution. (D) The topology of the regulatory sub-network during FLG stress condition. The TFs are represented as triangles and TGs are represented as ellipses. The size of the triangle reflects the number of connecting TGs to the corresponding TF. Previously reported key biotic stress regulators, such as WRKY22, MYB117, GNL, HSFA1A, MYB124, etc. are prominently visible as highly connected hubs in this network topology.

Figure 4.

Predicted activity profiles of 194 stress-regulated TFs as a heat map. Activity profiles of 194 stress-regulated TFs (x-axis) during eleven stress conditions (y-axis), in 10 Arabidopsis ecotypes, predicted using NCA. The activity values were scaled to +1 to −1. Hierarchical clustering of the TFs performed using Pearson correlation with average linkage method. The predicted profiles show the differential activity of the TFs in different stress conditions, as well as in various ecotypes.

Figure 5.

GSEA of the overlapping TGs. Network based visualization of GSEA results identified four clusters of statistically significant over-represented functional categories. Parental GO terms are shown in the centre. The differentially regulated primary processes were present in cluster 1, and secondary processes in cluster 2. The various stress response processes were distinctly clustered together in cluster 3. Processes that regulate the series of complex stress—response events were grouped together in cluster 4. These categories reflect the sequential modulation of stress-response processes. Upon exposure to stressed conditions, plants compromise with primary processes such as primary metabolism, photosynthesis, reproduction, to allocate extra resources to the activation of various defence and stress-response processes and secondary metabolic processes. Nodes were coloured according to their corrected P-values. Node sizes reflected the total number of genes in that category.