Genome-Wide Identification and Expression Analysis of bZIP Family Genes in Stevia rebaudiana Bertoni

1. Introduction

The basic leucine zipper (bZIP) is one of the most widely distributed and conserved transcription factor families among eukaryotes, playing a significant role in regulating plant growth and development. Extensive research has explored bZIP genes in various species, including Arabidopsis thaliana (78) [1], poplar (86) [2], Isatis indigotica (65) [3], licorice (Glycyrrhiza uralensis) (66) [4], safflower (Carthamus tinctorius) (52) [5]. The bZIP domain consists of basic amino acid region (N-X₇-R/K) and leucine zipper region that binds to the alkaline region [6]. The leucine (Leu) zipper region is located in the N-terminus, consists of highly conserved heptad repeats of leucine or other large hydrophobic amino acids, enabling the formation of homodimers and heterodimers among bZIP proteins through the leucine zipper. In plants, the bZIP proteins preferentially combined with ACGT core sequences [1]. The researchers initially divided the members of bZIP gene family in A. thaliana (Arabidopsis thaliana) into 10 subgroups (A, B, C, D, E, F, G, H, I, S) based on their conserved domain [6]. Then the AtbZIP genes further divided into 13 subgroups, added three groups (M, K and J) based on the original research [1].

In plants, members of the bZIP TF family participated in multiple biological processes, including tissues and organs development, responses to abiotic and biotic stresses, and secondary metabolism regulation. Co-expression of AtbZIP10, AtbZIP25 and ABI3 regulates seed specificity [7]. AtbZIP46 determines the number of floral organs and participates in stem and floral meristem expression [8]. HY5 and HYH mediates light response in A. thaliana as main regulators of photomorphogenesis [9]. AtbZIP17 functions as a transcriptional activator in the response to salt stress [10]. OsbZIP16 can reduce the sensitivity of overexpressed rice seedlings to abiotic stress during germination [11]. And, many studies have showed that bZIP transcription factors effectively regulate the biosynthesis of plant secondary metabolites, such as terpenoids, alkaloids, and flavonoids. AaHY5 regulates light-induced artemisinin biosynthesis by binding to G-box (CACGTG) site in the promoter of AaGSW1 [12]. CRISPR/Cas9-mediate knock-out of one allele of VvbZIP36 in grapevine promote anthocyanin accumulation [13]. MdHY5 promotes anthocyanin accumulation by binding to the G-box-2 site in the promoter of MdMYB10, thereby regulating the expression of downstream anthocyanin biosynthesis genes [14].

Stevia rebaudiana Bertoni, a perennial herb of the Asteraceae family, is renowned for being a valuable source of the tetracyclic diterpenoid derivative, natural sweetener steviol glycosides (SGs), known for their high sweetness and low calorie content. This plant holds not only high edible value but also significant medicinal importance, with steviol glycosides exhibiting various pharmacological effects such as anticancer, anticardiovascular disease, anti-inflammatory and antimicrobial. Despite its increasing application of, the production of S. rebaudiana (Stevia rebaudiana Bertoni) falls short in meeting the demand, necessitating urgent efforts to enhance production and ensure product quality by conducting insights into its developmental processes and responses to environmental factors.

Although there have been reports of studies sequencing the stevia genome [15], a systematic analysis of the complete genome information of S. rebaudiana is yet to be conducted. Hence, it becomes imperative to explore its molecular mechanism, network regulation and functional diversity. Specifically, the distribution, function and structure analysis of SrbZIP genes have not been reported. In this study, we identified the SrbZIP gene family members utilizing bioinformatics methods at the genome level. Additionally, we performed phylogenetic analysis, examined chromosomal locations and collinearity, analyzed motif compositions, cis-acting elements in the promoter, studied gene structures, conducted protein interaction network analysis and the explored the expression of SrbZIP genes in six S. rebaudiana cultivars. Moreover, we investigated the SrbZIP genes responses to various factors, including light, abiotic stress (low temperature, salt and drought), phytohormones (methyl jasmonate, gibberellic acid and salicylic acid), as well as in different tissues, utilizing qRT-PCR. The findings from this study offered valuable insights for further molecular genetic improvement of S. rebaudiana.

2. Results

2.1. Identification and Characterization of the SrbZIPs

A comprehensive search and structural domain determinatoin led to the identification of 105 bZIP family members in S. rebaudiana, named SrbZIP1-SrbZIP-

105. Analysis of the amino acid physicochemical properties showed that the amino acid length readings varied from 96 aa (SrbZIP13) to 704 aa (SrbZIP17), and the corresponding molecular weight size ranged from 10.88 KDa (SrbZIP13) to 77.31 KDa (SrbZIP17). The theoretical isoelectric point (PI) of SrbZIP proteins ranged from 4.64 (SrbZIP21) to 10.08 (SrbZIP101). The hydrophilic nature of the 105 sequences was indicated by their GRAVY values, with a maximum value of -0.24 (SrbZIP3) and a minimum value of -1.425 (SrbZIP27). The aliphatic index of SrbZIPs ranged from 37.60 (SrbZIP38) to 104.04 (SrbZIP88). Instability index analysis revealed that, except for five proteins (SrbZIP1, SrbZIP23, SrbZIP51, SrbZIP62, SrbZIP95), all SrbZIP proteins were predicted to be unstable, as their exceeded 40 [16]. The majority of SrbZIP proteins were located in the nucleus, while three SrbZIP proteins (SrbZIP4, SrbZIP25, SrbZIP27) were localized in cytoplasm, two SrbZIP proteins (SrbZIP12, SrbZIP17) in endoplasmic reticulum, and one SrbZIP protein (SrbZIP2) in chlorosomes through subcellular localization analysis. The proteins were neutral or alkaline accounted for 97.14% (102/105) of the total number (Table S1).

2.2. Classfication of SrbZIP Genes Based on Phylogram

