Identiﬁcation and Expression Analysis of bZIP Members under Abiotic Stress in Mung Bean ( Vigna radiata )

: The main aim of this study was to identify the bZIP family members in mung bean and explore their expression patterns under several abiotic stresses, with the overarching goal of eluci-dating their biological functions. Results identiﬁed 75 bZIP members in mung bean, which were unevenly distributed in the chromosomes (1–11), and all had a highly conserved bZIP domain. Phylogenetic analysis divided the members into 10 subgroups, with members in the same subgroup having similar structure and motif. The cis -acting elements in the promoter region revealed that most of the bZIP members might have the connection with abscisic acid, ethylene, and stress responsive elements. The transcriptome data demonstrated that bZIP members could respond to salt stress at different degrees in leaves, but the expression patterns could vary at different time points under stress. Differentially expressed genes (DEGs), such as VrbZIP12 , VrbZIP37 , and VrZIP45 , were annotated into the plant hormone signal transduction pathway, which might be regulated the expression of abiotic stress-related gene ( ABF ). Quantitative real-time polymerase chain reaction (qRT-PCR) was applied to determine the expression of bZIP members in roots and leaves under drought, alkali, and low-temperature stress. Results showed that bZIP members respond differently to diverse stresses, and their expression was tissue-speciﬁc, which suggests that they may have different regulatory mechanism in different tissues. Overall, this study will provide a reference for further research on the functions of bZIP members in mung bean.


Introduction
Mung bean (Vigna radiate), a leguminous crop, has been cultivated in China for more than two thousand years. Given its value in medicine, health care, and nutrition, the demand for mung bean has been rapidly increasing in recent years. However, mung bean production has been maintained at a relatively low level in some major producing areas in north China largely due to the limitation of drought, soil salinization, low temperature in spring, and other adverse environmental conditions. It is well-known that transcription factors (TFs) play important roles in almost all biological processes of plants and are key regulatory factors responding to various signals in plant growth and development as well as environmental stress. Mounting evidence has revealed that transcription factors activate or inhibit transcriptional expression of downstream target genes by binding to promoter or enhancer regions of regulatory genes [1][2][3]. Therefore, this calls for studies to explore the biological characteristics and stress response of transcription factors for breeding and variety improvement of mung bean. The basic leucine zippers (bZIP) family is one of the largest and most diverse TF families that only exists in eukaryotes [4,5]. All bZIP family members contain a highly conserved bZIP domain consisting of 60-80 amino acids, which is composed of two parts: a highly conserved alkaline region and a variable leucine zipper. The alkaline region, consisting of 16 amino acid residues with a constant N-X7-R/K motif, is responsible for specific recognition and binding to the DNA sequence on the promoter. The leucine zipper, consisting of seven-peptide repeats of leucine or other hydrophobic amino acids, plays an oligomerization function through which the bZIP transcription factor specifically recognizes and forms homologous or heterologous dimers to complete transcriptional activation or inhibition of downstream genes [6,7]. Numerous studies have shown that bZIP transcription factors play an important role in various aspects of plant biological processes, such as embryogenesis [8], seed maturation [9,10], and flower and vascular development [11,12]. On the other hand, bZIP protein is also involved in signal regulation and response to various biotic and abiotic stresses, including high salinity, drought, osmosis, cold stress, and pathogen infection [13]. Studies have shown that the OsbZIP73 gene can improve cold tolerance in rice [14], AtbZIP17 and AtbZIP28 genes can regulate root elongation during stress response [15], and the GmbZIP2 gene characterized from soybean can improve drought and salt tolerance of transgenic tobacco [16].
Over the years, members of the bZIP transcription factor family have been identified or predicted in many eukaryotic genomes [17][18][19][20]. However, no study has identified and characterized the bZIP genes in mung bean. Herein, given the completion of mung bean genome sequencing, we used mung bean genome and transcriptome data to analyze the bZIP family genes from the aspects of systematic evolution, gene structure, protein motif, and chromosome location and then analyzed the expression pattern of the genes under several stresses, including salt, alkali, drought, and low temperature. This study provides a theoretical basis for further exploring the function of bZIP genes in mung bean resistance to various stresses.

