*2.2. Phylogenetic Analysis of BSK Family Genes*

*2.2. Phylogenetic Analysis of BSK Family Genes*  To explore phylogenetic relationships among the 143 BSK proteins identified in different plant species, a neighbor-joining phylogenetic tree was constructed. This analysis divided all BSK proteins into six major groups according to the bootstrap values and phylogenetic topology (Figure 2). Only one protein from maize, named Zm00008a014861\_T01, was in group I, which contained protein with longer amino acid sequence (743 aa) and an additional MSP (Major Sperm Protein) domain along with a conserved kinase domain and tetratricopeptide repeats (TPRs) (Figure S1). Group II was composed of two putative BSK proteins from the monocot seagrass *Z. marina*. The other four groups consisted of the BSK proteins largely from both dicot and monocot species. In particular, the BSK proteins from bryophytes including *P. patens* and *M. polymorpha,* were integrated into the group III. *A. comosus* and *Z. marina* belong to the angiosperm species prior to the split of eudicots and To explore phylogenetic relationships among the 143 BSK proteins identified in different plant species, a neighbor-joining phylogenetic tree was constructed. This analysis divided all BSK proteins into six major groups according to the bootstrap values and phylogenetic topology (Figure 2). Only one protein from maize, named Zm00008a014861\_T01, was in group I, which contained protein with longer amino acid sequence (743 aa) and an additional MSP (Major Sperm Protein) domain along with a conserved kinase domain and tetratricopeptide repeats (TPRs) (Figure S1). Group II was composed of two putative BSK proteins from the monocot seagrass *Z. marina*. The other four groups consisted of the BSK proteins largely from both dicot and monocot species. In particular, the BSK proteins from bryophytes including *P. patens* and *M. polymorpha,* were integrated into the group III. *A. comosus* and *Z. marina* belong to the angiosperm species prior to the split of eudicots and monocots.

monocots. The phylogenetic analysis showed that the BSK proteins from *A. comosus* (Aco018845.1, Aco011823.1, Aco014133, Aco010223.1, and Aco000489.1) divided the BSK proteins from dicots and The phylogenetic analysis showed that the BSK proteins from *A. comosus* (Aco018845.1, Aco011823.1, Aco014133, Aco010223.1, and Aco000489.1) divided the BSK proteins from dicots and monocots in each group. Moreover, other five BSK proteins from *Z. marina* (Zosma313g00120, Zosma1g02160, Zosma37g01020, Zosma41g01020, and Zosma7g01140) further divided the BSK members from dicots into smaller groups. These results could be considered as evidence for lineage-specific expansion of the BSK genes after the divergence of dicots and monocots. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 4 of 17 members from dicots into smaller groups. These results could be considered as evidence for lineage-specific expansion of the BSK genes after the divergence of dicots and monocots.

**Figure 2.** Neighbor-joining phylogenetic analysis of BSK genes. The gene tree was constructed using 143 BSK proteins and visualized using Figure Tree v1.4.2. Representative BSK proteins from different groups were marked with various colors. BSK proteins from dicot plants were in dark, and those from monocot were in pink. Blue and green ones were from *Zostera marina* and *Ananas comosus*, **Figure 2.** Neighbor-joining phylogenetic analysis of BSK genes. The gene tree was constructed using 143 BSK proteins and visualized using Figure Tree v1.4.2. Representative BSK proteins from different groups were marked with various colors. BSK proteins from dicot plants were in dark, and those from monocot were in pink. Blue and green ones were from *Zostera marina* and *Ananas comosus*, respectively.

#### respectively. *2.3. Genomic Structure and Conserved Motif Analysis of the BSK Gene Family*

