The proteins sequences of the B. napus BZR genes were downloaded from the brassica genome browser GENOSCOPE¹ (Chalhoub et al., 2014) and BnPIR² (Song et al., 2020, 2021), respectively, by employing the six AtBZR genes from A. thaliana genome browser TAIR³ with corresponding Gene ID (At1G19350, At1G75080, At3G50750, At4G36780, At4G18890, and At1G78700). Furthermore, the amino acid sequences from B. oleracea and B. rapa were acquired from the Brassica database BRAD,⁴ using the AtBZRs as a query. The protein sequence with 80% sequence similarity was uploaded to ScanProsite⁵ (De Castro et al., 2006) and Pfam⁶ (Mistry et al., 2021) to search for the presence of the BZR1 domain in the sequence, protein sequences with lacking distinctive matching domains were eliminated. Moreover, multiple sequence alignment analysis was performed by aligning the 46 BZRs protein sequence from A. thaliana, B. oleracea, B. napus, and B. rapa to confirm the conserved domain. Physiochemical characteristics of the BnaBZRs proteins were investigated by the ProtParm⁷ tool (Gasteiger et al., 2005) and the cellular location was predicted by Plant-mPLoc⁸ (Chou and Shen, 2010).

The NJ (Neighbor-joining) method in MEGA X was utilized to plot the evolutionary tree. The tree’s validity was verified by performing 1,000 boot replications, and the Newick format of the phylogenetic tree was further annotated and visualized by the iTOL⁹ web server (Letunic and Bork, 2021). The genes structural features were recognized by uploading the coding and genomic sequences in the gene structure display server GSDS¹⁰ (Hu et al., 2015) and FGENESH (Solovyev et al., 2006). To further verify the intron/exon arrangement, we manually aligned the genomic and coding sequences of each BZR gene. Functional conserved motifs were examined by using the MEME¹¹ program (Bailey and Elkan, 1994) with the minimum motif search configured to 20, and the rest of the settings left to default.

The genome feature file (GFF) of the B. napus genome was employed to inspect the chromosomal distribution of the 21 BnaBZRs, and the map was drawn using the TBtools program (Chen et al., 2020). Moreover, Gene duplication events in BnaBZR genes were performed using the MCScanX toolkit (Wang et al., 2012b). The collinear relationship between BnaBZR genes and BZR genes from B. oleracea, A. thaliana, and B. rapa were acquired by Dual Synteny Plotter function in TBtools. Moreover, positive and purifying selection of the BZR gene family was predicted by utilizing a Selecton Server¹² (Stern et al., 2007). Besides this, the nonsynonymous (Ka) and synonymous mutation (Ks) of a single codon were also calculated by Ka/Ks Calculator 2.0 (Wang et al., 2010).

To determine the Type-I and Type-II functional divergence among clusters of the BnaBZRs, we used the DIVERGE 3.0 program (Gu et al., 2013), and Bayesian posterior probability Q_k was set to 0.9 to predict the specific amino acid site. Protein 3D structure was generated by submitting the BnaA06BZR1 protein sequence into the I-TASSER server¹³ with default parameters (Yang et al., 2015), and specific amino acid sites were marked and visualized using the PyMOL software.¹⁴

Seeds of B. napus cultivar Zhongshuang 11 (ZS11) were obtained from the Oil Crops Research Institute of the Chinese Academy of Agricultural Sciences, Wuhan. The seeds were sprouted on moistened filter paper under the listed conditions: 70% relative humidity and 20 ± 5°C, 16 h light/8 h dark at a light intensity of 50 μmol/m²/s. To investigate the expression profile of BnaBZRs in B. napus, tissues from the adult plant, such as mature-siliqua, root, stem, seed, leaf, flower, flower-bud, and shoot-apex, were collected and stored at −80°C for further experiment. Furthermore, for salt treatment, the three-week-old seedlings were immersed into 400 mM sodium chloride (NaCl; received 100 mM NaCl per day). After reaching a final concentration (400 mM), the leaf samples were collected at 2, 8, and 16 h post-treatment, while for drought treatment, 25% PEG solution was applied to simulate a drought stress-like condition. The leaf tissues were harvested at 2, 8, and 16 h post-treatment. The collected leaves were instantly frozen into the LN2 (liquid nitrogen) and kept at −80°C. To extract RNA and synthesize cDNA, we employed the same procedure as described in our prior study (Sarwar et al., 2021). The synthesized cDNA of the individual sample was utilized as a template in the Thermo Fisher QuantStudio Real-Time PCR system with three separate biological replicates for qRT-PCR analysis. For internal control, the B. napus Actin gene (geneBank ID: XM013858992) was used. The expression patterns of BnaBZRs were measured by the 2^−∆∆Ct method described by Livak and Schmittgen (2001), and a student t-test was performed to find the significant difference. The graph was generated using GraphPad Prism8.0 software (Swift, 1997). The qRT-PCR primers were designed according to the qPCR SYBR-Green Master Mix instruction and listed in Supplementary Table S12.

