Genome-Wide Identification of Brassinosteroid Signaling Downstream Genes in Nine Rosaceae Species and Analyses of Their Roles in Stem Growth and Stress Response in Apple

Brassinosteroid signaling downstream genes regulate many important agronomic traits in rice. However, information on such genes is limited in Arabidopsis and Rosaceae species. We identified these genes in Arabidopsis and nine Rosaceae species. They were, respectively, named based on chromosomal locations. Segmental duplication and whole-genome duplication under purifying selection, as determined by Ka/Ks analysis, likely contributed to Rosaceae gene expansion. Apple (Malus domestica), Arabidopsis, and rice genes were generally similar, while several Rosaceae genes differed from their rice homologs in various characteristics, such as gene length, subcellular localization, transmembrane topology, conserved domains, secondary structures, and responses to external signals. The brassinosteroid downstream genes in apple were, respectively, induced or repressed by five phytohormones. Furthermore, these apple downstream genes were differentially expressed in different apple grafting combinations (“Nagafu No. 2”/“Malling 9” and “Nagafu No. 2”/“Nagafu No. 2”) and long–short shoot varieties (“Yanfu No. 6” and “Nagafu No. 2”). Responses of the MdBZR genes to diverse stress signals were examined and candidate hub genes were identified. These findings indicated that several brassinosteroid signaling downstream genes in Rosaceae functionally differed from their rice homologs, and certain apple genes may play roles in plant height and stress responses. This study provided valuable information and presented enriched biological theories on brassinosteroid signaling downstream genes in apple. Identification of such genes serve to help expand apple breeding and growth. This study provides useful information for brassinosteroid signaling downstream genes.

Key BR signaling downstream genes and their activities have been widely reported in rice (Supplementary Table S1). The activity of BZR1-2 is regulated by BR signaling, and BZR1-2 can affect the transcriptional levels of genes involved in BR biosynthesis (Yin et al., 2005;Salazar-Henao et al., 2016;Wang et al., 2019). BZR1-2 is involved in stem elongation, cell elongation through cell wall modification, transportation of ions and water, cytoskeletal rearrangement, and gibberellin (GA) biosynthesis and signaling (Xie et al., 2011;Tong et al., 2014). BZR1-2 also participates in biological and abiotic stress responses. In rice, abscisic acid (ABA)-induced drought resistance, senility delay, and salt resistance are mediated by BZR1, ABA INSENSITIVE 5 (ABI5), and ABI3, respectively (Lopez-Molina et al., 2002;Yang et al., 2016). Furthermore, jasmonic acid methyl ester (MeJA) enhances plant resistance to insects through the BZR1 gene (Miyaji et al., 2014). Pathogen-associated molecular pattern (PAMP) receptors promote phosphorylation of BZR2, which controls plant immune responses to pathogens (Kang et al., 2015). DLT encodes a plant-specific GRAS (Tong et al., 2009) and participates in BR negative feedback regulation and GA metabolism in rice . As a CCCH-type zinc finger gene, LIC is induced by BR and regulates tiller angles, plant height, and production in rice (Wang et al., 2008). Acting downstream of BZR1, ILI1 controls shoot elongation and is upregulated by BR (Zhang et al., 2009). The OSH1 protein with a KNOX domain plays an indispensable role in cytokinin (CK)mediated shoot apical meristem (SAM) maintenance (Sentoku et al., 1999). Moreover, several BR catabolic genes are under the control of OSH1 during cell differentiation in SAM (Tsuda et al., 2014). BR treatment down-regulates AtRAV1, influencing lateral root and leaf growth (Hu et al., 2004). The RAV6 protein participates in BR homeostasis and affects plant architecture . SMOS1 (also named RLA1) can form a protein complex with BZR1 in regulating BR signaling and rice architecture (Qiao et al., 2017). OsCSA is a direct target of OsBZR1 and influences rice reproduction and yield (Zhu et al., 2015). OsSPY, an O-linked N-acetylglucosamine transferase, controls BR synthesis and GA signaling by interacting with the DELLA protein (Shimada et al., 2006;Wang et al., 2009). OsGSR1 is involved in crosstalk between GA and BR (Shimada et al., 2006;Wang et al., 2009).
In summary, BR signaling downstream genes in rice share various functions in regulating important agronomic traits, such as plant height, tillering, yield, and environmental adaptations. Rosaceae plants, as important economic crops, exhibit many similar agronomic traits to rice. To date, however, limited information on the BR signaling downstream mechanisms in Rosaceae is available. Previously published Rosaceae genomes (e.g., Malus domestica, Fragaria vesca, Pyrus communis, Prunus persica, Prunus avium, Prunus dulcis, Prunus mume, Rubus occidentalis, and Rosa chinensis) provide useful tools for genome-wide searching of BR downstream genes. Here, we systematically identified BR signaling downstream genes in the above Rosaceae genomes. In total, 196 BR signaling downstream genes were characterized in this study. We analyzed their transmembrane helices, chromosomal locations, conserved domains, evolutionary relationships, gene/protein structures, synteny, and Ka/Ks ratios. Apple, as an important Rosaceae species, occupies a decisive position in the global fruit market in terms of cultivation area and yield. As a perennial woody fruit tree, the apple tree exhibits complex growth performance mediated by multiple hormones. For example, ABA, MeJA, and salicylic acid (SA) play essential roles in stress adaptation in apples . Therefore, it is vital to identify apple BR signaling downstream genes and their functions. Promoter and protein-protein interaction analyses of the apple BR signaling downstream genes were performed, and their expression profiles were analyzed in various tissues. In addition, their responses to BR, brassinazole (BRZ), auxin (IAA), GA, and CK were explored using quantitative reverse transcription polymerase chain reaction (qRT-PCR). Their expression patterns in different shoot varieties ["Yanfu No. 6" (YF) and "Nagafu No. 2" (CF)] and grafting combinations [CF/"Malling 9" (M9) and CF/CF] were confirmed. In addition, the expression profiles of MdBZR genes in response to stress-related signals (ABA, MeJA, SA, drought, and salt) were identified. To provide a foundation for future studies on MdBZR genes, we cloned several candidate genes, i.e., MdBZR1, MdBZR4, MdBZR6, and MdBZR7, and identified the subcellular localization of MdBZR1. The data presented herein not only serve as a valuable resource for elucidating the functions of BR signaling downstream genes in nine Rosaceae species but also provide insights into their potential functions in apple.

