2.1 Phylogenetic relationship of BSKs from spinach and other plant species
To evaluate the molecular evolution relationship of BSKs from spinach with other plant species, an un- rooted neighbor- joining (NJ) phylogenetic tree was constructed using MEGA-X with the maximum likelihood method (Fig. 1). The protein sequences of seven BSKs from spinach, seven from sugar beet, 12 from Arabidopsis, 18 from common tobacco, 12 from woodland tobacco, five from rice, and 11 from maize were clustered in the phylogenetic tree (Fig. 1 and Additional file 1). Among them, 12 Arabidopsis AtBSKs [11, 17-32, 40, 41], five rice OsBSKs [15, 16], two woodland tobacco NsBSKs (i.e. NsBSK1 and NsBSK3) , a common tobacco NtBSK2 , and nine maize ZmBSKs  have been reported previously, and all the other sequences of BSKs were obtained by searching against the genomics database using BLASTp. Totally, 72 proteins were grouped into seven classes (Class I to Class VII), and seven SoBSKs distributed into four classes (Class I, Class II, Class V, and Class VII) (Fig. 1). In Class I, SoBSK1 and SoBSK7 were clustered into two sub-branches with their homologous sugar beet BvBSK1 and BvBSK7, respectively, and the BSKs from monocotyledonous rice and maize were grouped into another sub-branches. In Class II, SoBSK2, BvBSK2, and AtBSK2 were grouped into one subgroup, and the BSKs from rice (OsBSK2) and maize (ZmBSK2 and ZmBSK3) were grouped into another subgroups. In Class V, SoBSK3 was grouped into a subgroup with BvBSK3 and AtBSK5. Three NsBSKs and six NtBSKs were grouped into another three subgroups. Class VII included four main subgroups, Subgroup I (contained AtBSK8 and AtBSK7), Subgroup II (i.e. SoBSK5, SoBSK6, BvBSK5, and BvBSK6), Subgroup III (i.e. SoBSK4, BvBSK4, and two AtBSKs), and Subgroup IV contained two NsBSK2 and five NtBSKs. Besides, Class III, Class IV and Class VI did not contain any SoBSKs. Herein, AtBSK6, AtBSK9, and AtBSK10 were clustered into Class IV. Two OsBSK3 and five ZmBSKs were clustered into Class III (i.e. OsBSK4, ZmBSK8, and ZmBSK9) and Class VI (i.e. OsBSK3, ZmBSK5, ZmBSK6-1, and ZmBSK6-2). The phylogenetic analysis indicated that seven SoBSKs had close homology with BvBSKs from sugar beet, and all these reported BSKs were highly conserved among monocotyledons and dicotyledons.
2.2 Chromosomal distribution analysis of SoBSKs
The physical map positions of seven SoBSKs on chromosome were analyzed using MapGene2Chrom web v2 (http://mg2c.iask.in/mg2c_v2.0/) against spinach genome. Seven SoBSKs were distributed nonrandomly and unevenly across four out of six spinach chromosomes, except for chromosome 1 and 3 (Fig. 2A). SoBSK4, SoBSK7, and SoBSK1 were localized on chromosome 2, 4, and 5, respectively, while chromosome 6 harbored SoBSK2, SoBSK5, and SoBSK6. In addition, SoBSK3 was mapped onto an unanchored scaffold (Scf_01469) according to the draft spinach Sp75 genome, which was assembled with low coverage rate (47%) .
2.3 Collinearity analysis of BSKs between spinach and other plants
The collinearity relationship was constructed between spinach and three plant species (i.e. Arabidopsis, sugar beet, and maize) (Fig. 2B). Collinearity analysis showed that SoBSK genes had homologous genes in Arabidopsis, sugar beet, and maize, of which Arabidopsis had six homologous gene pairs located in chromosome 1, 4, and 5, followed by sugar beet (five gene pairs in chromosome 1, 7, and 8, as well as an unanchored scaffold chr5_nu.sca001), and maize (only one gene pair in chromosome 10) (Fig. 2B and Additional file 2). Interestingly, SoBSK4 (LOC110789384) had homologous gene pairs with all three plant species, especially three gene pairs in Arabidopsis, including AtBSK3 (AT4G00710), AtBSK4 (AT1G01740), and AtBSK7 (AT1G63500). SoBSK6 (LOC110778653) and SoBSK7 (LOC110788203) had homologous gene pairs with both sugar beet and Arabidopsis. In addition, both SoBSK4 and SoBSK6 had close relationship with AtBSK7. SoBSK6 was also homologous with AtBSK8 (AT5G41260). These results indicated that SoBSKs had close phylogenetic relationship with dicotyledons rather than monocotyledons.
2.4 Gene structure and protein function domain analysis of SoBSKs
Phylogenetic analysis indicated that SoBSKs were grouped into three groups according to their amino acid sequences. SoBSK5, SoBSK6, and SoBSK4 were included in Group I, SoBSK1, SoBSK7, and SoBSK2 were classed into Group II, while SoBSK3 was a distinct Group III (Fig. 3A).
To predict the evolutionary feature and functional diversification, exon-intron organization of seven SoBSKs was analyzed according to their gene sequences in SpinachBase (http://www.spinachbase.org/) using TBtools. SoBSK4 contained eleven introns, SoBSK2, SoBSK3, and SoBSK5 had nine introns, SoBSK1 and SoBSK7 had eight introns, while SoBSK6 had four introns (Fig. 3B). The diversity of various introns and exons implied the possible alternative splicing in SoBSKs was critical for spinach development and stress response, which were reported in Arabidopsis BSKs . In addition, SoBSK4 and SoBSK6 only had 3’ UTR, SoBSK5 only had 5’UTR, while the remaining four SoBSKs had both 5’UTR and 3’UTR (Fig. 3B).
The conserved motifs were detected using MEME online search software (https://meme-suite.org/meme/tools/meme) (Fig. 3C and Additional file 3). Thirteen conserved motifs were found in the six SoBSKs, except for SoBSK6 which only had seven motifs (Motif 1, 2, 4, 5, 8, 9, and 11) (Fig. 3C). Among 13 motifs, seven motifs (Motif 1, 3, 5, 6, 10, 11, and 12) were localized in the kinase domain, while three motifs (Motif 2, 4, and 9) were distributed in the TPR domain. Importantly, more than half amino acid sites (26 out of 50) in Motif 1 were highly conserved, among of which contained four conserved S sites (S22, S26, S36, and S43) (Fig. 3D). Besides, the 2nd serine in Motif 1 also conserved in SoBSKs, except for SoBSK6 (Fig. 3D and G). The Arabidopsis homologous site of the 2nd serine in Motif 1 (S230 in AtBSK1 and S210 in AtBSK3) had been reported to a phosphorylated site by AtBRI1 in the BR signaling pathway [11, 16]. Moreover, the glycine (G) 16, G28, and G37 in Motif 1 were also conserved in all SoBSKs (Fig. 3D and G). The homologous sites in AtBSK3 of G16 and G28 (G226 and G238 in AtBSK3) can interfere other BSKs in BR signaling and protein stability of AtBSK3, respectively . Motif 6 has 26 conserved amino acids out of 50 amino acids, including two lysines (K10 and K29 in Motif 6), and K29 was the key site for ATP binding and BSK kinase activity (Fig. 3E and G) . In addition, the 2nd glycine in Motif 13 was highly conserved myristoylation site for membrane localization of BSKs (Fig. 3F and G) [16, 18].
2.5 Analyses of conserved function domain and modification sites in SoBSKs
A common feature of BSKs from Arabidopsis and rice is the presence of a N-terminal kinase domain and two or three C-terminal TPR domains, which can interact directly with each other for autoregulating its kinase activity [11, 16]. In spinach, we also found the conserved N-terminal kinase domain and C-terminal TPR domains in seven SoBSKs using the NCBI web CD-search tool (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) with the default parameters (Fig. 3G). SoBSK2 and SoBSK3 contained two C-terminal TPR domains, but all the other SoBSKs had three C-terminal TPR domains (Fig. 3G).
We performed the prediction of some conserved amino acid sites for protein PTM. The sequence alignment with these reported functional sites in Arabidopsis [11, 18, 21, 26, 27, 29, 41] and rice , the GPS-SNO 1.0  with default parameters, and UbiComb (http://nsclbio.jbnu.ac.kr/tools/UbiComb/) with a threshold value of 0.5  were used for the predictions of phosphorylation, myristoylation, S-nitrosylation, and ubiquitination, respectively (Fig. 3G). Except for SoBSK6, all the other SoBSKs had a serine site and a lysine site in the kinase domain (Fig. 3G), which were the homologous sites with the conserved phosphorylated site by BRI1 and ATP- binding site for kinase activity determination in Arabidopsis, respectively [11, 26]. Similarly, most SoBSKs, except for SoBSK6, had a glycine site for myristoylation modulating its plasma membrane localization, and more than one cystines for S-nitrosylation in response to nitric oxide signal  (Fig. 3G). In addition, SoBSK1 and SoBSK5 had a predicted ubiquitination site (K100 in SoBSK1 and K104 in SoBSK5) in the kinase domain (Fig. 3G), which were supposed to be a critical site for BSK degradation through ubiquitin-proteasome system. This indicated that these conserved PTM sites determined subcellular localization, redox regulation, and protein turnover of SoBSKs.
2.6 Analysis of cis-acting elements in SoBSK promoters
In order to explore the potential function of SoBSKs, cis-acting elements in seven SoBSK promoter sequences of 2,000 bp upstream from the translation start site of individual gene were predicted using the PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), and the prevalence distribution of cis-acting elements for stress- and hormone- responses were schematically depicted (Fig. 4A). A total of 294 cis-acting elements were found in SoBSKs involved in diverse abiotic stresses, such as light (80 elements), drought (74), heat (29), anaerobism (22), wound (11), cold (2), and phytohormone (76) (Fig. 4B and Additional file 4). Among them, the elements in response to drought (i.e. DRE, MYB, MRS, MYC, and MBS), light (e.g., Box 4, G-box, and GT1-motif), heat (i.e. CCAAT-box, STRE, and AT-rich), and phytohormone accounted for the main parts in each SoBSK, and some of these elements were conserved across seven SoBSKs, which implied that SoBSKs could play important roles in response to drought, light, heat, and various hormones. Importantly, 76 hormone-responsive elements included 28 for ABA (i.e. ABRE), 23 for methyl jasmonate (MeJA) (i.e. TGACG-motif, CGTCA-motif, and JERE), nine for salicylic acid (SA) (i.e. TCA-element), eight for ethylene (i.e. ERE), three for gibberellin acid (GA) (i.e. P-box and TATC-box), and five for auxin (i.e. TGA-element and AuxRR-core) (Fig. 4B). Collectively, these conserved elements in the promoter region suggested that SoBSKs were pivotal for these abiotic stresses and phytohormone stimuli.
2.7 Localization of SoBSK1 and SoBSK6 proteins
The subcellular localization of SoBSKs was performed using tobacco leaves. SoBSK1 and SoBSK6 were fused respectively in frame to the 3’-terminus of the GFP reporter gene under the control of CaMV 35S promoter. The recombinant SoBSK1-GFP, SoBSK6-GFP, and GFP alone were transiently expressed in the epidermal cells of tobacco leaves. The 35S::GFP signals were distributed in the cytosol, nucleus, and PM in tobacco epidermal cells, while the signal from SoBSK1-GFP was concentrated in the PM (Fig. 5). This suggested that the SoBSK1-GFP recombinant protein localized on the cell membrane. Interestingly, unlike SoBSK1-GFP, SoBSK6-GFP recombinant protein localized not only in the cell membrane, but also in the cytoplasm and nucleus, because SoBSK6 was lack of the N-terminal myristoylation site for potential membrane localization . This indicated that the truncated N-terminus without myristoylation site in SoBSK6 led to the missing of exclusive localization to cell periphery.
2.8 Organ-specific expression analysis of SoBSKs
To better understand organ- specific expression of seven SoBSK genes, total RNA from roots, stems, and three pairs of alternate leaves (i.e. the 1st, 2nd, and 3rd leaves) from heat- sensitive variety Sp73 and heat- tolerant variety Sp75 were prepared. The transcription level of each SoBSK gene was evaluated for three replicates using reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) analysis. In heat- sensitive variety Sp73, SoBSK6 in stems, SoBSK3 in the 1st leaves, SoBSK4 in the 2nd leaves, and SoBSK2 in the 3rd leaves had the highest expression levels (Fig. 6). However, in heat- tolerant variety Sp75, SoBSK5 in stems, the 1st and 2nd leaves, as well as SoBSK4 in the 3rd leaves exhibited the highest expressions (Fig. 6). On the other hand, for each gene in variety Sp73, SoBSK1, SoBSK2, and SoBSK4 in the 3rd leaves, SoBSK3 in the 1st leaves, SoBSK5 and SoBSK7 in the roots, as well as SoBSK6 in the stems exhibited the highest expression (Fig. 6). While, in variety Sp75, it was SoBSK1 and SoBSK4 in the 3rd leaves, SoBSK2, SoBSK3, SoBSK5, and SoBSK6 in the 2nd leaves, as well as SoBSK7 in the roots that had the highest expression (Fig. 6). These suggested that, except of SoBSK7 in roots, the other SoBSKs mainly functioned in leaves rather than stems and roots in different varieties of spinach (Fig. 6).
2.9 ABA- and brassinolide (BL)- responses of SoBSKs in leaves
BSKs are critical components in BR signaling pathways . We found multiple SoBSKs had ABRE cis-acting elements in their promoters for ABA signaling (Fig. 4). This implied that the cross-talk of ABA and BR signaling would happen in regulating plant growth, development, and stress response [46, 47]. Therefore, we evaluated the expression profiles of seven SoBSKs in leaves under exogenous BL and ABA treatments using RT-qPCR analysis (Fig. 7).
For BL treatment, all the SoBSKs in variety Sp73 were induced more than 1.5- fold at 2 and 4 h. SoBSK1, SoBSK2, SoBSK3, SoBSK4, SoBSK5, and SoBSK7 reached their highest expression levels at 4 h, while SoBSK6 maintained a higher level with more than 2- fold during whole process (Fig. 7). Except of SoBSK1 and SoBSK3, the other SoBSKs in variety Sp73 were increased at 24 h (Fig. 7). However, in variety Sp75, SoBSK2, SoBSK3, and SoBSK5 were apparently induced at 2, 4, and 12 h, but reduced at 24 h (Fig. 7). Besides, SoBSK1 and SoBSK4 were increased, but SoBSK6 and SoBSK7 were decreased in variety Sp75 under most conditions (Fig. 7). These indicated that, although the seven SoBSKs were all involved in BL- response, most SoBSKs, except of SoBSK6 and SoBSK7, were BL- induced in both Sp73 and Sp75 at the early stages (2 and 4 h) of BL treatment, and most SoBSKs exhibited significantly different even opposite responses in Sp73 and Sp75 under 24 h of BL treatment (Fig. 7). Interestingly, no SoBSKs were obviously ABA- induced in Sp73, and most of them were reduced in leaves under various exogenous ABA treatment conditions (Fig. 7). On the contrary, in Sp75, most SoBSKs in leaves were ABA- induced at 2 and 4 h, but ABA- reduced at 24h (Fig. 7). All these suggested that heat- sensitive variety Sp73 and heat- tolerant variety Sp75 probably employed different BSK- mediated mechanisms in response to exogenous BR and ABA treatments.
2.10 Temperature stress- response of SoBSKs in leaves
Temperature extremes inhibit seed germination and reduce plant growth and reproduction . We have predicted that 29 heat- responsive cis-acting elements and two cold- responsive cis-acting elements existed in SoBSKs promoters (Fig. 4). To evaluate the function of these SoBSKs in temperature response, the expression profiles of SoBSK genes in the heat- sensitive variety Sp73 and heat- tolerant variety Sp75 under heat and cold stresses were validated using RT-qPCR analysis. Most SoBSKs were significantly altered in both varieties under certain temperature stress condition when compared with the normal condition (0 h) (Fig. 8). Importantly, the opposite temperature- responsive expression patterns of most SoBSKs were observed in varieties of Sp73 and Sp75. Five SoBSKs, except of SoBSK1 and SoBSK6, were reduced in Sp73, but induced in Sp75 at 4 and 12 h of cold and heat stresses (Fig. 8). Under heat conditions, SoBSK2 and SoBSK4 were obviously decreased in Sp73, but SoBSK5 and SoBSK7 were increased at 12 and 24 h. However, five SoBSKs (i.e. SoBSK2, SoBSK3, SoBSK4, SoBSK5, and SoBSK7) were heat- induced in Sp75 under each time point, which was represented by the significantly heat- induced SoBSK2 (Fig. 8). The similar expression patterns were found under cold conditions. Five SoBSKs (i.e. SoBSK2, SoBSK3, SoBSK4, SoBSK5, and SoBSK7) were cold- increased in Sp75, but cold- decreased in Sp73, while two SoBSKs (i.e. SoBSK1 and SoBSK6) were slightly induced in Sp73, but reduced in Sp75 (Fig. 8).
2.11 Heat- responsive genes in BR- and ABA- signaling in leaves
The expression levels of six BR signaling- related genes and 16 ABA signaling- related genes were detected in leaves from variety Sp73 and variety Sp75 under heat treatments for 0, 2, 4, 12, and 24 h using RNA-sequencing assay (Fig. 9). Two genes (BRL2 and BSK2) for BR signal perception were slightly heat- induced in variety Sp75 at 2 and 4 h, but significantly decreased in variety Sp73 at 12 and 24 h. However, the BKI1 for inhibition of BR binding with BRI1 was heat- reduced in Sp73 and Sp75 at 12 and 24 h (Fig. 9). Besides, three genes encoding transcription factors in the BR- mediated signaling pathway, including BZR1, BEH2, and BEH4, were induced in variety Sp75 under heat treatment (Fig. 9). This implied that BR signaling was probably employed in variety Sp75 for facilitating its heat tolerance.
The other 16 heat-altered genes were involved in ABA biosynthesis, transport, signal transduction, and transcriptional regulation (Fig. 9). Five genes for ABA biosynthesis, including two ZEPs, ABA4, NCED1, and NCED5, were reduced in variety Sp73, and four of them were also reduced in variety Sp75, except of NCED1 induced at 2 and 24 h (Fig. 9). Besides, ABCC2 for ABA transport from cytoplasm to vacuole was induced, but ABA efflux transporter ABCG25 and influx transporter ABCG22 for long-distance ABA transport were reduced in both Sp73 and Sp75 (Fig. 9). These indicated that the biosynthesis and transport of ABA were heat- disturbed in spinach, while some specific pathways probably were enhanced in heat- tolerant variety Sp75. Significantly, the expressions of four ABA receptor encoding genes (i.e. PYR1, PYL4, PYL5, and PYL9) were changed under heat treatment. Among them, PYL5 was significantly heat- induced in variety Sp75 and Sp73, and PYL4 was induced in Sp75 at 2 and 4 h (Fig. 9). Interestingly, PP2C24 and PP2C51, the negative regulation factor of ABA signaling, were obviously reduced in Sp75 under heat stress. Besides, the expressions of transcription factor ABI5 and AREB3 were altered, and AREB3 was increased in both Sp73 and Sp75 under heat treatments (Fig. 9). All these implied that variety Sp75 exhibited more enhanced ABA signaling to regulate heat- responsive process when compared with variety Sp73.
2.12 Heat-responsive phenotypes of SoBSK1-overexpressing Arabidopsis plants
The SoBSK1 expression pattern was distinct with other SoBSKs in heat- tolerant variety Sp75 in response to heat treatment (Fig. 8). To investigate the function of SoBSK1 on heat tolerance, we constructed overexpressing SoBSK1 transgenic Arabidopsis plants under the control of the strong constitutive CaMV 35S promoter, and compared the survival rates of wild-type (WT), bsk134678 mutant, and SoBSK1- overexpressed seedlings under heat treatment (Fig. 10). Three- day- old WT, bsk134678 mutant, and SoBSK1- overexpressed seedlings grew under 22 ℃ for five days as control, or treated under 43 ℃ for 4 h, and following under 22 ℃ for five days as heat treatments (Fig. 10A and B). There were no differences among WT, bsk134678 mutant, and SoBSK1- overexpressed seedlings under control condition, but they exhibited obvious different phenotypes under heat treatments (Fig. 10A). The bsk134678 mutant seedlings showed higher survival rate (92%), which was about 3.7- fold and 4.7- fold higher than WT and SoBSK1- overexpressed seedlings, respectively (Fig. 10A and C). While WT and SoBSK1- overexpressed seedlings exhibited similar phenotype, and no significantly different survival rate was observed (Fig. 10C). These implied that SoBSK1 probably regulated plant heat tolerance as a negative factor and the members of BSK family had function redundance.