Exploring transcription factors reveals crucial members and regulatory networks involved in different abiotic stresses in Brassica napus L.

Background Brassica napus (B. napus) encompasses diverse transcription factors (TFs), but thorough identification and characterization of TF families, as well as their transcriptional responsiveness to multifarious stresses are still not clear. Results Totally 2167 TFs belonging to five families were genome-widely identified in B. napus, including 518 BnAP2/EREBPs, 252 BnbZIPs, 721 BnMYBs, 398 BnNACs and 278 BnWRKYs, which contained some novel members in comparison with existing results. Sub-genome distributions of BnAP2/EREBPs and BnMYBs indicated that the two families might have suffered from duplication and divergence during evolution. Synteny analysis revealed strong co-linearity between B. napus and its two ancestors, although chromosomal rearrangements have occurred and 85 TFs were lost. About 7.6% and 9.4% TFs of the five families in B. napus were novel genes and conserved genes, which both showed preference on the C sub-genome. RNA-Seq revealed that more than 80% TFs were abiotic stress inducible and 315 crucial differentially expressed genes (DEGs) were screened out. Network analysis revealed that the 315 DEGs are highly co-expressed. The homologous gene network in A. thaliana revealed that a considerable amount of TFs could trigger the differential expression of targeted genes, resulting in a complex clustered network with clusters of genes responsible for targeted stress responsiveness. Conclusions We identified and characterized five TF families in B. napus. Some crucial members and regulatory networks involved in different abiotic stresses have been explored. The investigations deepen our understanding of TFs for stress tolerance in B. napus. Electronic supplementary material The online version of this article (10.1186/s12870-018-1417-z) contains supplementary material, which is available to authorized users.


MEME motif prediction for the identified TF families
A motif is a sequence pattern that occurs repeatedly in a group of sequences. Motifs in the MEME 7 Suite are represented as position-dependent letter-probability matrices that describe the probability of each possible letter at each position in the pattern. Individual MEME motifs do not contain gaps. Patterns with variable-length gaps are split by MEME into two or more separate motifs. The identified conserved motifs for each of the five TF families by using the MEME search tool 7 are shown in Figs.1-5. Sequence logos in each figure show the conserved residues in the associated family. The y − axis shows the information content 8 IC = 4.32 − H (bit), here, log 2 20 = 4.32 denotes the maximum entropy at a certain position of multiple sequences, which corresponds to uniformly distribution of the 20 amino acids at the position. During multiple sequences alignment, one can compute the frequency distribution of the 20 amino acids at each position, and obtain the entropy 8 H for the considered position according to Here, p i denotes the frequency of the i th amino acid, and ∑ 20 i=1 p i = 1. Each family is predicted with many conserved motifs. The bigger the amino acid letter, the higher the probability of the amino acid appeared at the associated position of protein sequences, and the corresponding amino acid is conserved.  Fig.6 shows that the identified BnAP2/EREBPs, BnbZIPs, BnNACs and BnWRKYs are all members of those in PlantTFDB. While, 50 of the identified 721 BnMYBs are not members of those as listed in the database. Moreover, for each family, we exclude several members in PlantTFDB or references in our work.

