The plant material P. grandiflorus was cultivated from the Medicinal Botanical Garden of Heilongjiang University of Chinese Medicine (HLJUCM), Harbin, China. There are two types of plant materials used in this study, one of which was the roots, stems, leaves and flowers of 1-year-old P. grandiflorus in the normal growth state without treatment. The other was a three-month-old P. grandiflorus seedling grown in a simulated abiotic stress environment. P. grandiflorus seedlings were watered with 200 mmol/L NaCl and 20% PEG6000 (Coolaber Science & Technology, Beijing, China) to simulate drought and salt stress, respectively, and roots were collected at 48 h. Following the sample process, all plant components were instantly frozen in liquid nitrogen immediately after sampling and kept in a medical refrigerator at -80°C for additional analysis.

The complete genome sequence and annotated file of P. grandiflorus can be downloaded from the Figshare database (https://doi.org/10.6084/m9.figshare.21511020). TBtools software was used to identify every member of the bZIP TFs in the genome of P. grandiflorus (Chen et al., 2020). TBtools software was used to extract all genes sequences in the P. grandiflorus genome. The hidden Markov model file in the bZIP domain was used as a template to identify each gene and obtain the candidate PgbZIP protein sequence. For validation, the candidate sequences were uploaded to the NCBI CD-search tools and the plant transcription factor database. After removing redundant and structurally incomplete domain sequences, the final PgbZIP protein was obtained. Use the web tool ExPASy ProtParam (https://web.expasy.org/protparam/) to search the molecular weight (MW), isoelectric point, grand average of hydropathicity (GRAVY) of PgbZIP proteins.

The genome annotation file of P. grandiflorus was analyzed and the chromosome location information of PgbZIP gene was obtained, and its corresponding chromosomal physical location was mapped using Tbtools software. The genomic information of Malus domestica (GCF_002114115.1), Cannabis sativa (GCA_900626175.1), Vitis vinifera (GCF_000003745.3), Solanum lycopersicum (GCF_000188115.4) and Oryza_sativa (GCA_001433935.1) were downloaded from the NCBI database. MCScanX software was used to analyze the gene duplication and collinearity of PgbZIP gene family (Wang et al., 2012). The PgbZIP gene replication events were visualized using the TBtools software’s “Advanced Circos” function, and the collinearity relationship between P. grandiflorus and other species was shown using the “Dual Systeny Plot” function.

The TAIR database (https://www.arabidopsis.org/) provided the sequence of the A. thaliana bZIP protein. MEGA X was used to perform multiple sequence alignment of AtbZIP, IibZIP (Jiang et al., 2022) and PgbZIP protein, respectively (Kumar et al., 2018). The phylogenetic tree was constructed by neighbor- joining (NJ) method, and the bootstrap was repeated 1000 times. PgbZIP proteins are divided into different subfamilies based on A. thaliana bZIP proteins (Jakoby et al., 2002).

The PgbZIP gene family was subjected to motif analysis using the Multiple Em for Motif Elicitation website. The parameters were set to Classic mode, and ten motifs were discovered (Bailey et al., 2009). Extract the upstream 2 kb sequence of the PgbZIP family gene using the “GXF Sequences Extract” tool in TBtools software, and then identify and retrieve the cis-elements of the gene family through the PlantCare database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for analysis and comparison (Lescot et al., 2002).

The NCBI Sequence Read Archive (SRA) database has the raw transcriptome data for the eight P. grandiflorus tissues that were used in this investigation. These accession numbers are as follows: SER912510-SER912517 (Kim et al., 2020). The expression of every PgbZIP gene used to be quantified in fragments per kilobase of exon model per million mapped fragments (FPKM). All FPKM values have been processed using row scale transformation, and heatmaps were generated using TBtools.

Total RNA was extracted from leaves, roots and stems of 1-year-old P. grandiflorus, drought, and salt stress treated P. grandiflorus seedlings, respectively, using the RNAprep Pure Plant Kit (TIANGEN Biotech Co.,Ltd, Beijing, China). Then, The first strand cDNA was synthesized using MS 1st Strand cDNA Synthesis SuperMix for qPCR (+gDNA wiper) kit (Msunflowers Biotech Co.,Ltd, Beijing, China) reverse transcribed RNA as a template, and finally, qRT-PCR experiments were performed using 2×SYBR Green qPCR Mix (With ROX) kit (Shandong SparkJade Biotechnology Co., Ltd, Shandong, China) kit. qRT-PCR experiments had been carried out the use of an AriaMx real-time PCR system (Agilent Technologies, USA). The PgGAPDH had been used as a reference gene (Ma et al., 2016), and three technical replicates had been set up for each experiment. The relative expression of genes was computed by the 2^{− ΔΔCT} Method (Livak and Schmittgen, 2001). The Primer3web (version 4.1.0) website (https://primer3.ut.ee/) was used to layout unique primers (Untergasser et al., 2012). GraphPad Prism (v8.0.2) software program was once used to draw the histogram of relative expression of genes, which was once analyzed by t-test.

PlantTFDB was used to identify TFs for the whole genome proteins of P. grandiflorus, and the online database of Pfam was used to annotate the conserved domains of genes. The correlation between transcriptomes was calculated using a Python script primarily based on the Pearson correlation coefficient (PCC) (Zhang et al., 2017). Visualization of correlation network diagrams using Cytoscape (v3.7.1) software (Shannon et al., 2003). The GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis using the R language package clusterProfiler (Zhang et al., 2017).

After removing redundant and incomplete domain sequences, we finally identified 47 bZIP genes in P.grandiflorus ( Supplementary Table 1 ). In subsequent analyses, we named the genes PgbZIP1 to PgbZIP47 based on their location on the chromosome or contig (from Pg_chr01_03120T to Pg_contig01856_00120T). Then, we analyzed the physicochemical parameters of the PgbZIP protein in detail. The results confirmed that the amino acid (AA) size of the PgbZIP proteins ranged from 132 AA (PgbZIP40) to 704 AA (PgbZIP13). The molecular weight (MW) of the PgbZIP proteins ranged from 14840.52 Da to 75818.79 Da, with an theoretical isoelectric point (pI) of 4.66 (PgbZIP28) to 9.77 (PgbZIP27). The GRAVY values of all PgbZIP proteins were less than 0, suggesting that all PgbZIP proteins may be hydrophilic.

The distribution of PgbZIP genes on the chromosomes of P.grandiflorus did not show any obvious pattern, with 47 PgbZIP genes unevenly distributed on 9 chromosomes and 2 contigs ( Figure 1 ). Chromosome 3 contained the most bZIP genes (9), and chromosome 9, Pg_contig00116 and Pg_contig01856 contained the fewest bZIP genes (only a single gene each). The bZIP genes were distributed on all chromosomes, but none of the genes showed a preferential distribution on a particular chromosome.

Investigations into the gene replication events in P. grandiflorus were conducted to understand the expansion of the PgbZIP genes. Gene duplication within the same chromosome or on different chromosomes, but without one following the other, is considered a segmental duplication event. In the PgbZIP gene family, four PgbZIP gene pairs were generated by segmental duplication. It is worth noting that 1 pair of segmental duplication events (PgbZIP5 and PgbZIP46) is not shown in the Figure 1B , because the PgbZIP46 gene was not mapped to the chromosome and was distributed on the contig.

To explore the evolutionary relationships and classification of the PgbZIP gene family, we constructed NJ phylogenetic trees using MEGA X software ( Supplementary Table 2 ). Based on the earlier classification results of the A. thaliana and I. indigotica bZIP gene family, 47 PgbZIP proteins were divided into 11 subfamilies ( Figure 2 ), comprising A (10), B (1), C (2), D (6), E (4), F (1), G (5), H (2), I (6), K (1) and S (9). The outcome of phylogenetic tree evaluation confirmed that the bZIP genes of P. grandiflorus, A. thaliana and I. indigotica were highly homologous in each cluster.

To explore the structural similarity of PgbZIP proteins, 10 conserved motifs have been recognized the use of the MEME website ( Figure 3A ). These motifs vary in size from 21–50 AA, with at least 1 to 5 motifs distributed across 47 PgbZIP proteins, and all sequences have motif 1. Members of gene families with close evolutionary relationships have similar motif compositions; for example, G subfamily members all have motifs 1 and 3. Furthermore, all PgbZIP proteins contain conserved domains of the bZIP genes, proving that the outcome of gene identification had been dependable ( Figure 3B ).

To research the potential transcriptional regulatory role of PgbZIP genes, in this study, cis-elements were extracted from the PlantCARE database using 2K bp upstream of the PgbZIP genes as the promoter region. A total of 929 valuable cis-elements were recognized in the promoter regions of 47 PgbZIP genes, these cis-elements could be broadly categorized into three groups ( Figure 3C ) ( Supplementary Table 3 ). Relevant to plant growth include light responsive elements, meristem expression, endosperm expression, anaerobically induced regulatory elements, seed-specific regulatory response elements, cell cycle regulatory elements, and circadian regulatory elements. Hormone response related elements include MeJA responsive elements, auxin responsive elements, abscisic acid responsive elements, and gibberellin responsive elements. Involved in abiotic stress related including stress responsive and defense elements and low-temperature responsive elements. Interestingly, the PgbZIP gene has cis-elements that bind to the MYB gene, participating in photoresponse, drought response, and flavonoid synthesis. This suggests that the PgMYB gene may regulate PgbZIP, forming regulatory networks that exercise different biological functions.

In order to further clarify the evolution of the PgbZIP genes, we compared the PgbZIP gene sequences and analyzed the coding regions and introns ( Figure 3D ). Except for seven members of the S family who do not have introns, the number of introns in other members is distributed between 1 and 12. As anticipated, members of the same subfamily have a relatively conserved number of introns, such as the H subfamily, which has 2 introns, and the E subfamily, which has 3 to 4 introns.

For investigating further the relationship of bZIP genes evolution among various plants, we compared interspecific synteny of P. grandiflorus with that of five dicotyledonous plants and one monocotyledonous plant (O. sativa) ( Figure 4 ). Dicotyledonous plants include three fruit (M. domestica, V. vinifera and S. lycopersicum), and two medicinal plants (C. sativa and C. lanceolata). The results showed that the collinearity relationship between PgbZIP genes and MdbZIP genes was the closest, with 79 pairs of genes, 40 pairs of PgbZIP and VvbZIP genes, 32 pairs of PgbZIP and SlbZIP genes, 23 pairs of PgbZIP and ClbZIP genes, and the same pairs of CsbZIP genes. PgbZIP and OsbZIP had the worst collinearity relationship, with only 7 pairs of genes. These results indicate that most of these homologous bZIP genes occur after the differentiation of dicotyledons and monocotyledons. It is worth noting that PgbZIP1 (Pg_chr01_03120T) and PgbZIP19 (Pg_chr03_29870T) have syntenic gene pairs with 6 other species, which may have existed before species divergence and participated in the evolution of these plants.

To investigate the expression profile of bZIP genes in P. grandiflorus tissues, the gene relative expression of 47 PgbZIP genes were analyzed based on RNA-seq data of P. grandiflorus root, leaf, seed, petal, stem, stamen, pistil, and sepal ( Figure 5 ) ( Supplementary Table 5 ). The findings indicated that 20 genes were expressed in leaf, 28 genes in petal, 34 genes in pistil, 32 genes in root, 38 genes in seed, 37 genes in sepal, 37 genes in stamen and 37 genes in stem (FPKM>0.5). A total of 10 genes (PgbZIP5, PgbZIP13, PgbZIP14, PgbZIP21, PgbZIP25, PgbZIP28, PgbZIP30, PgbZIP33, PgbZIP42 and PgbZIP45) were highly expressed in all tissues (FPKM>10), and these genes probably participating in the whole development processes of P. grandiflorus. PgbZIP16 and PgbZIP17 were not expressed in all tissues and may be pseudogenes or require specific conditions to activate expression. Interestingly, PgbZIP genes are also tissue-specific. For example, PgbZIP26 and PgbZIP34 are expressed only in stamens. PgbZIP12, PgbZIP35, PgbZIP37 and PgbZIP43 are only expressed in seeds. These genes involved in specific tissue expression may only be involved in the biological process of this tissue.

As mentioned above, subfamilies A (PgbZIP14), B (PgbZIP13), C (PgbZIP30), F (PgbZIP42), G (PgbZIP21), I (PgbZIP25/45), K (PgbZIP28) and S (PgbZIP5/33) have been expressed at very high levels in various tissue parts of P. grandiflorus, and are probably participated in the regulation of developmental process in all tissues. Therefore, those 10 PgbZIP genes were selected as candidate genes for qRT-PCR experiments to validate in this research ( Supplementary Table 6 ). The relative expression of the 10 candidate genes were basically consistent to the expression trends obtained from the RNA-seq data ( Figure 6 ). Notably, PgbZIP28 showed higher expression in petals, and PgbZIP30 and PgbZIP33 were expressed at higher levels in sepals than stems. All of the above suggest that these 10 PgbZIP genes possibly are closely involved to the developmental process of P. grandiflorus.

The root of P. grandiflorus is the main medicinal part. For studying the potential role of PgbZIP gene family under abiotic stress, the qRT-PCR experiments analysis was conducted using the roots of P. grandiflorus seedlings under drought and salt stress as templates ( Figure 7 ). Compared with the normally growth group (CK), the expression levels of 4 genes (PgbZIP5, PgbZIP21, PgbZIP25 and PgbZIP28) were increased under both drought and salt stress. PgbZIP33 was only highly expressed under salt stress, but its expression was reduced under drought stress. These PgbZIP genes with increased expression levels under drought and salt stress could help P. grandifloru to resist abiotic stress.

The bZIP genes often form an interaction networks with various TFs to participate in plant developmental processes. In order to mine the transcriptional regulatory network of PgbZIP in different tissues, this study characterized TFs in the RNA-seq data of 8 different tissues of P. grandiflorus and constructed a correlation network. A total of 1567 TFs were identified and categorized into 58 gene families, among which the top 5 were ERF (136), bHLH (118), C2H2 (115), MYB (105) and NAC (90). The expression profiles of all genes in the RNA-seq data of 8 different tissues of P.grandiflorus were analyzed using Python script to explore the TFs co-expressed with the PgbZIP gene. In the correlation network, a total of 19 PgbZIPs were co-expressed with 149 TFs, and 168 nodes with 361 network pairs were found. The positive (r > 0.7) (P < 0.05) and negative (r < -0.7) (P < 0.05) correlation network pairs are 232 and 129, respectively ( Figure 8 ). PgbZIP45, PgbZIP30, PgbZIP13 and PgbZIP32 were presented as hub genes (degree≧30). The highest correlations (degree≧7) among TFs were Pg_chr05_09390T (GATA), Pg_chr08_31480T (C3H), Pg_chr02_03150T (MYB-related), Pg_chr01_14460T (AP2), Pg_chr06_07740T (bHLH), and Pg_contig26001_00320T (NF-YA) ( Supplementary Table 7 ). These TFs has been shown to function in regulating plant growth and development in a large number of studies. Therefore, PgbZIP genes could be interacting with these TFs to form a regulatory network and participate in the development process of P.grandiflorus.

In order to further clarify the functions of TFs in the network in biological processes and the pathways involved, we performed GO annotation, GO and KEGG enrichment analyses of all TFs in the correlation network ( Figure 9 ). The results of the analysis of the GO annotation findings suggests that those genes function in the molecular function classification for recruitment of TFs and activation of transcriptional activity ( Figure 9A ). In the cellular component classification, these genes were found to be predominantly distributed in the telomeric region, chromosome, nuclear chromosome and spindle microtubule. In biological processes, these genes are mainly involved in response to ethylene, negative regulation of macromolecular biosynthesis processes, negative regulation of RNA metabolic processes and ethylene response ( Supplementary Table 8 ).

Similarly, we performed GO and KEGG enrichment analysis of TFs in the correlation network. The outcome of GO enrichment analysis indicated that the main functions of these genes are transcription factor binding transcription factor activity, cellular response to ethylene stimulus, meristem maintenance, response to ethylene, xylem development, phloem development, floral organ morphogenesis, meristem maintenance, and so on ( Figure 9B ). The outcome of KEGG enrichment indicated that these TFs were primarily enriched in circadian rhythm, phytohormone signal transduction, nucleotide excision repair, plant-pathogen interaction, pyrimidine metabolism, purine metabolism and so on ( Figure 9C ). These outcomes indicate that the TFs in the correlation network are participating in most bioprocesses in plants and play a certain role.