Genome-wide analysis of the bZIP gene family in Chinese jujube (Ziziphus jujuba Mill.)

Background Among several TF families unique to eukaryotes, the basic leucine zipper (bZIP) family is one of the most important. Chinese jujube (Ziziphus jujuba Mill.) is a popular fruit tree species in Asia, and its fruits are rich in sugar, vitamin C and so on. Analysis of the bZIP gene family of jujube has not yet been reported. In this study, ZjbZIPs were identified firstly, their expression patterns were further studied in different tissues and in response to various abiotic and phytoplasma stresses, and their protein-protein interactions were also analyzed. Results At the whole genome level, 45 ZjbZIPs were identified and classified into 14 classes. The members of each class of bZIP subfamily contain a specific conserved domain in addition to the core bZIP conserved domain, which may be related to its biological function. Relative Synonymous Codon Usage (RSCU) analysis displayed low values of NTA and NCG codons in ZjbZIPs, which would be beneficial to increase the protein production and also indicated that ZjbZIPs were at a relative high methylation level. The paralogous and orthologous events occurred during the evolutionary process of ZjbZIPs. Thirty-four ZjbZIPs were mapped to but not evenly distributed among 10 pseudo- chromosomes. 30 of ZjbZIP genes showed diverse tissue-specific expression in jujube and wild jujube trees, indicating that these genes may have multiple functions. Some ZjbZIP genes were specifically analyzed and found to play important roles in the early stage of fruit development. Moreover, some ZjbZIPs that respond to phytoplasma invasion and abiotic stress environmental conditions, such as salt and low temperature, were found. Based on homology comparisons, prediction analysis and yeast two-hybrid, a protein interaction network including 42 ZjbZIPs was constructed. Conclusions The bioinformatics analyses of 45 ZjbZIPs were implemented systematically, and their expression profiles in jujube and wild jujube showed that many genes might play crucial roles during fruit ripening and in the response to phytoplasma and abiotic stresses. The protein interaction networks among ZjbZIPs could provide useful information for further functional studies.


Background
Transcription factors (TFs) are an important component of regulatory networks. TFs bind to specific promoter sequences to activate or inhibit the expression of target genes [1]. Among several TF families unique to eukaryotes, the basic leucine zipper (bZIP) family is one of the largest and most diverse [2][3][4], containing two regions with different functions: the basic region and the leucine zipper region [5]. The basic region is highly conserved and consists of approximately 16 amino acid residues with a consistent N-X7-R/K motif for nuclear localization and sequentially specific DNA binding, while the leucine zipper region is less conserved and forms a helical structure with dimerization specificity [1,6].
Members of the bZIP family are involved in the regulation of plant resistance under biotic and abiotic stresses [7][8][9], and play some important roles during plant growth and development processes, such as hormone signal transduction [10], energy metabolism [11], seedling development [12] and flowering [13]. In plant, bZIPs can be combined with cis-acting elements such as G-box (CACGTG), A-box (TACGTA) and abscisic acid (ABA)responsive elements (ABRE) (CCACGTGG) to regulate the expression of downstream genes. In Arabidopsis, ABF1, ABF2, ABF3 and ABF4 could bind to the cis-acting element ABRE, and regulate many downstream salt and drought tolerances through the interaction between ABRE and bZIP proteins [14]. OsbZIP1 might enhance resistance to Magnaporthe grisea through salicylic acid (SA), jasmonic acid (JA) and ABA signal transduction pathways [15]. In addition, AtbZIPs also regulated the signal transduction of ABA-related pathways, thereby affecting seed germination and maturity [12].
The TGA (TGACG motif-binding factor) subfamily of bZIPs plays important roles in defense responses against pathogens [16]. As the target of SA signaling, the TGA factors can thus activate and connect the SA pathway with the JA/ET-dependent pathway [17][18][19][20]. Moreover, these WRKY transcription factors induced by SA could activate the promoter of pathogenesis-related (PR)-1 and involved in the regulation of the TGA/NPR1 complex [21,22]. And they also demonstrated some reverse functions, such as TGA2 suppressed the expression of PR1 whereas TGA6 was able to increase PR1 expression and could induce basal defense [23]. Thus, such interactions will be discovered in more plants when further exploring the function of TGA factors.
Chinese jujube (Ziziphus jujuba Mill.), a member of the Rhamnaceae family, is an important dry fruit and a traditional herbal medicine in Asia [24]. Both Chinese jujube and wild jujube, which are considered ideal fruit trees in arid and semiarid temperate regions, have strong tolerance to biotic and abiotic stress [25,26]. Although the bZIP gene family has been analyzed in many other plant species such as Arabidopsis [27], peach [28], apple [3] and so on [29,30], the analysis of bZIP family in jujube has not yet been reported. Based on the functions of bZIPs in other species, we thought the members of bZIP family should have multiple functions on jujube development, defense responses against pathogens and abiotic stress. Thus, the characteristics and expressions of bZIP members in Chinese jujube are identified and analyzed systematically. These results would provide a basis for future studies related to biological functions and the regulatory networks.