Identification, Chromosomal Location, Gene Duplication, Collinearity, and Ka/Ks Ratio Analyses of BR Signaling Downstream Genes in Nine Rosaceae Species
The amino acid sequences of rice BR signaling downstream genes (OsBZR1, OsBZR2, OsDLT, OsLIC, OsILI1, OsOSH1, OsRAVL1/OsRAV6, OsSMOS1, OsCSA, OsSPY, and OsGSR1) were used to blast the Arabidopsis genome 1 and nine Rosaceae databases using Blast v2.7.1 + software. Their distinctive domains were confirmed using pfam 2 , SMART 3 , and sequence alignment analysis. Protein sequences lacking corresponding characteristic domains were discarded from further analyses. Chromosomal locations were completed using the "gene location visualize" tool in TBtools (Chen et al., 2020). Genes were determined as duplicates in each of the genomes based on the following: (1) The aligned gene sequences were more than 70% identical, and the length of matching sequences was at least 70% of the longer gene (Yang et al., 2008). (2) On the same chromosome, duplicated genes separated by less than five genes within a 100 kb region were considered as tandem duplicates (Wang et al., 2010).
(3) Duplicates on different chromosomes were characterized as segmental duplications (Wei et al., 2007). Collinearity analysis and visualization were completed using the "one step MCScanX" and "Amazing Super Circos" tools in TBtools (Chen et al., 2020). The Ka/Ks ratio was calculated using the "simple Ka/Ks calculator" tool in TBtools (Chen et al., 2020).
Analyses of Phylogenetic Trees, Protein-Protein Interactions, Promoters, Gene/Protein Structures, Transmembrane Helices, Chemical Characters, and Sequence Alignments The maximum-likelihood approach in MEGA X was used to construct a phylogenetic tree, and protein sequences were aligned using the ClustalW program with default parameters. STRING 4 (option value >0.800) was used to construct an apple proteinprotein interaction network using rice homologous proteins. In the PlantCARE database, cis-elements were identified in promoter sequences (2 kb upstream of the transcriptional start site). Gene structures were identified using coding sequences and full-length gene sequences with the "gene structure view" tool in TBtools (Chen et al., 2020). PHYRE (v2.0) 5 and SOMPA 6 were used to analyze protein secondary structures, and transmembrane helices were identified in TMHMM (v2.0) 7 .
The plant materials and treatments are summarized in Supplementary Table S2. Gene Expression Analysis, Molecular Cloning, and Subcellular Localization of Apple MdBZR Genes RNA was isolated from the above apple samples (Gambino et al., 2008). A PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) was used to synthesize firststrand cDNA. Primers of apple genes were designed using Primer 3 software (Supplementary Table S3). qRT-PCR analysis was completed using a SYBR Green qPCR Kit (TaKaRa, Dalian, China) on a Bio-Rad CFX 134 Connect Real-Time PCR Detection System (United States). Relative gene expression was calculated with the 2 − Ct method (Livak and Schmittgen, 2001). Target genes were amplified using Phanta HS Super-Fidelity DNA Polymerase (Vazyme, Nanjing, China). After their introduction into the pMD-18T vector (Takara), target genes were sequenced. The coding sequence of the MdBZR1 gene without a termination codon was cloned in pCAMBIA1302 vector 9 and then transferred into tobacco (Nicotiana benthamiana) leaves by the Agrobacterium tumefaciens-mediated method (Yoo et al., 2007). The transgenic tobacco plants were observed using confocal microscopy (TCS SP8, Leica, Germany).