Materials and Treatment
Five mung bean seeds (Jilv 10 variety) were sown in a pot (diameter = 8 cm) filled with nutrient soil and cultured in a greenhouse at 22-24 • C and 14 h light and 10 h darkness. After the emergence of mung bean, the two needles were fully expanded, and then 200 mM NaCl, 200 mM mannitol, and 75 mM Na 2 CO 3 were, respectively, poured onto the root of seedlings. The other seedlings were put in an artificial climate chamber at 4 • C for cold stress, and other conditions were the same as in the greenhouse. Root and leaf samples were collected at 0 h, 6 h, and 24 h time points under 200 mM mannitol, 75 mM Na 2 CO 3 , and 4 • C stress and immediately immersed in liquid nitrogen, followed by storing in the refrigerator at −80 • C for RNA extraction and gene expression analysis. For transcriptome sequencing, leaf samples under 200 mM NaCl stress were collected at 0 h, 6 h, 12 h, and 24 h time points. Notably, all experiments were replicated three times.

Identification and Bioinformatics Analysis of bZIP Family Genes in Mung Bean
We first searched the bZIP genes in mung bean transcriptome data to obtain the gene ID and then searched in NCBI to obtain the gene and protein sequences. Next, the conserved domains of bZIP protein were evaluated using conservative structure domain tools in NCBI database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi, accessed on 1 October 2021). Physical and chemical properties of bZIP proteins were analyzed by ExPASy (http://web.expasy.org/protparam/, accessed on 3 October 2021), whereas the motif of bZIP protein was analyzed by MEME (http://meme-suite.org/tools/meme, accessed on 4 October 2021). TBtools was used to visualize the results of protein phylogeny, gene structure, and protein motif analysis [21]. Finally, 2000 bp promoter sequences upstream of the CDS were extracted to predict the cis-element using Plant CARE

Phylogeny of bZIP Proteins in Mung Bean and Arabidopsis Thaliana
The bZIP protein sequences of Arabidopsis thaliana were downloaded from the TAIR database (https://www.arabidopsis.org/, accessed on 12 September 2021). Phylogenetic analysis of bZIP proteins from mung bean and Arabidopsis thaliana was carried out by MEGA7.0 adjacency method, with the bootstrap value set to 1000 and other parameters set as default.

Identification and Bioinformatics Analysis of bZIP Family Genes in Mung BeanRNA Extraction and Quantitative Real-Time PCR (qRT-PCR) Analysis
Total RNA was extracted from mung bean tissues using the plant RNA extraction kit (TIANGEN), followed by reverse transcription into cDNA using TaKaRa's PrimeScriptTM RT Reagent Kit according to the manufacturers' instructions. Primers for bZIP and reference genes were designed by NCBI primer design tools (https://www.ncbi.nlm.nih.gov/tools/ primer-blast, accessed on 10 November 2021) ( Table 1). SYBR Green Realtime PCR Master Mix kit was used for quantitative analysis, and each sample was replicated three times. Reaction system 20 µL: 2 µL cDNA, upstream and downstream primers (10 µmol/L) 2 µL, 10 µL SYBR, 6 µL ddH 2 O. Reaction procedure: pre-denaturation at 95 • C for 30 s, denaturation at 95 • C for 10 s, annealing at 60 • C for 30 s, extension at 72 • C for 30 s, cycle 40 times. The genes expression were calculated by 2 -∆∆Ct method.

Expression Analysis of bZIP Genes in Mung Bean Leaves under Salt Stress
According to the transcriptome sequencing results, the expression patterns of bZIP family genes in mung bean leaves under salt stress for 0, 6, 12, and 24 h were analyzed using the FPKM value. Heat maps were then generated using R3.4.3 software. The differential gene screening criteria were FDR (false-discovery rate) < 0.01, FC (fold change) > 2.

Identification and Physicochemical Property Analysis of bZIP Family Members in Mung Bean
bZIP gene sequences were first obtained from the NCBI database and then named according to their position on the chromosome, ranging from VrbZIP1 to VrbZIP75. It was evident that the bZIP members were unevenly distributed on the chromosomes (1-11). Specifically, chromosome 5 contained the most bZIP members with ten; followed by chromosome 8 with nine bZIP members; chromosomes 2, 7, and 11 with six bZIP members each; chromosomes 4 and 6 each with five bZIP members; chromosome 3 with four bZIP members; chromosome 1 with three bZIP members; chromosome 10 with two bZIP members; and finally chromosome 9 with one bZIP gene. There were also 18 bZIP members that have not yet been mapped on the chromosomes (Figure 1). Specifically, chromosome 5 contained the most bZIP members with ten; followed by chromosome 8 with nine bZIP members; chromosomes 2, 7, and 11 with six bZIP members each; chromosomes 4 and 6 each with five bZIP members; chromosome 3 with four bZIP members; chromosome 1 with three bZIP members; chromosome 10 with two bZIP members; and finally chromosome 9 with one bZIP gene. There were also 18 bZIP members that have not yet been mapped on the chromosomes (Figure 1). The CDS for distribution of bZIP members ranges from 369 to 2337 bp and encodes 122~778 amino acids. The amino acid length of bZIP members in group D, E, and S is relatively concentrated. The length of amino acids in group S is shorter, distributed between 134 and 200, whereas the amino acid in group D is longer, distributed between 352 and 552. The molecular weight ranges between 13,714.35 and 84,495.52 Da, and the isoelectric point ranges between 2.42 and 10.23. In addition, the basic protein is mainly concentrated in group A, group S, and group H. Results obtained after conducting conserved domain analysis showed that all the 75 mung bean bZIP proteins had bZIP characteristic domains; however, the domain of VrbZIP37 was partially missing. The basic structure of the conserved domain of mung bean bZIP protein is N-x9-R/K-x7-L-x6-L-x6-L and is composed of basic domain and leucine zipper domain. The basic domain is made up of 20 amino acids, contains conserved structure N-x9-R/K, and is rich in basic amino acids, such as arginine and lysine residues. The N-terminal of the leucine zipper domain is linked to the basic domain, which is composed of three consecutive seven peptides (x6-L-x6-L-x6-L). There is a conserved basic amino acid (leucine) site at the seventh position of each peptide. It is worth noting that most of the conserved domains of mung bean bZIP proteins are 40~59 amino acids; the longest is VrbZIP32 with 72 amino acids, and VrbZIP20 is the shortest with only 35 amino acids ( Table 2).  The CDS for distribution of bZIP members ranges from 369 to 2337 bp and encodes 122~778 amino acids. The amino acid length of bZIP members in group D, E, and S is relatively concentrated. The length of amino acids in group S is shorter, distributed between 134 and 200, whereas the amino acid in group D is longer, distributed between 352 and 552. The molecular weight ranges between 13,714.35 and 84,495.52 Da, and the isoelectric point ranges between 2.42 and 10.23. In addition, the basic protein is mainly concentrated in group A, group S, and group H. Results obtained after conducting conserved domain analysis showed that all the 75 mung bean bZIP proteins had bZIP characteristic domains; however, the domain of VrbZIP37 was partially missing. The basic structure of the conserved domain of mung bean bZIP protein is N-x9-R/K-x7-L-x6-L-x6-L and is composed of basic domain and leucine zipper domain. The basic domain is made up of 20 amino acids, contains conserved structure N-x9-R/K, and is rich in basic amino acids, such as arginine and lysine residues. The N-terminal of the leucine zipper domain is linked to the basic domain, which is composed of three consecutive seven peptides (x6-L-x6-L-x6-L). There is a conserved basic amino acid (leucine) site at the seventh position of each peptide. It is worth noting that most of the conserved domains of mung bean bZIP proteins are 40~59 amino acids; the longest is VrbZIP32 with 72 amino acids, and VrbZIP20 is the shortest with only 35 amino acids ( Table 2).

Phylogenetic Analysis of bZIP Family Proteins in Mung Bean
Multiple sequence alignment and phylogenetic analysis were performed between the bZIP protein sequences of mung bean and A. thaliana (Figure 2).
According to the classification results of A. thaliana, mung bean bZIP family proteins were divided into 10 groups: A~I and S. Group A had the most members, with a total of Life 2022, 12, 938 6 of 14 17 accounting for 22.7%; followed by group D and group S, both with 14 bZIP proteins accounting for 18.7%; and group I, with nine members accounting for 12%. Group B had the least members, with only one bZIP protein accounting for 1.3%. The largest group in A. thaliana was group S accounting for 22.7%; followed by group A and group I, both accounting for 17.3%; and group D accounting for 13.3%. It was evident that the groups with the largest number of bZIP members were A, D, S, and I, with these four groups containing 72.1% and 70.6% of the total bZIP members in mung bean and A. Thaliana, respectively.

Phylogenetic Analysis of bZIP Family Proteins in Mung Bean
Multiple sequence alignment and phylogenetic analysis were performed between the bZIP protein sequences of mung bean and A. thaliana (Figure 2).

Analysis of Conserved Motif and Gene Structure of bZIP Proteins in Mung Bean
MEME software was applied to determine the conserved motifs of mung bean bZIP family proteins, with obtained results indicating 20 conserved motifs. The conserved motif, gene structure, and phylogenetic relationship of bZIP proteins were then visually analyzed by TBtools software (Figure 3). Figure 3A,B show that the motifs contained in the same group of bZIP proteins were similar in type and arrangement order. The members of group S contained only four kinds of motifs, namely motifs 1, 5, 15, and 20. Group C members contained four kinds of motifs, namely motifs 1, 5, 13, and 14, and motif 14 was unique to this group. Group G contained motifs 1, 5, 13, 17, and 20, with motif 13 and motif 17 being specific to this group, where five members contained both of them. In addition, all members of group G contained motif 20, whereas some of the members of group S had motif 20. Group A members had a total of six motifs, namely motifs 1, 5, 8, 9, 10, and 11, with motif 8, 9, and 10 being specific to this group. Group E had three motifs, namely motifs 1, 15, and 18, with motif 15 and motif 18 only existing in this group. With exception of VrbZIP46, all other members of group I contained motifs 1, 5, 7, and 12, with motifs 7 and 12 existing in all members of group I and remaining unique to this group. The two members in group F all contained motifs 1, 5, and 19, with motif 19 being unique to this group. Group B, which only had one member (VrbZIP70), contained three motifs, namely motifs 1, 5, and 13, with motif  13 being specific to this group. The members of group H only contained motifs 1 and 5 and had no specific motif. Group D had the largest number of motifs, which were motifs 1, 2, 3, 4, 5, 6, 11, and 16. Notably, motifs 2, 3, 4, 6, and 16 were specific to group D. Based on these results, it was evident that motif 1 was common to all members of the mung bean bZIP genes, and motif 5 existed in the vast majority of members. This is because motif 1 is the basic domain of the conserved domain of bZIP proteins, whereas motif 5 is the leucine zipper domain. The two domains are closely associated with each other and determine the binding of bZIP protein to DNA. Moreover, it was found that most of the motifs showed specific distribution in different subfamilies, and most groups had specific motifs different from other groups. Therefore, we speculate that the genes in the same subfamily may have similar functions, whereas the genes in different subfamilies contain specific conservative motif and have a different arrangement order, and thus, they may perform different biological functions. According to the classification results of A. thaliana, mung bean bZIP family proteins were divided into 10 groups: A~I and S. Group A had the most members, with a total of 17 accounting for 22.7%; followed by group D and group S, both with 14 bZIP proteins accounting for 18.7%; and group I, with nine members accounting for 12%. Group B had the least members, with only one bZIP protein accounting for 1.3%. The largest group in A. thaliana was group S accounting for 22.7%; followed by group A and group I, both accounting for 17.3%; and group D accounting for 13.3%. It was evident that the groups with the largest number of bZIP members were A, D, S, and I, with these four groups containing 72.1% and 70.6% of the total bZIP members in mung bean and A. Thaliana, respectively.

Analysis of Conserved Motif and Gene Structure of bZIP Proteins in Mung Bean
MEME software was applied to determine the conserved motifs of mung bean bZIP family proteins, with obtained results indicating 20 conserved motifs. The conserved motif, gene structure, and phylogenetic relationship of bZIP proteins were then visually analyzed by TBtools software (Figure 3).   Based on the analysis of the structure of bZIP family genes in mung bean, it was found that the bZIP genes clustered in the same group had similar gene structure ( Figure 3A,C). The number of exons of bZIP family gene in mung bean ranged between one and 18. Group S had the least number of exons, where each gene only had one exon and no intron, with the exception of VrbZIP39, which had two exons and one intron. VrbZIP26, a member of group F, also had only one exon and no intron. On the other hand, group D and G had the greatest number of exons, with group D exons ranging between 9 and 15 and group G exons ranging between 6 and 18.

Analysis of Promoter Cis-Acting Elements of bZIP Genes in Mung Bean
In this study, ten cis-acting elements were selected in the 2000 bp upstream CDS of bZIP genes, and a total of 633 elements were obtained (Figure 4). Among them, there were five kinds of hormone response elements, including ERE (ethylene response element), ABRE (abscisic acid response element), GARE-motif (gibberellin response element), TCAelement (salicylic acid response element), and AuxRR-core (auxin response element). There were five cis-acting elements associated with defense and stress response: MBS (MYB binding element), TC-rich repeats (defense response element), LTR (low temperature response element), Wbox (trauma and pathogen response element), and WUN-motif (trauma response element) ( Figure 4A). In addition, there were 410 hormone response elements accounting for 64.8% of all elements and 223 defense and stress response elements accounting for 35.2% of all elements. Results showed that VrbZIP38 contained the largest number of elements (18), whereas VrbZIP48 had the least number of elements (only one component). Among the 75 bZIP genes, 74 genes contained 1 to 16 hormone response elements, among which ethylene response elements and abscisic acid response elements were the most common. Moreover, 70 of the 75 genes contained defense and stress response elements, with the number ranging from 1 to 10 ( Figure 4B).

Analysis of Expression Patterns of bZIP Genes in Mung Bean under Salt Stress
To analyze the role of bZIP family genes in mung bean response to salt stress, we performed high-throughput RNA sequencing analysis on mung bean leaves under salt stress and then used the "pheatmap" package in R software to generate the expression heat map of bZIP genes ( Figure 5). Given that VrbZIP27 was not expressed under normal

Analysis of Expression Patterns of bZIP Genes in Mung Bean under Salt Stress
To analyze the role of bZIP family genes in mung bean response to salt stress, we performed high-throughput RNA sequencing analysis on mung bean leaves under salt stress and then used the "pheatmap" package in R software to generate the expression heat map of bZIP genes ( Figure 5). Given that VrbZIP27 was not expressed under normal conditions and salt stress, we analyzed the expression pattern of the other 74 mung bean genes. Figure 5 shows that bZIP genes of mung bean could respond to salt stress, and there were differences in the expression patterns of different bZIP genes under salt stress. It was found that the expression patterns of bZIP genes in the same group were similar, but there were differences among different groups. Results showed that VrbZIP1, VrbZIP 24, VrbZIP 35, and VrbZIP 53 genes in group A, which has 17 genes, were downregulated under salt stress. VrbZIP 23, VrbZIP 59, VrbZIP 68, and VrbZIP 74 had similar expression patterns, where their expression increased at 6 h and 12 h and then decreased at 24 h. There was no significant difference in the expression levels of VrbZIP 12, VrbZIP 42, VrbZIP 51, and VrbZIP 75 at 6 h and 12 h time points, but the levels increased significantly at 24 h. In addition, VrbZIP 14 and VrbZIP 71 as well as VrbZIP 37 and VrbZIP45 exhibited the same expression pattern, respectively. The three genes in group E had similar expression patterns. Therefore, considering that the expression patterns of bZIP genes in the same group could be the same, it was evident that bZIP genes in the same group may perform the same biological functions under salt stress.
ife 2022, 12, x FOR PEER REVIEW11 of 17 pression patterns. Therefore, considering that the expression patterns of bZIP genes in the same group could be the same, it was evident that bZIP genes in the same group may perform the same biological functions under salt stress.

Screening and Metabolic Pathway Analysis of Differentially Expressed bZIP Genes in Mung Bean under Salt Stress
According to the RNA-seq analysis results, genes with fold change ≥ 2 and FDR < 0.01 were regarded as differentially expressed genes. A total of 24 differentially expressed bZIP genes were identified ( Figures 6A,B), of which VrbZIP46, VrbZIP37, VrbZIP45, VrbZIP41, VrbZIP22, VrbZIP38, VrbZIP49, VrbZIP12, VrbZIP11, VrbZIP34, and VrbZIP70 were upregulated at different time points under stress. Enrichment analysis showed that VrbZIP12, VrbZIP37, and VrbZIP45 genes were enriched in the plant hormone signal transduction pathway mainly in response to the induction of abscisic acid ( Figure 6C), which activated the expression of abiotic stress-related gene ABF and then regulated the stomatal closure of leaves.

Screening and Metabolic Pathway Analysis of Differentially Expressed bZIP Genes in Mung Bean under Salt Stress
According to the RNA-seq analysis results, genes with fold change ≥ 2 and FDR < 0.01 were regarded as differentially expressed genes. A total of 24 differentially expressed bZIP genes were identified ( Figure 6A,B), of which VrbZIP46, VrbZIP37, VrbZIP45, VrbZIP41, VrbZIP22, VrbZIP38, VrbZIP49, VrbZIP12, VrbZIP11, VrbZIP34, and VrbZIP70 were upregulated at different time points under stress. Enrichment analysis showed that VrbZIP12, VrbZIP37, and VrbZIP45 genes were enriched in the plant hormone signal transduction pathway mainly in response to the induction of abscisic acid ( Figure 6C), which activated the expression of abiotic stress-related gene ABF and then regulated the stomatal closure of leaves.

Responses of bZIP Genes to Abiotic Stress
Eight bZIP genes were randomly selected and their expression patterns were analyzed in roots and leaves of mung bean seedlings under drought, low-temperature, alkali, and salt stress (Figure 7). Under drought stress, almost all of the eight bZIP genes were upregulated in different degrees in mung bean leaves compared to the control, and the expression level was the highest after 24 h under stress. However, most genes were not significantly upregulated or downregulated in roots. Under alkali stress, VrbZIP39, VrbZIP45, VrbZIP46, VrbZIP49, and VrbZIP68 were all upregulated in leaves, with the expression of VrbZIP46 increasing obviously at 6 h and 24 h time points. In the roots, VrbZIP11, VrbZIP39, VrbZIP46, and VrbZIP68 were first upregulated and then downregulated, whereas there was no significant change was detected in the expressions of VrbZIP34, VrbZIP45, and VrbZIP49. Under low-temperature stress, VrbZIP46 showed an upregulated expression pattern in leaves, whereas VrbZIP34, VrbZIP39, VrbZIP55, and VrbZIP68 were all downregulated. In roots, VrbZIP11 and VrbZIP39 were upregulated, VrbZIP46 was obviously downregulated, and there was no obvious change in the expressions of VrbZIP55 and VrbZIP68.

Responses of bZIP Genes to Abiotic Stress
Eight bZIP genes were randomly selected and their expression patterns were analyzed in roots and leaves of mung bean seedlings under drought, low-temperature, alkali, and salt stress (Figure 7). Under drought stress, almost all of the eight bZIP genes were upregulated in different degrees in mung bean leaves compared to the control, and the expression level was the highest after 24 h under stress. However, most genes were not significantly upregulated or downregulated in roots. Under alkali stress, VrbZIP39, VrbZIP45, VrbZIP46, VrbZIP49, and VrbZIP68 were all upregulated in leaves, with the expression of VrbZIP46 increasing obviously at 6 h and 24 h time points. In the roots, VrbZIP11, VrbZIP39, VrbZIP46, and VrbZIP68 were first upregulated and then downregulated, whereas there was no significant change was detected in the expressions of VrbZIP34, VrbZIP45, and VrbZIP49. Under low-temperature stress, VrbZIP46 showed an upregulated expression pattern in leaves, whereas VrbZIP34, VrbZIP39, VrbZIP55, and VrbZIP68 were all downregulated. In roots, VrbZIP11 and VrbZIP39 were upregulated, VrbZIP46 was obviously downregulated, and there was no obvious change in the expressions of VrbZIP55 and VrbZIP68.

Discussion
This study identified 75 bZIP genes in mung bean, which is the same number as in Arabidopsis thaliana (75) [6] but lower than in soybean (131) [20], rice (89) [17], corn (125) [4], wheat (156) [22], and millet (85) [23]. Notably, the number is greater than that in potatoes (56) [24]. According to systematic evolution analysis of bZIP genes in Arabidopsis and mung bean, it was found that bZIP family genes in mung bean could be divided into 10 groups, similar to Arabidopsis. The intron-exon structure of bZIP genes in mung bean was highly conserved in the same group and had similar intron region distribution. Previous studies reported that the structure of intron-exon is closely associated with the evolution of gene family [25][26][27][28]. In the present study, the number of exons of bZIP genes in mung bean ranged from 1 to 18. Among the 75 genes, 14 genes did not have introns, of which 13 were group S members, and 1 was in group F. Groups D and G had the largest number of exons, with 9-15 exons in group D and 6-18 exons in group G. Similar gene structural diversity has also been found in bZIP family genes in beans and potatoes [20,24], which suggests that different splicing states of exons and introns may be of great significance to the evolution of bZIP genes in mung beans, and there are similarities in the evolution of bZIP genes among different species.
This study also identified and classified 20 conserved motifs of mung bean bZIP proteins. Results showed that all mung bean bZIP proteins contained typical bZIP domain (motif 1), with each group having some common motifs with other groups. Most groups also contained some specific motifs. It is worth mentioning that the bZIP domain, the core of the bZIP family proteins, preferentially binds to the specific cis-acting elements of the downstream target gene promoter (such as ABRE). Liu and Chu (2015)

Discussion
This study identified 75 bZIP genes in mung bean, which is the same number as in Arabidopsis thaliana (75) [6] but lower than in soybean (131) [20], rice (89) [17], corn (125) [4], wheat (156) [22], and millet (85) [23]. Notably, the number is greater than that in potatoes (56) [24]. According to systematic evolution analysis of bZIP genes in Arabidopsis and mung bean, it was found that bZIP family genes in mung bean could be divided into 10 groups, similar to Arabidopsis. The intron-exon structure of bZIP genes in mung bean was highly conserved in the same group and had similar intron region distribution. Previous studies reported that the structure of intron-exon is closely associated with the evolution of gene family [25][26][27][28]. In the present study, the number of exons of bZIP genes in mung bean ranged from 1 to 18. Among the 75 genes, 14 genes did not have introns, of which 13 were group S members, and 1 was in group F. Groups D and G had the largest number of exons, with 9-15 exons in group D and 6-18 exons in group G. Similar gene structural diversity has also been found in bZIP family genes in beans and potatoes [20,24], which suggests that different splicing states of exons and introns may be of great significance to the evolution of bZIP genes in mung beans, and there are similarities in the evolution of bZIP genes among different species.
This study also identified and classified 20 conserved motifs of mung bean bZIP proteins. Results showed that all mung bean bZIP proteins contained typical bZIP domain (motif 1), with each group having some common motifs with other groups. Most groups also contained some specific motifs. It is worth mentioning that the bZIP domain, the core of the bZIP family proteins, preferentially binds to the specific cis-acting elements of the downstream target gene promoter (such as ABRE). Liu and Chu (2015) demonstrated that different protein motif composition may lead to functional diversity of mung bean bZIP proteins [19]. The gene structure and protein motif distribution of each phylogenetic group were highly conserved, which shows that the genes in the same group are closely associated with the process of evolution.
Furthermore, ten elements were identified in the promoter region of mung bean bZIP genes, and these elements could be divided into two categories: hormone response elements and stress response elements, accounting for 64.8% and 35.2%, respectively. Among the hormone response elements, the number of ERE and ABRE elements was the largest, accounting for 86.1% of the hormone response elements. In addition, 70 of the 75 bZIP genes contained stress response elements. These results suggest that mung bean bZIP genes can participate in the regulation pathway induced by ethylene and abscisic acid, and the existence of a large number of stress response elements lays a foundation for mung bean bZIP genes to participate in various stress responses.
At present, numerous studies have confirmed that bZIP genes play an important role in response to salt and other abiotic stresses in many species [16,[29][30][31][32][33]. For multi-gene families, such as bZIP, the analysis of gene expression can provide an important theoretical basis for gene function prediction. The transcriptome sequencing data indicated that mung bean bZIP genes could respond to salt stress in different degrees. Some bZIP genes in the same group had similar expression patterns, such as the three bZIP genes of group E. These results suggest that bZIP genes in the same group may perform the same biological function under salt stress. Studies have shown that PmbZIP genes with the same motifs exhibited similar expression patterns: PmbZIP31, PmbZIP36, and PmbZIP41 played more prominent roles in response to freezing stress [34]; 52 putative bZIP genes of P. Glaucum were identified, and 9 of them differentially expressed under water stress conditions. These are consistent with our conclusions [35]. It should be noted that some genes in the same group showed varying expression patterns. For example, four bZIP genes in group A were downregulated, and four genes were upregulated at 6 h and 12 h and then downregulated at 24 h. The bZIP gene family is a critical member of the ABA signaling pathway in plants, which plays an important role in plants response to drought [36]. Enrichment analysis of the 24 differentially expressed VrbZIP genes in RNA-seq analysis showed that VrbZIP12, VrbZIP37, and VrbZIP45 were annotated in the plant hormone signal transduction pathway, which mainly activates the expression of ABF, an abiotic stress-related gene, in response to abscisic acid induction, followed by regulation of stomatal closure in mung bean leaves. Studies have shown that group A bZIP genes in A. thaliana can activate the expression of abscisic-acid-dependent stress-resistance gene ABF under abiotic stress [37]. Accumulating evidence suggests that ABF is mainly involved in the response to salt, drought, hightemperature, and low-temperature stresses [38,39], which is consistent with our results. Importantly, two genes (VrbZIP12 and VrbZIP37) that regulate the expression of ABF gene belong to group A.
Finally, we analyzed the expression of bZIP genes in roots and leaves of mung bean seedlings under drought, low-temperature, and alkali stress. Results showed that the expression patterns of bZIP genes were different under different stress conditions. For example, bZIP68 was upregulated under drought and alkali stress but downregulated under low-temperature stress. The expression patterns of the same genes were also different in roots and leaves. Most of the bZIP genes selected in this study were upregulated in leaves under the three kinds of stress but downregulated in roots, indicating that the expression of mung bean bZIP genes under stress is tissue-specific, and the genes may have different regulatory mechanisms in different tissues.

Conclusions
Seventy-five mung bean bZIP genes were identified, the bZIP genes were unevenly distributed in mung bean chromosomes (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11), and all had highly conserved bZIP domains. Phylogenetic analysis divided the bZIP genes into 10 groups, with genes in the same group having similar gene structure and protein sequence. Results revealed that the promoter region contained more abscisic acid elements, ethylene response elements, and stress-response elements. Although the mung bean bZIP genes could respond to salt stress in different degrees, the expression patterns of the genes were different at different time points under stress. Moreover, it was evident that the bZIP genes could respond to drought, low-temperature, and alkali stress, and the expression patterns of the genes were tissue-specific. Metabolic pathway enrichment analysis of differentially expressed genes showed that VrbZIP12, VrbZIP37, and VrbZIP45 were mainly enriched in the plant hormone signal transduction pathway, which regulates the expression of abiotic stress-related gene ABF. Therefore, the three genes can be used as important candidate genes for further gene function verification.