To investigate the evolutionary connections and categorization of the bZIP family, we constructed a neighbor-joining (NJ) phylogenetic tree utilizing the MEGA-X software with 1000 bootstrap replicates The phylogenetic tree comprised 105 SrbZIPs and 73 AtbZIPs. Referring to the principles of classification in A. thaliana [1], S. rebaudiana bZIP genes were divided into 12 subgroups (SrbZIP-A, SrbZIP-B, SrbZIP-C, SrbZIP-D, SrbZIP-E, SrbZIP-F, SrbZIP-G, SrbZIP-H, SrbZIP-I, SrbZIP-J, SrbZIP-K and SrbZIP-S) (Figure 1). However, SrbZIP77 could not be aggregated into any subfamily, speculated that the structure of the SrbZIP77 has been greatly different in evolution or differed from the evolutionary direction from AtbZIPs. The SrbZIP genes classified in subgroup S contained the largest number of members, including 25 in S. rebaudiana and 15 in A. thaliana, while subgroup J and subgroup K contained the fewest members, both including only one in S. rebaudiana.

2.3. Chromosomal Location and Collinearity Analysis of SrbZIP Members

The distribution of the 105 SrbZIP genes across 11 chromosomes (Chr1-Chr11) and contigs in S. rebaudiana (Figure 2A). Notavly, Chromosome 2 carried the highest number of SrbZIP genes (21), whereas chromosome 5 and 10 had the fewest SrbZIP genes (4) (Figure 2A). Surprisingly, the number of SrbZIP genes on each chromosome was independent of the chromosome size.

During the process of gene evolution, gene duplication could effectively expand and obtain functional diversity via tandem duplication and segmental duplication [17]. To analyze the expansion pattern of SrbZIP genes, we analyzed their gene duplication events utilizing MCScanX. Thirty-two pairs of segmental duplication events were occurred on eleven chromosomes (Figure 2A; Table S2), as a major driver of bZIP genetic diversity and evolution of S. rebaudiana. Additionally, SrbZIP67 and SrbZIP68 genes experienced tandem duplication events (Table S3). .

In this study, interspecies collinearity analysis between SrbZIPs and AtbZIPs were constructed. The resuls showed conservative collinearity relationship with 62 orthologous pairs found between S. rebaudiana and A. thaliana, which on all the chromosomes (Figure 2B; Table S4).

2.4. Gene Structure and Conserved Motif Analysis of SrbZIPs

To investigate the gene structure and the motif composition of SrbZIPs, the phylogenetic tree, motif, protein domain and gene structure were built using TBtools (v1.047) (Figure 3). The analysis identified ten conserved motifs (Figure S1), with characteristic motif 1 regions found in all bZIP family members of S. rebaudiana (Figure 3B), indicating high conservation during the evolution of SrbZIP proteins. Motif 5 regions were detected in most bZIP subfamilies, except subfamily A and D. Gene structures within the same bZIP subfamily were mostly identical, although some exceptions were observed, for example, motif 9 regions were absent in seven bZIP members (SrbZIP47, SrbZIP23, SrbZIP24, SrbZIP71, SrbZIP72, SrbZIP67, SrbZIP68) of subgroup A compared to others. The C-terminal regions containing more motifs were crucial for bZIP dimerization and DNA binding . All SrbZIP proteins were found to have the bZIP superfamily domain (including BRLZ, MFMR and DOG1 domains), and six SrbZIP proteins contained the bZIP-HY5-like domain (SrbZIP5, SrbZIP10, SrbZIP13, SrbZIP70, SrbZIP54, SrbZIP96) (Figure 3C). Moreover, we investigated the exon-intron distribution patterns in SrbZIP genes (Figure 3D). Around 69.52% of the SrbZIP gene members contained 1-5 introns. The members of subgroup G, except SrbZIP16, SrbZIP17 and SrbZIP18, contained more than ten introns. A maximum of 13 introns were detected in SrbZIP87. Moreover, members of the same SrbZIP subgroup displayed similar intron/exon structures, for instance, most SrbZIP-S genes had no intron, while two members of subgroup H both contained two introns.

2.5. Cis-elements Analysis in SrbZIPs Promoter Regions

To understand the regulatory mechanisms underlying the response of SrbZIPs to abiotic or biotic stress, the promoter regions located within 2.0 kb upstream of SrbZIPs were submitted to PlantCARE database for predicting cis-elements. The cis- elements in the SrbZIP promoters fell into three main categories: development-related elements, defense and stress-responsive elements, and phytohormone-responsive elements (Table S5). Most of environmental stress-related elements were light-responsive elements, drought induction (MBS), low temperature responsive element (LTR) and antioxidant response element (ARE) (Figure 4). In addition, there were many phytohormone-responsive elements such as abscisic acid responsiveness (ABRE), MeJA-responsive elements (TGACG-motif), salicylic acid-responsive elements (TCA-element) were identified to exist in SrbZIP promoters. The members of subgroup G contained abundant light-responsive elements, SrbZIP93 promoter region contained 20 light-responsive elements. Subgroup A members, particularly SrbZIP61, exhibited an elevated number of abscisic acid responsive elements (ABREs),with 13 ABREs identified.

2.6. Protein Interaction Analysis Network of SrbZIP Proteins

To explore the interactions among SrbZIP proteins, we constructed a protein-protein interaction network (p-value: < 1.0e-16) relationship map of 105 SrbZIP proteins, the SrbZIPs were a complex family with 43 nodes and 172 edges. In the map of SrbZIP proteins interaction network, the results showed that 13 members (SrbZIP60, SrbZIP95, SrbZIP89, SrbZIP12, SrbZIP31, SrbZIP7, SrbZIP4, SrbZIP57, SrbZIP88, SrbZIP79, SrbZIP73, SrbZIP82 and SrbZIP69) were predicted to be involve in plant hormone signal transduction (ath04075), of which 3 members (SrbZIP7, SrbZIP57 and SrbZIP89) related to sugar and hormone signaling (WP3661). In addition, 4 SrbZIP proteins (SrbZIP83, SrbZIP88, SrbZIP7 and SrbZIP102) of the network of SrbZIP protein were predicted to be associated with seed development (WP2279). Base on their functions, SrbZIP proteins were categorized into four main groups: ABF2/ABF4/ABI5 (SrbZIP51/SrbZIP38/SrbZIP7), involved in phytohormon-e signaling: GBF1/GBF3/GBF4 (SrbZIP29/SrbZIP63/SrbZIP60), involved in environmental signaling; AREB3 (SrbZIP88), PAN (SrbZIP12), TGA1 (SrbZIP69), TGA4 (SrbZIP82), TGA7 (SrbZIP31), TGA9 (SrbZIP95), TGA10 (SrbZIP79) and HY5 (SrbZIP96), ivolved in cryptochrome signaling: and FD (SrbZIP23, SrbZIP24, SrbZIP47, SrbZIP67, SrbZIP68, SrbZIP71, SrbZIP72), promoting flowering (Figure 5; Table S6).

