Genome-Wide Identification of Gramineae Brassinosteroid-Related Genes and Their Roles in Plant Architecture and Salt Stress Adaptation

Brassinosteroid-related genes are involved in regulating plant growth and stress responses. However, systematic analysis is limited to Gramineae species, and their roles in plant architecture and salt stress remain unclear. In this study, we identified brassinosteroid-related genes in wheat, barley, maize, and sorghum and investigated their evolutionary relationships, conserved domains, transmembrane topologies, promoter sequences, syntenic relationships, and gene/protein structures. Gene and genome duplications led to considerable differences in gene numbers. Specific domains were revealed in several genes (i.e., HvSPY, HvSMOS1, and ZmLIC), indicating diverse functions. Protein-protein interactions suggested their synergistic functions. Their expression profiles were investigated in wheat and maize, which indicated involvement in adaptation to stress and regulation of plant architecture. Several candidate genes for plant architecture (ZmBZR1 and TaGSK1/2/3/4-3D) and salinity resistance (TaMADS22/47/55-4B, TaGRAS19-4B, and TaBRD1-2A.1) were identified. This study is the first to comprehensively investigate brassinosteroid-related plant architecture genes in four Gramineae species and should help elucidate the biological roles of brassinosteroid-related genes in crops.


Introduction
Leaf angle (LA), the angle between the leaf blade and stem, and plant height (PH) are important agronomic factors in determining plant architecture [1]. LA shapes plant architecture and influences light interception, which contribute to increasing crop yield and productivity [2]. PH is associated with grain yield and is involved in crop lodging resistance [3]. Soil salinization inhibits crop growth and development, which, in turn, reduces yield production. With the rapid increase in the global population, ever-increasing food demand has become a considerable challenge. Therefore, uncovering the mechanism regulating plant architecture and salinity resistance is critical.
Brassinosteroid (BR) phytohormones play various roles in plant growth and stress adaptation [4]. Abundant functional information on BRs has been obtained for Arabidopsis and rice [5][6][7], suggesting roles in regulating plant architecture and increasing LA and PH. Notably, both LA and PH are influenced by BR biosynthesis, signaling, and downstreamrelated genes.
The above BR-related genes also play pivotal roles in environmental stress responses. For example, GSKs and BES1/BZR1 have emerged as key factors antagonizing drought, cold, and heat responses [19]. ABA plays an important role in regulating responses to abiotic stress, and BR and ABA exhibit an antagonistic relationship mediated by BRI1, BAK1, and SnRK2s [20]. Furthermore, overexpression of AtDWARF4 in Brassica napus enhances drought tolerance [21].
While the above BR-related genes are known regulators of LA, PH, and stress responses in rice, their functions in plant architecture and salt stress tolerance in Triticum aestivum, Zea mays, Hordeum vulgare, and Sorghum bicolor remain unclear. In the current study, we conducted a genome-wide identification of BR-related genes in these species and determined chromosomal location, phylogenetic analysis, sequence alignment, gene/protein structure, transmembrane structure, synteny relationship, Ka/Ks values, promoter sequences, and protein-protein interactions. In addition, we investigated their potential involvement in regulating plant architecture and salt tolerance in wheat and maize using transcriptome data.

Multiple Sequence Alignment, Secondary Structure, and Phylogenetic Analyses of BR-Related Proteins
Conserved domains were identified from multiple protein sequence alignments ( Figure S2). For example, anchor regions and proline-rich domains were found in the N terminal of DWARF4/D11, while A-C and heme-binding domains were found in the C terminal ( Figure S2-1). Membrane anchor regions and proline-rich, dioxygen-binding, steroidbinding, and heme-binding domains were found in the D2/D3 proteins ( Figure S2-2). Furthermore, TPR, CD I, and C II domains were identified in the N to C terminals of SPY proteins ( Figure S2-8). However, several proteins shared specific motifs with their homologous ones. For example, there was an incomplete TPR domain in HvSPY, but no CD II domain ( Figure S2-8); the AKER and AP2-R1 domains were lost in HvSMOS1 ( Figure S2-11); and the CCCH domain was not in ZmLIC ( Figure S2-14).

