Genome-wide identification and expression analysis of the WRKY gene family in cabbage (Brassica oleracea var. capitata L.)

Abstract

WRKY proteins are a class of transcription factors involved in plant growth and development, biotic and abiotic stress, and other biological processes. In order to systematically study the structure and function of WRKY proteins in cabbage, we used the whole genome sequence of cabbage in the Brassica database (BRAD) and bioinformatics methods. A total of 150 WRKY genes were identified in cabbage, named BolWRKY1-BolWRKY150. These genes were unevenly distributed on 9 chromosomes. In this study, the WRKY proteins were classified into three classes: I, II and III. Cis-acting elements analysis showed that the BolWRKY gene was involved in multiple hormone and stress responses. The protein interaction network analysis suggested that BolWRKY19/100/105 genes might play similar functions in the same stress response. Based on transcriptome data analysis, we found that different BolWRKY genes showed tissue-specific expression. Quantitative real-time PCR analysis showed that several BolWRKY genes were strongly induced by drought stress and ABA (abscisic acid) treatment, indicating members of BolWRKY gene family were closely related to the stress responses such as drought. Our study lays a foundation for further understanding the regulatory mechanism and function of the WRKY gene in response to abiotic stress.

Keywords:

• Cabbage
• WRKY gene family
• abiotic stress
• expression analysis

Introduction

As one of the largest plant transcription factors (TFs) gene families, the WRKY family plays crucial roles in plant growth and development, metabolism, hormone crosstalk, biotic and abiotic stresses and other biological processes [1]. In the process of plant response to stress, TFs participating in resisting various environment stress by regulating the gene expression, and the upregulation of stress-responsive genes may contribute to plant survival [2]. SPF1 was the first WRKY TF to be identified and isolated from sweet potato [3]. Then, WRKYs (WRKY1/2/3) were identified in Petroselinum crispum, and Rushton et al. (1996) [4] categorized these genes as the WRKY gene family. The DNA binding ability of the WRKY proteins is mainly due to the hallmark heptapeptide at the N-terminus (WRKYGQK) and the zinc finger motif C2HC (C-X7-C-X23-H-X1-C) or C2H2 (CX4-5CX22-23HXHH) at the C-terminus. There are some variant forms of the WRKYGQK heptapeptide, such as WRKYGKK, WRKYGEK, WRRYGQK, WSKYGQK, WKRYGQK, WKKYGQK and WRICGQK [5].

WRKY TFs can be divided into three groups according to the number of WRKY domains (WDs) and the structure of zinc finger motifs. Proteins with two WDs and a C2H2 zinc finger pattern are classified as Group I. Although Group II has the same zinc finger pattern as Group I, it only has one WRKY domain. Group III has one WRKY domain and a C2HC motif. Group II could be further subdivided into five subgroups based on the primary amino acid sequence and their phylogenetic relationship (IIa-IIe) [1, 6, 7]. Cis-acting elements spatially and temporally regulate the gene expression at the transcriptional level. WRKY proteins regulate gene expression by preferentially binding to the cis-acting DNA element W-box (TTGACT/C) of the target genes promoter regions [4]. The ability to bind DNA is reduced significantly when the WRKY domain is mutated [8]. However, non-W-box elements, including PRE4 (TGCGCTT), WT-box (GGACTTTC), and WK-box (TTTTCCAC), can also be recognized and bound by WRKY protein [7].