2.7. Expression Pattern and qRT-PCR Validation of SrbZIP Genes

To investigate the expression pattern of SrbZIP genes in leaves from different cultivars, the fragments per kilobase million (FPKM) values of the SrbZIP genes were retrieved from the transcriptome data described in our previous study [18], and the heatmaps were generated based on these values (Figure 6; Table S7).

SrbZIP30 and SrbZIP54 showed higher expression in 6 varieties of S. rebaudiana leaves, while SrbZIP30 in leaves had the highest expression. SrbZIP30, a member of subgroup S, which were functionally interrelated in plant starvation signaling [1]. The highly expression of SrbZIP30 may provide the most basic conditions in S. rebaudiana growth. Likewise, SrbZIP54 was classified in subgroup H, which were composed of only two members, elongated hypocotyl 5 (HY5) and HY5 homolog (HYH). Meanwhile, HY5 has proved conclusively to be one of the most predominant transcriptional regulator involved in promoting photomorphogenesis, pigment accumulation downstream of phytochromes and chloroplast development [19, 20]. Performing as a master regulator in coordinating light, developmental and environmental signaling result in its highly expressed.

In order to further reveal the biological function of SrbZIP genes, the expression patterns of SrbZIP genes in different tissues and conditions was analyze by qRT-PCR. We preliminary selected 10 SrbZIP genes from different bZIP subfamilies that expressed higher of ‘023’ S. rebaudiana cultivar. As we can see from Figure 6, SrbZIP30 (from S subfamily), SrbZIP54 (from H subfamily), SrbZIP100 (from I subfamily), SrbZIP60 (from A subfamily) SrbZIP63 (from G subfamily), SrbZIP70 (from K subfamily), SrbZIP32 (from I subfamily), SrbZIP9 (from D subfamily), SrbZIP21 (from D subfamily) and SrbZIP10 (from B subfamily) have higher expression in 8 bZIP subfamilies of S. rebaudiana leaves, in descending order. Comparative DEG analysis revealed 5 bZIP genes (p＜ 0.01) significantly differed between the ‘023’ cultivar bZIP genes and other 5 cultivars. Among them, SrbZIP45 (log2(FC) = 4.14 p < 0.01, from C subfamily), SrbZIP93 (log2(FC) = 3.016 p < 0.001, from G subfamily), SrbZIP102 (log2(FC) = 2.13 p < 0.001, from G subfamily), SrbZIP19 (log2(FC) =5.3 p < 0.01, from A subfamily) and SrbZIP104 (log2(FC) = 5.14 p < 0.001, from F subfamily) were abundantly represented upregulated genes. Totally, we selected 15 upregulated and highly expressed SrbZIP genes from 10 subfamily, which were gathered to analyze the expression patterns via qRT-PCR in different tissues (root, stem, leaf and flower).

The expression levels of most of the 15 SrbZIP genes varied considerably among different tissues (Figure 7). Among the examined 15 SrbZIP genes, 3 SrbZIP genes (SrbZIP63, SrbZIP54, SrbZIP70) showed higher expression levels in leaves, 3 SrbZIP genes (SrbZIP60, SrbZIP104, SrbZIP10) presented high expression in roots, 5 SrbZIP genes (SrbZIP93, SrbZIP100, SrbZIP102, SrbZIP9, SrbZIP19) presented high expression in steams, SrbZIP30 presented high expression in flowers. These findings indicated that these screened SrbZIP genes had tissue-specific expression potentials in S. rebaudiana.

2.8. Expression Pattern of the SrbZIP Genes in Response to Phytohormones and Abiotic Stresses.

Based on the factors that affect the synthesis of stevia glycosides [21], we found that stevia glycosides varied from various abiotic stresses and phytohormones treatment. The distribution of cis-elements in the upstream regions of promoters for the selected 15 SrbZIP genes were investigated, the analysis revealed that the upstream regions of promoters of the 15 SrbZIP genes contained abundant cis-elements, including photoresponsive elements, phytohormone responsive elements and stress responsive elements (Figure 8; Table S5). We postulated that these 15 SrbZIP genes were influenced by light, abiotic stress and phytohormone. To further explore whether SrbZIP genes has the potential to participate in stevia glycoside biosynthesis, we analyzed the gene expressions of 15 SrbZIP genes under different treatment conditions to assess their potentials, which could provide reference for further exploration of the regulation mechanism of stevia glycoside synthesis pathway.

Phytohormones such as salicylic acid (SA) [22], methyl jasmonate (MeJA) [23] and gibberellic acid [24], enhanced the stevioside content in S. rebaudiana cultivated in vitro. Hereby, we evaluated the sensitivity of 15 SrbZIP genes to exogenous SA, MeJA and gibberellic acid. As shown in Figure 9, 12 SrbZIP genes (SrbZIP102, SrbZIP9, SrbZIP63, SrbZIP70, SrbZIP104, SrbZIP21, SrbZIP60, SrbZIP93, SrbZIP100, SrbZIP10, SrbZIP32, SrbZIP54) of analyzed genes were up-regulated under SA treatment. SrbZIP100 showed a substantial 21-fold increase in expression level after 1-h treatments, compared to 0-h. SrbZIP54 attained its highest expression level after 6-h, with a 9-fold increase, compared to 0-h. The expression levels of 9 SrbZIP genes (SrbZIP102, SrbZIP9, SrbZIP63, SrbZIP70, SrbZIP93, SrbZIP100, SrbZIP10, SrbZIP32, SrbZIP54) increased firstly and then decreased, showed the highest expression levels after 1-12h. Additionally, the expression levels of 3 SrbZIP genes (SrbZIP104, SrbZIP21 and SrbZIP60) showed an increasing trend and reached the maximum after 48-h.