For protein–protein interaction analysis, the AtBZRs orthologs were searched in the STRING server¹⁵ (Szklarczyk et al., 2019), and the interaction map was drawn using the Cytoscape program. Moreover, to confirm the interaction of miRNA with BnaBZRs, we used the publicly available B. napus miRNAs from miRbase¹⁶ (Kozomara et al., 2019) and PNRD¹⁷ (Yi et al., 2015), and searched against the identified 21 coding sequences of BnaBZRs in psRNAtarget server¹⁸ (Dai and Zhao, 2011). Furthermore, cis-elements were investigated by using the 1.5 bp genomic promoter sequence of the individual BnaBZR gene and uploaded to the plantCARE database¹⁹ (Lescot et al., 2002). Additionally, the functional properties of the BnaBZRs were evaluated using the AtBZRs orthologs in DAVID (Huang et al., 2009a,b) and PANTHER (Mi et al., 2019) program, GO terms with significant value p < 0.05 were considered to be enriched in BZR genes.

The expression data of BnaBZRs under drought, salinity, cold, and S. sclerotiorum stress treatment were acquired from the publically available RNA-seq data sets (CRA001775; Lin et al., 2020) and (SRP311601; Walker et al., 2021). DSEeq2 package in R-studio was applied to measure the genes differential expression, the detected values were adjusted by log2 fold change method, and heatmap graph was created by the TBtools program. To locate the sequence polymorphism of the B. napus BZR gene family, the SNPs from the BnaBZRs of 159 accessions were extracted from the B. napus genomic browser BnPIR² (Song et al., 2020, 2021), and the SNP distribution map was drawn using the TBtools.

The BZR gene family in A. thaliana and other plant species has been studied extensively. However, there are no studies available that have characterized the B. napus BZR genes in depth. Therefore, to examine the B. napus BZR genes, we obtained the peptide sequences of the AtBZR gene family (AtBZR1, AtBZR2, AtBEH1, AtBEH2, AtBEH3, and AtBEH4) from the A. thaliana genome browser and then used it as a query in BLASTP program in the B. napus genome browser GENOSCOPE (see footnote 1; Chalhoub et al., 2014). Moreover, the BnPIR (see footnote 2) B. napus genomic browser was also utilized to further verify the reliability of the BZR gene family, and redundant sequences were eliminated manually. A total of 21 BnaBZRs were recognized, including seven BnaBZRs and 14 BnaBEHs, which were named according to their loci (Supplementary Table S1). Simultaneously, seven and 12 BZRs were recognized in the genome of B. oleracea and B. rapa, respectively. This exhibits that the BnaBZR gene family might derive from their ancestor B. rapa and B. oleracea. Moreover, physicochemical properties of the BnaBZR proteins reveal that the instability index (II) was among 46.64–85.28, which shows that members of BnaBZRs were unstable. In addition, the isoelectric point (pI) values were between 5.14 and 9.81, indicating BnaBZRs were highly basic, and the negative GRAVY (grand average of hydrophobicity) value shows that these proteins are hydrophilic in nature. The amino acids encoded by the BnaBZR proteins are varied from 102aa to 400aa, in which the BnaA01BEH3 amino acid sequence was the shortest contained only 102 amino acids; however, the amino acid length of BnaC08BZR2-1 was the highest, around 400 amino acid residues. The molecular weight ranges from 10,990 KD to 44,156 KD, and the aliphatic index (AI) was between 48.12 and 73.3. Moreover, the subcellular localization signals of BnaBZRs were predicted by Plant-mPLoc. As transcriptional regulatory proteins, all BnaBZRs signals were detected in the nucleus (Supplementary Table S1).

