Genome-wide identification and evolutionary analyses of bZIP transcription factors in wheat and its relatives and expression profiles of anther development related TabZIP genes

Among the largest and most diverse transcription factor families in plants, basic leucine zipper (bZIP) family participate in regulating various processes, including floral induction and development, stress and hormone signaling, photomorphogenesis, seed maturation and germination, and pathogen defense. Although common wheat (Triticum aestivum L.) is one of the most widely cultivated and consumed food crops in the world, there is no comprehensive analysis of bZIPs in wheat, especially those involved in anther development. Previous studies have demonstrated wheat, T. urartu, Ae. tauschii, barley and Brachypodium are evolutionarily close in Gramineae family, however, the real evolutionary relationship still remains mysterious. In this study, 187 bZIP family genes were comprehensively identified from current wheat genome. 98, 96 and 107 members of bZIP family were also identified from the genomes of T.urartu, Ae.tauschii and barley, respectively. Orthology analyses suggested 69.4 % of TubZIPs were orthologous to 68.8 % of AetbZIPs and wheat had many more in-paralogs in the bZIP family than its relatives. It was deduced wheat had a closer phylogenetic relationship with barley and Brachypodium than T.urartu and Ae.tauschii. bZIP proteins in wheat, T.urartu and Ae.tauschii were divided into 14 subgroups based on phylogenetic analyses. Using Affymetrix microarray data, 48 differentially expressed TabZIP genes were identified to be related to anther development from comparison between the male sterility line and the restorer line. Genes with close evolutionary relationship tended to share similar gene structures. 15 of 23 selected TabZIP genes contained LTR elements in their promoter regions. Expression of 21 among these 23 TabZIP genes were obviously responsive to low temperature. These 23 TabZIP genes all exhibited distinct tissue-specific expression pattern. Among them, 11 TabZIP genes were predominantly expressed in anther and most of them showed over-dominance expression mode in the cross combination TY806 × BS366. The genome-wide identification provided an overall insight of bZIP gene family in wheat and its relatives. The evolutionary relationship of wheat and its relatives was proposed based on orthology analyses. Microarray and expression analyses suggested the potential involvement of bZIP genes in anther development and facilitated selection of anther development related gene for further functional characterization.