The expression levels of the 15 SrbZIP genes in S. rebaudiana leaf samples sprayed with MeJA was analyzed, the 15 SrbZIP genes were all affected and up-regulated by MeJA treatment (Figure 10). SrbZIP9 and SrbZIP63 exhibited highest induction after MeJA stress treatment, with a 48-fold and 56-fold increase, respectively. The SrbZIP54, SrbZIP100 and SrbZIP70 genes also showed significant increases, with peak expression levels observed after 6-h treatments, resulting in 18-fold, 30-fold and 56-fold increases, respectively. The expression of 6 SrbZIP genes (SrbZIP102, SrbZIP10, SrbZIP70, SrbZIP93, SrbZIP63 and SrbZIP60) generally remained up-regulated with the extension of treatment time, SrbZIP93 attained the maximum (19-fold) after 48-h treatment. Within 24 hours after MeJA treatment, the expression level of SrbZIP9 was remained 19-48 folds increase compared to the control. Meanwhile, the expression levels of SrbZIP45 and SrbZIP19 genes firstly increased and then decreased, attained the maximum after 1-h, resulting in a 17-fold and 15-fold increases, respectively, compared to 0-h.

The expression levels of SrbZIP30 and SrbZIP102 genes were initially down-regulated upon gibberellic acid treatment after 1-h and then recovered to the same level as the control (Figure 11). On the other hand, SrbZIP54 and SrbZIP32 genes showed moderate increases, reaching their highest expression levels after 6-h and 3-h treatments, with a 9-fold and 2-fold increase, respectively, compared to 0-h. The expression level of SrbZIP19 and SrbZIP63 genes began to be down-regulated after gibberellic acid treatment, and showed slight up-regulation after 24-h. In contrast, the expression level of SrbZIP70 and SrbZIP9 genes exhibited a down-regulated trend. SrbZIP45, SrbZIP100, SrbZIP104, SrbZIP60, SrbZIP21 and SrbZIP93 first decreased and then increased after gibberellic acid treatment, and showed significantly decreased after 24-h treatments. Overall, gibberellic acid treatment showed negative regulation to most of these 15 SrbZIP genes.

Previous research indicates that a long-day (16 h) photoperiod significantly increased S. rebaudiana leaf biomass and stevia glycoside (SGs) content [25]. Moreover, under varying light intensities, the biosynthetic genes of stevia glycoside exhibited differential expressions, leading to changes in stevioside (ST) and rebaudioside A (RA) contents [26]. Additionally, the maximum biomass accumulation in the callus culture of S. rebaudiana was induced by white light, compared to yellow, blue, green and red lights. We further explored the effect of the light on the 15 SrbZIP genes by respectively treating stevia leaves with white light and dark.

Except for SrbZIP70, SrbZIP63, SrbZIP21 and SrbZIP10 genes, other 11 SrbZIP genes were up-regulated to varying degrees after light treatment (Figure 12A). SrbZIP60 gene was rapidly induced and reached the maximum expression after 1-h light treatment, followed by a decline, but it remained upregulated compared to control. The expression levels of SrbZIP30, SrbZIP54 and SrbZIP19 genes increased firstly and then decreased, reached the maximum after 9-h and 6-h light treatment, resulting in 5-fold, 4-fold and 2-fold increases, respectively, compared to the control. Likewise, SrbZIP93, SrbZIP9 and SrbZIP32 genes displayed an initial up-regulation followed by down-regulation after light treatment, with down-regulation observed after 12-h treatments, then the expression levels increased slightly, indicating that these genes may be influenced by circadian rhythms. On the contrary, SrbZIP45 and SrbZIP102 genes first decreased and then increased after light treatment, in the early stage of light treatment, light negatively regulated the expression of the two SrbZIP genes. When the stevia leaves were shaded, some SrbZIP genes showed the opposite variation trend (Figure 12B). SrbZIP70, SrbZIP10, SrbZIP100, SrbZIP9, SrbZIP60 and SrbZIP63 genes were rapidly induced upon dark treatment, with the extension of treatment time, the up-regulated trend was maintained compared to 0-h, SrbZIP70 and SrbZIP10 genes showed highest expression with a 12-fold and 4-fold after 24-h dark treatment. Interestingly, whether SrbZIP60 gene was treated with light or dark, it was rapidly induced and maintained highly expressed after 1-h treatment. The expression levels of SrbZIP104 and SrbZIP30 continued to declining, while SrbZIP32 continued to increasing. SrbZIP19, SrbZIP93 and SrbZIP54 genes first decreased and then increased after dark treatment, and showed significantly decreased after 24-h treatments. Dark treatment negatively regulates the expression of these SrbZIP genes, which may trigger the defense mechanisms of stevia, it is also related to the circadian rhythm.

Inappropriate temperature conditions inhibit plant growth and development, which strongly influences secondary metabolism accumulation. And, a study has shown that all of fifteen genes of stevia glycoside (SGs) biosynthesis pathway were transcribed maximally at 25 ℃, while both low (15 ℃) and high temperatures (35 ℃) restrained their transcription [27]. As shown in Figure 13, the expression levels of all t 15 SrbZIP genes significantly down-regulated within a short time frame (1-3h). SrbZIP93, SrbZIP102, SrbZIP54, SrbZIP100 and SrbZIP32 genes have showed highest expression levels after 12-h low temperature treatment, and then decreased. Within a period of time, these genes may play a crucial role for the tolerance to low temperature stress in stevia, but inconspicuously recovered after sustained low temperature stress damage. Notably, the SrbZIP54 gene was dramatically increased and reach the maximum expression level after 12-h treatment, with a 34-fold increase, compared to the control, even when down-regulated, its expression level maintained higher than the control. The expression levels of the SrbZIP21 and SrbZIP45 gene were continuously down-regulated. The SrbZIP30, SrbZIP63, SrbZIP10, SrbZIP70, SrbZIP60, SrbZIP104, SrbZIP19 and SrbZIP9 genes firstly down-regulated and then up-regulated, reached maximum expression after 48-h low temperature treatment. Different SrbZIP genes played distinct physiological regulatory roles at different stress durations.