Literature and databases comparison analysis for the identified five TF families
Compared with 87 15 , 132 9 and 515 BnAP2/EREBPs 16 in B. napus that the references have reported, we identified 518 BnAP2/EREBPs, which include some novel members. For the 17 BnAP2/EREBPs that are not identified by us, we check their structures one by one. We find that 9 of the 17 genes (BnaC03g69940D, BnaC09g24990D, BnaA07g12190D,  12 , we identified 398 BnNACs, whose number is six times more than the existing references. WRKYs in various species have been extensively investigated 13 As an example, we show the highly conserved sections for the 278 WRKYs in Fig.7. Researchers have found that the WRKY amino acid sequence have been replaced by WRRY, WSKY, WKRY, WVKY, WKKY or WKNY in some WRKY proteins 23,24 . For the 278 BnWRKYs, 275 proteins are with the usual WRKY amino acid sequences. However, the WRKY amino acid sequences for three group IIc BnWRKYs have been replaced by WKKY or WKNY. Actually, the BnaCnng77260D has been replaced by WKNY, while BnaCnng28950D and BnaCnng39890D have been replaced by WKKY. The WRKY genes can be classified into three subfamilies 27 . Subfamily I WRKY proteins have two WRKY domains, the N-terminal (NT) domain and the C-terminal (CT) domain, as shown in Fig.7. Most of WRKY genes with one WRKY domain belong to subfamily II. Both subfamily I and II WRKY proteins have the same zinc motifs. WRKY proteins with only one WRKY domain but different patterns of zinc finger motifs are categorized into subfamily III. Subfamily III WRKY domains contain a C2-HC motif (C-X7-C-X23-H-X1-C). Subfamily II proteins can be further divided into IIa to IIe based on the primary amino acid sequences. Among the identified 278 BnWRKYs in B. napus, detailed analysis on their amino acid sequences reveal that 49 genes belong to subfamily I, 183 genes are subfamily II, while 46 genes are classified into subfamily III. For the four species, subfamily II WRKYs take up more than one half, especially, more than 25% WRKYs belong to group IIc.

Phylogenetic trees and subfamilies distributions of the identified five TF families
Phylogenetic tree analysis was performed using MEGA by neighbor-joining method. The unrooted phylogenetic trees for the five families are shown in Fig.8. In each panel, we show the cases for B. napus in comparison to A. thaliana, B. rapa, B. oleracea and B. napus. Subfamilies distributions of the identified five families of TFs are shown in Fig.9.
From the constructed phylogenetic trees and subfamilies distributions, we observe that TFs in each family can be classified into several different clades. Specifically, the BnAP2/EREBPs can be divided into five clades consisting of four subfamilies: AP2, ERF, DREB and RAV. The two super-families DREBs and ERFs in B. napus were further distributed into six subgroups, A1 to A6, B1 to B6, respectively. We found that more than 10% AP2/EREBPs in the three Brassica species belong to ERF B3, DREB A4 and ERF B1 respectively. However, more than 15% AtAP2/EREBPs in A. thaliana are classified into ERF b6. Specifically, for B. napus, we find that a total of 194 DREBs were distributed into groups A1-A6, containing 17, 30, 3, 66, 42 and 36 genes respectively. Additionally, 244 ERFs were distributed into groups B1-B6 containing 61, 18, 76, 21, 25, and 43 genes respectively, with 28 ERFs can not be assigned to any specific subgroups. Finally, 32 genes were classified into the AP2 subfamily and 20 genes are assigned into the RAV subfamily. The AtbZIPs in A. thaliana can be classified into 7 subfamilies, while the bZIPs in B. rapa, B. oleracea and B. napus can be divided into 9, 11 and 12 subfamilies. Almost 50% AtbZIPs (59) in A. thaliana belong to subfamily I, while only 12 of the 518 BnbZIPs are subfamily I in B. napus. In B. napus, subfamily III BnbZIPs is the biggest subfamily, which takes up about one quarter (63), while members in subfamilies XI and VI are both less than 5.

7/23
and 17 subfamilies. More than 20% BoNACs in B. oleracea belong to subfamilies X and IV respectively. The B. napus is most enriched in subfamily VI, while B. rapa is with more than 15% BrNACs in subfamily XIII. Among the identified 278 BnWRKYs in B. napus, detailed analysis on their amino acid sequences reveal that 49 genes belong to subfamily I, 183 genes are subfamily II, while 46 genes are classified into subfamily III. For the four species, subfamily II WRKYs takes up more than one half, where more than 25% WRKYs belong to group IIc. Interestingly, we have not found BnWRKY members in group IIa in B. napus. In fact, literature reported group IIa WRKYs are classified into group IIb in our work, and sequences of groups IIa and IIb proteins are with high similarity. The numbers of group specific WRKYs in each species can be found in Tab.3.

Heat maps for cis-acting elements of the five families
The five TF families are rich in various hormone or stress responsive cis-acting elements. In this paper, we consider 11 cis-acting elements, for details, see Tab.4. Heat maps of cis-acting elements for the five TF families can be found in Fig.10. Fig.10    stress response, ABA signaling pathways and heat stress response, respectively. The BnNACs are relatively rarely involved in W box, where no more than 40% genes are with W box. Chromosome distributions of TFs with various cis-acting elements show that different cis-acting elements prefer different chromosomes (Fig.?? (b) and Fig.4). For example, TC-rich repeats in BnAP2/EREBP are enriched in A03 and C07; Many elements in BnMYBs, BnNACs and BnWRKYs are enriched in A03 and C03. Elements in BnbZIPs are enriched in A09, C03 and C04. Our results also show that some chromosomes are always with rare cis-acting elements, such as C01.
Distributions of TFs with various cis-regulatory elements in each subgroup for each family show that different cis-acting elements prefer different subgroups (Fig.?? (c) Figure 4. Heat maps of cis-acting elements for the five TF families. Yellow represents with the associated regulatory elements, while green denotes no elements.

Synteny analysis
To investigate the evolutionary relatedness of the identified sequences together with AP2/ERF, bZIP, MYB, NAC, WRKY genes encoded by the fully sequenced Arabidopsis, we performed synteny analysis using the conserved AP2/ERF, bZIP, MYB, NAC, WRKY transcription factor domain, as shown in Fig.5. There are 66, 26, 54, 30, 24 single-copy, and 419, 204, 602, 339, 235 multi-copy AP2-EREBP, bZIP, MYB, NAC and WRKY genes in B. napus. Numbers of multi-copy genes are all more than six times higher than single-copy genes. Synteny analysis among the A subgenomes of B. rapa and B. napus, the C subgenomes of B. oleracea and B. napus show strong co-linearity, although chromosomal rearrangements have occurred between B. napus and B. rapa or B. oleracea. For example, 6 of the 66 single copied AP2/EREBP genes in B. napus locate on chrA08, genes on chrA08 are syntheny with genes on chromosomes A08 for B. rapa, C3 and C8 for B. oleracea. Similarly, the 25 multi-copied WRKY genes located on chrA03 in B. napus are synteny with chromosomes A03, C03, chrC07, chrC03, chrA03, chrA10 and so on.
Moreover, there may be several homologies of an Arabidopsis gene in each of the three Brassica species , which may be used to explain why numbers of TFs in the three Brassica species are higher than those in Arabidopsis, and why numbers of TFs in B. napus are higher than those in B. rapa and B. oleracea. At the same time, the above result also indicates the evolution of B. napus may incorporate with gene duplication events.

Expression profiles of five TF families under five stress conditions in B. napus
It is assumed that the five TF families play important roles during stress response. For B. napus, we consider five stress conditions including a: Cold 4 • C; b: heat treatment 40 • C; d: Drought treatment with 15% PEG; e: ABA treatment with 30um ABA; f: Salt treatment with 150mM NaCl. Figure 10. Heat maps of cis-acting elements for the five TF families. Yellow blocks denote that the corresponding genes have the associated cis-acting elements, while green blocks denote the corresponding genes do not have the corresponding cis-acting elements.

Synteny analysis of the five TF families in B. napus
Homologous genes were predicted by means of synteny-based methods in the CoGe database (https://genomevolution.org) 28 . Orthologous genes of B. napus were predicted in the target genomes of A. thaliana, B. rapa, and B. oleracea through the identification of a conserved syntenic region via SynFind with the following parameters: last comparison algorithm, minimum number of four genes, and unlimited syntenic depth. Putative paralogous genes in B. napus were identified via SynMap by the following parameters: last algorithm, 20 genes for the maximum distance between two genes, and five genes for the minimum number of aligned pairs. The algorithm for syntenic depth was defined as quota-align with a coverage depth of 1:1 for B. napus versus B. rapa and B. oleracea. The coverage depth of B. napus versus A. thaliana was determined at 6:1. Orthologous sets were visualized via McScanx (http://chibba.pgml.uga.edu/mcscan2/).
The identified five families of sequences in B. napus are evolutionary homologous with its diploid progenitors B. rapa and B. oleracea. Synteny analysis among the A sub-genomes of B. rapa and B. napus, the C sub-genomes of B. oleracea and B. napus show strong co-linearity, although chromosomal rearrangements and gene duplication events have occurred between B. napus and parental species after their divergence from the common ancestor (Additional file 3). Synteny analysis allows us to characterize four categories of genes: 1) Novel genes, which emerged during the evolution from parental species to B. napus. 2) Conserved genes, which have a unique copy from B. rapa, B. oleracea to B. napus. 3) Multi-copy genes. If two or more genes from the same species were detected in a homologous gene set of the four species 29 . 4) Lost genes 29 , which were missing during the evolution from parental genomes to B. napus. There are 33, 22, 65, 29, 19 novel, and 66, 26, 54, 30, 24 conserved BnAP2/EREBPs, BnbZIPs, BnMYBs, BnNACs and BnWRKYs among the identified TFs in B. napus respectively. While 419 BnAP2/EREBPs, 204 BnbZIPs, 602 BnMYBs, 339 BnNACs and 235 BnWRKYs are multi-copy genes. Novel genes averagely take up about 7.6% in the five families. BnMYBs posses the most novel genes (9.0%), followed by BnbZIPs (8.7%). BnAP2/EREBPs contain the lowest novel genes (6.4%), and BnWRKYs contain about 6.8% novel ones. About 9.4% TFs are conserved in B. napus in the five families, including 12.7% BnAP2/EREBPs, and about 7.5% BnMYBs and BnNACs. More than 80% TFs in B. napus are multi-copy ones. The BnNACs and BnWRKYs contain the most multi-copy genes (about 85%), while BnAP2/EREBPs and BnbZIPs contain about 80% multi-copy ones. For the A sub-genome, a total of 30 TFs lost during the evolution from B. rapa to B. napus, including 12 BrAP2/EREBPs, 7 BrbZIPs, 8 Figure 6. Heat maps of cis-acting elements for the five TF families. Yellow represents with the associated regulatory elements, while green denotes no elements.  . Synteny analysis for single-copy genes (UP) and multi-copy genes (Down) among Arabidopsis, B. rapa, B. oleracea and B. napus. The genes in the four species are plotted against their linked counterpart chromosomes. From left to right, we show the cases for AP2/EREBP, bZIP, MYB, NAC and WRKY. From outside to inside of each circos plot, the first layer shows the chromosome names, A* corresponds to B. rapa, chr* represnts B. napus; C* and S* represent B. oleracea and Scaffold, respectively; Chr* represents Arabidopsis. The second layer shows gene ID (partially); The third layer denotes length of chromosomes; The last layer links the gene IDs with synteny. Synteny block analysis reveals changes in transcription factors during the evolution from diploid to tetraploid.