WRKY TFs are considered to be important regulators of plant growth and development. In short-day conditions, AtWRKY12 and AtWRKY13 participated in late and early flowering of Arabidopsis, respectively [9]. The WRKY gene of Retama raetam was also related to seed dormancy process [10]. Many WRKY genes are involved in the response to drought, salt, heat and cold. PoWRKY1 overexpression in Arabidopsis improved seed germination and root growth in transgenic plants during cold stress and drought [11]. Overexpression of GhWRKY41 and SpWRKY1 resulted in elevated tolerance to salt and drought in Solanum pimpinellifolium [12, 13]. Silencing the GhWRKY21 resulted in elevated tolerance to heat and drought in cotton, reduced damage to cotton leaves, while overexpression of GhWRKY21 resulted in reduced tolerance to drought [14]. Additionally, transgenic expression of BoWRKY10 positively regulated tolerance to drought stresses in kale (Brassica oleracea var. acephala DC) [15]. GmWRKY27 interacts with GmMYB174 to suppress GmNAC29 expression in soybean, thereby improving cold tolerance in transgenic soybean hairy roots [16]. The physiological responses modulated by these WRKYs also include abscisic acid (ABA). In Arabidopsis, AtWRKY18 and AtWRKY60 controlled the ABI5 and ABI4 transcript levels through direct W-box interactions to regulate ABA related gene expression [17]. So far, numerous WRKY gene family members have been identified in many plants species, including 74 members of Arabidopsis [8], 182 members of soybean [16], 136 members of maize [18], 104 members of poplar [19], 97 members of rice [20] and 83 members in tomato [21].

Cabbage (Brassica oleracea var. capitata L.) is a biennial herb of the Brassicaceae family that is widely cultivated around the world for its leaves. The WRKY family has been extensively researched for its critical role in a variety of physiological processes in plants, especially its stress-resistance functions. So far, many studies have explored and analyzed the functions of WRKY genes in Brassicaceae crops. However, the function of the WRKY gene family in cabbage has received little attention. Therefore, the WRKY gene family was identified and analyzed in this study by analyzing the cabbage genome database. We also investigated the expression patterns of several WRKY genes using quantitative real-time polymerase chain reaction (qRT-PCR) under ABA and PEG (polyethylene glycol) treatments. The findings will serve as the foundation for future research into the roles of the BolWRKY genes and the molecular mechanisms of abiotic stress responses in cabbage.

Materials and methods

Plant materials and stress treatments

Seeds of ‘Jingfeng No. 1’ were germinated on filter paper wetted with distilled water at 26 °C for 2 days. They were sown in pots filled with a grass charcoal:vermiculite (2:1) mixture and grown in the greenhouse, maintained at 23 °C and 16 h light/8 h dark photoperiod to the 4-leaf stage. ABA treatment: cabbage seedlings with the same growth were selected and sprayed with 5 mL of ABA (20 mg/L) on the leaf surface of each seedling for 0 h, 3 h, 6 h, 12 h and 24 h. PEG treatment: cabbage seedlings with the same growth were selected. Each seedling was irrigated with 20 mL PEG 8000 (15%) and treated for 0 h, 1 h, 3 h, 6 h, 12 h and 24 h. Three plants were treated for each treatment. Then, all leaves were collected after treatment and frozen in liquid nitrogen immediately and stored at −80 °C for RNA extraction. All trials were performed with three biological replicates.

Genome-wide identification of WRKY genes in cabbage