Identification of ZjbZIPs in Chinese jujube
A total of 45 nonredundant putative bZIP transcripts (Table 1) within the jujube genome sequence were identified. They were named ZjbZIP1 to ZjbZIP45 according to their gene structure and motifs. The ORF length of the ZjbZIPs ranged from 384 bp (ZjbZIP44) to 2205 bp (ZjbZIP11), and they encoded proteins ranging from 127 (ZjbZIP44) to 749 (ZjbZIP11) amino acids (aa) in length, with predicted pIs ranging from 4.65 (ZjbZIP40) to 9.96 (ZjbZIP44) ( Table 1) and predicted molecular weights (MWs) of 14.67-82.15 kDa. The proteins with an isoelectric point greater than 7 accounted for 53% of the total number, which means that half of the ZjbZIP proteins were neutral or alkaline, and most of the proteins in the E, F and G subfamilies were weakly acidic.
The average GC content of 45 ZjbZIPs was 46.88%, and the contents of GC 1 , GC 2 and GC 3 were 53.45, 44.63 and 42.56%, respectively. Relative Synonymous Codon Usage (RSCU) analyses will help us to understand the patterns in ZjbZIPs, and the RSCU values greater than 1.5 was defined as high-frequency codons [31]. Among the 64 codons of 45 ZjbZIPs, seven highfrequency codons (AGG 1.80, AGA 1.79, GTT 1.73, GCT 1.69, TTG 1.65, TCT 1.59 and TTT 1.58) were investigated, and most of them were T-ended (Table 2). We also found that most ZjbZIPs prefer ATG as the stop codon. RSCU values of four NCG codons in ZjbZIPs were the lowest (TCG 0.62, CCG 0.61, ACG 0.45, GCG 0.36) ( Table 2), suggesting that ZjbZIPs were at a relative high methylation level [32]. Meanwhile, RSCU values of four NTA codons also displayed a lower level (ATA 0.75, TTA 0.73, GTA 0.66, CTA 0.55), which was beneficial to increase the protein production by inhibiting mRNA degradation [33].

Phylogenetic tree construction and conserved motifs of ZjbZIPs
Compared with bZIPs in Arabidopsis, ZjbZIPs were also divided into 10 subfamilies (A-I, S). In addition, we defined six newly discovered ZjbZIPs as four subfamilies of J, K, L and M (Fig. 1). The classification result was further supported by the phylogenetic tree of bZIPs between jujube and apple (Additional file 1). In jujube, the A subfamily is the largest subfamily, while in the Arabidopsis, the S subfamily is the largest subfamily [34]. As expected, bZIP proteins from the same group tended to cluster together and were named following the same scheme.
Using MEME software, a total of 10 conserved motifs among ZjbZIPs was identified (Fig. 2), of which motif 1 and motif 5 were identified as the core conserved domains and constituted the leucine zipper region of bZIP (Additional file 2). The proteins in each subfamily contain the same conserved motifs, which further support the above result of phylogenetic tree. However, they also have different conserved motifs among various subfamilies. For example, in addition to the core conserved domains, the A subfamily also contains three conserved motifs (Motif 6, 8 and 9), which may be related to their different biological functions.

