Genome-Wide Identification and Analysis of bZIP Family and their Expression in Response to Multiple Abiotic Stresses in Malus domestica

The basic leucine zipper proteins (bZIPs) are widely involved in the regulation of many biological processes by acting as transcriptional factors in plants [1]. bZIPs characteristically harbor a bZIP domain, which composed of a basic amino acids rich region (basic region) and a leucine zipper region. When bZIPs acting as transcription factor, basic region directly binds to negatively charged DNA and leucine zipper region formed homoand/or heterodimer to maintain the stability of the combination [2]. The basic region (also called DNA binding region) is highly conserved, often containing about 18 amino acid residues with a sequence like Asn-X7-Arg/Lys-X9, which is used as one distinct feature to identify bZIP proteins [3]. Leucine zipper region has a less conserved dimerization motif that made up of heptad(s) with repeated leucine or similar hydrophobic amino acids. The amino acids at the positions nearing interface of leucine zipper determine the dimerization stability and specificity [4,5].


Introduction
The basic leucine zipper proteins (bZIPs) are widely involved in the regulation of many biological processes by acting as transcriptional factors in plants [1]. bZIPs characteristically harbor a bZIP domain, which composed of a basic amino acids rich region (basic region) and a leucine zipper region. When bZIPs acting as transcription factor, basic region directly binds to negatively charged DNA and leucine zipper region formed homo-and/or hetero-dimer to maintain the stability of the combination [2]. The basic region (also called DNA binding region) is highly conserved, often containing about 18 amino acid residues with a sequence like Asn-X7-Arg/Lys-X9, which is used as one distinct feature to identify bZIP proteins [3]. Leucine zipper region has a less conserved dimerization motif that made up of heptad(s) with repeated leucine or similar hydrophobic amino acids. The amino acids at the positions nearing interface of leucine zipper determine the dimerization stability and specificity [4,5]. Apple (Malus domestica) is one of the most economically important woody plant and fruit crops in the world. It is cultivated worldwide and has great economic values and so there is considerable interest in identifying genes to improve apple horticultural characteristics, such as those that promote stress resistance. Recently, the divergence of the bZIP gene family in apple, strawberry and peach were reported. It showed that bZIP genes exist multiple modes of gene evolution after duplication [13]. Using iTAK-Plant Transcription factor and Protein Kinase Identifier and Classifier database, 112 apple bZIP members were obtained and divided into 11 groups [14]. Among 11 groups, one group have been investigated and found some of them were responsive to high salinity and drought, as well as to different phytohormones. Nevertheless, more questions need to be investigated further. For example, whether or not the apple bZIP members other than the above group respond to stress is still unknown. Besides, what is the relationship of the functional evolution of bZIP family with the gene structure and splicing patterns of intron? As the transcription factors of binding with DNA, whether do apple bZIP members possess dimerization properties?

Genome-Wide Identification and Analysis of bZIP Family and their Expression in Response to Multiple Abiotic Stresses in Malus domestica
To answer the above questions, in this investigation, 116 bZIP encoding genes were identified and classified on the basis of their evolutionary relationships. The intron position and splicing pattern in the basic and hinge region of apple bZIP gene family have been highly conserved during the course of evolution in each group. Detailed sequence analysis has been done to identify additional conserved motif outside the bZIP domain and predict their dimerization properties, which revealed the functional divergence and redundancy of apple bZIP family. Cis-acting elements related to abiotic stress tolerance and hormone signaling were massively found in the promoters of apple bZIP genes. The expressions of some randomly selected apple bZIP genes were quite responsive to the exogenous abscisic acid, salt, cold and osmotic treatments. Our results provide a prospective for the evolutionary history and general biological function, especially in tolerance to multiple abiotic stresses.

