Genome-Wide Identification of bZIP Family Genes Involved in Drought and Heat Stresses in Strawberry (Fragaria vesca)

Basic leucine zipper (bZIP) genes are known to play a crucial role in response to various processes in plant as well as abiotic or biotic stress challenges. We have performed an identification and characterization of 50 bZIP genes across the woodland strawberry (Fragaria vesca) genome, which were divided into 10 clades according to the phylogenetic relationship of the strawberry bZIP proteins with those in Arabidopsis and rice. Five categories of intron patterns were observed within basic and hinge regions of the bZIP domains. Some additional conserved motifs have been found with the group specificity. Further, we predicted DNA-binding specificity of the basic and hinge regions as well as dimerization properties of leucine zipper regions, which was consistent with our phylogenetic clade and classified into 20 subfamilies. Across the different developmental stages of 15 organs and two types of fruits, the clade A bZIP members showed different tissue-specific expression patterns and the duplicated genes were differentially regulated, indicating a functional diversification coupled with the expansion of this gene family in strawberry. Under normal growth conditions, mrna11837 and mrna30280 of clade A showed very weak expression levels in organs and fruits, respectively; but higher expression was observed with different set of genes following drought and heat treatment, which may be caused by the separate response pathway between drought and heat treatments.


Introduction
The transcription factors containing bZIP domain [1,2] constitute one of the largest gene families in plants. The common feature, bZIP domain with length ranging from 60 to 80 amino acids (aa), includes two unique structures, a highly conserved DNA-binding basic and hinge region and a relatively diversified leucine zipper region [3]. The basic and hinge region consists of 18 aa residues with an invariant motif (N-X 7 -R/K-X 9 ). In contrast, the downstream leucine zipper region is made up of heptad leucine repeats, a repeated pattern of several amino acids, or other hydrophobic amino acids, such as Ile, Phe, and Met [4]. The leucine zipper region is located exactly 9 amino acids downstream from the C-terminal of the basic region. It has been proved that the leucine zipper region is more favorable to constitute an amphipathic-helix with two helical turns within each heptad, because of its special amino acid composition [1,2].
Although the bZIP transcription factor family in strawberry (F. vesca), apple, and peach has been identified [33], only information on the evolutionary characterization of bZIP genes in Rosaceae was provided. Limited evidence is available regarding basic annotation of structural and functional features in the bZIP family among the economically important strawberry. The cultivated strawberry (F. × ananassa) is a model plant for nonclimacteric fruit species. F. × ananassa derived from 4 diploid ancestors has a relatively complex octaploid genome harboring 56 chromosomes (2n = 8x = 56). In contrast, our experiment system F. vesca, the diploid woodland strawberry, has been sequenced, which offers substantial advantages for molecular and physiological studies [34]. In this research, we performed comprehensive studies on expansion pattern, structural annotation, DNA-binding site specificity, dimerization properties, and expression levels in different stages and abiotic treatments to shed light on bZIP genes in F. vesca.