The phylogenetic tree and line charts for a lineage of gene groups for ZjbZIP genes
In order to further study the evolution pattern and direction of ZjbZIP genes, ZjbZIP26 and 29 were selected to perform evolutionary analysis. ZjbZIP26 and ZjbZIP29 are homologous genes of GBF3 and HY5, respectively. And GBF3 and HY5 were proved to participate in various biological processes [35,36]. As shown in the Fig. 3, two phylogenetic trees of ZjbZIP26 and − 29 with 20 other genes showing high homology indices (HIs) in various species were constructed, respectively. The values of HIs between all pairs in the two trees were all above 0.7, suggesting that they have similar amino acid sequences and might have conserved functions. To the phylogenetic tree of ZjbZIP29, three paralogous events were presumably occurred in a group of two genes (Cucumis sativus XP_004138731 and Cucumis melo NP_001284656), a group of two genes (Citrus clementina XP_006450470 and Citrus sinensis XP_006483336) and a group of two genes (Ziziphus jujuba XP_015885857 and XP_015868446) in Ziziphus jujuba. To the phylogenetic tree of ZjbZIP26, there were also five paralogous events.
Generally, along with evolutionary time the decrease of both the numbers of genes and species represents an orthologous event, and the number of genes decreases and the number of species is retained, which means a paralogous event [37]. For ZjbZIP26, along with evolutionary time, the numbers of genes (red line) and species (blue line) were both decreased in the timing between 0.619 and 0.632 of HIs, suggesting that an orthologous event happened; and only the numbers of genes (red line) were decreased in the timing between 0.570 and 0.582 of HIs, indicating that a paralogous event occurred. The two kinds of homologous events also found in ZjbZIP29 analysis. Therefore, the paralogous and orthologous duplication events should occurred during the evolutionary process of bZIPs.

The chromosomal location and gene structure of ZjbZIPs
Of the 45 ZjbZIP genes, 34 were mapped to 10 pseudochromosomes in the jujube genome (Fig. 4), and 11 genes were unanchored. ZjbZIPs were not evenly distributed across the 12 chromosomes: there were 6 ZjbZIPs on both Chr. 9 and 12, and no genes were located on Chr. 1 or 11. Further analysis found that ZjbZIP1 and ZjbZIP2, ZjbZIP6 and ZjbZIP19, ZjbZIP7 and ZjbZIP14 are tandem repeat genes, indicating that some ZjbZIPs underwent gene duplication during jujube evolution to increase the number of genes and enhance their biological functions.
Additionally, the gene structure within the same subfamily was highly conserved (Additional file 3). We found that the genes in the C, D and G subfamilies contained more introns than did the genes in the other groups.

Expression patterns of ZjbZIPs in various organs
To investigate the organ-specific expression of the ZjbZIP genes, the expression of 30 ZjbZIPs were analyzed in five organs of jujube and wild jujube by RT-PCR (Fig. 5, Additional file 4). It was shown that most ZjbZIPs were expressed in at least four organs, indicating that ZjbZIPs should involve in the development process of various organs in jujube and wild jujube. The expression of most ZjbZIP genes in the same subfamily showed similar patterns, suggesting that these genes in the same subfamily might have conserved functions. In wild jujube, some ZjbZIP genes display a special expression pattern, in which ZjbZIP10, − 25, − 31, − 36 were highly expressed in branch and leaf, indicating that these genes should play some roles in the development of the two organs. Compared with jujube, the expression levels of some genes, such as ZjbZIP3, − 5, − 12, − 34, − 35, − 36, − 38, and − 41, were significantly decreased in root of wild jujube, and only three genes, ZjbZIP15, − 31, − 40, showed lower expression (Fig. 5). These ZjbZIP genes may be related to their differing functions in root between wild jujube and jujube. We also found most of ZjbZIPs were expressed with varying degrees in fruit, and these genes should be candidate genes for further investigating in jujube fruit development. The broad and divergent expression patterns indicated that the ZjbZIPs should have multiple functions in jujube growth and development.