Multiple Sequence Alignment, Secondary Structure, and Phylogenetic Analyses of BR-Related Proteins
Conserved domains were identified from multiple protein sequence alignments (Figure S2). For example, anchor regions and proline-rich domains were found in the N terminal of DWARF4/D11, while A-C and heme-binding domains were found in the C terminal ( Figure S2-1). Membrane anchor regions and proline-rich, dioxygen-binding, steroid-binding, and heme-binding domains were found in the D2/D3 proteins ( Figure S2-2). Furthermore, TPR, CD I, and C II domains were identified in the N to C terminals of SPY proteins ( Figure S2-8). However, several proteins shared specific motifs with their homologous ones. For example, there was an incomplete TPR domain in HvSPY, but no CD II domain ( Figure S2-8); the AKER and AP2-R1 domains were lost in HvSMOS1 ( Figure S2 -11); and the CCCH domain was not in ZmLIC ( Figure S2-14).
In terms of secondary structures, the proteins were similar to their rice homologs ( Figure S3 and Table S2). For example, there were 50% alpha-helices, 35% random coils, 11% extended strands, and 4.5% beta-turns in the DWARF4 and D11 proteins. According to phylogenetic analysis, most T. aestivum, H. vulgare, Z. mays, and S. bicolor proteins were highly homologous to the corresponding rice proteins ( Figure 2). For example, all XIAOs were clustered in subgroup A; BRI1 proteins were divided into two classes in subgroup B; and SbBAK1 shared a relatively distant relationship with other genes in group C. In terms of secondary structures, the proteins were similar to their rice homologs ( Figure S3 and Table S2). For example, there were 50% alpha-helices, 35% random coils, 11% extended strands, and 4.5% beta-turns in the DWARF4 and D11 proteins. According to phylogenetic analysis, most T. aestivum, H. vulgare, Z. mays, and S. bicolor proteins were highly homologous to the corresponding rice proteins ( Figure 2). For example, all XIAOs were clustered in subgroup A; BRI1 proteins were divided into two classes in subgroup B; and SbBAK1 shared a relatively distant relationship with other genes in group C.

Transmembrane Structure, Synteny, Gene Structure, Promoter, and Protein-Protein Interaction Network Analyses
Transmembrane analysis showed that about half of the BR-related plant architecture proteins shared transmembrane structures (Figure 1c and Figure S4). Tandem duplications were identified in TaBRD1s (Figure 1d and Figure S5). In total, 37 pairs of segmentally duplicated genes were found in the four Gramineae species, including 33 in T. aestivum, one in H. vulgare, and three in Z. mays (Figure 1e). To analyze the evolutionary relationships and expansion patterns of the BR-related genes, we compared the four Gramineae species with rice and identified 14, 10, 15, and 19 pairs of syntenic orthologous genes, respectively ( Figure S6). The Ka, Ks, and Ka/Ks values were calculated for tandem and segmental gene pairs to identify selective pressure (Table S3). The Ka/Ks ratios of most tandem and segmental gene pairs were less than 1, indicating that they experienced negative selection during evolution.

Transmembrane Structure, Synteny, Gene Structure, Promoter, and Protein-Protein Interaction Network Analyses
Transmembrane analysis showed that about half of the BR-related plant architecture proteins shared transmembrane structures (Figures 1c and S4). Tandem duplications were identified in TaBRD1s (Figures 1d and S5). In total, 37 pairs of segmentally duplicated genes were found in the four Gramineae species, including 33 in T. aestivum, one in H. vulgare, and three in Z. mays (Figure 1e). To analyze the evolutionary relationships and expansion patterns of the BR-related genes, we compared the four Gramineae species with rice and identified 14, 10, 15, and 19 pairs of syntenic orthologous genes, respectively (Figure S6). The Ka, Ks, and Ka/Ks values were calculated for tandem and segmental gene pairs to identify selective pressure (Table S3). The Ka/Ks ratios of most tandem and segmental gene pairs were less than 1, indicating that they experienced negative selection during evolution.