The cabbage genome was obtained from the BRAD database (http://brassicadb.cn/#/). The Arabidopsis thaliana reference genome TAIR10.1 (GCF_000001735.41) was obtained from the National Centre for Biotechnology Information (NCBI; https://www.ncbi.nlm.nih.gov/). All the non-redundancy WRKY protein sequences of Arabidopsis were used to find the WRKY genes of cabbage in the local BLASTP program with an E-value cut-off <1e⁻¹⁰. Then, the conserved domains of WRKY proteins in cabbage were analyzed and annotated using the HMMER program (http://hmmer.org/) and the SMART (http://smart.embl.de/smart/batch.pl) databases [22]. The online ExPASy tool (https://web.expasy.org/protparam/) [23] was used to calculate the physical and chemical properties of the putative WRKY protein sequences, such as the number of amino acids, total average hydrophobic index and isoelectric point. The subcellular localization of these WRKYs was predicted using Cello v2.5 software (http://cello.life.nctu.edu.tw/) [24].

Phylogenetic analysis, gene structure, conserved motifs identification and cis-acting element analysis

Multiple alignments of WDs were conducted using Clustal X1.83 [25]. A neighbour-joining phylogenetic tree was constructed in MEGA5.0 with 1000 bootstrap replicates and displayed using EvolView (https://evolgenius.info//evolview-v2/). Then, BolWRKY genes were classified into different groups based on the domain characteristics and phylogenetic relationships. The structure of WRKY genes, including intron and exon information, was visualized using the online tools Gene Structure Display Server 2.0 (GSDS 2.0) (http://gsds.gao-lab.org/index.php). MEME (http://meme-suite.org/meme) was employed for the analysis of conserved motifs in BolWRKY proteins with the following parameters: maximum number of motifs, 10. Cis-acting elements in the promoter region of BolWRKY genes (2000 bp upstream of the start codon) were predicted and analyzed using PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [7]. We combined and eliminated the redundant replications among all the cis-elements, and visualized them with TBtools [26].

Chromosomal localization, duplication and synteny analysis

The chromosomal distribution of WRKY genes in cabbage was determined from the annotation file and plotted on cabbage chromosomes using TBtools. The analysis of gene duplication events was performed using TBtools, with the default parameters. We clarified the WRKY genes synteny relationship between cabbage and Arabidopsis using the TBtools software [27].

Interaction analysis

Based on the Arabidopsis thaliana association model, we used the STRING v11.0 (https://string-db.org/) database to construct functional protein association networks and draw proteins map to predict protein interaction pathways [28].

Expression profile analysis of WRKY genes in different tissues

The transcription data were obtained from GEO (NCBI GEO database: GSE42891) and then visualized using TBtools software. The transcription data covering the root, stem, leaf, flower, bud, silique and callus of cabbage were selected to analyse the tissue-specific expression.

RNA extraction and qRT-PCR analysis

The cabbage leaves were quickly transferred to a mortar precooled with liquid nitrogen. The tissue was ground with a pestle while being continuously sprayed with liquid nitrogen until it was ground into powder. Total RNA was isolated by using RNAiso Plus reagent Kit (TaKaRa, Japan). The integrity and concentration of RNA were evaluated by 1.0% agarose gel electrophoresis and Nanodrop 2000c (Thermo Scientific, USA). One microgram of total RNA was used to synthesize cDNA with the ReverTra Ace® qPCR RT Kit (TOYOBO, Japan). Real-time polymerase chain reaction (PCR) was performed using the ChamQ SYBR qPCR Master Mix Without ROX (Vazyme) on the Bio-Rad CFX96^TM Real Time PCR System (Bio-Rad, USA). Primer sequences for all the genes are listed in Supplemental Table S1. The BolActin2 was used as an internal control. All reactions were performed in triplicate, and the 2^−△△CT method was used to calculate the relative expression levels of WRKYs after cabbage leaves were treated with ABA and PEG [29].

Results

Identification and phylogenetic analysis of the WRKY gene family in cabbage

To identify the WRKY family genes, the amino acid sequences of Arabidopsis were used as queries in BLAST analyses, and 150 genes encoding WRKY proteins were identified in the cabbage genome database. Based on the positions of their corresponding genes on chromosomes from top to bottom, the identified WRKY genes are named as WRKY1-WRKY150 (Supplemental Table S2). The protein length of encoded potential proteins by these WRKY genes varied between 130 aa (BolWRKY142) and 1127 aa (BolWRKY6) with PI ranging from 4.59 (BolWRKY56) to 9.97 (BolWRKY3) (Supplemental Table S2).

Multiple alignments were performed to clarify the sequence features of BolWRKYs domain of 150 members in cabbage. We identified the WRKYGQK domain as well as 15 types of heptapeptide variants in the BolWRKY family. The WRKYGKK and WKKYGQK domains were the dominant types (Figure 1), and most mutants of heptapeptide occurred in Group IIc. In BolWRKY113 and BolWRKY114, the WRKYGQK motifs were replaced by WRKSYYK. In BolWRKY58, the WRKYGQK motifs were replaced by WKKYGQR, and a variant (WRKYGRK) was observed in BolWRKY59. Almost all WRKY transcription factors in all plant species can be divided according to the number of WDs and C-terminal zinc finger structures [7]. The neighbour-joining phylogenetic tree was constructed with 150 WRKYs protein sequences in cabbage (Figure 2). On the basis of a previous study of WRKYs [1, 6], the 150 WRKYs used in this study were divided into seven subgroups such as Subgroups I (28), IIa (6), IIb (17), IIc (36), IId (24), IIe (15) and III (24). Compared to Arabidopsis [1], the WRKY members for each Group in cabbage were found to expand nearly two times. The results suggest that the WRKY family has extensively expanded in its common ancestor after its divergence from the Arabidopsis lineage.

Figure 1. Alignment of WRKY conserved domain sequences in cabbage. (MEGA 5.0).

Note: The core sequences are shown in a red box.

Figure 2. Phylogenetic tree of WRKY proteins in cabbage (EvolView).

Gene structure, conserved motif and promoters analysis of BolWRKY genes

To further systematically identify the motifs and gene structure of cabbage, we conducted analysis using the online software GSDS and MEME. All BolWRKYs contained between 2 and 11 exons (Figure 3), and the number of motifs ranged from 2 to 11 (Figure 4A). BolWRKYs genes from the same Group always have similar structures and motifs. Motif 1 was found in the WRKY domain with the heptapeptide, and it existed in all BolWRKYs. Members of IIc, IIe and III had less exons, and the members of Group III all contained two introns. They also had simple structure, most of them only had 3 or 4 motifs. By contrast, members of Group I, IIa and IIb had more conserved motifs and exons. The C-terminal of group I BolWRKYs contained three consecutive motifs, and most IIc BolWRKYs are similar. Furthermore, Motif 4 was only found in Group I, while Motif 8 was found in Groups I and IId. In Group IIc, there were 10 BolWRKYs (BolWRKY25, BolWRKY30, BolWRKY74, BolWRKY92, BolWRKY142, BolWRKY143 and BolWRKY150) had only two exons, as well as Group IIe (BolWRKY48 and BolWRKY87) and IId (BolWRKY78). Most members of Groups I, IIa, and IIb, on the other hand, had complex gene structures with many introns and conserved motifs in the encoded products. The results of gene structure and motif analysis distinguished different classes of BolWRKYs and further corroborated the results of the phylogenetic analysis.

Figure 3. Gene structure and phylogenetic relationship of BolWRKYs. (GSDS 2.0).

Figure 4. Conserved motif and promoter analysis of BolWRKY genes. (A). The sequence information for each motif of BolWRKYs. (B). Conserved motif distribution of BolWRKY gene family. Different colors represent different motifs. (C). Cis-acting elements of the BolWRKY gene family. Different colors and shapes represent different cis-acting elements and their position in the BolWRKYs. (TBtools).

Transcription factor binding sites (TFBSs) are essential components of gene transcription regulation in promoters. To investigate the potential biological roles of BolWRKYs, the online tool PlantCARE was used to predict TFBSs in the 2-kb promoter regions of all BolWRKYs (Figure 4B) [30]. Eighteen cis-acting elements were identified in the promoter regions of all BolWRKY genes, including hormone response, stress response and plant growth and development. In hormone responses, cis-acting elements abscisic acid (ABRE), methyl jasmonate (CGTCA-motif), salicylic acid (TCA-element), gibberellin (TATC-box) and auxin (TGA-element) were identified in the promoter regions of BolWRKYs. For example, CGTCA-motif and ABRE motif were detected in 121 and 115 BolWRKYs, respectively. For stress response, defence and stress responses (TC-rich repeats), drought-inducible element (MBS), low temperature-responsive element and anaerobic induction element (ARE) were also detected. The presence of cis-acting elements indicate that some BolWRKY genes may respond to abiotic stress. In addition, MBS I, WUN-motif and HD-Zip 1 were present in the promoters of some BolWRKYs. WRKY TFs can specifically interact with (T) (T) TGAC (C/T) sequence (W-box) to regulate the expression of regulatory genes or functional genes containing W-box elements in promoters, thus participating in various defence responses and regulating plant growth and development. BolWRKY10 and BolWRKY117 had four W-box cis-elements, indicating that they play an important role in the whole plant growth process.

Chromosomal distribution, duplication and synteny of BolWRKY genes

The replication and distribution of BolWRKY genes in the cabbage genome were analyzed. On 9 cabbage chromosomes, all BolWRKYs were unevenly distributed, there being between 9 (C 03) and 30 (C 05, C 08) genes on each chromosome (Figure 5A). Gene duplication, which includes whole-genome duplication, segmental duplication and tandem duplication, is the most important factor in gene family expansion [7, 31]. About 142 pairs of segmental duplications were identified in the cabbage genome. Among them, C 03 and C 04 had the most collinear genes. BolWRKY genes with tandem duplication events were defined on a single chromosome region within 200 kb containing two or more than two genes [32]. In this study, we found tandem duplication event regions on chromosome 01 (BolWRKY2/3), chromosome 02 (BolWRKY10/20, BolWRKY26/27), chromosome 03 (BolWRKY32/33/34, 57/58/59), chromosome 04 (BolWRKY72/73, BolWRKY81/82), chromosome 06 (BolWRKY94/94, 98/99) and chromosome 07 (BolWRKY111/112, 124/125). These results suggest that the occurrence of both tandem and segmental duplication events might be accounted for the expansion of the WRKY gene family in cabbage. To further explore the phylogenetic mechanisms of the BolWRKY gene family, we constructed collinear relationship of cabbage associated with Arabidopsis (Figure 5B). The number of orthologous gene pairs between cabbage and Arabidopsis was 134.

Figure 5. Duplication and synteny of BolWRKY genes. (A). Synteny map of BolWRKY family genes. The red lines represent the collinear relationship between the BolWRKYs on different chromosomes. (B). Synteny map of WRKYs between Arabidopsis and cabbage. Lines of different colors represent the collinear relationship between BolWRKYs and AtWRKYs on different chromosomes. (TBtools).

Interaction network between BolWRKYs and AtWRKYs

STRING software was used to draw an interaction network map of BolWRKY proteins in cabbage and related proteins in Arabidopsis to identify their functional and physical interactions (Figure 6). The interaction network involved 150 BolWRKY genes and 74 AtWRKY genes. Three BolWRKYs (BolWRKY19/100/105) were homologous to AtWRKY40 and interacted with 6 AtWRKYs (AtWRKY18/33/38/46/53/70). Among them, AtWRKY40/46/70 [31, 32] were involved in drought stress response in Arabidopsis, AtWRKY40 was a positive regulator involved in the drought resistance process and AtWRKY46 and AtWRKY70 were negative regulators. BolWRKY19, 100 and 105 interacted with AtWRKY18 and AtWRKY18 can recognize the W-box in AtWRKY60 promoter and activate its expression, thus participating in plant response to abscisic acid and abiotic stress [33]. BolWRKY38/64/80/113/114, orthologous to AtWRKY33 [34], a modulator of the response to salt stress in Arabidopsis, interacted with AtWRKY46 and AtWRKY53 [35], the key regulatory factor in the response of salt. Frequent interactions were found for BolWRKYs. Even though the functions of most BolWRKYs have not been clarified, this interaction network lays the groundwork for future functional studies of the cabbage WRKY TF gene family [36].

Figure 6. The network of BolWRKY protein functional connections. (STRING v11.0).

Expression patterns of BolWRKY genes in different tissues

Based on available RNA-seq data (NCBI GEO database: GSE42891), the expression patterns of 44 BolWRKY genes in various tissues were visualized using TBtools (Figure 7). The 44 BolWRKY genes showed significantly different expression patterns among the seven cabbage tissues. BolWRKY104 and BolWRKY147 was widely expressed at high levels in all tissues, while BolWRKY72, BolWRKY77, BolWRKY119 and BolWRKY147 exhibited diverse expression patterns. BolWRKY72 was predominantly expressed in leaves, and BolWRKY77 was highly expressed in callus; BolWRKY119 was highly expressed in roots, and BolWRKY147 was predominantly expressed in flowers. The homologs in the class IIc family, BolWRKY1, BolWRKY47, BolWRKY93 and BolWRKY138 showed different expression patterns. The homologs in the class IIb family, BolWRKY18 and BolWRKY104 also showed different expression patterns. BolWRKY18 is specifically expressed in roots, and BolWRKY104 is expressed at low levels in roots. These results indicate that BolWRKY genes play an important role in plant growth and development.

Figure 7. RNA-Seq analysis of BolWRKY genes in different tissues. (TBtools).

Expression patterns of BolWRKY genes under abiotic stresses

To further verify the effects of different abiotic stresses on BolWRKY genes expression, the expression patterns of these genes were assessed by qRT-PCR under PEG and ABA conditions. Cabbage with the same growth trend was selected for stress treatment. As shown in Figure 8, 6 BolWRKYs (BolWRKY19/42/70/72/100/105) responded quickly to ABA treatment and their expression was upregulated. The expression level of BolWRKY19 peaked 3 h after treatment and then declined, but it remained significantly higher than that of the control throughout the treatment. At 3 h after treatment, the expression of BolWRKY42, BolWRKY72, BolWRKY100 and BolWRKY105 peaked, then decreased at 6 h and increased at 12 h. The expression level of BolWRKY70 increased with treatment time, and then decreased to levels lower than in the control Group.

Figure 8. Expression pattern analysis of BolWRKY genes in response to ABA stress. The X-axis represents the RNA samples from the whole plants in different treatment at five time points, from left to right: control (0 h), ABA (3 h, 6 h, 12 h, 24 h). The BolActin2 was used as an endogenous control for normalization of gene transcript level. The Y-axis represents the relative expression levels of BolWRKY genes using the 2^−ΔΔCT method. Data represent the mean ± standard deviation (p < 0.05). (GraphPad Prism 8.0).

As shown in Figure 9, the relative expression of 6 BolWRKYs was determined within 24 h after cabbage was treated with 15% PEG8000. Drought treatment enhanced the transcript levels of three BolWRKYs (BolWRKY19/100/105). The expression level of BolWRKY19 and BolWRKY100 increased at 6 h, then decreased. The expression level of BolWRKY105 increased with treatment time and reached its maximum at 24 h after PEG8000 treatment. Three BolWRKYs (BolWRKY19/100/105) during the entire treatment period were significantly higher than those in control. The gene expression levels of BolWRKY42, BolWRKY70 and BolWRKY72 during the process of treatment were significantly lower than those in the control. The results showed that BolWRKY19, BolWRKY100 and BolWRKY105 rapidly responded to PEG induction, and the expression levels of BolWRKY42, BolWRKY70 and BolWRKY72 were inhibited by PEG treatment.

Figure 9. Expression pattern analysis of BolWRKY genes in response to drought stress. The X-axis represents the RNA samples from the whole plants in different treatment at six time points, from left to right: control (0 h), PEG (1 h, 3 h, 6 h, 12 h, 24 h). The BolActin2 was used as an endogenous control for normalization of gene transcript level. The Y-axis represents the relative expression levels of BolWRKY genes using the 2^−ΔΔCT method. Data represent the mean ± standard deviation (p < 0.05). (GraphPad Prism 8.0).

Discussion

The WRKY gene family, which is one of the largest TF families, is mainly involved in plant growth and development, and various signal transduction process [6]. In this study, genome-wide identification, classification, characterization and phylogenetic analysis of WRKY gene family were conducted in cabbage. The expression of WRKY genes in cabbage was analyzed, and the BolWRKYs involved in ABA and PEG response were identified.

An essential step in clarifying a gene family's function is classification of the genes. According to the newly released BRAD V3.0 Genome database, 150 BolWRKY genes were identified in the cabbage genome, including 148 BolWRKYs previously reported by Yao et al. [37]. A widely accepted family classification system for the WRKY genes has been established in Arabidopsis [38, 39]. Based on the results of Rushton et al. (2010) [6] and Wu et al. (2010) [7], we constructed an integral phylogenetic tree of cabbage WRKY proteins [37].. Among 150 BolWRKY genes, 98 members belonged to Group II, which accounted for the majority of the proteins (65%). This situation is similar to that in Arabidopsis.

AtWRKY TFs with ‘WRKYGQK’ preferentially bind to the core sequence W-box element of the promoter region of the target gene [40]. Previous studies have reported that the ‘WRKYGQK’ sequence varies among plants, resulting in the binding of WRKY transcription factors to cis-regulatory elements other than the W-box [41, 42]. In this study, two variants of Group IIc, BolWRKY29 and BolWRKY47, possessed the ‘WRKYGKK’ motif, and all upstream sequences contained no W-box element, which may be consistent with the case in soybeans where binding involves other cis-regulatory elements due to core sequence changes [42]. It is reported that lacking the ability to recognize and bind W-box was due to core sequence variation [42]. Except for the WRKYGKK sequence, other heptapeptide variants, ‘WKKYGQK’, ‘WKKYGQR’, ‘WRKYGRK’ and ‘WRKSYYK’, appeared in the BolWDs. This may be the result of mutations in the original ancestor of the WRKY transcription factor in cabbage during long-term evolution [43].

Genetic structure analysis revealed that all BolWRKY genes contained between 2 to 11 exons, and BolWRKY genes in the same clades share similar exon/intron structures. However, there was a significant exception in that BolWRKY6 (Group IIe) contained 11 exons, while BolWRKY50, which is similar to BolWRKY6, possessed only three exons. Some of the BolWRKY members in IIc, IId and IIe are missing introns, and it is speculated that the rate of intron loss is higher than the rate of gain during the evolution of cabbage thus leading to the loss of introns [41]. We also determined the motif composition diversification in different Group members of the BolWRKY gene family under conserved motif analysis, and all BolWRKYs had the WRKY motif. Cis-acting elements are required for gene transcription and expression. The findings revealed that the WRKY promoters in cabbage contained a wide range of cis-acting elements, including hormone response elements, stress response elements, W-box elements and plant physiological metabolism related elements. Most WRKY promoters have a high frequency of CGTCA-motif, ABRE, ARE and W-box. Each WRKY promoter contains at least one hormone response element and one stress response element. Hence, these BolWRKYs are thought to be involved in hormonal and stress responses.

Gene duplication is one of the major pathways for genome evolution and gene family expansion, and whole genome duplication, segmental duplication and tandem duplication are the main pathways for plant gene family expansion [27]. In this study, we identified eleven pairs of tandem duplications on cabbage chromosomes 1, 2, 3, 4, 6 and 7, respectively, indicating a minor contribution of tandem duplication to BolWRKY gene family expansion. However, 142 pairs of segmental duplication pairs were identified in cabbage, indicating that fragment duplication events could be considered as a major mechanism driving the expansion of cabbage BolWRKY gene families. This result is similar to that of Arabidopsis [44]. Using comparative synteny mapping, we explored the orthologous and paralogous relationship among cabbage and Arabidopsis. There were 134 collinear pairs between cabbage and Arabidopsis. It can be speculated that gene replication doubled the number of BolWRKY genes, which may be the direct cause of the expansion of the BolWRKY gene family.

Network visualization and statistical analysis of the 150 BolWRKYs and 74 AtWRKYs using the STRING database showed that BolWRKYs (BolWRKY2/14/19/38/42/70/72/100/105) located at the centre of interactions play a critical role in the overall BolWRKY gene regulatory network. AtWRKY40 was a positive regulator involved in the drought resistance process in Arabidopsis. As the orthologues of AtWRKY40 in cabbage, BolWRKY19, BolWRKY100 and BolWRKY105, which belong to Group IIa, may have the same function as AtWRKY40. In addition, AtWRKY18 as the orthologue of BolWRKYs (BolWRKY2/51/125), can recognize the W-box in AtWRKY60 promoter and activate its expression, thereby participating in the response to ABA and abiotic stresses. By further analyzing the expression of BolWRKYs in different tissues of cabbage, some BolWRKY family members were found to be specifically expressed with high and no expression, and some genes were expressed at high levels in roots, flowers, leaves and callus. Among them, BolWRKY19, BolWRKY42, BolWRKY70, BolWRKY72, BolWRKY100 and BolWRKY105 had high expression levels in the roots and leaves. When plants are subjected to adversity stress, the roots are the first to respond to the stress. It is speculated that BolWRKY19, BolWRKY42, BolWRKY70, BolWRKY100 and BolWRKY105, which are mainly expressed in cabbage roots, may play a key role when subjected to abiotic stress, which is consistent with the results of the previous analysis. These ideas were further validated by the quantitative analysis of WRKY gene expression patterns under abiotic stress by quantitative real-time PCR. Our study found that the expression levels of BolWRKY19, BolWRKY42, BolWRKY70, BolWRKY72, BolWRKY100 and BolWRKY105 were upregulated under ABA. At the early stages of stress induction, they responded quickly to stress. BolWRKY100 expression rapidly increased and then sharply decreased, giving a rapid and transient response. Some BolWRKYs have also been instrumental in responding to late stress (e.g. BolWRKY70 in ABA stress; BolWRKY19, BolWRKY100 and BolWRKY105 in PEG stress). BolWRKY70 was continuously expressed during ABA stress, similar to Cheng et al. (2019) [45]. In contrast, BolWRKY42, BolWRKY70 and BolWRKY72 showed different degrees of down-regulation after PEG8000 stress, indicating that these three genes are sensitive to drought stress. These results indicate that the BolWRKY genes have a more complex expression pattern, which indicates that the BolWRKY transcription factor undergoes complex and variable regulation in order to accommodate this complex set of physiological changes during the process of adversity stress. Therefore, to gain insight into the regulation mechanism of the BolWRKY genes, a large amount of research is needed in the later stage.

Yao et al. [37] studied and identified the WRKY family in cabbage, and they found that BolWRKY genes were classified into seven subgroups that belong to three major groups. They found that all members of the WRKY family were unevenly distributed in the pseudo-chromosomes. Our results are consistent with this. They discovered a common variant (WRKYGKK) in 5 BolWRKY proteins and a common variant (WRKYGQR) in BolWRKY40 in the cabbage genome. We found another four heptapeptides, ‘WKKYGQK’, ‘WKKYGQR’, ‘WRKYGRK’ and ‘WRKSYYK’ variants, in the BolWDs. In addition, we predicted the function of WRKY genes by constructing a protein interaction network and analyzing cis-acting elements. The expression patterns of six tissue-specific genes were studied in response to abiotic stresses. These findings provide a theoretical foundation for further research into the mechanism and function of the WRKY gene family in cabbage in response to abiotic stress.

Conclusions

In our study, a total of 150 BolWRKYs were identified and further divided into Groups I, II (Groups IIa, b, c, d, and e) and III. We found similar structure and conserved motifs within the same groups and subgroups. Phylogenetic relationships and collinearity analysis of WRKY genes in cabbage can help to decipher the generation and evolutionary process of BolWRKY genes. Then, we analyzed the BolWRKY expression profiles in different tissues, and six tissue- specific BolWRKY genes were selected for response analysis of PEG and ABA stress. This study will aid in revealing the role of BolWRKY TFs in the regulatory mechanisms of the responses to abiotic stresses and will aid future research into the function of the WRKY genes.

Author contributions

Y.W. and X.Y. conceived and designed the experiments. X.Y., S.Z., W.G., T.W., and Z.F. performed the bioinformatics and data analysis; X.Y. and S.Z performed qRT-PCR measurement and analysis. X.Y wrote the paper.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability statement

The data that support the findings reported in this paper are available from the corresponding author [YW] upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Natural Science Foundation of China (NSFC) (No. 31401908); Natural Science Foundation of Heilongjiang Province (C2016056); Fundamental Research Funds in Heilongjiang provincial universities (135509221, 135409217); University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (UNPYSCT-2017155); Graduate Innovation Project of Qiqihar University, China (YJSCX2021031).