ZjbZIPs involvement in jujube fruit development
Based on organ-specific expression (Fig. 5), the highly expressed genes in fruit were further investigated at five development stages of jujube fruit. It is noteworthy that almost all of the ZjbZIPs tested were highly expressed at the first two stages and showed similar trends in both of the two cultivars (Fig. 6), indicating that these genes were positively involved in the fruit enlargement process.
In other word, ZjbZIPs should play some significant roles in jujube fruit development.
Especially, ZjbZIP28, − 29, − 30, − 36, − 38 showed the obvious increase in expression at the early white mature fruit stage (EWM), which period is just the fruit rapid expanding stage. These genes were identified as candidate genes for further study, in terms of their functions in fruit development.

ZjbZIP protein-protein interaction network prediction
Based on their orthologs in Arabidopsis, ZjbZIPs were predicted to interact with each other by STRING (Fig. 7a), which was in accordance with previous reports that the binding activity of bZIPs depends upon the formation of homodimers or heterodimers among bZIPs [38]. Several interactions including ZjbZIP27 and ZjbZIP28, ZjbZIP29 and ZjbZIP28, were further confirmed by yeast two-hybrid in Fig. 7b (Additional file 4). As shown in Fig. 7a, both TGA9 (homolog of ZjbZIP14) and TGA10 (homolog of ZjbZIP18) are involved in the regulation of anther development [39], B202H1 (homolog of ZjbZIP12) can interact with BZIP53 (homolog of ZjbZIP36, ZjbZIP38 and ZjbZIP44) to promote their expression in seeds and are involved in defense responses [40], and the interaction between GBF4 (homolog of ZjbZIP4) and bZIP68 (homolog of ZjbZIP25) is regulated by light or other hormones [41].

The expression of bZIP genes under abiotic and biotic stress
According to protein-protein interaction predictions and previous studies, 11 ZjbZIPs were selected for further validation under abiotic stress, including lowtemperature and salt stresses (Additional file 5). Under salt stress, the expression levels of ZjbZIP10, 11, 23 and 40 (Additional file 5) increased with prolonged treatment time. The expression of ZjbZIP10 increased significantly after 6 h of treatment and peaked at 48 h (Fig. 8a), and ZjbZIP40 expression peaked at 1 h of treatment but then decreased. The expression of ZjbZIP3 showed no difference under salt stress. We also found that the expression of ZjbZIP29, ZjbZIP36 and 38 increased markedly under low-temperature stress (Additional file 5 and Fig. 8b), and the expression of other three genes (ZjbZIP17, 22 and 24) also increased in varying degrees (Additional file 5). Jujube withes' broom (JWB), which is caused by phytoplasma, is a destructive disease that affects jujube production. Since bZIP genes have multiple functions in plants, whether they participate in jujube-phytoplasma interactions remains unclear. Among the 17 ZjbZIPs identified (Fig. 9), the expression of ZjbZIP3 and 12 in diseased leaves was significantly lower than that in healthy leaves. Moreover, the ZjbZIP11, 15,17,18,19,20 and 26 genes were highly expressed in diseased leaves. These results suggested that these ZjbZIPs participate in jujube-phytoplasma interactions and play different roles.