Expression Analysis of Maize BR-Related Genes in Regulating PH
Three different maize hybrids, with extremely distinct PH [i.e., low (L), middle (M), and high (H) groups], were used to conduct RNA sequencing (RNA-seq) [22]. These plants showed increased PH from stages 1 to 3 [i.e., jointing stage (1), big flare period (2) and tasseling stage (3)]. Cluster analysis of the maize BR-related genes was performed at the three key stem developmental stages. Maize BR-related genes were grouped into five clusters. The expression of three genes in cluster 1 decreased in the maize hybrids over time, showing the opposite trend to PH (Figure 3a,b). Among these genes, ZmBZR1 showed the highest expression level (Figure 3c). Coronatine (COR), a functional and structural analog of jasmonic acid (JA), can improve maize lodging resistance by inhibiting ninth internode elongation [23]. The meristem (M, basal 0-1 cm between internode) and elongation regions (E, basal 1-2 cm between internode) of the ninth internode were collected for RNA-seq. In the M and E regions, the expression patterns of the BR-related genes over time were identified in the control and COR-treated groups (Figure 4a,b,f,g). In the M region, nine BR-related genes (two in cluster 2, three in cluster 4, and four in cluster 5) showed a close relationship with PH and were induced or repressed by COR (Figure 4a). For example, ZmLIC and ZmLC2 in cluster 2 increased over time in the control group but were initially inhibited and then induced in the COR-treated group. Details on these genes are presented in Figure 4c-e, with highly expressed genes (ZmLC2, ZmXIAO, ZmBZR1, and ZmTUD1) marked with asterisks. In the E region, genes in clusters 2 and 6 sharing close links with PH showed different expression patterns between the control and COR-treated groups. In cluster 2, ZmBRI1-1, ZmGSK1/2/3/4-2, and ZmMADS22/47/55-1 gradually increased from E1C to E4C, whereas their expression levels became stable after COR treatment at the same time points (Figure 4f,h). The transcription patterns of the seven genes in cluster 6 decreased over time in the control group but fluctuated slightly in the COR-treated group (Figure 4f,i). In addition, 14 maize BR-related genes were either down-regulated (cluster A) or up-regulated (cluster B) at most time points in the COR-treated group ( Figure 4j). Based on gene co-expression network analysis, ZmGSK1 was identified as a likely key regulator of COR-mediated internode elongation (Figure 4k).

Transcriptome Analysis of Wheat BR-Related Genes in Regulating LA and PH
In wheat, BR treatment enhanced LA from days 2 to 8, especially at days 5, 7, and 8 ( Figure 5a). To investigate the potential functions of wheat BR-related genes in this process, their expression patterns were analyzed. In total, 12 genes were down-regulated by BR at different time points (Figure 5b). The expression profiles of BR-related genes over time were also identified in both the control and BR-treated groups (Figure 5c,d). These genes were clustered into eight and five clusters over time in the control and BR-treated wheat seedlings, respectively. In the control group, TaD1-1B, TaOFP8-3D, TaOFP8-3B, and TaARF11-2D showed consistent variation trends with LA with time ( Figure 5a,c,e). However, the expression levels of these four genes were very low (Figure 5f). In the BR-treated group, genes in cluster 5, especially TaGSK1/2/3/4-3A and TaGSK1/2/3/4-3D, may be negatively associated with BR-induced lamina inclination (Figure 5a,d,g,h).

Expression Analyses of Wheat BR-Related Genes in Response to Salt Stress
Qing Mai 6 (QM) and Chinese Spring (CS) are salt-tolerant and salt-sensitive wheat varieties, respectively. Here, QM and CS were used for salt stress treatment. Their roots were collected at 6, 12, 24, and 48 h after salt treatment for RNA-seq [25]. After salt treatment, the wheat BR-related genes in CS could be divided into four clusters over time (Figure 7a). The expression patterns of genes in cluster 4 were similar, with reduced chlorophyll content (Figure 7b,g), and seven genes (TaGSK1/2/3/4-3A, TaGSK1/2/3/4-3D, TaARF19-7B, TaARF19-7A, TaARF19-7D, TaIBH1-1A, and TaBRD1-2A.1) showed relatively high expression in this cluster (Figure 7c). In QM, genes in cluster 2 were highly correlated with the salt-related trait (Figure 7d,e,g), and TaGRAS19-4B showed relatively high expression. TaARF11-2D was lowly expressed in both control and salt-treated CS samples. Furthermore, TaDWARF4-4A expression was relatively consistent in CS and QM under control conditions but was induced in QM compared to CS after salt treatment. According to connectivity and logFC, TaDLT-4A and TaMADS22/44/57-4B were identified as core genes in the co-expression network (Figure 8c).