Statistical Analysis
The qRT-PCR was completed using three biological replicates for each sample, and three technical replicates for each reaction. EF-1α (GenBank accession no. DQ341381) served as the internal standard. Standard values, standard errors, and significance analyses of experimental data were subjected to an analysis of variance (ANOVA) at the 0.05 and 0.01 level with OriginPro 9.0 software (OriginLab Inc., Northampton, Massachusetts, United States).

Characterization, Synteny, and Gene Structure Analyses of Genes Encoding BR Signaling Downstream Components in Rosaceae Species
A genome-wide analysis of key BR signaling downstream genes (BZR1-2, DLT, LIC, ILI1, OSH1, RAVL1, RAV6, SMOS1, CSA, SPY, and GSR1) was completed for the nine Rosaceae species. We named them according to their chromosomal or scaffold positions (Supplementary Figure S1). A total of 196 BR signaling downstream genes were identified: 15-23 from nine Rosaceae species (Figure 1A), including 3-8 BZR1/2s and 1-3 other genes. All proteins, except for several SPYs, AtGSR1, and FvGSR1-3, shared instability indices above 40 ( Figure 1B and Supplementary Table S4) and were unstable according to the published criteria (Feng et al., 2019). Most proteins did not contain transmembrane structures, except for MdBZR8, PcBZR3, RcBZR4, and so on ( Figure 1C and Supplementary Figure S2).