Discussion
In this study, a total of 45 bZIP genes were identified in the jujube genome. A previous genome evolution study Fig. 6 Heat maps of the relative expression of ZjbZIPs during fruit ripening in 'Lizao' and 'Yazao'. Y, young fruit; EWM, early white mature fruit; WM, white mature fruit; HR, half-red fruit; FR, full red fruit. Scaled log2 expression values based on qRT-PCR data are shown from blue to yellow, indicating low to high expression showed that Chinese jujube is closely related to species of the Rosaceae family [24], so the numbers of bZIPs from two Rosaceae species (apple and peach) [42], grape and Arabidopsis were compared with the number of ZjbZIPs. Compared with the number in some other plant species, the smaller number of bZIP genes in jujube and peach should be caused by the occurrence of only one whole genome duplication (WGD) event during the species evolution. Meanwhile, the divergence of gene number in same gene family of various species also relate to their evolutionary differences or the genome size [28].
In Arabidopsis, various subfamilies of AtbZIPs contain different conserved domains and perform various biological functions. For example, most members of the class A subfamily are involved in ABA pathway and abiotic stress responses, the class D members mainly play some roles in pathogen resistance and plant development, and the class S subfamily is involved in carbohydrate metabolism [34]. Here, the ZjbZIPs are highly homologous to sequences of the corresponding subfamilies in Arabidopsis, indicating that ZjbZIPs also might play different functions as a result of different conserved domains.
There are 10 TGA family members in Arabidopsis (Additional file 6), which are divided into five categories.
A number of studies have shown that members of the TGA family play an important role in defense against biological and necrotizing pathogens [43,44]. One member of the TGA family was found to interact with an ankyrin-repeat protein, a nonexpressor of the pathogenassociated (PR) gene (NPR1), which is a key component of the SA defense signaling pathway; SA is a key signaling molecule involved in plant resistance [45,46]. Moreover, TGA members participate in mitotic reactions, regulate the growth and development of organisms, and play an important role in flower development [43]. Homology analysis and phylogenetic tree analysis indicated that subfamily D of ZjbZIP family is highly homologous to the TGA family in Arabidopsis (Additional file 6). It showed that ZjbZIP14 and TGA9 (At1g08320), ZjbZIP18 and TGA10 (At5g06839), ZjbZIP19 and TGA4 (At5g10030) were homologous (Additional file 6), suggesting that they might have similar biological functions and participate in pathogen defense reactions. Moreover, ZjbZIP17, − 18, − 19 and − 20, which are highly homologous to the TGA members in Arabidopsis, showed differential expression at the early stage of phytoplasma infection (Fig. 9). Phytoplasma causes one of serious disease in jujube production, named Jujube Witches' Broom. Thus, these genes should involve in jujubephytoplasma interaction. AtbZIPs can recognize and bind to ABREs within promoters, which are named ABRE-binding factors (ABFs)/ ABRE-binding proteins (AREBs) [9]. They play a key role in the regulation of the expression of downstream stress-responsive genes involved in ABA signaling. Genetic transformation of ABF/ABRE transcription factors has been suggested to be an effective approach for engineering stress-tolerant plants [47]. In this study, we also found that the expression level of ZjbZIP10 was increased with time under salt stress. Thus, ZjbZIP10 may have the function in terms of involvement in the resistance to salt stress. In addition, the expression of ZjbZIP36 was significantly induced under lowtemperature stress, indicating that it might be involved in the response to cold stress.
Overall, all the results imply that ZjbZIPs play multiple roles in jujube development and under biotic and abiotic stresses. Further function verifications are worthwhile and need to reveal their regulation metabolism in detail.

Conclusions
The bioinformatics analyses of 45 ZjbZIPs were firstly investigated in this study. Meanwhile, their expression patterns in jujube and wild jujube were measured by qPCR during fruit ripening and in the response to phytoplasma and abiotic stresses, and some candidate genes involved in these processes were screened out. The protein interactions among ZjbZIPs were predicted, which provide useful information for further functional studies to elucidate their regulation mechanism.