Transcriptome sequencing data and results
We consider five treatments, including cold treatment at 4 • C low temperature, heat treatment at 40 • C high temperature, drought treatment with 15% PEG-6000, salt treatment with 150mM NaCl and ABA stimulus with 30 mmol/L ABA. The 7 days old 11/23 seedlings were treated under each condition lasting 12 hours for transcriptional sequencing. Three repeats are considered under each treatment. Totally 18 samples (5 conditions × 3 repeats) are sequenced using Illumina Hiseq platform, including 15 treated samples and 3 control samples without any treatments. On average we generated about 6.56Gb bases from each sample. Averagely 72.09% clean reads are mapped to reference genome. After mapping sequenced reads to reference genome and reconstructing transcripts, we finally identify 53948 genes on average for each sample, including 50815 known genes and 3133 novel genes on average for each sample. For details of the RNA-Seq results, see Tab. 5. Table 5. Summary of sequencing reads after filtering, genome mapping and expressed genes. RL: Read Length(bp); TMR: Total Mapping Ratio; UMR: Uniquely Mapping Ratio; Uniquely Mapping: Reads that map to only one location of reference, called uniquely mapping. Q20: the rate of bases which quality is greater than 20.

Enrichment analysis of cis-acting elements for the identified 315 crucial DEGs
To explore the enriched cis-acting elements for the selected DEGs, we use Zscore 30 to quantify the relative abundance of each cis-acting element, defined as: Here, F i j denotes the actual number of curated DEGs in the i th (i =BnAP2/EREBP, BnbZIP, BnMYB, BnNAC, BnWRKY) family with the j th ( j =ABRE, GARE motif, P box,TCA element,TGA element, TGACG motif, HSE, LTR, MBS, TC-rich repeats, W box) cis-acting element. R i j is a vector with 5000 elements, with each element denotes the number of DEGs with the j th element, where the 5000 rounds of randomly sampled amounts of DEGs equal to the amount of curated DEGs in the i th family. std(R i j ) represents the standard deviation of R i j . Zscore for each cis-acting element of each family is shown in Tab.6. Hereinafter, Zscore > 0 represents enriched element, while Zscore < 0 indicates underrepresented element. Table 6 shows that the identified 93 BnAP2/EREBPs are especially enriched with TGA element with Zscore = 1.9025, but are underrepresented with TCA element, TC-rich repeats and W-box. Actually, 45 of the 93 BnAP2/EREBPs are with TGA element, however, for randomly sampled 5000 sets of BnAP2/EREBPs with 93 DEGs, the average number of genes with TGA element is around 32. The 42 BnbZIPs are enriched with ABREs, but are underrepresented with TCA element. Actually, 29 of the 42 BnbZIPs are with ABREs, which is higher than the average number 21.7 for 5000 randomly sampled sets. Though 21 of the 42 BnbZIPs are with TCA element, as compared with the average number 26.5 for 5000 sets of randomly sampled genes, it is relatively lower for the 42 crucial BnbZIPs. The 94 BnMYBs are enriched with HSE, P box, GARE motif, TGA element, however, they are shorted in TC-rich repeats and W box. The 48 BnNACs are enriched with GARE motif and TGACG motif, but they are underrepresented with four cis-acting elements, including TGA element, HSE, MBS and TC-rich repeats. For the 38 BnWRKYs, results reveal that they are enriched with MBS and TC-rich repeats, but are shorted in ABRE and HSE.

