Introduction

Transcription factors (TFs), a group of proteins that regulate expression of the target genes by binding specifically to its promoter, are involved in multiple biological processes. The bZIP transcription factor family is one of the largest and most diverse families in plants. The bZIP transcription factors have a conserved bZIP domain containing 60–80 residues with a basic region and a leucine zipper region (Nijhawan et al. 2008). The basic region, consisting of an invariant motif N-X7-R/K at the N terminus, is responsible for nuclear localization and sequence-specific DNA binding. The less-conserved leucine zipper region is composed of several heptad repeats of leucine (Leu) or other bulky hydrophobic amino acids toward the C terminus (Nijhawan et al. 2008). Through the interaction of the hydrophobic amino acids in the long α-helices of leucine zipper regions, the two subunits are tightly bound to each other to form a coiled-coil homo- and/or heterodimerization structure. Aside from the bZIP domain, additional motifs may contribute to the functional variety of the bZIP proteins. For example, two conserved motifs, R/KxxS/T and S/TxxD/E, were identified as potential phosphorylation sites for Ca2+-dependent protein kinase and casein kinase II (CKII), respectively (Jakoby et al. 2002; Nijhawan et al. 2008). Proline-rich, glutamine-rich and acidic domains have also been confirmed as transcriptional activators in bZIP proteins (Jakoby et al. 2002; Nijhawan et al. 2008).

To date, the bZIP gene family has been identified or predicted in several plant species, such as Arabidopsis thaliana (Jakoby et al. 2002), rice (Oryza sativa) (Nijhawan et al. 2008), maize (Zea mays) (Wei et al. 2012), tomato (Solanum lycopersicum L.) (Li et al. 2015), sorghum (Sorghum bicolor) (Wang et al. 2011), sesame (Sesamum indicum L.) (Wang et al. 2018b), barley (Hordeum vulgare L.) (Pourabed et al. 2015), soybean (Glycine max) (Zhang et al. 2018), apple (Malus domestica) (Zhao et al. 2016), Chinese cabbage (Brassica rapa) (Bai et al. 2016), rapeseed (Brassica napus) (Zhou et al. 2017), castor bean (Ricinus communis L.) (Jin et al. 2014), cucumber (Cucumis sativus L.) (Baloglu et al. 2014) and grapevine (Vitis vinifera) (Liu et al. 2014b). Previous studies had shown that bZIP genes are involved in regulating seed maturation and germination (Toh et al. 2012), somatic embryogenesis (Guan et al. 2009), cell elongation (Fukazawa et al. 2000) and photomorphogenesis (Huang et al. 2012). In addition to the vital roles of bZIPs in plant growth and development, they are also key components in response to stresses and hormone signaling, such as drought (Yoshida et al. 2010; Chen et al. 2012), high salinity (Hsieh et al. 2010; Hartmann et al. 2015), low-temperature (Liu et al. 2012), pathogens (Li et al. 2015), ABA (Yoshida et al. 2010), salicylic acid (SA) (Zhou et al. 2018) and gibberellin (GA) (Fukazawa et al. 2000). For instance, OsbZIP12/OsABF1, OsbZIP72, GmbZIP44, GmbZIP62 and GmbZIP78 regulate salt stress tolerance in plants (Liao et al. 2008; Lu et al. 2009; Amir-Hossain et al. 2010). Among them, OsbZIP12/OsABF1 and OsbZIP72 are positive regulators of ABA signaling (Lu et al. 2009; Amir-Hossain et al. 2010), while overexpression of GmbZIP44, GmbZIP62 and GmbZIP78 weakens ABA sensitivity in Arabidopsis (Liao et al. 2008). Furthermore, OsbZIP16, OsbZIP52/RISBZ5 and OsbZIP71 are involved in drought stress in rice (Chen et al. 2012; Liu et al. 2012, 2014a). Under cold conditions, the LIP19-OsOBF1 heterodimer induces the expression of cold-responsive genes (Shimizu et al. 2005), whereas OsbZIP52/RISBZ5 is a negative regulator of cold stress (Liu et al. 2012). However, little is known about the biological functions of the bZIP gene family in yellowhorn.

