Genome-wide analysis of bZIP transcription factors and their expression patterns in response to methyl jasmonate and low-temperature stresses in Platycodon grandiflorus

Plants and materials

The experimental materials were P. grandiflorus tissue culture seedlings grown in a growth chamber in 1/2 MS culture medium for 3 months (16 h/8 h light/dark cycle, 25 °C). The P. grandiflorus specimens were kept in the Plant Herbarium at Anhui University of Chinese Medicine (Hefei, Anhui, China). The P. grandiflorus seeds came from the School of Pharmacy, Anhui University of Chinese Medicine, specimen number 2020070 (Su et al., 2021). Healthy seedlings were randomly selected and divided into stress treatment groups in a climate chamber. In the MeJA treatment, MeJA (0.5 mM) was sprayed on P. grandiflorus seedlings with Tween 20 as a dispersant. P. grandiflorus samples were collected at 0, 6, 24, and 96 h after MeJA treatment. In the low-temperature treatment, the plants were stored in an incubator at 4 °C, and the P. grandiflorus samples were collected at 0, 1, 3, 6, 9, and 12 h after treatment. Each sample in this study had three independent biological duplicates, and the collected samples were immediately stored at −80 °C and used in the experiment within 12 h.

Identification of bZIPs in P. grandiflorus