To explain the phylogenetic relationship of the BnaBZRs, an evolutionary tree was generated by aligning the full-length coding sequence of six A. thaliana, 21 B. napus, 12 B. rapa, and seven B. oleracea BZR genes. As shown in Figure 1, all 47 genes were fell into four groups: group I carried BZR1 and BZR2 genes; group II contained BEH1; and group III contained BEH2 genes, while group IV contained BEH3 and BEH4 genes. Interestingly, the homologs of BEH3 and BEH4 were not detected in the B. oleracea related to those of A. thaliana, B. rapa, and B. napus indicating that the BZR gene family in B. oleracea might experience an evolutionary divergence in structural features. Additionally, the expansion of the BEH genes in A. thaliana, B. napus, and B. rapa take place in group II, III, and group IV due to whole-genome duplication. We also found that each BnaBZR gene shows higher homology with A. thaliana BZR genes. For instance, group I carried AtBZR1 and AtBZR2, along with BnaA06BZR1, BnaC05BZR1, and BnaC08BZR2-1, BnaC08BZR2-2, BnaC06BZR2, whereas BnaBEH genes were highly conserved with the A. thaliana BEH genes. Overall, in B. napus, group I, II, III, and IV had seven, two, four, and eight BZR family members, respectively. The genes cluster within the same subfamily displays a similar role in B. napus. Furthermore, to understand the diversity of BnaBZRs, Gene Structure Display (GSDS) analysis was performed using the coding and genomic sequence of the B. napus and A. thaliana. As shown in Figure 2, one to five exons were predicted in the BnaBZRs. Similar to the A. thaliana BZR gene family, most members in the BnaBZRs subfamily had the same number of exons and intron. For instance, most members of group I, group II, and group III have two to three exons and one to two introns, except for BnaA09BEH2, which contains four introns. Conversely, group IV contained one to three exons and only one intron, except for BnaA07BEH4, with no intron. Results from this analysis indicate that group I, II, and III might receive the original genes because, in the previous studies, it has been reported that the gain of introns is far slower than the intron loss during duplication (Nuruzzaman et al., 2010). Therefore, members in group IV might be derived from the original genes during evolutionary processes.

Our study observed that the B. napus BZR gene family was grouped into four groups. To determine the potential functional divergence between groups, we tested the Type I (ϴI) and Type II (ϴII) functional divergence utilizing the CodeML package in DIVERGE v2.0 (Gu and Vander Velden, 2002). For this analysis, we use group I, III, and IV members, while group II was excluded because it comprises only two members of the BnaBZRs. The results indicate that the all-divergence coefficients of Type I among group I/III, group I/IV, and group III/IV were statistically significant, with the ϴI and LRT values varying from 0.99 to 0.646 and 0.86 to 33.9, respectively, indicating that the amino acid sites among different groups of BnaBZR gene family have altered evolutionary constraints, which lead to group-specific functional evolution after divergence (Table 2). In contrast, the coefficients of Type II among group I/IV and group III/IV were also statistically significant with an ϴII value of 0.319 and 0.082, respectively, indicating the radical shift in amino acids properties of these pairs during evolution. However, the coefficient of ϴII among group I/III was not significant. We also found some crucial amino acid residues sites that may affect the functional divergence of the B. napus BZR gene family. To reduce false-positive results, posterior probability (Q_k > 0.9) was selected to detect specific amino acid residues among the BnaBZR family (Table 2). The results predict 21 similar Type-I sites among group I/III and group I/IV. In comparison, only two sites were examined in group III/IV, implying that significant functional divergence might exist in group I/III and group I/IV as compared to group III/IV. Conversely, no Type-II functional divergence sites were predicted in group I/III and group I/IV. However, three sites occurred in group III/IV, of which 99T, 100R, and 155S were examined in both ϴI and ϴII, indicating the altered selective constraints may happen at these sites.