Phylogenetic, Conserved Domain, and Secondary Structure Analyses of BR Signaling Downstream Proteins
To identify the characteristic domains of the BR signaling downstream proteins, multiple protein sequences were aligned (Supplementary Figure S6). In the N-terminus of most BZR1/2s, a nuclear localization signal (NLS) and atypical basic helixloop-helix (bHLH) DNA-binding motif with basic, helix1, loop, and helix2 regions were detected (Supplementary Figure S6-1). However, PaBZR1, PcBZR1, and PcBZR2 shared no NLS or bHLH DNA-binding motif, and PaBZR1 lost the S-rich phosphorylation site. There were less consistent PEST and C-terminal domains among all BZRs, except for MdBZR8 and PcBZR4 (Supplementary Figure S6-1). The DLTs contained conserved leucine heptad I domain with IA and IB, VHIID motif, leucine heptad II with A and II B, PFYRE, and SAW domains (Supplementary Figure S6-2). In rice, Arabidopsis, and Rosaceae RAVL1-RAV6s, a conserved B3 DNA-binding domain consisting of ∼90 aa residues was uncovered in the central part (Supplementary Figure S6-3). At the N-terminus, LIC proteins contained a single C-x8-C-x5-C-x3-H (CCCH)type domain. A conserved EELR domain was found in the middle part of the LICs and a less conserved serine-rich region was located at the C-terminus. However, PaLIC lacked the EELR and serine rich domains (Supplementary Figure S6-4). All OSH proteins contained a KNOX domain at their N-terminus, and there was a conserved ELK domain, including one conserved basic end and three conserved helix motifs at the C-terminus (Supplementary Figure S6-5). As AP2 TFs, the SMOS1 proteins contained five conserved domains, including AKER, AP2 (containing highly constant AP2-R1, linker region, and AP2-R2 motifs), EPY, ILS, and C-terminus WTNF domains; however, the EPY, ILS, and WTNF domains were not identified in PaSMOS1-2 (Supplementary Figure S6-6). Highly conserved R2 and R3 domains were detected in the CSA proteins (Supplementary Figure S6-7). There was a C-terminal region containing C, R, G, K, and Y residues in the GSR1 proteins (Supplementary Figure S6-8). Long conserved TRP, CD I, and CD II domains were identified in SPY proteins (Supplementary Figure S6-9). Two helix regions and one small loop region were identified in the ILI1 proteins (Supplementary Figure S6-10).
These protein secondary structures are shown in Supplementary Figure S7 and Supplementary Table S7. Most proteins contained an α helix, extended strand, random coil, and β-turn, and a mass of α helices or random coils were found in all proteins. Random coils and β-turns were, respectively, the most and least abundant structures in the BZR proteins, except for PmBZR3. In space, α helices occupied the main positions in most BZR proteins, except for AtBZR1, OsBZR2, PcBZR4, PaBZR1, RcBZR1-2, RcBZR4, and PmBZ1 (Supplementary  Table S7). The LIC proteins contained more than 70% random coils and less than 5% β-turns, with the proteins predominately made up of random coils in space (Supplementary Figure S7-3 and Supplementary Table S7). The ILI1, OSH1, and GSR1 proteins, which shared a few extended strands and β-turns, mainly consisted of α helices or random coils in space ( Supplementary  Figures S7-4 We next constructed a phylogenetic tree. As expected, Rosaceae BR signaling downstream genes were highly homologous to their corresponding rice and A. thaliana homologs (Figure 2). In group A, BZRs were divided into two subgroups. OsBZR1, AtBZR1-2, and MdBZR3 were clustered in subgroup a2, with other genes located in subgroup a1. The P. avium SMOS1 genes showed relatively distant phylogenetic relationships with other SMOS1s in group B. Rosaceae LIC genes showed a close evolutionary relationship in group C, and OsLIC and AtLIC were closely clustered together. Most RAV genes were located in subclass d1, although some were also identified in subclass d2. In group E, OsOSH1 separated with other OSH1 genes. AtDLT_OsDLT and Rosaceae DLT genes were closely clustered in group F. Class G was divided into subgroup g1, composed of FvILI1-2, RoILI1-2, RoILI1-3, AtILI1, PcILI1-1, PpILI1-2, and PmILI1-3, and subgroup g2, composed of the other ILI1s. All CSA genes clustered together in group H, while PaCSA1 shared a distant relationship with other CSA1 genes. All SPYs were clustered in subclass I, whereas OsSPY was relatively independent. MdGSR1 and RoGSR1-2 closely clustered together in group J.

Promoters and Protein-Protein Interaction Analyses of Apple BR Signaling Downstream Genes
To further investigate apple BR signaling downstream genes, we examined the cis-acting elements in their promoters. We identified ABA-, MeJA-, GA-, IAA-, SA-, and BR-related motifs in these genes (Supplementary Figure S8). For example, ABAresponse motifs (ABRE) and GA-related units (GARE, P-box, or TATC-box) were identified in all genes. In addition, MeJArelated elements (CGTCA-and TGACG-motifs), TGA element or AuxRR-core responding to auxin, SA-related SARE or TCA elements, BR-related elements (BRRE, G-box, E-box, or GATGTG), and stress-related elements (TC-rich repeats, MBS, and DRE) were found in most genes. In addition, elements (MSAlike or MBSI) referring to growth and development were found in MdBZR1, MdBZR4, MdCSA, MdDLT, and MdOSH1-1.
Expression Analysis of BR Signaling Downstream Genes in Different Tissues and Response to BR, BRZ, IAA, GA, and CK Treatments The expression levels of BR signaling downstream genes in STs, YX, YP, MX, MP, YLs, MLs, and Rs were investigated (Figure 3). Results showed that the genes exhibited different expression patterns in the above tissues. In STs, MdBZR4, MdBZR7, and MdBZR8 showed the highest levels. Twelve genes, such as MdBZR1, MdBZR5-6, MdILI1-1, and so on, were highly expressed in the stem tissues (YX, YP, MX, or MP). Nine genes (MdBZR2, MdBZR3, MdBZR4, and so on) were highly expressed in roots. MdBZR4 and MdBZR7-8 were found at the highest levels in STs. Other genes were mainly expressed in the YLs and MLs.
BR, IAA, GA, and CK are known to promote apple tree stem growth and development (Zheng et al., 2018). To identify the roles of BR signaling downstream genes in apple tree stem growth and development, their expression profiles were further investigated (Figures 4, 5). MdBZR1, MdBZR2, MdBZR3, and MdBZR6 were induced by BR, generally. MdBZR4, MdBZR5, and MdBZR8 were separately down-regulated in BR-treated STs. However, the expression pattern of MdBZR7 was irregular in response to BR treatment (Supplementary Figure S12). MdBZR1, MdBZR7 and MdBZR8 were upregulated by IAA, generally (Figure 4). The other MdBZR genes showed irregular expression patterns after IAA treatment (Supplementary Figure S12). GA up-regulated several MdBZR genes, including MdBZR1-4 and MdBZR7-8. For example, MdBZR1 showed a 2-6-fold increase in the GA-treated STs. From days 28 to 56, MdBZR2 expression was higher (p < 0.05) in GAtreated samples than in the controls (Figure 4). GA exposure irregularly influenced MdBZR5 and MdBZR6 (Supplementary Figure S12). The expression patterns of other BR signaling downstream genes were also explored after hormone treatments (Figure 5). At most time points, the expression levels of MdDLT, MdLIC, MdILI1-1, and MdILI1-2 decreased in the BR-treated plants (Figure 5). The expression levels of MdRAVL1s/MdRAV6s, MdSMOS1s, MdGSR1, and MdCSA decreased at 3-5 time points following BR exposure (Figure 5). MdSPY2 showed higher FIGURE 3 | Expression profiles of apple-specific BR downstream genes in different tissues. Each value represents mean ± standard error of three biological replicates. Different lowercase letters indicate significant differences at 0.05 level.

Expression Analysis of BR Signaling Downstream Genes in Two Branch Types and Grafting Combinations of Apple Trees
As several BR signaling downstream genes (BZR1-2, DLT, LIC, ILI1, and SMOS1) are involved in stem elongation Qiao et al., 2017), we analyzed their expression patterns in two branch-type apple trees (YF and CF) (Figure 6). YF exhibited a lower (p < 0.05) shoot elongation rate, internode number, and average internode length than CF . MdBZR1-2, MdBZR4, MdBZR6, MdDLT, MdILI1-1, and MdSMOS1s showed higher expression levels in CF than in YF at most time points. The transcriptional levels of MdLIC were lower (p < 0.05) in CF than in YF at 65, 85, 105, and 125 DABB. Regarding the grafting combinations, CF/CF grows more vigorously than CF/M9 . Several BR signaling downstream genes showed different expression patterns in CF/M9 and CF/CF (Figure 6). For example, MdBZR1 and MdBZR4 showed lower (p < 0.05) levels in CF/CF than in CF/M9 at 3-4 time points. MdBZR5-8, MdDLT, MdILI1s, and MdSMOS1-1 were highly expressed in CF/CF trees at most time points. The expression patterns of the remaining MdBZRs, MdDLT, MdLIC, MdILI1, and MdSMOS1 genes were irregular among YF-CF and CF/M9-CF/CF (Supplementary Figure S12).

