Genome-wide identification and expression analysis of the bZIP gene family in silver birch (Betula pendula Roth.)

YaNan Jing·Yue Yu·Hanzeng Wang·Ye Wu·Chenghao Li

Abstract Basic Leu zipper (bZIP) genes play pivotal，versatile roles in abiotic or biotic stress responses and in other biological processes.Knowledge on the evolutionary relationships and patterns of gene expression of bZIP family members in woody plants，however，has been limited.Here we identified and characterized 47 BpbZIP genes across the silver birch (Betula pendula Roth.) genome.With reference to bZIP classifications for Arabidopsis thaliana，all BpbZIP proteins clustered among 10 groups in phylogeny.The bZIP domains were divided into five patterns based on intron positions and splicing phases.A total of 24 conserved motifs were detected in BpbZIPs with high group specificity.We also analyzed the protein structure of the BpAREB/ABF/ABI5 subfamily，the most important subfamily in the bZIP family.Expression analyses demonstrated that BpbZIP genes were widely involved in abscisic acid，salt，drought，and heat stress responses，with BpbZIP07/ABF4 and BpbZIP21/ABF2 most highly expressed.Our results on genome-wide identification，evolutionary relationships，gene structure，and motif and promoter element identification for BpbZIP family members in silver birch provide a comprehensive understanding of bZIP transcription factors in birch and will lead to a deeper understanding of their evolution and potential biological functions.

Keywords Basic Leu zipper gene·Birch·Abscisic acid·Abiotic stress·Gene expression

Introduction

The basic Leu zipper (bZIP) transcription factors with the bZIP domain form one of the earliest transcription factor families and provide insight into the molecular structure of their conserved domain，which is characterized by a strictly conservative basic DNA binding region(N-× 7-R/K-× 9) containing 18 amino acids (aa) following an adjacent Leu zip，that comprises several heptad repeats of Leu or other hydrophobic amino acid (Landschulz et al.1988;Jakoby et al.2002).The basic region regulates nuclear localization and DNA binding，and the Leu zipper mainly takes part in dimerization (Kouzarides and Ziff1989;Vinson et al.1989).Foster et al.(1994)proposed thatbZIPgenes preferentially bind to motifs of DNA sequence that contain an ACGT core，including G-boxes (CAC GTG)，A-boxes (TAC GTA)，and C-boxes(GAC GTC).Besides the bZIP domain，additional conserved domains/sites act as transcriptional activators.For instance，ABA-dependent 42-kDa protein kinases are able to phosphorylate conserved domains of AtAREB1，and the domains contain two target sites，R/KxxS/T for calcium-dependent protein kinase (CDPK) and S/TxxD/E for casein kinase II (CKII) (Furihata et al.2006;Zhu et al.2007).

PlantbZIPgenes are involved in diverse biological processes in numerous species;75bZIPgenes have been found inArabidopsis thaliana(Jakoby et al.2002)，125 inZea mays (Wei et al.2012)，49 inRicinus communis(Jin et al.2014)，178 inPanicum virgatum(Wang et al.2020)，and 89 inOryza sativa(Nijhawan et al.2008).Considerable evidence shows thatbZIPgenes take part in the differentiation of various organs and tissues (Izawa et al.1994;Toh et al.2012)，transition to and initiation of flowering (Strathmann et al.2001;Wigge et al.2005)，and embryogenesis (Shiota et al.2008;Guan et al.2009).Growing evidence reveals that bZIP proteins respond to numerous biotic and abiotic stress conditions，including infection (Kaminaka et al.2006)，heat stress (Wang et al.2017，2020)，drought stress (Chen et al.2012;Tang et al.2012)，salinity stress (Uno et al.2000;Hsieh et al.2010)，light irradiation (Ulm et al.2004)，and osmotic stress(Balázs et al.2010).

Silver birch (Betula pendulaRoth.)，an important worldwide pioneer tree species that survived from the ice age，is highly adaptable and optimized for fiber and biomass production.It is also multi-purpose species for urban greening and afforestation to prevent soil erosion.The genome sequence ofB.pendulawas published in 2017 (Saloj?rvi et al.2017)，and that ofBetula platyphyllahas also been sequenced recently (Chen et al.2021)，thus facilitating genome-wide and functional analyses.Despite these advances，very few studies have focused on gene function in birch，except for the role ofBpCUC2in regulating leaf shape and internode development (Liu et al.2019a) and identification ofBpNACgene family (Chen et al.2019).