(Continued from previous page)
Conclusions: The genome-wide identification provided an overall insight of bZIP gene family in wheat and its relatives. The evolutionary relationship of wheat and its relatives was proposed based on orthology analyses. Microarray and expression analyses suggested the potential involvement of bZIP genes in anther development and facilitated selection of anther development related gene for further functional characterization.
Keywords: Basic leucine zipper (bZIP), Common wheat, Orthology analyses, Gene expression, Anther development, Over-dominance Background Transcription factors play vital roles in almost all plant biological processes [1]. Among the largest and most diverse dimerizing transcription factor families, the basic leucine zipper (bZIP) family of proteins in plants participates in regulating various processes including floral induction and development, stress and hormone signaling, photomorphogenesis, seed maturation and germination, and pathogen defense [1,2].
The bZIP proteins characteristically possess a bZIP domain composed of a basic region and a leucine zipper [3]. The highly conserved basic region consists of approximately 16 amino acid residues that contain an invariant N-x7-R/K motif and that are responsible for sequencespecific DNA binding [1,2]. The less conserved leucine zipper is composed of a heptad repeat of Leu or other bulky hydrophobic amino acids (i.e., Ile, Val, Phe or Met) positioned exactly nine amino acids toward the Cterminus to form an amphipathic helix, which confers homo-or hetero-dimerization specificity [2]. Plant bZIP proteins preferentially bind to DNA sequences contain an ACGT core, especially like G-box (CACGTG), C-box (GACGTC) and A-box (TACGTA) [4,5].
Members of the bZIP transcription factor family have been comprehensively identified or predicted in several plant species [1,2,[6][7][8][9][10][11][12]. In the model plants Arabidopsis and rice, 76 and 94 distinct members of the bZIP family have been identified and analyzed, respectively [1,2]. bZIP family members of Arabidopsis have been classified into ten groups based on a similar basic region and additional conserved motifs outside the bZIP domain [1]. OsbZIP proteins are divided into 11 groups based on their putative DNA-binding specificity and dimerization properties [2].
Until now, some bZIP transcription factors have been reported to function in the anther development. In Arabidopsis, AtbZIP1 overexpression is detrimental to pollen development and AtbZIP1 is thought to be closely related to water movement during anther development, when anthers first absorb water for growth but dehydrate before dehiscence [13][14][15]. AtbZIP34 controls pollen wall patterning and several metabolic pathways in developing pollen [16,17]. Among the ten members of TGACG (TGA) motif-binding proteins composing a distinct subgroup in the bZIP family of Arabidopsis, TGA 9 and TGA10 have been reported to play a role in anther patterning and dehiscence. TGA9/10, together with ROXY1/2, positively regulate a common set of genes that contribute to the development of the tapetal, which provides nutrients and secretes enzymes and structural components of the pollen coat [1,18]. In rice, OsABI5 is highly expressed in mature pollen and suppression of OsABI5 expression in transgenic rice lines causes abnormal mature pollen development and low fertility [19,20]. However, among the largest and most diverse transcription factor families, bZIP gene family have not been systematically identified from the wheat genome. And, to date, few wheat bZIP transcription factors have been identified to be involved in anther development.
The hexaploid bread wheat (Triticum aestivum; 2n = 6x = 42; AABBDD) derived from two hybridizations in the evolution history between three gramineous ancestors. Wheat arose as a result of hybridization between the cultivated tetraploid emmer wheat (T.dicoccoides, AABB) and diploid goat grass (Ae.tauschii, DD) approximately 8,000 years ago. The three sub-genomes of wheat evolved from T.urartu (A sub-genome), from a species that might be from the section Sitopsis (B sub-genome), and from Ae.tauschii (D sub-genome), respectively [21]. Wheat is one of the most widely cultivated and consumed food crops in the world [22,23]. Recently, an ordered draft sequence of the 17Gb hexaploid bread wheat genome has been generated by sequencing isolated chromosome arms and major efforts have been undertaken worldwide to sequence and annotate the wheat genome [21][22][23][24][25][26][27][28].
Wheat thermosensitive genic male sterile (TGMS) lines such as BS366 are of great importance for the utilization of heterosis in wheat breeding [29]. When exposed to low temperature conditions, BS366 becomes sterile, thus allowing the large-scale production of F1 hybrids via crossing with wheat restorer lines. The male sterility phenotype is heritable in wheat TGMS lines and controlled strictly by temperature [29,30]. Previous studies have demonstrated that male fertility in BS366 is controlled by temperature during the period from the pollen mother cell (PMC) stage to the meiosis stage [31].
The two-line cross combinations of BS-series are our exclusive plant materials for wheat breeding and several of them have been widely used in wheat breeding and production. As important parental materials used in twoline combined breeding, BS366 (male sterility line) and TY806 (restorer line) generate a fertility-restored F1 hybrid. Thus, the cross combination of BS366, TY806 and their hybrid F1 is the excellent candidate for studying anther development in wheat.
In this study, we report the identification of bZIP family genes in the current genomes of wheat and its relatives, including Triticum urartu, Aegilops tauschii and Hordeum vulgare. Then, we have identified the orthologs and in-paralogs between each pair of wheat and its relatives. Further, we have proposed the evolutionary relationship among wheat and its relatives based on the orthology analyses of bZIP family. bZIP transcription factors in wheat, T.urartu and Ae.tauschii have been classified on the basis of phylogenetic analyses. 23 anther development related bZIP genes are selected from differentially expressed genes by microarray data analyses. Their gene structure and cis-acting elements have been also analyzed. Expression analyses have been performed to reveal their expression profile in response to low temperature and tissue-specific expression pattern in wheat male sterility line BS366, restorer line TY806 and their F1 hybrid and an emphasis is laid on expression in anther. Our results deepen the understanding of evolutionary relationship of wheat and its relatives and provide a prospective for biological involvement of bZIP family genes in anther development and heterosis of wheat.

Results
Identification of bZIP transcription factors in wheat, T. urartu, Ae. tauschii and barley In the Gramineae family, except for T.urartu and Ae.tauschii, barley and Brachypodium are also close relatives of wheat. The evolutionary relationships between rice, Brachypodium, barley and wheat have been speculated based on the mean synonymous substitution rates (Ks) of orthologous gene pairs [32]. To identify bZIP transcription factors in wheat, T. urartu, Ae. tauschii, barley and Brachypodium, a hidden Markov model (HMM) search was performed using the HMM profiles of the bZIP domain (Pfam accession No.: PF00170 and PF07716) as queries against the respective protein databases. As a result, the amino acid sequences of 100 possible bZIP transcription factors were obtained using a Perl program from HMM search hits against the T. urartu and Ae. tauschii protein databases, respectively. Among these sequences, three contig sequences were removed from the analysis. Subsequently, two T. urartu proteins and one Ae. tauschii protein were discarded because they lacked the typical bZIP domain, finally resulting in the identification of 98 T. urartu and 96 Ae. tauschii bZIP transcription factors. Similarly, 182 nonredundant bZIP proteins were identified in the wheat genome and named. All these wheat bZIP genes were designated as TabZIP genes and given a number designation from 1 to 182 as an unique identifier as proposed for bZIP transcription factors in Arabidopsis [1]. The nomenclature was based on the positions of these genes on the wheat chromosomes, from the short arm to the long arm and in the order of 1A to 7A, 1B to 7B and 1D to 7D. Additionally, several TabZIP genes had been identified in previous studies [33][34][35][36][37][38][39]. TaABF1 [33] and Tab-ZIP60 [34] were identical to TabZIP23 and TabZIP51 in our 182 wheat bZIP family members, respectively. Wlip19a, Wlip19b and Wlip19d (three homoeologous loci of WLIP19 on chromosomes 1A, 1B, and 1D) [35] were identical to TabZIP2, TabZIP59 and TabZIP123, respectively. TaOBF1b and TaOBF1d [35] were identical to TabZIP98 and TabZIP149. However, TaOBF1a [35], WABI5 [36], TaABP1 [37], TaABI5 [38] and TaABL1 [39] was not within 182 bZIP family members identified by us. Finally, these 5 TabZIP genes were added to wheat bZIP gene family. So a total of 187 wheat bZIP genes were identified from the current wheat genome. Similarly, the 98 T. urartu and 96 Ae. tauschii bZIP genes were designated as TubZIP1-98 and AetbZIP1-96, respectively. Furthermore, 107 bZIP proteins were also seperated from the barley genome. The Brachypodium genome had been repoted to contain 96 bZIP genes [12]. The nomenclature and chromosomal location of all bZIP genes are shown in Additional file 1: Table S1.
Orthology and in-paralogy of T. urartu and Ae. tauschii bZIP proteins Orthologous genes in different species originate from a single gene in the last common ancestor of these species and are therefore likely to share the same function [40,41]. Homologs deriving from gene duplications are called paralogs. Because gene duplication events take place both before and after speciation, paralogy can also exist between genes in different species. Paralogs that originate following a gene duplication after speciation are termed as "in-paralogs" [42].
T. urartu and Ae. tauschii are the progenitors of the A and D genomes of hexaploid wheat (2 N = 6X = 42; AABBDD), respectively. A preliminary analyses have revealed that T. urartu and Ae. tauschii are evolutionarily close. To gain a better understanding of the evolutionary relationship of bZIP genes, InParanoid 7 (Center for Genomics and Bioinformatics, Stockholm, Sweden) was used to identify orthologs and in-paralogs between T. urartu and Ae. tauschii (Fig. 1). Among the bZIP proteins from T. urartu and Ae. tauschii, 65 groups of orthologs that were composed of 68 (69.4 %) T. urartu and 66 (68.8 %) Ae. tauschii proteins were detected using InParanoid 7 (Additional file 2: Figure S1). The largest group of orthologs contained 4 T. urartu proteins and 2 Ae. tauschii proteins and showed a many-to-many type of orthology. In this case, TubZIP18, TubZIP11, TubZIP12 and TubZIP15, as in-paralogs from the same cluster, were co-orthologous to AetbZIP10 and Aetb-ZIP5. According to the clustering algorithm [40], Tub-ZIP11, TubZIP12 and TubZIP15 originated from the gene duplication of TubZIP18, while AetbZIP5 was derived from AetbZIP10. The remaining groups of orthologs between T. urartu and Ae. tauschii represented straightforward one-to-one orthology. Results showed that there were few in-paralogs in whether T. urartu or Ae. tauschii for bZIP gene family. Furthermore, bZIP family members in T. urartu shared high similarity in protein sequence with those in Ae. tauschii, indicating gene functions of bZIP family might remain conserved between T. urartu and Ae. tauschii, which was consistent with the study by Li [43].
Evolutionary relationship of wheat, T. urartu, Ae. tauschii, barley and Brachypodium To further explore the evolutionary relationship of bZIP family genes among these five evolutionarily close species, we identified orthologs and in-paralogs between each pair of these five species ( Fig. 2 and Additional file 1: Table S2) using InParanoid 7, which allowed us to define the evolutionary point of the orthology precisely. To provide a simple explanation, the basic principle was that the number of orthologous groups between closely related species was expected to be greater than that between distantly related species [40]. The number of orthologs in each organism clustered with one other genome was summarized in Table 1. Notably, the table was not a symmetrical table because the gene duplication frequency in organism A generally differed from that in organism B since the speciation of organism A and B. Based on Table 1, wheat had a high average ortholog group size of 1.702 such that every TabZIP gene had an average of 0.702 paralogs. The average ortholog group size in T. urartu, Ae. tauschii, Barley and Brachypodium were found to be significantly smaller than that in wheat, which was consistent with the observation that these four genomes all had relatively fewer in-paralogs compared to the wheat genome [21-23, 32, 44]. Additional file 1: Table S2 also showed the intuitive observation that wheat had many more in-paralogs of bZIP genes than the other genomes.
To reflect the level of orthology of the bZIP family genes between these five species, we constructed a phylogenetic tree termed an orthophylogram [41]. The orthology distance from species A to B (dAB) was calculated using the following formula [41], and the average orthology distance (dAB + dBA)/2 was then used to generate a UPGMA tree [41], as shown in Fig. 3.

dAB ¼ proteins A−proteins A orthologous to B proteins A
This orthophylogram demonstrated quantitatively the level of orthology of the bZIP gene family between different species. It was observed from the orthophylogram that wheat had a closer phylogenetic relationship with barley and Brachypodium than T. uratu and Ae. tauschii.
In plants, the structure and function of most of the bZIP genes most likely remained conserved during angiosperm evolution [2], thus, the results of our study provided deeper insight into the possible evolutionary relationship of five close relatives, namely, wheat, T. urartu, Ae. tauschii, barley and B. distachyon.

Phylogenetic analysis of bZIP genes
The phylogenetic analysis was performed with all identified TabZIP, TubZIP and AetbZIP proteins as well as 76 Arabidopsis and 94 rice bZIP family members (Fig. 4). The TabZIP, TubZIP and AetbZIP proteins were classified primarily into 10 subgroups, which were named A, B, C, D, E, F, G, H, I and S according to a previously described classification method used for Arabidopsis bZIP family [1]. With the newly identified bZIP sequences, we found that subgroups S and I could be further divided into S1 and S2 and into I1 and I2, respectively. U1-U4 represented four previously unnamed clades, including AtbZIP60 and AtbZIP62, which did not fit into any of the 10 known subgroups [1]. Another ungrouped Arabidopsis bZIP member, AtbZIP72, was incorporated into the I subgroup based on our phylogenetic tree. The newly added Arabidopsis bZIP family (See figure on previous page.) Fig. 1 Orthologs between T.urartu and Ae.tauschii and in-paralogs within each organism of them. Each chromosome of T.urartu and Ae.tauschii was figured by actual proportion to others. The chromosomal scales of T.urartu was different from Ae.tauschii. The top and bottom of each stick standing for one chromosome were starting point and terminal point of the chromosome, respectively. The deletion bin maps of T. urartu 1A-7A chromosomes were noted on the left side, and bZIP genes in each deletion bin were averagely distributed. 'uA' and 'uD' mean genome regions not located to specific chromosomes, and genes in these regions were averagely distributed, thus the positions in the figure do not stand for the real locations in the genome members AT1G58110 and AT4G06598 were found to belong to the E subgroup based on the similar distribution of conserved motifs in the protein sequences with other E subgroup members. Similarly, the new member AT2G12940 was found to belong to the I subgroup.
The distribution of bZIP proteins in the subgroups differed among the five species (Additional file 2: Figure S2),  HvbZIP38 HvbZIP56 . "1A"-"7A", "1D"-"7D", "1 Bd"-"5 Bd" and "1H"-"7H" represented chromosomes of T. urartu, Ae. tauschii, B. distachyon and H. vulgare,respectively. The chromosomes for the same species were drawn to scale based on actual proportion, but chromosomes of different species were figured by different scales. The direction for each chromosome from starting point to terminal point was indicated by the red arrow. 'uA', 'uD' and 'uHv' mean genome regions not located to specific chromosomes, and genes in these regions were averagely distributed, thus the positions in the figure do not stand for the real locations in the genome. The light-grey dashed rectangles in each corner just mean linkers between species. Ortholog pairs for the same pairs of species were linked by straight lines of the same color as indicated by the colored block at the bottom. For example, each pair of orthologs between T. urantu and barley were linked with a green straight line suggesting that different plant species underwent different expansions of the bZIP family. The number of bZIP proteins in T. urartu was found to be exactly equal to that in Ae. tauschii for subgroups A, B, C, D, G, H, I, U1 and U3 and differed by one for subgroups F and U2 (Additional file 2: Figure S2). Furthermore, the ratio of total identified bZIP proteins in wheat to that in T. urartu (or Ae. tauschii) was 1.86 (or 1.90). For subgroups D, G, I and U1, the ratio of total identified bZIP proteins in wheat to that in T. urartu (or Ae. tauschii) was approximately 3. Thus, we estimate that approximately 100 bZIP proteins in wheat are not represented in the current draft genome.

Microarray analyses and identification of anther development-related TabZIP genes
To identify TabZIP genes related to anther development, we performed microarray analysis to identify differentially expressed TabZIP genes during the developmental period from the PMC stage to the meiosis stage of wheat anther. The microarray data included five groups of comparisons, namely, between 366A and 411A (CK), 366E and 411E (CK), 366AB and 366EF (CK), 366A and 366E (CK), and 366B and 366 F (CK). "366" and "411" referred to wheat varieties BS366 (thermosensitive genic male sterile line) and Jing411 (restorer line, male-fertile), respectively. The letters "A", "B", "E" and "F" represented the growth conditions of BS366 and Jing411 (A: 10°C with a 12 h light/ 12 h dark photoperiod; B: 10°C with a 14 h light/ 10 h dark photoperiod; E: 20°C with a 12 h light/ 12 h dark photoperiod; F: 20°C with a 14 h light/ 10 h dark photoperiod). "A" and "B" were male-sterile conditions while "E" and "F" were male-fertile conditions for BS366. For example, in comparisons groups 366A-411A (CK), the comparison was established between BS366 and Jing411 grown under the same conditions of "A" and Jing411 was regarded as control.
Altogether, among all TabZIP genes represented on the microarray, 48 TabZIP genes were differentially expressed ( Table 2 and Additional file 1: Table S4). It had been reported that 8 bZIP genes from Arabidopsis, rice and maize had an important role in floral development [16,18,20,[45][46][47][48]. These 8 genes included Atb-ZIP14 (FD) [45], AtbZIP21 (TGA9) [18], AtbZIP27 (FDP) [46], AtbZIP34 [16], AtbZIP46 (PAN) [47], Atb-ZIP65 (TGA10) [18], OsABI5 [20] and maize DLF1 [48]. So 23 candidate bZIP gene were further selected from 48 differentially expressed gene based on their phylogenetic relationship and conserved motifs similarity with the 8 known floral development-related bZIP genes. The phylogenetic relationship and conserved motifs of 23 TabZIP genes and 8 known floral development-related  The number of orthologs in an organism (y-axis) when clustered with another genome (x-axis) was shown on the left of the slash while the number on the right of the slash referred to the number of ortholog groups between two species. Thus, 76 wheat genes had orthologs in T.urartu, which were orthologous to a total of 52 T.urartu genes and wheat and T.urartu had totally 49 ortholog groups. 'bZIP genes' referred to the total number of bZIP genes in one organism. Only one protein for each bZIP gene in five genomes was used in In-paranoid clustering. 'Average size of ortholog groups' referred to the average number of in-paralogs in one ortholog group for an organism. For wheat, average size of ortholog group: bZIP genes were shown in Additional file 2: Figure S3 and S4.

Analysis of gene structures and cis-acting regulatory elements
First, we mapped the exon/intron organization of the 23 TabZIP genes (Fig. 6). We found that 3 genes were intronless and that the numbers of introns for the other TabZIP genes varied from 1 to 10. Based on the method described above, we have identified several pairs of inparalogs within these 23 TabZIP genes (Additional file 1: Table S5). As shown in Fig. 6, the in-paralogs tended to share an identical or similar exon-intron composition. Three in-paralog pairs (namely, TabZIP2 and TabZIP59, TabZIP14 and TabZIP134, as well as TabZIP26 and Tab-ZIP94) had preserved identical gene structure With the newly identified bZIP members, S subgroup were divided into two clades S1 and S2 by the members in C subgroup and the E subgroup seprated protein members in the I subgroup into two clades I1 and I2. U1-U4 were 4 unnamed clades for lack of grouped AtbZIPs organizations, while another in-paralog pair (TabZIP27 and TabZIP148) had a highly constant exon-intron composition except for differing by six nucleotides within the first exon. Gene structure analysis revealed that genes with a close evolutionary relationship tended to share similar gene structures. Divergent exon lengths within the coding sequences of several TabZIP genes could potentially lead to the generation of functionally distinct paralogs.
Expression patterns of TabZIP genes in response to low temperature treatment First, 23 genes were selected from differentially expressed genes based on microarray analyses, indicating these 23 genes were more likely to be involved in anther development. Additionally, the male fertility of BS366 was controlled by temperature. Exposing to 10°C for 5d (starting at the PMC stage and continuing to the meiosis stage) leaded to complete male sterility in BS366 [31]. Thus, we speculate that expression of genes involved in male fertility of BS366 is supposed to be affected by low temperature.
To investigate the expression patterns of TabZIP genes in response to low temperature (10°C) treatment, quantitative real-time PCR (qRT-PCR) was performed using BS366 anthers that were collected at different time points from the PMC to the meiosis stage as samples. Results showed that the expressions of 21 among the 23  TabZIP genes were obviously responsive to low temperature and these genes exhibited complicated expression changes after different durations of treatment (Fig. 7). Clear changes in expression values were observed for these 21 TabZIP genes post-treatment. The transcript levels of six of the 21 TabZIP genes were upregulated at each time point after treatment compared with the untreated control. TabZIP148 exhibited the most remarkable change in expression after treatment because the highest expression value was over 59 times higher than the lowest value. TabZIP153 and TabZIP148 shared a similar changing tendency that transcripts increased gradually and reached a maximum at 12 h, then declined to a relatively lower level at 24 h, followed by an increase throughout the following durations.
Tissue-specific expression patterns of TabZIP genes in the male-sterile line, restorer line and their F1 hybrid Among TY806, BS366 and F1 (TY806 × BS366) grown under local environmental conditions of Fuyang (Anhui Province, China), no differences in the morphology or structure of the stamens and pistils were observed before blooming (Additional file 2: Figure S5). However, after blooming, the anthers of BS366 failed to dehisce as those of TY806 and F1, and a few pollen grains formed and were devoid of starch (Additional file 2: Figure S5). The results of pollen iodine staining demonstrated that the pollen grains of BS366 were small, irregularly shaped, nearly transparent and easily broken compared with those of TY806 and F1. Statistical analysis showed that 95.9 % of the pollen grains in TY806 and 87.0 % in F1 were circular, opaque and dark brown-black after iodine staining, while BS366 exhibited a highly male-sterile type because no Type A pollen grains were present and because 80.6 % of the pollen grains that were present were circular, opaque or partially transparent and light brown-black after iodine staining (Additional file 2: Figure S5).
To explore tissue-specific and variety-specific expression pattern of the 23 TabZIP genes in the wheat cross combination TY806 × BS366, we obtained samples from the roots, stems, leaves and anthers of TY806, BS366 and their hybrid F1 at the meiosis stage to analyze their expression patterns. The results revealed that every TabZIP gene exhibited a distinct tissue-specific expression pattern in TY806, BS366 and their hybrid F1 at the meiosis stage (Fig. 8). Additionally, as in-paralogs, both the TabZIP2/TabZIP59 and TabZIP26/TabZIP94 pairs exhibited similar tissue expression patterns, while TabZIP14 showed a different expression profile compared to its in-paralog Tab-ZIP134. Furthermore, TabZIP27 and TabZIP148 were both expressed at low levels in anther tissue (Fig. 8).
Among the 23 TabZIP genes, 11 were predominantly expressed in the anther, which indicated that they were more likely to have a role in anther development during meiosis.
Furthermore, we selected the anther expression data of these 11 TabZIP genes to analyze the expression differences among the three wheat lines (Fig. 9). Significant differences between the F1 and TY806 or between the F1 and BS366 were assessed using a t-test. These 11 TabZIP genes, excluding TabZIP78, were further classified into three distinct gene expression modes based on their deviation from the mid-parent prediction: over-dominance, 10 TabZIP genes (for each of them, the highest expression level after treatment was more than 10 times over the untreated control); (Right) 11 TabZIP genes (for each of them, the highest expression level after treatment was less than 10 times over the untreated control) under-dominance and low-parent dominance [49]. The high-parent dominant genes and low-parent dominant genes were defined according to the following criterion that the expression level in F1 genotype was significantly different from that in one parent and no significantly different from that in another parent. And the genes were identified as over-dominance or under-dominance when the expression level in F1 genotype was significantly higher or lower than those in both inbred parents [49]. TabZIP134 displayed such extremely clear overdominance that the expression level in the F1 was approximately 6-fold greater than that in BS366, while the corresponding multiple for TabZIP37 was approximately 4-fold. Among the genes with over-dominance mode, Tab-ZIP26, TabZIP37, TabZIP86, TabZIP94 and TabZIP134 showed similar expression trend because they were expressed maximally in the F1 and had the lowest expression level in TY806. TabZIP76 showed an underdominance mode because the slight expression in the F1 was significantly lower than those in TY806 and BS366. TabZIP174 was expressed primarily in TY806 and exhibited slight expression in both BS366 and F1, which belonged to the low-parent dominance mode. In conclusion, the TabZIP genes showed multiple gene expression modes in TY806, BS366 and their F1 hybrid, and these results were consistent with relevant studies in maize and rice that supported the theory that multiple modes of gene action collaboratively contributed to heterosis [49,50].

Discussion
Evolutionary relationship of wheat, T. urartu, Ae. tauschii, barley and Brachypodium The evolutionary relationships between rice, Brachypodium, barley and wheat have been speculated based on the mean synonymous substitution rates (Ks) of orthologous gene pairs [32]. Wheat shared a common ancestor with rice approximately 40-54 million years BC, and the Fig. 9 Differential expression of 11 TabZIP genes in anther of TY806, BS366 and their hybrid F1. These 11 TabZIP genes were predominantly expressed in anther. The expression of each gene in anther of TY806 was regarded as a reference, and other values represented the expression levels of this gene in anther of BS366 and F1 relative to the reference. Single asterisk indicated the expression difference between TY806 and F1 or between BS366 and F1 was significant at p < 0.05 by the t-test and double asterisks meant a highly significant difference between TY806 and F1 or between BS366 and F1. These 11 TabZIP genes were sorted into three distinct gene expression modes based on their deviation from the mid-parent prediction: over-dominance, underdominance and low-parent dominance. TabZIP78 did not belong to any of four kinds of modes (over-dominance, under-dominance, low-parent dominance and high-parent dominance) divergence time of Brachypodium from other gramineous plants was approximately 32-39 million years BC. Subsequently, barley diverged from other gramineous plants approximately 3-4 million years BC [32].
Wheat derived from hybridization between the cultivated tetraploid emmer wheat (T. dicoccoides, AABB) and diploid goat grass (Ae. tauschii, DD) approximately 8,000 years ago [21]. The three diploid progenitor genomes of wheat, AA from T. urartu, BB from a species that might be from the section Sitopsis, and DD from Ae. tauschii, originated from a common Triticeae ancestor approximately 2.5-4.5 million years ago [21].
Notably, T. urartu, Ae. tauschii and barley, all of which were diploid, shared the same chromosome number (seven), suggesting that they had a relatively close evolutionary relationship. Furthermore, synteny block identification among Brachypodium, barley and Ae. Tauschii revealed that they were most likely to share a common ancestral set of five chromosomes [32]. The draft genome sequence of Ae. tauschii revealed 11,289 orthologous gene pairs between barley and Ae. tauschii and 14,675 between Brachypodium and Ae. tauschii [22]. Although the T. urartu genome was more than 18 times larger than that of B. distachyon, the average gene sizes of these two organisms were similar, and the predicted gene number (34,879) in the T. urartu genome was only approximately 1.37-fold that of B. distachyon (25,532) [23]. In conclusion, previous studies have demonstrated that wheat, T. urartu, Ae. tauschii, barley and Brachypodium are evolutionarily close in the Gramineae family. However, the real evolutionary relationship among them still remains mysterious.
bZIP family genes appeared before the divergence between monocots and dicots and the structure and function of most bZIP genes most likely remained conserved during angiosperm evolution [2]. bZIP transcription factors were associated with the evolution of plants [51]. Thus, in this study, we deduced the possible evolutionary relationship among them based on orthology analyses of bZIP transcription factors. In the case of the bZIP family, A sub-genome of wheat was more similar with T.urartu than other three genomes (34.5 % of bZIP family members in A sub-genome of wheat were orthologous to T.urartu). Similarly, among T.urartu, Ae.tauschii, barley and Brachypodium, D sub-genome of wheat was most similar with Ae.tauschii (42.6 % of bZIP family members in D sub-genome of wheat were orthologous to Ae.tauschii). However, B sub-genome of wheat was more similar with barley (49.2 %) and Brachypodium (47.6 %) than T.urartu (30.2 %) and Ae.tauschii  Figure S6). In addition, it was observed from Fig. 3 that wheat had a closer evolutionary relationship with barley and Brachypodium than T. urartu and Ae. tauschii. Approximately 8,000 years ago, wheat arose as a result of hybridization between the cultivated tetraploid emmer wheat (T. dicoccoides, AABB) and diploid goat grass (Ae. tauschii, DD) [21]. It could be deduced that much more than approximately 8,000 years ago, T. dicoccoides (AABB) was generated from the hybridization between T. urartu (donor of A sub-genome) and the donor of B sub-genome. So T. urartu and Ae. tauschii were very old species. From speciation, wheat had undergone at least 8,000 years of evolution course. And after long-term natural selection and artificial domestication and cultivation, A and D subgenomes of modern wheat were no more what they had been. If the self-evolution of T. urartu and Ae .tauschii was taken into consideration, there might exist bigger differences between A sub-genome (or D sub-genome) of modern wheat and modern T. urartu (or Ae. tauschii). However, it was a generally acknowledged fact that barley was one of close relatives of wheat, and Brachypodium had been regarded as a suitable model system for studies on temperate cereals, such as wheat, because of its small and similar genome [52]. So it was understandable that wheat had a closer phylogenetic relationship with barley and Brachypodium than T. uratu and Ae. tauschii that was concluded from orthology analyses. Our results provided deeper insight into evolutionary relationship of wheat and its relatives.

The evolution of bZIP transcription factor family
The phylogenetic analysis of bZIP transcription factors in rice and Arabidopsis indicated that these genes appeared preceding the divergence between monocots and dicots [2]. In plants, the structure and function of most bZIP genes most likely maintained conserved during angiosperm evolution [2]. Indeed, the majority of Osb-ZIP proteins were found to be orthologous to AtbZIP proteins [2]. According to the orthologs identified in wheat and its four relatives (Table 1 and Additional file 1:  Table S2), barley had 72 groups of orthologs with Brachypodium and 64 groups of orthologs with wheat. Notably, 69.4 % of the TubZIP proteins were orthologous to 68.8 % of the AetbZIP proteins ( Fig. 1 and Additional file 2: Figure S1).
During evolution process of bZIP transcription factor family, several new subgroups or clades had come into being in the T. urartu and Ae. tauschii genomes. Based on the phylogenetic tree (Fig. 4), Subgroups A, B, C, D, E, F, G, H, I, and S, as well as U1 and U3, clearly consisted of bZIP proteins that were present in all five species, indicating that the origin of all 10 subgroups and U1 and U3 might have predated the speciation events between eudicot and monocot plants. In contrast, U2 and U4 consisted of only bZIP proteins from T. urartu and Ae. tauschii but not from wheat, Arabidopsis or rice. Previous studies had deduced that wheat and rice share a common ancestor approximately 40 to 54 million years BC [32]. The observation that U2 and U4 bZIP proteins were only present in the T. urartu and Ae. tauschii genomes suggests that these genes likely arose after the evolutionary divergence of rice from other gramineous plants; however, these genes disappeared during the evolution or domestication of wheat, which was consistent with the observation that wheat underwent genetic loss as a consequence of domestication [21].
Additionally, bZIP family went through different degrees of expansion in different gramineous plants. From Table 1, it could be deduced that wheat had many more in-paralogs than its four relatives T. urartu, Ae. tauschii, barley and Brachypodium. On average, each TabZIP gene had 0.702 in-paralogs. Two explanations likely accounted for this high average: (i) wheat was derived from two hybridizations between three diploid progenitors, and as a consequence, several copies of homologous genes had been gathered together in the wheat genome; (ii) a vast number of repetitive genes were present in the wheat genome as a result of gene duplication. By contrast, both T. urartu and Ae. tauschii had few in-paralogs in the case of the bZIP family (Fig. 1, Table 1 and Additional file 1: Table S2), indicating that gene duplication events were likely to occur very rarely within both T. urartu and Ae. tauschii, or that the gene duplication frequency was low during the long evolutionary history of at least 8,000 years because gene duplication can accelerate the generation of in-paralogs within a species [32]. Few gene duplication events might limit the expansion of a gene family such as the bZIP transcription factor family, which might account for the similar sums of bZIP family members between T. urartu and Ae. tauschii.
Additionally, sequence analysis revealed that TabZIP88 and TabZIP175 were highly similar to Arabidopsis TGA9 and TGA10 in the protein sequence, which played a role in anther patterning and dehiscence [18]. And expression analysis showed that TabZIP88 is predominantly expressed in anther (Fig. 8). Therefore, TabZIP88 might share similar function with TGA9 and TGA10. Besides, TabZIP86, which was also predominantly expressed in anther (Fig. 8), had complete sequence and domain identities with OsABI5 of 61.93 % and 93 %, respectively. OsABI5 was required for pollen development and to maintain normal fertility in rice [20].
Expression in anther demonstrated that TabZIP genes exhibited multiple gene expression modes in the TY806 × BS366 heterotic cross (Fig. 9). There were mainly three possibilities accounting for the differential transcription of the same gene in TY806 × BS366: (i) the interplay between promoter INDEL and regulated expression of trans-acting factors (such as transcription factors); (ii) microRNA regulation; (iii) DNA methylation. Previous study have revealed that sequence polymorphism within promoter alleles between inbreds preferentially occurs in those differentially transcribed genes. Promoter INDELs and differential expression of transcription factors may result in all possible modes of gene action [49]. miRNAs generally act as posttranscriptional regulators of gene expression in plants [55][56][57]. DNA methylation significantly contributes to gene expression pattern and the promoter-methylated genes tend to show a greater degree of tissue-specific expression [58,59].
Compared with TY806 and F1, the most special fertility-related characters of BS366 under specific low temperature conditions were that the anthers fail to dehisce normally and that few pollen grains were produced and devoid of starch. As shown in Fig. 9, TabZIP78 was expressed predominantly in anthers of BS366 relative to those of TY806 and their F1 hybrid under low temperature conditions (10°C). Furthermore, TabZIP78 was highly homologous to AtbZIP34, which was involved in the control of several metabolic pathways of developing pollen in Arabidopsis, including pollen wall patterning [16]. Thus, TabZIP78 was probably involved in anther development, however, whether TabZIP78 had a correlation with the character of male sterility in BS366 required further functional verification.

Conclusions
In this study, 187 bZIP family genes were comprehensively identified from the current wheat genome. 98, 96 and 107 members of bZIP family were also identified from the genomes of T.urartu, Ae.tauschii and barley, respectively. The orthologs and in-paralogs between each pair of wheat and its relatives had been also identified. Results showed T.urartu and Ae.tauschii were evolutionarily close and wheat probably had a closer phylogenetic relationship with barley and Brachypodium than T.urartu and Ae.tauschii. Based on microarray data, 48 differentially expressed TabZIP genes were identified to be related to anther development from the comparison between the male-sterile line and the male-fertile line. Genes with close evolutionary relationship tended to share similar gene structures. Expression of 21 among 23 selected TabZIP genes were obviously responsive to low temperature. These 23 TabZIP genes all exhibited distinct tissue-specific expression pattern. Among them, 11 TabZIP genes were predominantly expressed in anther and most of them showed over-dominance expression mode in the cross combination TY806 × BS366. This study provided a comprehensive overview of the bZIP family in wheat and its relatives and laid a foundation for further evolution research and functional characterization of anther development related TabZIP gene.

Identification of bZIP transcription factors in wheat, T. urartu, Ae. tauschii and barley
To identify putative bZIP transcription factor genes in wheat, T. urartu, Ae. tauschii and barley, HMM searches of all annotated proteins in these five genomes were conducted using the HMM profiles as queries with the hmmsearch program, which is included in the HMMER software package [60]. Wheat genomic data were downloaded from the web site Gramene (Release note 41, http://www.gramene.org/) [61]. T. urartu and Ae. tauschii genomic data were downloaded from the web site http://gigadb.org/ [62]. Barley genomic data was downloaded from the web sites ftp://ftpmips.helmholtzmuenchen.de/plants/barley/. The HMM profiles of the bZIP domain (PF00170 and PF07716) were extracted from the Pfam database (http://pfam.sanger.ac.uk/) [63]. And the HMM search was performed using an E-value cutoff < 1.0. The contig sequences were removed, and the remaining protein sequences were verified for the presence of the typical bZIP domain using Pfam with an E-value cutoff < 1.0.

Identification of orthologs and in-paralogs between wheat and other species
Orthologs and in-paralogs between any pair of five species wheat, T. urartu, Ae. tauschii, barley and Brachypodium were identified using Inparanoid 7 software (http://inparanoid.sbc.su.se/cgi-bin/index.cgi) [40]. Two user-defined criteria were applied to each pairwise match: (1) an all-against-all blast score cutoff of 50 bits and (2) an overlap cutoff of 50 %, namely, the matching region must exceed 50 % of the longer sequence in length. Furthermore, a confidence value was defined to show the similarity of each in-paralog to the primary ortholog in the same species on a scale ranging from 0 to 100 % (Additional file 1: Table S2) [40]. A greater confidence value indicated that the in-paralog was more similar to the primary ortholog compared with the other in-paralogs. The orthology figures were generated using a self-writing Perl program employing the SVG module.

Multiple sequence alignment and phylogenetic analysis of bZIP proteins
Multiple sequence alignment of the identified TabZIP, TubZIP, AetbZIP proteins was conducted using the ClustalX program (Version 2.1) with default settings [64]. Meanwhile, 76 Arabidopsis and 94 rice bZIP family members were added to the phylogenetic analyses to determine the classification of TabZIP, TubZIP and Aetb-ZIP proteins because Arabidopsis and rice were model plants for eudicots and monocots, respectively [1,2]. The unrooted tree was generated using ClustalX and displayed by MEGA 5.0 [65].

Microarray analyses
The microarray analyses was primarily based on our previous original Affymetrix microarray data. Samples were collected from wheat cv. Jing411 and BS366 of different developmental stages and under different treatment conditions (10°C with 12 h light/12 h dark, 10°C with 14 h light/10 h dark, 20°C with 12 h light/12 h dark and 20°C with 14 h light/10 h dark). The operating procedures of the whole experimental process were detailedly described in our previous study [31]. For microarray data analysis, the image (cel) files were imported into the online tool Babelomics 4.3 (http://babelomics.bioinfo.cipf.es/) [66]. Subsequently, data normalization was performed using RMA methods, and expression differentiation analysis was conducted using limma methods. Five comparison groups were established: 366A-411A, 366E-411E, 366AB-366EF, 366A-366E and 366B-366 F. "366" represented wheat variety BS366 and "411" represented wheat variety Jing411. The letters "A", "B", "E" and "F" referred the growth conditions of BS366 and Jing411 (A: 10°C with a 12 h light/ 12 h dark photoperiod; B: 10°C with a 14 h light/ 10 h dark photoperiod; E: 20°C with a 12 h light/ 12 h dark photoperiod; F: 20°C with a 14 h light/ 10 h dark photoperiod). And the sampling timepoints were detailed in Additional file 1: Table S8.
Probe sets that were up-or downregulated with p < 0.05 and adjusted p < 0.05 in all comparison groups compared to the corresponding control were considered to be differentially expressed. Nucleic acid and amino acid sequence alignments were performed to establish the corresponding relationship between probe sets and 108,569 wheat genes from the wheat genomic data at the Gramene web site. A probe set was determined to represent one wheat gene using BLASTn and BLASTx tools in accordance with the following criteria: (1) an alignment parameter E-value ≤ 1e-5 between the target sequence of the probe set and CDS sequence of a wheat gene and (2) an alignment parameter Evalue ≤ 1e-3 for the protein sequences between the target of the probe set and this gene. The heatmaps were generated with Cluster 3.0 and displayed by Java Treeview.

Analysis of protein motifs, gene structures and cis-acting regulatory elements
The protein sequences of 23 TabZIP genes were searched for common conserved motifs using the MEME analysis tool (version 4.9.1, http://meme-suite.org/tools/meme) [67], and a limit of 15 motifs was specified with all other parameters set to default. To obtain information regarding the intron and exon structure, the CDS sequences of these 23 bZIP genes were aligned with their corresponding genomic sequences using GSDS (http://gsds.cbi.pku.edu.cn/index.php) [68]. To analyze their promoter regions, the 1.5-kb promoter regions were selected and screened against the Plant-CARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) [69].

Sample collection for morphological observation and pollen iodine staining
To observe the phenotype differences in floral organs among Jing411, TY806, BS366 and the F1 (TY806/ BS366) grown under natural male-sterile conditions of Fuyang, Anhui Province, approximately 60 reproductive organs (each including a pistil and three stamens) were dissected from 60 spikelets, which were collected from 20 primary stems for each wheat line 3 days before blooming. Additionally, approximately 150 anthers were excised from 50 spikelets that were collected from 15 primary stems with small tweezers for each wheat line 1 day after blooming. The representative reproductive organs and anthers were imaged using a Leica MZ16F stereomicroscope (Leica Microsystems, Wetzlar, Germany).
Thirty anthers from each wheat line were used for pollen iodine staining, and the procedures were as follows. (i) A mature anther was positioned on the slide, and a drop of distilled water was placed on the anther. (ii) The anther was squeezed using small tweezers to release the pollen grains. (iii) Then, 1 to 2 drops of I 2 -KI solution [6.7‰ (g/ml) KI and 3.3‰ (g/ml) I 2 ] were added, and a coverslip was mounted on the sample. (iv) The pollen grains were observed after iodine staining, and images were obtained using an Olympus BX41 laboratory microscope (Olympus, Tokoyo, Japan). Based on the morphology of the pollen grains after iodine staining, they were classified into 4 types: (A) circular, opaque and dark brown-black; (B) circular, opaque or partially transparent and light brown-black; (C) circular, transparent and light yellow; or (D) transparent with an irregular shape and light yellow. The number of pollen grains in each classification type was counted, and the percentage of each type for the 4 wheat lines was calculated.

Plant materials and treatments
For the low-temperature treatment experiment, the wheat TGMS line BS366 were grown in the soil in plastic pots buried in the field. The wheat plants underwent vernalization naturally in the field. Lowtemperature treatment was initiated when the flag leaf had half-emerged from the collar of the penultimate leaf (about 1.5 mm in anther length). Plants of uniform size were selected and transferred to the phytotron (Conviron, Winnipeg, Canada) and continue to grow at 10°C with a 12-h light/12-h dark photoperiod for 5d. During this period, anthers were harvested at 0 h, 5 h, 12 h, 24 h, 3d and 5d post-treatment with home-made dissecting needles under the dissecting microscope (Olympus SZX12, Tokoyo, Japan) and used for the expression analysis. To maximize the sampling consistency, about 30 anthers were separated from 10 spikelets, which were collected from 3 primary stems at each sampling timepoint.
For tissue-specific expression analyses of 23 TabZIP genes in the wheat cross combination TY806 × BS366, plant materials were collected from the wheat varieties TY806, BS366 and their F1 hybrid grown in the field in Fuyang, Anhui Province, where BS366 was a male-sterile variety, and TY806 and F1 were male-fertile varieties. The roots, stems, leaves and anthers of the three varieties were harvested 5d after the flag leaf had half-emerged from the collar of the penultimate leaf and three biological replicates were harvested for each sample.
After sampling, all plant materials used for the above studies were frozen immediately in liquid nitrogen and stored at −80°C.

RNA extraction and qRT-PCR
For the expression analysis of the bZIP genes, total RNA was extracted from plant tissues using TRIzol reagent (Ambion, USA) according to the manufacturer's instructions. First-strand cDNA synthesis was performed using a TaKaRa PrimeScript™ RT Reagent Kit with gDNA Eraser (TaKaRa, Dalian, China). qRT-PCR analysis was conducted using an Eco Real-Time PCR system (Illumina, San Diego, CA, USA) with TaKaRa SYBR® Premix Ex Taq™ (Tli RNase H Plus) (TaKaRa, Dalian, China). Wheat Actin (Gene ID: 542814) served as an internal control for the expression studies.
The primers for all genes were designed using the Primer Premier 5.0 program [70]. Each reaction was performed in triplicate in a reaction volume of 10 μl. The qRT-PCR parameters were as follows: 95°C for 30 s, 40 cycles of 95°C for 5 s and 55°C for 30 s. For the melting curve analysis, a program of 95°C for 15 s followed by a constant increase from 55°C to 95°C was included after the PCR cycles. The expression analysis of the bZIP genes and wheat actin gene was performed using the same PCR procedure as detailed above or with a slightly adjusted annealing temperature. The relative gene expression levels were analyzed according to the 2 −ΔΔC T method [71].

Availability of supporting data
All the supporting data is included within the article and its additional files.

Additional files
Additional file 1: Table S1. Nomenclature and chromosomal location of bZIP family genes in wheat, T.urartu, Ae.tauschii, barley and Brachypodium. Table S2. Orthologs and in-paralogs in bZIP family between each pair of wheat and its relatives. Table S3. Numbers of differentially expressed bZIP genes in five comparison groups. Table S4. Lists of differentially expressed bZIP genes in five comparison groups. Table S5. In-paralog pairs within 23 TabZIP proteins we studied. Table S6. Distribution of cis-acting elements in promoter regions of TabZIP genes. Table S7. Functions of cis-acting elements in promoter regions of TabZIP genes. Table  S8. Wheat lines, growth conditions, and sampling timepoints for microarray data analysis. The sample features were represented by the wheat line name and a letter (A, B, E or F) followed by a number referring to the anther length when samples were taken. The letters 'A', 'B', 'E' and 'F' indicated the growth conditions of wheat lines. 'A': 10°C with a 12 h-light/12 h-dark photoperiod; 'B': 10°C with a 14 h-light/10 h-dark photoperiod; 'E': 20°C with a 12 h-light/12 h-dark photoperiod; 'F': 20°C with a 14 h-light/10 h-dark photoperiod. (ZIP 146 kb) Additional file 2: Figure S1. Orthologous bZIP genes between T.urartu and Ae.tauschii. Figure S2. Distribution of bZIP family genes in subgroups and clades. bZIP family genes of wheat, T.urantu, Ae.tauschii, rice and Arabidopsis were divided into 10 subgroups and 4 clades. The number in the bracket indicated the total number of bZIP family genes in the corresponding genome. Figure S3. Phylogenetic relationship of 23 TabZIP proteins and 8 known floral development related bZIP proteins. Figure S4. Distribution of conserved motifs in 23 TabZIP and 8 known floral development related bZIP proteins. Figure S5. Morphological characteristics of floral organs and pollen iodine staining. The wheat materials were taken from Fuyang (Anhui Province), where Jing411, TY806 and F1 were male-fertile while BS366 was male-sterile. The first row showed the stamen and pistils before blooming; The second row showed anthers after blooming; The third row showed iodine staining results of pollen grains; The fourth row showed the statistical data of pollen iodine staining. The total pollen count used for data statistics was about 250 in four wheat lines. A, B, C and D meant 4 different kinds of shape-color types for pollen grains after iodine staining: (A) circular, opaque and dark brown-black; (B) circular, opaque or partially transparent and light brown-black; (C) circular, transparent and light yellow; (D) transparent in an irregular shape and light yellow. Figure S6. Number of bZIP family members in each genome (or sub-genome) and number of ortholog pairs between each pair of genomes and sub-genomes. Each box presented a genome or sub-genome, and the numbers in the boxes referred to the sums of bZIP family members in respective genomes. The number besides the straight line connecting each pair of genomes and sub-genomes indicated the sum of ortholog pairs between them. (ZIP 2480 kb)