Genome-Wide Identification of Brassicaceae Hormone-Related Transcription Factors and Their Roles in Stress Adaptation and Plant Height Regulation in Allotetraploid Rapeseed

Phytohormone-related transcription factors (TFs) are involved in regulating stress responses and plant growth. However, systematic analysis of these TFs in Brassicaceae is limited, and their functions in stress adaptation and plant height (PH) regulation remain unclear. In this study, 2115 hormone-related TFs were identified in nine Brassicaceae species. Specific domains were found in several Brassicaceae hormone-related TFs, which may be associated with diverse functions. Syntenic analysis indicated that expansion of these genes was mainly caused by segmental duplication, with whole-genome duplication occurring in some species. Differential expression analysis and gene co-expression network analysis identified seven phytohormone-related TFs (BnaWRKY7, 21, 32, 38, 52, BnaGL3-4, and BnaAREB2-5) as possible key genes for cadmium (Cd) toxicity, salinity stress, and potassium (K) and nitrogen (N) deficiencies. Furthermore, BnaWRKY42 and BnaARR21 may play essential roles in plant height. Weighted gene co-expression network analysis (WGCNA) identified 15 phytohormone-related TFs and their potential target genes regulating stress adaptation and plant height. Among the above genes, BnaWRKY56 and BnaWRKY60 responded to four different stresses simultaneously, and BnaWRKY42 was identified in two dwarf rapeseeds. In summary, several candidate genes for stress resistance (BnaWRKY56 and BnaWRKY60) and plant height (BnaWRKY42) were identified. These findings should help elucidate the biological roles of Brassicaceae hormone-related TFs, and the identified candidate genes should provide a genetic resource for the potential development of stress-tolerant and dwarf oilseed plants.


Phylogenetic and Syntenic Analyses of Brassicaceae Hormone-Related TFs
Phylogenetic trees showed that most B. napus, B. carinata, B. juncea, B. nigra, B. oleracea, B. rapa, C. sativa, and C. rubella hormone-related TFs were highly homologous to the Arabidopsis TFs ( Figure S4). In addition, the TFs were divided into distinct groups depending on phylogenetic analysis. For example, RGL3s were clustered in subgroup a, RGL2s in sub-group b, RGL1s in sub-group c, GAIs in sub-groups d and e, and RGAs in sub-group f ( Figure S4-4).

Expression Profiles of B. napus Hormone-Related TFs in Response to Abiotic Stress
To identify the potential roles of rapeseed hormone-related TFs in adapting to abiotic stress, their responses to four stresses were investigated. FPKM (fragments per kilobase of transcript sequence per million mapped reads) was calculated to assess gene expression levels, and p < 0.05 and |log2(fold-change)| ≥ 1 were set as the criteria for identifying differentially expressed genes (DEGs). The DEGs were grouped into two clusters (1 and 2) according to their expression level. We then identified hub genes using Cytoscape (v3.8.2). Deleted genes with FPKM values < 1 and correlation value > 0.95 were used as the threshold for screening interactions between genes.
In shoots, 46 and 28 hormone-related TFs were increased and decreased, respectively, after Cd treatment (Figure 1a). Among these DEGs, BnaWRKY7 was a core gene based on gene co-expression network analysis (GCNA) (Figure 1b). In roots, BnaWRKY49, BnaWRKY38, BnaBZR11, and BnaBZR23 in cluster 1 were up-regulated by Cd, while BnaARF33 and BnaARR75 in cluster 2 were significantly reduced after Cd treatment ( Figure 1c). BnaWRKY38 may be a hub gene responding to Cd in roots (Figure 1d). Salinity severely inhibited rapeseed growth and yield [36]. In shoots, BnaWRKY12 in cluster 1 was inhibited by salt, whereas BnaWRKY49 in cluster 2 was induced ( Figure 1e). Among them, BnaAREB2-5 was identified as a hub gene in the gene co-expression network (Figure 1f). In roots, 86 phytohormone-related TFs were regulated by salt, with down-regulated TFs K is an essential nutrient for plant growth and development [37]. Here, K shortage altered the expression of 121 phytohormone-related TFs (Figure 2a,c). In shoots, genes in cluster 1 were induced; in contrast, genes in cluster 2 were inhibited ( Figure 2a). Among them, BnaWRKY21 was identified as a key gene (Figure 2b). Fewer genes were regulated by K stress in the roots (Figure 2c), and BnaGL3-4 was identified as a likely core gene ( Figure 2d). N is another important nutrient required for crop growth [38]. Here, two phytohormone-related TFs were under the control of N in shoots ( Figure 2e). BnaWRKY49 expression was markedly down-regulated under low N stress in roots, while seven other genes were induced ( Figure 2f). Among the DEGs, BnaWRKY32 was identified as a hub gene ( Figure 2g).