Effects of Stress-Related Treatment on MdBZR Gene Expression Profiles
BZR1 and BZR2 play important roles in response to stress (Kang et al., 2015;Yang et al., 2016). We investigated the expression patterns of apple BZR genes after ABA, MeJA, SA, drought, and salt treatments in YLs (Figure 7). MdBZR1-3, MdBZR5-6, and MdBZR8 were generally induced by ABA and MeJA but was down-regulated in SA-treated samples. In ABA-treated samples, MdBZR4 was up-regulated ∼0. 9-, 0. 45-, 11-, and 2fold at 0, 6, 12, and 24 h, respectively. In MeJA-treated leaves, the expression of MdBZR4 was inhibited at 3, 6, and 24 h, but increased at 0 h. MdBZR4 was not affected by SA exposure. At all-time points, the transcript levels of MdBZR7 were higher (p < 0.05) in ABA-treated leaves than in the controls, but lower in MeJA-treated samples and repressed in SA-treated samples at p < 0.05 level.
The expression patterns of MdBZR1-8 were also identified after drought and salt treatments (Figure 7). The transcript levels of MdBZR1 were higher (p < 0.05) in the control group than in the drought-and salt-treated groups at most time points. In drought-treated leaves, MdBZR2 and MdBZR3 generally increased, but their expressions were down-regulated by salt. Drought and salt obviously repressed MdBZR4 at most time points. The expression of MdBZR5 was unaffected by drought but was down-regulated by salt on days 4, 6, and 8. The transcript levels of MdBZR6 and MdBZR7 in leaves were higher (p < 0.05) in the drought-treated group than in the control, whereas salt treatment decreased them. MdBZR8 expression increased on day 6 but decreased on day 8 in the drought-treated group.
In salt-treated trees, MdBZR8 was low at most time points, except for day 8.