Yellowhorn (Xanthoceras sorbifolium, NCBI: txid99658), the only species of the genus Xanthocera (Sapindaceae), is a woody oil tree (Bi et al. 2019). It is valuable for economy, pharmacology and ecology. The oil content of its seed kernels can be as high as 67% (Bi et al. 2019). Its numerous bioactive compounds include triterpenoid saponins and barringenol-like triterpenoids that have antitumor and anti-inflammatory activities and perhaps action against Alzheimer’s disease (Li et al. 2016; Wang et al. 2017). As an endemically significant species in northern China, yellowhorn is widely utilized for soil and water conservation owing to its capacity of resistance to cold, drought and salinity. It is a preferred species for greening and beautification and desertification control. It is also ideal for studying mechanisms of resistance to stresses. For that purpose, the yellowhorn genome has been completely sequenced and assembled (Bi et al. 2019). Since bZIP transcription factors function in stress responses, we identified and systematically analyzed these genes in the yellowhorn genome.

In this study, we identified 64 bZIP genes in the yellowhorn genome and analyzed phylogenetic relationships, group classification, gene structure, conserved domains, chromosomal distribution, evolution relationships and cis-element prediction of XsbZIP genes. Using RNA-seq and qRT-PCR, we assessed gene expression patterns of XsbZIPs in response to abiotic stresses and ABA to provide a foundation for deeper studies of XsbZIP biological functions in this woody species.

Materials and methods

Identification of the bZIP gene family in yellowhorn