Identification of bZIP Genes in Woodland Strawberry.
To generate a comprehensive list of F. vesca bZIP genes, F. vesca v1.1 genome was downloaded from Phytozome (https:// phytozome.jgi.doe.gov/pz/portal.html). Firstly, a HMM profile of the bZIP domain (PF00170) was downloaded from Pfam (http://pfam.xfam.org/). Secondly, HMMER analysis was performed to search all F. vesca protein sequences and extracted the corresponding coding DNA sequence (CDS) and amino acid sequences. Thirdly, all of the putative bZIP domains were manually verified (E value < 1.0) within Pfam and SMART (http://smart.embl-heidelberg.de/) to confirm their presence and completeness of bZIP domain. Molecular weight, instability index, isoelectric points, aliphatic index, and grand average of hydropathicity (GRAVY) of candidate FvbZIP proteins were detected in ExPASy ProtParam server (http://web.expasy.org/protparam/).

Phylogenetic Analysis and Structures of FvbZIP Gene.
Phylogenetic analysis was used for classifying FvbZIP genes into ten clades [33]. The amino acid sequences of bZIP domains located in bZIP proteins from strawberry, A. thaliana (AtbZIP) [8], and rice (OsbZIP) [2] were used to generate a phylogenetic tree through ClustalW alignment and the unrooted Neighbor-joining (N-J) method using MEGA 6.0 [35]. N-J analysis with pairwise deletion and the Jones-Taylor-Thornton (JTT) model was performed using the 1000 replicate bootstraps to support for each node in phylogenetic tree. The exon-intron organizations with intron phase in FvbZIPs were performed by alignment within CDS sequence and corresponding genomic sequence and presented by GSDS 2.0 (http://gsds.cbi.pku.edu.cn/).

Conserved Motif Identification of FvbZIP Proteins.
The motif distribution of FvbZIP proteins were extracted by MEME (http://meme.nbcr.net/meme/), and parameters are set as follows: the largest number of discovered and conserved motifs set 12 and motif width was from 30 to 50 aa. The remaining parameter settings were kept default.

Dimerization Properties of FvbZIP Proteins.
To describe the foundation of the dimerization stability and speculate the dimerization specificity of 50 FvbZIP members, we named special amino acid residue positions as g, a, b, c, d, e, and f to put in order of each heptad on the basis of standard nomenclature for the seven heptad repeats [36] (see Supplementary Figure S4 available online at https://doi.org/10.1155/2017/3981031). The first heptad was manually arranged from four amino acid residues (corresponding amino acid located at g position in the first heptad) before appearance of the first leucine in the bZIP domain. Both boundaries of N-terminal and C-terminal within the FvbZIP leucine zipper region have been the same as applied to bZIP proteins in Arabidopsis [37]. 2.6. Quantitative RT-PCR Analysis. Transcriptional changes of FvbZIP genes responding to drought and heat treatments were determined by quantitative RT-PCR analysis on Stratagene Mx3000P Real-Time PCR system using SYBR® Premix Ex Taq™ (TaKaRa, Japan). The PCR amplification conditions were as follows: 95°C 10 min, 40 cycles of 15 sec at 95°C, and 1 min at 60°C; at the end, the melting curve analysis was executed for verifying the specificity of the primer with the following stage: 95°C for 15 sec, 60°C for 1 min, and 95°C for 15 sec. The relative gene expression was determined using an interspacer 18S-26S strawberry gene as an endogenous control gene [38]. The relative expression levels of the target genes were calculated by the 2 −ΔΔCt method [39]. Each specific primer pair was listed in Supplementary Table S1. Three independent biological replications were performed at each treatment point.

Plant
RPKM (reads per kilobase of exon per million mapped reads) of FvbZIP genes in different tissues and early stage fruit development were directly downloaded from available database (http://bioinformatics.towson.edu/strawberry/) (Table S5), which reflects mRNA levels. The expression levels (RPKM) for FvbZIP genes in yellow and red fruits of strawberry were from Tianbao Yang's lab in USDA-ARS. To compare the differential expression of bZIP genes (Table S5), the expression level of each gene in the different stages was read and normalized using multiexperiment viewer (MeV) software [40].

Results and Discussion
3.1. Identification and Characterization of Strawberry bZIP Family. A total of 50 FvbZIP genes were identified in strawberry genome in Phytozome. The characteristic parameters of all predicted FvbZIP proteins are listed in Table 1, including chromosome location, protein length, molecular weight, theoretical pI, instability index, aliphatic index, and grand average of hydropathicity (GRAVY). The length of amino acids of FvbZIPs is from 114 (mrna26148) to 1106 (mrna03778) aa with an average of 402.54 aa. Compared to those in Arabidopsis [8] (harboring 75 members with an average 321 aa in length) and rice [2] (harboring 89 members with an average 311 aa in length), strawberry has 50 bZIP members with longer average length than these two plants.
After divergence of monocots and eudicots, lower frequency of bZIP candidates' evolution occurred in eudicots than monocots [9], which may explain the fact that strawberry and Arabidopsis have fewer bZIP members than rice. The predicted molecular weights of FvbZIP proteins range from 13.51 (mrna26148) to 122.71 (mrna03778) kDa, with pI from 4.73 (mrna29159) to 9.56 (mrna11666) ( Table 1).

Expansion Pattern of bZIP Family in Strawberry,
Arabidopsis, and Rice. To survey the extent of species-specific expansion of the bZIP genes in strawberry, Arabidopsis, and rice, we performed a joint phylogenetic analysis of strawberry, Arabidopsis, and rice bZIP proteins. As shown in Arabidopsis and rice analysis [2,8], the joint phylogenetic tree of three species grouped all bZIPs in ten distinct clades (A, B, C, D, E, F, G, H, I, and S) (Supplemental Figure S1A). Six bZIP genes in strawberry, Arabidopsis, and rice could not be conclusively mapped to any clade and was renamed U. We investigated some specific nodes which led to strawberry-, Arabidopsis-, and rice-specific classes within clades (red circle on node in Supplemental Figure S1A). As these nodes located in specific classes among strawberry, Arabidopsis, and rice represent the divergence point in the evolution and indicate the most recent common ancestors (MRCAs) prior to division, we found that some bZIP genes possibly have been present in the MRCAs of strawberry, Arabidopsis, and rice, but one or two species-specific bZIPs were lost in strawberry, Arabidopsis, or rice. Six classes contained only strawberry bZIP genes (Supplemental Figure S1A, red arrows), fifteen classes contained only Arabidopsis bZIP genes (Supplemental Figure S1A, green arrows), and thirty classes only contained rice bZIP genes (Supplemental Figure S1A, yellow arrows), suggesting that gene loss event has occurred in those classes in other species. The interspecies classes revealed that a parallel evolution of bZIP genes happened in three plant species and the orthologous bZIP proteins were observed to play a similar role [41]. The number of classes illustrated a fact that at least 85 ancestral bZIP genes existed before the strawberry, Arabidopsis, and rice split. After splitting from the MRCA, strawberry has lost the largest number of bZIP genes, followed by Arabidopsis and rice (Supplemental Figure S1B). A wide range interspecies comparison of strawberry, Arabidopsis, and rice revealed the species specific and nonspecific proteins.
3.3. Gene Structure of FvbZIP Genes. In order to observe the structure of the FvbZIP genes, we performed an analysis of the number and distribution of exon-intron, which are recorded to play key roles in gene family evolution and further promote to comprehend the emergence and evolution of a given gene [42,43]. As shown in Figure 1, 10 (20%) of total FvbZIP genes exhibit intronless, which occurs exclusively in group S, accounting for 88.9%. About 20% of peach, 19.1% of rice, 15.3% of maize, and 20.4% of castor bean bZIP genes were predicted to be intronless [2,4,10,44]. Among FvbZIP genes containing introns, the number of introns in the open reading frame (ORF) varied from one (mrna29159) to 20 (mrna03778), indicating a great diversity in the FvbZIP members ( Figure 1). According to number, position, and splicing phase of the introns located in FvbZIP genes, the considerably conserved basic and hinge regions within the bZIP domain were divided into five patterns (a, b, c, d, and e) ( Figure 2; Supplemental Figure S2). A previous study has showed that patterns of intron positions and splicing phases located at the bZIP domain regions were widely regarded as a promoter to understand homology and evolution of bZIP genes [2]. Based on intron location relative to the codon, introns can be classified into three patterns. Introns that do not interrupt the codons, between the first and second bases of codon and between the second and third bases of the codon, are termed phase 0 (P0), phase 1 (P1), and phase 2 (P2), respectively [45].
Patterns a (including 14 FvbZIPs) and e (including 14 FvbZIPs) were the most widespread ( Figure 2). Pattern a had an intron in P0 within the hinge region behind the position −6 (Q). Both patterns b and c had one intron each followed the position −22 (R) located at the basic region, but the phase was different. Pattern b had an intron in P0, whereas pattern c had an intron in P2. Pattern d with two introns (both located at P0), one followed position −26 (K) of the basic region and the other one located at the hinge region at position −5, was interrupted by Lys and Ala. Pattern e had no intron distribution neither in the basic nor hinge region. Ten (71.2%) FvbZIPs within pattern g were intronless, while the other 4 had introns located at the outside of their basic and hinge regions (Figures 1 and 2). In particular, the position and phase were observed to be highly conserved when there was an intron interrupting in the hinge region. However, positions of introns located in the basic region were variably distributed in P0 or P2 ( Figure 2; Supplemental Figure S2). The similar intron distribution and splicing patterns were observed in the bZIP domain regions of rice, castor bean, and maize, which declares that the intron pattern within the basic and hinge regions of bZIP genes were highly conserved in plants [2, 4, 10]. Though the lengths of intron were variable (Figure 1), most members clustered in the same clade classified above (Supplemental Figure S1) exhibited the same or similar intron number and splicing phase. For instance, same splicing phases and gene structures were observed in four FvbZIP members within clade A and eight members within clade S; five members in clade I also showed a similar splicing phase and gene structure. The intron pattern a was shared in the whole members in clade A, clade S consisted of the most intronless members showed the intron pattern e, members in clade C showed intron pattern c, and clade D (except for mrna31321) shared pattern d. In conclusion, the conserved splicing phases and gene structures in each clade were in agreement with the classification of FvbZIPs in phylogenetic analysis.

Conserved Motifs of FvbZIP Proteins.
Besides the bZIP domain, some of the other multiple conserved motifs had been detected in Arabidopsis, rice, maize, caster bean, and peach [2,4,8,10,44]. As shown in Figure 3, a total of 12 conserved motifs were identified in FvbZIP genes. The consensus sequence and the amino acid width of those conserved motifs are given in Supplemental Table S2. It can be observed that some motifs are shared by several clades, such as motif 9 which was located in clades C and S, motif 4 shared by clades E and I, and motif 12 which was present in clades A and D. Furthermore, majority of conserved motifs appear in a specific clade, for instance, motifs 2, 3, 5, 6, and 11 in clade D; motifs 7 and 8 in clade A; and motif 10 in clade G ( Figure 3). It is a proved hypothesis that the group-specific conserved motifs promote to decide specific functions of FvbZIP genes in each clade and foundation for the functional divergences between members from different clades [46][47][48]. The conserved motif analysis was further proved to be consistent with phylogenetic relationship and classification of FvbZIP genes.
Notably, motifs 2, 9, and 10 were identified as DOG1 (delay of germination), bZIP C, and MFMR (multifunctional mosaic region) according to Pfam databases, respectively. The previous study has found DOG1-corresponding gene to be related to seed dormancy in a form of quantitative trait locus [49]. It has been suggested that N-terminal half of MFMR motif is rather rich in proline residues and has been termed the PRD [50], and some of these motifs may play a key role in regulating protein-protein interactions [51]. Motif 10 with a Pro-rich domain was supposed to be involved in mediating protein-protein interactions [51]. Motifs [54]. We can speculate a similar function of bZIP proteins due to conserved motifs shared by Arabidopsis, rice, castor bean, and maize, while the roles of other conserved motifs found in FvbZIP proteins are not yet clear.

DNA-Binding Site Specificity of FvbZIP
Proteins. DNAbinding specificity owed by bZIP proteins is determined by certain key amino acids located at the most conserved basic and hinge region of the bZIP domain, and these two regions have a direct interaction with DNA cis-elements [55,56]. The annotation and classification of amino acids located at the basic and hinge region were characterized by diagnostic sequence structure in Arabidopsis and rice [2,8], which indicates that some highly conserved amino acids are present in each clade (Supplementary Figure S3). The first leucine of the leucine zipper region was designed as +1, and C-terminal amino acid of the hinge region was numbered −1 [55]. It is proved that some functional replacements on invariant sites (−18 and −10) by other amino acids will lead to new DNAbinding specificities [57]. In the strawberry bZIP family, these replacements infrequently happened and just occurred in one group. mrna30252 in clade U with a hydrophobic Ile (Isoleucine) residue at position −10 replaced arginine/lysine, indicating that perhaps they are not able to bind DNA or might hold a unique DNA-binding specificity. Furthermore, we can speculate the structural characteristics and DNAbinding specificity of FvbZIP proteins in a clade model, as listed in Supplementary Table S3. Expression level of a particular gene is sometimes determined by binding affinity to specific sites of promoters of genes, which is resulted of variations in their target sites in promoters of different genes [2]. Computational discovery and further experimental verification of DNA binding ability would promote to identify which genes are selectively activated by various FvbZIP genes.
3.6. Dimerization Properties of FvbZIP Proteins. The function of bZIP proteins requires dimerization between parallel coiled-coil constructions [58]. Amino acids present at positions a, d, e, and g near the leucine zipper interface play an important role in regulating oligomerization of leucine zipper domain as well as specificity and stability of dimerization [4]. Thus, we performed detailed characterization of the amino acids located at positions a, d, e, and g to predict their  The numbers "0" and "2" denote different splicing phases, "0" means splicing occurred after the third nucleotide of the codon, and "2" means splicing occurred after the second nucleotide.
functions. The length of the FvbZIP leucine zipper is variable to range from two to nine heptads. About 18% of amino acids located at position a were Asns (asparagine) in 41 FvbZIP proteins (Figure 4(a)), which is similar to both 16% in humans and H. sapiens, but lesser than observed in Arabidopsis (22%) [37], rice (23%) [2], maize (22%) [4], and castor bean (26%) [10]. The highest frequency of Asn at position a was displayed in the second heptad (accounting for 44%) followed by the fifth heptad (accounting for 36%) (Figure 4(b)), which is also observed earlier for rice and Arabidopsis [2,37]. Since Asns can form a more stable N-N interaction at the a ↔ a′ position than the other a position amino acids, high frequency of Asn at the a position promote to create homodimerizing Leu zippers among the bZIP family [59]. This result raises the possibility that a majority of FvbZIPs prefer to homodimerization. It is observed that minority of charged amino acids (R, K, E, and D) displayed in position a, which may contribute to form hetero-dimerization.
Amino acids located at positions a and d are generally hydrophobic residues and pack on the surface of "knobs and holes" pattern to establish a hydrophobic core area that is important to homo-/hetero-dimerization [60]. The frequency of the conserved Leu responsible for dimer stability at d position, was found in 68% of 50 FvbZIP proteins. The abundance of other hydrophobic amino acid (including I, V, and M) at d position was comparable in FvbZIPs (10%, Figure 4(a)), Arabidopsis (19%), and rice (17%). The Leu at position d is a type of highly conserved and stable aliphatic amino acid residue [60], which is critical for keeping the dimerization stability. Compared with the AtbZIP proteins (the proportion of Leu at d position accounts for 56%) [37] and RcbZIP proteins (56%) [10], FvbZIPs showed more abundant Leu located at d position, which suggests that strawberry genes probably hold a stronger stability dimerization on Leu zipper than Arabidopsis.
It is observed that the abundance of charged amino acids displayed at e and g positions, such as acidic amino acids (including E and D) and basic amino acids (including R and K), accounts for 50% and 50%, respectively (Figure 4(a)). Compared to Arabidopsis (charged amino acids at e and g positions were, respectively, accounting for 41% and 53%), the frequencies of charged amino acids located at e position was higher but was lower in strawberry at g position. The charged amino acids located at positions g and e were replaced by alanine can result in creating tetramers instead of dimmers [61]. A greater number of charged amino acids (containing E, D, R, and K) located at e and g positions were found to flank the dimerization interface, which are thought to interact electrostatically to form salt bridges mediating the dimerization and determining the dimerization specificity as well as stability [62].
As dimerization properties of FvbZIP proteins were governed by different contributions of charged amino acids located at e and g positions, we calculated the presence of attractive and repulsive g ↔ e ′ pairs in each heptad of strawberry Leu zippers (Supplementary Figure S4) and displayed corresponding frequency histogram in Figure 4(c). Four groups (basic repulsive, acidic repulsive, +/− attractive, and −/+ attractive) were found in complete g ↔ e′pairs according to both the g and following e position amino acids are charged, while incomplete g ↔ e ′ pairs were defined as only one charged position in g or e [10]. The highest frequency of the complete g ↔ e ′ pairs was observed to locate at the first heptad for 46%, majority of which was found to be the basic repulsive g ↔ e ′ pair. In following three heptads, the frequencies of the complete g ↔ e ′ pairs showed a sharp decline. The frequency of complete g ↔ e ′ pairs appeared to be increased in the fifth heptad, especially on the −/+ attractive gene pairs. Moreover, attractive g ↔ e′ (including +/− attractive and −/+ attractive) pairs were found most frequently within the second, fifth, and ninth heptad. It is observed that the attractive g ↔ e ′ pairs had no tend to form homo-or hetero-dimerization, while the repulsive g ↔ e ′ pairs preferred to make a conformation hetero- dimerization [62]. Furthermore, biophysical measurements also reveal that repulsive g ↔ e′ pairs present a more primary role in driving dimerization specificity than attractive g ↔ e ′ pairs [4]. According to the observation of the dimerization properties in our study, the 50 FvbZIP proteins could be classified into 20 types (Supplementary Figure S4 and Table S4).

Expression of FvbZIP Genes in Different
Organs. Many bZIP proteins have been shown to be involved in the processes of plant growth and development [63]. To explore the tissue-specific expression patterns of genomic scale FvbZIP genes, the transcriptome data were analyzed to determine those genes expressed in different organs and tissues ( Figure 5). FvbZIP genes were showed to rank from highest (mrna21832 in style and yellow fruit) to lowest (mrna31321 in embryo) according to their differential expression across main tissues. The FvbZIP genes clustered in the same clade showed different tissue-specific expression patterns. mrna15193, mrna21832, mrna14942, and mrna02284 within clade S showed particularly high expression in most organisms, while the others of clade S presented a relatively low expression in all tissues except for mrna18282 in ghost (refers to the entire seed with its embryo removed) and ovule. The duplicated genes in clades (predicted in the Plant Genome Duplication Database) showed different transcript abundance in specific tissues, such as mrna32022 and mrna18928 in flowered and perianth ( Figure 5), which is also supported by the duplicated bZIP genes in tomato [64]. The divergences in expression profiles between paralogs revealed that some of them may acquire new functions after duplication in the evolutionary process.
As an activator of downstream gene expression, members of clade A bZIP genes share function as important regulatory   roles in the ABA signaling pathway, abiotic stresses, and tissues development in Arabidopsis [8]. MdbZIP26 was abundantly expressed in roots, stems, and apical buds [65]. OsABF2 showed expression in different tissues of rice and induction of drought, salinity, cold, and ABA. In fact, clade A bZIPs often have played key roles not only in seed development but also in fruit maturation [66]. We further analyzed the expression profiles of the 8 clade A FvbZIP genes in various organs and fruits using transcriptome data.
Firstly, the members show different tissue-specific expression. For example, mrna00393 and mrna09110 are moderately expressed in most tissues, while mrna30280 is more enriched only in flower ( Figure 5(a)). Secondly, the respective members of the duplicated FvbZIPs pair (mrna14556 and mrna08566, predicted in the Plant Genome Duplication Database) have different pattern specificities. The duplicated pair has a similar expression pattern in each tissue but a different pattern in red and yellow fruits (Figures 5(a) and 5(b)). 39    It is evident that although highly conserved on amino acids sequence, those duplicated genes are differently regulated [67]. Thirdly, the same gene show different expression profiles in different organs. mrna11837 is expressed at an extremely high level in seedling but very low in majority of other tissues (Figure 5(a)). Overall, the expression of clade A FvbZIP genes has different organ specificities, including a functional diversification coupled with the expansion of this gene family in strawberry. At least one clade A bZIP gene (mrna09110) is expressed in most investigated strawberry organs. It is consistent with previous studies in Arabidopsis and apple [8,65].

Expression Patterns of FvbZIP Genes in Leaves under
Drought and Heat Treatments. We next performed an expression investigation of the clade A FvbZIP genes following exposure to drought and heat. A range of expression patterns was observed and some of the genes clearly upregulated, while others were downregulated (Figures 6(a) and 6(b)). As shown in Figure 6(a), for example, mran00393, mrna08566, mrna30280, and mrna11837 showed a significant upregulation after the plants were treated by drought for 2 days, while all genes in clade A were significantly downregulated at day 4 ( Figure 6(a)). These four genes showed a higher expression level than the previous point after rewatering (Figure 6(a)). When exposed to drought, the expression of mrna14556 and mrna28250 decreased gradually (Figure 6(a)). Heat treatment caused mrna30280 transcription levels to increase in all points, and the mrna30280 transcription level peaked 12 h after the heat treatment was started (Figure 6(a)). Similarly, mrna00393 showed an obvious increase in expression at 12 h (Figure 6(a)). mrna14556 and mrna28250 kept low expression during heat treatment (Figure 6(a)). It is worth noting that the expression of mrna11837 in most organs in Figure 5(a) and mrna30280 in fruit in Figure 5(b) were weak under normal growth conditions but both higher following drought and heat. We conclude that these two genes (mrna11837 and mrna30280) have a different pathway to respond to drought and heat treatments. Drought and heat responses in plants may be a separate pathway, which is consistent with the finding in grape [11]. However, mrna14556 expression was not detected in any organ, even in response to drought or heat. A similar observation was made in grape, apple, and sorghum with some bZIP genes [65,68,69], and we speculate that those genes may be degenerated into a nonfunctional gene, or it is involved in other processes that were not investigated here, such as other abiotic and biotic stresses [68,69], or they may be expressed at other organs and time points after treatment [65]. This can be supported by the fact that the bZIP ABF genes, ABF1 and ABF4, are induced by cold as well as ABF2 and ABF3 expression are induced by high salt in vegetative tissues [26,53]. The constitutive overexpression of ABF3 in rice and Arabidopsis consequently made the transgenic plants more tolerant to stress conditions [70,71]. Moreover, OsbZIP23 and OsbZIP72 play an important role in training drought tolerant rice and resistance to drought treatment by activating ABA signaling [21,22]. The genes upregulated by stresses was suppressed in mutants of OsbABF1 under ABA treatment, which indicates that OsbABF1 is related to abiotic stress responses through ABA pathway in rice. We concluded that the clade A FvbZIP genes showed a range of response patterns when exposed to drought and heat stresses, revealing that they may respond to stress signaling, which is also observed in Arabidopsis, apple, and grape [8,65,68].

Conclusion
We performed genome-wide analysis of the strawberry bZIP transcription factor family and conducted a detailed investigation of their evolutionary relationship, structural feature, and organization. In addition, we characterized the expression of some clade A genes following drought and heat treatment. These results may prove useful in developing strategies for the further improvement of stress tolerance in strawberry.

Conflicts of Interest
The authors declare no conflict of interest.