Genome-Wide Analysis of Basic Helix–Loop–Helix Superfamily Members Reveals Organization and Chilling-Responsive Patterns in Cabbage (Brassica oleracea var. capitata L.)

Basic helix–loop–helix (bHLH) transcription factor (TF) family is commonly found in eukaryotes, which is one of the largest families of regulator proteins. It plays an important role in plant growth and development, as well as various biotic and abiotic stresses. However, a comprehensive analysis of the bHLH family has not been reported in Brassica oleracea. In this study, we systematically describe the BobHLHs in the phylogenetic relationships, expression patterns in different organs/tissues, and in response to chilling stress, and gene and protein characteristics. A total of 234 BobHLH genes were identified in the B. oleracea genome and were further clustered into twenty-three subfamilies based on the phylogenetic analyses. A large number of BobHLH genes were unevenly located on nine chromosomes of B. oleracea. Analysis of RNA-Seq expression profiles revealed that 21 BobHLH genes exhibited organ/tissue-specific expression. Additionally, the expression of six BobHLHs (BobHLH003, -048, -059, -093, -109, and -148) were significantly down-regulated in chilling-sensitive cabbage (CS-D9) and chilling-tolerant cabbage (CT-923). At 24 h chilling stress, BobHLH054 was significantly down-regulated and up-regulated in chilling-treated CS-D9 and CT-923. Conserved motif characterization and exon/intron structural patterns showed that BobHLH genes had similar structures in the same subfamily. This study provides a comprehensive analysis of BobHLH genes and reveals several candidate genes involved in chilling tolerance of B. oleracea, which may be helpful to clarify the roles of bHLH family members and understand the regulatory mechanisms of BobHLH genes in response to the chilling stress of cabbage.


Introduction
Brassica oleracea is widely grown all over the world, including cabbage (B. oleracea var. capitata), kale (B. oleracea var. acephala), Brussels sprout (B. oleracea var. gemmifera), cauliflower (B. oleracea var. botrytis), and broccoli (B. oleracea var. italica). B. oleracea is one of the most important leafy vegetable crops for studying morphological diversity and polyploidy evolution because of its unique U's triangle theory and the Brassica lineage-specific whole-genome triplication (WGT) event [1]. The available genome sequences of B. oleracea, B. rapa, and Arabidopsis thaliana provided valuable information for plant improvement and biological research [1,2]. A large number of studies indicate that the regulations of plant growth and stress responses are coordinated by the combination of different network regulators, such as transcription factors (TFs). TFs can specifically recognize and bind the gene promoters in the nucleus, regulating the temporal and spatial expression of genes [3][4][5]. TFs respond to plant physiological, biochemical, and stress conditions by activating or inhibiting the associated downstream genes [5][6][7]. The basic helix-loop-helix (bHLH) TF is one of the largest families of TFs in plants [4]. The bHLH regions contain two functionally distinct regions, including the basic region and the HLH region [8,9]. The N-terminal is a basic region [10], which is mainly associated with the binding ability of the E-box (CANNTG) DNA cis-acting element. The C-terminal is a helix-loop-helix region, which includes two amphipathic α helices and one ring of variable length hydrophobic amino acids that are connected to form homodimers or heterodimers and to interact with other proteins. Nowadays, plenty of plant genomes are sequenced, a large number of bHLH proteins have been identified [4,6,7,[11][12][13][14][15][16]. The bHLH proteins were divided into six major groups (A to F) through the phylogenetic relationship and specific protein motifs [17][18][19][20].
Chilling stress has a great effect on the metabolism and transcriptomes of plants. Chilling directly inhibits the enzymes of metabolism and reprograms the expression of genes [38]. ICE1 (inducer of CBF3 expression1)-CBF (CRT binding factor3)-COR (cold-regulated) chilling response pathway is one of the main signal transduction pathways for plant response and resistant to chilling stress. In A. thaliana, AtICE1 encodes a bHLH TF, AtICE1 can activate AtCBF3 transcription by combining with MYC sites on the promoter of AtCBF3 and regulate the response to chilling stress in A. thaliana [39][40][41]. ICE1, a CBF expression inducer, is an upstream regulator of chilling-responsive genes in various plants. ICE1-CBF pathway in diverse plants is conserved [42]. ICE1 can integrate different signals to adjust the chilling resistance of A. thaliana [43,44]. It has been shown that the protein kinase OST1 can phosphorylate and stabilize ICE1 under chilling stress, and increase the expression of the downstream chilling-tolerant genes [45]. However, mitogen-activated protein kinases (MPK3 and PMK6) reduce ICE1 stability and transcriptional activity by phosphorylate ICE1 protein, thus, negatively regulating the expression of CBF and chilling-tolerance of plants [40]. The transgenic analysis of different plants showed that the genetic engineering CBF pathway could significantly improve the chilling tolerance of plants. Overexpression of ICE1 in cucumber [46] and rice [47] has been shown to enhance the chilling tolerance of the plants. It has been found that the MdCIbHLH1 genes can be used to improve the chilling tolerance of the apple [48]. The bHLH family had been broadly studied in various plants. However, no systematic investigations on the bHLH TF family in B. oleracea have been reported. In this study, we comprehensively described the bHLH TFs in B. oleracea by comparative genomic analysis.