Cloning and Subcellular Localization of Selected Apple BZR Genes
MdBZR genes (especially MdBZR1) may play important functions in both stem growth and abiotic stress responses. To further clarify their functions, MdBZR1, MdBZR4, MdBZR6, and MdBZR7 were cloned. Compared with its reference sequence, the cloned MdBZR1 was missing a 21 and 17 bp sequence at the 5 end and 3 end, respectively (Supplementary Figure S13A). There was a 36 bp insertion and one mismatch at the 5 end of the cloned MdBZR4 compared to its reference sequence (Supplementary Figure S13B). There were sequence deletions and mismatches at the 3 end of the MdBZR6 clone compared with its reference (Supplementary Figure S13C). Several sequence mismatches and a small insertion were identified in the MdBZR7 clone (Supplementary Figure S13D). The cloned MdBZR1 and MdBZR4 proteins were shorter and longer, respectively, than their references (Supplementary  Figures S14A,B). Like their nucleotide sequences, the protein FIGURE 5 | Effects of BR, IAA, and GA on apple BR downstream genes (except for MdBZRs). qRT-PCR data are shown relative to 0 day. Bars show mean ± standard error (n = 3). *, ** indicate significant differences at 0.05 and 0.01 level, respectively. sequences of cloned and reference MdBZR6 shared variation at the C-terminal (Supplementary Figure S14C). There were only nine mismatches in the protein sequence between the cloned and reference MdBZR7 (Supplementary Figure S14D). The accuracies of the identified genes in this study were confirmed from the above findings.
As TFs, MdBZRs were predicted to be located in the nucleus (Supplementary Table S4). To validate the results from bioinformatics analysis and lay a foundation for further identifying MdBZR1's function, 35:MdBZR1-GFP and 35:GFP were transiently transformed into tobacco leaf. As a positive control, GFP fusion protein was targeted in both the nucleus and cytoplasm. We found that MdBZR1-GFP was characteristically located in the nucleus (Figure 8).

DISCUSSION
BR influences various important agronomic traits in plants, such as agricultural and environmental adaptations, through downstream components (Xie et al., 2011;Miyaji et al., 2014). The functions of BR downstream components have been revealed in rice, thus demonstrating the feasibility of acquiring desirable traits in the breeding of other crops. However, information regarding corresponding Arabidopsis and Rosaceae genes is not available. Here, Arabidopsis and nine Rosaceae BR downstream genes were characterized, with information on chromosomal localizations, transmembrane helices, physicochemical properties, synteny analyses, phylogenetic relationships, Ka/Ks ratios, gene/protein structures, promoter analyses, and protein-protein interactions further clarified. Moreover, their expression patterns in various tissues and responses to different growth-related factors as well as stress signals were investigated in apple. These findings should provide rich information on Rosaceae BR downstream genes and support the utilization of these genes to improve the agronomic traits of apple trees.

