Abstract
Cabbages (Brassica oleracea L.) are an important vegetable crop around world, and cold temperature is among the most significant abiotic stresses causing agricultural losses, especially in cabbage crops. Plant bZIP transcription factors play diverse roles in biotic/abiotic stress responses. In this study, 119 putative BolbZIP transcription factors were identified using amino acid sequences from several bZIP domain consensus sequences. The BolbZIP members were classified into 63 categories based on amino acid sequence similarity and were also compared with BrbZIP and AtbZIP transcription factors. Based on this BolbZIP identification and classification, cold stress-responsive BolbZIP genes were screened in inbred lines, BN106 and BN107, using RNA sequencing data and qRT-PCR. The expression level of the 3 genes, Bol008071, Bol033132, and Bol042729, was significantly increased in BN107 under cold conditions and was unchanged in BN106. The upregulation of these genes in BN107, a cold-susceptible inbred line, suggests that they might be significant components in the cold response. Among three identified genes, Bol033132 has 97% sequence similarity to Bra020735, which was identified in a screen for cold-related genes in B. rapa and a protein containing N-rich regions in LCRs. The results obtained in this study provide valuable information for understanding the potential function of BolbZIP transcription factors in cold stress responses.
1. Introduction
Cabbage (Brassica oleracea L.) plants represent one of the major vegetable crops grown worldwide. Most crops of B. oleracea and its sister species Brassica rapa produce a range of phytochemicals with diverse functions for plant defense such as polyphenolic compounds, carotenoids, and glucosinolates [1, 2]. The draft genome sequences of B. oleracea (with the CC genome) and B. rapa (with the AA genome) were recently published [3, 4]. A total of 66.5% (34,237) of B. oleracea genes and 74.9% (34,324) of B. rapa genes were clustered. In total, 5,735 B. rapa-specific genes and 9,832 B. oleracea-specific genes among 45,758 protein coding genes were identified. The availability of published genome sequence for these crop plants facilitates studies of structural and functional genomics in agronomically important species.
Plant bZIP transcription factors play diverse roles in developmental and physiological processes and biotic/abiotic stress responses such as ABA signaling for osmotic stress responses during vegetative growth [5], seed germination and flowering time [6], glucose-ABA signaling [7], sugar signaling during metabolism [8], lipid stress responses [9], response to zinc deficiency [10], salicylic acid- (SA-) dependent plant systemic defense responses and the activation of jasmonic acid- (JA-) and ethylene (ET-) dependent defense mechanisms [11], anthocyanin accumulation during photo morphogenesis [12], floral patterning [13], auxin-mediated histone acetylation related AtbZIP11 [14], and ABA signaling related to stress tolerance [15]. As the focus of recent studies due to their importance as regulator of responses to the biotic and abiotic stresses, bZIP transcription factors have been identified in diverse plants. Based on the presence of the UARR and LCRs, 136 bZIPs were identified in B. rapa; 64 were found in cucumber based on predicted structural features, 92 in sorghum through genome-wide identification and characterization, 89 in rice according to their DNA binding specificity and amino acid sequences in basic and hinge regions, 131 in soybean based on the basic region of the bZIP domain and the presence of additional conserved motifs, 75 in Arabidopsis according to sequence similarities of their basic region and additional conserved motifs, and 141 in Hordeum vulgare [16–22]. However, little is known about the genome-wide survey and expression patterns of bZIP transcription factors in B. oleracea. Among the BolbZIPs, the function of only one gene related with drought stress and ABA has been reported. Expression of BolABI5 was dramatically induced by drought stress and exogenous ABA [23]. Heterogeneous expression of BolABI5 rescued the ABA-insensitive phenotype of the Arabidopsis abi5-1 mutant during seed germination, suggesting that BolABI5 likely functions in positive regulation of plant ABA responses.
The bZIP domain includes a basic region and a leucine zipper located on a contiguous α-helix. An N-x7-R/K motif comprising ~16 amino acids constitutes the basic region, which binds DNA containing a nuclear localization signal. The leucine zipper is composed of leucine residue repeat and is positioned precisely at nine amino acids towards the C-terminus from the arginine in the basic region, creating an amphipathic helix. To bind DNA, two subunits adhere via interactions between the hydrophobic sides of their helices, which create a superimposed coiled-coil structure for homo- or/and heterodimerization. Plant bZIPs preferentially bind to specific sequences, namely, the A-box (TACGTA), C-box (GACGTC), and G-box (CACGTG), but there are also examples of nonpalindromic binding sites [21].
In this study, we identified 119 BolbZIP proteins using the consensus sequence of several bZIP proteins and classified them based on specific amino acid sequence, unique amino acid repeat regions (UARRs), and low complexity regions (LCRs). Additionally, transcriptome analysis related to cold stress responses using RNA sequencing provided valuable information for research into stress tolerance and molecular breeding in B. oleracea.
2. Materials and Methods
2.1. Database Searches for bZIP Transcription Factors in B. oleracea
The AtbZIP, BrbZIP, and BolbZIP amino acid sequences obtained from TAIR (http://www.arabidopsis.org/), BRAD (http://brassicadb.org/brad/), and Bolbase (http://ocri-genomics.org/bolbase/). To confirm the presence of bZIP domain, UARR and LCRs in putative AtbZIP, and BrbZIP and BolbZIP proteins, the Motif scan tool (http://myhits.isb-sib.ch/cgi-bin/motif_scan), SMART tool (http://smart.embl-heidelberg.de/), and Batch CD-search tool (http://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) were used. bZIP proteins that showed the presence of a bZIP domain, UARR, and LCRs with confidence (-value < 0.1) in the Motif scan tool and Batch CD-search tool were used for further analyses. Next, LCRs were identified using the SMART tool.
2.2. Plant Material and Cold Treatment
Seeds of B. oleracea (inbred lines “BN106” and “BN107”) were germinated in soil and then grown for approximately 3 weeks in a growth chamber at 25°C under long day condition (16 h day/8 h night). For cold treatment, the 5-week-old plants were transferred to a 4°C growth chamber under continuous light conditions. The plants were then treated with cold temperature at 4°C for 6 h, followed by 0°C for 2 h. Further, the plants were subjected to freezing treatment at −2°C for 2 h followed by 4°C for 6 h.
2.3. RNA Extraction and cDNA Synthesis
Total RNA was isolated from plant tissues using an RNA extraction kit (Qiagen, USA) according to the manufacturer’s protocol. Total RNA was treated with RNase-free DNase (Promega, USA) to remove the genomic DNA contamination. The quality of total RNA was checked using a nanoDrop Spectrometer (nD-1000 Spectrophotometer, Peqlab) and agarose gel electrophoresis. cDNA was then synthesized using Superscript II reverse-transcriptase (Invitrogen), after which 5 μL (about 2 μg) total RNA and 1 μL of oligo dT (500 μg/mL) were mixed in the reaction tube and then heated at 65°C for 10 min. The enzyme was then added into the tube and incubated at 42°C for 50 min. Finally, the reaction tube was incubated at 70°C for 15 min to inactivate the enzyme.
2.4. RNA Sequencing
Two cabbage lines, BN106 and BN107 which exhibit different sensitivity to cold stress, were used for RNA sequencing. Total RNA was extracted from leaves of BN106 and BN107 at 2 h in 0°C. The total RNA was isolated using TRIzol reagent (Invitrogen, USA) following the manufacturer’s instructions. Total RNA (20 μg) from each sample, BN106_22°C and BN107_22°C (control) and BN106_0°C and BN107_0°C (treated), were used for Illumina sequencing (33 G 101 bp paired-end reads; Seeders, Republic of Korea). Transcripts of unigenes assembled from the total reads were validated by direct comparison with gene sequences in the Phytozome 15 (https://phytozome.jgi.doe.gov/pz/portal.html) using BLASTx (threshold -value ≤ ). The number of mapped clean reads for each unigene was counted and normalized using the DESeq package in R on two independent biological replicates. From the differentially expressed gene dataset, the transcripts of bZIP transcription factors were analyzed for up- and downregulated differentially expressed genes. BolbZIP sequence and RNAseq database sequences were aligned to each other using ClustalW with default parameters (http://www.genome.jp/tools/clustalw/).
2.5. RT-PCR and qRT-PCR
Quantitative real-time PCR (qRT-PCR) and reverse transcription PCR (RT-PCR) were conducted using cDNA from cold treated plants using primers specific for the BolbZIP gene (see Table S1 in Supplementary Material available online at http://dx.doi.org/10.1155/2016/4376598). RT-PCR was conducted using cDNA of plants exposed to cold and freezing temperatures (22°C, 4°C, 0°C, and −2°C). The PCR procedure involved predenaturation at 95°C for 5 min followed by cycles of denaturation at 95°C for 30 s, annealing at 60°C for 30 s, extension at 72°C for 30 min, and then a final extension for 5 min at 72°C. qRT-PCR was conducted by subjecting the samples to initial denaturation at 95°C for 10 min followed by 40 cycles of 95°C for 20 s, 60°C for 20 s, 72°C for 30 s, and final extension at 72°C for 2 min. An actin primer set for B. oleracea was used for normalization of RT-PCR and qRT-PCR.
3. Results
3.1. Identification of bZIP Transcription Factors in B. oleracea
To search for bZIP transcription factors in B. oleracea, we used the conserved bZIP domain consensus sequences (Table S2) of several proteins as BLASTP queries against the Brassica database (http://brassicadb.org/brad/). In addition, homology searches using 136 BrbZIP proteins were performed [16]. A total of 126 BolbZIP candidates were initially obtained with a probability -value threshold of 0.05. To confirm the presence of a bZIP domain in the selected bZIP proteins, domain searches were performed with several tools (see Section 2). After exclusion of the proteins lacking a bZIP domain, 119 putative BolbZIP transcription factors were identified. The position of each candidate BolbZIP gene in B. oleracea chromosome data available at Bolbase (Version 1.0) was then determined.
Among 119 candidate BolbZIP genes, 112 were mapped on chromosomes C01–C09 (Figure 1). 14 genes of BolbZIP were mapped on C01, 12 genes on C02, 15 genes on C03, 23 genes on C04, 8 genes on C05, 7 genes on C06, 10 genes on C07, 12 genes on C08, and 11 genes on C09. In particular, 20% of the BolbZIP genes mapped to chromosome 4 (Table S3). In addition, 7 genes were found in scaffolds that have yet been mapped to chromosomes. Bol024237 was anchored on Scaffold000093, Bol019052 on Scaffold000133, Bol016607 on Scaffold000153, Bol004200 on Scaffold000329, Bol003614 on Scaffold000345, Bol001886 on Scaffold000417, and Bol000879 on Scaffold000492.
3.2. Classification of BolbZIP Transcription Factors
We have classified the BolbZIP transcription factors based on amino acid sequence similarity to 136 BrbZIP and 75 AtbZIP proteins previously reported (Table 1) [16]. For the majority of bZIP proteins, we found orthologous groups including counterparts from each species, although occasionally no BrbZIP or AtbZIP homologs were found. AtbZIP and BrbZIP homologs of the BolbZIP proteins are summarized in Table 1. The proteins were divided into 63 categories based on the amino acid sequence similarity (Table 1). Most categories included two or three BolbZIP and BrbZIP proteins but a single AtbZIP. Analysis of the amino acid sequences revealed that the similarity between BolbZIP, BrbZIP, and AtbZIPs ranged from 50% to 90%. Several BolbZIP proteins showed over 90% similarity to the corresponding AtbZIP. For example, the similarity among Bol010308, At3g12250, and At5g06950 was 91–94%. For other genes, the closest homologs (with over 90% amino acid homology) were between the BolbZIP and the BrbZIP such as Bol004832 and Bra004689. BolbZIP proteins were also classified according to the method by Hwang et al. [16] based on UARRs and LCRs, which were further divided into 9 groups: glutamine (Q), aspartic acid (D), proline (P), asparagine (N), serine (S), glycine (G) rich domain, transmembrane (Tm) domain, LCRs only, and no LCRs except bZIP domain (Table 2, Tables S4 and S5). BolbZIP proteins and their orthologs from B. rapa and A. thaliana were found in the same groups. For example, BolbZIP of category 1 and its homologs Bra004550 and At2g46270 were classified into group 3A. LCRs of group 11 (only LCRs present) bZIP proteins composed single and mixed repeat natural amino acids. Group 12 contained bZIP proteins with no LCRs or specific amino acid-rich regions.
3.3. Candidate BolbZIP Genes for Responses to Cold Stress
To identify BolbZIP genes that might function in responses to cold stress, we carried out comparative analysis of the expression of BolbZIP gene in two B. oleracea inbred lines, cold-tolerant BN106 and cold-susceptible BN107. BolbZIP genes were selected from an RNA sequencing dataset based on their annotations and their expression profiles were analyzed (data not shown). Among the 119 BolbZIP genes, the expression of 41 genes was remarkably changed in responses to cold treatment, whereas 78 genes of them showed no significant changes in their expression. BolbZIP genes with significantly different expression were determined in 4°C-treated sample base on fold change (FC) ≥2 and ≤0.5 relative to 22°C-treated sample. Cold treatment at this temperature caused the upregulation of 18 genes in BN106 and of 7 genes in BN107, whereas 15 genes were downregulated in BN106 and 8 genes were in BN107 by cold treatment. In total, the expression of 21 genes was upregulated and 20 genes downregulated by cold treatment (Table 3). In addition, 6 genes were not showing any expression within BN106 lines and therefore not calculated (Table 3). Finally, 47 BolbZIP genes’ expression level was confirmed using quantitative real-time PCR (qRT-PCR) (Table 3). To obtain detailed expression for the putative cold-response BolbZIP genes thus identified, qRT-PCR was carried out using samples from plants treated at several temperatures (22°C, 4°C, 0°C, or −2°C). Totally, 25 BolbZIP genes with significantly different expression were selected based on fold-changes (FC) ≥3 and ≤0.5 relative to the control sample (22°C). Most of the tested genes were significantly upregulated by cold treatment except Bol021255. Among 25 tested genes, 22 genes are displayed in Figure 2 and three genes by RT-PCR in Figure 3. We were not able to determine the analogous relative expression for the latter three genes because they were not expressed in the 22°C treated sample. The expression levels of several BolbZIP genes were comparable between the two lines. However, no significant change in the expression of Bol008071, Bol033132, and Bol042729 was observed in response to cold treatment in BN106, whereas these genes were upregulated at all temperatures in BN107 (Figure 2(a)). By contrast, Bol009713, Bol013712, Bol016432, and Bol022925 were upregulated in BN106, but not in BN107 (Figure 2(b)). The increased expression of 17 BolbZIP genes was more pronounced after severe cold treatment at 4°C, 0°C, and −2°C (Figure 2(c)) and one gene was downregulated by cold treatment in both BN106 and BN107 (Figure 2(d)). Homologs of cold stress-response BrbZIP genes were included in the qRT-PCR [16]. These expression patterns are summarized in Figure 4. Moreover, several genes including Bol016432, Bol022925, Bol026864, Bol027732, and Bol028975 displayed differential expression between cold (4°C) and freezing (−2°C) temperature. The expression level of the 3 genes, Bol008071, Bol033132, and Bol042729, was significantly increased in BN107 under cold conditions and was unchanged in BN106. Among three genes, Bol033132 has 97% sequence similarity to Bra020735 which was previously reported gene. Two proteins, Bol033132 and Bra020735, contained N-rich regions in LCRs (Figure 5(a)). Moreover, Bol042729 included the N-containing LCR (Figure 5(b)). We suggest the possibility that BolbZIP proteins as well as BrbZIP proteins containing N-rich regions might be involved in cold stress response.
(a)
(b)
4. Discussion
It was known that B. rapa and B. oleracea genomes are highly similar in their gene structure, but there still exist species-specific genes in two species. Hence this study was carried out in B. oleracea and identified 119 BolbZIP proteins and placed them into 63 categories according to sequence similarity (Table 1). To identify the bZIP proteins in B. oleracea, a few bZIP domain consensus sequences of several species were used (Table S2). It is possible that this approach could lead us to underestimate the number of bZIP proteins present, despite the high number of BolbZIP proteins we identified. To address this, other search methods or more detailed consensus sequences for bZIP proteins in plants could be examined. In Arabidopsis, bZIP proteins were classified into different groups and subfamilies according to sequence similarities in their basic region and additional conserved motifs in order to elucidate the likely function of the proteins [21]. In rice, Nijhawan et al. [19] published 89 bZIP transcription factor-encoding genes based on DNA binding specificity and amino acid sequences in basic and hinge regions. Recently BrbZIP and AtbZIP proteins were divided into 9 groups based on their UARR and LCRs, which are highly enriched in one or a few amino acids [16]. In this study, 119 BolbZIP proteins were categorized into 63 groups and also classified according to UARR and LCRs based on the classification method of Hwang et al. [16]. In addition, the sequence similarity of the bZIP proteins of B. oleracea, B. rapa, and A. thaliana was analyzed. Most of homologs were found to have the same UARR and LCRs. UARRs were composed of 6 amino acids including Q, D, P, N, S, and G in the B. oleracea (Tables 2 and S4). This conservation of amino acid composition suggests that these 6 amino acids are important for biological functions and formation of protein structures in bZIP proteins.
BolbZIP gene family members were physically mapped to all the nine chromosomes of B. oleracea. Among them, chromosome 04 was found to contain the highest number of BolbZIP genes (21%), while chromosomes 05 and 06 harbored the fewest (6-7%) (Figure 1, Table S3). In B. rapa, the highest number of BrbZIP genes was detected in chromosome 09 (21%) [16]. Additionally, most BolbZIP genes were distributed in the arm end of each chromosome. The clustered distribution pattern of the BolbZIP genes on some chromosomes might be indicated in significant regions evolutionarily. For example, BolbZIP genes located on chromosomes 01, 02, 04, 07, and 08, and chromosomes 09 appear to be clustered at the arm end in those chromosomes (Figure 1).
To screen for cold stress-responsive BolbZIP genes, we tested the transcription patterns of BolbZIP genes enhanced or decreased by cold treatment in two B. oleracea lines that showed different cold tolerance [16]. Based on their expression patterns, the cold-responsive BolbZIP transcription factors were divided into four groups (Figure 2). We found that the expression of three genes, Bol008071, Bol033132, and Bol042729, was upregulated in cold-susceptible BN107 but not changed in cold-tolerant BN106. Additionally, when compared with 6 genes published for significant BrbZIP factors involved in the cold response, 4 BolbZIP genes (Bol004832, homologous to Bra000256, Bra004689, and Bra003320; Bol033132, homologous to Bra020735; Bol018688, homologous to Bra011648; and Bol021255, homologous to Bra023540) showed similar patterns of expression in response to cold treatment. For example, Bol033132 showed an expression pattern like that of its homolog Bra020735, indicating that these genes might be conserved key regulator in cold stress responses. Moreover, Bol033132 and Bol042729 encode bZIP proteins that include the LCR containing amino acid N or N-rich region (Figure 5, Tables S4 and S5). These results indicated that the N-containing region of BolbZIP proteins might be involved in cold stress responses. Although the functions of the N-containing region are largely unknown, the regions might be biologically active [24, 25]. This genome-wide identification and expression profiling of bZIP proteins from B. oleracea provides new opportunities for functional analyses, which may be used in further studies for improving stress tolerance in plants.
Competing Interests
The authors declare that there are no competing interests regarding the publication of this paper.
Acknowledgments
This research was supported by Golden Seed Project (Center for Horticultural Seed Development, no. 213003-04-4-SB110), Ministry of Agriculture, Food and Rural Affairs (MAFRA), Ministry of Oceans and Fisheries (MOF), Rural Development Administration (RDA), and Korea Forest Service (KFS).
Supplementary Materials
The supplementary materials contain 5 files, they are some important data supported to the methods and results of the presented study. These data make paper easier to read and understand.