The genome database, the protein database, and related annotation files for P. grandiflorus were downloaded from the National Genomics Data Center (https://ngdc.cncb.ac.cn), ID: PRJCA003843 (Jia et al., 2022), and the hidden Markov model (HMM) files for bZIP transcription factor conserved domains (PF07716, PF00170, PF03131, and PF16326) were downloaded from the Pfam database (http://pfam.xfam.org/) (Li et al., 2021). AtbZIP family transcription factors were downloaded from the Plant Transcription Factor Database (http://planttfdb.gao-lab.org/). In the BioEdit software, the P. grandiflorus protein database was used to blast with AtbZIP transcription factors and HMM models, and the E-value was set to 0.000001. The obtained sequences were merged, and the redundant sequences were removed using CD-HIT online analysis (http://cd-hit.org/) (Li & Godzik, 2006). Candidate PgbZIP proteins were subjected to multiple sequence alignments in the MEGAX software to find conserved residues specific to the bZIP family, and proteins that did not possess conserved residues were deleted. Then, the physical and chemical properties of these protein sequences were evaluated using the ExPASy5 ProtParam tool (Savojardo et al., 2018), and the protein subcellular localizations were predicted using the Busca online website (http://busca.biocomp.unibo.it/) (Gasteiger et al., 2003).

Constructing phylogenetic trees of bZIPs in P. grandiflorus

To analyze the evolutionary relationships of conserved protein sequences, we constructed a phylogenetic tree using the MEGAX software for bZIP transcription factor protein sequences obtained from P. grandiflorus and A. thaliana (Kumar et al., 2018). Sequences were aligned using Clustal W method, and evolutionary trees were generated using the neighbor-joining method (Bootstrap values set to 1,000, pairwise deletion of gaps, and the p-distance method were used) in MEGAX.

Chromosomal localization and gene duplication events of PgbZIPs

The inferred position information and chromosome lengths of bZIP genes were obtained from the general feature format file of the P. grandiflorus genome. According to the related literature (Wang et al., 2012), the chromosome location and the collinearity analysis were carried out using MG2C (http://mg2c.iask.in/mg2c_v2.1/) and the TB-tool software. We use the Quick Genome gene dot plot tool to describe gene duplication events and the Simple Ka/Ks Calculator tool for the Ka/Ks analysis.

Identification of the PgbZIP conserved motifs

The conserved PgbZIP motifs were analyzed using the online website MEME (http://meme-suite.org/tools/meme) (Bailey et al., 2015). The predicted value was set to 20, and the conserved PgbZIP motifs were obtained. The parameters are set as follows: motif sites are distributed sequentially, and each sequence has a contributing motif site; the maximum number that needs to be found is 20; the minimum and maximum motif widths are between 10 and 50; the motif value is retained <1 × 10⁻²⁰; the default value was used for other options for further analysis.

Identification of cis-regulatory elements in PgbZIP gene promoters

The bZIP gene regulatory regions at 2 kb upstream were extracted from the P. grandiflorus coding sequence using the TB-tool software (Chen et al., 2020). The promoter sequence’s cis-regulatory elements were detected using the TBtools software, and the same software was used to display these elements.

Real-time quantitative PCR

Using the TransZol kit (TransGen Biotech, Inc., Beijing, China), the total RNA of P. grandiflorus was extracted according to the manufacturer’s instructions. cDNA was synthesized using a FastKing RT Kit (Tiangen, Beijing, China) and reverse transcription reaction. The primers were designed using the Primer 5.0 software (Table S1). Gel electrophoresis showed the RNA integrity, and the RNA concentration was detected using a Nanodrop microspectrophotometer. The RNA purity was good, RIN > 7, and the A260/280 ratio was between 1.9 and 2.1. Optimus Biotechnology Co., Ltd. synthesized the primers, and the 18s rRNA gene was used as the internal reference gene (Su et al., 2021). The RT-qPCR analysis was performed using a SuperReal PreMix Plus SyBr Green PCR kit (Qiagen, Hilden, Germany) on a Cobas z480 Real-Time PCR System using the method described by Su et al. (2021). The thermal program of the PCR is as follows: pre-incubation at 95 °C for 15 min; three steps for amplification over 50 cycles, 95 °C for 10 s, 58 °C for 20 s, and 72 °C for 20 s. Each reaction was repeated three times, and a negative control group was set up at 0 h. The relative gene expression was calculated using the 2^−ΔΔCt method. One-way analysis of variance (ANOVA) was applied for data processing, and statistical differences were compared using t-tests based on IBM SPSS Statistics 23. Finally, visual analysis was performed using GraphPad Prism 8.0.1 (GraphPad Software, La Jolla, CA, USA).

Analysis of PgbZIP transcription factors

According to the similarity to AtbZIPs and the HMM of the bZIP–DNA binding domain, a total of 75 proteins were found in P. grandiflorus. The bZIP basic region contains ~16 amino acid residues. The most striking characteristic of this basic region is the presence of invariant N-x7-R/K conserved residues (Lai et al., 2024). Through multiple sequence alignments of the 75 proteins, we found that only 46 proteins contain the typical N-x7-R/K conserved residues (Fig. S1). These proteins are named PgbZIP001 to PgZIP046 (Table S2). According to the analysis of protein sequence characteristics, the lengths of these 46 PgbZIP proteins varied between 139 and 730 amino acids; the relative molecular weights were between 16.46 and 78.77 kDa; the average relative molecular weight was 36.38 kDa; the isoelectric points were between 4.65 and 10.23; the average isoelectric point was 7.17. All PgbZIPs have negative grand average of hydropathy (GRAVY) values, which indicates that PgbZIP proteins are hydrophilic proteins. The subcellular localization prediction showed that PgbZIPs are located in the nucleus (45) and cytoplasm (one) (Table S3).

Phylogenetic evolution and classification of PgbZIPs

A.thaliana is one of the earliest plants studied for the bZIP transcription factor family, holds significant reference value, and is commonly employed as a model genetic plant for plant genetic research (Droge-Laser et al., 2018; Jakoby et al., 2002). To analyze the phylogenetic relationship of these PgbZIPs, we constructed a phylogenetic tree using 127 AtbZIP proteins and 46 PgbZIP proteins (Fig. 1). According to the reports of A. thaliana, the PgbZIP phylogenetic tree was divided into nine groups (A-E, G-I, and S) (Droge-Laser et al., 2018). In the PgbZIP family, the largest groups are Groups S and A (each with ten members), and the smallest group is Group H (one member).

Chromosomal location and collinearity analysis of PgbZIPs

These PgbZIP genes were physically located on the chromosome in P. grandiflorus (Pgchr) using the MG2C online tool (Fig. 2A). The PgbZIP lengths and positions on the chromosomes are shown in Table S3. A total of nine chromosomes were shown in P. grandiflorus genomic data. The results showed that these forty-six PgbZIPs were distributed on eight Pgchrs. Among them, most PgbZIPs were found on Pgchr1, Pgchr2, and Pgchr3 (each with eight members), and the smallest number of PgbZIPs (two members) were found on Pgchr9. This PgbZIP gene distribution pattern on chromosomes may be related to differences in chromosome structure and size.

Gene duplication events are important in promoting plant evolution, and new gene functions are derived through this process. In this study, we analyzed gene duplication events during P. grandiflorus evolution (Figs. 2B and 2C). Gene duplication event analysis uncovered 1,476 tandem duplication events in P. grandiflorus, but among them, only one tandem duplication event was identified on Pgchr2 among the PgbZIPs studied (Fig. 2C). Collinearity analysis showed 1,913 suspected gene duplication events during P. grandiflorus evolution. After further analysis, 1,689 duplication events were found to be WGDs or segmental duplications. The dot map shows that most WGD events occurred on Pgchr1, Pgchr2, Pgchr5, Pgchr6, and Pgchr7 (Fig. 2B). In addition, 11 gene segmental duplication events were identified between PgbZIPs (Fig. 2C). The frequency of PgbZIP gene segmental duplication events was the highest on Pgchr1, with seven gene segmental duplication events, followed by Pgchr2 with five gene segmental duplication events. As shown in the Fig. 2C, the PgbZIP genes between Pgchr1 and Pgchr2 segmental were duplicated most frequently and have the strongest homology. For a gene that is not under pressure from natural selection, the value of Ka/Ks is close to 1. The Ka/Ks analysis showed that the Ka/Ks of segmental duplication gene pairs was much less than 1 (Table 1). Therefore, we can judge that they have a strong purifying selection.

Conserved motifs, conserved domains, and gene structures of PgbZIPs

The existence of motifs in the same phylogenetic group of PgbZIPs was similar (Figs. 3A and 3B). Twenty motifs were identified in PgbZIPs, ranging from six to 50 amino acids in length (Figs. 3B and 3E). Motif 2, Motif 3, Motif 4, Motif 5, Motif 10, Motif 11, and Motif 16 were found mainly in PgbZIPs from Group D. Motif 6, Motif 15, and Motif 20 were found mainly in PgbZIPs from Groups E and I. Motif 13, Motif 14, Motif 17, and Motif 19 were found mainly in PgbZIPs from Group G. Motif 8, Motif 9, and Motif 12 are found mainly in PgbZIPs from Group A, and Motif 18 is found mainly in PgbZIPs from Group S. Motif 7 is shared by Groups B, H, G, A, C, and S. Moreover, Motif 1 is widely found in all PgbZIP groups. Motif 1 included N-x7-R/K conserved residues so that we can find them in PgbZIP001–PgbZIP045. Motif 1 is not present in PgbZIP046, but in the conserved domain analysis, we found that PgbZIP046 has the multifunctional mosaic region (MFMR) superfamily domain. The MFMR region contains a nuclear localization signal in bZIP opaque, and this region is found in the N-terminus of the bZIP transcription factor domain.

Conserved domain analysis reveals that the existence of domains in the same phylogenetic group of PgbZIPs was similar as well (Fig. 3C). The bZIP_HBP1b-like domain, shared by most PgbZIPs in Group D, contains a C-terminal DOG1 domain that is a specific seed dormancy plant regulator (Mikami, Sakamoto & Iwabuchi, 1994). The bZIP_plant_RF2 domains are owned by the PgbZIPs in both Groups E and I. Transcription factors with bZIP_plant_RF2 domains are similar to Oryza sativa RF2a and RF2b, and transgenic rice with increased RF2a and RF2b display increased resistance to rice tungro disease with no impact on plant development (Dai et al., 2004). The bZIP_plant_BZIP46 domain, present in most PgbZIPs in Group A, is similar to GmBZIP46, which may be a drought-response gene (Liao et al., 2008). The MFMR domain is possessed by the PgbZIPs in Group G, and it has been suggested that some of these domains may be involved in mediating protein–protein interactions (Sibéril, Doireau & Gantet, 2001). The bZIP_plant_GBF1 domain of PgbZIPs in Group S is a DNA-binding and dimerization domain (Norén Lindbäck et al., 2023). The bZIP_plant_GBF1 domain contains an N-terminal proline-rich domain in addition to the bZIP domain. The transcription factors with this domain are involved in developmental and physiological processes in response to stimuli such as light or hormones and can also affect lysine content and carbohydrate metabolism, acting indirectly on the starch/amino acid ratio.

To reveal the PgbZIP gene family’s structural evolution, the PgbZIP exon–intron structures were analyzed (Figs. 3D and 3F). The PgbZIP gene structures displayed significant variability and diversity in terms of the relative positions and numbers of introns and exons. PgbZIPs in Group B contained 1–3 introns, PgbZIPs in Group D contained 7–11 introns, PgbZIPs in Groups A and E contained 2–4 introns, PgbZIPs in both Groups H and I contained three introns, PgbZIPs in Group G contained 7–11 introns, PgbZIPs in Group A contained 0–4 introns, PgbZIPs in Group C contained five introns, and PgbZIPs in Group S contained 0–1 introns. PgbZIPs with close phylogenetic relationships may possess similar gene structures.

The similar domains, motifs, and gene structures in the same PgbZIP subgroups may lead to their similar functions, and the above results may provide some theoretical basis for future PgbZIP functional studies.

Cis-regulatory elements of PgbZIPs

To explore PgbZIP gene responses to various stresses, a 2,000 bp DNA sequence region upstream of PgbZIPs was retrieved as the promoter region and used to analyze the PgbZIPs’ cis-regulatory elements (Fig. 4). The results showed that among the PgbZIP gene promoters, the number of light responsive elements was the largest (518), followed by MeJA responsive elements (126). In addition, the numbers of abscisic acid responsive elements (90), anaerobic induction element (88), and MYB binding sites (57) are also extremely high. Notably, the low-temperature responsive elements (33) were the second most numerous in addition to the light responsive elements for natural environmental stresses.

Expression pattern of PgZIP genes upon exposure to low-temperature and MeJA stresses

There is usually a correlation between gene expression pattern and function, and plant bZIP transcription factors are widely involved in various stresses. To investigate the response patterns of PgbZIPs to MeJA and low temperatures, PgbZIP003 (eight), PgbZIP005 (six), PgbZIP02 (four), PgbZIP021 (two), and PgbZIP011 (zero) with different numbers of MeJA responsive elements and PgbZIP046 (five), PgbZIP019 (three), PgbZIP041 (two), PgbZIP021 (one), and PgbZIP032 (zero) with different numbers of low-temperature responsive elements were selected to study their gene expression patterns upon exposure to low-temperature and MeJA stresses. RT-qPCR was conducted at 0, 6, 24, and 96 h after MeJA treatment and 0, 1, 3, 6, 9, and 12 h after low-temperature treatment (Fig. 5). The expression levels of nine PgbZIPs (PgbZIP003, PgbZIP005, PgbZIP011, PgbZIP019, PgbZIP021, PgbZIP022, PgbZIP024, PgbZIP032, PgbZIP042, and PgbZIP046) were significantly decreased up to 6 h of MeJA treatment and gradually increased again after 6 h. The expression levels of seven PgbZIPs (PgbZIP003, PgbZIP005, PgbZIP019, PgbZIP021, PgbZIP024, PgbZIP042, and PgbZIP046) remained significantly lower than that of untreated (0 h) when treated with MeJA for 24 h. However, when the treatment time reached 96 h, most PgbZIPs (PgbZIP011, PgbZIP019, PgbZIP021, PgbZIP024, PgbZIP042, and PgbZIP046) returned to the untreated expression level or higher than that of the untreated expression level. In addition, PgbZIP032 expression increased gradually after MeJA treatment. The expression levels of most PgbZIPs may decrease significantly (up to four-fold) within 6 h and then gradually increase to normal or above normal levels within 90 h after MeJA treatment. This result indicated that MeJA may have an inhibitory effect on the PgbZIP gene in the early response stage but increased PgbZIP gene expression after some time.

After low-temperature treatment, the expression trends of PgbZIPs was generally decreasing from 0 to 12 h. PgbZIP003, PgbZIP005, PgbZIP011, and PgbZIP032 were all significantly lower than untreated at 12 h after the low-temperature treatment. At the same time, the expression levels of the other six PgbZIPs were also significantly lower than untreated at most time points after treatment. Notably, the expression levels of all PgbZIPs were significantly reduced after 3 h of the low-temperature treatment, with the largest decrease of about 35-fold for PgbZIP032. This shows that PgbZIP gene expression levels may be inhibited by low-temperature stress, and the occasional increase in expression level may be caused by P. grandiflorus adaptation to the low-temperature environment.

PgbZIPs response to low-temperature and MeJA stresses

The bZIP genes are involved in a variety of abiotic stresses in plants, including low-temperature, drought, and high-temperature stresses, and biotic stresses, such as diseases and pathogens, as well as in a variety of hormone-induced processes (Amorim et al., 2017; Banerjee & Roychoudhury, 2017; Lakra et al., 2015; Liu, Wu & Wang, 2012; Zhang et al., 2015, 2018). The RT-qPCR results demonstrated that PgbZIPs showed substantial changes in gene expression levels upon exposure to low-temperature and MeJA stresses. PgbZIP gene expression may be significantly inhibited by low-temperature stress. Transgenic A. thaliana plants heterologously expressing Triticum aestivum bZIP6 (TabZIP6) and Camellia sinensis bZIP6 (CsbZIP6) upon exposure to low-temperature treatments showed reduced survival, increased relative conductivity, increased malondialdehyde content, and reduced soluble sugar content (Cai et al., 2018). This indicated that the overexpression of TabZIP6 or CsbZIP6 reduced the freezing tolerance of transgenic A. thaliana plants. Rice plants overexpressing OsABF2 showed enhanced cold, salt, and drought tolerances (Banerjee & Roychoudhury, 2017). However, in the rice cold tolerance test, it was found that the OsbZIP52 transgenic plant survival rate was 18%, and the survival rate of the wild type was 78.03%, which shows that OsbZIP52 expression was negatively correlated with the low-temperature tolerance in rice (Liu, Wu & Wang, 2012). This is consistent with our findings that high PgbZIP expression levels may lead to a decrease in cold tolerance in P. grandiflorus plants so that a decrease in PgbZIP expression levels was observed after the low-temperature treatment. The main active constituents in P. grandiflorus are platycodins, and platycodin synthesis may be triggered by MeJA (Huang et al., 2021; Kim et al., 2004; Scholz et al., 2009). A significant increase in the final expression level of most PgbZIPs occurred in P. grandiflorus after MeJA treatment. This suggests that PgbZIPs respond to MeJA induction, but whether PgbZIPs are involved in saponin content regulation in P. grandiflorus requires further study.