The genome annotations of yellowhorn were downloaded from the yellowhorn genome database (Bi et al. 2019). The bZIP protein sequences of Arabidopsis, maize and soybean were obtained from the PlantTFDB database (https://planttfdb.cbi.pku.edu.cn/; Jin et al. 2017). The sequences of all the annotated bZIP proteins in yellowhorn genome were aligned with AtbZIPs, ZmbZIPs and GmbZIPs using blastp (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Sequences with E-value < e−5 were identified by HMMER_ v3.2.1 based on the hidden Markov model (HMM) profile of the bZIP domain (PF00170) downloaded from the Pfam database (https://pfam.xfam.org/; Mistry et al. 2013; El-Gebali et al. 2019). We selected the candidate XsbZIPs with E-value ≤ 0.01. Next, the online software CD-Search (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and SMART (https://smart.embl-heidelberg.de/) were used to confirm the integrity of the bZIP domains (Marchler-Bauer et al. 2017; Letunic et al. 2018). Finally, the validated genes were assigned as yellowhorn bZIP genes. Isoelectric point (PI) and molecular weight (MV) of the identified bZIP proteins were obtained using the online software ExPASy (https://web.expasy.org/compute_pi/; Gasteiger et al. 2003), with subcellular location predicted using the online software MBC (https://cello.life.nctu.edu.tw/; Yu et al. 2006).

Phylogenetic analysis and classification of the yellowhorn bZIP family

Based on the bZIP domain sequences of AtbZIP, ZmbZIP, GmbZIP and XsbZIP proteins, the phylogenetic tree was constructed using the unweighted pair group method with arithmetic mean (UPGMA) method in MEGA6 (https://www.megasoftware.net/) program with the following optimized parameters: mode, p-distance; rates among sites, uniform rates; pattern among lineages, same (homogeneous); gaps, complete deletion (Tamura et al. 2013). The bootstrap analysis was carried out with 1000 replicates. The unrooted phylogenetic tree was visualized using iTOL (https://itol.embl.de/; Letunic et al. 2019).

Gene structure and conserved motif analyses

Detailed information on the XsbZIP gene structure was derived from the general feature format (GFF) files in the yellowhorn genome database. TBtools (https://github.com/CJ-Chen/TBtools) was exploited to analyze and display the intron/exon structural patterns of XsbZIP genes (Chen et al. 2018). Intron patterns within the basic and hinge regions of the bZIP domains were characterized and divided into different types according to the illuminated sequences of all XsbZIP genes. The online program MEME-5.0.4 (https://meme-suite.org/) was used to identify conserved motifs of XsbZIP proteins including the bZIP domain with the maximum number of motifs specified as 20 and all other parameters set to the default option (Bailey et al. 2009). The conserved motifs were annotated by the Pfam database (https://pfam.xfam.org/; El-Gebali et al. 2019).

Chromosomal distribution and gene duplication analyses

Gene duplication events were analyzed using MCScanX (https://chibba.agtec.uga.edu/duplication/mcscan) with default parameters (Wang et al. 2012). TBtools (https://github.com/CJ-Chen/TBtools) was used to determine the location of XsbZIP genes and segmental duplication regions on yellowhorn linkage groups (Chen et al. 2018). Syntenic maps of yellowhorn associated with six representative species were visualized by MCScanX (Wang et al. 2012). The Ka (nonsynonymous)/Ks (synonymous) ratio for each bZIP gene pair was calculated using KaKs_Calculator 2.0 (Wang et al. 2010).

Prediction of cis-elements in XsbZIP promoters

DNA sequences of the 2 kb upstream initiation codons of the XsbZIPs were submitted to the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/; Lescot et al. 2002). The cis-elements related to stress and hormone treatments were selected and showed in a figure corresponding to the phylogenetic tree of the XsbZIPs.

Plant preparation and treatment

Yellowhorn seedlings were grown in a 1:3 mixture of vermiculite and humus in a culture room at 22–24 °C with 60% relative humidity and 12 h light/12 h dark for 28 days before treatment. For salt stress, four-week-old seedlings were transplanted to the mixture which had been thoroughly watered with 150 mM NaCl solution before transplant. And these seedlings weren’t watered again during two-day treatment. For phytohormone treatment, seedlings were treated by foliar spraying with 100 μM ABA. For cold stress, seedlings were placed in a 4 °C growth chamber, and the other conditions were the same as untreated. For each treatment, three independent plants (biological replications) were sampled with three corresponding untreated control, and the leaves were collected at 3, 6, 9, 12, 24 and 48 h. All leaf samples were frozen in liquid nitrogen immediately and stored at − 80 °C until analyzed.

Quantitative real-time PCR analysis

Total RNA was extracted using OminiPlant RNA Kit (DNaseI) (CWBIO, Beijing, China) and quality and quantity tested using agarose gel electrophoresis and a DHS NanoPro 2010 Spectrophotometer (Dinghaoyuan, Beijing, China). For first-strand DNA synthesis, 500 ng high quality total RNA was reverse transcribed using TaKaRa PrimeScriptTMRT reagent Kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China) and the manual. Quantitative real-time PCR (qRT-PCR) primers of the XsbZIP genes were designed using the predicted cDNA sequences (Table S1). All primers were tested by PCR amplification, electrophoresis gels and melting curve analyses. qRT-PCR was carried out using a Bio-Rad/CFX Connect Real-Time PCR Detection System and the UltraSYBR Mixture (CWBIO, Beijing, China). The cycling conditions were: 95 °C for 10 min, and 40 cycles of 95 °C for 30 s, 60 °C for 30 s, and 72 °C for 30 s. Additionally, relative expression of XsbZIP genes was analyzed by the 2−ΔΔCt method and normalized against the average transcript quantity of XsACTIN (EVM0013888.1) gene. The qRT-PCR experiments were conducted with three replicates.

Analysis of transcriptome data

A total of 6 sequence libraries, including for 24 h and 48 h after treatment with NaCl, cold or ABA, were constructed to analyze the expression pattern of the XsbZIP genes. RNACocktail (https://bioinform.github.io/rnacocktail/) was utilized to analyze the transcriptome data (unpublished) from our laboratory (Sahraeian et al. 2017). Transcript abundance of XsbZIP genes was calculated as transcripts per million (TPM). Heatmaps and hierarchical cluster analyses were constructed by pheatmap (https://cran.r-project.org/web/packages/pheatmap/index.html) in the R environment (R Foundation for Statistical Computing, Vienna, Austria).

Results

Identification and nomenclature of bZIP transcription factors in yellowhorn

To identify XsbZIP genes, we used blastp with an E-value < e−5 and the program HMMER3.0 with E-value threshold of 0.01 to search bZIP genes in the yellowhorn genome (Mistry et al. 2013). A total of 64 bZIP genes were identified in the yellowhorn genome, and 56 nonredundant genes were named XsbZIP1 to XsbZIP56 in accordance with their physical location on the linkage groups. Four redundant genes (EVM0006705.1, EVM0019189.1, EVM0021491.1 and EVM0020123) were named as XsbZIP8.1, XsbZIP47.1, XsbZIP48.1 and XsbZIP56.1, respectively. Five bZIP genes (EVM0019189.1, EVM0023400.1, EVM0001297.1, EVM0011266.1 and EVM0020934.1) that could not be mapped to any linkage groups were named XsbZIP47.1 and XsbZIP57 to XsbZIP60, respectively.

Domain analysis revealed that 57 of the potential XsbZIP proteins have the typical bZIP domain with the basic region containing an invariant N- × 7-R/K motif and the leucine zipper composed of several aheptad repeat of Leu or other bulky hydrophobic amino acids positioned exactly nine amino acids upstream of R/K toward the C terminus. For the remaining seven, in XsbZIP28, the heptad repeat of Leu is positioned five amino acids toward from C terminus instead of the usual nine. In XsbZIP50, there is an isoleucine (Ile) in place of the conserved arginine (Arg)/lysine (Lys) in the basic region, which is similar to those reported in rice previously (Nijhawan et al. 2008). In the sequences of XsbZIP8, XsbZIP8.1, XsbZIP56 and XsbZIP56.1, the conserved asparagine (Asn) has been replaced by Lys; such unusual replacement in the bZIP basic region has been observed in the yeast MET4 and rice bZIP21 and bZIP82 (Thomas et al. 1992; Nijhawan et al. 2008). XsbZIP9 also has a substitution of aspartic acid (Asp) for the conserved Asn.

Gene characteristics, including the full-length coding sequence (CDS) and the amino acid residues, protein molecular weight (MW), isoelectric point (pI) and subcellular localization, were analyzed (Table S2). The smallest of the 64 XsbZIP proteins was XsbZIP16 with 135 amino acid (aa) and 408-bp CDS, whereas the largest one was XsbZIP13 (1317 aa and 3954-bp CDS). The MW of proteins ranged from 15,699.02 (XsbZIP16) to 144,949.3 (XsbZIP13) Da, and the pI varied from 4.48 (XsbZIP29) to 10.37 (XsbZIP57). Remarkably, all of the XsbZIPs were located in the nucleus, consistent with functioning as regulators of target genes in nucleus.

Phylogenetic analysis and classification of the XsbZIP family

To explore the phylogenetic relationships of the XsbZIP transcription factor family, the bZIP domain amino acid sequences of yellowhorn, Arabidopsis, soybean and maize bZIP proteins were applied to establish the unrooted phylogenetic tree (Fig. 1). The phylogenetic tree consists of 64 XsbZIPs, 75 AtbZIPs, 138 GmbZIPs and 125 ZmbZIPs subdivided into 10 groups, supported by classifications in Arabidopsis, soybean and maize (Jakoby et al. 2002; Wei et al. 2012; Wang et al. 2015). Ten groups of XsbZIPs were named as A to I and S, and consisting of 13, 1, 3, 8, 7, 2, 3, 3, 8 and 13 proteins, respectively. However, four small unique clades, AtbZIP60/XsbZIP29/GmbZIP12, AtbZIP72/XsbZIP43, AtbZIP62/XsbZIP50/GmbZIP21/GmbZIP128/ZmbZIP51/ZmbZIP118 and GmbZIP75/GmbZIP102, might have individual evolutionary trajectories that differ from other clades and were called unclassified group (group U). Furthermore, the amino acid residues in bZIP domains of XsbZIP proteins are highly conserved within each group (Fig. S1).

Fig. 1
figure 1

Phylogenetic tree of yellowhorn, Arabidopsis, soybean and maize bZIP proteins based on bZIP domain in protein sequences of XsbZIPs, AtbZIPs, GmbZIPs and ZmbZIPs. The phylogenetic tree was constructed using the UPGMA method in the MEGA 6 program. The bZIPs clustered into 11 groups (A to I, S and U) as indicated by different-colored arcs. Group names are shown after genes within the corresponding groups. Purple triangles: XsbZIPs; ochre yellow stars: AtbZIPs; brown stars: GmbZIPs; light blue stars: ZmbZIPs. Bootstrap values from 1,000 replicates are shown at branches

Full size image

Gene structure of XsbZIP genes

The structure of each gene, including the intron position, phase and length, likely presents a record of key evolutionary events and thus offers insight on the emergence and evolution of a given gene family (Betts et al. 2001). Therefore, intron numbers, positions and splicing phases of XsbZIP genes were examined to understand their evolutionary trajectory (Fig. S2). Ten of the 64 XsbZIPs lack introns, a state exclusive to group S. Among intron-containing XsbZIP genes, the number of introns within the open reading frame (ORF) varies from 1 to 14, and intron position and splicing sites within the ORF of XsbZIPs differ from each other. However, the intron patterns in the bZIP domains are highly conserved. XsbZIP genes presented five types of patterns (a, b, c, d and e) on the basis of the intron number, position and splicing phase within the basic and hinge regions (Fig. 2 and Fig. S3). Pattern a has a single intron in phase 0 (P0) at the − 5 position within the hinge region. Two introns both having P0 splicing phase models, one at the − 5 position insertion between codons for Lys and alanine (Ala) in the hinge region and the other at the − 25 position in the basic region, were observed in pattern b. Both pattern c and d have just one intron within the basic region in phase 2 (P2), an interruption of the codon at the − 22 and − 17 position. Pattern e, devoid of any intron in the basic or hinge regions, comprised 16 members, 10 of which lack introns in the ORF. Note that the five intron patterns in yellowhorn corresponded one-to-one with those observed in castor bean, and one or more of them were also consistent with certain patterns in rice, maize, barley and tomato, owing to the same positions and splicing phases of introns (Nijhawan et al. 2008; Wei et al. 2012; Jin et al. 2014; Li et al. 2015; Pourabed et al. 2015). Thus, intron patterns remain considerably conserved in the basic and hinge regions of the bZIP genes in plants. All members of most groups classified above share the same intron/exon structural patterns: group A and G sharing pattern a; group D sharing pattern b; group C, H and I sharing pattern c; group B, F and S sharing pattern e. Patterns c and d, having similar structure with identical splicing phases and different positions of introns, contain all genes in group E. In summary, the intron patterns, acting as an index for the phylogenetic relationships in a gene family, have been well conserved during evolution and play an important role in the functional divergence of XsbZIP genes.

Fig. 2
figure 2

Intron patterns within the basic and hinge regions of encoding bZIP domains in XsbZIP genes. The five patterns (a–e) are based on intron numbers, positions and splicing phases. Phase 0 (P0) and phase 2 (P2) mean the intron splicing sites are located between codons and between the second nucleotide and the third nucleotide in one codon, respectively. Horizontal black bars represent codons for amino acid sequences within the basic and hinge regions; vertical lines represent intron positions. Details on intron positions within the region encoding the bZIP domains of the XsbZIP genes are given in Supplementary Fig. S3

Full size image

Conserved domain analysis of the XsbZIP proteins

The bZIP domain is the core of the bZIP proteins and can bind to some specific cis-elements in the promoters of their downstream target genes. Nonetheless, additional domains may contribute to the functional diversity of the bZIP proteins. In total, 20 conserved motifs were detected in 64 XsbZIP proteins by the MEME analysis tool (Fig. 3 and Table S3). Motif 1 is the basic region of the bZIP domain; motif 3 and 9 are leucine zipper regions. In addition, the overall distribution of these conserved motifs is similar within the same group but differs from that in other groups, which indicates similar function of members within the same group but divergent from other groups and also implies that the group classification is reliable. Most of the motifs are shared by several groups, such as motif 6, 12, 13 and 17 in two groups; motif 7, 9, 16, 19 and 20 in three groups; and motif 18 in four groups. However, motif 5, 8, 10 and 11 are exclusive to group A, and motif 2, 4, 14 and 15 are only present in group D.

Fig. 3
figure 3

Phyolgenetic tree based on conserved motifs of XsbZIP proteins. Twenty conserved motifs, numbered 1–20 and coded with different colors in the key on the right, were detected among the 64 XsbZIP proteins using the MEME tool. Structures of XsbZIPs based on the position of conserved motifs are shown to the right of the tree. Nonconserved sequences are depicted by grey lines. Details on predicted motifs are given in Supplementary Table S3. Scale at the bottom shows estimated length of proteins

Full size image

Previous studies suggested that several motifs might be linked to specific biological functions. For example, motif 4 and 17 contain potential Ca2+-dependent protein kinase phosphorylation sites (R/KxxS/T), shown as R(QHK)(QE)T and RSHS, respectively (Jakoby et al. 2002; Nijhawan et al. 2008). And motif 5 and 11 also have potential phosphorylation sites (S/TxxD/E) presented as TLE(DE) and S(SN)(AG)D, respectively, for casein kinase II (CKII) (Jakoby et al. 2002; Nijhawan et al. 2008). Motif 16, characterized by a proline-rich domain having transcriptional activation potential, is present in XsbZIP11 and XsbZIP59 (Nijhawan et al. 2008).

Chromosomal distribution and evolution analysis of XsbZIP genes

As shown in Fig. 4a, 59 XsbZIP genes were disproportionately mapped on the yellowhorn linkage groups except LG5. LG4 with the largest number of bZIP genes presents 10 members, whereas a single XsbZIP is distributed on LG6.

Fig. 4
figure 4

Chromosomal distribution and analysis of evolution of XsbZIP genes. a XsbZIPs genomic distribution and collinear relationships. A total of 59 XsbZIPs were disproportionately mapped on the yellowhorn linkage groups. Segmentally duplicated genes are purple; the others are black. b Syntenic relationships of bZIP genes between yellowhorn and six representative species. Gray lines indicate all syntenic blocks among yellowhorn linkage groups or between yellowhorn and the other species; collinear pairs of bZIPs are connected by red lines. Species: Xanthoceras sorbifolium, Arabidopsis thaliana, Glycine max, Malus domestica, Oryza sativa, Populus trichocarpa and Zea mays

Full size image

During the evolution of plants, gene families generally have undergone gene duplication, thus driving their quantitative expansion and functional diversification. To explore the evolutionary mechanisms of XsbZIP genes, we analyzed the contribution of tandem duplication and segmental duplication. According to Holub’s (2001) perspective, tandem duplication was defined as a chromosomal region within 200 kb containing two or more collinear genes. No tandem duplication was found between XsbZIP genes. In addition, there are 14 pairs of segmentally duplicated genes, comprising 19 XsbZIPs in the segmental duplication regions of yellowhorn linkage groups, which may be attributed to genome fusion or whole-genome duplication (Fig. 4a; Table S4). Furthermore, these segmentally duplicated XsbZIP pairs were present either within a linkage group or between linkage groups and might have been duplicated once, twice or three times. The above results indicate that segmental duplication events predominantly contributed to the evolution of the XsbZIP gene family.

To further elucidate the evolution mechanisms of yellowhorn bZIP family, we constructed six comparative syntenic maps of yellowhorn associated with representative species containing four dicots (Arabidopsis, soybean, apple and Populus trichocarpa) and two monocots (rice and maize) (Fig. 4b, Table S5). A total of 51 (Arabidopsis), 108 (soybean), 74 (apple), 89 (Populus trichocarpa), 36 (rice) and 30 (maize) orthologous pairs were found between yellowhorn and the other species. Interestingly, among these gene pairs, eight XsbZIP genes (XsbZIP3, XsbZIP11, XsbZIP30, XsbZIP32, XsbZIP36, XsbZIP37, XsbZIP38 and XsbZIP54) were shown to have collinear relationships with all of the above six species, showing that these orthologous pairs might already have been present before the divergence of dicotyledonous and monocotyledonous plants. Some collinear gene pairs with 10 XsbZIP genes (XsbZIP1, XsbZIP2, XsbZIP7, XsbZIP8, XsbZIP12, XsbZIP23, XsbZIP26, XsbZIP41, XsbZIP42 and XsbZIP53), however, were exclusively identified between yellowhorn and dicots (Arabidopsis/soybean/apple/populus trichocarpa), indicating that these pairs might have formed after ancestral divergence. Besides, the Ka/Ks values of the bZIP pairs were calculated to better understand the selective pressure on bZIP family genes. For most orthologous bZIP gene pairs, the Ka/Ks values were less than 1, which suggested that the XsbZIPs had mainly undergone purifying selection.

Prediction of cis-elements in XsbZIP promoters

To explore the potential regulatory mechanisms of XsbZIPs in response to various stimuli, cis-elements associated with stress and hormone treatments in the 2-kb upstream initiation codons of the XsbZIP genes were analyzed using the PlantCARE database (Lescot et al. 2002) (Fig. 5; Table S6). Except for XsbZIP38, XsbZIP57 and XsbZIP59, the promoters of all XsbZIPs cover at least two types of cis-elements, such as methyl-jasmonate (MeJA) responsive element CGTCA-motif and TGACG-motif, abscisic acid response element ABRE, drought stress response element MBS, anaerobic induction element ARE, etc. Thus, the transcription of XsbZIP genes might be regulated by diverse cis-elements during stress and hormone signalings. Owing to the cis-elements in gene promoters, 51, 38 and 32 genes among 64 XsbZIP genes are probably involved in the response to ABA, MeJA and SA, respectively. In addition, 48, 37 and 30 genes may be related to anaerobic, drought and low temperature responses, respectively, and 58 have at least one cis-element related to stress responses in their promoter regions. Among individual genes, the promoter of XsbZIP16 contains seven AREs, implying that it possibly functions in anaerobic induction. XsbZIP27, having five MBS in promoter region, may take part in drought stress. Therefore, XsbZIP gene family probably plays a vital role in the regulation of stresses in yellowhorn.

Fig. 5
figure 5

The cis-elements in the promoter regions of XsbZIP genes. The cis-elements, in the 2000 bp upstream of initiation codons and related to stress and hormone treatments, were found in 61 XsbZIP genes using the PlantCARE tool. And color-coded based on the phylogenetic tree of the XsbZIPs. Grey lines indicate nucleotide sequences without the cis-elements. Positions of cis-elements are shown by the scale at the bottom

Full size image

Expression of XsbZIP genes in response to abiotic stresses and ABA

Increasing evidence has shown that bZIP genes play a crucial role in response to abiotic stress and phytohormone treatments in higher plants. To gain more insight into the possible functions of bZIP genes in yellowhorn, we used RNA-seq data to document expression levels of XsbZIP genes at 24 h and 48 h after treatment with abiotic stresses (NaCl and Cold) and ABA (Table S7). A heatmap was constructed using the TPM values of each XsbZIP genes (Fig. 6a). Most XsbZIP genes had different expression patterns in response to at least one of the treatments. Next, we selected 18 genes that were up- or downregulated (fold change > 2, P < 0.05) under one or more of the treatments for further study. Seven were validated by qRT-PCR, and these results were highly consistent with the RNA-seq (Fig. 6b). As for the eighteen genes with prominent differential expression, four were upregulated, and three were downregulated in response to NaCl; one was upregulated, and one was downregulated in response to ABA; ten were upregulated, and six were downregulated in response to cold. Interestingly, five (XsbZIP10, XsbZIP21, XsbZIP34, XsbZIP42 and XsbZIP54) responded to both NaCl and cold, while two (XsbZIP3 and XsbZIP17) responded to both ABA and cold, indicating that these genes are probably involved in multiple environmental adaptations in yellowhorn.

Fig. 6
figure 6

Expression of XsbZIP genes over time after salt, cold or ABA treatment. a Heatmap based on TPM values of the relative expression levels of 64 XsbZIP genes in yellowhorn leaves after treatments. b Expression of representative XsbZIP genes validated by qRT-PCR. Data was analyzed by the 2−ΔΔCt method and normalized to XsACTIN (EVM0013888.1). Error bars indicate standard deviations. Relative expression at 0 h was set as 1

Full size image

Discussion

Evolution of XsbZIP genes

We identified 60 non-redundant bZIP genes in yellowhorn (Table S2), fewer than in other plants, except castor bean (49, Jin et al. 2014) and grapevine (55, Liu et al. 2014b) (Arabidopsis thaliana, 75 members (Jakoby et al. 2002); rice, 89 (Nijhawan et al. 2008); maize, 125 (Wei et al. 2012); tomato, 69 (Li et al. 2015); sorghum, 92 (Wang et al. 2011); sesame, 63 (Wang et al. 2018b); barley, 89 (Pourabed et al. 2015); soybean, 131 (Zhang et al. 2018); apple, 112 (Zhao et al. 2016); Chinese cabbage, 136 (Bai et al. 2016); rapeseed, 247 (Zhou et al. 2017); cucumber, 64, Baloglu et al. 2014). It also seems likely that the monocots share more bZIP genes than the dicots, possibly owing to the higher number of bZIP genes that evolved in monocots than in dicots after the ancestral divergence. According to the prediction, all XsbZIPs were located in nucleus (Table S2). Most of the XsbZIPs share highly conserved bZIP domain, with the basic region having an invariant N- × 7-R/K motif and the leucine zipper having a heptad repeat of Leu or other bulky hydrophobic amino acids, while some XsbZIPs (XsbZIP8, XsbZIP8.1, XsbZIP9, XsbZIP28, XsbZIP50, XsbZIP56 and XsbZIP56.1) have substitutions (Fig. S1).

Among the 11 groups identified in the phylogenetic analysis of the 64 yellowhorn bZIP proteins (Fig. 1), five intron patterns were found within the bZIP domains and additional conserved motifs identified, both being highly group-specific and possibly linked to functional diversity of the XsbZIP genes (Figs. 2, 3). XsbZIP genes were not evenly distributed on the yellowhorn linkage groups except for LG5, and the number of these genes ranged from 1 to 10 on every LG (Fig. 4a). This disproportionate distribution of bZIP genes is also present in other species, such as barley (Pourabed et al. 2015), grapevine (Liu et al. 2014b) and cucumber (Baloglu et al. 2014). In addition, expansion of XsbZIP genes might be the consequence of segmental duplication, consistent with the results in Arabidopsis thaliana (Jakoby et al. 2002), rice (Nijhawan et al. 2008) and grapevine (Liu et al. 2014b). The synteny analysis and phylogenetic comparison of bZIP genes also provided valuable clues about the evolutionary characteristics of yellowhorn bZIP genes.

Expression and potential functions of XsbZIPs in response to abiotic stresses and ABA

In group A of bZIP transcription factors are the ABA-responsive element binding proteins (AREBs)/ABA-responsive element binding factors (ABFs) which have been suggested to function in ABA and stress signalings. For instance, AtbZIP36/ABF2/AREB1, AtbZIP37/ABF3 and AtbZIP38/ABF4/AREB2 have overlapping roles in response to water stress and cooperatively regulate the expression of the ABRE-dependent genes in ABA signaling under water stress (Yoshida et al. 2010). Overexpression of SIAREB induced the expression of AtRd29A and AtCOR47 in Arabidopsis and induced the expression of SICI7-like dehydrin in tomato, which increased their tolerance to drought and salt (Hsieh et al. 2010). Also, OsbZIP72, an ABRE binding transcription factor, plays a positive part in drought resistance through ABA signaling (Lu et al. 2009). In addition, OsbZIP71 participates in ABA-mediated drought and salt tolerance by directly binding to the promoters of the abiotic stress-related genes, OsNHX1 and COR413-TM1 (Liu et al. 2014a). In japonica rice, bZIP73 interacts with bZIP71 and remarkably enhances resistance to cold by regulating ABA and ROS homeostasis (Liu et al. 2018). ZmbZIP72 is induced by ABA, high salinity and drought treatments in maize and its heterologous overexpression improved resistance under drought and salt stresses in Arabidopsis (Ying et al. 2012). The group S1 bZIP TFs (AtbZIP1 and AtbZIP53), functioning in reprogramming carbohydrate and amino acid metabolism, play important roles in the root-specific response to salt (Hartmann et al. 2015). VIP genes also belong to the bZIP gene family, and pathogenesis-related protein 1 (PR1) gene is upregulated and resistance to Brenneria salicis infection increases in transgenic Populus overexpressing PtVIP1 or AtVIP1 (Wang et al. 2018a).

The cis-elements and expression pattern analyses suggested that XsbZIPs are widely involved in abiotic tolerance. Determining orthologous relationships between the AtbZIP genes, which are the best studied, and XsbZIP genes will allow researchers to predict XsbZIP functions. XsbZIP11 and XsbZIP25, orthologues of AREB/ABFs in Arabidopsis, were highly induced by cold and salt stresses, respectively (Fig. 6). These two XsbZIPs, members of group A, have the conserved motif 5 and 11, which contain potential casein kinase II (CKII) phosphorylation sites and ABRE in their promoters. AtVIP1, orthologues of XsbZIP21, acts as a regulator of osmosensory signaling in Arabidopsis (Tsugama et al. 2014). XsbZIP21 was highly induced by salt and cold treatments in this study, indicating its possible roles in osmosensory response (Fig. 6). Relative expression of XsbZIP35 increased with NaCl treatment (Fig. 6), consistent with its ortholog AtbZIP24, which is induced by salt stress and acts as a negative transcriptional regulator of salt tolerance by repressing the maintenance of osmotic and ionic balance (Yang et al. 2009). However, XsbZIP35 was downregulated by cold stress (Fig. 6), suggesting that it maybe play different roles under various abiotic stresses. These results lay a solid foundation for further study on the evolutionary history and biological functions of bZIP gene family in yellowhorn.

Conclusions

Among the 64 XsbZIP genes identified here for yellowhorn, 59 were disproportionately localized on all LGs except LG5, and they had been named according to their order on the LGs. XsbZIP proteins clustered into 11 groups based on their phylogenetic relationships, which was supported by intron patterns in the bZIP domains and additional conserved motifs. The evolution of XsbZIPs was largely attributed to segmental duplication events. The expression patterns of the bZIP genes in response to abiotic stresses revealed that XsbZIP transcription factors play significant roles in stress responses. Analyses of the conserved motifs and cis-elements and overview of the XsbZIP gene family and their potential functions in stress responses provide a theoretical basis for research on the biological roles of bZIP genes in yellowhorn and breeding programs for stress tolerance.