To predict the three-dimensional structure of the BnaBZR proteins, we utilized the representative BnaA06BZR1 protein sequence in the I-TASSER database (Zhang, 2008). The normalized C-score for the predicted model was > −3, suggesting the best template for the BnaBZR1 protein 3-D structure prediction. As shown in Figure 3, we found that the 3D structure of the BnaA06BZR1 contains nine α-helicases and one β-strands, in which the BZR1 domain has only four α-helicases. Furthermore, important amino acid sites, which take part in the functional divergence, were also displayed on the three-dimensional structure of the BnaA06BZR1. Out of 21 amino acid residues, only 14 were dispersed on the BZR1 domain, and the remaining eight sites were spread on the C-terminal region of the BnaA06BZR1, which indicates the BZR1 domain might be more vulnerable to positive selection during the evolution of the BnaBZR gene family. Additionally, amino acids residues 19T and 20R that are present in both ϴI and ϴII functional divergence were located on the N-terminal BZR domain, indicating the essential role of these sites in the evolution of the B. napus BZR gene family (Figure 3).

The full-length 42 coding sequences of the BZR proteins from B. oleracea, B. rapa, A. thaliana, and B. napus were employed to perform the multiple sequence alignment analysis to predict the highly conserved amino acids sites and examine the structural features of the BnaBZR gene family (Supplementary Figure S1). Based on the alignment, we found that six AtBZRs showed the highest percent amino acid homology with BnaBZRs (Supplementary Table S2). Furthermore, a highly conserved domain BRASSINAZOLE RESISTANT1 (BZR1) was identified in the N-terminal region of the BnaBZRs and other denoted plant species (Supplementary Figure S1). In the previous studies, it has been reported that the N-terminal BZR1 domain contains the serine-rich phosphorylation site and participates in modulating the expression of targeted genes to promote various plant physiological processes, such as reproduction, senescence, cell elongation, and cell division (Clouse et al., 1996; Ye et al., 2010). Overall, the domain distribution in B. napus and A. thaliana was similar (Figure 4). To further classify the diversification of the B. napus BZR gene family, we built a graph presenting 20 putative motifs and their distribution on BnaBZRs (Figure 4). Our results showed that all BnaBZRs contain five to 15 highly conserved motifs within a length range from eight to 50 amino acids. In which motif, 1, 3, 4, 5, and 11 encoded a BZR1 domain (Figure 4). The consensus sequences of these motifs were also shown in (Supplementary Figure S2). All BnaBZRs contain motif 5, except for BnaC03BEH3 and BnaC06BEH4, while other remaining motifs existed in some but not in all members (Figure 4). For instance, motif 8 and 10 were only present in group I, II, and III, whereas motif 7 and 18 are only found in group IV, except for BnaC03BEH3 and BnaA01BEH3 lacking motif 7. Furthermore, motif 19 was only detected in BnaC08BZR2-1 and BnaA09BZR2. In contrast, motif 20 was only found in AtBEH4, BnaA07BEH4, and BnaC06BEH4. Motif 12 and 16 were only present in group I, except for BnaA09BZR2 and BnaA07BZR2. However, the structure and motif distribution within the same group were similar, suggesting similar biological functions among the same group (Figure 4).