Plant materials
The Chinese jujube and wild jujube trees used in this study were cultivated at the Experimental Station of Chinese Jujube, Hebei Agricultural University, Baoding, China. They are very common fruit trees in China and are not endangered species. No specific permits were required for the sample collection.
Five organs (roots, branches, leaves, flowers and fruits) collected from three jujube trees and three wild jujube trees were used for organ-specific expression analysis. In general, jujube fruit development can be divided into five typical stages: the young fruit stage (Y), early white mature fruit stage (EWM), white mature fruit stage (WM), half-red fruit stage (HR) and full-red fruit stage (FR). The first two stages of jujube fruit development (young fruit stage (Y) and early white mature fruit stage (EWM)) are crucial to fruit enlargement, and the later three stages are important to fruit maturity. Z. jujuba Mill. 'Lizao' and 'Yazao' are two jujube fresh cultivars, which fruit display five typical developmental stages and suitable to study the fruit development. Thus, fruit samples of the two cultivars were taken at the five developmental stages and used to investigate the expression patterns of ZjbZIPs. Each treatment consisted of three biological replicates.
Three healthy trees and three JWB diseased trees of Z. jujuba Mill. 'Dongzao' were used to collect the leave samples at four stages (June, July, August and September). There are four kinds of leaves, i.e., healthy leaves (HL), apparently normal leaves (ANL), phyllody leaves (PL), and witches' broom leaves (WBL). All treatments consisted of three biological replicates. All fresh samples were frozen in liquid nitrogen immediately, and then stored at − 80°C for RNA extraction.
Cold-and salt-stress treatments were performed on the callus tissues of Z. jujuba Mill. 'Guanyangchangzao'. For the cold treatment, calli at same growth stage and condition were transferred at 4°C and then collected within 0, 1, 6, 12, 24, and 48 h. Calli incubated at 25°C were collected and used as negative controls [26]. For the salinity treatment, callus tissues were subjected to a 150 mM NaCl and NaHCO 3 -NaOH solution (pH 9.5) for 0, 1, 6, 12, 24, and 48 h. Samples subjected to sterile water treatment rather than the NaCl and NaHCO 3 -NaOH solution were used as negative controls [48]. Identification and protein structure analysis of ZjbZIPs A hidden Markov model (HMM) file of the bZIP domain (PF00170) was downloaded (Pfam 31.0, https:// pfam.xfam.org/) and used as the template to identify ZjbZIP sequences in the jujube genome [24]. Above sequences were further confirmed by using the Plant Transcription Factor Database (PlantTFDB, http:// planttfdb.cbi.pku.edu.cn/) and the online CD search tool in NCBI database (https://www.ncbi.nlm.nih.gov/). Some characters of the ZjbZIP genes were also predicted by the bioinformatical tools [49]. GC content and RSCU were performed by software codonW1.4.4 (http:// codonw.sourceforge.net/). The conserved motifs of ZjbZIP proteins were detected by MEME (http://memesuite.org/).

Multiple sequence alignment and phylogenetic tree construction
The multiple sequence alignment results were analyzed by using DNAMAN 8.0 software. A phylogenetic tree consisting of 45 ZjbZIPs was constructed based on their conserved domains. bZIP proteins from Arabidopsis thaliana and Malus were downloaded from the NCBI (https://www.ncbi.nlm.nih.gov/) and PlantTFDB (http:// planttfdb.cbi.pku.edu.cn/) (Additional file 7) and used to construct the phylogenetic tree [52,53]. In addition, the phylogenetic trees and line charts for a lineage of gene groups for ZjbZIP26 and 29 were analyzed in the Gcorn plant database (http://www.plant.osakafu-u.ac.jp/~kagiana/gcorn/p/) [37].

RNA isolation and expression analysis
The procedures of total RNA isolation, detection and cDNA synthesis were performed according to previous methods [49]. Both semiquantitative PCR and qRT-PCR were used to measure the expression of ZjbZIPs. The primer sequences were listed in Additional file 8. The PCR reaction system and procedures in the study were implemented as the method described [49,54,55]. Three biological replicates were analyzed for each treatment.
All statistical analyses were performed with SPSS software 17.0. Duncan's multiple range tests were used to assess differences among different treatment times. Different letters mean significant difference at 0.05 levels between the corresponding treatments.

Protein-protein interaction predictions
45 ZjbZIP protein sequences were used as the targets and their orthologs of Arabidopsis thaliana were appointed as references, and the STRING website (https://string-db.org/) was employed to predict proteinprotein interactions [49]. Yeast two-hybrid screening (Y2H) [56] was further applied to verify several interactions.