It was found that salt stress (NaCl) promoted the accumulation of stevia glycoside (SGs) significantly [28], and up-regulated several stevia glycoside (SGs) biosynthesis pathway genes (CMS, CMK, HDR and UGT76G1) [29]. To investigate the expression pattern of the 15 SrbZIP genes in stevia leaves under abiotic stress induced by the salts for different durations, we constructed heatmaps to evaluate their responsiveness to salt stress (Figure 14). Evidently, ten genes (SrbZIP60, SrbZIP70, SrbZIP19, SrbZIP9, SrbZIP93, SrbZIP21, SrbZIP10, SrbZIP30, SrbZIP45 and SrbZIP104) were significantly down-regulated after 1-h salt treatment. SrbZIP60, SrbZIP70 and SrbZIP19 genes were continuously down-regulated after salt stress. Following salt treatment for 6 h, two genes (SrbZIP54 and SrbZIP102) were significantly up-regulated by approximately 6 and 11 folds, respectively, and then decreased. These positively regulated SrbZIP genes under salt stress may be conducive to S. rebaudiana of improving the tolerance to abiotic stress and promoting plant growth and development. The expression of five genes (SrbZIP9, SrbZIP93, SrbZIP21, SrbZIP10 and SrbZIP30) firstly decreased and then increased, showed highest expression upon 6-h salt treatment. SrbZIP104 gene was significantly down-regulated, and then slightly recovered expression but lower than the control, while SrbZIP45 gradually reached to the 0-h expression level. In summary, most SrbZIP genes were early responsed to salt stress, while a few SrbZIP genes showed a slightly delayed response process.

A study has shown that the production of steviol glycosides (SGs) in callus and suspension culture of S. rebaudiana are enhanced [30]. The 15 SrbZIP genes showed different expression patterns under drought stresses caused by PEG treatment (Figure 15). With the prolonged stress time, the expression level of 10 SrbZIP genes (SrbZIP21, SrbZIP10, SrbZIP32, SrbZIP9, SrbZIP45, SrbZIP102, SrbZIP19, SrbZIP70, SrbZIP60 and SrbZIP104) showed an initial increase, followed by a decline, and reached significantly lower levels after 24-h treatment, finally increased again after 48-h treatment. The expression levels of SrbZIP21, SrbZIP63, SrbZIP30, SrbZIP60 and SrbZIP104 genes showed highest expression after 48-h, with a 13-fold, 9-fold, 10-fold, 33-fold and 4-fold increase, respectively, compared to the control. SrbZIP60 exhibited a rapid response to drought stress and maintained higher expression level than the control, except for 24-h. After 24-h treatment, other factors may come into play, with SrbZIP genes working in concert with other genes, enhancing the drought tolerance. The expression levels of SrbZIP100, SrbZIP54 and SrbZIP93 generally maintained down-regulated after PEG treatment. In summary, most SrbZIP genes were showed a fuctuating state for responsing to drought stress, while the drought stress negatively regulated few SrbZIP genes.

There were 3 and 12 SrbZIP genes that were differentially expressed under SA stress, respectively. 12 SrbZIP genes (SrbZIP102, SrbZIP9, SrbZIP63, SrbZIP70, SrbZIP104, SrbZIP21, SrbZIP60, SrbZIP93, SrbZIP100, SrbZIP10, SrbZIP32 and SrbZIP54) of analyzed genes were up-regulated upon SA treatment, while other 3 SrbZIP genes were down-regulated under SA stress condition. As for MeJA treatment, all 15 SrbZIP genes were affected and up-regulated, the speed of response to MeJA was different. SrbZIP70 and SrbZIP9 genes were dramatically increased and reach the maximum expression levels after 6-h MeJA treatment. The gibberellic acid treatment showed negative regulation to most genes of the 15 SrbZIP genes, but SrbZIP54, SrbZIP100 and SrbZIP60 genes were modestly increased and highly expressed upon 6-h and 24-h treatments, with 9-fold, 14-fold and 8-fold increases, respectively. Except for SrbZIP70, SrbZIP63, SrbZIP21 and SrbZIP10 genes, other 11 SrbZIP genes were up-regulated to varying degrees after light treatment, and SrbZIP54, SrbZIP30, SrbZIP60, SrbZIP93 and SrbZIP9 genes were rapidly induced. Among them, whether SrbZIP60 gene was treated with light or dark, it was rapidly induced and maintained a higher expression level. Additionally, SrbZIP54, SrbZIP70, SrbZIP60 and SrbZIP19 were up-regulated from the control, at highest expression levels with approximately 34 folds, 23 folds, 10 folds and 6 folds, respectively, after low temperature treatment. After 6-h salt treatment, two genes (SrbZIP54 and SrbZIP102 ) were dramatically up-regulated with a 6-fold and 11-fold, respectively, compared to 0-h. With prolonged salt stress, the expression level of 10 SrbZIP genes (SrbZIP104, SrbZIP60, SrbZIP70, SrbZIP30, SrbZIP19, SrbZIP102, SrbZIP9, SrbZIP32, SrbZIP10 and SrbZIP21) showed a total trend of increased, except for 24-h treatment, and SrbZIP60 extremely fast response to drought stress. These results elucidated the different responsive mechanisms of SrbZIPs under light-treatment, phytohormones treatment and abiotic stresses. The characteristics of SrbZIPs can be more effectively explored to tap their potential in the future.

2.9. Analysis of Terpenoid Synthesis-Related SrbZIP Genes Which Responsed to Light-treatment, Phytohormones Treatment and Abiotic Stresses