Topological analysis on subnetworks of the constructed gene co-expression network
The constructed gene co-expression network consists of five families of genes, we withdraw subnetworks for each family and draw the corresponding co-expression networks, as shown in Fig.12. For the constructed co-expression networks, we perform detailed topological analysis. Table 7 shows the statistics information 31 on the constructed gene co-expression network and the 12/23 corresponding five subnetworks as shown in Fig.12. Specifically, average degree of undirected and unweighted network is defined as E is the number of edges, and N is the number of nodes in the complex network respectively. The link density of the network is calculated as: .
Here, 1 2 N(N − 1) represents all possible edges between any two of the N nodes in the complex network. Higher degree or higher link density indicates each node may have more neighbors.
If high degree nodes tend to connect with high degree ones, then the network is assortative; whereas, if high degree nodes tend to connect with low degree ones, then the network is disassortative. The degree-degree correlation coefficient DCC 32 can act as an indicator of assortativity and disassortativity, which is defined as Here, M is the total number of edges, k i , j i are the degrees of the nodes at the ends of the i th(i = 1, 2, ...M) edge. DCC < 0 indicates the disassortativity of the network, and DCC > 0 indicates the assortativity, while DCC = 0 indicates no degree correlations. The obtained statistical results for the co-expression networks are summarized in Tab.7. The average degrees and link densities of the co-expression networks are very low. The average degrees and link densities of the BnNAC and BnMYB subnetworks are higher than the other subnetworks, which indicate genes in the two TF families may have relatively stronger co-expression relationships. The co-expression network is assortative, indicating that highly co-expressed genes tend to connect with each other in the network.
Moreover, to evaluate the relative connection density within each TF family in comparison to those among families, we define the following index, called relative density (RD): Here, N maximum within represents all possible pairs of edges within each of the five TF families.
The RD reveals that the relative connection density among families is lower than that within families, and it indicates that TFs in the same family tends to exhibit similar expression profiles under various stresses or stimuli. Furthermore, we declare that the degrees of the 315 DEGs in the gene co-expression network have no apparent correlation with their expression levels. To verify this declaration, we compute the Spearman correlation coefficient between degree sequence and the expression levels under the five treatments. Scatter plots for each treatments can be found in Fig.13. The Spearman correlation coefficients between degree sequence and the expression levels under the five treatments are all lower than 0.17.   Topological analysis subnetworks for the constructed process-gene network    Figure 13. Scatter plots of node degree versus relative expression levels under cold, heat, drought, salt and ABA treatments for the gene co-expression network. r represents the corresponding Spearman correlation coefficient between degree sequence and the relative expression level.