Expression Profiles of Hormone-Related TFs in Rapeseed with Different PH and Stem Breaking Resistance (SBR)
Both PH and SBR are crucial agronomic traits [39]. Comparative transcriptome analysis between dwarf (df59, ed1, Ldt, and dwf) and wild-type (WT) rapeseeds was performed to define candidate genes. Compared with df59, nine hormone-related TFs were down-regulated in WT rapeseeds, while 18 were up-regulated ( Figure 3a). As shown in Figure 3b, BnaWRKY42 was identified as a hub gene among the above DEGs. In ed1, BnaMYC2-2, BnaARR21, and BnaARR71 were down-regulated, while 12 other TFs (especially BnaWRKY36 and BnaWRKY59) were up-regulated ( Figure 3c). The key gene BnaARR21 was appraised according to the gene co-expression network ( Figure 3d). In Ldt, five rapeseed hormone-related TFs exhibited low expression ( Figure S7a). BnaWRKY29, BnaAREB2-6, and BnaWRKY38 were clearly increased in dwf ( Figure S7b). Compared with rapeseeds with low SBR during flowering (FL), the expression levels of six phytohormonerelated TFs were significantly altered with high SBR during flowering (FH) ( Figure S7c). Seven genes were down-regulated in rapeseeds with high SBR during silique development (SH) compared to those with low SBR (SL) ( Figure S7d).

Identification of Weighted Gene Co-Expression Network Analysis (WGCNA) Modules and Hub Genes Associated with Target Traits
WGCNA uses data from all genes to identify gene sets of interest, rather than genes only showing differential expression, and to analyze significant associations with phenotypes. WGCNA has two main advantages: i.e., loss of fewer genes and the ability to collate many genes into gene sets and identify their association with phenotypes without multiple hypothesis testing [40]. Therefore, we used WGCNA to analyze the RNA-sequencing data of Cd, salt, K, and N treatments and six dwarf mutants to identify hub genes and their target genes involved in stress adaptation and PH regulation.