To investigate the potential SrbZIP genes associated with steviol glycosides (SGs) biosynthesis and focused on specific genes, we utilized the STRING software to analyze the protein interaction network between the 15 SrbZIP proteins and their homologous AtbZIP proteins in A. thaliana (Figure 16A; Table S8). Additionally, by searching relevant references, we selected candidate related proteins involved terpenoid synthesis, including AabZIP1 (GenBank: PWA69369.1) [31], AabZIP9 (GenBank: MG584701) [32], AaTGA6 (GenBank: MH201467) [33]and AaABF3 (GenBank: MH734935) [34], OsbZIP79 (Os11g0152700) [35] and OsTGAP1 (Os04g0637000) [36] (Figure 16B). We constructed a protein interaction network (p-value: < 1.0e-16) relationship map with 14 nodes and 25 edges (Figure 16A), it was found that there were 6 members (SrbZIP54 (HY5), SrbZIP60 (GBF4), SrbZIP9 (AHBP), SrbZIP70 (bZIP60), SrbZIP30 (bZIP2) and SrbZIP10 (bZIP17)) were predicted to be involved in positive regulation of metabolic process (GO: 0009893), of which 3 members (SrbZIP70, SrbZIP10 and SrbZIP9) and other 2 members (SrbZIP32 (VIP1) and SrbZIP45 (BZO2H3)) related to cellular response to stress (GO:0033554). Through the prediction of protein-protein interaction between SrbZIP proteins and 6 reported bZIP proteins involved in terpenoid synthesis, it was found that SrbZIP54 (HY5), SrbZIP63 (GBF3), SrbZIP32 (VIP1) and SrbZIP45 (BZO2H3) had interaction with terpenoid synthesis-related bZIP proteins (Figure 16B). SrbZIP54 (HY5) was predicted to interact with other proteins (COP1, PIF3, PHY and SPA1), functioning as the centre of a transcriptional network hub and a master regulator of light signal (Figure 16C). SrbZIP63 (GBF3) was predicted to be associated with most terpenoid synthesis-related proteins, indicating its potential involvement in terpenoid synthesis. Furthermore, SrbZIP63 (GBF3) was predicted to interact with abscisic acid responsive element binding factor 3 (ABF3) /abscisic acid insensitive (ABI5) participated in ABA signaling pathway, with FT and AP1 which involved in regulation of flowering and with HY5 which specifically bind G-box (Figure 16D). SrbZIP32 (VIP1) was predicted to interact with (Mitogen-activated protein kinases) MPKs (Figure 16E) which induced by stress, cytokines, plant hormones, growth factors and participated in plant signaling [37]. The results showed the potential function of SrbZIP32 (VIP1) participated in responding to stress and resisting disease. SrbZIP45 (BZO2H3) was predicted to interact with AtbZIP53 and AtbZIP1 (Figure 16F) which were pivotal regulatory factors involved in energy deficiency, sucrose starvation, and senescence-induced nutrient translocation [38] and with AKINBETA1 which involved in regulating nitrogen and sugar metabolism [39]. In addition, based on the expression and response speed of 15 SrbZIPs after various treatments, we selected other two SrbZIP proteins (SrbZIP60 and SrbZIP9) to constructed a protein-protein interaction network. SrbZIP60 (GBF4) was predicted to interact with (sucrose non-fermenting-1-related protein kinase) SNRK (Figure 16G) which were involved in different stress signal transduction pathways and participated in resisting to adverse environments [40], and with (open stomata 1) OST1 which involved in resisting to low temperature stress [41]. SrbZIP9 (AHBP-1B) was predicted to interact with (Arabidopsis nonexpresser of pr genes) NPRs (Figure 16H) which were positively regulated SA-dependent signaling pathway while negatively regulated JA-dependent signaling pathway [42,43]. These results suggested that SrbZIP60 (GBF4) and SrbZIP9 (AHBP-1B) may play a joint role in tolerance to abiotic stress.

3. Discussion

3.1. Classification and Gene Duplication of SrbZIPs

Stevia rebaudiana Bertoni, a perennial herb renowned for its high sweetness and low-calorie characteristics. Its leaves contain natural sweetener steviol glycosides (SGs) with numerous health benefts, which is a tetracyclic diterpenoid derivative. The plant bZIP transcription factors have proved to involved in multiple biological processes such as tissues and organs development, responses to abiotic and biotic stress, and secondary metabolism [44]. Although there have been reports of studies sequencing the stevia genome [15], bZIP family genes in S. rebaudiana have not been comprehensively identified and their roles in terpenoid synthesis are unclear. A total of 105 SrbZIP genes were discovered in the S. rebaudiana genome by a homology search in this study. Comparing with other reported species, such as Arabidopsis (78) [1], rice (Oryza sativa) (89) [45], sorghum (92) [46], soybean (Glycine max) (160) [47], poplar (86) [2], licorice (Glycyrrhiza uralensis) (66) [4], safflower (Carthamus tinctorius) (52) [5] and Marijuana (Cannabis sativa) (58) [48], S. rebaudiana had the second highest number of bZIP genes. Phylogenetic relationship of SrbZIPs showed that bZIP family genes in S. rebaudiana were classified into 12 subgroups, which were similarly observed in those of C. tinctorius (Carthamus tinctorius) [5]. The SrbZIP genes classified in subgroup S contained the most members, while subgroup J and K contained the fewest members. However, SrbZIP77 could not be aggregated into any subfamily, speculated that the structure of the SrbZIP77 has been greatly different in evolution or differed from the evolutionary direction from A. thaliana bZIP family. Additionally, we observed that the 105 SrbZIP genes were distributed across 11 chromosomes in S. rebaudiana with no apparent pattern (Figure 2A). SrbZIP genes distributed on chromosome 2 shared the biggest number of members (21). This uneven distribution may be attributed to differences in chromosome size and structure. The variation in the bZIP gene family among different plant species could be influenced by the genome size or the gene duplication events during evolution [49]. Gene duplication generally occurs through two major mechanisms: segmental and tandem duplications, representing distinct evolutionary patterns in plants. Chromosome rearrangement generates numerous duplicated chromosomal blocks that result in segmental duplication events occur [50]. It has been revealed that segmental duplication during evolution resulted in the expansion of various gene families. In our study, the results suggested that the expansion of bZIP gene family in S. rebaudiana was mainly attributed to segmental duplication. Thirty-two pairs of gene segmental duplications were detected on eleven chromosomes and one pair (SrbZIP67 and SrbZIP68) with evidence for tandem duplication (Figure 2A; Table S2). Furthermore, that there was great collinearity between S. rebaudiana and A. thaliana, with sixty-two orthologous pairs were identified to exist which on all the chromosomes (Figure 2B; Table S4).

3.2. Structure Characteristic and Function Prediction of SrbZIPs

The SrbZIP gene family exhibited highly conservation, with all members sharing characteristic motif 1 regions, and the conserved motifs within the same subfamily were closely related to the gene structure (Figure 3). Ten motifs were identified in SrbZIP gene family through the motif analysis, named motif 1 to motif 10 (Figure S1), the overall compositions of motifs were similar within the same subfamily. The numbers and lengths of exons varied considerably among different SrbZIP subfamilies, resulting in gene lengths ranging from 96 to 704 amino acids (Figure 3D; Table S1). In eukaryotes, the variation in intron locations and numbers may account for some specific gene functions and evolutionary trajectories, and intron evolution is often associated with gene segmental duplication [51]. Our study revealed that 19% of the total SrbZIP genes were intronless (Figure 3D), which were similar to those of rice (15.3%) [45], tomato (17.4%) [52] and poplar (22%) [2]. Subgroup S contained 25 members of SrbZIP genes, among them, 18 members of S subgroup were intronless, similarly, 19 members classified in S subgroup in poplar, with18 members containing no introns [2]. These results indicated that the gene structures of SrbZIP genes which belonged to the S subgroup were highly conserved across different species.

