BACKGROUND
a defining feature of transcription factors is that they contain one or more sequence-specific dna-binding domains that bind to the promoter and/or enhancer regions of target genes to regulate gene expression  <cit> .
the basic leucine  zipper  transcription factor family, one of the most diverse transcription factors families, is characterized by a highly conserved bzip domain which is 60– <dig> amino acids in length and composed of two parts: a basic region and a leu zipper  <cit> .
the basic and leu zipper regions are structurally and functionally distinct.
the basic region comprises approximately  <dig> amino acid residues with the invariant motif n-x7-r/k-x <dig> and is responsible for nuclear localization and dna binding, whereas the leu zipper is composed of heptad repeats of leu or other bulky hydrophobic amino acids positioned exactly nine amino acids towards the c-terminus and mediates homo- and/or heterodimerization of bzip proteins  <cit> .
bzip transcription factor encoding genes have been identified extensively in plants including arabidopsis  <cit> , rice  <cit> , sorghum  <cit> , maize  <cit> , grapevine  <cit> , cucumber  <cit> , castor bean  <cit>  and barley  <cit>  with the availability of their whole genome sequences.
like other transcription factors, members of the bzip transcription factor family are expressed constitutively or in an organ-specific  <cit> , stimulus-responsive  <cit> , development-dependent  <cit>  or cell cycle-specific  <cit>  manner in plants.
it has been reported that bzip transcription factors are involved in various biological processes including organ and tissue differentiation, embryogenesis, seed maturation and storage protein gene regulation, floral transition and initiation, and vascular development .
moreover, bzip transcription factors are also regarded as important regulators in signaling and responses to abiotic/biotic stimuli, including abscisic acid  signaling, hypoxia, drought, high salinity, cold stress, hormone and sugar signaling, light responses, osmotic stresses and pathogen defense .
in soybean, three bzip genes were found to function as negative regulators of aba signaling and confer salt and freezing tolerance in transgenic arabidopsis  <cit> .
in recent years, legume genome sequencing projects have been initiated and completed in lotus japonicus   <cit> , soybean   <cit> , medicago truncatula   <cit> , pigeonpea   <cit> , chickpea   <cit>  and common bean   <cit> .
these six legumes belong to the large papilionoideae subfamily and fall within two sub-clades of the papilionoid legumes: the phaseoloids  and hologalegina .
the phaseoloids are mostly tropical and include the crops soybean , pigeonpea  and common bean , while the hologalegina are mostly temperate and include medicago , chickpea  and lotus .
among the widespread genome duplications throughout the history of flowering plants  <cit> , two recent whole genome duplication  events have affected the evolution of legume genomes.
the older polyploidy event, shared by all legumes, is estimated to have occurred 56– <dig> million years ago   <cit> , while the more recent genome duplication event occurred up to  <dig> mya only in the lineage leading to glycine  <cit> .
genome duplication and subsequent fractionation have played key roles in shaping present-day legume genomes and also the sizes of gene families  <cit> .
with the availability of these legume genome sequences, the members of the bzip transcription factor family were systematically investigated and analyzed in this study.
we identified all legume bzip genes and analyzed their bzip domain sequences, gene structure and additional meme motifs, which was in agreement with and supported the phylogenetic classification.
then, we predicted the dna-binding-site specificity and dimerization properties of the legume bzip proteins.
we also investigated the impact of the two legume-lineage wgds and tandem duplication on the expansion of the legume bzip gene family.
by analyzing their expression profiles, legume bzip genes constitutively or specifically expressed in different tissues and seed developmental stages were identified, as well as candidate legume bzips responsive to drought and salt stresses.
methods
identification of bzip transcription factors in six legume genomes
all genomic sequences and annotated proteins of the six legumes were downloaded from ftp://ftp.jgi-psf.org/pub/compgen/phytozome/v <dig> /gmax/ , http://jcvi.org/medicago/display.php?pagename=general&section=download , http://genome.jgi.doe.gov/pages/dynamicorganismdownload.jsf?
organism = phytozomev <dig> , http://cicar.comparative-legumes.org/ , http://www.icrisat.org/gt-bt/iipg/genome_manuscript.html  and ftp://ftp.kazusa.or.jp/pub/lotus/lotus_r <dig> / .
to identify all the possible bzip proteins in the six legume genomes, both local blast and hidden markov model  searches were performed.
for blastp, the known bzip proteins from arabidposis  <cit> , rice  <cit>  and maize  <cit>  were used as queries and the e-value was set to 1e- <dig>  for the hmm search, the profile of bzip domain was used and the e-value threshold was set at  <dig>  the sequences were further analyzed to confirm the presence and integrity of the bzip domain through the expasy proteomics server   <cit>  and interpro   <cit> .
all bzip domain sequences were aligned using mafft  <dig>  <cit>  to manually check and remove sequences with incomplete domains.
the nomenclature was based on the exact positions of the bzip genes on the chromosomes/scaffolds from top to bottom.
distinct transcripts encoded by the same gene locus shared the same gene number with an additional decimal part, such as point  <dig> or  <dig> .
sequence alignment and phylogenetic analysis
the bzip amino acid sequences from a. thaliana and the six legume genomes were aligned using clustalx  <dig>   <cit>  with gap opening and gap extension penalties of  <dig> and  <dig> , respectively.
the phylogenetic tree was reconstructed by the maximum likelihood  method using the phyml  <dig>  software  <cit> .
jtt + g was selected as the best model for constructing the phylogenetic tree by the akaike information criterion implemented in prottest  <dig>   <cit> .
bootstrap values from  <dig> replicates were indicated at each node.
mega <dig>  <cit>  was used to show the tree.
structure of bzip genes
the positional information of both the gene sequence and the corresponding coding sequence were loaded into the gene structure display server v <dig>    <cit>  to obtain information on the intron/exon structure.
the coordinates of the bzip domain in each protein were recalculated into the coordinates in gene sequence and featured in gene structure.
we used genewise  <cit>  to analyze the intron distribution pattern and intron splicing phase within the basic and hinge regions of the bzip domains in the six legumes.
detection of additional conserved motifs
to identify additional conserved motifs outside the bzip domain of legume bzip transcription factors, we used the multiple em  for the motif elicitation tool   <cit> .
the limits for maximum width, minimum width and maximum number of motifs were specified as  <dig>   <dig> and  <dig>  respectively.
fifty motifs were finally confirmed because of their low e-values .
the motifs were numbered according to the order displayed in meme and were considered as group-specific signatures for their presence of high frequency in the given groups.
detection of duplicated genes and estimation of synonymous  and nonsynonymous  substitutions per site and their ratio
the duplicated gene pairs derived from segmental duplication were identified in the legume genomes based on the method from the plant genome duplication database  <cit> .
an all-against-all blastp comparison  provided the gene pairs for syntenic clustering determined by mcscan  .
tandem duplication arrays were identified using blastp with a threshold of e < 10− <dig>  and one unrelated gene among cluster members was tolerated, as described in a. thaliana  <cit> .
pairs from segmental or tandem duplications were used to estimate ka, ks and their ratio.
amino acid sequences from segmentally or tandemly duplicated pairs were first aligned and then guided and transferred into a cdna sequences alignment using in-house perl scripts.
then, the software kaks_calculator was used to compute ka and ks values for each pair following the yn model  <cit> .
expression analysis of legume bzip genes
for different tissues/organ and seed developmental stages, the normalized counts for bzip genes from rna-seq were obtained from the soyseq   <cit>  and pvgea  databases  <cit> .
microarray expression values for medicago and lotus were downloaded from http://mtgea.noble.org/v3/  <cit>  and http://ljgea.noble.org/ v2/  <cit> .
to identify candidate bzip genes responsive to drought and/or salt stresses, microarray gene expression data in medicago were downloaded from zhang et al.
<cit>  and li et al.
<cit> , respectively.
the gene expression changes in drought-stressed roots and shoots  were obtained by comparing with levels in the watered control .
similarly, the fold changes in gene expression were calculated when comparing the salt-stressed roots  and control .
the corresponding relationships between microarray probes and legume bzip genes were built using blast .
the expression values or normalized counts were log10-transformed and the gplots package was used to make a heatmap in r.
RESULTS
identification and nomenclature of the legume bzip transcription factor family
through sequence similarity and domain searches,  <dig> gm,  <dig> mt,  <dig> pv,  <dig> ca,  <dig> cc and  <dig> lj bzip genes, encoding  <dig>   <dig>   <dig>   <dig>   <dig> and  <dig> distinct proteins, respectively, were identified in the six legume genomes.
based on their exact positions on chromosomes/scaffolds , we gave a unique name to each bzip protein.
the related information on bzip transcription factors are listed in the additional file  <dig>
domain analysis showed that all  <dig> of the bzip transcription factors except eight had a typical bzip domain with an invariant n-x7-r/k motif in the basic region and a heptad repeat of leu or other bulky hydrophobic amino acids positioned exactly nine amino acids upstream of r/k toward the c-terminus .
of the remaining eight, in glyma12g <dig>   and phvul.011g <dig>  , the conserved asn  in the basic region was replaced by lys .
in glyma03g <dig>  , glyma19g <dig>  , medtr7g <dig>   and ca_ <dig> , the conserved arg/lys  in the basic region was substituted by an ile , whereas in glyma11g <dig>  , the conserved arg/lys  in the basic region was substituted by a trp .
in c. cajan <dig> , the heptad repeat of leu was positioned at  <dig> amino acids toward the c terminus instead of the usual nine.
all of these unusual changes in the bzip domain have been observed in rice bzip sequences previously  <cit> .
phylogenetic analysis and classification of legume bzip genes
to investigate the phylogenetic relationships of the bzip transcription factors,  <dig> protein sequences from the six legumes and  <dig> protein sequences from arabidopsis   <cit>  were analyzed .
in accordance with the bzip classification in arabidopsis  <cit> , the phylogenetic tree was subdivided into  <dig> clades with well-supported bootstrap values.
all groups contained legume bzip proteins clustered together with atbzip proteins in the same clade, except group f, which included members from only three genomes: g. max, p. vulgaris and c. arietinum.
in addition,  <dig> legume bzip proteins and  <dig> atbzip genes  formed two small unique and several other individual clades , which were classified into the unclassified group  based on their possible independent evolutionary trajectories from other clades.
group classification was supported by the group-specific sequence characteristics identified in the following analyses of gene structures, intron phases, additional conserved motifs outside the bzip domain and dna-binding site specificity in each group.
it is evident that the group-specific sequence characteristics of the bzip members formed before the divergence of arabidopsis and legumes since conserved sequence characteristics were present in the same group containing both arabidopsis and legumes bzip genes.
nevertheless, it seems that intra-species duplication and parallel evolution of the bzip family in each legume has occurred afterward and contributed to the member variation in each group.fig.
<dig> phylogenetic tree of legume and arabidopsis bzip genes.
bzip protein sequences were aligned by clustal x and the phylogenetic tree was constructed using phyml by the maximum likelihood method.
bootstrap values are based on  <dig> replicates.
genes on branch ends from different species are denoted by different colored circles.
the legume bzip proteins were grouped into  <dig> distinct clades , which are indicated by colored branches
gene structure of legume bzip genes
the intron-exon organization can reflect the evolutionary trajectory of gene families  <cit> .
we examined the gene structures of all  <dig> legume bzips and found that the structural patterns were similar among members within each group but distinct between different groups .
the number of introns in each group was uneven but relatively concentrated .
we detected  <dig> gmbzip genes ,  <dig> mtbzip genes ,  <dig> pvbzip genes ,  <dig> cabzip genes ,  <dig> ccbzip genes  and  <dig> ljbzip genes  with no introns.
most of these intronless genes were clustered into groups s and f .
among the intron-containing bzip genes, the number of introns within the open reading frame  varied from  <dig> to  <dig> in gmbzip,  <dig> in mtbzip,  <dig> in pvbzip,  <dig> in cabzip,  <dig> in ccbzip, and  <dig> in ljbzip, close to the highest number of introns within the orf reported in arabidopsis   <cit> , rice   <cit> , sorghum   <cit> , maize   <cit> , castor bean   <cit> , barely   <cit>  and cucumber   <cit> .
the bzip genes with the most introns were commonly found in groups d and g .
the intron positions within the orf were diverse.
the phases of the splicing sites within the orfs also differed, but the positions and phases of introns in the basic and hinge regions of the bzip domain were highly conserved.
the legume bzip genes showed four intron patterns  based on the intron positions, presence or number and splicing phases within the basic and hinge regions .
pattern a, having one intron in phase  <dig>  within the hinge region at the − <dig> position, was seen in all members of groups a and g in the legumes.
pattern b, having two introns each in phase 0: one in the basic region at the − <dig> position and the other in the hinge region at the − <dig> position, was seen in all members of group d in the legumes.
pattern c, having a single intron in phase  <dig>  at the − <dig> position in the basic region, was seen in all members of groups c, e and h in the legumes.
pattern c also was seen in all members of group i in g. max, m. truncatula, c. arietinum and c. cajan, and most members  of group i in l. japonicus.
pattern d was no introns in the basic and hinge regions and was seen in groups b, s and f, which contained most of the intronless genes.
among  <dig> genes showing pattern d,  <dig> were intronless, while the remaining  <dig> had introns outside the basic and hinge regions.
in summary, the overall pattern of intron position acts as an index for the group classification and phylogenetic relationships in the legume bzip gene family.
the splicing phase has remained well conserved during the course of evolution in legume bzip genes.fig.
<dig> intron patterns within the basic and hinge regions of the legume bzip domain.
the primary structure of the bzip domain is shown at the top of the picture.
p <dig> indicates the intron splicing site is between codons, and p <dig> means the intron splicing site is located between the second and third nucleotides of one codon.
based on the intron incidence and positions, as well as the splicing phase, the legume bzip genes were divided into four patterns .
details of the intron positions within the bzip domain in the legume bzip proteins are shown in additional file 2
identification of additional structural features in the legume bzip genes
all legume bzip protein sequences were loaded into the meme analysis tool and a total of  <dig> additional conserved motifs outside the bzip domain were identified.
the multi-level consensus sequences and the amino acid lengths of these conserved motifs are given in additional file  <dig>  the legume bzip proteins within the same group had similar motif compositions, suggesting conserved evolution and supporting the group classification .
additionally, some motifs were shared by different groups, such as motifs  <dig> and  <dig> in two groups, motif  <dig> in three groups, motifs  <dig> and  <dig> in four groups, and motifs  <dig>  and  <dig> in five groups.
nevertheless, most of the conserved motifs appeared specific to each group  and therefore the group-specific motifs could help determine the specific functions of members in each group.
notably, we did not detect any additional conserved motifs outside the bzip domain in groups b or f.fig.
<dig> distribution of additional conserved motifs identified by meme.
motif compositions based on the position of the bzip domain and additional conserved motifs outside the bzip domain in representative legume bzip proteins for each group are shown.
the bzip domains are shown in pink and different motifs are highlighted in different colored boxes with numbers  <dig> to  <dig>  the motifs shared by different groups are marked with black stars.
details for the predicted conserved motifs are given in additional file 6
a few of these motifs in legume bzips have been studied for their possible biological functions.
a part of motifs  <dig>   <dig> and  <dig> represent potential casein kinase ii  phosphorylation sites , indicated by the motif patterns e, tlge and tvde.
motif  <dig> also contained a phosphorylation site for the ca2+-dependent protein kinase , presented as rqs.
motif  <dig> in group h was group-specific, and contained potential ck ii phosphorylation sites , indicated by cye.
such motifs were identified in group a exclusively.
in addition, all members in group d shared motif  <dig>  which was a dog <dig> domain, and motif  <dig>  the function of which was unclear .
motifs  <dig>   <dig> and  <dig> were observed in group g exclusively and were characterized by a part of the proline-rich domain, which has been shown to have transcriptional activation potential.
interestingly, there were some common motifs among the six legumes, maize and rice  <cit> .
for example, motif  <dig> in group a was the same as motif  <dig> in group a in maize.
motifs  <dig>   <dig> and  <dig> in group d are in common with motifs  <dig>   <dig> and  <dig> in group d of maize and motifs  <dig>   <dig> and  <dig> of rice, respectively.
prediction of dna-binding-site specificity in legume bzip proteins
experiments of mutant proteins demonstrated that the bzip tfs binding specificity is independently determined by the core basic region and the hinge region, and the two regions have an additive effect on dna-binding specificity  <cit> .
to predict the dna-binding-site specificity of the legume bzip proteins, the amino acid sequences of the basic and hinge regions of  <dig> legume bzip proteins were aligned, revealing some highly conserved amino acid residues within each group .
we can predict the dna-binding specificity in a group manner, as described in additional file  <dig>  furthermore, the amino acids were numbered as previously reported  <cit> , and the first leu in the leu heptad repeats was numbered + <dig>  for the two invariant sites asparagine  and arginine , numbered − <dig> and − <dig>  respectively, new dna-binding specificities will occur if other amino acids functionally replace these two invariant sites   <cit> .
nevertheless, such replacements were infrequent and occurred only in groups g and u in the legume bzip proteins.
at position − <dig>  the conserved asparagine  was replaced with lysine  in gmbzip <dig> and pvbzip <dig> in group g .
the same replacement  was also observed in four bzips in maize  <cit> , two bzips in barley  <cit> , and two bzips in castor bean  <cit> .
at position − <dig>  nine members in group u had a hydrophobic isoleucine  residue instead of arginine  or lysine .
an identical substitution pattern  was observed in two bzips in maize  <cit> , two in barley  <cit> , one bzip in grapevine  <cit> , and one bzip in castor bean  <cit> .
it was demonstrated that an arginine to isoleucine mutation in the basic domain of the yeast bzip factor gcn <dig> completely inhibited its affinity for the ap <dig> site  <cit> .
additionally, it was reported that oszip-2a belonging to group u in rice, because of this replacement, does not recognize g-boxes  <cit> .
this evidence suggests that unusual substitutions in the dna-binding domain affect the dna-binding specificity.
these predictions were made to facilitate further studies on the dna-binding patterns of the legume bzip transcription factors.
prediction of dimerization properties in legume bzip proteins
studies have demonstrated that the leu zipper region of the bzip domain, arranged in the form of heptad repeats, mediates homo- and/or heterodimerization between the parallel coiled-coil structures .
within each heptad, the amino acid positions are recognized in order as g, a, b, c, d, e, and f  <cit>  .
leu zipper oligomerization, dimerization stability and specificity are determined mainly by the four amino acids present at the a, d, e and g positions because of their special positions near the leu zipper interface.
the a and d residues are typically hydrophobic on the surface of the helix and create a hydrophobic core that promotes the interaction between two monomers and is essential for dimer stability  <cit> .
the a position contains asparagine , which can form a polar pocket in the hydrophobic interface that limits oligomerization in interhelical interactions  <cit>  and produce more stable n–n interactions at a↔aʹ  than other amino acids  <cit> .
the conserved leu at the d position, one of the most stabilizing aliphatic amino acids  <cit> , is important to maintain dimer stability.
however, the e and g positions that flank the dimerization interface frequently contain charged amino acids including the acidic amino acids glutamic acid  and aspartic acid , and the basic amino acids arginine  and lysine , which are thought to form salt bridges between helices in electrostatic interactions  <cit> .
to predict the dimerization properties of the legume bzip transcription factors, a detailed analysis was carried out to characterize the amino acids present at the a, d, e and g positions .
figure 4a shows the composition of different kinds of amino acids found in the a, d, e and g positions in the six legume bzips, respectively.
at the a position, about 25 % of the residues were asparagine , suggesting that there will be a greater number of homodimerizing leu zippers through stable n–n interhelical interaction at the a↔aʹ position among legume bzip proteins.
the frequency of asparagine  residues in the a position was highest in the second heptad followed by the fifth heptad, accounting for  <dig>  % and  <dig>  %, respectively , similar to earlier observations for atbzip proteins  <cit> .
at the d position, the frequency of leu, responsible for dimer stability, was about 68 % in the legume bzips , which is a little less than in osbzips   <cit>  and zmbzips   <cit> , but significantly greater than in atbzips   <cit> .
at the e and g positions , the frequencies of charged amino acids  were  <dig> and 59 %, respectively.fig.
<dig> prediction of dimerization properties of the legume bzip proteins.
a pie charts presenting the frequency of amino acids in all the a, d, e and g positions of the leu zipper in legume bzip domains.
b histogram of the frequency of asn  in the a positions of the leu zippers among legume bzip proteins.
c histogram of the frequency of attractive or repulsive g↔eʹ pairs per heptad for legume bzip proteins
in electrostatic interactions, the amino acids at position g and oppositely charged amino acids at position eʹ  can form interhelical salt bridges that determine the dimerization specificity and stability  <cit> .
the dimerization specificity of human, drosophila, arabidopsis, rice, maize, barley and castor bean bzip tfs has been predicted on the basis of attractive or repulsive interhelical g↔eʹ electrostatic interactions  <cit> .
to analyze the contribution of charged residues to the dimerization properties of the legume bzip proteins, we calculated the frequency of attractive and repulsive g↔eʹ pairs in each heptad of the bzip leu zippers; the corresponding histograms are shown in fig.
4c.
if both the g and corresponding e position amino acids are charged, they are referred to as complete g↔eʹ pairs.
this analysis was carried out on the basis of the frequency of attractive and repulsive g↔eʹ pairs, which were classified into four groups, attractive basic-acidic pairs , attractive acidic-basic pairs , repulsive basic pairs  and repulsive acidic pairs  in the heptads.
attractive g↔eʹ pairs were predominant in the second, fifth and sixth heptads, thereby suggesting the chances of heterodimerization .
in contrast, few complete g↔eʹ pairs were observed in the eighth heptads except two attractive acidic-basic  pairs from gmbzip <dig> and pvbzip <dig> .
moreover, only attractive acidic-basic  pairs were present in ninth heptads , which was similar to observations in hvbzips  <cit> , osbzips  <cit>  and zmbzips  <cit> .
we divided the  <dig> legume bzip proteins into  <dig> subfamilies  on the basis of the defining properties of dimerization specificity  <cit> .
these subfamilies were divided into three general groups:  those that strongly favor homodimerization within eight subfamilies ,  those with both homo- and heterodimerization properties  and  those that strongly favor heterodimerization within five subfamilies .
the results indicated the complexity and diversity of dimerization patterns in legume bzip proteins, with the potential to homodimerize with themselves or with members in the same subfamily as well as heterodimerize with other subfamily members, which has been shown in maize  <cit> , rice  <cit> , and arabidopsis  <cit> .
based on the criteria used to define the boundaries and natural c-terminus, we observed that the length of the leu zipper in the bzip transcription factor family was variable, ranging from two to nine heptads.
among the bzip proteins,  <dig>  % had only two short zippers and belonged to bz <dig> and group g, about 28 % had only three short zippers , and about 9 % had no α-helix breakers for  <dig> or more heptads, were mainly distributed in bz37–bz <dig> and included members of group i.
the impact of whole genome duplication and tandem duplication on the expansion of the legume bzip gene family
to explore the expansion mechanism in detail, we analyzed the contribution of wgd –derived segmental duplication and tandem duplication to the expansion of the legume bzip gene family.
first, a phylogenetic tree was constructed based on the concentrated orthologous sequences of rbc and matk genes from the six legume genomes, which supported the classification of two sub-clades of papilionoideae legumes: phaseoloids  and hologalegina  .
we carried out a genome-wide identification of collinear duplicated blocks derived from segmental duplication in each species, and then examined the pairwise synonymous distances  of paralogs within duplicated collinear blocks.
by plotting their distribution, two distinct peaks were found : one was specific to gm-gm paralogues, while the other was observed in all six legumes.fig.
<dig> whole genome duplication -derived legume bzips genes.
a the phylogenetic relationships of the six legumes based on the concentrated orthologous sequences of rbc and matk genes.
two wgd events, the glycine-lineage-specific and early-legume wgd events, are indicated on the corresponding clades.
b the ks distribution of paralogs derived from wgd-derived duplicated genomic blocks in each legume.
two obvious peaks correspond to the old and recent legume-lineage wgd events.
c, d the black lines in the ideogram show the chromosomal positions of all identified bzip genes, and the duplicated bzip pairs are linked by lines in soybean  and common bean .
e patterns of gene retention or loss indicated by two different informative tree topologies using common bean and medicago as outgroups, respectively.
the numbers for each pattern are shown in parentheses.
f the ks  and ka  distribution for each duplicated legume bzip gene pair.
the red dashed line means the ka/ks ratio is equal to  <dig>  particularly, eight gm duplicated pairs  and one mt pair , whose ka/ks values were greater than  <dig> , are indicated in the inner figure .
the red dashed line in inner figure means the ka/ks ratio is equal to  <dig> .
the numbers  correspond to duplicated pairs: 1: gmbzip23-gmbzip131; 2: gmbzip49-gmbzip109; 3: gmbzip55-gmbzip107; 4: gmbzip8-gmbzip113; 5: gmbzip4-gmbzip69; 6: gmbzip26-gmbzip42; 7: gmbzip87-gmbzip101; 8: gmbzip65-gmbzip113; 9: mtbzip2-mtbzip26
second, the chromosomal distribution of the legume bzip genes was plotted and the bzip gene pairs on duplicated chromosomal collinear segments were connected by lines .
among these genes, some were segmentally duplicated once and some were duplicated twice or thrice.
the duplication occurred within a chromosome or between chromosomes.
we detected  <dig> gm,  <dig> mt,  <dig> pv,  <dig> ca,  <dig> cc and  <dig> lj bzip genes involved in segmental duplication, accounting for around  <dig>  % ,  <dig>  % ,  <dig>  % ,  <dig>  % ,  <dig>  % , and  <dig>  %  of the bzip genes in each species.
the higher ratio in soybean reflects the preferential gene retention after multiple rounds of wgd, while the different ratios among the other five legumes may be mainly attributable to the genome assembly quality  or species-specific evolution in each genome.
moreover, we roughly identified different wgd event origins for duplicated bzip gene pairs according to their pairwise synonymous distances using criteria from soybean  <cit> : ks values of  <dig> – <dig>  correspond to the 13-mya glycine-lineage-specific wgd, and ks values of  <dig> – <dig>  in soybean and 0– <dig>  in the other five legumes correspond to the 59-mya early-legume wgd, while larger ks values correspond to more ancient wgd events like the ‘gamma’ event  <cit> .
by ordering the ks values,  <dig> gm bzip gene pairs were associated with the  <dig> mya glycine-lineage-specific wgd and  <dig> gm,  <dig> mt,  <dig> pv,  <dig> ca,  <dig> cc, and  <dig> lj bzip gene pairs were associated with the  <dig> mya early-legume wgd; the others were associated with more ancient wgds .
because only soybean has undergone the glycine-lineage-specific wgd, the other five legume genomes could be considered putative ancestors for investigating gene retention and loss after the recent genome duplication in soybean.
mt and pv were chosen, because of their good genome assembly and annotation, to construct phylogenetic trees with gm genes.
we identified informative tree models of gene retention:  and , and gene loss  and , showing that 83 %  and 87 %  ancestral loci were retained as paralogs in soybean after the recent wgd .
compared with the average  <dig>  % of genes retained as paralogs at the whole genome level  <cit> , genes of the bzip family were highly retained.
to better understand the evolutionary constraints acting on the legume bzip genes, we calculated the ka/ks ratios for all legume duplicated bzip gene pairs.
the resulting pairwise comparison data showed that the ka/ks values of only eight gm duplicated pairs and one mt pair were larger than  <dig>   while all of the remaining ka/ks ratios were less than  <dig>  , suggesting that the bzip family has mainly undergone strong purifying selection and the legume bzip genes are slowly evolving at the protein level.
we further compared the strength of selection on bzip paralogs between the recent and old whole-genome duplications in soybean.
the average ka/ks ratio for the recently duplicated bzip genes  was higher than that of the early duplicated bzip genes , and there was a significant difference between these ratios .
this indicated that the younger bzip proteins may be under stronger evolutionary constraints than older proteins, but supported the notion that the legume bzip gene family is essential for the regulation of cellular processes.
third, we determined members with tandem duplication in each legume .
we detected only  <dig> gm,  <dig> mt,  <dig> pv,  <dig> ca,  <dig> cc and  <dig> lj tandem gene pairs, indicating the limited contribution of tandem duplication to the expansion of the gene family.
therefore, it seems that segmental duplication rather than tandem duplication is the major mechanism driving the expansion of this gene family.
expression analysis of bzip transcription factors
expression data from different tissues  and seed developmental stages in soybean  <cit> , medicago  <cit> , common bean  <cit>  and lotus  <cit>  were downloaded, and hierarchical clustering was performed to visualize a global transcription profile of the legume bzip genes.
as shown in fig.
<dig> , the heatmaps were always divided into three recognized clusters, which was similar to the results in rice  <cit> , maize  <cit>  and other plants  <cit> .
the different clusters corresponded to overall differences in expression patterns including expression values and specificity across tissues.
the latter was indicated by cv values , which were calculated for each gene, and helped to recognize genes expressed in specific tissues or stages .
cluster i included genes with relatively high expression levels and the least expression variability , indicating an extensive and stable expression pattern relative to the other legume bzip genes.
cluster ii included genes with variable and moderate expression.
cluster iii included genes with inconsistent  but low expression in tissues.
the broad expression pattern across various tissues indicated that members of the legume bzip transcription factor family are either expressed constitutively or in an organ-specific, development-dependent manner and may be involved in organ and tissue differentiation and seed developmental processes.
among the genes that were highly expressed during seed developmental stages, some were identified as homologs/orthologs of well-studied bzip genes from arabidopsis  <cit> , rice  <cit>  and maize  <cit> .
the atbzip39/abi <dig> gene, a homolog of mtbzip <dig>  pvbzip <dig> and pvbzip <dig> , is functionally involved in aba signaling and mediating embryogenesis in late embryo development  <cit> .
atbzip66/areb3/dpbf <dig> and atbzip67/dpbf <dig>  which were homologs of pvbzip <dig>  ljbzip <dig> and ljbzip <dig>  have been confirmed to play important roles in aba-mediated seed development, germination, and embryo maturation  <cit> .
in addition, gmbzip <dig> and mtbzip <dig> were homologous to the maize bzip factor opaque <dig>  which controls the transcription of a-zein, b- <dig> and b-prolamin genes and regulates protein accumulation, and amino acid and sugar metabolism in maize seeds .fig.
<dig> expression profiles of legume bzip genes.
a–d clustering of legume bzip genes according to their expression profiles in tissues including nodules, roots, stems, flowers, and pods and seeds at different developmental stages in gm , pv , mt  and lj .
the color scale represents log <dig> of the average signal values.
e, f the up- and down- regulated bzip genes identified in drought-stressed roots and shoots  and salt-stressed roots .
the color scale represents the fold change in the gene expression value compared with the control
as a kind of ubiquitous transcription factor, bzip proteins regulate the expression of a wide spectrum of stress-related genes.
we analyzed the expression values measured in drought-stressed roots and shoots corresponding to  <dig>   <dig>   <dig>   <dig>  and 14 days of drought, and in salt-stressed roots upon 180 mm nacl treatment.
we identified the obvious up- and down-regulated genes  and the log <dig>  ratio values were illustrated by heatmaps .
we detected five genes  up-regulated both in roots and shoots under drought stress, one gene  up-regulated only in roots, and three genes  up-regulated only in shoots.
among the genes down-regulated under drought stress, mtbzip18_i, mtbzip41_a and mtbzip34_a were down-regulated in both roots and shoots, mtbzip56_c was only down-regulated in roots, and mtbzip30_d was only down-regulated in shoots.
in roots under salt stress, mtbzip60_d, mtbzip32_i, mtbzip26_s, mtbzip2_s and mtbzip4_i were up-regulated and mtbzip54_s, mtbzip34_a and mtbzip28_s down-regulated.
most of the genes showing a response to drought and salt stress were concentrated in groups a and s. bzip genes from these groups have been frequently reported to be involved in sugar signaling and abiotic stress regulation  <cit> .
notably, the ckii and ca2+-dependent protein kinase phosphorylation site motifs  confined to group a in this study have been proposed to be involved in stress and/or aba signaling, which plays an important the role in the adaptation of plants to various abiotic environmental stress conditions like drought, high salinity, and cold stress  <cit> .
abi <dig>  and abfs /areb  have been shown to be key aba-dependent signal transduction factors involved in abiotic stress tolerance  <cit> .
among the stress-responsive legume bzip genes , mtbzip <dig> and mtbzip <dig> were homologous to abi <dig> and abfs/areb, respectively.
in addition, two genes, mtbzip <dig> and mtbzip <dig>  were responsive to salt stress and homologous to atbzip <dig>  which has been reported to be transcriptionally induced by salt treatment  <cit>  and leads to enhanced or reduced tolerance to salt stress when overexpressed or knocked out, respectively  <cit> .
overall, the expression analysis presented here improves our understanding of plant responses to stress at the molecular level and provides candidate legume bzip genes for future research.
CONCLUSIONS
bzip transcription factors have been extensively characterized in eukaryotic genomes and have been shown to play crucial roles in plant development, physiological processes, and biotic/abiotic stress responses.
using the six legume genomes available, we performed an extensive study of legume bzip genes including structure, phylogeny, sequence, and expression analyses.
the group classification of legume bzip genes based on their phylogenetic relationships was supported by subsequent analyses of gene structure, intron phases in the bzip domain, meme motif composition, dna-binding specificity and dimerization patterns, which showed group-specificity.
the group-specific sequence characteristics of the bzip members should have formed before the divergence of arabidopsis and legumes since conserved sequence characteristics were present in the same group containing both arabidopsis and legume bzip genes.
the global expression profile supports the role of legume bzip proteins in performing diverse developmental and physiological functions during tissue differentiation and seed development, as well as drought and salt stresses.
availability of supporting data
the data set supporting the results of this article is available in the dryad digital repository, http://dx.doi.org/ <dig> /dryad.m0qc <dig>
additional files
additional file 1:
the identified legume bzip proteins and their related information.
additional file 2:
position and pattern of introns within the basic and hinge regions of bzip domains of the legume bzip transcription factors.
additional file 3:
the phylogenetic analysis of bzip proteins based on the bzip protein sequences from arabidopsis  and six legumes  corresponding to fig.
1
.
additional file 4:
the map of intron-exon arrangement of legume bzip genes.
additional file 5:
histogram of intron number of legume bzip genes in each group.
additional file 6:
the detailed information of additional conserved motifs in bzip proteins as predicted by meme.
additional file 7:
alignment of basic and hinge regions of  <dig> legume bzip proteins.
additional file 8:
dna binding specificity prediction of each group.
additional file 9:
amino acid sequence alignment of the leucine zipper region of  <dig> legume bzip proteins.
additional file 10:
chromosomal distributions of legume bzip genes.
additional file 11:
bzip genes present on duplicated chromosomal segments.
additional file 12:
expression values of legume bzips in multiple tissues.
abbreviations
wgdwhole genome duplication
myamillion years ago
hmmhidden markov model
mlmaximum likelihood
orfopen reading frame
mememultiple expectation maximization for motif elicitation
abi5aba-insensitive 5
abreaba response element
abfabre binding factor
arebaba-responsive element binding protein
competing interests
the authors declare that they have no competing interests.
authors’ contributions
zw, yl and bl conceived and designed the experiments, zw, kc, lw, ly, hj, sl, yl and bl analyzed the data, zw, yl and bl wrote the paper.
all authors have read and approved the final version of the manuscript.

 a defining feature of transcription factors is that they contain one or more sequence-specific dna-binding domains that bind to the promoter and/or enhancer regions of target genes to regulate gene expression. the basic leucine zipper transcription factor family, one of the most diverse transcription factors families, is characterized by a highly conserved bzip domain which is 60– amino acids in length and composed of two parts: a basic region and a leu zipper. the basic and leu zipper regions are structurally and functionally distinct. the basic region comprises approximately amino acid residues with the invariant motif n-x7-r/k-x and is responsible for nuclear localization and dna binding, whereas the leu zipper is composed of heptad repeats of leu or other bulky hydrophobic amino acids positioned exactly nine amino acids towards the c-terminus and mediates homo- and/or heterodimerization of bzip proteins. bzip transcription factor encoding genes have been identified extensively in plants including arabidopsis, rice, sorghum, maize, grapevine, cucumber, castor bean and barley with the availability of their whole genome sequences. like other transcription factors, members of the bzip transcription factor family are expressed constitutively or in an organ-specific, stimulus-responsive, development-dependent or cell cycle-specific manner in plants. it has been reported that bzip transcription factors are involved in various biological processes including organ and tissue differentiation, embryogenesis, seed maturation and storage protein gene regulation, floral transition and initiation, and vascular development. moreover, bzip transcription factors are also regarded as important regulators in signaling and responses to abiotic/biotic stimuli, including abscisic acid signaling, hypoxia, drought, high salinity, cold stress, hormone and sugar signaling, light responses, osmotic stresses and pathogen defense. in soybean, three bzip genes were found to function as negative regulators of aba signaling and confer salt and freezing tolerance in transgenic arabidopsis. in recent years, legume genome sequencing projects have been initiated and completed in lotus japonicus, soybean, medicago truncatula, pigeonpea, chickpea and common bean. these six legumes belong to the large papilionoideae subfamily and fall within two sub-clades of the papilionoid legumes: the phaseoloids and hologalegina. the phaseoloids are mostly tropical and include the crops soybean, pigeonpea and common bean, while the hologalegina are mostly temperate and include medicago, chickpea and lotus. among the widespread genome duplications throughout the history of flowering plants, two recent whole genome duplication events have affected the evolution of legume genomes. the older polyploidy event, shared by all legumes, is estimated to have occurred 56– million years ago, while the more recent genome duplication event occurred up to mya only in the lineage leading to glycine. genome duplication and subsequent fractionation have played key roles in shaping present-day legume genomes and also the sizes of gene families. with the availability of these legume genome sequences, the members of the bzip transcription factor family were systematically investigated and analyzed in this study. we identified all legume bzip genes and analyzed their bzip domain sequences, gene structure and additional meme motifs, which was in agreement with and supported the phylogenetic classification. then, we predicted the dna-binding-site specificity and dimerization properties of the legume bzip proteins. we also investigated the impact of the two legume-lineage wgds and tandem duplication on the expansion of the legume bzip gene family. by analyzing their expression profiles, legume bzip genes constitutively or specifically expressed in different tissues and seed developmental stages were identified, as well as candidate legume bzips responsive to drought and salt stresses. methods identification of bzip transcription factors in six legume genomes all genomic sequences and annotated proteins of the six legumes were downloaded from ftp://ftp.jgi-psf.org/pub/compgen/phytozome/v /gmax/, http://jcvi.org/medicago/display.php?pagename=general&section=download, http://genome.jgi.doe.gov/pages/dynamicorganismdownload.jsf? organism = phytozomev, http://cicar.comparative-legumes.org/, http://www.icrisat.org/gt-bt/iipg/genome_manuscript.html and ftp://ftp.kazusa.or.jp/pub/lotus/lotus_r /. to identify all the possible bzip proteins in the six legume genomes, both local blast and hidden markov model searches were performed. for blastp, the known bzip proteins from arabidposis, rice and maize were used as queries and the e-value was set to 1e- for the hmm search, the profile of bzip domain was used and the e-value threshold was set at the sequences were further analyzed to confirm the presence and integrity of the bzip domain through the expasy proteomics server and interpro. all bzip domain sequences were aligned using mafft to manually check and remove sequences with incomplete domains. the nomenclature was based on the exact positions of the bzip genes on the chromosomes/scaffolds from top to bottom. distinct transcripts encoded by the same gene locus shared the same gene number with an additional decimal part, such as point or. sequence alignment and phylogenetic analysis the bzip amino acid sequences from a. thaliana and the six legume genomes were aligned using clustalx with gap opening and gap extension penalties of and, respectively. the phylogenetic tree was reconstructed by the maximum likelihood method using the phyml software. jtt + g was selected as the best model for constructing the phylogenetic tree by the akaike information criterion implemented in prottest. bootstrap values from replicates were indicated at each node. mega was used to show the tree. structure of bzip genes the positional information of both the gene sequence and the corresponding coding sequence were loaded into the gene structure display server v to obtain information on the intron/exon structure. the coordinates of the bzip domain in each protein were recalculated into the coordinates in gene sequence and featured in gene structure. we used genewise to analyze the intron distribution pattern and intron splicing phase within the basic and hinge regions of the bzip domains in the six legumes. detection of additional conserved motifs to identify additional conserved motifs outside the bzip domain of legume bzip transcription factors, we used the multiple em for the motif elicitation tool. the limits for maximum width, minimum width and maximum number of motifs were specified as and respectively. fifty motifs were finally confirmed because of their low e-values. the motifs were numbered according to the order displayed in meme and were considered as group-specific signatures for their presence of high frequency in the given groups. detection of duplicated genes and estimation of synonymous and nonsynonymous substitutions per site and their ratio the duplicated gene pairs derived from segmental duplication were identified in the legume genomes based on the method from the plant genome duplication database. an all-against-all blastp comparison provided the gene pairs for syntenic clustering determined by mcscan. tandem duplication arrays were identified using blastp with a threshold of e < 10− and one unrelated gene among cluster members was tolerated, as described in a. thaliana. pairs from segmental or tandem duplications were used to estimate ka, ks and their ratio. amino acid sequences from segmentally or tandemly duplicated pairs were first aligned and then guided and transferred into a cdna sequences alignment using in-house perl scripts. then, the software kaks_calculator was used to compute ka and ks values for each pair following the yn model. expression analysis of legume bzip genes for different tissues/organ and seed developmental stages, the normalized counts for bzip genes from rna-seq were obtained from the soyseq and pvgea databases. microarray expression values for medicago and lotus were downloaded from http://mtgea.noble.org/v3/ and http://ljgea.noble.org/ v2/. to identify candidate bzip genes responsive to drought and/or salt stresses, microarray gene expression data in medicago were downloaded from zhang et al. and li et al., respectively. the gene expression changes in drought-stressed roots and shoots were obtained by comparing with levels in the watered control. similarly, the fold changes in gene expression were calculated when comparing the salt-stressed roots and control. the corresponding relationships between microarray probes and legume bzip genes were built using blast. the expression values or normalized counts were log10-transformed and the gplots package was used to make a heatmap in r.  identification and nomenclature of the legume bzip transcription factor family through sequence similarity and domain searches, gm, mt, pv, ca, cc and lj bzip genes, encoding and distinct proteins, respectively, were identified in the six legume genomes. based on their exact positions on chromosomes/scaffolds, we gave a unique name to each bzip protein. the related information on bzip transcription factors are listed in the additional file domain analysis showed that all of the bzip transcription factors except eight had a typical bzip domain with an invariant n-x7-r/k motif in the basic region and a heptad repeat of leu or other bulky hydrophobic amino acids positioned exactly nine amino acids upstream of r/k toward the c-terminus. of the remaining eight, in glyma12g and phvul.011g, the conserved asn in the basic region was replaced by lys. in glyma03g, glyma19g, medtr7g and ca_, the conserved arg/lys in the basic region was substituted by an ile, whereas in glyma11g, the conserved arg/lys in the basic region was substituted by a trp. in c. cajan, the heptad repeat of leu was positioned at amino acids toward the c terminus instead of the usual nine. all of these unusual changes in the bzip domain have been observed in rice bzip sequences previously. phylogenetic analysis and classification of legume bzip genes to investigate the phylogenetic relationships of the bzip transcription factors, protein sequences from the six legumes and protein sequences from arabidopsis were analyzed. in accordance with the bzip classification in arabidopsis, the phylogenetic tree was subdivided into clades with well-supported bootstrap values. all groups contained legume bzip proteins clustered together with atbzip proteins in the same clade, except group f, which included members from only three genomes: g. max, p. vulgaris and c. arietinum. in addition, legume bzip proteins and atbzip genes formed two small unique and several other individual clades, which were classified into the unclassified group based on their possible independent evolutionary trajectories from other clades. group classification was supported by the group-specific sequence characteristics identified in the following analyses of gene structures, intron phases, additional conserved motifs outside the bzip domain and dna-binding site specificity in each group. it is evident that the group-specific sequence characteristics of the bzip members formed before the divergence of arabidopsis and legumes since conserved sequence characteristics were present in the same group containing both arabidopsis and legumes bzip genes. nevertheless, it seems that intra-species duplication and parallel evolution of the bzip family in each legume has occurred afterward and contributed to the member variation in each group.fig. phylogenetic tree of legume and arabidopsis bzip genes. bzip protein sequences were aligned by clustal x and the phylogenetic tree was constructed using phyml by the maximum likelihood method. bootstrap values are based on replicates. genes on branch ends from different species are denoted by different colored circles. the legume bzip proteins were grouped into distinct clades, which are indicated by colored branches gene structure of legume bzip genes the intron-exon organization can reflect the evolutionary trajectory of gene families. we examined the gene structures of all legume bzips and found that the structural patterns were similar among members within each group but distinct between different groups. the number of introns in each group was uneven but relatively concentrated. we detected gmbzip genes, mtbzip genes, pvbzip genes, cabzip genes, ccbzip genes and ljbzip genes with no introns. most of these intronless genes were clustered into groups s and f. among the intron-containing bzip genes, the number of introns within the open reading frame varied from to in gmbzip, in mtbzip, in pvbzip, in cabzip, in ccbzip, and in ljbzip, close to the highest number of introns within the orf reported in arabidopsis, rice, sorghum, maize, castor bean, barely and cucumber. the bzip genes with the most introns were commonly found in groups d and g. the intron positions within the orf were diverse. the phases of the splicing sites within the orfs also differed, but the positions and phases of introns in the basic and hinge regions of the bzip domain were highly conserved. the legume bzip genes showed four intron patterns based on the intron positions, presence or number and splicing phases within the basic and hinge regions. pattern a, having one intron in phase within the hinge region at the − position, was seen in all members of groups a and g in the legumes. pattern b, having two introns each in phase 0: one in the basic region at the − position and the other in the hinge region at the − position, was seen in all members of group d in the legumes. pattern c, having a single intron in phase at the − position in the basic region, was seen in all members of groups c, e and h in the legumes. pattern c also was seen in all members of group i in g. max, m. truncatula, c. arietinum and c. cajan, and most members of group i in l. japonicus. pattern d was no introns in the basic and hinge regions and was seen in groups b, s and f, which contained most of the intronless genes. among genes showing pattern d, were intronless, while the remaining had introns outside the basic and hinge regions. in summary, the overall pattern of intron position acts as an index for the group classification and phylogenetic relationships in the legume bzip gene family. the splicing phase has remained well conserved during the course of evolution in legume bzip genes.fig. intron patterns within the basic and hinge regions of the legume bzip domain. the primary structure of the bzip domain is shown at the top of the picture. p indicates the intron splicing site is between codons, and p means the intron splicing site is located between the second and third nucleotides of one codon. based on the intron incidence and positions, as well as the splicing phase, the legume bzip genes were divided into four patterns. details of the intron positions within the bzip domain in the legume bzip proteins are shown in additional file 2 identification of additional structural features in the legume bzip genes all legume bzip protein sequences were loaded into the meme analysis tool and a total of additional conserved motifs outside the bzip domain were identified. the multi-level consensus sequences and the amino acid lengths of these conserved motifs are given in additional file the legume bzip proteins within the same group had similar motif compositions, suggesting conserved evolution and supporting the group classification. additionally, some motifs were shared by different groups, such as motifs and in two groups, motif in three groups, motifs and in four groups, and motifs and in five groups. nevertheless, most of the conserved motifs appeared specific to each group and therefore the group-specific motifs could help determine the specific functions of members in each group. notably, we did not detect any additional conserved motifs outside the bzip domain in groups b or f.fig. distribution of additional conserved motifs identified by meme. motif compositions based on the position of the bzip domain and additional conserved motifs outside the bzip domain in representative legume bzip proteins for each group are shown. the bzip domains are shown in pink and different motifs are highlighted in different colored boxes with numbers to the motifs shared by different groups are marked with black stars. details for the predicted conserved motifs are given in additional file 6 a few of these motifs in legume bzips have been studied for their possible biological functions. a part of motifs and represent potential casein kinase ii phosphorylation sites, indicated by the motif patterns e, tlge and tvde. motif also contained a phosphorylation site for the ca2+-dependent protein kinase, presented as rqs. motif in group h was group-specific, and contained potential ck ii phosphorylation sites, indicated by cye. such motifs were identified in group a exclusively. in addition, all members in group d shared motif which was a dog domain, and motif the function of which was unclear. motifs and were observed in group g exclusively and were characterized by a part of the proline-rich domain, which has been shown to have transcriptional activation potential. interestingly, there were some common motifs among the six legumes, maize and rice. for example, motif in group a was the same as motif in group a in maize. motifs and in group d are in common with motifs and in group d of maize and motifs and of rice, respectively. prediction of dna-binding-site specificity in legume bzip proteins experiments of mutant proteins demonstrated that the bzip tfs binding specificity is independently determined by the core basic region and the hinge region, and the two regions have an additive effect on dna-binding specificity. to predict the dna-binding-site specificity of the legume bzip proteins, the amino acid sequences of the basic and hinge regions of legume bzip proteins were aligned, revealing some highly conserved amino acid residues within each group. we can predict the dna-binding specificity in a group manner, as described in additional file furthermore, the amino acids were numbered as previously reported, and the first leu in the leu heptad repeats was numbered + for the two invariant sites asparagine and arginine, numbered − and − respectively, new dna-binding specificities will occur if other amino acids functionally replace these two invariant sites. nevertheless, such replacements were infrequent and occurred only in groups g and u in the legume bzip proteins. at position − the conserved asparagine was replaced with lysine in gmbzip and pvbzip in group g. the same replacement was also observed in four bzips in maize, two bzips in barley, and two bzips in castor bean. at position − nine members in group u had a hydrophobic isoleucine residue instead of arginine or lysine. an identical substitution pattern was observed in two bzips in maize, two in barley, one bzip in grapevine, and one bzip in castor bean. it was demonstrated that an arginine to isoleucine mutation in the basic domain of the yeast bzip factor gcn completely inhibited its affinity for the ap site. additionally, it was reported that oszip-2a belonging to group u in rice, because of this replacement, does not recognize g-boxes. this evidence suggests that unusual substitutions in the dna-binding domain affect the dna-binding specificity. these predictions were made to facilitate further studies on the dna-binding patterns of the legume bzip transcription factors. prediction of dimerization properties in legume bzip proteins studies have demonstrated that the leu zipper region of the bzip domain, arranged in the form of heptad repeats, mediates homo- and/or heterodimerization between the parallel coiled-coil structures. within each heptad, the amino acid positions are recognized in order as g, a, b, c, d, e, and f. leu zipper oligomerization, dimerization stability and specificity are determined mainly by the four amino acids present at the a, d, e and g positions because of their special positions near the leu zipper interface. the a and d residues are typically hydrophobic on the surface of the helix and create a hydrophobic core that promotes the interaction between two monomers and is essential for dimer stability. the a position contains asparagine, which can form a polar pocket in the hydrophobic interface that limits oligomerization in interhelical interactions and produce more stable n–n interactions at a↔aʹ than other amino acids. the conserved leu at the d position, one of the most stabilizing aliphatic amino acids, is important to maintain dimer stability. however, the e and g positions that flank the dimerization interface frequently contain charged amino acids including the acidic amino acids glutamic acid and aspartic acid, and the basic amino acids arginine and lysine, which are thought to form salt bridges between helices in electrostatic interactions. to predict the dimerization properties of the legume bzip transcription factors, a detailed analysis was carried out to characterize the amino acids present at the a, d, e and g positions. figure 4a shows the composition of different kinds of amino acids found in the a, d, e and g positions in the six legume bzips, respectively. at the a position, about 25 % of the residues were asparagine, suggesting that there will be a greater number of homodimerizing leu zippers through stable n–n interhelical interaction at the a↔aʹ position among legume bzip proteins. the frequency of asparagine residues in the a position was highest in the second heptad followed by the fifth heptad, accounting for  % and  %, respectively, similar to earlier observations for atbzip proteins. at the d position, the frequency of leu, responsible for dimer stability, was about 68 % in the legume bzips, which is a little less than in osbzips and zmbzips, but significantly greater than in atbzips. at the e and g positions, the frequencies of charged amino acids were and 59 %, respectively.fig. prediction of dimerization properties of the legume bzip proteins. a pie charts presenting the frequency of amino acids in all the a, d, e and g positions of the leu zipper in legume bzip domains. b histogram of the frequency of asn in the a positions of the leu zippers among legume bzip proteins. c histogram of the frequency of attractive or repulsive g↔eʹ pairs per heptad for legume bzip proteins in electrostatic interactions, the amino acids at position g and oppositely charged amino acids at position eʹ can form interhelical salt bridges that determine the dimerization specificity and stability. the dimerization specificity of human, drosophila, arabidopsis, rice, maize, barley and castor bean bzip tfs has been predicted on the basis of attractive or repulsive interhelical g↔eʹ electrostatic interactions. to analyze the contribution of charged residues to the dimerization properties of the legume bzip proteins, we calculated the frequency of attractive and repulsive g↔eʹ pairs in each heptad of the bzip leu zippers; the corresponding histograms are shown in fig. 4c. if both the g and corresponding e position amino acids are charged, they are referred to as complete g↔eʹ pairs. this analysis was carried out on the basis of the frequency of attractive and repulsive g↔eʹ pairs, which were classified into four groups, attractive basic-acidic pairs, attractive acidic-basic pairs, repulsive basic pairs and repulsive acidic pairs in the heptads. attractive g↔eʹ pairs were predominant in the second, fifth and sixth heptads, thereby suggesting the chances of heterodimerization. in contrast, few complete g↔eʹ pairs were observed in the eighth heptads except two attractive acidic-basic pairs from gmbzip and pvbzip. moreover, only attractive acidic-basic pairs were present in ninth heptads, which was similar to observations in hvbzips, osbzips and zmbzips. we divided the legume bzip proteins into subfamilies on the basis of the defining properties of dimerization specificity. these subfamilies were divided into three general groups: those that strongly favor homodimerization within eight subfamilies, those with both homo- and heterodimerization properties and those that strongly favor heterodimerization within five subfamilies. the results indicated the complexity and diversity of dimerization patterns in legume bzip proteins, with the potential to homodimerize with themselves or with members in the same subfamily as well as heterodimerize with other subfamily members, which has been shown in maize, rice, and arabidopsis. based on the criteria used to define the boundaries and natural c-terminus, we observed that the length of the leu zipper in the bzip transcription factor family was variable, ranging from two to nine heptads. among the bzip proteins,  % had only two short zippers and belonged to bz and group g, about 28 % had only three short zippers, and about 9 % had no α-helix breakers for or more heptads, were mainly distributed in bz37–bz and included members of group i. the impact of whole genome duplication and tandem duplication on the expansion of the legume bzip gene family to explore the expansion mechanism in detail, we analyzed the contribution of wgd –derived segmental duplication and tandem duplication to the expansion of the legume bzip gene family. first, a phylogenetic tree was constructed based on the concentrated orthologous sequences of rbc and matk genes from the six legume genomes, which supported the classification of two sub-clades of papilionoideae legumes: phaseoloids and hologalegina. we carried out a genome-wide identification of collinear duplicated blocks derived from segmental duplication in each species, and then examined the pairwise synonymous distances of paralogs within duplicated collinear blocks. by plotting their distribution, two distinct peaks were found: one was specific to gm-gm paralogues, while the other was observed in all six legumes.fig. whole genome duplication -derived legume bzips genes. a the phylogenetic relationships of the six legumes based on the concentrated orthologous sequences of rbc and matk genes. two wgd events, the glycine-lineage-specific and early-legume wgd events, are indicated on the corresponding clades. b the ks distribution of paralogs derived from wgd-derived duplicated genomic blocks in each legume. two obvious peaks correspond to the old and recent legume-lineage wgd events. c, d the black lines in the ideogram show the chromosomal positions of all identified bzip genes, and the duplicated bzip pairs are linked by lines in soybean and common bean. e patterns of gene retention or loss indicated by two different informative tree topologies using common bean and medicago as outgroups, respectively. the numbers for each pattern are shown in parentheses. f the ks and ka distribution for each duplicated legume bzip gene pair. the red dashed line means the ka/ks ratio is equal to particularly, eight gm duplicated pairs and one mt pair, whose ka/ks values were greater than, are indicated in the inner figure. the red dashed line in inner figure means the ka/ks ratio is equal to. the numbers correspond to duplicated pairs: 1: gmbzip23-gmbzip131; 2: gmbzip49-gmbzip109; 3: gmbzip55-gmbzip107; 4: gmbzip8-gmbzip113; 5: gmbzip4-gmbzip69; 6: gmbzip26-gmbzip42; 7: gmbzip87-gmbzip101; 8: gmbzip65-gmbzip113; 9: mtbzip2-mtbzip26 second, the chromosomal distribution of the legume bzip genes was plotted and the bzip gene pairs on duplicated chromosomal collinear segments were connected by lines. among these genes, some were segmentally duplicated once and some were duplicated twice or thrice. the duplication occurred within a chromosome or between chromosomes. we detected gm, mt, pv, ca, cc and lj bzip genes involved in segmental duplication, accounting for around  %,  %,  %,  %,  %, and  % of the bzip genes in each species. the higher ratio in soybean reflects the preferential gene retention after multiple rounds of wgd, while the different ratios among the other five legumes may be mainly attributable to the genome assembly quality or species-specific evolution in each genome. moreover, we roughly identified different wgd event origins for duplicated bzip gene pairs according to their pairwise synonymous distances using criteria from soybean: ks values of – correspond to the 13-mya glycine-lineage-specific wgd, and ks values of – in soybean and 0– in the other five legumes correspond to the 59-mya early-legume wgd, while larger ks values correspond to more ancient wgd events like the ‘gamma’ event. by ordering the ks values, gm bzip gene pairs were associated with the mya glycine-lineage-specific wgd and gm, mt, pv, ca, cc, and lj bzip gene pairs were associated with the mya early-legume wgd; the others were associated with more ancient wgds. because only soybean has undergone the glycine-lineage-specific wgd, the other five legume genomes could be considered putative ancestors for investigating gene retention and loss after the recent genome duplication in soybean. mt and pv were chosen, because of their good genome assembly and annotation, to construct phylogenetic trees with gm genes. we identified informative tree models of gene retention: and, and gene loss and, showing that 83 % and 87 % ancestral loci were retained as paralogs in soybean after the recent wgd. compared with the average  % of genes retained as paralogs at the whole genome level, genes of the bzip family were highly retained. to better understand the evolutionary constraints acting on the legume bzip genes, we calculated the ka/ks ratios for all legume duplicated bzip gene pairs. the resulting pairwise comparison data showed that the ka/ks values of only eight gm duplicated pairs and one mt pair were larger than while all of the remaining ka/ks ratios were less than, suggesting that the bzip family has mainly undergone strong purifying selection and the legume bzip genes are slowly evolving at the protein level. we further compared the strength of selection on bzip paralogs between the recent and old whole-genome duplications in soybean. the average ka/ks ratio for the recently duplicated bzip genes was higher than that of the early duplicated bzip genes, and there was a significant difference between these ratios. this indicated that the younger bzip proteins may be under stronger evolutionary constraints than older proteins, but supported the notion that the legume bzip gene family is essential for the regulation of cellular processes. third, we determined members with tandem duplication in each legume. we detected only gm, mt, pv, ca, cc and lj tandem gene pairs, indicating the limited contribution of tandem duplication to the expansion of the gene family. therefore, it seems that segmental duplication rather than tandem duplication is the major mechanism driving the expansion of this gene family. expression analysis of bzip transcription factors expression data from different tissues and seed developmental stages in soybean, medicago, common bean and lotus were downloaded, and hierarchical clustering was performed to visualize a global transcription profile of the legume bzip genes. as shown in fig. , the heatmaps were always divided into three recognized clusters, which was similar to the results in rice, maize and other plants. the different clusters corresponded to overall differences in expression patterns including expression values and specificity across tissues. the latter was indicated by cv values, which were calculated for each gene, and helped to recognize genes expressed in specific tissues or stages. cluster i included genes with relatively high expression levels and the least expression variability, indicating an extensive and stable expression pattern relative to the other legume bzip genes. cluster ii included genes with variable and moderate expression. cluster iii included genes with inconsistent but low expression in tissues. the broad expression pattern across various tissues indicated that members of the legume bzip transcription factor family are either expressed constitutively or in an organ-specific, development-dependent manner and may be involved in organ and tissue differentiation and seed developmental processes. among the genes that were highly expressed during seed developmental stages, some were identified as homologs/orthologs of well-studied bzip genes from arabidopsis, rice and maize. the atbzip39/abi gene, a homolog of mtbzip pvbzip and pvbzip, is functionally involved in aba signaling and mediating embryogenesis in late embryo development. atbzip66/areb3/dpbf and atbzip67/dpbf which were homologs of pvbzip ljbzip and ljbzip have been confirmed to play important roles in aba-mediated seed development, germination, and embryo maturation. in addition, gmbzip and mtbzip were homologous to the maize bzip factor opaque which controls the transcription of a-zein, b- and b-prolamin genes and regulates protein accumulation, and amino acid and sugar metabolism in maize seeds.fig. expression profiles of legume bzip genes. a–d clustering of legume bzip genes according to their expression profiles in tissues including nodules, roots, stems, flowers, and pods and seeds at different developmental stages in gm, pv, mt and lj. the color scale represents log of the average signal values. e, f the up- and down- regulated bzip genes identified in drought-stressed roots and shoots and salt-stressed roots. the color scale represents the fold change in the gene expression value compared with the control as a kind of ubiquitous transcription factor, bzip proteins regulate the expression of a wide spectrum of stress-related genes. we analyzed the expression values measured in drought-stressed roots and shoots corresponding to and 14 days of drought, and in salt-stressed roots upon 180 mm nacl treatment. we identified the obvious up- and down-regulated genes and the log ratio values were illustrated by heatmaps. we detected five genes up-regulated both in roots and shoots under drought stress, one gene up-regulated only in roots, and three genes up-regulated only in shoots. among the genes down-regulated under drought stress, mtbzip18_i, mtbzip41_a and mtbzip34_a were down-regulated in both roots and shoots, mtbzip56_c was only down-regulated in roots, and mtbzip30_d was only down-regulated in shoots. in roots under salt stress, mtbzip60_d, mtbzip32_i, mtbzip26_s, mtbzip2_s and mtbzip4_i were up-regulated and mtbzip54_s, mtbzip34_a and mtbzip28_s down-regulated. most of the genes showing a response to drought and salt stress were concentrated in groups a and s. bzip genes from these groups have been frequently reported to be involved in sugar signaling and abiotic stress regulation. notably, the ckii and ca2+-dependent protein kinase phosphorylation site motifs confined to group a in this study have been proposed to be involved in stress and/or aba signaling, which plays an important the role in the adaptation of plants to various abiotic environmental stress conditions like drought, high salinity, and cold stress. abi and abfs /areb have been shown to be key aba-dependent signal transduction factors involved in abiotic stress tolerance. among the stress-responsive legume bzip genes, mtbzip and mtbzip were homologous to abi and abfs/areb, respectively. in addition, two genes, mtbzip and mtbzip were responsive to salt stress and homologous to atbzip which has been reported to be transcriptionally induced by salt treatment and leads to enhanced or reduced tolerance to salt stress when overexpressed or knocked out, respectively. overall, the expression analysis presented here improves our understanding of plant responses to stress at the molecular level and provides candidate legume bzip genes for future research.  bzip transcription factors have been extensively characterized in eukaryotic genomes and have been shown to play crucial roles in plant development, physiological processes, and biotic/abiotic stress responses. using the six legume genomes available, we performed an extensive study of legume bzip genes including structure, phylogeny, sequence, and expression analyses. the group classification of legume bzip genes based on their phylogenetic relationships was supported by subsequent analyses of gene structure, intron phases in the bzip domain, meme motif composition, dna-binding specificity and dimerization patterns, which showed group-specificity. the group-specific sequence characteristics of the bzip members should have formed before the divergence of arabidopsis and legumes since conserved sequence characteristics were present in the same group containing both arabidopsis and legume bzip genes. the global expression profile supports the role of legume bzip proteins in performing diverse developmental and physiological functions during tissue differentiation and seed development, as well as drought and salt stresses. availability of supporting data the data set supporting the results of this article is available in the dryad digital repository, http://dx.doi.org/ /dryad.m0qc additional files additional file 1: the identified legume bzip proteins and their related information. additional file 2: position and pattern of introns within the basic and hinge regions of bzip domains of the legume bzip transcription factors. additional file 3: the phylogenetic analysis of bzip proteins based on the bzip protein sequences from arabidopsis and six legumes corresponding to fig. 1. additional file 4: the map of intron-exon arrangement of legume bzip genes. additional file 5: histogram of intron number of legume bzip genes in each group. additional file 6: the detailed information of additional conserved motifs in bzip proteins as predicted by meme. additional file 7: alignment of basic and hinge regions of legume bzip proteins. additional file 8: dna binding specificity prediction of each group. additional file 9: amino acid sequence alignment of the leucine zipper region of legume bzip proteins. additional file 10: chromosomal distributions of legume bzip genes. additional file 11: bzip genes present on duplicated chromosomal segments. additional file 12: expression values of legume bzips in multiple tissues. abbreviations wgdwhole genome duplication myamillion years ago hmmhidden markov model mlmaximum likelihood orfopen reading frame mememultiple expectation maximization for motif elicitation abi5aba-insensitive 5 abreaba response element abfabre binding factor arebaba-responsive element binding protein competing interests the authors declare that they have no competing interests. authors’ contributions zw, yl and bl conceived and designed the experiments, zw, kc, lw, ly, hj, sl, yl and bl analyzed the data, zw, yl and bl wrote the paper. all authors have read and approved the final version of the manuscript. 