To predict the distribution of BnaBZRs in B. napus chromosomes, we mapped the 21 BnaBZR members into the 19 B. napus chromosomes (Figure 5A). Based on their physical position, we found that members of the BnaBZR have distributed unevenly between B. napus chromosomes. In which chromosome, A07, A01, CO1, C06, and C08 contain two members of BnaBZR. Chromosomes A03 and C07 hold three members of the BnaBZR, while chromosomes A02, A06, A09, C03, and C05 have only one BnaBZR gene (Figure 5A). Furthermore, the MCScanX tool was also implemented to reveal the duplication of BnaBZRs in B. napus; 29 segmental duplication pairs were identified using this method. However, no tandem duplication pair was detected, implying that segmental duplication events played the central driving role in the functional divergence of BnaBZRs during B. napus genome duplication (Figure 5B). Moreover, genome duplication variation between four Brassicaceae species A. thaliana, B. napus, B. oleracea, and B. rapa was also investigated (Figure 6). The results show that B. napus contains both orthologous and paralogous copies of BZR genes. For instance, six AtBZR genes (AT1G19350, AT1G75080, AT3G50750, AT4G36780, AT4G18890, and AT1G78700) have two, five, two, four, five, and three orthologous copies in the B. napus genome, respectively (Supplementary Table S3). Moreover, the orthologous copies of BZR genes from the B. napus ancestor were also predicted. A total of 12 and 13 BZR genes from B. rapa and B. oleracea showed a syntenic relationship with BnaBZRs, respectively. These orthologous copies might be the reason for the distribution and functional diversification of the B. napus BZR gene family. In addition, we also measured the nonsynonymous (Ka) and synonymous (Ks) ratios for each pair of duplicated BnaBZR proteins (Supplementary Table S4). The results showed that all BnaBZR pairs contain the Ks/Ka ratio < 1, suggesting that BnaBZRs practiced a robust purifying selection (Supplementary Figure S3). However, there are some limitations in calculating the Ka/Ks ratio of the BnaBZRs, because the members of the BnaBZR have minute variations in their sequences, which may bring Ka/Ks value to lower than 1 (Nielsen and Yang, 1998; Nekrutenko et al., 2002). Therefore, we used the site models to predict the selection pressure on single amino acid codons (Yang, 2007; Supplementary Table S5). Model MO is one ratio model, which presumes the same ratio ɯ at all amino acid sites. In our study, we predicted the ratio value ɯ = 0.228 in M0, suggesting the strong purifying selection was the main reason for the evolution of the BnaBZRs. Furthermore, we also compared the model M0 and M3 to predict dN/dS ratio differentiation among codon sites. Interestingly, the log-likelihood 2∆InL was statistically significant (p < 0.01), indicating BnaBZR experienced immense selective pressure across different sites. Additionally, the models M7 (beta) and M8 (beta + ɯ > 1) were also compared, which considered a very stringent test of positive selection (Anisimova et al., 2001), to predict whether positive selection involves promoting divergence of BZR gene family in B. napus. However, no positive selection sites were identified. The results from this analysis indicate that the BnaBZR gene family underwent divergent selective pressure during duplication (Supplementary Table S5).

Protein interaction study provides an essential insight to examine gene function in a systematic study (Pellegrini et al., 2004). Our analysis predicted that the BnaBZR proteins are highly homologous to the A. thaliana BZR proteins. Therefore, to further examine the interacting relationship of the BnaBZRs with its targeting partners, we constructed a network using the AtBZRs orthologs. As shown in Figure 7, we found that the members in the group I interact with more targeted proteins than others. BES1-INTERACTING MYC-LIKE1 (BIM1), BRI1 KINASE INHIBITOR 1 (BKI1), ARABIDOPSIS MYB-LIKE 2 (MYBL2), and PHYTOCHROME INTERACTING FACTOR 4 (PIF4) were predicted to bind with BnaBZR1, BnaBZR2, and BnaBEH1, while PHYTOCHROME B (PHYB) was only found to interact with BnaBZR2 and BnaBEH1 (Supplementary Table S6). However, the biological functions associated with their interaction are still unclear. In contrast, BRASSINOSTEROID-INSENSITIVE 2 (BIN2), BRASSINOSTEROID-INSENSITIVE 1 (BRI1), and BRI1 SUPPRESSOR 1 (BSU1) only interact with BnaBZR1 BnaBZR2, BnaBEH1, BnaBEH2, and BnaBEH4, whereas PHYTOCHROME INTERACTING FACTOR 3 (PIF3), REPRESSOR OF GA (RGA1), RGA-LIKE PROTEIN 3 (RGL3), and AUXIN RESPONSE FACTOR 6 (ARF6) are found to interacts with BnaBZR1 and BnaBZR2. In previous studies, it has been reported that the interaction of BZR1 with ARF6 and PIF mediates auxin-induced plant growth in response to environmental stresses. However, under constant stress conditions, gibberellin negative regulator RGA binds with the BZR1 to inhibit its transcriptional activity (Oh et al., 2014). Additionally, SLEEPY1 (SLY1) and TOPLESS (TPL) are only found to interact with BnaBZR2. Similarly, GA INSENSITIVE DWARF1B (GID1B) shows only interaction with BnaBEH2, which indicates the possible role of BnaBEH2 in regulating BR and GAs crosstalk in response to stresses. It is also worth noting that BnaBEH3 displays no interaction with other proteins (Figure 7). This analysis provides significant clues to understand the relation of the BnaBZRs with unknown proteins in different signaling pathways.