Cis-elements are essential for the transcriptional regulation of gene expression in response to abiotic stresses [53]. Photoresponsive elements, followed by ABA, MeJA, gibberellin and SA responsive elements, were among the major cis-elements in the SrbZIP promoter regions (Figure 4). In addition, the defense and stress-responsive elements (low-temperature, salt and drought stress, etc) were frequently discovered in SrbZIP promoters. These cis-regulatory elements in SrbZIP promoters provided foundational evidences for the functional relevence of SrbZIP genes in regulating the growth and response to various abiotic/biotic stresses.

Proteins interact with each other to participate in intracellular/ intercellular signal transduction energy and material metabolism, regulation of gene expression, and regulation of cell cycle. Through the homology search of AtbZIPs, some orthologous pairs were identified to exist between SrbZIPs and AtbZIPs, which revealed the potential functional relevence of most SrbZIP proteins. To further tap their potential, we carried out an interaction prediction network analysis utilizing the STRING database, with searching multiple sequences (Figure 5). SrbZIP proteins were mainly classified into four major categories: involved in phytohormone signaling, involved in environmental signaling, involved in cryptochrome signaling and promoted flowering.

3.3. Expression Patterns of SrbZIP Genes and Light, Phytohormone and Abiotic Stress Response

Each SrbZIP gene’s profile of expression in six S. rebaudiana cultivars was elucidated. The SrbZIP genes were also verified the differential expression among stevia cultivars. Comparative DEG analysis revealed 5 bZIP genes (SrbZIP45, SrbZIP93, SrbZIP102, SrbZIP19 and SrbZIP104) (p＜0.01) significantly differed between the ‘023’ cultivar bZIP genes and other 5 cultivars. In addition, 10 SrbZIP genes from different bZIP subfamilies that expressed higher of ‘023’ S. rebaudiana cultivar (Figure 6). Most genes of the 15 SrbZIP members were differently expressed among various tissues (Figure 7). SrbZIP63, SrbZIP54 and SrbZIP70 showed higher expression levels in leaves, which may be related to leaf development. SrbZIP30 presented high expression in flowers, which may relate to regulation of flowering.

We found that stevia glycosides varied among various abiotic stresses and phytohormones treatment [21]. The expression of 15 SrbZIP genes under different treatment conditions was analyzed, which provided reference for further exploration of the regulation mechanism of stevia glycoside synthesis pathway. Except for drought stress and dark treatment, SrbZIP54 exhibited up-regulation in all phytohormone treatments and abiotic stresses, with expression levels consistently higher than the control (Figure 9, 10, 11, 12, 13, 14, 15), attributed to the abundant cis-elements in the upstream regions of promoters of SrbZIP54 gene in S. rebaudiana, including stress response, phytohormone response and light response elements (Figure 8). SrbZIP54 (HY5) was predicted to interact with other proteins (PIF3, COP1, PHY and SPA1), suggesting it functions as a central transcriptional network hub and a master regulator of light signaling (Figure 16C). ELONGATED HYPOCOTYL5 (HY5) inhibits hypocotyl growth and lateral root development, and promotes pigment accumulation in a light dependant manner in A. thaliana [54]. HY5 is a major regulator regulating plant growth and development, including cell elongation and proliferation, chloroplast development, pigment accumulation and nutrient assimilation [54,55]. Recently, the role of HY5 in other aspects has also been reported, such as hormone signal transduction, plant defense and temperature response [56,57]. Most importantly, HY5 also regulated the terpenoid synthesis, reported in Artemisia (AaHY5) [12] and in A. thaliana (AtHY5) [58]. It provided evidence for HY5 may act as a potential transcriptional regulator for stevia biosynthesis. SrbZIP63 (GBF3) was up-regulated in drought treatments (Figure 15), the result can also be verified in A. thaliana. Overexpression of AtGBF3 in A. thaliana enhanced the tolerance to osmotic stress, salinity and drought stress [59]. SrbZIP63 was up-regulated after dark treatment, on the contrary, was down-regulated after light treatment, which consistent with the results reported in the study of GBF3 mRNA predominantly in dark grown leaves and in roots [60]. In the same subfamily (G-subfamily), SrbZIP93 and SrbZIP102 showed different expression trends after light treatment, with SrbZIP102 preferred dark treatment. Meanwhile, SrbZIP102 and SrbZIP63 showed same expression trends after phytohormone treatment (SA, MeJA and GA). It indicated bZIP genes from the same subfamily showed similar functions with similar structure. SrbZIP63 (GBF3) was also predicted the potential function in terpenoid synthesis, abscisic acid response and regulation of flowering (Figure 16D). SrbZIP32 (VIP1) was upregulated with approximately 10-fold higher expression under drought stress, compared with the control (Figure 15). It was predicted that interacted with mitogen-activated protein kinase 3 (MPK3) (Figure 16E), indeed it is, VIP1 protein was reported not only participated in agrobacterium-mediated plant transformation, but also related to plant immune signal transduction pathways with phosphorylated by MPK3 [61]. The I-subfamily bZIP genes in A. thaliana were reported the potential value of disease resistance and stress resistance [1]. Although only a subset of I-subfamily bZIP genes has been analyzed yet, providing insight into potential functional relevence of I-subfamily bZIPs to defense and stress response, regulation of cell cycle and various developmental aspects. SrbZIP45 (BZO2H3) was predicted to interact with AabZIP9 (Figure 16B), however, it lackluster responsed to most stress treatments, with the expression level was down-regulated under treatments. Above all, SrbZIP54, SrbZIP63 and SrbZIP32 showed promise as potential candidates for enhancing abiotic stress tolerance and secondary metabolite production in S. rebaudiana through genetic improvement.

4. Materials and methods

4.1. Plant Materials

Six S. rebaudiana cultivars ( ‘GX’, ‘GP’, ‘1188’, ‘110’, ‘11-14’ and ‘023’ ) were obtained from the growth room of Sichuan Agricultural University (Chengdu, China) and propagated utilizing a 1:1 mixture of nutrient soil and vermiculite. The ‘023’ S. rebaudiana cultivar contains more types of stevia glycosides, we chose it as the experimental material. To explore the discrepancy of expression patterns of SrbZIP genes in various tissues, different S. rebaudiana tissue samples (root, stem, leaf and flower) were collected in triplicate. The samples weres snap frozen in liquid nitrogen and then stored at - 80 ℃.