14/23
Topological analysis on subnetworks of the constructed process-gene network The constructed process-gene network and its subnetworks are shown in Fig.14. The obtained statistical indexes are summarized in Tab.8. From Fig.14 and Tab.8, we conclude that the process-gene network and its subnetworks are all densely connected. Especially, the link density of the BnMYB subnetwork equals 0.9, which indicates the 315 TFs involve in some closely related biological processes. The BnMYB subnetwork and the BnWRKY subnetwork are disassortative, which indicate highly connected genes in the process-gene network tend to connect with low degree ones. In the process-gene network, highly connected genes correspond to those genes who involve some popular biological processes, while low degree genes may involve some unpopular biological processes. Thus, the disassortativity of the BnMYB subnetwork and the BnWRKY subnetwork further indicate that some genes in the two families may simultaneously involve some popular and unpopular biological processes.
For the process-gene network, similar to the co-expression network, we can calculate the index RD. For the process-gene network, N within = 7652, N among = 29466 − 7652 = 21814. N maximum within = 11341, N maximum among = 38114. Therefore, This also reveals that the relative connection density among families is lower than that within families, and it indicates genes in the same TF families prefer similar biological processes.  The scatter plots of unweighted node degrees versus relative expression levels under cold, heat, drought, salt and ABA treatments for the process-gene network are shown in Fig.15. Fig.15 also reveals that there are weak correlations between node degree and relative expression levels. DEGs with high relative expression levels are unnecessarily to involve many popular biological processes.