*2.3. Genomic Structure and Conserved Motif Analysis of the BSK Gene Family*  The structure of the BSK genes was analyzed by comparing their full-length coding sequences (CDS) and the corresponding genomic DNA sequences using the software program GSDS 2.0 (http://gsds.cbi.pku.edu.cn/index.php). The exon/intron structure analysis showed that the structures of all BSK genes exhibited diversity both in their intron numbers and in their lengths (Figure 3). In particular, the intron numbers of the *Arabidopsis* BSK genes ranged between 6 and 9 (Figure 4A). Furthermore, we analyzed the different domain architectures, motif compositions, and gene structures of the BSK from *Arabidopsis*. The 12 BSK genes from *Arabidopsis* were located in all The structure of the BSK genes was analyzed by comparing their full-length coding sequences (CDS) and the corresponding genomic DNA sequences using the software program GSDS 2.0 (http: //gsds.cbi.pku.edu.cn/index.php). The exon/intron structure analysis showed that the structures of all BSK genes exhibited diversity both in their intron numbers and in their lengths (Figure 3). In particular, the intron numbers of the *Arabidopsis* BSK genes ranged between 6 and 9 (Figure 4A). Furthermore, we analyzed the different domain architectures, motif compositions, and gene structures of the BSK from *Arabidopsis*. The 12 BSK genes from *Arabidopsis* were located in all five chromosomes (Figure 4B). All BSK proteins had a putative kinase catalytic domain at the N-terminus and TPR

five chromosomes (Figure 4B). All BSK proteins had a putative kinase catalytic domain at the N-terminus and TPR repeats at the C-terminus, whereas BSK2 and BSK12 contained only two TPR repeats, one less than that from the other BSK members (Figure 4C). In addition, a glycine-rich repeats at the C-terminus, whereas *BSK2* and *BSK12* contained only two TPR repeats, one less than that from the other BSK members (Figure 4C). In addition, a glycine-rich region and unknown motif were found at the N-terminus in *BSK1* and *BSK11*, respectively. Previous studies showed that *BSK1* was phosphorylated by BRI1 at S-230, and a mutation in the ATP-binding site K-104 led to enzyme inactivity [10,17]. These conserved sites were also consistent with other BSK family members, except for *BSK12*. Moreover, most of the BSK proteins, except for *BSK9*, had a myristoylation site at the N-terminus (Figure S2). Other kinds of protein post-translational modifications, including sumoylation and ubiquitination, were also analyzed. It was found that all BSK family proteins possessed at least one sumoylation site and multiple ubiquitination sites, some of which were conserved (Figure S2). For instance, most of the ubiquitination sites were located in the kinase catalytic or TPR domain. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 5 of 17 site K-104 led to enzyme inactivity [10,17]. These conserved sites were also consistent with other BSK family members, except for BSK12. Moreover, most of the BSK proteins, except for BSK9, had a myristoylation site at the N-terminus (Figure S2). Other kinds of protein post-translational modifications, including sumoylation and ubiquitination, were also analyzed. It was found that all BSK family proteins possessed at least one sumoylation site and multiple ubiquitination sites, some of which were conserved (Figure S2). For instance, most of the ubiquitination sites were located in the kinase catalytic or TPR domain.

**Figure 3.** Phylogenetic relationship and exon-intron structure of BSK genes in representative plants. The BSKs were divided into five subgroups, which are indicated by different colors. The black boxes represent exons, solid lines represent introns, and green boxes represent untranslated regions (UTRs). The lengths of the BSK genes are indicated by horizontal lines (kb). **Figure 3.** Phylogenetic relationship and exon-intron structure of BSK genes in representative plants. The BSKs were divided into five subgroups, which are indicated by different colors. The black boxes represent exons, solid lines represent introns, and green boxes represent untranslated regions (UTRs). The lengths of the BSK genes are indicated by horizontal lines (kb).

*Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 6 of 17

**Figure 4.** The features of the BSK genes in *Arabidopsis*. (**A**) Exon-intron structures of the BSK family genes. (**B**) Chromosomal location of BSK genes. (**C**) Schematic diagram of *Arabidopsis* BSK genes. The putative domains or motifs were identified using the Conserved Domain Database (CDD), Pfam and Simple Modular Architecture Research Tool (SMART) databases with the default parameters. The kinase domain was labeled in pink. MY, myrislation; TPR: tetratricopeptide repeat; GYR: glycine rich domain; Un: unknown domain. **Figure 4.** The features of the BSK genes in *Arabidopsis*. (**A**) Exon-intron structures of the BSK family genes. (**B**) Chromosomal location of BSK genes. (**C**) Schematic diagram of *Arabidopsis* BSK genes. The putative domains or motifs were identified using the Conserved Domain Database (CDD), Pfam and Simple Modular Architecture Research Tool (SMART) databases with the default parameters. The kinase domain was labeled in pink. MY, myrislation; TPR: tetratricopeptide repeat; GYR: glycine rich domain; Un: unknown domain.

#### *2.4. Expression Profiles Analysis of BSK Family Genes from Arabidopsis 2.4. Expression Profiles Analysis of BSK Family Genes from Arabidopsis*

To compare the expression patterns of *Arabidopsis* BSK genes in different tissues or in responses to various abiotic stresses and hormones, the publicly available transcriptome sequencing and microarray datasets were analyzed. The results showed that most BSK members were constitutively expressed in different parts at different stages of plant development, while tissue-specific expression patterns were also observed (Figure 5A). For instance, *BSK1*, *BSK2*, *BSK3*, and *BSK5* were expressed at much lower levels in mature pollen, whereas *BSK11* and *BSK12* were more highly expressed. Among the 12 BSK genes, *BSK9* was specifically highly expressed in root, while *BSK5* was To compare the expression patterns of *Arabidopsis* BSK genes in different tissues or in responses to various abiotic stresses and hormones, the publicly available transcriptome sequencing and microarray datasets were analyzed. The results showed that most BSK members were constitutively expressed in different parts at different stages of plant development, while tissue-specific expression patterns were also observed (Figure 5A). For instance, *BSK1*, *BSK2*, *BSK3*, and *BSK5* were expressed at much lower levels in mature pollen, whereas *BSK11* and *BSK12* were more highly expressed. Among the 12 BSK genes, *BSK9* was specifically highly expressed in root, while *BSK5* was specifically expressed at very

specifically expressed at very low levels during seed development. A cell-specific expression pattern analysis for the BSK family genes showed that most BSK family members possessed different

one light:dark cycle.

low levels during seed development. A cell-specific expression pattern analysis for the BSK family genes showed that most BSK family members possessed different expression levels within three root zones (Figure 5B). *BSK11* and *BSK12* were undetectable, whereas *BSK9* and *BSK10* were found only in the differentiation zones. Moreover, the latest RNA-seq transcriptome data analysis showed that the expression of some BSK family genes were light-dependent (Figure 5C). *BSK1*, *BSK3*, *BSK5,* and *BSK8* were found to be highly expressed at night, while their expression was significantly depressed upon exposure to light. By contrast, *BSK6* was expressed in a rhythmic manner independent of light. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 7 of 17 expression levels within three root zones (Figure 5B). *BSK11* and *BSK12* were undetectable, whereas *BSK9* and *BSK10* were found only in the differentiation zones. Moreover, the latest RNA-seq transcriptome data analysis showed that the expression of some BSK family genes were light-dependent (Figure 5C). *BSK1*, *BSK3*, *BSK5,* and *BSK8* were found to be highly expressed at night, while their expression was significantly depressed upon exposure to light. By contrast, *BSK6* was expressed in a rhythmic manner independent of light.

Under various abiotic and biotic stress treatments, significant alterations in the transcript amounts of some BSK family members were observed at certain time points (Figure 6A). For instance, *BSK5* and *BSK9* showed significant responses to multiple stimuli, while at most times, these transcripts were down-regulated. *BSK6* was upregulated after one hour of exposure to salt or oxidative stress, whereas *BSK2* was strongly repressed after one hour of heat stress. Under phytohormone treatments, the transcript levels of *BSK5*, *BSK9,* and *BSK12* showed rapid responses with either up- or down-regulation (Figure 6B). Surprisingly, only *BSK12* was rapidly induced by BL, while the other three members (*BSK1*, *BSK2,* and *BSK3*) reported to be involved in BR signaling were insensitive to BL for three hours: this, however, remains to be further elucidated. Under various abiotic and biotic stress treatments, significant alterations in the transcript amounts of some BSK family members were observed at certain time points (Figure 6A). For instance, *BSK5* and *BSK9* showed significant responses to multiple stimuli, while at most times, these transcripts were down-regulated. *BSK6* was upregulated after one hour of exposure to salt or oxidative stress, whereas *BSK2* was strongly repressed after one hour of heat stress. Under phytohormone treatments, the transcript levels of *BSK5*, *BSK9,* and *BSK12* showed rapid responses with either up- or down-regulation (Figure 6B). Surprisingly, only *BSK12* was rapidly induced by BL, while the other three members (*BSK1*, *BSK2,* and *BSK3*) reported to be involved in BR signaling were insensitive to BL for three hours: this, however, remains to be further elucidated.

**Figure 5.** Development- and tissue-specific expression profiles of *Arabidopsis* BSK genes. (**A**) Heat map of the relative expression levels of all identified *Arabidopsis* BSK genes during development process. The colour bar represents the expression values. Samples from specific developing stages of *Arabidopsis* for gene expression analysis were highlighted with inverted triangles, such as roots from seedlings (1) and mature plants (2). (**B**) Expression profiles of BSK genes in specific root regions. DZ: division zone; EZ: elongation zone; MZ: mature zone. (C) Expression profiles of BSK genes across **Figure 5.** Development- and tissue-specific expression profiles of *Arabidopsis* BSK genes. (**A**) Heat map of the relative expression levels of all identified *Arabidopsis* BSK genes during development process. The colour bar represents the expression values. Samples from specific developing stages of *Arabidopsis* for gene expression analysis were highlighted with inverted triangles, such as roots from seedlings (1) and mature plants (2). (**B**) Expression profiles of BSK genes in specific root regions. DZ: division zone; EZ: elongation zone; MZ: mature zone. (**C**) Expression profiles of BSK genes across one light:dark cycle.

**Figure 6.** Expression patterns of all identified *Arabidopsis* BSK genes among hormone and abiotic stress. Public microarray data sets were obtained from BAR. Heat map representation the BSK genes under (**A**) six stress treatments, namely, cold, salt, drought, oxidative, wounding, and heat and (**B**) seven hormones, namely, ACC, zeatin, IAA, abscisic acid (ABA), methyl jasmonate (MeJA), gibberellic acid (GA-3), and brassinolide(BL). The colour bar represents the expression values. **Figure 6.** Expression patterns of all identified *Arabidopsis* BSK genes among hormone and abiotic stress. Public microarray data sets were obtained from BAR. Heat map representation the BSK genes under (**A**) six stress treatments, namely, cold, salt, drought, oxidative, wounding, and heat and (**B**) seven hormones, namely, ACC, zeatin, IAA, abscisic acid (ABA), methyl jasmonate (MeJA), gibberellic acid (GA-3), and brassinolide(BL). The colour bar represents the expression values.

#### *2.5. Diverse Gene Expression Patterns of BSK Family from* Arabidopsis *Using Reverse 2.5. Diverse Gene Expression Patterns of BSK Family from Arabidopsis Using Reverse Transcription-Polymerase Chain Reaction (RT-PCR)*

stresses in previous studies, but was observed to be strongly induced in our study.

*Transcription-Polymerase Chain Reaction (RT-PCR)*  The expression profiles of the representative BSK family members that demonstrated a diverse pattern of induced gene expression were validated by RT-PCR. The results showed that the expressions of some BSK genes were, indeed, regulated by stress or hormone treatments (Figure 7A). For instance, *BSK1* and *BSK3* were shown to significantly respond to salt, cold, and heat stresses after 6 or 24 hours. In response to heat stress, *BSK2* was first sharply inhibited (within 30 minutes of the stress treatment) and then highly accumulated within one hour, and then kept constant after incubation at the normal growth condition for an additional three hours. However, it was also found that stress responses in some BSK genes were inconstant to the trends observed with the microarray data. For example, *BSK2* was found to be significantly depressed by cold and salt The expression profiles of the representative BSK family members that demonstrated a diverse pattern of induced gene expression were validated by RT-PCR. The results showed that the expressions of some BSK genes were, indeed, regulated by stress or hormone treatments (Figure 7A). For instance, *BSK1* and *BSK3* were shown to significantly respond to salt, cold, and heat stresses after 6 or 24 h. In response to heat stress, *BSK2* was first sharply inhibited (within 30 min of the stress treatment) and then highly accumulated within one hour, and then kept constant after incubation at the normal growth condition for an additional three hours. However, it was also found that stress responses in some BSK genes were inconstant to the trends observed with the microarray data. For example, *BSK2* was found to be significantly depressed by cold and salt stresses in previous studies, but was observed to be strongly induced in our study.

The tissue-specific expressional patterns of these BSK family genes were also investigated (Figure 7B). The results showed that promising expressional levels of *BSK2*, *BSK3*, *BSK5*, *BSK7*, and

The tissue-specific expressional patterns of these BSK family genes were also investigated (Figure 7B). The results showed that promising expressional levels of *BSK2*, *BSK3*, *BSK5*, *BSK7*, and *BSK8* were found in leaves, flowers, fluorescents and siliques, which were consistent with the results from previous analysis. It also showed that all seven BSK genes were pronounced during reproductive growth. Our results further provided the information on the relative expression of these selected BSK genes. For instance, *BSK1* showed the strongest expression in flower and was dominantly expressed. Taken together, these results indicated the BSK genes are involved in multiple stress responses and developmental processes. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 9 of 17 *BSK8* were found in leaves, flowers, fluorescents and siliques, which were consistent with the results from previous analysis. It also showed that all seven BSK genes were pronounced during reproductive growth. Our results further provided the information on the relative expression of these selected BSK genes. For instance, *BSK1* showed the strongest expression in flower and was dominantly expressed. Taken together, these results indicated the BSK genes are involved in multiple stress responses and developmental processes.

**Figure 7.** Reverse transcription-polymerase chain reaction (RT-PCR) validation of gene expression patterns of representative BSK genes under stresses and organs. Expression profiles of selected BSK genes under stresses (**A**) and in different tissues (**B**). For various stresses treatment, 7-day-old seedlings were subjected to 10 μM ABA, 150 mM NaCl and 4 °C on plates, respectively. Heat stress treatment was monitored with four-week-old seedlings under 42 °C for one hour and then recovered at 22 °C for three hours. The intensity corresponding to the individual PCR band was calculated by Image J with reactions from untreated samples as 1.00. L: leaves; F: flowers; B: buds; S: siliques. **Figure 7.** Reverse transcription-polymerase chain reaction (RT-PCR) validation of gene expression patterns of representative BSK genes under stresses and organs. Expression profiles of selected BSKgenes under stresses (**A**) and in different tissues (**B**). For various stresses treatment, 7-day-old seedlingswere subjected to 10 <sup>µ</sup>M ABA, 150 mM NaCl and 4 ◦C on plates, respectively. Heat stress treatment was monitored with four-week-old seedlings under 42 ◦C for one hour and then recovered at 22 ◦C forthree hours. The intensity corresponding to the individual PCR band was calculated by Image J with reactions from untreated samples as 1.00. L: leaves; F: flowers; B: buds; S: siliques.

#### *2.6. Alternative Splicing (AS) Analysis for BSK Family Genes in Arabidopsis 2.6. Alternative Splicing (AS) Analysis for BSK Family Genes in Arabidopsis*

 The precursor mRNAs (pre-mRNAs) with introns can be spliced by alternative splicing (AS) to generate multiple mRNAs that are translated into different proteins. By analyzing the genome annotation obtained from the Phytozome v12.1 database (https://phytozome.jgi.doe.gov), we found that many BSK genes from different species possessed several transcripts resulting from AS (Figure S3). In particular, the BSK genes from *P. patens* underwent extensive AS events, resulting in an increase in the complexity of the BSK gene family transcriptome (Figure S3K). On the contrary, only one BSK gene was found to exhibit AS in *Arabidopsis* and rice. The precursor mRNAs (pre-mRNAs) with introns can be spliced by alternative splicing (AS) to generate multiple mRNAs that are translated into different proteins. By analyzing the genome annotation obtained from the Phytozome v12.1 database (https://phytozome.jgi.doe.gov), we found that many BSK genes from different species possessed several transcripts resulting from AS (Figure S3). In particular, the BSK genes from *P. patens* underwent extensive AS events, resulting in an increase in the complexity of the BSK gene family transcriptome (Figure S3K). On the contrary, only one BSK gene was found to exhibit AS in *Arabidopsis* and rice.

As annotated in the Phytozome database, *BSK5* had two spliced transcripts in *Arabidopsis* (Figure S3A). *BSK5.1* encodes a complete BSK protein, while the spliced transcript *BSK5.2* encodes an N-terminal-truncated BSK with the first intron retained from *BSK5.1*. To investigate the expression patterns of the two transcripts, specific primers were designed for RT-PCR (Figure 8A). Our results showed that *BSK5.1* was the dominant transcript, which was expressed in different tissues and responded to various environmental stimuli (Figure 8B,C). To our surprise, a third As annotated in the Phytozome database, *BSK5* had two spliced transcripts in *Arabidopsis* (Figure S3A). *BSK5.1* encodes a complete BSK protein, while the spliced transcript *BSK5.2* encodes an N-terminal-truncated BSK with the first intron retained from *BSK5.1*. To investigate the expression patterns of the two transcripts, specific primers were designed for RT-PCR (Figure 8A). Our results showed that *BSK5.1* was the dominant transcript, which was expressed in different tissues and responded to various environmental stimuli (Figure 8B,C). To our surprise, a third isoform *BSK5.3*

isoform *BSK5.3* was found in siliques (Figure 8B). A sequence analysis indicated two introns

In addition to *BSK5*, we also detected a novel splicing transcript of *BSK9* in flowers and buds (Figure 7B). Further analysis showed that the intron between exons 7 and 8 was retained in *BSK9*

retained in *BSK5.3*, implying that *BSK5* spliced in a tissue-dependent manner.

was found in siliques (Figure 8B). A sequence analysis indicated two introns retained in *BSK5.3*, implying that *BSK5* spliced in a tissue-dependent manner.

In addition to *BSK5*, we also detected a novel splicing transcript of *BSK9* in flowers and buds (Figure 7B). Further analysis showed that the intron between exons 7 and 8 was retained in *BSK9* (Figure 9). We further screened the ASIP (Alternative Splicing in Plants) database (http://www. plantgdb.org/ASIP/) for intron retention (IR) events among the BSK family genes in *Arabidopsis*, and there was only one record for *BSK7*. However, no IR event was detected in *BSK7* from the leaves under normal growth conditions (Figure 9). Instead, this IR event within *BSK7* was strongly induced by cold stress, and this same phenomenon was also found in *BSK5* and *BSK9*. Notably, *BSK5.3* was also found to be specifically induced by heat stress. These results indicated a novel post-transcriptional regulation pattern in the BSK genes. *Int. J. Mol. Sci.* **2019**, *20*, x FOR PEER REVIEW 10 of 17 (http://www.plantgdb.org/ASIP/) for intron retention (IR) events among the BSK family genes in *Arabidopsis,* and there was only one record for *BSK7*. However, no IR event was detected in *BSK7* from the leaves under normal growth conditions (Figure 9). Instead, this IR event within *BSK7* was strongly induced by cold stress, and this same phenomenon was also found in *BSK5* and *BSK9*. Notably, BSK5.3 was also found to be specifically induced by heat stress. These results indicated a novel post-transcriptional regulation pattern in the BSK genes.

**Figure 8.** Alternative splicing analysis of *BSK5* in *Arabidopsis*. (**A**) The scheme shows the primers (black arrow) used in the PCR. The black boxes represent exons, solid lines represent introns, and green boxes represent untranslated regions (UTRs). BSK5-ASF, BSK5.3-ASF: forward primers; BSK5-ASR, BSK5.2-ASR: reverse primers. (**B**) The mRNA corresponding to the three variants of *BSK5* in different tissues was detected by RT-PCR. The numbers indicate the lengths (bp) of the spliced transcripts and schematic diagram outlining the organization of the transcript variants were on the right. L: leaves; F: flowers; B: buds; S: siliques. (**C**) RT-PCR analysis of the expression profiles of *BSK5* variants under stress treatment. **Figure 8.** Alternative splicing analysis of *BSK5* in *Arabidopsis*. (**A**) The scheme shows the primers (black arrow) used in the PCR. The black boxes represent exons, solid lines represent introns, and green boxes represent untranslated regions (UTRs). BSK5-ASF, BSK5.3-ASF: forward primers; BSK5-ASR, BSK5.2-ASR: reverse primers. (**B**) The mRNA corresponding to the three variants of *BSK5* in different tissues was detected by RT-PCR. The numbers indicate the lengths (bp) of the spliced transcripts and schematic diagram outlining the organization of the transcript variants were on the right. L: leaves; F: flowers; B: buds; S: siliques. (**C**) RT-PCR analysis of the expression profiles of *BSK5* variants under stress treatment.

**Figure 8.** Alternative splicing analysis of *BSK5* in *Arabidopsis*. (**A**) The scheme shows the primers (black arrow) used in the PCR. The black boxes represent exons, solid lines represent introns, and green boxes represent untranslated regions (UTRs). BSK5-ASF, BSK5.3-ASF: forward primers; BSK5-ASR, BSK5.2-ASR: reverse primers. (**B**) The mRNA corresponding to the three variants of *BSK5* in different tissues was detected by RT-PCR. The numbers indicate the lengths (bp) of the spliced transcripts and schematic diagram outlining the organization of the transcript variants were

(http://www.plantgdb.org/ASIP/) for intron retention (IR) events among the BSK family genes in *Arabidopsis,* and there was only one record for *BSK7*. However, no IR event was detected in *BSK7* from the leaves under normal growth conditions (Figure 9). Instead, this IR event within *BSK7* was strongly induced by cold stress, and this same phenomenon was also found in *BSK5* and *BSK9*. Notably, BSK5.3 was also found to be specifically induced by heat stress. These results indicated a

novel post-transcriptional regulation pattern in the BSK genes.

**Figure 9.** Intron retention analysis of BSK genes under stresses. Plants were stressed for 6 h under ABA/salt/cold and one hour under heat. Samples were then collected for mRNA isolation. Specific primers in exons were designed for RT-PCR in BSK genes. The detected PCR products were outlined on the right.

#### **3. Discussion**

Several studies have reported the roles of BSK proteins from *Arabidopsis* and rice in BR signaling and immunity, as well as in abiotic stress responses [10–19]. However, the knowledge of BSK proteins in other plant species is still quite limited. In the present study, we identified a total of 143 BSK proteins from 17 plant species. Our analysis showed the evolutionary origin of BSK genes and BR receptor BRI1 or receptor-like genes in embryophytes, indicating the origin of BR signaling in plants (Figure 1A). A previous study reveals that plant BRI1 is highly conserved across taxa [20]. Taken together, the origin and development of the BR signaling system seems to be highly relevant with the evolution from aquatic to terrestrial plants, which have been observed in the ABA signaling system [21]. Bryophytes have an intermediate position between aquatic and terrestrial plants, and the establishment of the BR signaling pathway might have a great impact on the explanation of the movement of plants from water to land as an adaptation to environmental conditions. We also found evidence that genome duplication likely contributed the most to the expansion of the BSK gene family in many plant lineages. Prior to the split of eudicots and monocots, there was no evidence of whole-genome duplication (WGD). However, several rounds of WGD or triple genome duplication (WTD) have been reported in the 12 angiosperms studied here [22,23]. This might result in the expansion of the BSK genes from angiosperms. Indeed, we found that approximately two-fold expansion in the BSK genes was identified in the plants with WGD, including *A. thaliana*, *G.*, *P. trichocarpa*, and *Z. mays*. Furthermore, in *B. rapa*, we observed over three-fold BSK expansion, which was likely due to WTD.

In *Arabidopsis*, the BSK gene family includes 12 members, which share many similarities with each other with some exceptions. For instance, additional motifs were found in *BSK1* and *BSK11* (Figure 4C). The glycine-rich (GXXXG) motif, which provides conformational flexibility in the unstructured region, has been identified in many proteins, but has not yet been characterized in the BSK gene family [24]. It is thus, possible that this motif in *BSK1* might confer flexibility in the interaction of BSK with *BRI1* or in homo-dimmer conformation. Post-translational protein modifications, such as phosphorylation and myristoylation, have been reported in the BSK gene family. In the current study, no myristoylation site, which is critical for membrane association of BSK was found in *BSK9* [17]. A prediction of the subcellular localization of all BSK proteins from *Arabidopsis* showed that *BSK9*, as well as *BSK10* and *BSK11*, were absent from the plasma membrane, indicating that myristoylation might be sufficient, but not necessary, to anchor a protein to the membrane (Table S1). Further validation of the subcellular localization is required to support the accurate role of the myristoylation site. Other kinds

of post-translational protein modification in this gene family were also investigated. The prediction results showed that conserved sumoylation and ubiquitination sites were found in all BSK family proteins (Figure S2). Overall, these post-translational protein modification sites were mainly located in the N- and C-terminals, where the residues were probably the targets of a SUMO or ubiquitin complex. Many proteins are both sumoylated and ubiquitinated, and often at the same lysine residues. In most BSK proteins, some conserved residues were found to serve as both ubiquitination and sumoylation sites. These two protein modifications can act either synergistically or antagonistically, indicating the existence of a fine-tuned post-translational regulatory mechanism of the BSK family, which remains to be elucidated.

In the present study, comprehensive microarray and RNA-seq analyses of the *Arabidopsis* BSK proteins at several developmental stages showed that some of these BSK proteins might extensively function in many tissues (Figure 5). Our RT-PCR data also showed that the selected BSK genes, including *BSK3*, *BSK7*, and *BSK8*, were expressed in leaves, buds, wide open flowers, and siliques; this result was consistent with the previous findings (Figure 7B). Some BSK proteins seemed to function in specific tissues, likely resulting in their sub-functionalization. For instance, *BSK12* was upregulated only at the early stages of see development, whereas the expression of *BSK5* was found to be reduced gradually; this was also verified by our RT-PCR data. BRs were reported to regulate seed size and shape in *Arabidopsis* by BZR1-mediated activation and repression of a number of known regulators of seed size [25]. AtBZR1-like genes are highly conserved in angiosperms, and both soybean GmBZL2 rice OsBZR1 play an important role in seed set/size [26–29]. Thus, it could be possible that *BSK5* and *BSK12* function upstream of BZR1 in seed development. BRs were also reported to play important roles in regulating root meristem maintenance and root elongation [30]. Here, significantly high expression of *BSK9* was observed in the meristematic zone of roots (Figure 6C). It could be exciting to examine the function of *BSK9* in root meristem development. Moreover, the expression of *BSK1* was the lowest in mature pollen among the BSK family genes (Figure 5A). On the contrary, *BSK11* and *BSK12*, which were in the same sub-group, showed much stronger expression. BRs have been reported to promote pollen germination and growth [31]. Some BSK downstream genes, including *Arabidopsis* BES1 and rice OsBZR1, have been found to be involved in anther and pollen development, indicating that these BSK members might play a role upstream of them in male gametophyte development [29,32]. In addition, BRs also regulate photomorphogenesis and link with light signals through BES1 [33,34]. Here, we found that several BSK family members showed light-dependent transcriptional regulation. For instance, the expressions of *BSK1*, *BSK3*, *BSK5*, and *BSK8* were significantly upregulated in the dark and repressed by light (Figure 5C). It will be interesting to investigate the function of these BSK members upstream of BZR1 in photomorphogenesis, broadening the functions of the BSK gene family.

A systematical microarray analysis of the BSK proteins under distinct abiotic stresses (cold, salt, drought, et al) and phytohormones (ABA, 1-aminocyclopropane-1-carboxylic acid (ACC), indole-3-acetic acid (IAA), BR, etc)) based on the *Arabidopsis* eFP Browser database indicated that most BSK proteins were involved in these abiotic stress and hormone responses (Figure 6). For instance, *BSK2*, *BSK3*, *BSK5*, *BSK6*, and *BSK12* were involved in the salt-stress response. Among these BSK genes, most were repressed by salt stress, whereas *BSK6* was constantly upregulated. In addition, *BSK5*, *BSK9,* and *BSK12* showed broad responses to multiple stresses. In contrast, only a small number of BSK genes were regulated by hormones. To our surprise, most BSK members, including *BSK1*, showed no response to BRs in 3 hours, whereas *BSK12* was sharply induced in 30 minutes, indicating an unexpected role of *BSK12* through BR perception (Figure 6B). Our RT-PCR results further confirmed that several selected BSK genes indeed showed a response to various abiotic stresses and hormones, although there was some inconsistency with the microarray data (Figure 7A). For instance, *BSK1* was upregulated under salt and cold stresses after 24 hours. Under the same situations, *BSK2* and *BSK9* were also highly induced. By contrast, all BSK genes were repressed by heat stress and ABA the inhibited expression of *BSK3*, *BSK5* and *BSK8*. Further experiments are required to elucidate the important roles of these BSK genes in various stress responses.

Alternative splicing, an important modulator of gene expression, can increase proteome diversity and regulate mRNA levels by generating multiple transcripts from the same gene [35,36]. In plants, this post-transcriptional mechanism is noticeably induced by environmental stimuli and plays an important role in the regulation of gene expression for biotic/abiotic stresse responses as well as for plant growth and development [37–42]. The RNA-seq analysis has revealed that in *Arabidopsis* more than 61% of intron-containing genes and even more undergo alternative splicing under normal growth and stress conditions, respectively [31]. Notably, IR is the most prevalent form of AS in higher plants. For instance, the AS forms of *Arabidopsis* contain an unusually high fraction of retained introns (above 30%), and more than half of the AS events belong to IR in rice [43–45]. IR is mostly accepted on account of mis-splicing and thought to be non-functional because they are likely to result in nonsense-mediated decay [46]. However, several studies have highlighted the functional importance of intron-retaining mRNAs in plants. IR in the *Arabidopsis* INDETERMINATE DOMAIN 14 (IDD14) transcription factor gene generates a competitive inhibitor that modulates starch accumulation in response to cold stress [47]. IR in the 50 UTR of the Zinc-Induced Facilitator 2 gene (ZIF2) improves zinc tolerance in *Arabidopsis* [48]. *HAB1*, a group A protein type 2C phosphatase (PP2C), undergoes IR to produce a splice variant that plays opposing roles in fine-tuning ABA signaling [49]. Thus, IR should no longer be underestimated and exploring its roles in the development and/or stress response is of increasing importance. Here we found that the IR event in some *Arabidopsis* BSK genes was tissue-specific and responded to stress/hormone. The intron-containing form of *BSK9* preferred to express in development stages from buds to flowers, but was spliced completely in siliques. We also found that cold and heat stresses promoted the IR rate within the first two introns of *BSK5*, whereas ABA operated only in the first intron to generate and accumulate differently spliced transcripts. IR in *BSK7* and *BSK9* were also specifically induced by cold stress. Although the role of an individual transcript is largely unknown currently, these observations indicate the post-transcriptional regulation of these BSK genes in response to environmental stimuli. Future work will be performed to elucidate the mechanism and physiological significance of tissue-specific and stress-induced IR in these BSK genes.

Overall, our study represented diverse features of the BSK gene family in *Arabidopsis*, providing valuable insights into their important functions.

## **4. Materials and Methods**

## *4.1. Identification of BSK Family Genes and Phylogenetic Analysis of BSKs in Plants*

To identify ortholog(s) of AtBSK1 from *Arabidopsis* in the plant genome, the amino acid sequence of AtBSK1 was used as the seed sequence to perform a BLASTp search of the database whole genome sequences in the Phytozome database (https://phytozome.jgi.doe.gov). A similar method was applied on AtBRI1. Multiple protein sequence alignment for the deduced amino acid sequences of BSK proteins from different plant species was performed using the ClustalW software (http://www. clustal.org/clustal2/) with default settings. Then, the phylogenetic tree was constructed using MEGA software (Tokyo Metropolitan University, Tokyo, Japan) using the neighbor-joining (NJ) method and the bootstrap test carried out with 1000 replicates. The phylogenetic trees were visualized using FigTree (http://tree.bio.ed.ac.uk/software/figtree/).

#### *4.2. Conserved Domain Recognition and Gene Structure Analysis*

The CDS and the corresponding genomic sequences of the BSKs from the selected plant species were retrieved from the databases. Analysis of the exon-intron structures for the BSKs was carried out using the Gene Structure Display Server (GSDS v2.0, http://gsds.cbi.pku.edu.cn/index.php). Conserved domain for BSK proteins from *Arabidopsis* were identified using Conserved Domain Search Service (CD Search) from NCBI. The protein post-modification analysis was performed using SUMOgo (http://predictor.nchu.edu.tw/SUMOgo/) and UbiSite (http://csb.cse.yzu.edu.tw/UbiSite) with high

specificity level of 95%. Subcellular prediction for BSK family genes were performed using Cell eFP Browser (http://bar.utoronto.ca/cell\_efp/cgi-bin/cell\_efp.cgi).