Elucidation of the expression patterns of members of the plantbZIPgene family under abiotic stress will be instrumental to deepening our understanding of genetic and molecular mechanisms underlying stress responses and thus improve stress tolerance in plants.In the present study of theB.pendulagenome，47BpbZIPgenes were found and analyzed using bioinformatics methods determine their chromosomal location，phylogenetic relationships，conserved motifs，exon-intron organization，andcis-acting promoter elements.Quantitative real-time reverse transcription-polymerase chain reaction (qRTPCR) assays were also used to analyze gene expression ofBpbZIPgenes when plants were exposed to stress treatments using abscisic acid (ABA)，salt，drought，and heat stress.These fundamental findings on the function ofBpbZIPgenes during stress treatment provide a basis for improving wood tolerance to adverse environments.

Materials and methods

Identification of bZIP genes in silver birch

The silver birch genomic sequence was obtained from the Comparative Genomics (CoGe) database ( https://genom evolu tion.org/CoGe/Genom eInfo.pl? gid=35079) (Saloj?rvi et al.2017).The hidden Markov model (HMM) (ID:PF01700，PF07716，PF03131，and PF12498) and the web server HMMER 3.3 (Finn et al.2011) were used with a cutoffE-value of 0.01 for a sequence similarity search for the bZIP domain in the genomic sequence.Candidate BpbZIPs were further analyzed using the SMART database ( http://smart.embl-heidelberg.de) (Schultz et al.1998).The subcellular locations of BpbZIP proteins were then predicted using the WoLF PSORT database ( https://wolfpsort.hgc.jp)(Hortonetal.2007).

Phylogenetic and conserved motif analyses

All full-length bZIP amino acid sequences for the model plantPopulus trichocarpawere obtained from Phytozome(www.phyto zome.net) and forA.thalianafrom The Arabidopsis Information Resource (TAIR，https://www.arabidopsis.org).

Clustal X and BioEdit were used for multiple alignment and manual editing with a gap opening penalty and gap extension penalty of 10 and 0.1，respectively (Thompson et al.1997;Tippmann 2004).Then，an unrooted evolutionary tree was generated by MEGA X software using the neighbor-joining (NJ) algorithm based on 1000 bootstrap replicates and the Poisson model to support each node in the tree (Kumar et al.2018).The web-based tool MEME 5.1.1(https://meme-suite.org/tools/meme) was used to search for additional conserved motifs using default options except for optimum motif width of 6-200 aa and number of motifs of 50.

Gene structure and chromosomal locations analyses

The CoGe database was used to obtainB.pendulagenomic structural information (Saloj?rvi et al.2017).We visualized theB.pendulagenomic gene density profile and locations ofBpbZIPgenes on chromosomes using TBtools v1.082 (Chen et al.2020).The Gene Structure Display Server (GSDS，http://gsds.cbi.pku.edu.cn/) was used for graphical representation of the exon-intron organization ofBpbZIPgenes (Hu et al.2015).

Promoter analysis

As regulatory promoter regions ofB.pendula，2000-bp upstream regions of ATG were fetched by TBtools v1.082.The PlantCARE database was used for analysis ofcis-elements ( http://bioin forma tics.psb.ugent.be/webto ols/plant care/html/) (Lescot et al.2002).

Identification of the BpAREB/ABF/ABI5 subfamily

Using nineA.thalianaAREB/ABF/ABI5 protein sequences(Suzuki et al.2003) as references，the wholeB.pendulagenome was searched for BpAREB/ABF/ABI5 subfamily members using BioEdit.Jalview 2.11 (http://www.jalvi ew.org) was used for protein sequence alignment and analysis.

Plant material and stress treatments

One mature seed embryo of wild typeB.pendulawas used as the source of the clonally propagated plants and the collection and differentiation of seed were as previously described (Liu et al.2019a).Clonally propagated birch in vitro plantlets were cultured on solid rooting medium [1/2 Murashige and Skoog (MS) plus 0.02 mg/L of naphthalene acetic acid (NAA)] with 16 h light/8 h dark photoperiod at 25°C.Under the culture condition，one-month-old birch in vitro plantlets in the tissue culture flasks were then transferred to the solid rooting medium containing the following chemicals respectively for treatments:200 μM ABA，7% v/v polyethylene glycol 6000 (PEG6000)，150 mM sodium chloride (NaCl)，or a heat stress in a lit growth chamber at 42°C for 0 h，1.5 h，and 3 h.After 0 h，3 h，6 h，12 h，24 h，and 7 d，apical buds were collected from three independent birch plants for each treatment for qRT-PCR assays.

RNA extraction，primer design，and qRT-PCR

Total RNA was isolated from each apical bud separately using a Universal Plant Total RNA Extraction Kit，and gDNA Wiper Mix was used to remove the genomic DNA(BioTeke，Beijing，China).Agarose gel electrophoresis was used for quality control of the isolated RNA.HiScript II Q RT Supermix Kit (Vazyme，Nanjing，China) was used to synthesize first-strand cDNA.Primer pairs ofBpbZIPgenes for qRT-PCR (Supplementary Table S1) were designed using the IDT Realtime PCR tool ( https://sg.idtdna.com/Scito ols/Appli catio ns/RealT imePCR)，and NCBI Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/)was used to assess specificity.Before qRT-PCR，all primer pairs were tested for amplification of targets using PCR and gel electrophoresis.

ChamQ Universal SYBR qPCR Master Mix (Vazyme)was used for the qRT-PCR.BpUbiquitin(BpUBQ) andBpTubulinwere chosen as internal birch reference genes for normalization of gene expression.Expression levels were calculated using the 2-ΔΔCtmethod (Livak and Schmittgen 2001).

Statistical analyses

The mean-expression levels after a stress treatment at each time point were compared to the control group (0 h) using an analysis of variance (ANOVA) in SPSS v26.0 (IBM，Armonk，NY，USA) withp＜0.05 (*) andp＜0.01 (**).

Results

Predicted BpbZIP proteins and genes in genome

In silver birch，49 putative BpbZIP proteins were identified by HMMER 3.3.Analysis of conserved domains using the SMART database showed that all but two of the putative BpbZIPs had the bZIP domain (bZIP_1，bZIP_2，or bZIP_C).Among all predicted BpbZIP proteins，45 had either the bZIP_1 or bZIP_2 domain，demonstrating the reliability of the results.Two putative BpbZIPs，BpbZIP13 and BpbZIP14，had the bZIP_C domain，as found in AtbZIP25 and AtbZIP63 (Jakoby et al.2002).In the remaining two BpbZIP proteins (Supplementary Table S2)，no conserved bZIP domain was found，which may indicate a truncated BpbZIP;both were excluded from subsequent studies.

The 47BpbZIPgenes were renamedBpbZIP01toBpbZIP44based on the coordinates on theB.pendulachromosomes (Table 1;Fig.1).Chromosome mapping showed that 44BpbZIPgenes were located on 12 chromosomes，and their distributions were uneven.Chromosome 6 and 12 had nobZIPgenes，and chromosome 10 had just one.HigherbZIPgene density appeared at the upper end of chromosome 9 and the lower end of chromosomes 7 and 8.Three genes (BpbZIP45-BpbZIP47) were mapped to unanchored scaffolds.

Fig.1 Gene density profile and chromosome locations of BpbZIP genes from Betula pendula.Gene density and BpbZIP genes varied in different scaffold regions across the B.pendula genome.The scale bar is 5 Mb

The 47 BpbZIPs ranged from 145 (BpbZIP15) to 770(BpbZIP23) aa long，with a mean of 341 aa，similar to that forA.thaliana(mean of 321 aa for 75 members) (Jakoby et al.2002)，strawberry (Fragaria ananassa，402.54 aa for 50 members) (Wang et al.2017)，grapevine (Vitis vinifera，323.5 aa for 55 members) (Liu et al.2014)，rice (Oryza saliva，311 aa for 89 members) (Nijhawan et al.2008)，and sorghum (Sorghum bicolor，306.1 aa for 112 members)(Wang et al.2011).Of the 47 BpbZIP proteins，44 were predicted to be in the nucleus，BpbZIP23 in the endoplasmic reticulum，and BpbZIP42 and BpbZIP45 in peroxisomes.Details on the genes are in Table 1.

Table 1 Basic information of the bZIP genes in Betula pendula

Phylogenetic analysis and classification

To analyze the phylogenetic relationships among all bZIPs fromB.pendulaand the two model plant species，a NJ phylogenetic tree was generated with 218 protein sequences，including 47 protein sequences fromB.pendula，97 fromP.trichocarpa，and 75 fromA.thaliana.As shown in Fig.2，based on the branches of phylogenetic tree and classification method ofA.thalianabZIP proteins，the 47 BpbZIP proteins were classified into 10 distinct groups，except for BpbZIP27 in a separate group，named group U.

Fig.2 Phylogenetic analysis of bZIP proteins from Betula pendula，Arabidopsis thaliana，and Populus trichocarpa.Colored rings indicate the different groups (A-I，S，and U) of bZIP proteins.Green circles represent BpbZIPs;red stars indicate BpAREB/ABF/ABI5 subfamily members

Each group of BpbZIP proteins contained the same type of bZIP domain as that ofA.thaliana(Jakoby et al.2002).For example，group C of BpbZIPs and AtbZIPs had the bZIP_C domain.Group F of BpbZIPs and AtbZIPs had the bZIP_2 domain.Group S of BpbZIPs and AtbZIPs had the bZIP_1 domain.BpbZIP proteins were distributed across all groups，although the distribution was uneven.Group A was the largest group with nine BpbZIP proteins，followed by groups S with eight and D with seven.However，groups B，E and U contained only one BpbZIP protein.

Sequence-structure features of the BpbZIP gene family

In some gene families，exon-intron organization and intron patterns are considered to be typical signatures of evolution (Li et al.2015;Wang et al.2020).Therefore，gene structure ofBpbZIPgenes was analyzed to get insights to the emergence and evolutionary history.As for the same phylogenetic groups，the structure ofBpbZIPgenes had several conserved characters (Fig.3 c).NineBpbZIPgenes (19.1%) were intronless，exclusively in group S andBpbZIP37of group F.Similarly，about 20% ofPrunus persica bZIPgenes are intronless (Sun et al.2016).For introns in otherBpbZIPgenes，the quantity varied from one (BpbZIP09，BpbZIP23，andBpbZIP34/AREB3) to 11(BpbZIP16，BpbZIP39，andBpbZIP44) in open reading frames (ORFs)，indicating diverse structures ofBpbZIPgenes.In group D，the average number of introns was the highest (9.3).

Fig.3 Phylogenetic analysis and conserved motifs of BpbZIP proteins and sequence structure of the BpbZIP genes from Betula pendula.Protein structures (b) with conserved motifs and exon-intron organization (c) of BpbZIP genes are rearranged according to the phylogenetic tree (a) of BpbZIP proteins.Boxes in different colors，numbered from 1 to 24，represent different conserved motifs.Exons，introns，bZIP domains and untranslated regions (UTRs) are represented by different color boxes and lines.Numbers above lines represent the splicing phase

Further，within the bZIP domain，the patterns of intron positions and splicing phases related to codons were also analyzed (Supplementary Table S3).The boundaries of the bZIP domain were demarcated based on those ofA.thalianabZIP proteins (Jakoby et al.2002).By analyzing the phase，position and number of introns of the bZIP domain，five patterns，nameda，b，c，d，ande，were found(Fig.4).The exception is patterne，containing all of group C，because the bZIP_C domain did not contain the N-× 7-R/K-× 9 basic region.

Fig.4 Diagram of intron patterns within the bZIP domain of BpbZIP proteins from Betula pendula.At the top，a typical structure of the bZIP domain(bZIP_1 or bZIP_2) is shown，and an example reference sequence is given below.The amino acids in red and green are highly conserved residues.Black vertical lines labeled with P0 or P2 above indicate intron positions and phase.Details of intron positions with BpbZIP domain are given in Supplementary Table S4

Patternawas the most widespread，occurring in 19 of theBpbZIPgenes.Patternboccurred in the fewest members，only three.Across different species，intron patterns corresponded to each other because they shared the same position.Patternsa，b，andcofBpbZIPgenes corresponded to patternsa，d，andcofF.ananassaand patternsa，b，andcofZ.mays，respectively.There were 13BpbZIPgenes in patternd，which had no introns in the bZIP domain，similar to patternfofZ.mays(Wei et al.2012;Wang et al.2017).

Most of theBpbZIPgenes in each group had the same intron pattern.Patternahad only one phase 0 (P0) intron behind the position of-6 (Q or K) in the hinge and occurred in group A，group G withoutBpbZIP43，and half of the members of group D.The other members of group D formed patternb，having two P0 introns;one was inserted after position-26 (K) in the basic region，and the other had the same phase and position with patterna.In addition，patterncwas divided into three subcategories:c-I，c-II，andc-III.The three subcategories shared an identical insertion of phase 2(P2) in front of position-22 (R) in the basic region but had distinct insertion cases in the Leu zippers.Patterncoccurred in groups I，E，and H except inBpbZIP02.Patterndlacked introns in the bZIP domain，and this pattern occurred in all members of groups S，B，and F;BpbZIP43in group G;andBpbZIP02in group H.None of the patterns contained an intron of phase 1 (P1).

Identification of conserved motifs

In BpbZIP proteins，24 distinct conserved motifs were detected (Fig.3 b)，and the logos of multilevel consensus sequences are illustrated in Supplementary Table S4.Overall，all BpbZIPs contained motif 1，which represented the basic region and hinge.Motifs 4，5，10，and 20 represented the diversity of the Leu zipper structure.All BpbZIPs had motif 1 and one of motifs 4，5，10，except BpbZIP02，BpbZIP23，BpbZIP28，and BpbZIP32.Notably，the Leu zippers were conserved without exception in the different groups.

Apart from the bZIP domain，additional motifs that were distributed across different groups were conserved to a certain extent.For instance，motifs 8，9，11，and 16 were only present in group A;motifs 18，19，and 21 in group I;motif 15 in group S;and motif 14 in group F.These results indicated that BpbZIP proteins from one group had similar motif characteristics，suggesting functional similarities.The conserved motif results were supported by the phylogenetic relationships and the classification of theBpbZIPgenes.

There is evidence that mostbZIPgenes in group A inA.thalianatake part in ABA signaling transduction pathway (Finkelstein and Lynch 2000;Yoshida et al.2010).Typically，in group A ofAtbZIPgenes，AREB/ABF/ABI5genes have two sites，S/TxxD/E and R/KxxS/T，that can be activated by CKII and CDPK，respectively.(Uno et al.2000;Furihata et al.2006).Therefore，some additional conserved motifs in BpbZIPs may modulate protein-protein interactions.Motifs 8，9，and 16 all contained putative targets for CKII and CDPK.For the CKII putative targets，motifs 8，9，and 16 had multilevel consensus sequences TVDE，TLE[DE]，and T[FL]DE，respectively;for CDPK，these three motifs were RQ[GA]S，[RG]Q[PS][TE]，and RQSS，respectively.Comparisons with additional conserved motifs fromA.thalianamay show thatbZIPgenes inB.pendulahave similar functions.

Promoter cis-elements inBpbZIP genes

In the 2000-bp promoter regions ofBpbZIPgenes，a variety ofcis-elements were identified，containing a few TATAboxes and CAAT-boxes，which are necessary to initiate transcription (Supplementary Table S5).There were four broad categories of response elements:light，hormone，stress，and growth regulation.Twenty forms of light-responsive elements were found (e.g.，AE-box，Lamp-element，GTGGCmotif，MRE，and Sp1).Numerous hormone-responsive elements were also detected such as the CGTCA/TGACGmotif belonging to methyl jasmonate (MeJA)-responsiveness elements;P-box，TATC-box，and GARE-motif classified as gibberellin-responsive elements;ABRE and TCA-elements related to ABA-responsiveness and salicylic acid (SA)-responsiveness，respectively;AuxRR-core and TGA-element related to auxin-responsiveness.For growth regulation，the CAT-box is related to meristem development.Stress-responsive elements，such as drought-related elements (MBS) and anaerobic responsive elements (ARE)were observed.A small number of regulatory elements of other stress responses were also detected，such as LTR for freezing responsiveness.

Our study focused on seven of the most widely distributed elements associated with phytohormone and environmental-stress signal responsiveness (Fig.5，Supplementary Table S6).Most of theBpbZIPpromoters contained multiple identifiedcis-elements.Among these stress-related elements，AREB，MBS，and CGTCA-motif were the most abundant.BpbZIP27had the mostcis-elements (11).BpbZIP03andBpbZIP29contained six of these sevenciselements;BpbZIP11，BpbZIP30，andBpbZIP43had none.The variety ofcis-elements further supported the diversity of functions ofBpbZIPgenes.

Fig.5 The cis-elements of BpbZIP gene promoters that are related to responses to abiotic stress and phytohormones

Protein structure in the BpAREB/ABF/ABI5 subfamily

Since the group A sequence analysis showed obvious correlations between the sequences and the ABA-responsive signaling pathway，we further identified AREB/ABF/ABI5 members ofB.pendulausing nineA.thalianaAREB/ABF/ABI5 members (Suzuki et al.2003) as the reference sequences to search theB.pendulagenome.Five proteins were thus identified as members of the BpAREB/ABF/ABI5 subfamily and renamed according to the genome annotation ofB.pendulaasBpbZIP07/ABF4，BpbZIP21/ABF2，BpbZIP26/ABI5，BpbZIP34/AREB3，andBpbZIP38/AREB1.They are homologous to the genes in group A in the phylogenetic tree (Fig.3 a).

As shown in Fig.6，all five BpAREB/ABF/ABI5 proteins had four regions containing potential phosphorylation domains，which play a key role in regulatingAREB1inA.thaliana(Furihata et al.2006).The bZIP basic region of BpAREB/ABF/ABI5 members was highly conserved.In particular，the sequence included four additional conserved domains (C1，C2，C3，and C4)，which were first found inA.thalianaand later inSolanum tuberosumandP.trichocarpa(Furihata et al.2006;Ji et al.2013;Liu et al.2019b).The analysis of the protein sequences showed C1，C2，and C3 domains，corresponding to motifs 16，8，and 9，respectively，near the N-terminal，and the C4 domain at the C-terminal.Multiple S/TxxD/E (C1-2，C2-3，and C3-2) and R/KxxS/T(C1-1，C2-1，C2-2，C3-1，and C4-1) sites were distributed in conserved domains and verified the previous motif analysis.

Fig.6 Multiple sequence alignment of BpAREB/ABF/ABI5 subfamily members.Blue shading marks identical residues;light blue shading marks conserved residues.Positions of the basic region of bZIP domain and conserved domains C1-C4 are demarcated by lines above sequences.Sites of S/TxxD/E(C1-2，C2-3，and C3-2) and R/KxxS/T (C1-1，C2-1，C2-2，C3-1，and C4-1) are framed in red.Three black triangles indicate conserved Leu residues in the Leu zipper

Expression analysis of BpbZIP genes under ABA treatment

In this first study on the ABA-stress response inB.pendula，the expression analysis of the 47BpbZIPgenes after exposure to 200 μM ABA within 7 days showed that expression of 45 of the genes was significantly upregulated or down regulated to varying degrees at various times (p＜0.05).

None of the two members of group C (BpbZIP29andBpbZIP33) responded to ABA treatment.Six (BpbZIP14，BpbZIP18，BpbZIP22，BpbZIP25，BpbZIP35，andBpbZIP37) suppressed significantly to various extents after different durations (Fig.7).For the 24-h ABA treatment，15BpbZIPgenes，including all genes from group A，were significantly upregulated (p＜0.01) compared to the control group (0 h);while five genes (BpbZIP13，BpbZIP18，BpbZIP22，BpbZIP35，andBpbZIP37) were downregulated(p＜0.01).Among these five genes，BpbZIP11，BpbZIP13，andBpbZIP18were homologous genes in group S，but they had contrary trends.The same trends were observed in group D betweenBpbZIP22andBpbZIP24.

Fig.7 (continued)

Notably，throughout the 7-day ABA treatment，at all sampling times，BpbZIP07/ABF4，BpbZIP10，BpbZIP20，BpbZIP21/ABF2，andBpbZIP34/AREB3were significantly upregulated (p＜0.01);andBpbZIP18andBpbZIP22were significantly downregulated (p＜0.01).Expression ofBpbZIP21/ABF2was the highest with 28-fold greater expression compared to the control group (0 h).Further，all members ofBpAREB/ABF/ABI5subfamily were strongly induced (p＜0.01)，which showed that they were closely related to ABA signaling.

Expression analysis of BpbZIP genes under salt，drought，and heat stresses

ABA signaling pathways are known to be closely involved in abiotic stress responses (Suzuki et al.2 016;Zong et al.2016).Subsequently，expression patterns of the 47BpbZIPgenes were again analyzed during salt (150 mM NaCl)，drought (7% v/v PEG6000)，and heat (42°C) stress.

Under salt stress，the expression of allBpbZIPgenes was upregulated or down regulated (p＜0.05)，except forBpbZIP37(group F)，and four genes were suppressed(BpbZIP06，BpbZIP13，BpbZIP17，andBpbZIP25) at various times (Fig.8).Moreover，at all sampling times during the salt treatment，BpbZIP07/ABF4，BpbZIP12，BpbZIP21/ABF2，andBpbZIP24were significantly upregulated，andBpbZIP06，BpbZIP13，andBpbZIP25were downregulated.The expression ofBpbZIP24and its homologous geneBpbZIP25in group D had opposite trends.Similar to the ABA treatment，all genes in theBpAREB/ABF/ABI5subfamily were significantly upregulated.In particular，the expression ofBpbZIP07/ABF4，the gene with the highest expression of all theBpbZIPgenes，was 33-fold higher than that of the control group (0 h).

Fig.8 Relative expression of 47 BpbZIP genes from B etula pendula at different sampling times during treatment with 150 mM NaCl.Error bars represent±SD of the means of three biological replicates ((ANOVA，*p＜0.05，**p＜0.01)

Fig.8 (continued)

During the drought treatment，the expression of all butBpbZIP32(group H) was upregulated or down regulated(p＜0.05)，and fvie genes (BpbZIP13，BpbZIP14，BpbZIP18，BpbZIP24，andBpbZIP37) were suppressed at various times(Fig.9).At six sampling times within 7 days，BpbZIP02，BpbZIP11，BpbZIP15，BpbZIP21/ABF2，andBpbZIP29were significantly upregulated (p＜0.01)，andBpbZIP13andBpbZIP37were significantly downregulated (p＜0.01).The expression level ofBpbZIP20，the most-upregulated gene，was 28-fold higher than that of the control group (0 h)，whileBpbZIP28had the lowest expression level at 24 h，only 7.8%higher than at 0 h.Different from other stresses，at day 7 of the drought，among the 41 upregulated genes，31 had maintained the upregulated trend (p＜0.05)，suggesting that theseBpbZIPgenes take part in regulating the response to the longer drought stress.

Fig.9 Relative expression of 47 BpbZIP genes from Betula pendula at different times during drought stress (7% v/v polyethylene glycol 6000).Error bars represent±SD of the means of three biological replicates ((ANOVA，*p ＜0.05，**p＜0.01)

Fig.9 (continued)

Heat treatment up regulated the expression of 37bZIPgenes (p＜0.05);eight genes were suppressed at different times (Fig.10).Significantly，at 1.5 h，the expression of 39 genes was rapidly and strongly upregulated or down regulate(p＜0.01).BpbZIP02(group H) andBpbZIP10(group I)did not change in their expression.BpbZIP05was rapidly and strongly upregulated at 1.5 h，16-fold higher than at 0 h.

Fig.10 Relative expression of 47 BpbZIP genes from Betula pendula at different times during heat stress (42°C).Error bars represent±SD of the means of three biological replicates ((ANOVA，*p＜0.05，**p＜0.01)

Fig.10 (continued)

These results showed that manyBpbZIPgenes responded positively to multiple stresses.The expression of allBpAREB/ABF/ABI5members were induced by all stress treatments.BpbZIP34/AREB3increased the most under heat stress，BpbZIP38/AREB1the most under drought stress.BpbZIP07/ABF4andBpbZIP21/ABF2were simultaneously and significantly upregulated by the drought and the salt treatments.

Discussion

InB.pendula，theBpbZIPgene family includes 47 fulllength genes，each of which contained a typical conserved bZIP domain.Sorghum and rice have morebZIPgenes than in birch，probably because，after the divergence of monocots and eudicots，the frequency ofbZIPgenes in eudicots was lower than that in monocotyledon (Wang et al.2011).

The length ofBpbZIPsequences significantly varied，but to a certain extent the length was conserved within the same group，indicating the high complexity inBpbZIPgenes.WoLF PSORT analyses showed that all but three of theBpbZIPgenes were predicted to be located in the nucleus，agreeing with previous studies showing that the basic region regulates nuclear localization.Interestingly，the three exceptionalBpbZIPgenes，BpbZIP23，BpbZIP42，andBpbZIP45，differed from their homologous genes.These results indicate that subcellular localization is irrelevant to the phylogenetic relationships and similarity of sequence，and homologous genes may be involved in different signal transduction pathways.

Comparative analyses between well-characterized species such asA.thalianaand those less-studied such asB.pendulacan be instrumental in the elucidating functional information on the less-studied.Therefore，determining phylogenetically orthologous relationships among bZIP proteins ofA.thaliana，P.trichocarpa，andB.pendulawill help birch researchers to predict the functions of protein more accurately.BpbZIPtranscription factors were found here to be phylogenetically closely related to the model forest treeP.trichocarpa.BpbZIPs were classified based on their phylogenetic relationship with AtbZIPs，and members that clustered in the same group also shared common motifs and conserved domains.Regarding sequence conservation of bZIP domains in BpbZIP proteins，only a few structural differences were found in the Leu zippers，while the basic region sequences were highly conserved，similar to the case for the bZIP domain in grapevine (Liu et al.2014).For example，the Leu zipper of BpbZIP06 was much shorter than that of BpbZIP46.The structural variations within these conserved regions may point to their functional specialization.By comparing additional conserved motifs of AtbZIPs and BpbZIPs in the same group，we found several similar motifs.For example，in group A，motifs 8，9，and 16 in birch corresponded to motifs 1，2，and 3 inA.thaliana，respectively.In group D，motif 3 in birch corresponded to motif 1 inA.thaliana.In group G，motifs 1，2，and 3 inA.thalianacorresponded to motif 6 in birch.In group E，motif 14 in birch corresponded to motif 2 inA.thaliana.In addition to its conformity to the grouping forA.thaliana，the classification grouping of BpbZIP proteins conformed to that found for maize (Wei et al.2012).In the present study，the findings for the conserved motifs，gene structure and intron patterns further supported the classification ofBpbZIPgenes (Figs.3，4).For the intron patterns in the bZIP domains，according to a previous study，the positions and phases can be considered be informative of the evolution of bZIP proteins (Nijhawan et al.2008).

To understand the roles ofBpbZIPgenes in various abiotic stresses，cis-elements (Fig.5) as well as expression levels under abiotic stresses (Figs.7，8，9，10) were analyzed.Various stress-relatedcis-elements were found，and almost allBpbZIPgenes were induced by stresses，confirming the potential transcriptional regulatory function ofBpbZIPgenes.Notably，BpbZIP09andBpbZIP12were simultaneously upregulated by ABA，salt，and drought，but they decreased significantly during the heat stress.This result indicates that there are two distinct sets of pathways by whichBpbZIPgenes respond to stresses，consistent with findings from grapevine and strawberry (Liu et al.2014;Wang et al.2017).However，BpbZIPgenes were coexpressed under ABA，salt，drought，and heat stress，suggesting that these genes may share similar roles under multiple stresses.We only investigated the effects of four stresses on apical buds of birch，but based on our results，we speculate that theBpbZIPgenes may take part in other biological processes not explored here，such as the responses of roots and stems to freezing or heavy metal exposure.Similarly，AtABF1andAtABF4ofA.thalianaare induced by cold(Choi et al.2000) and molybdenum application significantly increases the expression ofbZIPgenes in winter wheat (Sun et al.2009).

TheB.pendulagenome (435 Mb，14 chromosomes)is much larger than that ofA.thaliana(125 Mb，5 chromosomes) (Swarbreck et al.2007;Saloj?rvi et al.2017).Therefore，a smaller cutoffE-value of 1e-10 was selected to identify BpAREB/ABF/ABI5 subfamily members to avoid false positives.Subsequently，we confirmed the reliability of the by analyzing the conserved domains.It will also be meaningful to investigate the regulation of BpAREB/ABF/ABI5 protein phosphorylation in the future.

InA.thaliana，AtbZIP36/AREB1，AtbZIP37/ABF3，andAtbZIP38/AREB2are critical genes in the drought response via synergistic regulation of ABA signaling transduction(Yoshida et al.2010).InB.pendula，BpbZIP07/ABF4，BpbZIP21/ABF2，BpbZIP34/AREB3，andBpbZIP38/AREB1were also upregulated strongly in response to the drought stress (Fig.9).We also analyzed the somewhat differing patterns ofBpAREB/ABF/ABI5members in response to multiple stresses.Except forBpbZIP34/AREB3，allBpAREB/ABF/ABI5members were only slightly upregulated under heat stress，indicating that they are not very sensitive to heat，probably because heat stress-relatedcis-elements，such as heat-shock factor (HSF) elements，were absent from their promoters (Storozhenko et al.1998;Nishizawa-Yokoi et al.2009)，and the regulation of transcript abundance may rely only on the ABA signaling pathway activated by heat stress(Suzuki et al.2016).Considering these results，AREB/ABF/ABI5genes might serve as an ABA-response indicator forB.pendula.In addition to the familiar ABA-induced phosphorylation pathway，AtAREB/ABFgenes are rapidly induced under exogenous ABA treatment，thus offering an additional role in the ABA response for theAREB/ABF/ABI5subfamily(Furihata et al.2006;Yoshida et al.2015).Our results fromB.pendulashowed that exogenous ABA also increased the expression ofBpAREB/ABF/ABI5(Fig.6).

AtAREB/ABF proteins can form functional hetero-or homodimers，and AREB/ABF/ABI5 proteins in potato have different transcriptional activities than those in yeast cells(Yoshida et al.2010;Liu et al.2019b).Therefore，various dimer combinations of AREB/ABF/ABI5 proteins may be generated inB.pendula，so dimerization may provide another way that the transcriptional activity of AREB/ABF/ABI5 is altered and thus regulate the expression of target genes.

Conclusions

In summary，in this first comprehensive analysis of 47BpbZIPmembers in theB.pendulagenome，intron splicing phase，gene structure，and conserved motifs showed that theBpbZIPgenes were highly conserved in each clade and conform to the classification groups of BpbZIP proteins.Furthermore，our qPCR assays and analysis of promotercis-elements demonstrated that theBpbZIPgene family was widely involved in responses to ABA，salt，drought，and heat stress.These results will provide a platform for the study of other woody species responding to harsh environments.