Genome-Wide Identification and Characteristics Analysis of BR Downstream Genes
In rice, most BR downstream genes have been found with only one copy, with four OsBZRs and two OsRAVL1/OsRAV6s reported ( Figure 1A). In the current study, we identified several additional BR downstream genes in Rosaceae species (Figure 1B). For example, we identified eight MdBZR1/2, two MdILI1, two MdOSH1, two MdSPY, and three MdSMOS1 genes in apple. Six PcBZR1s or PcBZR2s, two PcILI1s, two PcOSH1s, two PcSMOS1s, four PcRAVL1s, or PcRAV6s, and three PcGSR1s were also identified. It is possible that gene duplications may result in gene expansions . Here, tandem and segmental duplications were found in almost all Rosaceae BR downstream genes ( Figure 1E and Supplementary Figure S3). Wholegenome duplication (WGD) has also been reported in many Rosaceae species during their domestication (Ou et al., 2019;FIGURE 6 | Expression analysis of apple BR downstream genes in YF-CF and CF/M9-CF/CF. DABB: days after bud break. Each value represents mean ± standard error of three biological replicates. Asterisks indicate significant differences at 0.05 level (*) and 0.01 level (**). Feng et al., 2020). Therefore, BR downstream gene expansion may have resulted from gene duplication and WGD. In addition, the Ka/Ks values of most duplicated or homologous gene pairs in each Rosaceae species were < 1 (Supplementary Tables S5, S6), indicating that these genes likely expanded under purifying selection. Based on their instability indices, most BR downstream genes were unstable ( Figure 1B). Introns play important roles in regulating gene expression (Wang et al., 2005). We found that one Rosaceae DLT or four Rosaceae RAVL1/RAV6 did not contain introns (Figure 1F), potentially highlighting their unique expression patterns.
Orthologous genes in pairs were identified between rice and each Rosaceae species (Supplementary Figure S4). The functions of the Rosaceae BR signaling downstream genes may be predicted according to their corresponding rice genes. For example, OsBZR1/2 participates in regulating culm length and lamina joint inclination (Bai et al., 2007), OsGSR1 is involved in controlling BR biosynthesis and plant size , and OsOSH1 maintains SAM, with osh1 mutants being dwarfs (Tsuda et al., 2011). Most BR signaling downstream proteins contained similar characteristic domains, although PaBZR1, PcBZR1, PcBZR2, PaLIC, and PaSMOS1-2 proteins all lacked one domain each (Supplementary Figure S6), which may result in altered function. For example, PaBZR1, PcBZR1, and PcBZR2 proteins all lacked a bHLH DNA-binding motif, which may lead to a loss in their DNA binding capacity. PaBZR1 protein lacked a serine-rich phosphorylation site (Supplementary Figure S6-1), which may result in its non-phosphorylated degradation in the cytoplasm (Yin et al., 2002). PaLIC protein lacked a serine-rich domain (Supplementary Figure S6-4), which could affect its phosphorylation status (Zhang et al., 2012). PaSMOS1-2 protein exhibited a shortened WTNF domain, which is a transcriptional activation domain (Supplementary Figure S6-6), which may reduce its transactivation activity (Hirano et al., 2017). In terms of secondary structures, most of these proteins were similar, except for PmBZR3, RcDLT, and RcRAVL1-1_RcRAV6-1 proteins (Supplementary Figure S7 and Supplementary Table S7), indicating that these three proteins may share different protein spatial configuration. PcBZR3 and PaLIC proteins were different from their homologous proteins in terms of length and subcellular localization (Supplementary Table S4), which may be indicative of their special characters. The parallel biological functions of these homologous genes can be predicted from their similarities in evolutionary relationships, synteny,   conserved domains, gene/protein structures, and chemical characteristics. However, several special genes should be studied in detail in the future.
Considering the pivotal roles of BZRs, DLT, LIC, ILI1s, and SMOS1s in stem elongation (Zhang et al., 2009;Qiao et al., 2017), we analyzed their expression patterns in YF-CF and CF/M9-CF/CF (Figure 6). With short nodes and fewer internodes, YF is a spur type of CF . Dwarf rootstock (M9) has been widely used to inhibit apple tree elongation . Here, the expression levels of MdBZR1-2, MdBZR4, MdBZR6, MdDLT, MdILI1-1, and MdSMOS1s were higher in CF than in YF at p < 0.05, whereas MdLIC was strongly expressed in YF (Figure 6). MdBZR1 and MdBZR4 were down-regulated but MdBZR5-8, MdDLT, MdILI1s, MdSMOS1-1, and MdSMOS1-3 were up-regulated in CF/CF. These results imply that these genes may have various effects on apple shoot growth, which should be confirmed in the future.
Previous studies have reported that BZR1-2 plays a role in stress resistance (Yang et al., 2016;Xie et al., 2019). Thus, we investigated the expression patterns of MdBZR genes under various stress treatments (Figure 7). ABA, MeJA, and SA are known for their roles in stress responses. We found that ABA induced MdBZRs (Figure 7), and ABA-responsive elements (ABREs) were recognized in all MdBZR genes (Supplementary Figure S8). AtBZR1 binds to the promoter of ABI5 and AtBZR2 down-regulates AtABI3 expression to repress ABA responses (Lopez-Molina et al., 2002;Ryu et al., 2014). In turn, ABA induces AtBZR1/2 and phosphorylates AtBZR2 protein (Zhang et al., 2009). In maize, ZmBZRs are in response to exogenous ABA, functionally (Yu et al., 2018). Furthermore, BvBZR genes respond to ABA-mediated processes . Here, in MeJA-treated leaves, MdBZR1-3, MdBZR5-6, and MdBZR8 were strongly expressed, while the expression levels of MdBZR4 and MdBZR7 were low. In addition, we identified at least one MeJA-responsive element in MdBZR1-8 (Supplementary Figure S8), which may explain the possible effects of MeJA on MdBZR1-8 expression. In sugar beet, certain BvBZR genes are up-and down-regulated by MeJA ). In the current study, the expression levels of all MdBZR genes, except MdBZR4, were low in the SAtreated samples (Figure 7). In the MdBZR2, MdBZR4-5, and MdBZR7-8 genes, several SA-related elements (TCA-element and SARE) were identified (Supplementary Figure S8). In Eucalyptus grandis, SA down-regulates EgBZR genes . Our results suggest that all MdBZRs, like AtBZR1/2, may potentially participate in ABA responses, several may positively or negatively regulate MeJA signals, and most may repress SA-induced defense reactions. We also found that drought stress repressed MdBZR1 and MdBZR4 expression levels, but up-regulated MdBZR2-3 and MdBZR6-7 (Figure 7). Salt treatment generally decreased the transcription levels of MdBZR1-8 (Figure 7). In Arabidopsis, rapeseed, and eucalyptus, BZR1/2 genes are induced or inhibited by drought or salt . Legume BZR genes show different expression patterns in response to drought treatment (Li et al., 2018). After salt treatment, expression levels of legume BZR and EgBZR genes are down-regulated (Li et al., 2018). In our study, several defense and stress elements, including drought-induced motifs, were detected in the MdBZR genes (Supplementary Figure S8). These results indicate that MdBZR genes play roles in drought and salt stress, like other BZR genes. The information on MdBZRs involving in stress adaption should be discussed in transgenic plants.