8/11
BnNAC BnWRKY Figure 9. Subnetworks for each of the five family for the process-gene network.  and ABA 1, 35, 36 , and we found its homologous genes in B. napus are ABA, cold, heat and salt responsiveness. Therefore, we speculate that the identified TFs can actually trigger the differentially expression of their down-stream genes.
On the other hand, we find 25, 36, 33, 121 and 13 up-stream genes that can regulate TFs of the five families respectively, and their homologous genes in B. napus are differentially expressed under at least one of the five treatments. There are hundreds of up-stream genes for AtNACs. Most of the up-stream genes are heat responsiveness. Except AtWRKYs, about 50% up-stream genes for each TF families are cold responsiveness. About 25% up-stream genes of AtAP2/EREBPs and AtbZIPs are salt responsiveness, and more than 30% up-stream AtAP2/EREBPs, AtMYBs and AtWRKYs are ABA responsiveness. However, only one up-stream gene for AtWRKYs is cold responsiveness. The above results indicate the identified TFs can trigger the differentially expression of tens of down-stream genes, and they can also be regulated by tens to hundreds of up-stream genes. To see whether the response of a gene could possibly be triggered by up-stream genes and whether the gene can trigger differentially expression of down-stream genes, we take two subnetworks for AtWRKYs as concrete examples, which are shown in Fig.17 (a). The first subnetwork consists of three AtWRKY genes 17/23  Table S7). AT1G80840, AT2G30250 and AT1G62300, and the second subnetwork consists of another three AtWRKY genes, including AT4G24240, AT2G23320 and AT4G31550.
In the first subnetwork, the homologous genes of the AtWRKY gene AT1G80840 are cold and heat responsive in B. napus, and these homologous genes are up-regulated under the two stresses. The homologous genes of AT2G30250 is heat inducible, and the homologous genes of AT1G62300 are cold and heat responsive. AT1G80840 can mutually regulate gene AT5G13630 and regulated by gene AT1G78300. The homologous genes for AT5G13630 all can respond to heat stress and down-regulated. Five of the seven homologous genes for the up-stream gene AT1G78300 are heat responsive, and two of the seven homologous genes are ABA inducible. Therefore, we speculate that the stress responsiveness of AT1G80840 is closely related to its up-stream and down-stream genes. AT1G78300 can regulate another AtWRKY gene AT1G62300, while one homologous gene of AT1G62300 can respond to cold and one can respond to heat stress. Homologous genes of another three up-stream genes for AT1G62300 are all heat responsiveness, and homologous gene of AT2G46020 in B. napus can further respond to drought stress, and AT5G38480 can further respond to ABA stimulus. Thus, we speculate that the cold and heat responsiveness of AT1G62300 is possible regulated by its up-stream genes. The third AtWRKY gene AT2G30250 can regulate two genes, and homologous genes for one of the two genes can respond to salt, ABA, and another one can respond to ABA and drought. Since drought and salt stresses can be induced by heat stress, thus, we speculate that the heat responsive AtWRKY AT2G30250 is responsible for the stress responsiveness of its two down-stream genes.
Similarly, in the second subnetwork, homologous genes of the three AtWRKY genes AT4G24240, AT2G23320 and AT4G31550 can respond to heat stress, and homologous genes of AT2G23320 can also respond to ABA stimulus. By analyzing the expression profiles of up-stream or down stream genes for the three AtWRKYs, we can also declare that the stress responsiveness of the three AtWRKYs correlate with their up-stream genes, and the three AtWRKYs are responsible for the stress responsiveness of their down-stream genes.