plant bzip proteins characteristically harbor a highly conserved bzip domain with two structural features: a dna-binding basic region and a leucine zipper dimerization region. four kinds of intron patterns within the basic and hinge regions were defined and additional conserved motifs were identified, both presenting high group specificity and supporting the group classification. we also detected several key legume bzip genes involved in drought- and salt-responses by comparing fold changes of expression values in drought-stressed or salt-stressed roots and leaves. the chromosome distribution and analysis for wgd-derived duplicated blocks revealed that the legume bzip genes have expanded mainly by segmental duplication rather than tandem duplication. the legume bzip genes were categorized into groups according to their phylogenetic relationships with genes from arabidopsis. expression data further revealed that the legume bzip genes were expressed constitutively or in an organ-specific, development-dependent manner playing roles in multiple seed developmental stages and tissues. in summary, this genome-wide identification, characterization and expression analysis of legume bzip genes provides valuable information for understanding the molecular functions and evolution of the legume bzip transcription factor family, and highlights potential legume bzip genes involved in regulating tissue development and abiotic stress responses. we predicted the dna-binding patterns and the dimerization properties, based on the characteristic features in the basic and hinge regions and the leu zipper, respectively, which indicated that some highly conserved amino acid residues existed across each major group. despite the availability of six completely sequenced legume genomes, a comprehensive investigation of bzip family members in legumes has yet to be presented. bzip gene familylegume genomesevolutionexpression analysishttp://dx.doi.org/ /501100001809national natural science foundation of china 31461143022liao boshou http://dx.doi.org/ /501100001809national natural science foundation of china 31301256liao boshou http://dx.doi.org/ /501100001809national natural science foundation of china 31071456lei yong http://dx.doi.org/ /501100001809national natural science foundation of china 31471536liu shengyi national high technology research and development program of china 2013aa102602wang zhihui national basic research and development program 2011cb109304jiang huifang national basic research and development program 2011cb109305liu shengyi issue-copyright-statement© the author 2015 they have been shown to be diverse transcriptional regulators, playing crucial roles in plant development, physiological processes, and biotic/abiotic stress responses.
 plant bzip proteins characteristically harbor a highly conserved bzip domain with two structural features: a dna-binding basic region and a leucine zipper dimerization region. they have been shown to be diverse transcriptional regulators, playing crucial roles in plant development, physiological processes, and biotic/abiotic stress responses. despite the availability of six completely sequenced legume genomes, a comprehensive investigation of bzip family members in legumes has yet to be presented.  in this study, we identified bzip genes encoding distinct proteins in six legumes, glycine max, medicago truncatula, phaseolus vulgaris, cicer arietinum, cajanus cajan, and lotus japonicus. the legume bzip genes were categorized into groups according to their phylogenetic relationships with genes from arabidopsis. four kinds of intron patterns within the basic and hinge regions were defined and additional conserved motifs were identified, both presenting high group specificity and supporting the group classification. we predicted the dna-binding patterns and the dimerization properties, based on the characteristic features in the basic and hinge regions and the leu zipper, respectively, which indicated that some highly conserved amino acid residues existed across each major group. the chromosome distribution and analysis for wgd-derived duplicated blocks revealed that the legume bzip genes have expanded mainly by segmental duplication rather than tandem duplication. expression data further revealed that the legume bzip genes were expressed constitutively or in an organ-specific, development-dependent manner playing roles in multiple seed developmental stages and tissues. we also detected several key legume bzip genes involved in drought- and salt-responses by comparing fold changes of expression values in drought-stressed or salt-stressed roots and leaves.  in summary, this genome-wide identification, characterization and expression analysis of legume bzip genes provides valuable information for understanding the molecular functions and evolution of the legume bzip transcription factor family, and highlights potential legume bzip genes involved in regulating tissue development and abiotic stress responses. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords bzip gene familylegume genomesevolutionexpression analysishttp://dx.doi.org/ /501100001809national natural science foundation of china 31461143022liao boshou http://dx.doi.org/ /501100001809national natural science foundation of china 31301256liao boshou http://dx.doi.org/ /501100001809national natural science foundation of china 31071456lei yong http://dx.doi.org/ /501100001809national natural science foundation of china 31471536liu shengyi national high technology research and development program of china 2013aa102602wang zhihui national basic research and development program 2011cb109304jiang huifang national basic research and development program 2011cb109305liu shengyi issue-copyright-statement© the author 2015 