Genome-wide identification of bZIP proteins in apple
The genome and proteome sequences of Arabidopsis and rice were downloaded from the NCBI website. Two methods were used. For the first method, BLASTp was performed by the Bioedit program with an E value cutoff of 0.001 to search the Golden Delicious peptide database (GDR http://www.rosaceae.org/) using bZIP protein sequences in Arabidopsis and rice as queries. For the second method, Hmmer3.0 software was downloaded from the HMMER website (http://hmmer. janelia.org/) to perform a global search of the apple proteome. The HMM profile PF00170) of the bZIP domain was download from the pfam website http://pfam.xfam.orf/) [15]. All of the proteins that were obtained by the two methods were submitted to the InterPro Database http://www.ebi.ac.uk/interpro/) and Smart Database http://www.smart. embl-heidelberg.de/) to ensure the presence of bZIP domain.

Phylogenetic analysis of M.bZIP genes
116 M.bZIP full length protein sequences, 89 OsbZIP full length protein sequences and 72 AtbZIP full length protein sequences were executed multiple alignments by Clustal X 1.83 program [16]. Phylogenetic trees were constructed using phylip-3.695 by the Neighbor-Joining (NJ) method, and the bootstrap test carried out with 1000 replications. The picture of those trees was drawn using FigTree (Verson1.4.2).

Chromosomal locations and gene duplications for M.bZIP genes
Location data were retrieved from genome annotations downloaded from the Genome Database for Rosaceae (GDR http://www.rosaceae. org/). The chromosome map showing the physical location of all M.bZIP genes was generated using revised version of MapDraw [17]. Tandem duplications were characterized as adjacent genes located within 10 predicted genes apart or within 30 kbp of each other. The intragenomic homologous regions were obtained from SyMAP Synteny Browser (http://www.symapdb.org/).

Gene structures of M.bZIP transcription factors
To explore the diverse exon/intron organizations of M.bZIP genes, we compared the predicted coding sequences of M.bZIP genes with their corresponding genomic sequences using GSDS software (http:// gsds.cbi.pku.edu.cn) [18]. The intron distribution and splicing phase within the genomic sequences were further predicted based on the alignments of the genomic sequences and cDNA full-length sequences obtained by Spidey website (http://www.ncbi.nlm.nih.gov/spidey/). The isoelectric point (PI) of each protein was calculated by online ExPASy programs (http://www.expasy.org/tool/). The conserved motifs were identified by the online MEME analysis tool (http:// meme-suite. org/). The InterPro Database and Smart Database were used to identify conserved functional domains in apple M.bZIP proteins.

Promoter analysis of M.bZIP genes
To investigate cis-acting elements in the promoter sequences of M.bZIP genes, 2000 bp of genomic DNA sequences upstream of the transcriptional start site were obtained from the apple genome. Promoter sequences were submitted to the PLACE database (http://www.dna.affrc.go.jp/PLACE/) and cis-acting elements were obtained [19].

Sample preparation and total RNA extraction
Golden Delicious plants which were cultivated in the fields were selected as the experimental material for qRT-PCR analysis. The apple (Golden Delicious) seedlings were separately treated in 100 μM ABA and 20% PEG6000 for 0 h, 1 h, 3 h, 6 h and 9 h. 1 g leaf samples from three individual plants of stressed and controlled seedlings were harvested for total RNA extraction. The total RNA was isolated from the leaves using the CTAB procedure [20]. RNA concentrations and A 260 /A 280 ratios were determined using a NanoDrop Spectrometer (ND-1000 Spectrophotometer, peqlab). The integrity of the RNA was detected by agarose gel electrophoresis. Qualified RNA was used for cDNA synthesis and quantitative real-time PCR. All of materials were stored at -80°C after freezing in liquid nitrogen.

Quantitative Real-time PCR (qRT-PCR) analysis
cDNA fragments were synthesized from total RNA using the TransScripTMone-step gDNA Removal and cDNA Synthesis SuperMix TransGen Biotech, Beijing china). The gene-specific primer pairs were designed based on target gene sequences using the software Beacon Designers 8.10. The apple actin gene which from our laboratory was used as an internal control. The primer sequences were listed in Table S6. The qRT-PCR was carried out with a Stratagene Mx3000P thermocycler (Agilent) in a final volume of 20 μl that contained 1.4 μl cDNA, 10 μl 2*SYBR premix Ex TaqTM (Takara, Shiga, Japan), 7.8 μl H 2 O and 0.8 μl (10 μM) primers. The thermal cycling conditions were as follows: 44 cycles of 95°C denaturation for 15 s, 55°C annealing for 30 s and 72°C extension for 15 s. The real-time PCR experiment was carried out at least three times under identical conditions. The relative expression levels were calculated as 2 -(ΔCt of treatment−ΔCt of control) .

Identification, chromosomal location and expansion patterns of bZIP genes in apple
To identify bZIP proteins in apple genome, two methods were used. For the first method, all of the identified sequences of bZIP proteins from Arabidopsis and rice were used as queries to search the Golden Delicious peptide database (GDR http://www.rosaceae.org/) using software BLASTp [21]. Trying not to miss potential proteins and taking account of credible results, the E value cutoff was set to 0.001 as in references of similar investigations [22]. From this method, about 700 proteins were screened out. For the second method, Hidden Markov Model (HMM) profile of the bZIP domain from characteristics of the amino acid sequences established by the pfam website (http://pfam. xfam.org/) were used as a query to search in the Golden Delicious peptide database. Total 199 proteins were obtained from this method. To further identify whether all above proteins contains the bZIP domain, the sequences were submitted to the Interpro Database to be verified. Finally, 116 proteins in apple have been demonstrated to have bZIP domain were named as M.bZIP. Detailed information about each M.bZIP genes was showed in Table 1. Presumably these nine genes could be localized to unassembled genomic sequence scaffolds [22]. The number of M.bZIP genes on different chromosomes was variable. There were four chromosomes harboring low density of M.bZIP genes, chromosome 1, 5, 6 and 17 with only two M.bZIP genes on each chromosome. Other chromosomes had relatively high density and M.bZIP genes usually emerged with gene clusters on each chromosome. Further analyzed the distance between genes which located in the cluster, within 30 kb of each other was defined as tandem duplication gene [23]. These tandem duplicates genes were highlighted with green block as shown in Figure 1 and listed in Table  S1. Based on this result, we hypothesized that the tandem duplication of genes was one reason for the expansion of the M.bZIP genes. In addition to gene tandem duplication, a genome-wide duplication (GWD) event also occurred in apple and formed several intra-genomic homologous regions [22]. Large numbers of M.bZIP genes were located in homologous regions (gray shadow region in Figure 1) on different chromosomes, which suggest that genome duplication was another important way for expansion of M.bZIP genes.

Phylogenetic analysis and classfication of the M.bZIP family
Previous studies have shown that the evolution of bZIP genes was conservative in plants [24]. To obtain the credible phylogenetic relationship of the M.bZIP proteins, bZIP members from apple (116), Arabidopsis (72) and rice (89) were used together to construct a phylogenetic tree ( Figure 2). Comparing evolutionary grouping of Arabidopsis and rice from previous reports, the phylogenetic tree constructed in this investigation was credible. All M.bZIP members can be divided into six groups, designated as group I to VI, with 25, 35, 13, 24, 8 and 11 members in each group, respectively. The group II was the largest group with 35 members, and the group V was the least with only 8 members. Most of the M.bZIP genes exhibited closer relationship with AtbZIPs than OsbZIPs and they were hierarchical clustered in the same clade in the phylogenetic tree, which suggested that dicots of apple and Arabidopsis might had the similar patterns of their bZIP families in evolutionary.

Gene structure of M.bZIP encoding genes
As a result of key events in evolution, gene structure including the number and distribution of exon/intron, offers insights into understanding the emergence and evolution of a given gene [25]. The exon/intron organization of 116 M.bZIP genes was analyzed based on the number, distribution and splicing patterns of intron. In total, 23 of 116 M.bZIP genes contained no intron and 93 of 116 M.bZIP genes contained varied introns, the number from one to twenty-two ( Table  1). The members of M.bZIP genes in each group possessed similar gene structure in introns number and distribution. The number of introns in each group was uneven, with intronless genes only found in group

Group
Gene name Gene ID Chromosome location  IV and VI, dramatically varied from one to twenty-two introns in the group I. In the group II, most genes had less than six introns, and in the group V, most genes had more than ten introns (Table 1).

Number of introns
Further analyzing splicing patterns of DNA sequence encoding the conservative DNA-binding region (basic and hinge region) of bZIPs, five splicing patterns of intron were characterized and designated as pattern S187 I to S-V ( Figure 3 and Table S2). The pattern S-I was the most prevalent and lacked intron in this regions with 43 M.bZIPs from group IV and VI. The pattern S-II had one intron in phase 2 (splicing occurred between the second and the third nucleotide in one codon) in the basic region, while the pattern S-III had one intron in phase 0 (splicing occurred between codons) located between amino acid Gln and Ala. The members in pattern S-III was completely corresponded to that from group III and V except M.bZIP100. Pattern S-IV had two introns (each in phase 0), one located between amino acid Lys and Arg within the basic region and the other located between amino acid Gln and Ala within the hinge region, which member was exclusive in group I. The pattern S-V had only one intron in phase 1 (splicing occurred between the first and the second nucleotide in one codon) in the hinge region. Gene structure showed that each group of M.bZIPs mostly shared the same intron/exon structural pattern and remained conserved during the course of evolution of M.bZIP family genes, which were very consistent with those in other plants.

Conserved structural features of M.bZIP proteins
The bZIP domain is the core of the bZIP proteins, which preferentially binds to the promoter of their downstream target genes on specific cis206 elements. To analyze the characteristics of the bZIP domain of M.bZIPs, the amino acid sequence logo of bZIP domain was generated (Figure 4). bZIP domain from most M.bZIPs contained 39 amino acids and residues at position -27, -17, 1, 8 were much conserved as Asn, Arg/Lys, Leu and Leu. A consensus sequence was identified as Asn-X7-Arg/Lys-X9-Leu-X6-Leu started froom residue Asn (position-17 of bZIP domain) (Table S3), which was accord with that of bZIP members in other plants. Besides, some amino acids especially charged (Asp, Glu, Arg and Lys) were also conservative in this region (gray shadow in Figure 4). Besides the bZIP domain, bZIP proteins usually contain additional conserved motifs which might be related to the function of bZIP proteins as potential sites. A total of 25 motifs, including the conserved bZIP domain (be named motif 1) were identified in 116 M.bZIPs (Table S4) and their distribution was shown in Figure 5. Except motif 1, motif 7 was the most widely distributed in 69 M.bZIP members in each group. Some motif shared by several groups, motif 9 presented in four groups (group II, III, IV and V), motif 12 presented in three groups (group II, IV and VI), and motif 11 presented in two groups (group IV and V). Most of the conserved motifs, however, were presented in specific group, for example, motif 2, 3, 6, 14, 15 and 19 only existed in the group I, motif 4, 5, 13, 16, 20, 21, 23 and 24 only presented in the group II. Group-specific motifs may indicate their specific functions in each group. Protein phosphorylation, catalyzed by protein kinase is an important modification and directly related to protein function. Two conserved sequences, R/KxxS/T and S/TxxD/E, were verified as phosphorylated sites by a Ca 2+ independent protein kinase and casein kinase II, respectively [8,26,27]. Six motifs (motif 6, 8, 9, 12, 14 and 24) were verified harbouring phosphorylation sites of casein kinase II. Nine motifs (motif 2, 3, 4, 11, 15, 16, 19, 20 and 23) harbored function site of Ca 2+ independent protein kinase. Motif 3 and motif 18 contained both sites of above two kinases. Other functional sites were also identified in motifs. Motif 4 and motif 7 were extension of Leucine zipper region, function for the dimerization of M.bZIPs. Motif 6 only presented in the members of group I, was part of the DOG1 domain, which was related to the dormancy of the [28]. Motif 17 was found in six members of group V, with the presence of small stretches of Pro and aromatic amino acids Tyr, Phe and Trp, and is a part of the Pro-rich domain, similar motifs have been identified in Arabidopsis, which have been shown to have transcription activation potential [29].

Predicted dimerization of M.bZIP proteins
bZIPs bind to DNA as dimer when acting as transcription factors [30]. The leucine zipper region in bZIPs, arranged in the form of heptad repeats, is responsible for dimerization. Within each heptad, the amino acid positions are named as a, b, c, d, e, f and g in order in human and Arabidopsis [5,30]. Following this criteria, two to nine heptads were identified in different M.bZIPs (Table S5). Stability and specificity of leucine zipper dimmer are closely related to the amino acids present at the position a, d, e and g, since they are near the interface of leucine zipper [31]. Amino acids at the position a and d are typically hydrophobic, forming a hydrophobic core is essential for interaction between two monomers. The position e and g usually present charged amino acids, which mediate the specificity of dimerization as well as stability [8]. By analyzing all the 116 M.bZIPs, less than one tenth of charge amino acids (K, R, E and D) presented at position a, while other amino acids especially hydrophobic residues were predominant at the position a and d ( Figure 6A). About 22% of amino acids present at the position a were Asn, equivalent amount with that of bZIP from maize (22%). Further analyzing the composition of each heptad, the highest    Table S4. frequency of Asn at the position a occurred both in the second and fifth heptads, accounting for 59% and 58% respectively ( Figure 6B). The high frequency of Asn at the position a implied that it will form homo-dimerizing leucine zippers among the bZIP members, because Asn produce more stable N-N interactions than other amino acids [8]. Leucine present in position d was also import for dimerization stability. The frequency of Leu was high to 70% at the position d in M.bZIPs ( Figure 6A), a similar level with that of bZIPs from rice (71%) and greater than that from Arabidopsis (56%).
Different from the position a and d, the charged amino acids occupied near half in the position e and g, with frequencies of 45% and 60% respectively ( Figure 6A). The attractive and repulsive g and e pairs governs dimerization properties of bZIPs [30]. The amino acids present in each heptad of leucine zippers of M.bZIPs were characterized (Table S5) and the histogram of their frequency was presented in Figure  6C. Maximum frequency of the complete g and e' pairs appeared in the first heptad, showing basic repulsive in blue block. High frequency of attractive pairs was found in the second, fifth and ninth heptad, moreover, only -/+ attractive pairs were found in the ninth heptad.
Attractive pairs and presence of Asn in position a contribute to homo285 dimerization, while repulsive and incomplete pairs, and charged residues in position a may tend to hetero-dimerization [24,31]. According to these principles, total 116 M.bZIPs were classified into 29 sub-families (Table S5). One subfamily may tend to homo289 dimerization (subfamily BZ1 in Table S5), since attractive pairs presented in the first and second heptads and lacking any repulsive interactions. Two subfamilies could be considered having heterodimerizing specificity because they contained only repulsive interhelical interactions (subfamily BZ2 and BZ3). There was a subfamily (subfamily BZ29) was thought not to form dimer because these proteins lacked leucine zipper region. Other twenty-five subfamilies had both homoand hetero-dimerization properties.

Analysis of promoter sequences of M.bZIP genes
Transcriptional control of stress-responsive genes is a crucial part of the plant response to a range of abiotic and biotic stresses. Transcription factors have the potential to activate or repress genes through cis-acting sequences in promoter regions that respond to specific stresses (Y). Previous studies suggested that some bZIP proteins are involved in signaling and responses to abiotic/biotic stimuli [32]. To investigate the role of M.bZIP in stress response, the promoter sequences of seventeen randomly selected M.bZIP genes in different phylogenetic groups were analyzed.
Total 165 cis-acting elements were detected, which distribution and function annotation were shown in Table S7. Important cis-acting elements related stress-response including ABA-response element (ABRE), dehydration-response element (DRE), low-temperature responsive element (LTRE) were widely existed in the promoters of bZIP genes as shown in Table 2. These elements have been demonstrated to bind the relative proteins to modulate stress response in other plants [33][34][35][36]. Besides, some elements found in the promoters of M.bZIPs function for binding of important transcription factors, which factors were closely related to various abiotic stresses such as transcription factor MYB and MYC ( Table 2). High frequency of above elements in the promoters suggested that M.bZIP genes were involved in the response of multiple stresses in apple.

Expression of M.bZIP genes under multiple abiotic stress
To further examine whether the expressions of M.bZIP gene were induced by abiotic stresses, 20% PEG, 300 mM NaCl, 4°C and 100 μM exogenous ABA were treated to 1-year-old apple seedlings. The expression levels of seventeen genes, randomly selected in different groups, were checked by qRT-PCR. To ensure the results obtained are credible, two stress-related genes were selected as marker genes and their expressions were significantly up-regulated as the previous studies reported under stresses ( Figure S1). Checking the expression of M.bZIP genes selected, most of them were up-regulated at some extent in both leaves and roots under different stresses and exogenous ABA. Only several genes, their expressions did not change significantly or decreased at some extent. Under PEG treatment, all of seventeen M.bZIP genes were up-regulated more that 2-fold and nine genes were up-regulated more that 5-fold in roots ( Figure 7A). Most of them increased obviously in leaves ( Figure 7B). Under NaCl treatment, twelve of seventeen genes were up-regulated more than two folds in roots. Two of these (M.bZIP46 and M.bZIP87) changed most, reaching about 22-fold (at 6 h) and 20-fold (at 9 h), respectively Figure 7A. In leaves, four members (M.bZIP46, 49, 87 and 93) changed more than ten  folds respectively ( Figure 7B). In response to cold treatment, fourteen M.bZIP genes were up339 regulated in roots and thirteen of these were up-regulated more than two folds and seven were in the leaves ( Figure  7A). As to exogenous ABA treatment more than half genes selected showed up-regulated at some extent, with six genes increased more than 2-fold in roots and three in leaves (Figure 7). All the results above suggested that substantial amount of M.bZIP genes are sensitive to abiotic stresses and ABA induction.

Discussion
In recent years, gene family analysis has become an important measure to understand gene structure characteristic and evolution relationship which were used to predict the function of genes. The bZIP genes, encoding a big family of transcription factors, are ubiquitous present in eukaryote, from unicellular saccharomycetes [37] to vascular plants. The structure, evolution and function of bZIP genes have been thoroughly studied in Arabidopsis [7], rice [8], maize [31], cucumber [3], soybean [10], sorghum [9] and Brachypodium distachyon [24]. Most of the studies were focused on Gramineae and Cucurbitaceae, while little is known about the bZIP family in Rosaceae. In this study, total 116 genes encoded bZIP proteins were identified in the 'Golden Delicious' apple genome, which was more than that from other plants, such as 76 in Arabidopsis [7] and 89 in rice [8]. Multiplication of M.bZIP genes might be linked with the formation of sophisticated transcriptional control in apple. Nevertheless, the genome size of apple is 700 Mb, almost 7 times that of Arabidopsis (115 Mb) and 1.7 times that of rice (420 Mb). The difference of densities of bZIP genes among the three species may mean that in the evolutionary process bZIPs have different degrees of gene duplication and/or gene loss under the action of natural selection. Similar results were obtained by investigation of the divergence of the bZIP gene family in strawberry, peach and apple [13]. DNA-binding regions of bZIPs are relatively conservative. By analyzing DNA sequence encoding this region, different splicing patterns were found. Among 116 M.bZIP genes, 42 members have no intron in this region while other members with introns have four splicing patterns. In other plants, OsbZIPs contained seven splicing patterns and HvbZIPs contained six splicing patterns [8,38].
Comprehensive considering the splicing patterns of bZIP genes family of apple, rice and barley, we find that although there are some different of splicing phase and the interrupted amino acids, the positions of intron were highly conserved. Moreover, most of the intronless or intron-containing genes were clustered into the same group. Similar cases have also been observed in that of Arabidopsis [7] and rice [8]. The above results indicated that the positions of intron may be important for the evolution and divergence of the function in the M.bZIP family.
Protein structure of M.bZIPs is directly related to the function of transcription factors. bZIP domain is the core domain of M.bZIPs protein, which determines the characteristics of bZIP protein binding to DNA. However, the function diversity of the M.bZIPs protein may also be contributed by additional regions. Apart from bZIP domain, additional 24 motifs were identified in M.bZIPs. The majority of the M.bZIPs in the same group shared similar motif composition, suggesting that these conserved motifs play important roles in group-specific functions. However, among different groups, high divergence was found in their compositions of motifs. Duplication genes of M.bZIP possessed similar motif composition suggesting the ancestor genes with various motif composition have been maintained throughout evolution. More detailed research showed that seventeen motifs appeared to contain potential phosphorylation sites R/KxxS/T and S/TxxD/E). These phosphorylation sites were closely related to ABRE-dependent ABA signaling pathways [39]. Furthermore, motif 4 and motif 7 were found to be extension of leucine zipper regions, which suggested that these two kinds of motifs may function in dimerization of M.bZIPs. When function as transcription factors by binding to DNA, leucine zipper region of bZIPs should be dimerized to keep the stability of binding. Based on analysis of the leucine zipper dimerization specificity of bZIPs from H. sapiens [5], Arabidopsis [40], rice [8] and maize [31], the same structural rules of dimerization specificity of M.bZIP proteins were investigated. Base on the standard nomenclature for the seven unique amino acid positions (a, d, e and g) to arrange in the leucine zipper region, 116 M.bZIPs can be categorized into 29 subfamilies (BZ1-BZ29). Most of the subfamilies have both homo-and heterodimerization properties (BZ4-BZ28). Two subfamilies (BZ2 and BZ3) can be considered as having apparent heterodimerization specificity because they contain only repulsive interhelical interactions, and only one subfamily (BZ1) apparent homodimerization. All these go to prove diversity and complexity of the dimerization patterns in the leucine zipper of the M.bZIP protein family, with the potential to homodimerize with themselves or with members in the same subfamily as well as heterodimerize with other subfamily members. Moreover, the subfamily (BZ29) was thought not to form dimer, suggested that the five members of this subfamily may not have the activity of transcription factor.
Previous studies have shown that some definite bZIP genes play pivotal roles in developmental processes and responses of multiple stresses [41][42][43]. In this study, the analysis of 17 M.bZIP promoter sequences revealed that promoter regions of these genes contained variety of cis-acting elements associated with stresses, including ABRE, DRE, LTRE, MYB and MYC. To further study the response of M.bZIP genes to abiotic stress, qRT-PCR was used to determine the expression of M.bZIP genes under abiotic stresses. A range of expression patterns was observed and Most M.bZIPs were obviously up regulated under different treatments, including PEG, salt, cold and exogenous ABA (Figure 7). Detailed observations shown that M.bZIP genes were more sensitive to drought and salt stress, followed by cold stress. It is noteworthy that most M.bZIP genes response to stresses were more sensitive in root than that in leaves. For example, M.bZIP87 was obviously up-regulated in roots under PEG treatment while the expression level in leaves did not change obviously. The reason for this phenomenon may be explained that roots directly accepted stresses stimuli, and there were genes in the leaves may require a longer response time. Moreover, some M.bZIP genes have higher expression changes under different treatments both in roots and leaves. For example, M.bZIP46 has higher expression change under PEG, salt and cold treatment. This phenomenon suggested that these genes may be involved in variety pathways to respond to stresses. In conclusion, we performed a genome-wide analysis of the family of M.bZIPs in apple (Malus domestica) and proved that this family has the ability to respond to multiple stresses. Our investigation may lead to further understanding the relationship of function and structure of bZIP family members.