Diverse Roles of Apple BR Signaling Downstream Genes
Based on the above results, a Venn diagram was constructed to investigate the transcriptional responses of apple BR signaling downstream genes to BR, IAA, GA, CK, YF-CF, CF/M9-CF/CF, and stress (Figure 9). MdSPY2, MdRAVL1-1/MdRAV6-1, MdRAVL1-2/MdRAV6-2, MdSMOS1-2, and MdCSA were specifically controlled by BR. Six apple BR signaling downstream genes (MdSMOS1-1, MdSMOS1-3, MdGSR1, MdLIC, MdOSH1-2, and MdOSH1-3) were regulated by two factors. For example, MdSMOS1-1 and MdSMOS1-3 were regulated by BR and IAA; MdGSR1 was influenced by BR and GA; MdLIC was altered by BR and showed differential expression between YF and CF. MdILI1-2, MdBZR5, MdDLT, MdBZR3, and MdOSH1-1 were altered in response to three stimulations at the same time. MdBZR2, MdBZR6, and MdBZR7 were differentially expressed in four groups. MdBZR4, MdBZR8, and MdILI1-1 were up-or down-regulated by five treatments, simultaneously. MdBZR1 was consistently regulated by BR, IAA, GA, and stress treatments, and was also differentially expressed in YF-CF and CF/M9-CF/CF. The responses of apple downstream genes to different clues suggest their various functions, but their functional information should be obtained in a detailed study.

CONCLUSION
In this study, a total of 196 BR downstream genes were systematically identified from Arabidopsis and nine Rosaceae species. Their chromosomal locations, transmembrane helices, conserved domains, chemical characterizations, phylogenetic tree, gene/protein structures, synteny relationships, Ka/Ks ratios, promoter sequences, protein-protein interactions, and potential functions were analyzed. We found that gene duplication events (mainly segmental duplication and WGD) resulted in BR downstream gene expansion through purifying selection and diverse evolutionary patterns. Their functions were similar to their rice homologous genes. However, several genes (PaBZR1, PcBZR1-3, PmBZR3, RcBZR4, PaLIC, PaSMOS1-2, and all DLTs and RAVL1/RAV6s except for MdRAVL1/MdRAV6) may be species-specific, as predicted from the instability, transmembrane structure, subcellular localization, gene structure, conserved domain analyses, and responses to exogenous factors. Moreover, many apple BR downstream genes were predicted to control apple tree architecture via various hormones (BR, IAA, GA, and CK) and may play vital functions in shoot elongation of different varieties and grafting combinations. Furthermore, several MdBZR genes appeared to respond to stress-related signals. MdBZR1 was identified as a core gene due to its various roles, and MdBZR2-8, MdILI1-1, MdILI1-2, MdDLT, and MdOSH1-1 should also be noted.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JH designed the experiment. LZ, YY, SM, JC, and CY performed the experiments. LZ and YW analyzed the data. LZ drafted the manuscript. JH revised the final version of the manuscript. WW, MS, and XH assisted with revisions to the manuscript. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank WW, MS, and XH for critical revision of the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.