Furthermore, we evaluated the expression levels of SrbZIP genes under light, abiotic stress (low temperature (4 ℃), salt (NaCl) and drought (PEG) and phytohormones treatment (salicylic acid, Methyl jasmonate and gibberellic acid), some stevia plants were treated by foliar spraying with 200 mM NaCl [28], 5% PEG 4000 [30], 2 mM SA [22], 0.1 mM MeJA [23], 2 mg·L^-1 GA [24], and samples were collected at indicated time points after treatment. For light treatment, partial stevia plants were grown under normal light illumination (150 µm^-2s^-1), while another set of plants under dark environment. For low temperature treatment, stevia plants were grown at 4 ℃, and other growth conditions were the same [27]. Three biological replicates were used for each treatment, the 0-h treatment was set as the control. The samples weres snap frozen in liquid nitrogen and then stored at - 80 ℃.

4.2. Data Sources

The complete genome sequences for Stevia rebaudiana Bertoni were retrieved from NCBI database (http://www.ncbi.nlm.nih.gov/, accessed on 14 February 2023). The RNA transcriptome databasebdownloaded from our previous study [18].The AtbZIP protein sequences were obtained from PlantTFDB database (https://planttfdb.cbi.pku.edu.cn, accessed on 1 March 2023).

4.3. Identification of S. rebaudiana bZIP Gene Family

First, to classify the SrbZIP gene family members, the known AtbZIP protein sequences were regarded as queriers for performing BLASTp searches in the S. rebaudiana genome (E-value < 10⁻¹⁴). Second, the profile hidden Markov models (HMMs) of the bZIP domain (PF00170, PF07716, PF12498 and PF03131) were downloaded from the Pfam database (http://pfam.xfam.org/, accessed on 2 March 2023). Then, we further identified and screened the conserved domains utilizing NCBI Conserved Domain Search (CD-Search) Tool (https://www.ncbi.nlm.nih.gov/S-tructure/cdd/wrpsb.cgi, accessed on 5 March 2023). After removing incomplete and redundant sequences, members of the SrbZIP gene family were obtained. The identified amino acid sequence were submitted to ProtParam tool of the ExPASy website (http://web.expasy.org/protparam/, accessed on 5 March 2023) for preditionsof physicochemical properties, and to Wolf Psort website (https://www.genscript.com/wolf-psort.html, accessed on 5 March 2023) for protein subcellular localization prediction.

4.4. Stevia rebaudiana bZIP Chromosomal Location and Collinearity Analysis

By analyzing Stevia rebaudiana genome annotation information and using TBtools (https://github.com/CJ-Chen/TBtools-Manual) [62,63], the sequence of the SrbZIP genes was mapped to individual chromosomes. To investigate gene duplication events, MCScanX software was employed [64]. The collinearity relationships between the SrbZIP genes and related genes were visually analyze by TBtools software.

4.5. Phylogenetic Analysis of SrbZIP Genes

To explore the phylogenetic relationships between the bZIP proteins of S. rebaudiana and A. thaliana, we used the neighbor-joining (NJ) method constructed in MEGA-X software [65]. We constructed a multiple alignment of the amino acid sequences by using ClustalW. The phylogenetic analysis included 1000 bootstrap replicates and visualized using the Interactive Tree of Life (iTol) server (https://itol.e-mbl.de/, accessed on 20 March 2023).

4.6. Gene Structure, Conserved Motif, and Cis-acting Element Analysis

The conserved motifs of bZIP proteins were analyzed utilizing the MEME Suite tools (https://meme-suite.org/meme/tools/meme, accessed on 6 May 2023) [66].The upstream 2.0 kb sequence of each bZIP gene was defined as the promoter region, and the promoter sequence was retrieved and submitted to the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools//plantcare/html/, accessed on 6 May 2023). Visualization of the conserved domains, exons, introns, motifs and cis- elements in the bZIP genes was performed utilizing the TBtools software package [63].

4.7. Protein–Protein Interaction Analysis Network of SrbZIPs

SrbZIP proteins with highest homology to AtbZIP proteins were retrieved by OrthoVenn3 software (E-value: 1×10^-2, inflation value: 1.5) (https://orthovenn3.bioin-fotoolkits.net, accessed on 7 May 2023). And we carried out an interaction prediction network map utilizing the STRING database (https://string-db.org, accessed on 7 May 2023) [67].

4.8. Stevia rebaudiana bZIP Expression Pattern Analysis

The fragments per kilobase million (FPKM) value of each SrbZIP gene was retrieved from the RNA-seq data described in our previous study [18]. TBtools software was used to analyze the expression pattern by constructing heatmaps of SrbZIP gene expression. The candidate SrbZIP genes were selected for qRT-PCR validation analysis, and the stevia β-actin [AF548026] gene was applied as a reference sequence. Relative transcript abundances were calculated via the 2^-ΔΔCt method [68]. The primers were listed in Table S9. All data were normalized based on setting up the relative expression level, the expression level of 0-point treatments for phytohormones and abiotic stress was set as 1.

5. Conclusion

Based on the complete genome data of Stevia rebaudiana Bertoni, a total of 105 SrbZIP genes were aggregated into 12 subfamilies bioinformatics method. Conserved domains, exons, introns and motifs indicated similarities within bZIP clusters. Segmental duplications were predominately responsible for the expansion of SrbZIP gene family. The cis-element analysis evaluated the potential multiple roles of SrbZIPs in light, phytohormone treatment and abiotic stresses. Utilizing the RNA-seq data from leaves of six S. rebaudiana cultivars, real-time fluorescence PCR data of SrbZIP gene family in different tissues and expression patterns of the selected 15 SrbZIP genes in response to light, phytohormones and abiotic stresses, 6 SrbZIP (SrbZIP54, SrbZIP63, SrbZIP32, SrbZIP45, SrbZIP60 and SrbZIP9) genes were further screened and 3 SrbZIP genes (SrbZIP54, SrbZIP63 and SrbZIP32) were finally identified as highly induced by one or more of these factors, potentially influencing stress-responsive terpenoid synthesis in S. rebaudiana. The protein-protein interaction network further expand the potential functions of SrbZIP genes. Our results paved the way for future functional studies to explore the roles of SrbZIP genes in stress-responsive terpenoid synthesis in S. rebaudiana.