Identification of Weighted Gene Co-Expression Network Analysis (WGCNA) Modules and Hub Genes Associated with Target Traits
WGCNA uses data from all genes to identify gene sets of interest, rather than genes only showing differential expression, and to analyze significant associations with phenotypes. WGCNA has two main advantages: i.e., loss of fewer genes and the ability to collate many genes into gene sets and identify their association with phenotypes without multiple hypothesis testing [40]. Therefore, we used WGCNA to analyze the RNA-sequencing    WGCNA was used to identify hub genes in response to K stress. As shown in Figure  6a, the "lightcyan" module (r = 0.92 and p < 0.05) was positively correlated with biomass, while the "turquoise" module (r = −0.92 and p < 0.05) was negatively correlated with SPAD. Gene interaction networks were constructed for these modules (Figure 6b,c and Tables S4-5 and S4-6). Three key genes (BnaEIL1-2, BnaWRKY56, and BnaARR14) were selected for their high connectivity (Figure 6b,c). Four K + transport genes (BnaKUP5, two BnaKUP6s and BnaPCP) [43,44] were found in the "turquoise" module, and W-boxes and RR binding cis-elements were present in their promoters, suggesting they may be the targets of BnaWRKY56 and BnaARR14 (Figure 6d). Under N stress conditions, the "green" module (r = −0.83 and p < 0.05) was negatively correlated with SPAD (Figure 7a). BnaBZR1 WGCNA was used to identify hub genes in response to K stress. As shown in Figure 6a, the "lightcyan" module (r = 0.92 and p < 0.05) was positively correlated with biomass, while the "turquoise" module (r = −0.92 and p < 0.05) was negatively correlated with SPAD. Gene interaction networks were constructed for these modules (Figure 6b,c and Tables S4-5 and S4-6). Three key genes (BnaEIL1-2, BnaWRKY56, and BnaARR14) were selected for their high connectivity (Figure 6b,c). Four K + transport genes (BnaKUP5, two BnaKUP6s and BnaPCP) [43,44] were found in the "turquoise" module, and W-boxes and RR binding ciselements were present in their promoters, suggesting they may be the targets of BnaWRKY56 and BnaARR14 (Figure 6d). Under N stress conditions, the "green" module (r = −0.83 and p < 0.05) was negatively correlated with SPAD (Figure 7a). BnaBZR1 and BnaBZR14 were identified as critical genes in the "green" module ( Figure 7b and Table S4-7). and BnaBZR14 were identified as critical genes in the "green" module ( Figure 7b and Table S4-7).    The relationships between WGCNA modules and PH were also investigated. A total of 38 modules were obtained, with the "darkmagenta" module (r = −0.84 and p < 0.05) showing a negative correlation with PH (Figure 8a). We found several phytohormonerelated TFs in the "darkmagenta" module (Table S4- 8). BnaARF42 and BnaARF26 with high connectivity were identified as hub TFs (Figure 8b). ARFs are reported to regulate PH [45].
Here, two ARF-binding cis-elements were found in BnaARF10, BnaARF18, and BnaARF54, indicating they were the targets of BnaARF42 and BnaARF26 (Figure 8c and Table S5).
The relationships between WGCNA modules and PH were also investigated. A total of 38 modules were obtained, with the "darkmagenta" module (r = −0.84 and p < 0.05) showing a negative correlation with PH (Figure 8a). We found several phytohormonerelated TFs in the "darkmagenta" module (Table S4- 8). BnaARF42 and BnaARF26 with high connectivity were identified as hub TFs (Figure 8b). ARFs are reported to regulate PH [45]. Here, two ARF-binding cis-elements were found in BnaARF10, BnaARF18, and BnaARF54, indicating they were the targets of BnaARF42 and BnaARF26 (Figure 8c and Table S5).

Venn Analyses of Phytohormone-Related TFs Mediating Stress Adaptation and PH Regulation
Based on differential expression analysis, DEG co-expression analysis, and WGCNA, two Venn diagrams were constructed to investigate the diverse roles of rapeseed hormone-related TFs in regulating Cd, salt, K, and N stress (Figure 9a), and PH (Figure 9b). Many TFs responded

Comparison of Plant Hormone-Related TFs among Brassicaceae
Genes within a family usually exhibit obvious variations during their evolutionar history, which contribute to gene family division, expansion, and functional divergenc [46]. Brassicaceae separated from Arabidopsis approximately 43.2 million years ago an underwent genome triplication, with B. napus (2n = 38), B. carinata (2n = 34), and B. junce
Gene functions are always associated with conserved domains [63]. Several distinct domains, including C2, HSF, and RALF, were identified in some Brassicaceae hormone-related TFs ( Figure S3). Notably, C2, which is associated with calcium-binding [64], was found in BolGL3-1; HSF, which participates in the regulation of biotic and abiotic stress [65], was found in BjuARF4; and RALF, which is involved in biotic and abiotic stress responses [66], was found in BniBZR1. These results indicate that most Brassicaceae hormone-related TFs shared conserved functions, while several genes may vary among Brassicaceae species.
A total of 1291 gene pairs were identified between A. thaliana and non-model Brassicaceae species (Figures S1c and S6), and the functions of Brassicaceae hormone-related TFs were predicted based on homologous Arabidopsis TFs. For instance, AtARF6 and At-BZR1 regulate hypocotyl and stem elongation in Arabidopsis [67] and homologous genes (BcaARF20, CruARF4, and BolBZR12) may share similar roles; CsaARR63 and BnaARR49 may control defense response in Brassicaceae based on the homologous gene AtARR6 [68].