To further examine the individual role of BnaBZR members and how these genes respond to different stresses, we investigated the cis-regulatory elements in the promoter sequences of the BnaBZRs. It has been proposed that the homologs genes with similar roles might consist of the same cis-regulatory elements in their promoter regions. As shown in Supplementary Figure S4, a set of cis-elements belonging to the plant development, stress and defense, and hormonal response-related elements were detected (Supplementary Table S7), but their numbers and distribution were uneven in all BnaBZRs (Supplementary Figure S4A). For instance, BnaC01BEH3, BnaA03BEH2, and BnaC07BEH2 contained the highest number of development-related, hormone-responsive, and light-responsive cis-elements, respectively. In contrast, BnaA06BZR1 contains many stress and defense-related cis-regulatory elements (Supplementary Figure S4B). Nearly, all of the BnaBZRs hold the basic promoter elements including, CAAT-box and TATA-box. However, several core cis-elements were only being seen in a few BnaBZRs. For instance, Sp1 (light-responsive cis-element) and TGA-box (cis-element take part in ABA responsiveness) were only found in BnaA06BZR1 and BnaC05BZR1. A-Box (cis-element involved in the development process), CCGTCC-Box (cis-acting element activates the meristem development), and DRE-core (regulates drought and cold-responsive gene expression) were detected in BnaC05BZR2, BnaC08BZR2-1, and BnaC06BZR2. Furthermore, OCT (cis-element involved in cell proliferation), GCN4-motif (cis-element involved in nitrogen-related response), chs-CMA2a (light-responsive element), and DRE1 (regulates drought and cold-responsive gene expression) were detected in BnaC06BZR2, BnaC07BEH1, and BnaA01BEH2, respectively. GATT-motif, SARE (cis-element involved in salicylic acid response), 3-AF1 binding site (light-responsive cis-element), CARE (involved in ABA response), and LAMP-element were absent in all BnaBZRs except in BnaC08BZR2-2, BnaA07BZR2, BnaA03BEH1, BnaC07BEH2, and BnaC07BEH3, respectively. Similarly, AT1-motif, Box-III, motif-I, and CTAG-motif were only found in BnaC01BHE3. In addition, CAG-motif, ATC-motif, and dOCT (cis-element involved in cell cycle responses) were detected in BnaC01BEH2 and BnaC03BEH3. Overall, light-responsive cis-elements were the most abundant in all BnaBZRs promoters, which indicates the transcriptional activity of the BnaBZRs might be induced by light. Additionally, cis-elements related to hormone responses, such as auxin (AUX), gibberellins (GA), ABA, and salicylic acid (SA), were predicted, in which cis-element involved in ABA responses were the most common in all BnaBZRs. Moreover, stress and defense responsive cis-elements, including ARE (anaerobic regulatory element), WUN-motif (regulates wound-related response), LTR (associated with cold response), and STRE, which activated by the several stress condition, especially heat shock, nutrient deficiency, and osmotic stress were also predicted. Results from this analysis showed that the BnaBZRs contain various kinds of development and stress-related responsive elements, indicating that the BnaBZRs might control plant physiology in response to several phytohormones and environmental stresses.

