Mining expressed sequence tags of rapeseed (<em>Brassica napus</em> L.) to predict the drought responsive regulatory network
Abstract
It is of great significance to understand the regulatory mechanisms by which plants deal with drought stress. Two EST libraries derived from rapeseed (Brassica napus) leaves in non-stressed and drought stress conditions were analyzed in order to obtain the transcriptomic landscape of drought-exposed B. napus plants, and also to identify and characterize significant drought responsive regulatory genes and microRNAs. The functional ontology analysis revealed a substantial shift in the B. napus transcriptome to govern cellular drought responsiveness via different stress-activated mechanisms. The activity of transcription factor and protein kinase modules generally increased in response to drought stress. The 26 regulatory genes consisting of 17 transcription factor genes, eight protein kinase genes and one protein phosphatase gene were identified showing significant alterations in their expressions in response to drought stress. We also found the six microRNAs which were differentially expressed during drought stress supporting the involvement of a post-transcriptional level of regulation for B. napus drought response. The drought responsive regulatory network shed light on the significance of some regulatory components involved in biosynthesis and signaling of various plant hormones (abscisic acid, auxin and brassinosteroids), ubiquitin proteasome system, and signaling through Reactive Oxygen Species (ROS). Our findings suggested a complex and multi-level regulatory system modulating response to drought stress in B. napus.
Introduction
Growth and development of crop plants are severely affected by several abiotic stresses (Seki et al. 2003; Shinozaki and Yamaguchi-Shinozaki 1996). Among the abiotic stresses, drought has a crucial devastating impact on agricultural production (Farooq et al. 2009; Farooq et al. 2012). Plants recruit complex overlapping mechanisms to reprogram their gene expression to cope with drought stress (Ahuja et al. 2010; Golldack et al. 2011). The orchestrated regulatory networks at transcriptional, post transcriptional and post translational levels are involved in forming the appropriate transcriptome and proteome of stress-exposed plants to show an adaptive response (Duque et al. 2013; Mazzucotelli et al. 2008). Understanding the role of regulatory components will be of great benefit for breeding drought tolerant crops.
Various regulatory proteins including transcription factors, protein kinases and protein phosphatases have been identified to function as the core regulatory elements of stress responses (Ashraf 2010; Hadiarto and Tran 2011; Seki et al. 2007; Shinozaki and Yamaguchi-Shinozaki 1996; Umezawa et al. 2006; Xoconostle-Cazares et al. 2010; Yamaguchi-Shinozaki and Shinozaki 2006). Also, it has been revealed that the expression of some plant microRNAs (miRNAs), as negative regulator of gene expression, alters during stress conditions such as drought, salinity and cold (Khraiwesh et al. 2012; Lu and Huang 2008; Sunkar et al. 2007; Sunkar et al. 2012). Notably, most of these miRNAs target genes encoding transcription factors, which place miRNAs at the center of stress regulatory networks (Lu and Huang 2008; Sunkar et al. 2012).
Expressed sequence tag (EST) analysis is known as a powerful tool to discover and characterize genes involved in different stress responses and developmental cues in several plant species (Deokar et al. 2011; Gao et al. 2008; Gorantla et al. 2007; Gruber et al. 2012; Lata et al. 2010; Pratt et al. 2005; Shamloo-Dashtpagerdi et al. 2013; Zhuang and Zhu 2014). ESTs are also a valuable source for computational identification of miRNAs and study their possible roles in plant stress responses (Guo et al. 2014; Han et al. 2014; Hwang et al. 2011; Panda et al. 2014; Patanun et al. 2013; Ye et al. 2013).
Rapeseed (Brassica napus L.) is globally an important oilseed crop whose quality and quantity is often limited by drought stress in arid and semi-arid regions. Development of rapeseed cultivars with improved tolerance to drought is of great priority in rapeseed breeding programs (Diepenbrock 2000; Raymer 2002; Xie et al. 2007). A high level of homology exists in the protein-coding sequences of rapeseed and Arabidopsis thaliana making possible to use Arabidopsis genome information for functional interpretation in rapeseed (Fourmann et al. 2002; Parkin et al. 2002). EST-based approaches have been frequently employed to elucidate the genetic basis of regulation of rapeseed transcriptome in response to developmental cues, biotic and abiotic stresses, however little is known about the interconnections among different types of regulatory genes in the context of a network system. Generation and analysis of EST sequences from rapeseed developing seeds led to characterize a number of seed specific transcripts which had significant similarities to their Arabidopsis counterparts (Dong et al. 2004). The multiple regulatory factors involved in fatty acid metabolism were identified in rapeseed using EST chip hybridization (Niu et al. 2009). A comparative EST study revealed that the expression pattern of the rapeseed cotyledon transcriptome shifted towards a defense state after damage by flea beetles (Gruber et al. 2012). Indeed, the activity of the transcription factors known to promote or repress genes controlling stress responses, as well as genes involved in primary and secondary metabolism pathways, cell wall synthesis, and transport altered during flea beetle attack (Gruber et al. 2012). Analysis of rapeseed ESTs also discovered some mitogen-activated protein kinases and AP2/ERF transcription factors involved in abiotic stress signaling and regulatory mechanisms (Zhuang and Zhu 2014).
In this study, two EST libraries derived from rapeseed leaves in non-stress and drought stress conditions were mined in order to gain an insight into the transcriptomic landscape changes of drought-exposed rapeseed plants. Afterwards, the differentially expressed regulatory genes encoding transcription factors, protein kinases and protein phosphatases were identified. A computational analysis was also performed to identify and characterize rapeseed drought-responsive miRNAs and their possible target(s). Finally, a regulatory network was depicted which represented a complex and multi-level regulatory system to respond to drought stress in rapeseed.
Materials and methods
Source of B. napus expressed sequence tags data
Sequences of two 5 EST libraries from leaves of B. napus line DH12075, generated by Agriculture and Agri-food Canada (http://www.agr.gc.ca/), were retrieved from the TIGR Gene Indices of DFCI (http://compbio.dfci.harvard.edu/tgi/). One library (Lib ID: LIBEST 021477; http://www.ncbi.nlm.nih.gov/nucest/?term=LIBEST_021477) contained 6040 sequences from non-stress condition and the other (Lib ID: LIBEST_021491; http://www.ncbi.nlm.nih.gov/nucest/?term=LIBEST_021491) contained 13,063 sequences from drought stress condition. DH12075 is a double haploid line derived from an F1 cross between the French variety Cresor and the Canadian variety Westar (Young et al. 2004). Plants at the four-leaf stage were subjected to drought stress by applying 200 mM mannitol for 4 h and the cDNA libraries were constructed from small to medium sized healthy leaves using SuperScript Plasmid System with Gateway Technology for cDNA Synthesis and Cloning kit (Invitrogen) (http://www.ncbi.nlm.nih.gov/nucest/EV208948). EST sequences were the direct results of Base calling software Phred, and then they were processed using Lucy software. The analysis of these EST data has not been published so far.
Sequence processing, similarity search and annotation
The EST sequences were checked for vector contamination, sequence length and complexity using EGassembler, an online bioinformatics service (http://egassembler.hgc.jp) (Masoudi-Nejad et al. 2006). The vector, repetitive, chloroplast, and mitochondrial sequences were removed and the trimmed sequences were excluded from further analysis when they were less than 100 bp or they had greater than 4 % ambiguous bases (Carlson et al. 2006; Masoudi-Nejad et al. 2006). The remaining high quality ESTs of each library were subsequently clustered and assembled into unigenes (contigs and singletons) using EGassembler with Overlap percent identity cutoff >80 %. All the unigenes were searched against the Arabidopsis protein database (ftp://ftp.arabidopsis.org/) by local BLASTX using CLC Genomics Workbench software (version: 3.6.5) with the cut-off E-value ≤10. To categorize the unigenes with significant BLASTX scores into functional modules, modular enrichment analysis (MEA) was performed in each of the non-stress and drought stress conditions using Gene Annotations Co-occurrence Discovery; GeneCodis (http://genecodis.cnb.csic.es/) (Tabas-Madrid et al. 2012). GeneCodis is a web tool for functional enrichment analysis that integrates different information sources (GO, KEGG and Swiss-Prot gene accession databases) to find concurrent annotation of genes and form groups of genes with similar biological functions according to their statistical relevance (Carmona-Saez et al. 2007). In GeneCodis, Arabidopsis was selected as the reference organism and hypergeometric and FDR statistical tests were applied to compute and correct p-values, respectively. Radar charts were used to illustrate the percentage of functional modules.
Identification of drought-responsive regulatory components
In order to identify significant differentially expressed genes between the non-stress and drought EST libraries, all the EST sequences of the two libraries were clustered and assembled together using EGassembler with default parameters. The contigs were subsequently subjected to Audic and Claverie test (Audic and Claverie 1997) with α = 0.05, available at IDEG6 web tool (Romualdi et al. 2003), to check whether the number of ESTs of each library contributing to each contig were significantly different. The total number of ESTs in each library and the number of library specific ESTs in each contig were the inputs to perform Audic and Claverie test.
Of those differentially expressed contigs, the regulatory genes encoding transcription factors, protein kinases and protein phosphatases were specified by two online bioinformatic tools. The contigs were annotated and classified as transcription factors using PlantTFact (http://plantgrn.noble.org/PlantTFcat/) (Dai et al. 2013). The genes encoding protein kinases and protein phosphatases were identified by Plant Transcription factor and Protein Kinase Identifier and Classifier (iTAK) program (http://bioinfo.bti.cornell.edu/cgi-bin/itak/index.cgi).
We also performed a computational analysis to pinpoint and characterize significant differentially expressed miRNAs and their target genes. The Contigs matched with the regulatory genes in the previous step were excluded from the set of differentially expressed contigs. The remaining contigs were searched against all the available mature miRNAs from miRBase (http://www.mirbase.org/) (Griffiths-Jones et al. 2006) using CLC genomics workbench (version: 3.6.5) with expected values ≤10 to raise the hit chance for more potential sequences. The Contig sequences with the maximum four mismatches were subjected to similarity searches against non-redundant protein database using BLASTX (E-value ≤10) by CLC Genomics Workbench software (version: 3.6.5). Also, in order to distinguish between potential miRNAs and other small RNAs such as tRNA, rRNA, snRNA and snoRNA, a BLASTN search was carried out against Rfam database (rfam.sanger.ac.uk/). After that, the secondary structures of the remaining contigs were predicted using Mfold (Zuker 2003), The parameters were set as RNA sequence (linear), folding temperature (37 °C), maximum interior/ bulges loop size (30), and all others with default values. Also, the minimal folding free energy index (MFEI) was calculated for the predicted miRNAs (Zhang et al. 2006). The potential target genes of the predicted differentially expressed miRNAs were identified using psRNATarget with default parameters (http://plantgrn.noble.org/psRNATarget/?function=1) (Dai and Zhao 2011). All the differentially expressed genes found in this study were used as the database of psRNATarget. In each EST library, the normalized expression values of the regulatory genes and miRNAs were obtained using IDEG6 web tool as the ratio of the number of ESTs for the gene or miRNA of interest to the total number of ESTs. The fold change of expression for each of the genes and miRNAs was calculated by dividing the expression level in drought stress conditions by that of in non-stress conditions.
The identified regulatory components were networked using Pathway Studio software version 9 based on Fisher test with p-values ≤ 0.01 (Nikitin et al. 2003; Subramanian et al. 2005). The database of software could provide a number of relevant components, which were not identified among the differentially expressed genes, to make the regulatory network more complete.
Source of B. napus expressed sequence tags data
Sequences of two 5 EST libraries from leaves of B. napus line DH12075, generated by Agriculture and Agri-food Canada (http://www.agr.gc.ca/), were retrieved from the TIGR Gene Indices of DFCI (http://compbio.dfci.harvard.edu/tgi/). One library (Lib ID: LIBEST 021477; http://www.ncbi.nlm.nih.gov/nucest/?term=LIBEST_021477) contained 6040 sequences from non-stress condition and the other (Lib ID: LIBEST_021491; http://www.ncbi.nlm.nih.gov/nucest/?term=LIBEST_021491) contained 13,063 sequences from drought stress condition. DH12075 is a double haploid line derived from an F1 cross between the French variety Cresor and the Canadian variety Westar (Young et al. 2004). Plants at the four-leaf stage were subjected to drought stress by applying 200 mM mannitol for 4 h and the cDNA libraries were constructed from small to medium sized healthy leaves using SuperScript Plasmid System with Gateway Technology for cDNA Synthesis and Cloning kit (Invitrogen) (http://www.ncbi.nlm.nih.gov/nucest/EV208948). EST sequences were the direct results of Base calling software Phred, and then they were processed using Lucy software. The analysis of these EST data has not been published so far.
Sequence processing, similarity search and annotation
The EST sequences were checked for vector contamination, sequence length and complexity using EGassembler, an online bioinformatics service (http://egassembler.hgc.jp) (Masoudi-Nejad et al. 2006). The vector, repetitive, chloroplast, and mitochondrial sequences were removed and the trimmed sequences were excluded from further analysis when they were less than 100 bp or they had greater than 4 % ambiguous bases (Carlson et al. 2006; Masoudi-Nejad et al. 2006). The remaining high quality ESTs of each library were subsequently clustered and assembled into unigenes (contigs and singletons) using EGassembler with Overlap percent identity cutoff >80 %. All the unigenes were searched against the Arabidopsis protein database (ftp://ftp.arabidopsis.org/) by local BLASTX using CLC Genomics Workbench software (version: 3.6.5) with the cut-off E-value ≤10. To categorize the unigenes with significant BLASTX scores into functional modules, modular enrichment analysis (MEA) was performed in each of the non-stress and drought stress conditions using Gene Annotations Co-occurrence Discovery; GeneCodis (http://genecodis.cnb.csic.es/) (Tabas-Madrid et al. 2012). GeneCodis is a web tool for functional enrichment analysis that integrates different information sources (GO, KEGG and Swiss-Prot gene accession databases) to find concurrent annotation of genes and form groups of genes with similar biological functions according to their statistical relevance (Carmona-Saez et al. 2007). In GeneCodis, Arabidopsis was selected as the reference organism and hypergeometric and FDR statistical tests were applied to compute and correct p-values, respectively. Radar charts were used to illustrate the percentage of functional modules.
Identification of drought-responsive regulatory components
In order to identify significant differentially expressed genes between the non-stress and drought EST libraries, all the EST sequences of the two libraries were clustered and assembled together using EGassembler with default parameters. The contigs were subsequently subjected to Audic and Claverie test (Audic and Claverie 1997) with α = 0.05, available at IDEG6 web tool (Romualdi et al. 2003), to check whether the number of ESTs of each library contributing to each contig were significantly different. The total number of ESTs in each library and the number of library specific ESTs in each contig were the inputs to perform Audic and Claverie test.
Of those differentially expressed contigs, the regulatory genes encoding transcription factors, protein kinases and protein phosphatases were specified by two online bioinformatic tools. The contigs were annotated and classified as transcription factors using PlantTFact (http://plantgrn.noble.org/PlantTFcat/) (Dai et al. 2013). The genes encoding protein kinases and protein phosphatases were identified by Plant Transcription factor and Protein Kinase Identifier and Classifier (iTAK) program (http://bioinfo.bti.cornell.edu/cgi-bin/itak/index.cgi).
We also performed a computational analysis to pinpoint and characterize significant differentially expressed miRNAs and their target genes. The Contigs matched with the regulatory genes in the previous step were excluded from the set of differentially expressed contigs. The remaining contigs were searched against all the available mature miRNAs from miRBase (http://www.mirbase.org/) (Griffiths-Jones et al. 2006) using CLC genomics workbench (version: 3.6.5) with expected values ≤10 to raise the hit chance for more potential sequences. The Contig sequences with the maximum four mismatches were subjected to similarity searches against non-redundant protein database using BLASTX (E-value ≤10) by CLC Genomics Workbench software (version: 3.6.5). Also, in order to distinguish between potential miRNAs and other small RNAs such as tRNA, rRNA, snRNA and snoRNA, a BLASTN search was carried out against Rfam database (rfam.sanger.ac.uk/). After that, the secondary structures of the remaining contigs were predicted using Mfold (Zuker 2003), The parameters were set as RNA sequence (linear), folding temperature (37 °C), maximum interior/ bulges loop size (30), and all others with default values. Also, the minimal folding free energy index (MFEI) was calculated for the predicted miRNAs (Zhang et al. 2006). The potential target genes of the predicted differentially expressed miRNAs were identified using psRNATarget with default parameters (http://plantgrn.noble.org/psRNATarget/?function=1) (Dai and Zhao 2011). All the differentially expressed genes found in this study were used as the database of psRNATarget. In each EST library, the normalized expression values of the regulatory genes and miRNAs were obtained using IDEG6 web tool as the ratio of the number of ESTs for the gene or miRNA of interest to the total number of ESTs. The fold change of expression for each of the genes and miRNAs was calculated by dividing the expression level in drought stress conditions by that of in non-stress conditions.
The identified regulatory components were networked using Pathway Studio software version 9 based on Fisher test with p-values ≤ 0.01 (Nikitin et al. 2003; Subramanian et al. 2005). The database of software could provide a number of relevant components, which were not identified among the differentially expressed genes, to make the regulatory network more complete.
Results and discussion
Sequence analysis, annotation and functional classification
Pre-processing of 6040 and 13,063 EST sequences derived from B. napus leaves in non-stress and drought stress conditions resulted in 6037 and 13,059 high quality sequences with average length of 863 and 741 bp, respectively. ESTs obtained from the cDNA library of non-stressed plants were clustered and assembled into 4149 unigenes (794 contigs and 3355 singletons), in which the contigs encompassed 2682 (44.43 %) of EST sequences. Also, assembly of ESTs obtained from the cDNA library of drought-stressed plants yielded 5817 unigenes (1945 contigs and 3872 singletons) that 9191 (70.36 %) of ESTs fell into the contigs. The average contig lengths were 898 and 791 bp for the non-stress and drought EST libraries, respectively. The number of ESTs ranged from 2 (465 contigs) to 97 (one contig) in the non-stress contigs and 2 (1140 contigs) to 179 (one contig) in the drought stress contigs representing genes with different levels of expression. Consequently, total lengths of approximately 3395.7 and 3873.5 kb of the non-stress and drought EST libraries were analyzed (Table 1).
Table 1
Summary statistics of the B. napus ESTs obtained from non-stress and drought stress conditions
| Library | Number of high-quality ESTs | Average length of high-quality ESTs (bp) | Total number of unigenes | Number of contigs | Average length of contigs (bp) | Number of Singletons | Approximate total length of unigenes (Kbp) |
|---|---|---|---|---|---|---|---|
| Non-stress | 6037 | 863 | 4149 | 794 | 898 | 3355 | 3395.7 |
| Drought stress | 13059 | 741 | 5817 | 1945 | 791 | 3872 | 3873.5 |
Three thousand six hundred and four thousand seven hundred sixty-eight unigenes of the non-stress and drought stress B. napus transcriptomes were annotated based on BLASTX homology search against the Arabidopsis protein database. 80.82 and 75.45 % of the unigenes respective to non-stress and drought stress conditions showed significant similarity with genes of known or putative function (Table 2). 5.95 and 6.52 % of the unigenes from the non-stress and drought stress EST libraries were matched to genes with unknown function. The remaining proportion of the unigenes, i.e., 13.23 and 18.03 % corresponding to the non-stress and drought stress libraries respectively, showed no significant match with the Arabidopsis protein database implying that they may be specific genes to B. napus or genes whose biological functions have not yet been reported in Arabidopsis. Therefore, they can be considered as a source for gene discovery.
Table 2
Summary of BLASTX search for the B. napus unigenes obtained from non-stress and drought stress EST libraries against the Arabidopsis proteins database with expect value cut-off at 10
| Library | Unigenes with known function | Unigenes with putative function | Unigenes with unknown function | Unigenes with no significant matches |
|---|---|---|---|---|
| Non-stress | 3236 (78 %) | 117 (2.82 %) | 247 (5.95 %) | 549 (13.23 %) |
| Drought stress | 4165 (71.6 %) | 224 (3.85 %) | 379 (6.52 %) | 1049 (18.03 %) |
GeneCodis was utilized for functional enrichment of the unigenes with significant hits derived from the two EST libraries in the context of biological process and molecular function ontologies. Based on GeneCodis database, the 1573 and 2676 unigenes of the non-stress and drought stress libraries were ontologically informative and selected for enrichment analysis. Of those, the 988 and 2101 sequences were specific to the non-stress and drought stress EST libraries, respectively, while the 575 unigenes were common between them. The results are shown by Venn diagram in Fig. 1.
Venn diagram showing the common and specific informative unigenes obtained from the non-stress and drought EST libraries
To demonstrate how the genome functions altered in response to drought stress, modular enrichment analysis was conducted for all the informative unigenes obtained from each non-stress and drought stress conditions, and also for the specific informative unigenes of each condition. As a result, all the unigenes of the non-stress and drought stress libraries were enriched in 11 and 10 annotation modules based on biological process ontology (Fig. 2a). The two transcriptomes shared some common annotation modules namely “Photosynthesis”, “Response to cold”, “Response to cadmium ion”, “Defense response to bacterium”, “Response to salt stress” and “electron transport in photosystem Ι”. Over 50 % of the unigenes derived from non-stress conditions were enriched in non-stress responsive groups such as translation and photosynthesis. The results showed that drought stress caused significant changes in the abundance of unigenes assigned to various functional modules. As shown in Fig. 2a, the activity of the enriched annotation modules namely response to “salt stress”, “abscisic acid (ABA)”, “water deprivation”, “light stimulus”, “cadmium ion” and “glycolysis” significantly increased under drought stress conditions, indicating the dramatic change in B. napus genome functions toward response to environmental cues. On the other hand, the percentage of unigenes assigned to “translation”, “response to cold”, “photosynthesis”, “response to bacterium” and “response to far red, red and blue light” modules were significantly underrepresented under drought stress conditions. The results of modular enrichment analysis for the unigenes specific to each condition interestingly revealed that all the unigenes specifically expressed in drought stress fell into the functional modules involved in response to various stresses including cadmium ion, salt, cold, ABA, water deprivation and wounding (Fig. 2b).
The percentages of functional modules of (a) all and (b) specific informative unigenes derived from non-stress and drought stress B. napus transcriptomes according to biological process ontology
From point of view of molecular function ontology, all the informative unigenes derived from non-stress and drought stress conditions were enriched in 15 and 14 annotation modules, respectively; some of which were common for both conditions but with different abundances of the corresponding unigenes (Fig. 3a). The protein binding module encompassed the highest percentage of the unigenes in non-stress (24.82 %) and drought stress (22 %) conditions. The modules of identical protein binding (2.37 %) and Protein kinase activity (1.37 %) contained the lowest percentage of the unigenes obtained from non-stress and drought stress conditions, respectively. According to molecular function ontology, the activity of ATP binding, zinc ion binding, oxidoreductase, ubiquitin protein ligase, GTP binding, protein kinase, nucleotide binding and transcription factor functional modules enhanced in response to drought. Conversely, the number of transcripts assigned to the functional groups of protein binding, ribosomal constituent, RNA binding, copper ion binding, chlorophyll binding, calcium ion binding and ATPase activity declined under drought stress. Moreover, Fig. 3b illustrated that over 60 % of the unigenes specifically expressed in drought stress were assigned to the three modules including protein binding, ATP binding and catalytic activity. On the other hand, the specific unigenes of non-stress conditions showed different distribution among the functional modules and were divided into more distinct modules; however the module of protein binding still had the maximum percentage (30.45 %) of the unigenes.
The percentages of functional modules of (a) all and (b) specific informative unigenes derived from non-stress and drought stress B. napus transcriptomes according to molecular function ontology
The functional ontology analysis revealed a substantial shift in the B. napus transcriptome to govern cellular drought responsiveness via different stress-activated mechanisms including biosynthesis and signaling of various hormones, signaling and neutralization of Reactive Oxygen Species (ROS) and regulation of cellular energy homeostasis. The results were consistent with the fact that plant stress responses cause reduction in photosynthesis and cellular energy thereby plants need to constantly adjust their energy-associated transcriptome (Avin‐Wittenberg et al. 2012). Functional enrichment analysis in Arabidopsis and Rice also demonstrated significant changes in the functional modules related to response and regulation of various plant hormones and photosynthesis under drought and bacterial stresses (Shaik and Ramakrishna 2013).
Identification of the regulatory components
In order to discover the differentially expressed genes between the B. napus transcriptomes in non-stress and drought stress conditions, the total of 19,096 high quality ESTs generated from both EST libraries were assembled and consequently 2345 highly represented genes (contigs) comprising 9534 ESTs (49.93 %) were further analyzed. The number of ESTs in the contigs ranged from 2 (1322 contigs) to 180 (1 contig). The contigs showed the average length of 824 bp.
The 814 contigs with significantly differential expression between non-stress and drought stress conditions were identified using Audic and Claverie test (p-value ≤ 0.05). The contigs are from hereon named as genes. Of those, 441 genes which showed equal or greater than three differences in the number of ESTs between the two EST libraries were used to identify the highly differential regulatory genes and subsequently to infer a simplified plausible drought-responsive regulatory network. As a result, the 26 regulatory genes of the three major groups encoding transcription factors, protein kinases and protein phosphatases were annotated among the differentially expressed genes. The log2 fold change of expression of the genes was calculated and shown in Table 3.
Table 3
The list of proteins encoded by the regulatory genes with significant differential expressions (p-value ≤ 0.05) between the B. napus transcriptomes derived from non- and drought stress conditions. The proteins were annotated based on homology to Arabidopsis proteins
| Regulatory protein | Family | TAIR ID | Annotation | Fold change (Log2) |
|---|---|---|---|---|
| Transcription factor | C3H | AT4G29190 | Zinc finger C-x8-C-x5-C-x3-H type | 2.54 |
| FHA-SMAD | AT5G67030 | ABA1, zeaxanthin epoxidase (ZEP) | 2.18 | |
| WD40-like | AT1G21680 | DPP6 N-terminal domain-like | 1.87 | |
| CCHC(Zn) | AT2G24590 | RSZ22a | 1.67 | |
| AUX-IAA | AT3G04730 | IAA16, indoleacetic acid-induced protein 16 | 1.45 | |
| C2H2 | AT5G14420 | RGLG2, RING domain ligase2 | 1.44 | |
| GRAS | AT2G37650 | GRAS family transcription factor | 1.44 | |
| Homobox-WOX | AT2G22430 | ATHB6, homeobox protein 6 | 1.44 | |
| MADS-MIKC | AT2G45650 | AGL6, AGAMOUS-like 6 | 1.44 | |
| MYB-HB-like | AT5G05090 | Homeodomain-like | 1.44 | |
| WD40-like | AT1G21680 | DPP6 N-terminal domain-like | 1.44 | |
| bHLH | AT4G36540 | BEE2, BR enhanced expression 2 | −2.86 | |
| C2C2-CO-like | AT5G24930 | ATCOL4, CONSTANS-like 4 | −2.86 | |
| Hap2/NF-YA | AT5G12840 | HAP2A, nuclear factor Y subunit A1 | −2.86 | |
| MADS-MIKC | AT2G45660 | AGL20, SOC1,AGAMOUS-like 20 | −2.86 | |
| Znf-B | AT1G68190 | B-box zinc finger | −2.86 | |
| LIM | AT1G10200 | WLIM1 | −3.38 | |
| Protein kinase | Casein Kinase II | AT3G50000 | CKA2, casein kinase II, alpha chain 2 | 1.87 |
| Leucine Rich Repeat | AT4G33430 | BAK1, BRI1-associated receptor kinase | 1.67 | |
| MAPK | AT2G43790 | MPK6, MAP kinase 6 | 1.44 | |
| SNF1 Related Protein | AT4G30960 | CIPK6, SNRK3.14 | 1.44 | |
| CDC2 Like | AT1G18670 | IBS1, Protein kinase | −2.86 | |
| Receptor Like Kinase I | AT3G56050 | Protein kinase | −2.86 | |
| S Domain (Type 1) | AT1G56145 | Leucine-rich repeat protein kinase | −2.86 | |
| Serine/threonine | AT2G42960 | Protein kinase | −2.86 | |
| Protein phosphatase | Protein phosphatase 2CA | AT4G26080 | ABA insensitive 1; ABI1 | 2.02 |
The 17 transcription factor genes belonging to 15 different families were found indicating that B. napus has evolved a complex regulatory machinery to cope with drought stress. These genes showed different magnitudes and directions of expression changes in response to drought stress (Table 3). The expression of eleven transcription factors significantly increased whereas the six transcription factors showed decreased expression. The gene (AT4G29190) encoding a member of C3H transcription factor family showed the highest level of upregulation in response to drought stress, whilst WLIM1 transcription factor gene had the maximum downregulation following drought stress. In addition, eight protein kinase genes and one protein phosphatase gene were identified showing significant alterations in their expressions during drought stress (Table 3). A significant rise was observed in the expression of the four protein kinase genes (CKA2, BAK1, MPK6 and CIPK6) as well as the protein phosphatase gene, ABI1; while the other four protein kinase genes were sharply downregulated in response to drought stress.
We also found the six miRNAs which were significantly differentially expressed during drought stress revealing the involvement of a post-transcriptional level of regulation for B. napus drought response (Table 4). Five miRNAs were upregulated and one miRNA was downregulated under drought stress conditions. Interestingly, only one miRNAs namely miR396a has already been reported in B. napus and the other five miRNAs are newly reported in this study. There were four miRNAs (miR5658, miR396b, miR1536 and miR1509) with 21 nucleotides in length, while the other two miRNAs (miR2919 and miR435) had the size of 19 and 20 nucleotides, respectively. Further analysis predicted that the seven identified regulatory genes with differential expression between non-stress and drought conditions, including six transcription factor genes and one protein kinase gene were the direct targets of the five differentially expressed miRNAs (Table 4). No target gene was found for miR1509b among the differential regulatory genes. It is noteworthy that miR1536 targeted the three differentially expressed transcription factors (AGL6, IAA16 and RGLG2) showing its importance in regulating drought responses of B. napus.
Table 4
The miRNAs with significant differential expressions (p-value ≤ 0.05) between the B. napus transcriptomes derived from non- and drought stress conditions. The best matched miRNAs obtained from miRBase database and their potential target genes are also represented
| B. napus mature miRNA sequence | The best matched miRNA | Fold change (Log2) | Target gene | UPE | Inhibition mechanism |
|---|---|---|---|---|---|
| AUGAUUAUGAUUAUGAUGAGG | ath-miR5658 | 2.37 | BAK1 | 8.098 | Translation |
| AAAGGGGGGGGGGGAAUUU | osa-miR2919 | 0.91 | ABA1 | 22.646 | Cleavage |
| UUCCACAACUUUCUUCAAGUC | Bna-miR396a | 0.91 | B-box zinc finger | 17.6 | Translation |
| UAGAAAAGACAAAUCUGUUUA | gma-miR1536 | 0.91 | AGL6, IAA16, RGLG2 | 13.749 | Cleavage |
| UUUCCCGGUAUUGGACUUGG | osa-miR435 | 0.91 | COL4 | 21.26 | Translation |
| UAAAUCAAGGAAAGAAGGGUU | gma-miR1509b | −2.85 | – | – | – |
Maximum energy to unpair the target site. The less energy means the more possibility that small RNA is able to contact (and cleave) target mRNA
The regulatory network was predicted among the identified components (Fig. 4). Eleven differentially expressed regulatory genes including seven transcription factors, three protein kinases and one protein phosphatases were placed in the network. Also, two genes namely ABI2 and BRI1 were added to the network by Pathway Studio software. Furthermore, the four identified drought responsive miRNAs (miR5658, miR1536, miR435 and miR2919) were connected to their targets. The network uncovered the relationships among the differentially expressed regulatory genes and the interactions between various regulatory pathways towards an adaptive response to drought in B. napus. These are detailed as follows.
Regulation of ABA biosynthesis
Notably, all the differential regulatory components within the predicted network were, directly or indirectly, connected to the phytohormone ABA, highlighting the significant roles of ABA in B. napus drought responsiveness. Alteration of cellular ABA levels is essential for several developmental and stress responsive processes (Fujita et al. 2011). The results showed that drought stress directly and significantly activated the expression of B. napus ortholog of ABA1 involved in the first step of ABA biosynthesis. The expression level of ABA1 was clearly upregulated by ABA, drought and salt in Arabidopsis shoots (Xiong et al. 2002), whereas its expression was induced by drought in tomato roots but not leaves (Thompson et al. 2000). The regulatory network presented a mutual connection between ABA and ABA1 that is consistent with that ABA synthesis is under a positive feedback regulation as the expression of some genes involved in ABA metabolism are regulated by ABA (Xiong and Zhu 2003). Nevertheless, B. napus ABA1 was known as the potential target of miR2919 which was induced by drought stress. This finding suggested a possible negative regulatory mechanism which antagonistically acts against positive feedback mechanism to precisely regulate ABA1 expression and subsequently ABA level.
Stress signaling through ROS regulation
In this study, we identified a protein phosphatase gene (ABI1), which its expression was significantly enhanced by drought stress. In Arabidopsis, ABI1 and ABI2 are the members of protein phosphatase 2C (PP2C) group and act as the negative regulators of ABA signaling (Leung et al. 1997). As shown in Fig. 4, ABI1 expression is under the direct influence of drought stress, and ABA hormone. ABI1 mediates signal transduction between ABA reception and ROS production under stress conditions (Murata et al. 2001). The Previous studies indicated that MAP kinase signaling is activated following ROS accumulation, as testified to the induction of MPK3 and MPK6 by H2O2 (Jammes et al. 2009; Xing et al. 2008). Interestingly, in this study, the expression of a B. napus homolog of MPK6 was significantly upregulated during drought stress. Arabidopsis MPK6 promotes stomatal closure under osmotic stress (Liu 2012). MPK6 recruits ethylene signaling pathway via phosphorylation and thereby activation of the Arabidopsis transcription factor ERF6 (ETHYLEN RESPONSE FACTOR 6) (Meng et al. 2013). ERF6 specifically binds to the ROS-responsive elements (ROSEs) and interacts with MPK6 to regulate ROS-responsive gene expression at the transcriptional level (Wang et al. 2013). As a result, the identification of the B. napus ABI1 and MPK6 provided evidence to support the role of ROS-mediated-signaling in coping with drought stress.
Response to drought stress via the ubiquitin 26S proteasome proteolytic pathway
The expression of RGLG2, a central component of the ubiquitin 26S proteasome pathway, was significantly downregulated in the drought-stressed B. napus transcriptome. Arabidopsis RGLG2 is a RING domain E3 ligase which mediates ubiquitination of the transcription factor Ethylene Response Factor 53 (ERF53) for proteasome degradation (Cheng et al. 2012). The expression of some drought-responsive transcription factors genes such as DREB2A, CBF4 and RAP2.4 is controlled by ERF53 (Haake et al. 2002; Lin et al. 2008; Sakuma et al. 2006). Therefore, downregulation of RGLG2 expression is required for induction of ERF53 and thereby activation of its downstream genes under drought stress. This study revealed that the B. napus RGLG2 is one of the possible targets of miR1536 which was identified for the first time in the B. napus transcriptome. Similar to ERF53, the involvement of a ubiquitin E3 ligase was known to regulate the abundance of the drought responsive homeobox-leucine zipper transcription factor, AtHB6 (Lechner et al. 2011). AtHB6 is a target of ABI1 and a negative regulator of ABA signaling pathway in Arabidopsis (Himmelbach et al. 2002). The significant upregulation of the B. napus HB6 expression along with the downregulation of the B. napus RGLG2 expression under drought stress resembled the regulatory mechanism, as occurred in Arabidopsis.
Drought-induced Ca signaling and auxin response
The transcript level of B. napus CIPK6 gene was significantly upregulated during drought stress. Arabidopsis CIPK6 encodes a protein kinase involved in sensing Ca signal and induction of Ca-dependent ABA signaling (Batistič and Kudla 2012; Halford and Hey 2009; Kim and Maik 2010; Luan 2009; Weinl and Kudla 2009). The overexpression of GhCIPK6 from cotton significantly enhanced tolerance to salt, drought and ABA stresses in transgenic Arabidopsis (He et al. 2013). CIPK6 is also known to function in auxin transport (Tripathi et al. 2009). It has been demonstrated that auxin might participate in positive regulation of drought tolerance through regulation of root architecture, expression of ABA-responsive genes, ROS metabolism and metabolic homeostasis (Shi et al. 2014). In accordance with the activation of B. napus CIPK6 and auxin in response to drought, the results also showed significant downregulation of the B. napus IAA16, a negative regulator of auxin signaling (Park et al. 2011), under drought stress by the increased activity of miR1536. Taken together, these findings uncovered the significance of ABA-mediated-auxin signaling in the modulation of B. napus drought stress responses.
The role of brassinosteroids
The results presented the contribution of brassinosteroids (BRs) to drought response in B. napus. BRs play essential roles in regulating various biological processes including responses to the environment (Wang et al. 2014). BR signaling is induced by binding of BR to a membrane receptor kinase, BRI1. BRI1 makes hetero oligomers with another receptor kinase, BAK1 to shape a phosphorylation string that ultimately causes functional changes of target genes (Nam and Li 2002; Wang et al. 2008). The BRI1 receptor complex induces the redundant genes encoding basic helix-loop-helix (bHLH) namely BEE1, BEE2 and BEE3 (Friedrichsen et al. 2002). The B. napus transcriptome showed a significant increase in the expression of BAK1 under drought stress, although the transcripts of this gene were the potential targets of the drought responsive miR5658 indicating the existence of a negative post-transcriptional regulation for BAK1. On the other hand, the expression of BEE2 was significantly downregulated under drought stress in B. napus leaves. As illustrated in the regulatory network, BEE2 was under the negative regulation of ABA (Friedrichsen et al. 2002). These results showed that ABA and BR signaling pathways act antagonistically in response to drought stress in B. napus. BR and ABA antagonize each other in many biological processes (Wang et al. 2014).
Developmental regulation during drought stress
The expression of the B. napus homologs of the transcription factor genes CONSTANS-like 4 (COL4) and AGAMOUS-LIKE 20 (AGL20) (also known as SOC1) were significantly downregulated under drought stress. The drought responsive regulatory network (Fig. 4) highlighted the induction of miR435 which potentially targets the expression of the B. napus COL4. In agreement with our results, the expression of a number of CONSTANS-like family members was downregulated by drought stress in rice (Shaik and Ramakrishna 2013). The expression of Arabidopsis COL4 was strongly induced by ABA, salt, and osmotic stress (Min et al. 2014). It has been demonstrated that the Arabidopsis SOC1, as a positive regulator of flowering, was upregulated following 4 to 5 days drought treatment indicating that plants reprogram flowering process to shorten the life cycle and promote survival under extreme environments (Su et al. 2013). The apparently contrary responses of these genes to drought may be due to different timing and intensity of stress, and even the existence of various strategies to set time to flowering under stress conditions in plants. Further investigations are required to determine the regulatory roles of B. napus COL4 and SOC1 under stress conditions.
In conclusion, this EST-based transcriptome study provided an insight into how B. napus genome responds to drought stress. The results clearly showed a substantial change in the B. napus genome expression in response to drought in which the corresponding regulatory mechanisms were significantly activated. We identified a number of significant drought responsive regulatory genes and miRNAs with various functions and characterized their interactions using the regulatory network. Our findings indicated that B. napus recruits a complex and overlapping regulatory system to modulate response to drought stress. The results highlighted the importance of biosynthesis and signaling of various plant hormones and their crosstalks to coordinate physiological events under drought stress in B. napus. The obtained drought responsive regulatory network can be utilized to select appropriate genes for transgenic breeding to develop B. napus drought tolerant cultivars.
Sequence analysis, annotation and functional classification
Pre-processing of 6040 and 13,063 EST sequences derived from B. napus leaves in non-stress and drought stress conditions resulted in 6037 and 13,059 high quality sequences with average length of 863 and 741 bp, respectively. ESTs obtained from the cDNA library of non-stressed plants were clustered and assembled into 4149 unigenes (794 contigs and 3355 singletons), in which the contigs encompassed 2682 (44.43 %) of EST sequences. Also, assembly of ESTs obtained from the cDNA library of drought-stressed plants yielded 5817 unigenes (1945 contigs and 3872 singletons) that 9191 (70.36 %) of ESTs fell into the contigs. The average contig lengths were 898 and 791 bp for the non-stress and drought EST libraries, respectively. The number of ESTs ranged from 2 (465 contigs) to 97 (one contig) in the non-stress contigs and 2 (1140 contigs) to 179 (one contig) in the drought stress contigs representing genes with different levels of expression. Consequently, total lengths of approximately 3395.7 and 3873.5 kb of the non-stress and drought EST libraries were analyzed (Table 1).
Table 1
Summary statistics of the B. napus ESTs obtained from non-stress and drought stress conditions
| Library | Number of high-quality ESTs | Average length of high-quality ESTs (bp) | Total number of unigenes | Number of contigs | Average length of contigs (bp) | Number of Singletons | Approximate total length of unigenes (Kbp) |
|---|---|---|---|---|---|---|---|
| Non-stress | 6037 | 863 | 4149 | 794 | 898 | 3355 | 3395.7 |
| Drought stress | 13059 | 741 | 5817 | 1945 | 791 | 3872 | 3873.5 |
Three thousand six hundred and four thousand seven hundred sixty-eight unigenes of the non-stress and drought stress B. napus transcriptomes were annotated based on BLASTX homology search against the Arabidopsis protein database. 80.82 and 75.45 % of the unigenes respective to non-stress and drought stress conditions showed significant similarity with genes of known or putative function (Table 2). 5.95 and 6.52 % of the unigenes from the non-stress and drought stress EST libraries were matched to genes with unknown function. The remaining proportion of the unigenes, i.e., 13.23 and 18.03 % corresponding to the non-stress and drought stress libraries respectively, showed no significant match with the Arabidopsis protein database implying that they may be specific genes to B. napus or genes whose biological functions have not yet been reported in Arabidopsis. Therefore, they can be considered as a source for gene discovery.
Table 2
Summary of BLASTX search for the B. napus unigenes obtained from non-stress and drought stress EST libraries against the Arabidopsis proteins database with expect value cut-off at 10
| Library | Unigenes with known function | Unigenes with putative function | Unigenes with unknown function | Unigenes with no significant matches |
|---|---|---|---|---|
| Non-stress | 3236 (78 %) | 117 (2.82 %) | 247 (5.95 %) | 549 (13.23 %) |
| Drought stress | 4165 (71.6 %) | 224 (3.85 %) | 379 (6.52 %) | 1049 (18.03 %) |
GeneCodis was utilized for functional enrichment of the unigenes with significant hits derived from the two EST libraries in the context of biological process and molecular function ontologies. Based on GeneCodis database, the 1573 and 2676 unigenes of the non-stress and drought stress libraries were ontologically informative and selected for enrichment analysis. Of those, the 988 and 2101 sequences were specific to the non-stress and drought stress EST libraries, respectively, while the 575 unigenes were common between them. The results are shown by Venn diagram in Fig. 1.
Venn diagram showing the common and specific informative unigenes obtained from the non-stress and drought EST libraries
To demonstrate how the genome functions altered in response to drought stress, modular enrichment analysis was conducted for all the informative unigenes obtained from each non-stress and drought stress conditions, and also for the specific informative unigenes of each condition. As a result, all the unigenes of the non-stress and drought stress libraries were enriched in 11 and 10 annotation modules based on biological process ontology (Fig. 2a). The two transcriptomes shared some common annotation modules namely “Photosynthesis”, “Response to cold”, “Response to cadmium ion”, “Defense response to bacterium”, “Response to salt stress” and “electron transport in photosystem Ι”. Over 50 % of the unigenes derived from non-stress conditions were enriched in non-stress responsive groups such as translation and photosynthesis. The results showed that drought stress caused significant changes in the abundance of unigenes assigned to various functional modules. As shown in Fig. 2a, the activity of the enriched annotation modules namely response to “salt stress”, “abscisic acid (ABA)”, “water deprivation”, “light stimulus”, “cadmium ion” and “glycolysis” significantly increased under drought stress conditions, indicating the dramatic change in B. napus genome functions toward response to environmental cues. On the other hand, the percentage of unigenes assigned to “translation”, “response to cold”, “photosynthesis”, “response to bacterium” and “response to far red, red and blue light” modules were significantly underrepresented under drought stress conditions. The results of modular enrichment analysis for the unigenes specific to each condition interestingly revealed that all the unigenes specifically expressed in drought stress fell into the functional modules involved in response to various stresses including cadmium ion, salt, cold, ABA, water deprivation and wounding (Fig. 2b).
The percentages of functional modules of (a) all and (b) specific informative unigenes derived from non-stress and drought stress B. napus transcriptomes according to biological process ontology
From point of view of molecular function ontology, all the informative unigenes derived from non-stress and drought stress conditions were enriched in 15 and 14 annotation modules, respectively; some of which were common for both conditions but with different abundances of the corresponding unigenes (Fig. 3a). The protein binding module encompassed the highest percentage of the unigenes in non-stress (24.82 %) and drought stress (22 %) conditions. The modules of identical protein binding (2.37 %) and Protein kinase activity (1.37 %) contained the lowest percentage of the unigenes obtained from non-stress and drought stress conditions, respectively. According to molecular function ontology, the activity of ATP binding, zinc ion binding, oxidoreductase, ubiquitin protein ligase, GTP binding, protein kinase, nucleotide binding and transcription factor functional modules enhanced in response to drought. Conversely, the number of transcripts assigned to the functional groups of protein binding, ribosomal constituent, RNA binding, copper ion binding, chlorophyll binding, calcium ion binding and ATPase activity declined under drought stress. Moreover, Fig. 3b illustrated that over 60 % of the unigenes specifically expressed in drought stress were assigned to the three modules including protein binding, ATP binding and catalytic activity. On the other hand, the specific unigenes of non-stress conditions showed different distribution among the functional modules and were divided into more distinct modules; however the module of protein binding still had the maximum percentage (30.45 %) of the unigenes.
The percentages of functional modules of (a) all and (b) specific informative unigenes derived from non-stress and drought stress B. napus transcriptomes according to molecular function ontology
The functional ontology analysis revealed a substantial shift in the B. napus transcriptome to govern cellular drought responsiveness via different stress-activated mechanisms including biosynthesis and signaling of various hormones, signaling and neutralization of Reactive Oxygen Species (ROS) and regulation of cellular energy homeostasis. The results were consistent with the fact that plant stress responses cause reduction in photosynthesis and cellular energy thereby plants need to constantly adjust their energy-associated transcriptome (Avin‐Wittenberg et al. 2012). Functional enrichment analysis in Arabidopsis and Rice also demonstrated significant changes in the functional modules related to response and regulation of various plant hormones and photosynthesis under drought and bacterial stresses (Shaik and Ramakrishna 2013).
Identification of the regulatory components
In order to discover the differentially expressed genes between the B. napus transcriptomes in non-stress and drought stress conditions, the total of 19,096 high quality ESTs generated from both EST libraries were assembled and consequently 2345 highly represented genes (contigs) comprising 9534 ESTs (49.93 %) were further analyzed. The number of ESTs in the contigs ranged from 2 (1322 contigs) to 180 (1 contig). The contigs showed the average length of 824 bp.
The 814 contigs with significantly differential expression between non-stress and drought stress conditions were identified using Audic and Claverie test (p-value ≤ 0.05). The contigs are from hereon named as genes. Of those, 441 genes which showed equal or greater than three differences in the number of ESTs between the two EST libraries were used to identify the highly differential regulatory genes and subsequently to infer a simplified plausible drought-responsive regulatory network. As a result, the 26 regulatory genes of the three major groups encoding transcription factors, protein kinases and protein phosphatases were annotated among the differentially expressed genes. The log2 fold change of expression of the genes was calculated and shown in Table 3.
Table 3
The list of proteins encoded by the regulatory genes with significant differential expressions (p-value ≤ 0.05) between the B. napus transcriptomes derived from non- and drought stress conditions. The proteins were annotated based on homology to Arabidopsis proteins
| Regulatory protein | Family | TAIR ID | Annotation | Fold change (Log2) |
|---|---|---|---|---|
| Transcription factor | C3H | AT4G29190 | Zinc finger C-x8-C-x5-C-x3-H type | 2.54 |
| FHA-SMAD | AT5G67030 | ABA1, zeaxanthin epoxidase (ZEP) | 2.18 | |
| WD40-like | AT1G21680 | DPP6 N-terminal domain-like | 1.87 | |
| CCHC(Zn) | AT2G24590 | RSZ22a | 1.67 | |
| AUX-IAA | AT3G04730 | IAA16, indoleacetic acid-induced protein 16 | 1.45 | |
| C2H2 | AT5G14420 | RGLG2, RING domain ligase2 | 1.44 | |
| GRAS | AT2G37650 | GRAS family transcription factor | 1.44 | |
| Homobox-WOX | AT2G22430 | ATHB6, homeobox protein 6 | 1.44 | |
| MADS-MIKC | AT2G45650 | AGL6, AGAMOUS-like 6 | 1.44 | |
| MYB-HB-like | AT5G05090 | Homeodomain-like | 1.44 | |
| WD40-like | AT1G21680 | DPP6 N-terminal domain-like | 1.44 | |
| bHLH | AT4G36540 | BEE2, BR enhanced expression 2 | −2.86 | |
| C2C2-CO-like | AT5G24930 | ATCOL4, CONSTANS-like 4 | −2.86 | |
| Hap2/NF-YA | AT5G12840 | HAP2A, nuclear factor Y subunit A1 | −2.86 | |
| MADS-MIKC | AT2G45660 | AGL20, SOC1,AGAMOUS-like 20 | −2.86 | |
| Znf-B | AT1G68190 | B-box zinc finger | −2.86 | |
| LIM | AT1G10200 | WLIM1 | −3.38 | |
| Protein kinase | Casein Kinase II | AT3G50000 | CKA2, casein kinase II, alpha chain 2 | 1.87 |
| Leucine Rich Repeat | AT4G33430 | BAK1, BRI1-associated receptor kinase | 1.67 | |
| MAPK | AT2G43790 | MPK6, MAP kinase 6 | 1.44 | |
| SNF1 Related Protein | AT4G30960 | CIPK6, SNRK3.14 | 1.44 | |
| CDC2 Like | AT1G18670 | IBS1, Protein kinase | −2.86 | |
| Receptor Like Kinase I | AT3G56050 | Protein kinase | −2.86 | |
| S Domain (Type 1) | AT1G56145 | Leucine-rich repeat protein kinase | −2.86 | |
| Serine/threonine | AT2G42960 | Protein kinase | −2.86 | |
| Protein phosphatase | Protein phosphatase 2CA | AT4G26080 | ABA insensitive 1; ABI1 | 2.02 |
The 17 transcription factor genes belonging to 15 different families were found indicating that B. napus has evolved a complex regulatory machinery to cope with drought stress. These genes showed different magnitudes and directions of expression changes in response to drought stress (Table 3). The expression of eleven transcription factors significantly increased whereas the six transcription factors showed decreased expression. The gene (AT4G29190) encoding a member of C3H transcription factor family showed the highest level of upregulation in response to drought stress, whilst WLIM1 transcription factor gene had the maximum downregulation following drought stress. In addition, eight protein kinase genes and one protein phosphatase gene were identified showing significant alterations in their expressions during drought stress (Table 3). A significant rise was observed in the expression of the four protein kinase genes (CKA2, BAK1, MPK6 and CIPK6) as well as the protein phosphatase gene, ABI1; while the other four protein kinase genes were sharply downregulated in response to drought stress.
We also found the six miRNAs which were significantly differentially expressed during drought stress revealing the involvement of a post-transcriptional level of regulation for B. napus drought response (Table 4). Five miRNAs were upregulated and one miRNA was downregulated under drought stress conditions. Interestingly, only one miRNAs namely miR396a has already been reported in B. napus and the other five miRNAs are newly reported in this study. There were four miRNAs (miR5658, miR396b, miR1536 and miR1509) with 21 nucleotides in length, while the other two miRNAs (miR2919 and miR435) had the size of 19 and 20 nucleotides, respectively. Further analysis predicted that the seven identified regulatory genes with differential expression between non-stress and drought conditions, including six transcription factor genes and one protein kinase gene were the direct targets of the five differentially expressed miRNAs (Table 4). No target gene was found for miR1509b among the differential regulatory genes. It is noteworthy that miR1536 targeted the three differentially expressed transcription factors (AGL6, IAA16 and RGLG2) showing its importance in regulating drought responses of B. napus.
Table 4
The miRNAs with significant differential expressions (p-value ≤ 0.05) between the B. napus transcriptomes derived from non- and drought stress conditions. The best matched miRNAs obtained from miRBase database and their potential target genes are also represented
| B. napus mature miRNA sequence | The best matched miRNA | Fold change (Log2) | Target gene | UPE | Inhibition mechanism |
|---|---|---|---|---|---|
| AUGAUUAUGAUUAUGAUGAGG | ath-miR5658 | 2.37 | BAK1 | 8.098 | Translation |
| AAAGGGGGGGGGGGAAUUU | osa-miR2919 | 0.91 | ABA1 | 22.646 | Cleavage |
| UUCCACAACUUUCUUCAAGUC | Bna-miR396a | 0.91 | B-box zinc finger | 17.6 | Translation |
| UAGAAAAGACAAAUCUGUUUA | gma-miR1536 | 0.91 | AGL6, IAA16, RGLG2 | 13.749 | Cleavage |
| UUUCCCGGUAUUGGACUUGG | osa-miR435 | 0.91 | COL4 | 21.26 | Translation |
| UAAAUCAAGGAAAGAAGGGUU | gma-miR1509b | −2.85 | – | – | – |
Maximum energy to unpair the target site. The less energy means the more possibility that small RNA is able to contact (and cleave) target mRNA
The regulatory network was predicted among the identified components (Fig. 4). Eleven differentially expressed regulatory genes including seven transcription factors, three protein kinases and one protein phosphatases were placed in the network. Also, two genes namely ABI2 and BRI1 were added to the network by Pathway Studio software. Furthermore, the four identified drought responsive miRNAs (miR5658, miR1536, miR435 and miR2919) were connected to their targets. The network uncovered the relationships among the differentially expressed regulatory genes and the interactions between various regulatory pathways towards an adaptive response to drought in B. napus. These are detailed as follows.
The network showing interactions among the drought responsive identified in B. napus
Regulation of ABA biosynthesis
Notably, all the differential regulatory components within the predicted network were, directly or indirectly, connected to the phytohormone ABA, highlighting the significant roles of ABA in B. napus drought responsiveness. Alteration of cellular ABA levels is essential for several developmental and stress responsive processes (Fujita et al. 2011). The results showed that drought stress directly and significantly activated the expression of B. napus ortholog of ABA1 involved in the first step of ABA biosynthesis. The expression level of ABA1 was clearly upregulated by ABA, drought and salt in Arabidopsis shoots (Xiong et al. 2002), whereas its expression was induced by drought in tomato roots but not leaves (Thompson et al. 2000). The regulatory network presented a mutual connection between ABA and ABA1 that is consistent with that ABA synthesis is under a positive feedback regulation as the expression of some genes involved in ABA metabolism are regulated by ABA (Xiong and Zhu 2003). Nevertheless, B. napus ABA1 was known as the potential target of miR2919 which was induced by drought stress. This finding suggested a possible negative regulatory mechanism which antagonistically acts against positive feedback mechanism to precisely regulate ABA1 expression and subsequently ABA level.
Stress signaling through ROS regulation
In this study, we identified a protein phosphatase gene (ABI1), which its expression was significantly enhanced by drought stress. In Arabidopsis, ABI1 and ABI2 are the members of protein phosphatase 2C (PP2C) group and act as the negative regulators of ABA signaling (Leung et al. 1997). As shown in Fig. 4, ABI1 expression is under the direct influence of drought stress, and ABA hormone. ABI1 mediates signal transduction between ABA reception and ROS production under stress conditions (Murata et al. 2001). The Previous studies indicated that MAP kinase signaling is activated following ROS accumulation, as testified to the induction of MPK3 and MPK6 by H2O2 (Jammes et al. 2009; Xing et al. 2008). Interestingly, in this study, the expression of a B. napus homolog of MPK6 was significantly upregulated during drought stress. Arabidopsis MPK6 promotes stomatal closure under osmotic stress (Liu 2012). MPK6 recruits ethylene signaling pathway via phosphorylation and thereby activation of the Arabidopsis transcription factor ERF6 (ETHYLEN RESPONSE FACTOR 6) (Meng et al. 2013). ERF6 specifically binds to the ROS-responsive elements (ROSEs) and interacts with MPK6 to regulate ROS-responsive gene expression at the transcriptional level (Wang et al. 2013). As a result, the identification of the B. napus ABI1 and MPK6 provided evidence to support the role of ROS-mediated-signaling in coping with drought stress.
Response to drought stress via the ubiquitin 26S proteasome proteolytic pathway
The expression of RGLG2, a central component of the ubiquitin 26S proteasome pathway, was significantly downregulated in the drought-stressed B. napus transcriptome. Arabidopsis RGLG2 is a RING domain E3 ligase which mediates ubiquitination of the transcription factor Ethylene Response Factor 53 (ERF53) for proteasome degradation (Cheng et al. 2012). The expression of some drought-responsive transcription factors genes such as DREB2A, CBF4 and RAP2.4 is controlled by ERF53 (Haake et al. 2002; Lin et al. 2008; Sakuma et al. 2006). Therefore, downregulation of RGLG2 expression is required for induction of ERF53 and thereby activation of its downstream genes under drought stress. This study revealed that the B. napus RGLG2 is one of the possible targets of miR1536 which was identified for the first time in the B. napus transcriptome. Similar to ERF53, the involvement of a ubiquitin E3 ligase was known to regulate the abundance of the drought responsive homeobox-leucine zipper transcription factor, AtHB6 (Lechner et al. 2011). AtHB6 is a target of ABI1 and a negative regulator of ABA signaling pathway in Arabidopsis (Himmelbach et al. 2002). The significant upregulation of the B. napus HB6 expression along with the downregulation of the B. napus RGLG2 expression under drought stress resembled the regulatory mechanism, as occurred in Arabidopsis.
Drought-induced Ca signaling and auxin response
The transcript level of B. napus CIPK6 gene was significantly upregulated during drought stress. Arabidopsis CIPK6 encodes a protein kinase involved in sensing Ca signal and induction of Ca-dependent ABA signaling (Batistič and Kudla 2012; Halford and Hey 2009; Kim and Maik 2010; Luan 2009; Weinl and Kudla 2009). The overexpression of GhCIPK6 from cotton significantly enhanced tolerance to salt, drought and ABA stresses in transgenic Arabidopsis (He et al. 2013). CIPK6 is also known to function in auxin transport (Tripathi et al. 2009). It has been demonstrated that auxin might participate in positive regulation of drought tolerance through regulation of root architecture, expression of ABA-responsive genes, ROS metabolism and metabolic homeostasis (Shi et al. 2014). In accordance with the activation of B. napus CIPK6 and auxin in response to drought, the results also showed significant downregulation of the B. napus IAA16, a negative regulator of auxin signaling (Park et al. 2011), under drought stress by the increased activity of miR1536. Taken together, these findings uncovered the significance of ABA-mediated-auxin signaling in the modulation of B. napus drought stress responses.
The role of brassinosteroids
The results presented the contribution of brassinosteroids (BRs) to drought response in B. napus. BRs play essential roles in regulating various biological processes including responses to the environment (Wang et al. 2014). BR signaling is induced by binding of BR to a membrane receptor kinase, BRI1. BRI1 makes hetero oligomers with another receptor kinase, BAK1 to shape a phosphorylation string that ultimately causes functional changes of target genes (Nam and Li 2002; Wang et al. 2008). The BRI1 receptor complex induces the redundant genes encoding basic helix-loop-helix (bHLH) namely BEE1, BEE2 and BEE3 (Friedrichsen et al. 2002). The B. napus transcriptome showed a significant increase in the expression of BAK1 under drought stress, although the transcripts of this gene were the potential targets of the drought responsive miR5658 indicating the existence of a negative post-transcriptional regulation for BAK1. On the other hand, the expression of BEE2 was significantly downregulated under drought stress in B. napus leaves. As illustrated in the regulatory network, BEE2 was under the negative regulation of ABA (Friedrichsen et al. 2002). These results showed that ABA and BR signaling pathways act antagonistically in response to drought stress in B. napus. BR and ABA antagonize each other in many biological processes (Wang et al. 2014).
Developmental regulation during drought stress
The expression of the B. napus homologs of the transcription factor genes CONSTANS-like 4 (COL4) and AGAMOUS-LIKE 20 (AGL20) (also known as SOC1) were significantly downregulated under drought stress. The drought responsive regulatory network (Fig. 4) highlighted the induction of miR435 which potentially targets the expression of the B. napus COL4. In agreement with our results, the expression of a number of CONSTANS-like family members was downregulated by drought stress in rice (Shaik and Ramakrishna 2013). The expression of Arabidopsis COL4 was strongly induced by ABA, salt, and osmotic stress (Min et al. 2014). It has been demonstrated that the Arabidopsis SOC1, as a positive regulator of flowering, was upregulated following 4 to 5 days drought treatment indicating that plants reprogram flowering process to shorten the life cycle and promote survival under extreme environments (Su et al. 2013). The apparently contrary responses of these genes to drought may be due to different timing and intensity of stress, and even the existence of various strategies to set time to flowering under stress conditions in plants. Further investigations are required to determine the regulatory roles of B. napus COL4 and SOC1 under stress conditions.
In conclusion, this EST-based transcriptome study provided an insight into how B. napus genome responds to drought stress. The results clearly showed a substantial change in the B. napus genome expression in response to drought in which the corresponding regulatory mechanisms were significantly activated. We identified a number of significant drought responsive regulatory genes and miRNAs with various functions and characterized their interactions using the regulatory network. Our findings indicated that B. napus recruits a complex and overlapping regulatory system to modulate response to drought stress. The results highlighted the importance of biosynthesis and signaling of various plant hormones and their crosstalks to coordinate physiological events under drought stress in B. napus. The obtained drought responsive regulatory network can be utilized to select appropriate genes for transgenic breeding to develop B. napus drought tolerant cultivars.