Comparison of BR-Related Genes among Gramineae Species
A total of 35 BR-related genes have been identified in rice. Here, 68 corresponding genes were identified in wheat, showing a one-fold increase compared to rice (Figure 1a). We identified two tandem and 33 segmental BR-related plant architecture gene duplications in wheat (Figure 1d,e). Gene duplications can lead to the evolution of species [28]. Moreover, the wheat genome has undergone genome-wide duplications [29]. Therefore, the expansion and evolution of BR-related genes in wheat may have been induced by gene and genomewide duplications.
Phylogenetic and syntenic analyses of rice-wheat, rice-maize, rice-barley, and ricesorghum were determined to predict the functions of the wheat, maize, barley, and sorghum BR-related genes based on rice orthologs (Figure 2 and Figure S6). The wheat, maize, barley, and sorghum BR-related genes in pairs might originate from common ancestors with rice BR-related genes, indicating their similar roles. In rice, BR-related genes are involved in regulating LA and PH [1,12,[14][15][16]18,[30][31][32]. Therefore, their homologous genes are likely to patriciate in regulating plant architecture.

Putative Functions of Maize and Wheat BR-Related Genes in Plant Architecture and Salt Adaptation
PH is closely associated with maize yield and lodging resistance. Here, three maize hybrids with distinct PH were selected to identify PH-related candidate genes. We found that the expression levels of ZmBZR1, ZmD3, and ZmILI1 decreased over time, in contrast to the increase in PH (Figure 3a,b), indicating that these genes may play negative roles in regulating PH. ZmBZR1 showed the highest expression level (Figure 3c), and therefore may be a key gene. As COR can effectively reduce maize height, we studied the potential roles of the maize BR-related genes in COR-mediated internode shortening and identified their time-series expression patterns (Figure 4). Results showed that several genes, including ZmLC2, ZmXIAO, ZmBZR1, ZmTUD1, ZmBRI1-1, ZmGSK1/2/3/4-2, ZmMADS22/47/55-1, ZmELT1, and ZmBAK1, were likely to participate in this process. BZR1 is a positive regulator of BR signaling, and bzr1 loss-of-function produces a dwarf mutant [33]. In rice, BZR1 directly binds to the promoters of GA biosynthesis genes to affect GA biosynthesis, cell elongation, and PH [34]. The tomato bri1 mutant exhibits short stem and internodes [35]. Rice bri1 is insensitive to BR and severe dwarfism [36]. Zmbri1-RNAi plants show short internodes and dwarf stature [37]. Down-regulating GSK2 to release its repression on BZR1 could increase PH in rice [15]. The gain-of-function elt1-D mutant reduces PH and leaf inclination in rice [18]. In Arabidopsis, bak1 mutation shows a semi-dwarf phenotype [38]. OsBAK1 is also involved in BR signal transduction and plant architecture [39]. Our study indicated that these genes (especially ZmBZR1) are important for maize architecture, and future studies should focus on their exact roles.

Materials and Methods
Step MCScanS' and 'Visualize Gene Structure (Basic)' in TBtools [44]. BR-responsive elements, including BRRE (CGTGT/CG), G-box (CACGTG), E-box (CACTTG), GATGTG, and CTCGC, were identified in the promoter region (-3-kb upstream of the transcription start site). STRING (15 May 2022, v10, http://string-db.org/) (option value > 0.800) was used to construct a protein interaction network according to the interolog proteins from rice. Using Genesis (15 May 2022, v1.8.1, http://genome.tugraz.at/projects.shtml), time-series expression data were analyzed to identify the best-fit model and phase of expression. Gene expression correlations were calculated with the cor. test function in R (v4.1). Using high-quality genes, WGCNA was performed in R (WGCNA package v1.51). To reduce hierarchical clustering, the dynamic tree cut algorithm, with a coefficient of variation cut-off of 0.25, power β of 18, and minimum module size of 30 genes, was used to filter genes with low variation among the samples. Significant module-trait relationships with root length and ROS accumulation were identified by calculating the module eigengene value. Gene co-expression networks were visualized with Cytoscape (v3.8.2, https://cytoscape.org/download.html 13 April 2022).

Plant Materials, BR Treatment, and Transcriptome Analysis
Using water in a 50-mL conical flask, wheat seeds (Bainong 321) were shaken until the seeds with a radicle length of approximately 5 mm. The germinated seeds were immersed in wet gauze above water for 5 d at room temperature. The wheat seedlings were then hydroponically grown with 1/4 Hoagland nutrient solution in an illuminated growth chamber using published culture conditions [45]. After approximately 2 days of growth, the second leaf of the wheat seedlings became visible, and 3 µM 24-epitestosterone (BR) was added to the nutrient solution. Lamina joints for transcriptome (RNA-seq) analysis were sampled on days 2, 3, 5, 7, and 8. Three replicates containing at least 18 plants each were used for each time point.
Total RNA was extracted from the lamina joints using TRIzol reagent (Cwbio, Beijing, China) according to the manufacturer's instructions. The cDNA was synthesized using a reverse transcription kit (Cwbio, Beijing, China). Using Poly-A Purification TruSeq library reagents (Illumina), barcoded cDNA libraries were constructed and then sequenced on the NovaSeq 6000 platform (Illumina). For each library,~10 Gb of high-quality 150-bp pairedend reads were generated. The raw RNA-seq data used for identifying BR-mediated LA can be found in the NCBI Sequence Read Archive database (accession number SRP157960). The fastp software (v0.20.1) was used to evaluate the overall sequencing quality of the reads and to remove low-quality reads in each sample. Hisat2 (v2.1.0) and SAMtools (v1.6) were used to align the high-quality reads to the wheat reference genome sequence IWGSC RefSeq (v2.1) [46]. Stringtie (v1.3.3b) was used to calculate the expression levels of high-confidence genes in each sample. The R package "edgeR" was used to identify DEGs between the control and BR-treated groups with the parameters "|log2 (fold change)| >= 1 and p < 0.05".

Conclusions
In summary, we performed a genome-wide analysis of BR-related genes in four Gramineae species. Their chromosomal location, conserved domain, secondary structure, transmembrane topology, gene structure, promoter sequence, synteny, and protein-protein interaction networks were analyzed. The expression profiles of BR-related genes in different maize and wheat lines with highly distinct PH and responses to COR and BR were delineated (Figures 3-6). As shown in Figures 3-6, ZmBZR1 and TaGSK1/2/3/4-3D appear to be the main regulators of maize and wheat height, respectively. Based on comprehensive analysis (Figures 7-10), TaBRD1-2A.1, TaGRAS19-4B, and TaMADS22/47/55-4B may participate in wheat resistance to salt stress. Future studies should validate the functions of the above-mentioned candidate genes.

Data Availability Statement:
The datasets used in this study were deposited in the NCBI Gene Expression Omnibus (15 May 2022, http://www.ncbi.nlm.nig.gov/geo) with accession number GSE115796, NCBI Sequence Read Archive (SRA) (15 May 2022, http://www.ncbi.nlm.nih.gov/sra) with accession number PRJNA633707, NCBI under BioProject ID PRJNA417210 with SRA submission ID SUB3198378, and NCBI SRA database under SRP062745. The raw RNA-seq data used to identify BR-mediated LA were deposited in the NCBI SRA database under accession number SRP157960. The datasets used and/or analyzed in the current study are available from the corresponding author upon reasonable request.