Analysis of bZIP gene family in lotus (Nelumbo) and functional study of NnbZIP36 in regulating anthocyanin synthesis

Background The basic leucine zipper (bZIP) family is a predominant group of transcription factors in plants, involved in regulating plant growth, development, and response to stressors. Additionally, the bZIP gene family has a key role in anthocyanin production. Despite the significant role of bZIP genes in plants, their potential contribution in lotus remains understudied. Results A total of 124 bZIP genes (59 NnbZIPs and 65 NlbZIPs) were identified from genomes of two lotus species. These genes were classified into 13 groups according to the grouping principle of the Arabidopsis bZIP gene family. Analysis of promoter cis-acting elements indicated that most bZIP gene family members in lotus are associated with response to abiotic stresses. The promoters of some bZIP genes contain MYB binding sites that regulate anthocyanin synthesis. We examined the anthocyanin content of the petals from three different colored lotus, combined with transcriptome data analysis and qRT-PCR results, showing that the expression trends of NnbZIP36 and the homologous gene NlbZIP38 were significantly correlated with the anthocyanin content in lotus petals. Furthermore, we found that overexpression of NnbZIP36 in Arabidopsis promoted anthocyanin accumulation by upregulating the expression of genes (4CL, CHI, CHS, F3H, F3'H, DFR, ANS and UF3GT) related to anthocyanin synthesis. Conclusions Our study enhances the understanding of the bZIP gene family in lotus and provides evidence for the role of NnbZIP36 in regulating anthocyanin synthesis. This study also sets the stage for future investigations into the mechanism by which the bZIP gene family regulates anthocyanin biosynthesis in lotus. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-023-04425-2.


Background
The basic region/leucine zipper motif (bZIP) transcription factor is one of the largest gene families of transcription factors, found in almost eukaryotes, and is highly conserved [1,2].The bZIP transcription factor consists of a basic region and a leucine zipper structure [3,4], where the basic region is relatively conserved and consists of approximately 20 amino acid residues.The residues are able to bind to specific DNA sequences with an ACGT core such as a-box (TAC GTA ), C-box (GAC GTC ) and G-box (CAC GTG ) [3,5].The bZIP transcription factors play important roles in plant growth, development and the response to external environment, such as regulation of plant growth [6], flower development [7][8][9], seed maturation and dormancy [10], senescence [11], light signaling [12,13], damage [14] and response to various environmental stresses [15][16][17].
With the development of high-throughput sequencing technology, a large amount of plant genomic data has been published and the bZIP gene family has been identified and studied in many plants, such as Arabidopsis [3,18], rice (Oryza sativa) [19], soybean (Glycine max) [20], tomato (Solanum lycopersicum) [21], apple (Malus domestica) [22], etc. Involvement in response to various abiotic stresses/biotic stresses is a very typical and also the most reported function of bZIP transcription factors, such as salt stress, drought, high temperature, cold stress and pathogen infection [12,19,23].For example, soybean contains 131 members of the bZIP gene family.More than 1/3 of them extensively involved in defense response including ABA, salt, drought, and cold stresses [24].In Arabidopsis, AtbZIP17 and AtbZIP60 can enhance salt tolerance [25,26].AtbZIP1 can be promoted in Arabidopsis in response to abiotic stresses such as drought and high salt stresses [27].Preeti et al. identified 191 bZIP transcription factors in wheat and found that TabZIPs were may play a role in various stress (high temperature, drought, salt stress) relief mechanisms.Overexpression of TabZIPs in Arabidopsis enhanced the ability of transgenic Arabidopsis to tolerate salinity, drought, high temperature, and oxidative stresses [28].In recent years, with the in-depth study of bZIP transcription factors, it has been found that bZIPs are involved in regulating the synthesis of plant secondary metabolites in addition to responding to abiotic stresses and participating in the regulation of plant growth and development.It has been shown that some bZIP genes are involved in the regulation of anthocyanin synthesis, such as HY5 (Elongated hypocotyl 5) and HYH (HY5 homolog) [29][30][31].In Arabidopsis, HY5 can promote anthocyanin accumulation by interacting with the promoters of some MYB transcription factors or directly with MYB transcription factors to promote the expression of genes related to anthocyanin synthesis [30].Similarly, in tomato and apple, HY5 has a similar function in the regulating anthocyanin biosynthesis [29,31].Tu et al. demonstrated that VvbZIP36 is a negative regulator of anthocyanin synthesis using CRISPR/Cas9 technology.VvbZIP36 plays a role in balancing the synthesis of stilbene (α-glucosin), lignans, flavonols and anthocyanins [32].This evidence suggests that bZIP plays an important role in the regulation of the plant flavonoid synthesis pathway.
Lotus (Nelumbo) is one of the ten traditional Chinese flowers with high ornamental value [33].Lotus are rich in flavonoid substances, which provide them with rich floral color and high medicinal value [34].It is of great significance to study of flavonoids biosynthesis in lotus to improve the medicinal value and enrich the flower color of lotus.However, most of the reports related to anthocyanin biosynthesis of lotus have focused on MYB-bHlH-WD40 complex, while there is no report on the regulation of anthocyanins in lotus by bZIP.In this study, 59 NnbZIP genes and 65 NlbZIP genes were identified from the genomes of Asian lotus and American lotus, respectively.The chromosomal distribution, gene structure, conserved motifs and evolution of bZIP gene family were further investigated.The regulation of lotus bZIP genes on anthocyanin synthesis in lotus was also analyzed.The NnbZIP36 gene was cloned as a representative bZIP gene, the regulation of anthocyanins by NnbZIP36 was also studied by ectopic expression in Arabidopsis.These results provide a reference for studying the regulation of anthocyanin synthesis by the bZIP gene family in lotus, and provide a basis for further elucidating the mechanism of lotus flower color formation.

Acquisition of sequencing data
Genomic data for the lotus of Nelumbo nucifera (N.nucifera) and Nelumbo lutea (N.lutea) used in this study are available in the databases (http:// nelum bo.cngb.org/ nelum bo/ home) and NCBI (https:// www.ncbi.nlm.nih.gov/) respectively under BioProject number PRJNA747731.The transcriptome data associated with this study were downloaded from the NCBI's SRA database (https:// www.ncbi.nlm.nih.gov/ sra/?term=) and the transcriptome data numbering information is shown in Table S1.

Calculation of transcriptome gene expression
The transcriptome cleaned data were aligned to the lotus reference genome using hisat2 software [35][36][37].Gene expression levels were calculated using the R package DESeq2 [38].

Identification of bZIP genes in the genomes of N. nucifera and N. lutea
The Hidden Markov Models (HMM) files containing sequences of conserved domains of bZIP_1 (PF00170), bZIP_2 (PF07716), bZIP_C (PF12498), bZIP_Maf (PF03131), and HLH (PF00010) were downloaded from the Pfam database (http:// pfam.xfam.org).Based on the Hidden Markov Model (HMM) files, the bZIP gene sequences of N. nucifera and N. lutea were queried using HMMER software.The bZIP genes in both genomes of lotus were further identified based on 79 AtbZIP protein sequences from Arabidopsis using BLASTp software.Conserved structural domain analysis of the bZIP protein sequences obtained was performed using the Simple Modular Architecture Research Tool (SMART: http:// smart.embl-heide lberg.de) to ensure the validity of the identified bZIP genes.

Characterization of the bZIP TFs in lotus
According to the gff file information of N. nucifera and N. lutea genomes, the length and position information of the lotus bZIP genes were extracted.Molecular weight (MW) and isoelectric point (pI) of the lotus bIZP genes were calculated using the ProtParam tool in ExPASy Server (https:// web.expasy.org/ protp aram/).The information of the structure, location and conserved components of the NnbZIPs and NlbZIPs genes were examined using the MEME online tool, identifying conserved motifs shared by the bZIP proteins.The MAST xml files were downloaded and visualized using TBtools software.The chromosomal distribution and gene structures (including introns and utr) of NnbZIPs and NlbZIPs were analyzed using TBtools [39].

Phylogenetic analysis of the NnbZIPs and NlbZIPs protein
The amino acid full-length sequences of bZIP proteins from seven species were aligned using muscle software with default parameters, and a neighbour-joining (NJ) phylogenetic tree was constructed using MEGA 7.0 software with a bootstrap replication number of 1000 [40].The bZIP transcription factors were classified into different groups according to the topology of the phylogenetic tree.
In order to study the expression characteristics of bZIP genes in N. nucifera and N. lutea, the 2000 bp sequence upstream of the start codon (ATG) of bZIP genes were obtained as the promoter region, and the phylogenetic tree of the promoter sequence was constructed using MEGA 7.0 software with performed 1000 bootstrap replications.

Analysis of the cis-acting element of the promoter of the bZIP genes in lotus
The obtained bZIP gene promoter sequences were uploaded to the PlantCARE database (http:// bioin forma tics.psb.ugent) for cis-acting regulatory element prediction.Visualized using TBtools based on the database analysis results.

Identification of gene duplication patterns and covariance analysis
The MCScan (https:// github.com/ tangh aibao/ jcvi/ wiki/ MCscan-Python) software was used to analyze the gene duplication events of NnbZIP and NlbZIP members [41], and the gene duplication relationship was visualized using TBtools software.The TBtools software was used to estimate the non-synonymous substitution rate (Ka), the synonymous substitution rate (Ks) and their ratio (Ka/Ks) [39].
The python version of MCscan (JCVI v1.1.7)[42] was used to perform a comparative analysis of the genomes of six plants.Co-linear modules with a co-linear relationship to the lotus bZIP genes were highlighted.

Experiment material
The experimental materials used in this study were collected at Baima Base for Teaching and Scientific Research of Nanjing Agricultural University.The experimental material cultivars were: the ancient Chinese lotus and the American yellow lotus.The samples were frozen in liquid nitrogen and stored in a -80 °C refrigerator.Wild-type Arabidopsis plants were grown in incubators (Arabidopsis seeds were kept in our laboratory).

RNA isolation, cDNA synthesis and quantitative real-time PCR analysis
The total RNA was extracted using the FastPure Universal Plant Total RNA Isolation Kit (Novazyme Biotechnology, China), and the experimental steps were referred to the kit instructions.The integrity of total RNA was detected by 1% agarose gel electrophoresis, and the concentration was determined using NanoDrop 1000 (Thermo, USA).The first-strand cDNA was synthesized using HiScript III RT SuperMix for qPCR (gDNA wiper) kit (Novazyme Biotechnology, China), and 1 μg of total RNA was used for each 20 μL reaction.The cDNA product was diluted tenfold with deionized water before use.
qRT-PCR experiments were performed using ChamQ SYBR qPCR Master Mix kit (Novozymes Biotechnology, China) on a CFX96 Touch ™ Real Time PCR Detection System (BIORAD, USA). 2 μl of diluted cDNA was used for each reaction and other reaction components and conditions were carried out according to the manufacturer's instructions.Specific primers were designed for the qPCR (Table S2).The Actin gene of lotus [43] was used as an internal reference gene and was performed in three biological and three technical replicates of each treatment.Relative expression was calculated using the 2 −ΔΔCt method [44,45].

Gene cloning of NnbZIP36 and its heterologous expression in Arabidopsis
The complete transcript sequence of the NnbZIP36 gene was cloned using 2 × Hieff Canace ® Plus PCR Master Mix (With Dye) (Yeasen Biotech Co.).PCR reactions were carried out in a volume of 50 μL under the following conditions: denaturation at 98 °C for 2 min, 40 thermal cycles (98 °C / 10 s; 60 °C / 20 s; 72 °C / 30 s) and a final extension at 72 °C for 5 min.The PCR products were separated by 1% agarose gel electrophoresis and the gel was cut to recover the target fragments.After recovery, ligation, transformation and sequencing, the final CDS sequence of the NnbZIP36 gene was obtained.
The CDS sequence of the NnbZIP36 gene was assembled in the overexpression vector pFAST-R05 using homologous recombination approach.The correct recombinant plasmid obtained by sequencing was transferred into Agrobacterium tumefaciens GV3101 using freeze-thawing method.It was inoculated in 50 mL of LB liquid medium containing 100 mg/L kanamycine and 50 mg / L rifampicin, incubated at 200 r / min for 48 h at 28 °C, and then centrifuged at 4000 rpm for 10 min to collect the bacteria.Bacteria were resuspended with an osmotic agent (0.05% sliwet77 + 5% sucrose + 1/2 MS liquid medium).The Arabidopsis plants were subsequently transformed according to the flower dip method and incubated in the dark for 48 h.After three infestations, Arabidopsis seeds were collected and screened for NnbZIP36 transgenic positive plants [46].The positive plants were identified using PCR and qRT-PCR methods.Seeds of the identified positive plants were grown in 1/2 MS Petri dishes and positive seedlings with two true leaves were photographed after two weeks.The primers used for cloning and vector construction are listed in Table S2.

Statistical analysis
Data in this study are shown as mean ± standard error (SE) of 3 or 6 independent biological replicates.Statistical differences between samples were analyzed by LSD and Duncan (D) (p < 0.05).Data analysis and visualization were processed using SPSS 20.0 and GraphPad Prism 8.0.

Phylogenetic classification and analysis of bZIP TFs in lotus
The bZIP protein sequences of two species of lotus, Ginkgo, Amborella, Nymphaea colorata, Vitis vinifera and Arabidopsis were aligned and a phylogenetic tree was constructed using MEGA 7.0 software.The bZIP proteins were classified into 13 groups (A, B, C, D, E, F, G, H, I, J, K, M  and S) according to the grouping rules of Arabidopsis bZIP gene family (Fig. 1).The number of bZIP genes distributed among the 13 groups are differed considerably, with groups A, D, I and S possessing a higher number of bZIP genes.Group J had the fewest members with only six members, and the bZIP members of N. nucifera and Amborella were present.In addition, the bZIP gene families of these seven species were distributed in the remaining 12 groups.The results showed that most of the NnbZIP genes and NlbZIP genes have a high degree of similarity.However, there are still NnbZIP genes and NlbZIP genes that specific to each other.It indicated that the N. nucifera and N. lutea have evolved new bZIP genes after a long period of isolation in order to adapt to the local natural environment.

Structural analysis of the bZIP gene family in lotus
Gene structure is an important reference information for measuring the evolution of gene families which can further support the results of phylogenetic trees [50].Genetic structure analysis showed that most of the genes of NnbZIP and NlbZIP were structurally similar in the phylogenetic tree.Only a small number of NnbZIP and NlbZIP genes had large differences between them, and these genes are specific to the presence of N. nucifera and N. lutea, respectively (Fig. 2).In addition, bZIPs can be divided into two types according to the presence or absence of introns.Here, a total of 18 bZIP genes have no introns.15 of the 18 bZIP genes are classified within group S, while NnbZIP17, NnbZIP19 and NnbZIP54 are classified in group E. The protein conserved motif analysis showed that there were 9 conserved motifs in the bZIP gene family in lotus.Motif5 is the most conserved motif and is present in all bZIP genes.Furthermore, Motif1 was present in all but NnbZIP13, NnbZIP17, NnbZIP27, NlbZIP16 and NlbZIP30.It is speculated that Motif1 and Motif5 may be associated with conserved structural domains of the bZIP gene family.The conserved structural domains of the lotus bZIP gene family were identified using the SMART database.Model analysis showed that the bZIP gene family has three types of conserved structural domains, with the largest number of bZIP gene members possessing the BRLZ structural domain.The bZIP_C structural domain contains 4 members (NnbZIP13, NnbZIP27, NlbZIP16 and NlbZIP30).The bZIP_2 structural domain only has one member, NnbZIP17 (Fig. 2).

The cis-acting elements of the bZIP gene promoter in lotus
Prediction and analysis of cis-acting elements in gene promoters can help predict their function [51].To understand the function of the bZIP gene family of lotus, the 2000 bp sequences upstream of the ATG of 124 bZIP genes (59 NnbZIPs and 65 NlbZIPs) were analyzed and a large number of cis-acting elements were identified (Fig. S2).In order to facilitate further analysis, these cisacting elements were divided into four categories: light response, stress response, hormone response, growth and development-related elements (Fig. 3).Among them, light response elements were widely distributed in the promoters of all bZIP genes.With the exception of NlbZIP02 and NnbZIP09, stress response-related elements are present in the promoters of all remaining bZIP genes.And they are generally more numerous compared to growth and development-related elements.This indicated that the bZIP gene family of lotus may play an important role in stress resistance of lotus.In addition, hormone response-related elements are also widely present in the promoter of the lotus bZIP genes.Interestingly, 9 NlbZIP genes and 8 NnbZIP genes have a large number of ABA-responsive elements in their promoters, so we speculate that these genes have similar functions to ABA-responsive factors (e.g.ABI3, ABI5, etc.).The large distribution of MYB-binding sites and MYC-binding sites also enhances the role of the bZIP genes in regulating the growth and development of the lotus and coping with the external natural environment.In conclusion, our results suggested that the lotus bZIP genes may respond to different signaling pathways through different types of cis-acting elements within their promoter regions.

Collinearity analysis of bZIP genes in lotus
MCScanX software was used to analyze the gene duplication events of NlbZIP genes and NnbZIP genes.A total of 24 gene duplication events were identified in N. lutea, and 23 gene duplication events in N. nucifera (Fig. 4AB).These gene duplication events were all segmental duplication events, and no tandem duplication events occurred.The Ka and Ks of these gene duplication events were analyzed.The synonymous substitution rate distribution of the NlbZIP and NnbZIP duplication gene pairs were calculated.It was found that the Ks of the NlbZIP genes had a peak between 0.5 and 0.6.In contrast, the NnbZIP genes showed two peaks distributed at Ks between 0.4-0.5 and 0.6-0.7,respectively (Fig. 4C).The Ka/Ks ratios of all bZIP-replicated gene pairs were less than 1 (Fig. 4D).These results indicated that these genes underwent strong purification selection after fragment replication, and the function of these genes did not significantly differentiate.
To further understand the evolution of the lotus bZIP gene family, the covariance of two lotus species Fig. 1 Phylogenetic tree of bZIP genes in seven species.The protein sequences of the bZIP TFs were compared using the maximum likelihood (ML) method.These proteins were divided into 12 groups and each group was assigned a different color with four other species (Nymphaea colorata (N.colorata), Vitis vinifera (V.vinifera), Glycine max (G.max) and Arabidopsis) were aligned using a python version of the MCScanX software (Fig. 5).Each of these four species represents an important point in the evolution of angiosperms.The results show that there is little difference in the number of co-linear pairs between the two lotus species and the N.colorata and V.vinifera.The number of homologous genes for bZIP was significantly increased on the G.max side of the lotus compared to G.max.The number of homologous genes on the Arabidopsis side was the lowest for these species.These results indicate that bZIP genes have experienced multiple expansion and contraction processes during the evolution of angiosperms, which may be the result of plants adapting to the changes in their natural environment.

Alignment analysis of the NnbZIP genes with the NlbZIP genes
The protein sequences of NnbZIPs and NlbZIPs were aligned using BLATp software.The results showed that the protein sequences of 53 bZIP genes were highly similar between N. nucifera and N. lutea (Fig. S3).There are 10 bZIP genes specifically present in N. lutea and 6 bZIP genes specific to N. nucifera (Fig. 6).This indicates that most of the functions of the bZIP gene family are similar between N. nucifera and N. lutea.The few unique bZIP genes may promote the formation of the characteristic traits of N. nucifera and N. lutea.This may be the result of the evolution of N. nucifera and N. lutea to adapt to the local environment.The gene information of unique to N. nucifera and N. lutea is represented in Table 2.

The bZIP gene family involved in the regulation of anthocyanin biosynthesis in lotus
In order to understand the regulatory effect of bZIP gene family on lotus anthocyanins, the transcriptome data from three lotus cultivars with different flower colors were downloaded and reassembled.The anthocyanin content of three lotus cultivars with different flower colors were detected.The results showed that anthocyanin was only accumulated in red lotus (Fig. 7B), and the content of anthocyanin in red lotus was highest in late bud stage (Fig. 7D).Further transcriptome data showed that the individual bZIP gene members of N. nucifera and N. lutea did not seem to be directly related to the anthocyanin content of lotus.Interestingly, NlbZIP38, NlbZIP50 and NnbZIP36 were highly expressed in the red lotus cultivar 'Jinlinghuodu' , while lower expressed in the white lotus 'Baiyinlian' and yellow lotus 'Jinsenianhua' (Fig. 7A).Further research found that the expression levels of NlbZIP38 and NnbZIP36 were consistent with the change trend of anthocyanin content in red lotus cultivar (Fig. 7C).This indicated that NlbZIP38 and NnbZIP36 may be involved in regulating the synthesis of lotus anthocyanins.The protein sequence comparison analysis showed that the protein sequences of NlbZIP38 and NnbZIP36 had a very high similarity (Fig. S4).NlbZIP38 and NnbZIP36 were considered to have the same function.Therefore, NnbZIP36 was selected as a candidate gene for further verification in this paper.

NnbZIP36 promotes anthocyanin accumulation in Arabidopsis
The result of qRT-PCR verified that NnbZIP36 was only highly expressed in red lotus cultivar (Fig. 8A).To understand the regulation mechanism of NnbZIP36 on anthocyanin synthesis, the open reading frame (ORF) of NnbZIP36 gene was cloned (Fig. 8B).The pFAST-R05-NnbZIP36 overexpression vector was constructed.Three NnbZIP36-OE transgenic Arabidopsis lines were obtained by PCR and qRT-PCR experiments (Fig. 9BC).It was found that the petioles and leaves of the transgenic Arabidopsis were red (Fig. 9A).The anthocyanin content of NnbZIP36-OE Arabidopsis plants was significantly up-regulated (Fig. 9D).The expression of most structural genes (4CL, CHI, CHS, F3H, F3'H, DFR, ANS and UF3GT) on the anthocyanin biosynthetic pathway were significantly up-regulated in transgenic Arabidopsis compared with wild-type Arabidopsis (Fig. 10).These results showed that the NnbZIP36 gene has a positive regulatory effect on the accumulation of anthocyanins in transgenic Arabidopsis.

Discussion
bZIP TFs are a family of transcription factors that widely exist in eukaryotes and highly conserved in evolution [18].As one of the largest transcription factor families, bZIP TFs play important roles in regulating plant growth, development and abiotic stress [7,10,30].The bZIP family genes on a genome-wide scale have been systematically analyzed in many plant species, including some key  crops, such as: Arabidopsis (79) [18], Oryza sativa (86) [19], Vitis vinifera (55) [52], Solanum lycopersicum (69) [21], Actinidia chinensis (81) [53], Punica granatum (65) [54], etc.However, until now, although the entire lotus genome has been sequenced [35,36], the studies of bZIP gene family in lotus are limited.In this study, the bZIP gene family from genomes of two lotus species were identified and systematically analyzed.The potential functions of bZIP gene family were explored in the regulation of anthocyanin biosynthesis in lotus.
In this study, a total of 124 bZIP genes (59 NnbZIPs and 65 NlbZIPs) were identified in the genomes of two species of lotus (Table 1).It is similar to the number of bZIP gene family members in Vitis vinifera (55) [52], Punica granatum (65) [55] and other species.The bZIP genes in lotus were classified into 13 groups by constructing phylogenetic trees of the bZIP gene family in two species of lotus, Ginkgo, Amborella, Nymphaea colorata, Vitis vinifera and Arabidopsis, which is similar to the grouping of Arabidopsis bZIP gene family [18] (Fig. 1).In addition to the absence of bZIP gene members of N. nucifera and Amborella in group J.The bZIP gene members of seven species were distributed in the other 12 subgroups.This suggests that the number of bZIP genes is more conserved during evolution.Due to the highly conserved nature of bZIPs sequences, the bZIP genes in the same group have the same or similar functions, which provides a reference for studying the functions of this gene family [55].
Analysis of gene structure and conserved motifs can help us predict the function of genes [56].It is also an important basis for studying gene evolution and gene replication [55].The results of structural analysis of lotus bZIP genes showed that the bZIP gene structure of lotus is relatively simple, with the number of introns ranging from 0 to 11.It can be divided into two categories according to the presence or absence of introns.A total of 18 bZIP genes without introns were found in the genomes of two lotus species.Of these 18 bZIP genes, 15 were classified in group S and three bZIP genes (NnbZIP17, NnbZIP19 and NnbZIP54) were classified in group E. A similar situation exists in species such as Arabidopsis, strawberry and gourd [18,52,57].bZIP genes located in different subgroups may play different roles in plant growth and development.For example, most members of group A are involved in the ABA biological pathway and regulate plant responses to abiotic stresses [58][59][60].Genes located in group D can be involved in plant defense against pathogens [3].Genes located in group G and their homologs are mainly involved in blue-violet signaling [61].The lack of introns contributes to an accelerated post-transcriptional response to abiotic stresses [62].Similar findings were observed in soybean and watermelon, revealing evolutionary imprinting and functional The number of conserved motifs of each bZIP protein also varies greatly.There are many bZIP proteins with only 1 protein conserved motif, and the most have 6 protein conserved motifs (Fig. 2).It has been shown that the bZIP gene family has a variety of different conserved structural domains.In addition to the specific bZIP_1 (PF00170) and bZIP_2 (PF07716) domains, bZIP proteins in plants also have additional functionally conserved domains that participate in various biological processes [67].In lotus, the identified bZIP members have three different conserved domains: BRLZ, bZIP_C and bZIP_2.Among them, the BRLZ conserved domain is widely present in the bZIP gene family of lotus.bZIP_C is found in 6 members (NnbZIP26, NlbZIP29, NlbZIP02, NlbZIP18, NnbZIP27 and NlbZIP30), while bZIP_2 is only present in NnbZIP17 (Fig. 2).
The cis-acting element of the promoter plays an important role in the transcription of the gene, its activity can directly affect the gene expression and function [68].The function of the gene can be predicted according to the type of the cis-acting element in the gene promoter.Among the promoters of the 124 bZIP gene family members in lotus, there are generally more cis-acting elements related to stress resistance than those related to growth and development (Fig. 3).Implying that the bZIP gene family of lotus plays an important role in response to stress.A large number of response hormone-related elements (especially ABA response elements and MeJA response elements) also widely exist in the promoters of bZIP genes.Studies have shown that the transcription and regulation of many bZIP genes are induced by ABA  The expansion and contraction of gene families are the result of plant adaptation to different environmental conditions, while the tandem and fragment duplication events of gene families are crucial for the expansion of gene families and the diversification of gene functions [70].The results of the MCScanX software analysis showed that the 24 pairs and 23 pairs of bZIP fragment duplication genes were detected in N. lutea and N. nucifera, respectively, and no tandem duplication events were detected in two lotus species (Fig. 4AB).The Ka/Ks ratios calculated for all gene pairs were less than 1 (Fig. 4D).This implies that these genes may have undergone strong purifying selection pressures during evolution and did not significantly alter the function of the bZIP genes.Furthermore, we compared the collinearity of the two lotus species with four other species (N.colorata, V.vinifera, G.max and Arabidopsis).The number of collinear pairs between the two kinds of lotus and N.colorata or V.vinifera is not much different.There are more collinear pairs between lotus and G.max, while the number of collinear pairs between lotus and Arabidopsis is greatly reduced (Fig. 5).This shows that lotus has a relatively close evolutionary relationship with N.colorata and V.vinifera.However, there is a relatively distant evolutionary relationship with G.max and Arabidopsis.It also indicated that the bZIP gene family experienced multiple expansion and contraction events during the evolution of angiosperms.
As a water-soluble pigment involved in the secondary metabolic pathway of plants, anthocyanins play an important role in providing bright colors to plants and helping them cope with the external environment [71].In recent years, members of the bZIP gene family in many plants have been reported to be involved in the regulation of anthocyanins.For example, HY5 in Arabidopsis can promote the expression of genes related to anthocyanin synthesis.The accumulation of anthocyanins by binding to the promoter of MYB transcription factors or directly interacting with MYB transcription factors [10,13,30].In apple, the bZIP transcription factor gene MdHY5 can actively regulate the accumulation of anthocyanins by directly promoting the expression of MdMYB10 and MdMYB1 genes [29,72,73].RsbZIP011 and RsbZIP102 were reported to be associated with anthocyanin content in radish [74].PgbZIP16 and PgbZIP34 of pomegranate were overexpressed in Arabidopsis to promote anthocyanin content accumulation in Arabidopsis leaves [55].VvbZIP36 was a negative regulator of anthocyanin synthesis in grapes and played a role in balancing the synthesis of stilbene (α-glucosin), lignans, flavonols, and anthocyanins [32].These clues indicated that the bZIP gene family plays an important role in regulating anthocyanin synthesis and promoting anthocyanin accumulation in plants.
Anthocyanin is one of the main pigment substances that affect the flower color of lotus.Understanding the regulation mechanism of anthocyanin synthesis in lotus is of great significance for improving ornamental value and stress resistance [34].In the present study, we found that NnbZIP36 and its homologous gene NlbZIP38 in lotus were significantly correlated with the anthocyanin content of lotus.The results of qRT-PCR further showed that NnbZIP36 was highly expressed in red lotus cultivar of 'Jinlinghuodu' .This implies that NnbZIP36 has the ability to regulate the synthesis of anthocyanins in lotus.In order to study the function of NnbZIP36, we constructed the pFAST-R05-NnbZIP36 overexpression  vector and transformed Arabidopsis.The results showed that the leaves and petioles of Arabidopsis overexpressed with NnbZIP36 appeared red and the anthocyanin content was significantly increased (Fig. 9).Most genes in the anthocyanin synthesis pathway (4CL, CHI, CHS, F3H, F3'H, DFR, ANS and UF3GT) were significantly up-regulated in NnbZIP36 transgenic Arabidopsis.The expression level of PAL in Arabidopsis overexpressing NnbZIP36 was lower than that of wild-type Arabidopsis (Fig. 10).This may be because PAL is a structural gene involved in the synthesis of several secondary metabolites but not directly involved in the synthesis of anthocyanins [75].According to the basis of existing experimental studies, we have found that NnbZIP36 can promote the anthocyanins synthesis in overexpressed Arabidopsis.However, the mechanism of NnbZIP36 gene regulating anthocyanin synthesis in lotus needs further study.And the NnbZIP36 interacts with other gene to co-regulate in flower color formation of lotus will be a new direction for our research.

Conclusions
In this study, we identified a total of 124 bZIP genes (59 NnbZIPs and 65 NlbZIPs) from the genomes of two lotus species.The bZIP genes of lotus were divided into 13 groups by constructing the phylogenetic tree of seven species.Due to the high degree of conservation of the bZIP genes, proteins with similar functions were clustered into a group.The combination of gene structure, conserved motifs and promoter cis-element analysis provided a reliable basis for studying the functions of related genes in the lotus bZIP gene family.Candidate genes affecting anthocyanin biosynthesis were identified using transcriptome data analysis.The gene cloning, qRT-PCR analysis, and Arabidopsis transformation were performed.Our results suggested that NnbZIP36 has a promoting role in anthocyanin accumulation.In our future studies, we will focus on the regulatory network comprising NnbZIP36 and anthocyanin synthesis-related genes, and elucidate the molecular mechanism by which NnbZIP36 interacts with other transcription factors to promote anthocyanin accumulation.

(
See figure on next page.)Fig. 2 Gene structure analysis of the bZIP gene family in the lotus.A: Structure of the bZIP gene family in N. nucifera.B: Structure of the bZIP gene family in N. lutea.Exons and introns are shown using green bars and grey lines.bZIP members in protein motifs.The colored boxes depict the different patterns.Clustering is based on the results of phylogenetic analysis

Fig. 3
Fig. 3 Cis-acting elements of the promoter region of bZIP genes in lotus.A: Number of different classes of cis-acting elements in the promoter region of the N. lutea NlbZIP genes.B: Number of different classes of cis-acting elements in the promoter of the N. nucifera NnbZIP genes [69].Through ABA response elements (ABRE), bZIP genes are involved in coordinating responses to drought and other environmental factors as well as seed development processes [69].It proves the role of bZIP gene family in regulating plant stress resistance.It is further shown that the cis-acting element of the promoter plays an indispensable role in the response of the bZIP gene to external environmental changes.In addition, a large number of MYB and MYC binding sites expanded the regulatory network of the bZIP gene family.

Fig. 4
Fig. 4 Covariance analysis and Ks distribution of bZIP replication genes in lotus.A: Covariance analysis of bZIP replication genes in N. lutea.B: Covariance analysis of bZIP replication genes in N. nucifera.C: Ks distribution of bZIP replication genes in two lotus species.D: Ka/Ks ratios of the two lotus species

Fig. 5
Fig. 5 Synteny analysis of bZIP genes between two lotus species with N.colorata, V. vinifera, G. Max, and Arabidopsis.The green lines indicate covariance between N. luctea and the other four species, the blue lines indicate covariance between N. nucifera and the other four species

Fig. 7 Fig. 8
Fig. 7 Potential ability of the bZIP genes to regulate anthocyanin synthesis in lotus.A: Expression profiling of the bZIP genes in lotus cultivars of different flower colors.B: Detection of anthocyanin content in different lotus cultivars.C: Expression profiling of the bZIP genes in different opening periods of red lotus cultivar.D: Detection of anthocyanin content in different opening periods of red lotus cultivar.(In the pictures, a indicates the lotus cultivar 'Jinsenianhua' , b indicates the lotus cultivar 'Jinlinghuodu' , c indicates 'Baiyinlian' , A indicates the bud period of the ancient lotus, B indicates the late bud period of the ancient lotus, C indicates the early blooming period of the ancient lotus, and D indicates the full blooming period of the ancient lotus)

Table 1
The detailed characteristics of bZIP genes identified in lotus

Table 2
The specifical bZIP genes in two lotus species