Phylogenetic Analysis and Genomic Location of BobHLH Genes
Complete amino acid sequences of bHLH proteins were screened by the Pfam databases and identified the bHLH domains. MEGA 6.0 software [51] was used to construct neighbor-joining (NJ) distance trees with the bHLH proteins of B. oleracea (234 BobHLH proteins), A. thaliana (162 AtbHLH proteins) and B. rapa (230 BrbHLH proteins), the bootstrap values set 1000 replicates, which provided the reliability of the statistic. To localize the BobHLH genes on the chromosomes, the position information (starting and ending positions) of BobHLHs were collected from the B. oleracea genome, and the MapChart [52] software was used to show BobHLH genes on nine chromosomes.

The Expression Patterns of BobHLH Genes in Different Organs/tissues
To analyze the expression profile of BobHLH genes in different organs/tissues, RNA-Seq was downloaded from the NCBI GEO database (GSE42891), which contains the expression levels in seven organs/tissues (callus, root, stem, leaf, bud, flower, and silique). The expression abundance of BobHLH genes in different organs/tissues were calculated using the fragments per kilobase of transcript per million fragments mapped (FPKM) values. The hierarchical clustering heat map of the BobHLH genes was generated by using the pheatmap package (https://cran.r-project.org/web/packages/pheatmap/) based on the log 2 (FPKM +1) values.

Plant Materials, Growth Conditions, and Chilling Treat
The seeds of chilling-sensitive cabbage (CS-D9) and chilling-tolerant cabbage (CT-923) were germinated and grown in pots containing sterilized soil, the environmental conditions were set: 14 h light/10 h dark, 25 • C day/18 • C night. When the seedlings grew up to five true leaves, the seedlings were selected and transferred to the vernalization chambers for 6 h and 24 h under 4 • C chilling stress, while the non-treated seedlings (mock) were still under normal conditions. The young leaves from the chilling treated and non-treated seedlings were collected at 6 h and 24 h, respectively, three biological replicates were set at each time point, and immediately frozen in liquid nitrogen and used for RNA-Seq analysis. Finally, a total of 24 RNA-Seq libraries were constructed and performed on the Illumina HiSeq TM 2500 platform.

The Characteristics Analysis, Conserved Motifs, and Gene Structure of the BobHLH Genes
Chilling-responsive BobHLH genes were selected to further study the characteristics analysis, conserved motifs, and gene structure. The theoretical isoelectric point (pI), amino acid number, molecular weight (MW), and hydrophobicity analysis of BobHLH proteins were analyzed by the ProtParam tool (https://web.expasy.org/protparam/). The conserved motifs of BobHLH protein sequences were identified by the Multiple Expectation-maximization for Motif Elicitation program (MEME, http://meme-suite.org/tools/meme) [53]. The parameters of the optimum motif width were set to 6 to 50 amino acids and the maximum motif number was set to 10. The Gene Structure Display Server (GSDS, http://gsds.cbi.pku.edu.cn/) [54,55] was used to draw the gene structure diagram according to the coding sequences and the corresponding genomic sequences of BobHLH genes.

Characterization of Putative Cis-Acting Elements in the Promoter Regions of BobHLH Genes
The 2000 bp upstream of the translation start codon of chilling-responsive BobHLH genes were obtained from the B. oleracea genome. The cis-acting elements in the 2000 bp upstream regions of the chilling-responsive BobHLH genes were predicted by using the PlantCare database (http: //bioinformatics.psb.ugent.be/webtools/plantcare/html/) [56].

Construction and Localization of Subcellular Localization Vector
WoLF PSORT (http://wolfpsort.org/) [57] and BacelLo (http://gpcr.biocomp.unibo.it/ bacello/pred.htm) were used to predict the sub-cellular localization of BobHLH proteins. In order to investigate the sub-cellular localization of BobHLH proteins, we used the transient transformed GPF fusion expression vectors to transform the epidermal cells of tobacco (Nicotiana benthamiana) to observe the position of the fluorescent protein. Forward primers with Sca I restriction sites and reverse primers with Xba I restriction sites were designed, respectively (the specific primers were in the Table S10), cDNAs were used as templates for PCR amplifications, to obtain the full-length coding sequences of the BobHLH054, BobHLH119, and BobHLH134 genes. The PCR amplifications were digested with Sca I and Xba I, then the fragments were ligated to the pCAMBIA2300-GFP vector. The constructed vector, as well as the empty vector (control), were transformed into Agrobacterium tumefaciens strain GV3101. A. tumefaciens were infiltrated into the tobacco leaves. After two days of infiltration, the GFP fluorescence was observed using a fluorescence microscope (OLYMPUS).

Identification and Phylogenetic Analysis of BobHLH Genes in B. oleracea
For the genome-wide identification of BobHLH genes, 234 available protein sequences were identified from the B. oleracea genome. It was confirmed that the 234 BobHLH proteins contain HLH domain. The 234 BobHLH genes were named from BobHLH001 to BobHLH234 based on their chromosomal positions (Table S4).
To study the classification and evolutionary relationships of BobHLH proteins, and to gain the potential function of BobHLH proteins for further investigation, the protein sequences of B. oleracea, A. thaliana, and B. rapa bHLH were used to generate a neighbor-joining (NJ) phylogenetic tree ( Figure 1). AT4G49770 (AtbHLH95), AT4G38071 (AtbHLH131), and AT2G20095 (AtbHLH133) were not found in A. thaliana. According to the clade support values and the classification of A. thaliana [11], the 234 BobHLH proteins were clustered into twenty-three subfamilies, while most of these subfamilies are common and have been reported in the phylogenetic trees of other species [6,12,13]. The results showed that the number of Ia subfamily member was the biggest and it contained 23 BobHLH proteins ( Figure 1). III(d+e) subfamilies contained only two BobHLH proteins (BobHLH151 and BobHLH178). There were twenty-five BobHLH proteins that differed highly from the other subfamilies members and these were classified as "orphans."

Chromosomal Localization of the BobHLH Genes in Cabbage
Based on the genomic position information, only 187 of the 234 BobHLH genes were unevenly distributed onto the nine chromosomes. The exact position of each BobHLH gene on chromosome (C01-C09) is shown in Figure 2. The rest of the 47 BobHLH genes were not anchored onto any of the B. oleracea chromosomes. Overall, the distribution of BobHLH genes was relatively dispersed, but some gene clusters also had relatively higher accumulation. Most of BobHLH genes were found on chromosome C03 (27,14.4%) and chromosome C08 (24, 12.8%), while chromosome C06 (7.5%), chromosome C05 (8.0%), chromosome C09 (9.1%) had 14, 15, and 17 members, respectively. Several BobHLH genes were densely distributed on chromosomes regions, several chromosomes regions had no BobHLH genes distribution ( Figure 2).

Expression Patterns of BobHLH Genes in Different Organs/Tissues
To investigate the expression pattern of BobHLH genes, we used the FPKM values normalized data from RNA-Seq (GSE42891) and evaluated their expression patterns among seven organs/tissues. As showed in Figure 3 (and Table S5), the expression profiles of BobHLHs in callus, root, stem, leaf, bud, flower, and silique had significant differences. PA total of 197 BobHLHs expressed in more than two organs/tissues, the relative transcript abundance of three BobHLHs (BobHLH081, BobHLH177, BobHLH218) were 0, and 13 BobHLHs did not have transcriptional abundance data. Moreover, some BobHLHs showed unique tissue expression patterns, we found seven BobHLHs (BobHLH022, -061, -076, -082, -083, -191, -228) expressed only in the callus, BobHLH154 was expressed only in the root, four BobHLHs (BobHLH014, -055, -152, -222) were specifically expressed in the leaves, six BobHLHs (BobHLH 058, -089, -141, -157, -158, -219) were expressed in the bud, and only BobHLH011 was expressed in the flower. Two BobHLHs (BobHLH062, BobHLH078) were specifically expressed in silique. These results suggested that the specifically expressed BobHLH genes had special functions in specific organs/tissues.

Expression Analysis of BobHLH Genes in Response to Chilling Stress
We also analyzed the expression patterns of BobHLH genes in CS-D9 and CT-923 under chilling stress (4 • C) and observed variation between the two different cultivars under different chilling stress (chilling for 6 h and 24 h; Figure 4; Table S6). A total of 181 BobHLHs were detectable in both of the two cultivars under chilling stress. There were several BobHLHs that have undetectable relative transcript abundance across the two cultivars under chilling stress (6 h and 24 h) and the mock-treated plants. When compared with mock-treated plants, three BobHLHs (BobHLH033, -156, and -221) were significantly up-regulated and seven BobHLHs (BobHLH003, -048, -059, 093, -109, -129, and -148 Table S6). In CT-923, BobHLH126 was significantly up-regulated in chilling-treated plants at 6 h under chilling stress, but the expression level was significantly down-regulated at 24 h. At the 24 h chilling stress, BobHLH054 was significantly down-regulated and up-regulated in chilling-treated CS-D9 and CT-923.

Characteristics, Protein Motifs, and Gene Structures of BobHLH Genes
Sequence analysis showed (Table 1; Table S7) that the length of the BobHLH proteins varied widely. The molecular weight ranged from 10.25 (BobHLH112) to 109.54 kDa (BobHLH008). The protein length of BobHLH proteins ranged from 90 (BobHLH112) to 994 (BobHLH008) amino acids and the theoretical isoelectric point (pI) ranged from 4.49 (BobHLH185) to 10.26 (BobHLH188), a total of 91 BobHLH proteins was more than 7, and they were in the alkaline range, the BobHLH proteins were rich in basic amino acid. The average instability index of all BobHLH proteins was 56.93 (> 40.00), indicating that most BobHLH proteins were unstable. The analysis of the grand average of hydropathicity (GRAVY) showed that the value for BobHLH proteins was negative, which indicated that almost all the BobHLH proteins belonged to hydrophilic proteins. To further analyze the diversification of BobHLH protein structure, the MEME program was used to predict the conserved motifs of BobHLH proteins ( Figure 5, Figure S1). A total of ten conserved motifs were characterized by chilling-responsive BobHLH proteins ( Figure 5), those of others were in Figure S1, and their logos were also generated ( Figure S2). The species and numbers of these conserved motifs tended to be consistent with our phylogenetic tree results ( Figure S1), several BobHLH proteins in the same subfamilies shared similar motifs, which means they might have similar biological functions [58]. Motif 1 and motif 2 were conserved in each BobHLH protein, except for BobHLH050 (only contains motif 1) and BobHLH143 (only contains motif 2; Figure S1). There were no motif 8 and motif 9 in the chilling-responsive BobHLHs ( Figure 5). Exon-intron structural diversity is considered to play an important role in the evolution of bHLH genes [59]. We analyzed the genomic DNA sequences and the corresponding coding sequences of 234 BobHLH genes to gain their exon-intron organization information ( Figure S3). While only the gene structures of chilling-responsive BobHLH genes were shown in Figure 6. The others were shown in Figure S3. The same subfamilies shared a more similar gene structure including the number and length of exons and introns ( Figure S3). The results indicated that the close evolutionary relationship in the same subfamilies [60]. Among the BobHLH genes, the exon numbers ranged from one to eleven. A total of twenty-five BobHLH genes had no intron, which accounted for 9.7% of the total BobHLH genes. The exon number of subfamily III(a+c) ranged from four to ten, while in subfamily Ia, it ranged from one to nine ( Figure S3). The gain or loss of these exons may lead to the functional diversity of BobHLH genes.

Subcellular Localization Analysis of the BobHLH Proteins
WoLFPSORT was used to predict subcellular localization, the results showed that most of BobHLH proteins are present in the nucleus (except for the nineteen BobHLH proteins), while BaceLO predicted that 219 BobHLH proteins exist in the nucleus (Table S9), these predictions indicated that BobHLH054, BobHLH119, and BobHLH134 proteins are distributed in the nucleus. We observed that the fluorescence signal of the GFP signal was present on both cell membrane and nucleus in tobacco leaf epidermal cells (Figure 8), while the green fluorescent signal in epidermal cells transferred with BobHLH054-GFP, BobHLH119-GFP, and BobHLH134-GFP was concentrated in the nucleus, indicating that the three BobHLHs are nuclear proteins, this is consistent with the predicted results.

Discussion
The rapid development of plant genome sequencing has accelerated the studies of genomics and functional genomics, making it more convenient to identify and explore important genes on the whole genome level, while comparative genomic analyses have become more significant. The whole genome information is also used to explore the mechanisms of plant growth and development, the response to biotic and abiotic stresses. B. oleracea belongs to Brassica species and it is one of the most important leafy vegetables in the world. The completion of B. oleracea genome sequencing has given researchers more opportunities to characterize more gene families and to identify their phylogenetic relationships. bHLH TF is a large class of plant TF family, which play important roles in plant physiological metabolism and response to various stresses.
In this study, we identified 234 bHLH genes in B. oleracea genome. In recent studies, the number of bHLH genes was varied; Chinese cabbage contained 230 BrbHLHs [12], 188 MdbHLHs had been identified in the apple "Golden Delicious" [64], and 113 FvbHLHs were detected in woodland strawberry [65]. The bHLH proteins classifications of several plants had been reported but the exact number is unknown. A. thaliana and rice are divided into twenty-six subfamilies and Chinese cabbage int twenty-four subfamilies [12]; according to the A. thaliana bHLH group nomenclature, bHLH proteins were separated into twenty-three subfamily in B. oleracea (Figure 1). These subfamilies were common in most species, suggesting that the bHLH proteins in conserved subfamilies might play an important role in plant evolution, while some non-conserved bHLH subfamilies may have evolved from plant-specific development or stress resistance [66]. Combined with the results of chromosome localization of BobHLH genes (Figure 4), several adjacent BobHLHs were clustered together on the same chromosome and they belonged to the same subfamily in the phylogenetic tree. BobHLH013, BobHLH014, and BobHLH015 belonged to Ia (Figure 1) and they were clustered into chromosome C01 (Figure 4) while BobHLH076 and BobHLH077 belonged to Ib(2) and were clustered into chromosome C04. These results can also be observed in other plants [66], these clustered genes may share similar conserved molecular functions. The conserved motifs ( Figure 5; Figure S1) and gene structures ( Figure 6; Figure S3) of these clustered BobHLHs were similar. The results of predicted cis-acting elements of BobHLHs showed the Ia subfamily member BobHLH013 contained LTR, while BobHLH014 and BobHLH015 did not. Usually, these clustered genes have some overlapping functions which might make them partially or totally redundant and lost [11].
According to the previous reports, bHLH genes were involved in stress responses. In general, the members of the same subfamilies may have similar molecular functions because of the specific conserved motifs, which could be used to predict the functions of unknown proteins. B. oleracea, A. thaliana, and B. rapa belong to Brassica and have close relationships, and the homologous genes among them were used to predict the functions in B. oleracea. It has been reported that BrICE1 encodes a protein with the bHLH domain, which accumulates rapidly in response to chilling stress in non-heading Chinese cabbage [67]. BrbHLH212 was expressed highly under chilling stress (12 h) in Chinese cabbage [12]. In this study, BobHLH093 and BrbHLH212 were clustered within subfamily VII(a+b) but it was significantly down-regulated expressed under chilling treat (6 h and 24 h), which may be due to its interaction with other bHLH proteins or non-bHLH proteins [11]. The homodimers or heterodimers formed between bHLH proteins or between bHLH and non-bHLH proteins determined the regulatory functions of bHLH proteins [62,68]. In the study of A. thaliana, ICE1 and ICE2 encode MYC-type bHLH TFs, and their orthologs were found to cluster into the same clade of subfamily VII(a+b) and III(d+e) with two corresponding members in B. oleracea (BobHLH138 and BobHLH151), respectively. They increased the chilling tolerance by inducing CBF/DREB1 regulon and regulating stomatal formation [69], chilling-induction activated its post-translational modifications, including ubiquitination, phosphorylation, and sumoylation [43,70], so that the whole plant showed chilling resistance. The results of cis-acting element analysis showed that there was a number of defense-related cis-acting elements in the promoter region of AtICE2, which indicated that AtICE2 might play an important role in plant biotic stress defense [71]. The promoter region of BobHLH054 contained LTR cis-acting elements, which were involved in the chilling stress response, and it was significantly down-regulated in CS-D9, while it was significantly up-regulated in CT-923 at 24 h under chilling stress, indicating that an important role of BobHLH054 gene in enhancing resistance to chilling stress. In the results of RNA-Seq, BobHLH156 and BobHLH003, -048, -059, -093, -109, and -148 were significantly upand down-regulated, respectively. These BobHLH genes showed similar expression patterns. Genes with similar expression trends might interact at the protein level [72] and coordinate the regulation of downstream genes [4]. These mechanisms depend on the interaction of these genes and it is possible that the expression of one bHLH is affected by or regulated by another one, resulting in similar changes. Similar changes in their expression requires further study to explore the interrelationships between these genes.

Conclusions
In this study, a total of 234 BobHLH genes were identified in the B. oleracea genome and were clustered into twenty-three subfamilies. Among them, 187 BobHLH genes were mapped to nine chromosomes, the rest were located on the scaffolds. Organ/tissue expression pattern analysis showed the organ/tissue-specific expression of BobHLH genes. RNA-Seq data analysis of chilling treat in B. oleracea indicated that compared with mock-treated plants, BobHLH054 was significantly down-regulated in CS-D9 and up-regulated in CT-923 at 24 h under chilling stress, which means that it may be involved in the chilling-response. The molecular weight of the BobHLH proteins ranged from 10.25 to 109.54 kDa. The characterization of the proteins and genes, including conserved motifs and gene structures, were similar in the same subfamilies. The stress-responsive cis-acting elements analysis showed that several of chilling-responsive BobHLHs contained LTR in the upstream regions. These findings in this study provide a foundation to investigate the information of BobHLH genes including the candidate genes of B. oleracea against chilling stress for further functional studies.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/11/914/s1: Table S1. The genomic sequences of the 234 BobHLH genes; Table S2. The coding sequences of the 234 BobHLH genes; Table S3. The protein sequences of the 234 BobHLH proteins; Table S4. The characteristics of BobHLH family in B. oleracea; Table S5. Organ/tissue-expression FPKM values of BobHLH genes extracted from RNA-Seq data; Table S6. Chilling-expression FPKM values of BobHLH genes extracted from RNA-Seq data; Table S7. The characteristics of 234 BobHLH proteins in B. oleracea; Table S8. The details of cis-acting elements in the 2000 bp upstream of the BobHLH genes; Table S9. Subcellular localization prediction; Table S10. Primers of vector construction. Figure S1. The distribution of conserved motifs in each BobHLH protein; Figure S2. Sequence logos of BobHLH proteins domains; Figure S3. The distribution of exon-intron in BobHLH genes.
Author Contributions: Z.D. conceived and supervised the work. X.S. performed the bioinformatics analysis and drafted the manuscript. F.Y. and W.Z. performed the experiment. J.T., W.Z., S.W., and J.L. provided guidance and manuscript reviews. All authors read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.