Putative Functions of B. napus Hormone-Related TFs in Regulating Stress Adaptation and PH
We found that most BnaWRKY DEGs were up-regulated, while BnaARRs were downregulated in both shoots and roots under Cd and K treatment (Figures 1a,c and 2a,c). We speculated that Cd poisoning and K starvation may induce SA expression and inhibit CTK expression in rapeseed, thereby enhancing plant Cd resistance and K import and translocation. We also found that salt stress and N deficiency led to a significant decrease in the expression of most BnaWRKYs in the shoots and roots (Figures 1e,f and 2e,f). Therefore, we concluded that salinity and limited N supply may restrict transcription of SA, thereby enhancing plant salt resistance and N uptake. Taken together, our results showed that most phytohormone TFs were responsive to diverse abiotic stresses, implying essential roles in the resistance or adaptation of rapeseed plants to stress.
Various genes and gene families involved in PH regulation have been characterized in plants. Different from our study, however, previous research has primarily focused on specific TFs. In A. thaliana, WRKY46, WRKY54, WRKY70, ARF6, and BZR1 regulate cell elongation and PH [75,76]. To further analyze the regulation of phytohormone TFs on PH, we performed DEG co-expression analysis and WGCNA. Results revealed that PH was mediated by BnaWRKY42, BnaARR21, BnaARF26, and BnaARF42 (Figures 3b,d and 7). Therefore, these genes are likely important for oilseed architecture and stress adaptation.
K-and Cd-related genes (BnaWRKY56 and BnaWRKY60) identified through WGCNA (Figures 4 and 6) responded to all four stress conditions. Furthermore, most key genes revealed through DEG co-expression analysis and WGCNA may play core roles in regulating oilseed resistance to three stresses ( Figure 9a). BnaWRKY38, BnaWRKY42, and BnaWRKY59, which may regulate PH, were identified in two dwarf mutants (Figure 9b). In addition, BnaWRKY38 was identified as a key gene in response to Cd stress, and BnaWRKY42 was differentially expressed under all four stress conditions. This suggests that the phytohormone TFs can simultaneously regulate PH and stress adaptation. Therefore, future studies should focus on the potential functions of these key genes.

Identification and Chromosome Locations of Hormone-Related TFs in Brassicaceae
The protein sequences of Arabidopsis hormone-related TFs were used as queries to BLAST the B. napus, B. oleracea, B. rapa, B. nigra, B. juncea, B. carinata, C. rubella, and C. sativa genomes. We retrieved the phytohormone TF gene sequences using the following databases: Arabidopsis  and visualized using the "Gene Structure View (Advanced)" tool in TBtools (v1.09876) to confirm their highly conserved segments. The "Gene Location Visualize from GTF/GFF" tool in TBtools (v1.09876) was used to visualize chromosomal locations, named according to their chromosome order [77].

Phylogeny, Gene Structure, and Syntenic Analyses of Brassicaceae Hormone-Related TFs
After aligning the full-length protein sequences using ClustalW with default parameters, MEGA X (v10.1.6, University of Pennsylvania, Philadelphia, PA, USA) was used to construct the phylogenetic tree with the maximum-likelihood method. Using generic feature format version 3 files of TFs, gene structure analyses were completed using 'Visualize Gene Structure (Advanced)' in TBtools (v1.09876). We performed syntenic analysis, with results visualized using the "One Step MCScanX" and "Advanced Circos" functions in TBtools (v1.09876) [77]. Genes were determined as duplicates in each genome based on the following [78]: (1) Aligned gene sequences were more than 70% identical, and the length of matching sequences was at least 70% of the longer gene; (2) Duplicates on different chro-mosomes were characterized as segmental duplications. The "Simple Ka/Ks Calculator (NG)" function in TBtools (v1.09876) was used to calculate Ka, Ks, and Ka/Ks values.

Conclusions
In this study, 2115 phytohormone-related TFs were systematically identified in nine Brassicaceae species. Their chromosome locations, gene/protein structures, and phylogenetic and syntenic relationships were characterized. Genes responding to Cd, salt, K, and N adaptation in B. napus were investigated through differential expression analysis and DEG co-expression network analysis. In addition, WGCNA indicated that 15 and two phytohormone-related TFs and their potential target genes responded to stress and PH regulation, respectively. Taken together, BnaWRKY56 and BnaWRKY60 were identified as potential hub genes of rapeseed resistance to stress. BnaWRKY42 may play an essential role in regulating PH. Our results showed that SA-related BnaWRKY TFs may be crucial for regulating PH and stress adaptation in rapeseed. The above-mentioned candidate genes should be validated in future studies.

Data Availability Statement:
The datasets used in this study can be found in published papers. The datasets used and/or analyzed in the current study are available from the corresponding author upon reasonable request.

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