To further explore the potential functions of the BnaBZR gene family, we study their expression patterns in cotyledon, leaf, silique, root, petal, bud, lower stem, sepal, vegetative rosette, pollen, upper stem, middle stem, seed, and filament of the B. napus variant Zhongshuang 11 (ZS11), by using the transcriptomic data available in BnTIR database http://yanglab.hzau.edu.cn/BnTIR. As shown in Supplementary Figure S5, 17 members of the BnaBZR gene family exhibited higher expression in sepal, petal, seed, lower stem, middle stem, upper stem, root, vegetative rosette, cotyledon, filament, and siliqua, while no or mild expression was observed in pollen and bud (Supplementary Table S8.1). In detail, pairs of homologs genes BnaC08BZR2-1/BnaC08BZR2-2, BnaA09BZR2, and BnaA07BZR2/BnaC06BZR2 showed higher expression in sepal, petal, and seed. In contrast, BnaA01BEH2, BnaA03BEH3, BnaC07BEH3, and BnaA02BEH4 higher expressions were observed in cotyledon and siliqua. Nevertheless, BnaA06BZR1 and BnaC05BZR1 display relatively lower expression in four tissues (root, sepal, upper stem, and vegetative rosette; Supplementary Figure S5). To further validate the organ-specific expression of the BnaBZRs, the relative expression pattern of the six putative BnaBZRs was determined in eight tissues of B. napus cultivar Zhongshuang 11 (ZS11) including mature-silique, shoot-apex, seed, leaf, flower-bud, flower, stem, and root by qRT-PCR. As shown in Figure 8, two of six BnaBZRs (BnaA06BZR1 and BnaC06BZR2) were consistently expressed in all selected tissues, except for BnaA06BZR1, which shows relatively lower expression in flower-bud. Additionally, the increased expression level of the BnaCO6BZR2, BnaC07BEH1, BnaA03BEH3, and BnaA02BEH4 was mainly observed in the mature siliqua (Figure 8, Supplementary Table S8.2). In contrast, BnaC07BEH2 was expressed highly in the leaf. Combined with transcriptomic data, the expression pattern of the BnaC07BEH1, BnaA03BEH3, and BnaA02BEH4 was also observed in the siliqua. However, the correlation test between the two datasets showed a non-significant association (r = 0.435, p = 0.713). The non-significant association between transcriptome data and qRT-PCR might be related to harvesting of siliqua at 38 and 44 days after flowering, respectively. Therefore, the fold change should not be expected to have same for both datasets. Overall, most members of the BnaBZRs expressed differently in the plant tissues, suggesting the key roles of BnaBZRs in various tissues during plant development.

To gain further insights into the BnaBZR gene family in response to the environmental stresses, we investigated their expression levels in drought, salinity, cold, and S. sclerotiorum stresses using publicly available RNA-seq data (CRA001775 and SRP311601). As shown in Supplementary Figure S6A, we found that the expression level of group I and group II was downregulated after 2 h of the drought treatment (Supplementary Table S9.1). Nevertheless, the expression pattern of the BnaA03BEH3 and BnaC06BEH4 were induced. In contrast, the expression pattern of the four BnaBZRs (BnaA07BZR2, BnaC06BZR2, BnaA03BEH3, and BnaC07BEH3) was peaked at 8 h of drought treatment. Under salinity treatment, the expression pattern of the BnaA03BEH3 and BnaC07BEH3 was first induced at 4 h but then reduced in 24 h of salinity treatment. Furthermore, under cold treatment, the RNA transcripts of the group IV members BnaC03BEH3, BnaA02BEH4, and BnaA07BEH4 were upregulated during 24 h of cold treatment. However, BnaA07BEH4 and BnaC06BEH4-induced expression was only observed in S. sclerotiorum infection (Supplementary Figure S6B). To elucidate the functions of the BnaBZRs in abiotic stresses, we monitored their expression changes by qRT-PCR in leaves of B. napus cultivar (ZS11) treated with different stresses, such as drought and salinity (Figure 9). On the whole, the expression patterns of the BnaBZRs were significantly changed during 8–16 h of drought and salt treatment (Supplementary Table S9.2). During salt treatment, BnaA06BZR1, BnaC06BZR2, and BnaA02BEH4 were dramatically induced at 8 h and reduced at 16 h, while BnaC07BEH1, BnaA03BEH3, and BnaC07BEH2 showed downregulation in their expression patterns, except for BnaC07BEH2, which showed upregulated expression at 16 h of drought treatment. Furthermore, expression patterns of the BnaC07BEH2 and BnaA02BEH4 were significantly increased at 8 h but then reduced at 16 h of NaCl treatment, indicating the expression patterns of the BnaC07BEH2 and BnaA02BEH4 were induced in more than one stress. In our, organ-specific expression analysis, we observed the different expression patterns of the BnaBZRs in various tissues (Supplementary Figure S5, Figure 8). Therefore, it is possible that the BnaBZR gene family in B. napus responded to drought, salinity, cold, and S. sclerotiorum stress might play a diverse regulatory function in various tissues.

MicroRNAs regulate gene expression by cleaving the mRNA of the targeted genes in response to stresses (Song et al., 2019). Identifying the miRNA-targeted genes employing the computational approaches is vital for predicting the transcriptional regulatory role of the gene family. Therefore, to understand the transcription network and post-transcriptional regulation of the BnaBZRs, we examine the potential miRNA target sites using the psRNATarget database. Out of the 21 BnaBZRs, only nine genes were targeted by 12 novel bna-miRNA, with an average of one to two targets per bna-miRNA (Table 3). The two members of group I were predicted to be targeted by three novel bna-miRNA (bna-miR395, bna-miR171, and bna-miR170), whereas six members of group IV were targeted by the bna-miR158, bna-miR164, bna-miR159, bna-miRN273, and bna-miRN292. In addition, only one member of group II (BnaC07BEH1) was targeted by the bna-miRN291 and bna-miRN285. However, except for the bna-miRN285/BnaC07BEH1 and bna-miR164/BnaA02BEH4, the majority of the bna-miRNA-targeted genes was expected to be silenced by the cleavage inhibition. Furthermore, we have also predicted the UPE value, which is within 8.94–21.43, needed to unpair the secondary structure around the target site. Nevertheless, most of the bna-miRNA identified in this study, such as bna-miRNA395 (Huang et al., 2010; Zhang et al., 2013a), bna-miR171 (Jian et al., 2016), and bna-miR159, are stimulated by the environmental cues, suggesting that different bna-miRNAs direct post-transcriptional regulation of some of the members of the BnaBZR gene family during biotic and abiotic stresses.

It is well reported that a family of transcriptional factors BZR interacts with several sequences to regulate plant tolerance to environmental cues. To further examine the biological role of BnaBZRs, we conducted the gene ontology (GO) enrichment analysis. Our results showed that a total of 43 GO terms were categorized into three groups, BP (biological process), CC (cellular component), and MF (molecular function; Supplementary Table S10). Among BP category, all BnaBZRs were within the “metabolic processes (GO:1901360, GO:0006725, GO:0034641, GO:0046483, GO:0006807),” “biosynthetic processes (GO:1901576, GO:0034645, GO:0009059, GO:0009058),” “gene expression (GO:0010467),” and “cellular process (GO:0050794),” whereas BnaC05G0163000ZS, BnaA06G0135100ZS, BnaA02G0233100ZS, BnaA01G0099000ZS, and BnaC03G0707200ZS might be involved in plant tolerance to stresses (GO:0006950, GO:0006952, GO:0006955, GO:0009607, and GO:0009617). Under the CC category, all BnaBZRs were in the nucleus (GO:0005634), while some of the genes, such as BnaA09G0609500ZS, BnaC08G0463100ZS, BnaC05G0163000ZS, BnaA06G0135100ZS, BnaC06G0273700ZS, BnaA07G0248700ZS, BnaC08G0253000ZS, BnaA03G0416600ZS, and BnaC07G0389100ZS, were also predicted in the cytosol (GO:0005829) and cytoplasm (GO:0005737) category, respectively. Additionally, in the MF category, all BnaBZRs were significantly enriched in binding (GO:0003677, GO:0003700, and GO:1901363) together with “organic cyclic compound binding (GO:0097159)” (Supplementary Figure S7). Results from this analysis displayed the diverse roles of the BnaBZRs in various biological processes, mainly in cellular and metabolic processes and biological process regulation.

To locate the sequence polymorphism of the BnaBZRs, we used the sequencing data of 159 accessions of rape seed from the B. napus genomic browser [BnPIR (Song et al., 2020, 2021)]. On the whole, 55.5% of SNPs were predicted in the BnaBZR gene family (Supplementary Table S11). However, the SNP percentage frequency on each subfamily was dissimilar. For instance, group I, II, III, and IV hold an average of 49, 18, 39.2, and 64.8 SNPs, respectively. For BnaA07BEH4 and BnaC06BEH, no SNP was found. Meanwhile, we also observed that the SNPs density of the BnaBZRs on CC chromosomes was relatively higher than the AA chromosomes. Furthermore, a detailed SNPs distribution on the coding and non-coding regions of the BnaBZR gene family was also shown (Figure 10), which display that the distribution of SNPs in introns is varied between 9.5–52.8%, in which the highest number of SNPs (224) was found in the intron region of BnaC03BEH3. However, in exons, the distribution of SNPs was 8.5–20.14%, whereas the lowest number of SNPs (4) was found in the exon region of BnaA03BEH1. Results from this analysis indicate that the sequence variation in the members of the BnaBZRs might be associated with their differential expression pattern under different conditions.