Genome-Wide Identification of bZIP Gene Family in Lilium davidii var. unicolor and Expression Analysis during Dormancy Release and Regeneration

1 Introduction

Lilium spp., belonging to the family Liliaceae, comprises more than 100 wild species of perennial bulbous plants that are widely distributed across the Northern Hemisphere [1,2]. China is one of the major centers of Lilium diversity, harboring about 55 species and 18 varieties [2]. Owing to their diverse floral morphology, attractive fragrance, and edible or medicinal bulbs, lilies have considerable economic importance in the global floriculture industry, as well as in traditional medicine and food production. Among them, L. davidii var. unicolor, commonly known as Lanzhou lily, is a botanical variety of L. davidii and is the only traditionally cultivated sweet edible lily in China [2,3]. With a cultivation history of more than 150 years, Lanzhou lily has become an important local crop [2].

In commercial production, lily yield and quality are strongly constrained by two physiological bottlenecks: bulb dormancy and asexual regeneration. As the central organ for both nutrient storage and vegetative reproduction [4], the bulb is subject to tight regulation by endogenous signals and environmental cues. The transition from dormancy maintenance to release strictly dictates phenological phases and year-round production schedules [5]. Concurrently, since commercial propagation relies heavily on asexual reproduction, reproductive efficiency is intrinsically linked to the regeneration capacity of plants [4]. At the molecular level, both dormancy release and asexual regeneration are characterized by profound transcriptional reprogramming [6,7]. Previous studies have highlighted phytohormones and carbohydrates were known to play an important role in regulating the bulb dormancy and asexual regeneration. In L. brownii var. viridulum, a staged temperature regime of 25°C–15°C–4°C promoted dormancy release and subsequent bulb enlargement more effectively than direct transfer from 25°C to 4°C. Integrated metabolomic and proteomic analyses associated this response with the early accumulation of soluble sugars and increased abundance of glycolytic enzymes, indicating that temperature history may prime carbon remobilization before visible bud growth [8]. CYCLING DOF FACTOR 2 (LdCDF2) functioned as a central hub in a gibberellic acid (GA) self-amplifying regulatory network, where it drived GA biosynthesis but was dually restricted by the DELLA protein LdSLR1, forming a dynamic feedback loop that precisely controlled lily bulb germination [9]. Progress in bulblet regeneration has provided more direct functional evidence. Genome-enabled transcriptomic analyses identified expanded gene families associated with starch and sucrose metabolism in lily and demonstrated that manipulating XYLOGLUCAN ENDOTRANSGLUCOSYLASE/HYDROLASE (LaXTH) expression significantly altered the induction of scale-derived bulblets [1]. More recently, spatial transcriptomics resolved nine major cell layers during bulblet formation and identified the type-B cytokinin RESPONSE REGULATORS (LiRR12/14) as positive regulators of both regeneration and expansion. LiRR14 activated CELL WALL INVERTASE 4 (LiCWIN4) to promote sucrose hydrolysis during regeneration, whereas LiRR12 and LiRR14 subsequently activated SUCROSE SYNTHESIS (LiSPS3) and STARCH BIOSYNTHESIS (LiSSS1), respectively, to support sucrose biosynthesis and starch accumulation during bulblet expansion [10].

The basic leucine zipper (bZIP) transcription factor family is one of the most widespread and evolutionarily conserved transcription factor families in plants [11]. Members of this family are characterized by a conserved bZIP domain that generally consists of two functional regions: a basic DNA-binding region and a leucine zipper dimerization region. The basic region typically contains a conserved N-x7-R/K motif and mediates the recognition of specific cis-regulatory elements in target promoters. In many bZIP proteins, this region is also associated with nuclear localization. The leucine zipper region contains regularly spaced hydrophobic leucine residues, which facilitate homo- or heterodimer formation and are required for the transcriptional regulatory activity of bZIP proteins [12].

The bZIP transcription factors participate in diverse biological processes of plants, including seed germination [13], organ differentiation, embryogenesis [14], photomorphogenesis, cell elongation [15], floral development [16], hormone signaling [13], secondary metabolite biosynthesis [17], and responses to biotic [18] and abiotic stresses [19]. With the rapid development of genome sequencing and annotation, genome-wide identification and functional analysis of bZIP gene families have been reported in many plant species, including Arabidopsis thaliana [12], rice [20], papaya [21], Chinese cabbage [19], cotton [22], and Eutrema salsugineum [23]. These studies have provided important insights into the evolutionary conservation, structural diversification, and functional specialization of bZIP genes in plants.

Accumulating evidence suggests that bZIP transcription factors act as key regulatory nodes in plant growth and developmental signaling pathways. In Arabidopsis, ELONGATED HYPOCOTYL 5 (HY5), a central positive regulator of photomorphogenesis whose loss of function results in elongated hypocotyls, interacted with HISTONE DEACETYLASE 15 (HDA15) to repressing hypocotyl cell elongation, thereby promoting photomorphogenesis [15]. AtbZIP1 has been implicated in sugar signaling and growth regulation, indicating a connection between bZIP-mediated transcriptional regulation and carbohydrate status [24]. In addition, AtbZIP59 formed a complex with LATERAL ORGAN BOUNDARIES DOMAIN 16 (LBD16) to directly activate FAD-BINDING BERBERINE (FAD-BD) genes, thereby promoting callus initiation and cell fate transition [25]. FLOWERING LOCUS D (FD) and FLOWERING LOCUS T (FT) promoted floral transition and development via APETALA 1 (AP1) [16]. These findings indicate that bZIP transcription factors are involved not only in environmental signal responses but also in developmental reprogramming, organogenesis, and plant regeneration.

The bZIP transcription factors are also important components of hormone- and stress-responsive regulatory networks. In rice, OsbZIP23, OsbZIP66, and OsbZIP72 interacted with MOTHER OF FT AND TFL 2 (OsMFT2) and positively regulated the expression of ABA-responsive genes, thereby affecting seed germination [13]. In Fagopyrum tataricum, FtbZIP5 acted as positive regulators of drought resistance through ABA signaling, directly binding to SNF1-RELATED PROTEIN KINASES 2.6 (FtSnRK2.6) [26]. In non-heading Chinese cabbage, BcbZIP72 binded to the G-box element in the promoter of C-REPEAT BINDING FACTOR 1 (BcCBF1) and activated its transcription, thereby improving cold adaptability [19]. These studies suggested that bZIP proteins integrated environmental stimuli with endogenous hormonal pathways, particularly ABA-dependent signaling, to coordinate plant developmental and stress responses. In addition to their roles in growth regulation and stress adaptation, bZIP transcription factors are increasingly recognized as regulators of plant secondary metabolism. In Artemisia annua, AabZIP19 regulated the biosynthesis of phenolic acids and flavonoids by binding to the promoter of PHENYLALANINE AMMONIA-LYASE 1 (AaPAL1) [17]. In grapes, VvibZIPC22 boosted flavanol production by activating the transcription of flavonoid biosynthetic genes like CHALCONE SYNTHASE (VviCHS3), CHALCONE ISOMERASE (VviCHI), FLAVONOL SYNTHASE 1 (VviFLS1), and ANTHOCYANIDIN REDUCTASE (VviANR) [27].

Despite this broad functional relevance, understanding of bZIP genes in lily remains limited and is largely based on studies of individual genes involved in floral pigmentation, trichome development, and bulblet formation. For instance, in L. leichtlinii, LlbZIP11 was found to bind to the LlMYB19S promoter to activate its expression, and the two proteins cooperatively promote anthocyanin accumulation in the raised spots on tepals [28]. During trichome development in L. pumilum, LpNAC48, a natural 259-bp variation in the LpNAC48 with an ABA-responsive element (ABRE), which could be bound and activated by LpbZIP29 to promoting trichome formation [29]. Conversely, LlbZIP11 acted as a transcriptional repressor by binding to the ACGTT cis-regulatory element in the WUSCHEL-RELATED HOMEOBOX 11 (LlWOX11) promoter to suppress LlWOX11 expression and inhibit bulblet formation [30]. Although these studies demonstrate the functional importance of specific lily bZIP genes, they did not provide a systematic genome-wide view of bZIPs.

High-quality genome assembly in lily has long been impeded by its enormous genome size, high heterozygosity, and high proportion of repetitive sequences. The release of the L. davidii var. unicolor reference genome in 2024 has now enabled genome-wide investigations of genes associated with important agronomic traits [2]. Despite the established roles of bZIP transcription factors in plant growth and development, the composition and evolutionary history of the lily bZIP family remain poorly understood. It is also unclear which family members may participate in bulb dormancy release and asexual regeneration. This genomic resource now enables a systematic investigation of the bZIP family in lily, provided an unprecedented opportunity for systematic annotation and functional dissection of its transcription factor families. In the present study, we performed a genome-wide identification of LdbZIP genes in L. davidii var. unicolor and characterized their phylogenetic relationships, gene structures, conserved motifs, and chromosomal distribution. We further examined their expression dynamics during bulb dormancy and under different regeneration conditions. These findings establish a framework for dissecting the functional roles of LdbZIP transcription factors in bulb dormancy and asexual regeneration and provide candidate targets for molecular breeding in lily.

2 Materials and Methods

2.1 Plant Materials and Sample Collection

Three-year-old bulbs of the lily cultivar ‘Chengse Yangguang’ were obtained from the greenhouse at Shanghai Academy of Agricultural Sciences. Bulbs, with uniform maturity and without visible mechanical damage, pests, or disease symptoms were stored at 4°C for long-term cold treatment (LTCT). Healthy middle scales were collected after 0, 30, and 60 d (S1–S3), with three biological replicates for each stage. For in vitro bulblet induction, scales from healthy tissue-cultured plantlets of L. davidii var. unicolor were excised and cultured on water agar medium (6 g L−1) under a 16 h light/8 h dark photoperiod at a light intensity of 200 μmol m−2 s−1. The cultures were maintained for 3 weeks to induce bulblet formation at the basal ends of the scales. Basal scale tissue was sampled at 0, 5, 10, and 15 d. All samples were frozen immediately in liquid nitrogen and stored at −80°C.

2.2 Identification and Characterization of LdbZIPs

78 annotated Arabidopsis bZIP protein sequences were downloaded from the TAIR database (https://www.arabidopsis.org/) according to previous annotations and used as query sequences [31]. The genome assembly and annotation files of L. davidii var. unicolor, including genomic sequences, coding sequences, protein sequences, and GFF3 annotation files, were downloaded from the China National GeneBank Database (CNP0005511) [2]. To identify putative bZIP proteins in L. davidii var. unicolor, two complementary approaches were employed. A local protein database was constructed using TBtools-II (v2.363), and Blastp searches were performed against L. davidii var. unicolor protein sequences with an E-value threshold of <10−5 [32]. In parallel, the Pfam HMM profile of the bZIP domain (PF00170.27; http://pfam.xfam.org/) was searched against the same proteome using hmmsearch in HMMER v3.4. Model-specific trusted cutoff (TC) bit scores provided with PF00170.27 were applied using the -cut_tc option. The results from both approaches were merged, and redundant sequences were removed to generate a preliminary candidate list. Furthermore, to minimize the likelihood of overlooking bZIP genes, we conducted a supplementary search using tBlastn against the L. davidii var. unicolor genome, utilizing the 78 Arabidopsis bZIP protein and 89 rice bZIP protein, with an E-value threshold of <10−5 [20]. To verify the presence of the conserved bZIP domain, all candidate sequences were examined using the NCBI Batch CD-Search tool (https://www.ncbi.nlm.nih.gov/cdd/), and sequences lacking a complete bZIP domain were excluded. The physicochemical properties of the identified LdbZIP proteins, including amino acid length, molecular weight (MW), theoretical isoelectric point (pI), instability index (II), aliphatic index, and grand average of hydropathy (GRAVY), were computed using the ProtParam module in BioPython. Subcellular localizations were predicted utilizing the DeepLoc 2.0 server (https://services.healthtech.dtu.dk/services/DeepLoc-2.0/).

2.3 Phylogenetic Analysis and Classification of LdbZIPs

To elucidate the evolutionary relationships and classify the LdbZIP genes, a combined set of LdbZIP, AtbZIP, and OsbZIP full-length protein sequences was aligned with MAFFT v7 using the auto strategy [33]. No additional automated trimming was applied after alignment. The phylogenetic tree was subsequently constructed by IQ-TREE (v3.0.1) under the maximum likelihood (ML) framework. The optimal substitution model (LG+I+G) was automatically determined by ModelFinderPlus (MFP). Branch support was assessed using 1000 ultrafast bootstrap replicates [34]. To assess the potential effect of variable regions outside the bZIP domain on phylogenetic inference, an additional ML tree was constructed using only the conserved bZIP domain sequences. LdbZIP proteins were assigned to subfamilies according to the classification established for AtbZIPs [31].

2.4 Gene Structure and Conserved Motif Analysis of LdbZIPs

The exon-intron organizational features of LdbZIP genes were extracted from the L. davidii var. unicolor genome annotation file (gff3 file). Gene length, exon number, total exon length, intron length, and coding-sequence length were calculated for each LdbZIPs. Conserved motifs of the LdbZIP proteins were predicted using the MEME Suite (v5.5.7) under the following parameters: maximum number of motifs = 10, motif width = 6–50 amino acids, and site distribution set to ZOOPS (zero or one occurrence per sequence). The gene structures, conserved motifs, and the phylogenetic tree, were visualized collectively using TBtools-II and the ggtree and ggplot2 packages in R. The LdbZIPs subfamilies were delineated based on the phylogenetic clustering results.

2.5 Cis-Regulatory Element and Enrichment Analysis in LdbZIPs Promoters

To investigate potential regulatory mechanisms, the 4000 bp genomic sequences upstream of the start codon (ATG) of the LdbZIP genes were extracted as putative promoter regions using the GTF/GFF3 sequence extraction tool in TBtools-II. These sequences were submitted to the PlantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to identify cis-regulatory elements. The predicted functional elements were filtered, categorized, and visualized alongside the phylogenetic tree using the ggtree and ggplot2 packages in R. To perform a genome-background-based enrichment analysis of cis-regulatory elements, 4000 bp promoter regions of the 64 LdbZIP genes were used as the foreground set, whereas the corresponding promoter regions of all annotated non-LdbZIP genes located on the primary lily chromosomes served as the background set. Using the same cis-regulatory elements defined, each cis-regulatory element was scored for presence or absence at the promoter level. Enrichment significance was determined using Fisher’s exact test, followed by Benjamini-Hochberg false discovery rate (FDR) correction for multiple comparisons.

2.6 Chromosomal Distribution, Synteny Analysis and Gene Duplications of LdbZIPs

Physical chromosomal locations of the LdbZIP genes were retrieved from the L. davidii var. unicolor GFF3 file and visualized by ggplot2 in R. For gene duplication analysis, MCScanX was utilized to compute collinearity blocks and classify duplication events (singleton, dispersed, proximal, tandem, and segmental/WGD). These intra-species syntenic networks were visualized using the R package Circlize. To estimate the selective constraints on the duplicated genes, nonsynonymous (Ka) and synonymous (Ks) values were calculated using the ParaAT2.0.

2.7 Expression Pattern Analysis of LdbZIPs

Tissue-specific transcriptome data for L. davidii var. unicolor were obtained from the Liliales Genome Database (LGD), developed by the Lily Laboratory at Nanjing Agricultural University (https://lgd.njau.edu.cn/lily/home) [35]. These data were used to analyze the expression profiles of LdbZIP genes across different tissues. To examine the dynamic expression patterns of LdbZIP genes during dormancy release and regeneration, publicly available transcriptome datasets were downloaded from LGD and integrated with an in-house RNA-seq dataset. RNA-seq libraries of ‘Chengse Yangguang’ were constructed and sequenced using an Illumina high-throughput sequencing platform. In the present study, three developmental stages (S1–S3) were selected to analyze the expression patterns of LdbZIP genes during dormancy release. All values were transformed as log2(TPM + 1) or log2(FPKM + 1), and heatmaps were row-scaled by Z score. Two-group comparisons were evaluated with two-sided t test. Heatmaps of LdbZIP expression profiles were generated using TBtools-II and ClusterGVis in R [36].

2.8 RNA Extraction and qRT-PCR Analysis of LdbZIPs Expression Patterns

Total RNA from samples was extracted using Eastep® Super Total RNA Extraction Kit (LS1040, Promega). First-strand cDNA was then synthesized according to the manual of a PrimeScriptTM RT reagent Kit with gDNA Eraser (RR047, TaKaRa). The sequences of primers used in this study are listed in Supplementary Table S1. qRT-PCR was performed using TB Green® Premix Ex TaqTM (RR420A, TaKaRa) in a Roche Light Cycler 480 II System. The thermal cycling program was: 95°C for 30 s, followed by 40 cycles of 95°C for 5 s and 60°C for 30 s. A melting curve analysis was performed to confirm amplification specificity. Three biological replicates and three technical replicates were tested for each gene. F-BOX FAMILY PROTEIN (FP) and ELONGATION FACTOR 1-α (EF1-α) were used as the internal reference gene. Relative gene expression levels were calculated using the 2−ΔΔCt method.

3 Results

3.1 Identification and Characterization of LdbZIPs

Through comprehensive Blastp, HMMER homology searches and tBlastn against the L. davidii var. unicolor genome assembly (Supplementary Table S2–S4), we identified 64 non-redundant proteins harboring the conserved bZIP domain. These genes were systematically renamed from LdbZIP1 to LdbZIP64 according to their chromosomal locations, following the order of chromosome number and their physical positions from the top to the bottom of each chromosome. More detailed information of LdbZIPs was provided in Supplementary Table S4. To explore their structural and functional diversification, we conducted an in-depth statistical characterization of encoded proteins, quantifying six critical physicochemical parameters: protein length, MW, pI, GRAVY, II, and aliphatic index. Notably, the physicochemical properties of LdbZIPs displayed marked heterogeneity. Sequence lengths spanned 95 residues (LdbZIP63; MW = 11.12 kDa) to 628 residues (LdbZIP25; MW = 67.55 kDa). The pI values ranged from 4.49 (acidic LdbZIP16) to 11.39 (basic LdbZIP37), with 45.31% (29/64) of proteins classified as acidic properties (pI < 7) and the remaining as basic (pI > 7). Hydropathic analysis revealed predominantly hydrophilic characteristics, as evidenced by negative GRAVY values (−1.21 to 0.202), with a single exception (LdbZIP29; GRAVY = 0.202). Instability indices (27.73–78.40) suggested inherent structural lability, with 96.88% (62/64) of proteins exceeding the stability threshold (index > 40). Aliphatic indices demonstrated substantial variability (45.86–104.10), correlating with differential thermostability potentials. Subcellular localization prediction supported the putative transcriptional regulatory roles of LdbZIP proteins, as all members were predicted to localize to the nucleus. Notably, the predictions identified members with additional localizations, two paralogs (LdbZIP25 and 59) showed dual localization to the endoplasmic reticulum, while three members (LdbZIP32, 50, and 54) displayed additional cytoplasmic distribution patterns.

3.2 Phylogenetic Analysis and Classification of LdbZIPs

To elucidate the evolutionary relationships and subgroup classification of the bZIP family in L. davidii var. unicolor, we constructed a ML phylogenetic tree from L. davidii var. unicolor, A. thaliana, and Oryza sativa using full-length protein sequences and conserved bZIP domain sequences, respectively (Fig. 1 and Fig. S1). The two ML phylogenetic trees were largely congruent. Based on their clustering with AtbZIP homologs, the 64 LdbZIP proteins were assigned to 11 of the 13 recognized subfamilies (I-XIII). Subfamily VI formed the largest clade with 15 members, followed by subfamilies XIII (13), I (9), and V (9), whereas the remaining subfamilies contained between one and six members. No LdbZIP protein was assigned to subfamily IV or VIII in either the full-length or domain-based phylogenetic trees. Consistent with this result, the additional tBLASTn search did not identify any unannotated genomic locus encoding an intact bZIP protein assignable to either subfamily under the search criteria used (Supplementary Table S3). Thus, member of subfamilies IV and VIII were not detected in the current L. davidii var. unicolor genome assembly. This phylogenetic framework provides useful clues for predicting the potential biological roles of LdbZIPs based on characterized Arabidopsis homologs. Furthermore, several LdbZIP proteins clustered closely with OsbZIPs, indicating high sequence homology and potentially conserved functions among these orthologs.

Figure 1: Phylogenetic analysis of full-length bZIP proteins in L. davidii var. unicolor (Ld), A. thaliana (At), and O. sativa (Os). The red star, blue tick, and green triangle represent proteins from L. davidii var. unicolor, rice, and A. thaliana, respectively. The different colors in the outer circle indicate different bZIPs subfamilies. Node symbols indicate ultrafast bootstrap support.

3.3 Exon-Intron Structure and Conserved Motif Analyses of LdbZIPs

To further explore the structural characteristics and potential functional divergence of the LdbZIP gene family, we systematically analyzed their exon-intron structures and conserved motif compositions. Considerable variation was observed in the genomic lengths of LdbZIP genes. The majority of LdbZIPs were shorter than 100 kb, several members showed pronounced genomic length expansion. Specifically, LdbZIP41 exceeded 300 kb (352.6 kb), and four others (LdbZIP21, 29, 46, and 57) spanned over 100 kb (Fig. 2). The number of exons per gene ranged from 1 to 12, with an average of 4.92 exons. Notably, all nine members of subfamily I contained only a single exon and lacked introns, suggesting a highly conserved gene structure within this subfamily. Compared with bZIP genes reported in other plant species, the exon lengths of LdbZIP genes did not show obvious expansion (Supplementary Table S4). By contrast, their introns were markedly longer, which may be associated with the large genome size of the Lilium.

Figure 2: Phylogenetic relationships and exon-intron structure of LdbZIPs. The phylogenetic tree was generated from full-length sequences using ggtree by R. Different colored backgrounds indicate various subfamilies of LdbZIPs. The location of the exon and intron of LdbZIPs was visualized using ggplot2, with green boxes indicating exons and black lines indicating introns.

To identify conserved motifs in the LdbZIPs, we analyzed the protein sequences using MEME. A total of 10 conserved motifs (Motif 1–10) were detected among LdbZIP transcription factors (Fig. 3 and Supplementary Table S5). Notably, Motif 1 was present in almost all LdbZIP members, highlighting its high degree of conservation. In contrast, none of 10 motifs was detected in LdbZIP59, suggesting that it has undergone considerable sequence divergence during evolution. Furthermore, motif distribution was closely associated with phylogenetic classification. For instance, motif 2 and 9 were specific to subfamily XIII, especially motif 9. Within subfamily V, the most members possessed motif 3, 4, 6, and 7. However, LdbZIP46 lacked motif 3, 4, and 6, implying potential functional differentiation relative to other subfamily V proteins. In subfamily VI, motif 5 and 10 were present in most members, and motif 10 appeared to be unique to this subfamily. Overall, while motif compositions vary considerably among different subfamilies, members within the same subfamily exhibit highly similar motif patterns. This structural consistency suggests that LdbZIP proteins within the same subfamily likely share conserved biological functions.

Figure 3: Phylogenetic relationships and conserved motif distributions of LdbZIP proteins. The phylogenetic tree was generated from full-length sequences using ggtree by R. Different colors indicated various subfamilies of LdbZIPs. Conserved motifs in LdbZIP proteins were analyzed by MEME and were visualized by ggplot2. Motif 1–10 are represented by boxes in different colors, and their sequences are listed in Supplementary Table S5. The scale bar indicates protein length.

3.4 Cis-Regulatory Element Analysis and Functional Prediction of LdbZIPs

Analysis of cis-regulatory elements in the promoter regions of LdbZIPs may provide valuable clues to their potential biological functions and regulatory mechanisms. To characterize the regulatory landscape of this gene family, promoter sequences of 64 LdbZIPs were extracted and analyzed for cis-regulatory element composition. The results showed that these promoter regions were contained a wide variety of cis-regulatory elements (Fig. 4). In addition to the core promoter elements CAAT-box and TATA-box, which were essential for transcription initiation, as well as enhancer-related elements, we also identified various elements associated with light responsiveness (12 types), plant growth and development (9 types), plant hormone (10 types) and stress responses (13 types) (Supplementary Table S6). Stress-responsive elements were the most abundant in absolute frequency, dominated by MYB, MYC, and stress response elements (STRE). Notably, 11 LdbZIPs harbored 8 or more MYB elements, with LdbZIP47 containing the highest number (16), followed by LdbZIP18 (12). LdbZIP26 possessed the greatest number of MYC elements (8), whereas LdbZIP42 contained the highest number of STRE elements (11). Plant-hormone-responsive elements were also widely distributed. These included ABRE, AuxRR-core, TGA-element, CGTCA-motif, TGACG-motif, ERE, P-box, TATC-box, and TCA-element. Among these, LdbZIP4 contained the largest number of ABRE motifs (17), followed by LdbZIP22 (11), suggesting that these genes may be particularly responsive to ABA-related regulation. Many of these phytohormones were closely linked to stress signaling, the enrichment of hormone-responsive elements and stress responsive elements further supported the potential involvement of LdbZIPs in stress responsiveness. Although LdbZIP genes promoter regions contain large numbers of cis-regulatory elements, enrichment analysis comparing the presence or absence of elements between the LdbZIPs set and the genome-wide promoter background revealed that these cis-regulatory elements were not significantly enriched in LdbZIPs promoters (Supplementary Table S6). Overall, the PlantCARE annotations provide testable candidates for subsequent studies.

Figure 4: Phylogenetic relationships and cis-regulatory elements analysis of LdbZIPs promoters. The phylogenetic tree was generated from full-length sequences using ggtree by R. Different colors indicate various LdbZIPs subfamilies. The cis-regulatory elements in the LdbZIPs promoters were classified into four functional categories for statistical analysis. The proportion of each category of cis-regulatory elements is illustrated in a stacked bar chart. Green, blue, dark gray, and red represent light responsiveness, plant growth and development, plant hormone-related, and stress-related elements, respectively. The heatmap summarize the number of different cis-regulatory elements within each functional category of LdbZIPs promoters, with a color gradient from white to red indicating an increase in the number of elements.

3.5 Chromosomal Distribution, Synteny Analysis and Gene Duplications of LdbZIPs

To explore the chromosomal distribution of the LdbZIP gene family, we anchored all 64 identified members onto the 12 chromosomes of L. davidii var. unicolor using the annotated genome assembly (Fig. 5). The genes were clearly unevenly distributed. Chr 3, 5, and 9 carried the highest numbers, with 11, 10, and 10 members respectively, while Chr 7 contained only a single gene (LdbZIP37). Notably, gene density showed no obvious correlation with chromosome length. Genome-wide synteny analysis uncovered an intricate network of homologous relationships. Within this network, 55 duplicated LdbZIP gene pairs were identified, including 54 segmentally duplicated pairs and only one tandemly duplicated pair (Supplementary Table S7). The predominance of segmental duplication indicates that it was the principal duplication mode associated with the retention and expansion of the LdbZIP family. All duplicated gene pairs had Ka/Ks ratios below 1, indicating that their protein-coding sequences evolved predominantly under purifying selection to preserving the structural and functional integrity. Together, these results suggest that segmental duplication contributed substantially to LdbZIP family evolution.

Figure 5: Chromosomal distribution and intraspecies synteny analysis of LdbZIPs. From the outermost to the innermost layer, the figure shows gene names and locations on the chromosome, followed by the chromosomal scaffold, gene density on the chromosomes, the GC content, and collinearity relationships, respectively. Gray lines represent all syntenic blocks within the genome, and red lines correspond to syntenic LdbZIP gene pairs.

3.6 Expression Patterns of LdbZIPs in Different Tissues

To explore the potential functional diversification of the LdbZIP gene family, we analyzed their expression profiles across 13 tissues of L. davidii var. unicolor using transcriptomic data from the LGD. The tissues examined included bulb, bulb root, stem root, ovary, main stem, lateral stem, anther, filament, inner petal, outer petal, style, leaf, and bulblet. The analysis revealed highly divergent expression patterns (Fig. 6). 10 LdbZIPs (LdbZIP29, 46, 47, 48, 50, 54, 55, 57, 59, and 63) were undetectable across all tissues, implying that their transcription may be induced by specific environmental conditions or restricted to unexamined developmental stages. In contrast, several LdbZIPs, including LdbZIP9, 10, 17, 20, and 64, showed relatively high transcript abundance across most tissues, implying possible roles in basic developmental and physiological processes.

A subset of LdbZIPs displayed clear tissue specificity, particularly in reproductive organs. Notably, LdbZIP49 and 56 were exclusively expressed in the ovary and were undetectable in the other tissues, suggesting potential roles in ovary development or function. Likewise, LdbZIP30, 39, 41, and 58 were markedly upregulated in anthers, with LdbZIP30 exhibiting especially strong anther-specific expression. In addition, several genes showed differential expression between inner and outer petals. For example, LdbZIP22 and 37 were more highly expressed in inner petals, whereas LdbZIP13 was preferentially expressed in outer petals, implying possible involvement in petal differentiation and morphogenesis. Collectively, these patterns suggest that distinct LdbZIP members participated in the development and functional regulation of different floral organs.

In addition to floral tissues, several LdbZIPs were localized accumulation in vegetative tissues. LdbZIP5, 15, 43, and 61 showed their highest expression levels in roots (bulb roots and stem roots), indicating potential roles in root development or environmental adaptation. Meanwhile, LdbZIP2, 7, 23, and 25 were highly enriched in bulbs, suggesting critical roles in storage organ formation or nutrient accumulation.

Overall, the LdbZIP gene family exhibited extensive tissue-specific expression diversity in L. davidii var. unicolor. These findings support the view that LdbZIP members had undergone substantial functional diversification and played distinct roles in plant growth, reproductive development, and the formation of underground storage organs. This expression atlas provides a valuable basis for future functional studies of candidate LdbZIPs involved in vegetative growth, floral organ development, and organ-specific regulatory processes in L. davidii var. unicolor.

Figure 6: Tissue-specific expression profiles of LdbZIP genes in L. davidii var. unicolor. Heatmap illustrates the relative transcript abundance of LdbZIP genes in different tissues. The horizontal axis shows different tissues (including bulb, bulb root, stem, leaf, etc.). The vertical axis lists individual LdbZIPs. The color scale (ranging from blue to red) represents the normalized log₂(TPM + 1) values, with red indicating higher expression levels. For heatmap visualization, expression values were scaled by row using Z-score transformation.

3.7 Expression Patterns of LdbZIPs during LTCT-Induced Dormancy Release

To further clarify the potential involvement of the LdbZIP gene family in LTCT-induced dormancy release in lily, we systematically examined two transcriptome datasets derived from different lily cultivars and tissues.

Previous work has shown that LTCT at 4°C effectively promotes sprouting, bud development after bulb dormancy release, and floral transition (vernalization) in the lily cultivar ‘Siberia’ [6]. To examine the potential involvement of LdbZIP genes in these processes, we analyzed transcriptional dynamics from the shoot apical meristem (SAM) of dormant bulbs (DB) and growth-phase-transited bulbs (GTB) treated by LTCT. In total, 48 LdbZIP genes were detected. Between DB and GTB, 26 genes (54.17%) showed decreased transcript abundance, whereas remaining 22 showed increased abundance (Fig. 7a). Among them, 24 genes were significant changed (p < 0.05) in transcript abundance between DB and GTB, with 12 showing highly significant changes (p < 0.01; Supplementary Table S8). Notably, these contrasting responses indicate substantial transcriptional reprogramming of the LdbZIP genes during the transition from dormancy to growth.

Because lily bulbs depend heavily on nutrient reserves stored in scales during dormancy and subsequent reactivation, we further generated a time-course transcriptome dataset using bulb scales of the newly bred cultivar ‘Chengse Yangguang’ during LTCT at 4°C. 37 LdbZIP genes were expressed in during scales during dormancy release in ‘Chengse Yangguang’. The expression profiles of these LdbZIPs were classified into eight clusters, which could be further grouped into four major patterns (Fig. 7b). The first pattern comprised genes showing continuous upregulation during LTCT (Clusters 5, 6 and 8). The second included genes that were initially induced and subsequently repressed (Clusters 3 and 4). The third consisted of genes showing continuous downregulation during cold storage (Clusters 2). The fourth included genes that were initially downregulated and then upregulated (Clusters 1 and 7). Notably, 11 genes showed significant differential expression in both S2 vs. S1 and S3 vs. S2 (Supplementary Table S9). Consequently, the expression patterns of 7 selected genes in scales were further examined during LTCT by qRT-PCR (Fig. 7c–i). LdbZIP25 and LdbZIP11 were induced during LTCT, whereas the expression levels of LdbZIP42, 41, and 17 declined. LdbZIP62 exhibited a transient, sharp increase at S2, followed by a significant decrease at S3. Although LdbZIP34 showed a brief rise at S2, its expression level at S3 remained lower than at the start of the treatment. All 7 candidate LdbZIPs responsed significantly to low temperature.

Taken together, these results demonstrate that the LdbZIPs are broadly responsive to LTCT and is likely involved in the regulatory network underlying bulb dormancy release and subsequent growth recovery in lily. This expression profiling analysis provides important evidence for the involvement of LdbZIP genes in LTCT-induced dormancy release in lily and lays a foundation for the identification of key candidate genes for future functional studies.

Figure 7: Expression profiles of LdbZIP genes during LTCT-induced dormancy release in lily. (a) Expression heatmap of LdbZIPs in the SAM from DB (LTCT for 0 weeks) and GTB (LTCT for 8 weeks) of ‘Siberia’. (b) Temporal expression profile of LdbZIPs in scales of ‘Chengse Yangguang’ at three developmental stages, which LTCT for 0 d (S1), 30 d (S2), and 60 d (S3). Genes are grouped into eight clusters based on expression profiles, with line graphs showing the consensus trajectory and gene count per cluster. (c–i) Relative expression levels of 7 candidate LdbZIP genes in scales during LTCT by qRT-PCR. In (a,b), color scales represent row-scaled Z-scores derived from log₂(TPM + 1) and log₂(FPKM + 1), respectively, ranging from blue (down-regulated) to red (up-regulated). Heatmap data are derived from three biological replicates per stage. Error bars in (c–i) represent mean ± SEM of three technical replicates. Asterisks indicate statistically significant differences by one-way ANOVA with Tukey’s multiple comparisons test (*p < 0.05, **p < 0.01; ***p < 0.001, ****p < 0.0001, ns: no significant difference (p > 0.05)).

3.8 Expression Patterns of LdbZIPs during Bulbil Formation

To elucidate the potential regulatory roles of the LdbZIPs in asexual regeneration, we first investigated the expression profiles of LdbZIPs during bulbil initiation and development. Previous studies have shown that bulbils formed spontaneously in the leaf axils of lily [37]. We therefore retrieved LdbZIP transcript abundances of four distinct developmental stages, Ll_DN, Ll_UN, Ll_UT and Ll_UE. 59 members were detectable during this process, and 43 of them differed significantly in at least one pairwise comparison (p < 0.05), indicating that bulbil formation was accompanied by reprogramming of LdbZIPs (Fig. 8a and Supplementary Table S10). A subset of genes, including LdbZIP1, 6, 27, 34, 58, and 60, showed progressive transcript accumulation as axillary buds developed from stage Ll_DN to Ll_UE. Notably, LdbZIP60 expression was near the limit of detection at Ll_DN (mean TPM = 0.02) but increased to 27.27 at Ll_UE, indicating strong transcriptional activation during primordium outgrowth. In contrast, genes such as LdbZIP8, 11, 14, 22, 28, and 45 exhibited highest transcript levels in Ll_DN and underwent down-regulation during bulbil differentiation. A third group, represented by LdbZIP9 and 10, displayed transient expression peaks specifically at the Ll_UT stage, coinciding with the morphological initiation of bulbil.

Phytohormones, particularly auxin, have been implicated in lily bulbil formation, and inhibition of polar auxin transport by N-1-naphthylphthalamic acid (NPA) promoted this process [37,38]. We therefore compared LdbZIPs expression in leaf axils treated with IAA or NPA with that in the corresponding control (Fig. 8b and Supplementary Table S11). 5 LdbZIPs (LdbZIP2, 18, 31, 34 and 42) differed significantly between the IAA and control. Among these responses, the decrease in LdbZIP42 was the largest. NPA treatment affected a broader set of genes. 11 LdbZIPs differed significantly from the control, of which two (LdbZIP15 and 17) were induced and nine (LdbZIP12, 20, 22, 23, 28, 31, 34, 42, and 58) were repressed.

Combining the transcriptome data, 10 candidated LdbZIPs were selected for further validation during bulbil formation by qRT-PCR. Scales of L. davidii var. unicolor were cultured on water agar medium to induce the bulbil in vitro. During the process, 5 LdbZIPs accumulated markedly. LdbZIP60, 58, and 28 shared a similar temporal expression profile, with expression level increasing progressively during bulblet emergence and reaching its highest after bulblet formation (Fig. 8c–e). LdbZIP2 and LdbZIP12 exhibited a more dynamic, fluctuating transcriptional profile throughout this developmental process (Fig. 8f,g). While LdbZIP42 and 34 maintained robust expression during the initial induction phase, their transcription was significantly attenuated after the bulblets were formed (Fig. 8h,i). LdbZIP17, 45, and 22 were repressed the activity during the organogenesis process, especially LdbZIP22, which expression was almost detectable during the later developmental stages (Fig. 8j–l). These results provide a basis for further functional validation of key LdbZIPs in the regulatory network controlling bulbil formation in lily.

Figure 8: Expression dynamics of LdbZIP genes involved in bulbil formation of lily. (a) Heatmap illustrating the expression profiles of LdbZIP genes at four stages of axillary bulbil formation. Samples were collected from four development stages, which were no bulbils at 10 cm down from the bottom (Ll_DN), no upper bulbils (Ll_UN), transparent bulges in upper leaf axils (Ll_UT), and white bulges emerging in upper leaf axils (Ll_UE), with three biological replicates per stage. (b) Heatmap displaying the transcriptional responses of LdbZIP genes in the leaf axils under IAA and NPA treatment, with three biological replicates per stage. The leaf axils treated with distilled water served as the control, with two biological replicates. Color scales in (a,b) represent row-scaled Z-scores derived from log₂(TPM + 1), ranging from blue (down-regulated) to red (up-regulated). (c–l) Relative expression levels of 10 candidated LdbZIPs during in vitro bulbil formation from scales cultured on water-agar medium. Samples were collected at 0, 5, 10, and 15 d after culture. Expression levels were determined by qRT-PCR. Error bars indicate mean of ± SEM of three technical replicates. Asterisks indicate statistically significant differences by one-way ANOVA with Tukey’s multiple comparisons test (*p < 0.05, **p < 0.01; ***p < 0.001, ****p < 0.0001, ns: no significant difference (p > 0.05)).

3.9 Expression Patterns of LdbZIPs during Callus Induction

To elucidate the potential function of LdbZIP genes during cell dedifferentiation and callus formation, we analyzed time-course transcriptome data from root explants treated with picloram (PIC) [39]. The LdbZIPs underwent marked transcriptional reprogramming during this developmental transition. Based on their expression dynamics, these genes were assigned to five distinct co-expression modules (Clusters 1–5) (Fig. 9).

A substantial portion of the genes, comprising Cluster 2 (13) and Cluster 4 (13) alongside members of Cluster 5 (10), exhibited immediate downregulation by treated with PIC. Specifically, transcripts in Clusters 4 and 5 peaked strictly in CK and plummeted significantly upon entry into ECI1. Several members of Cluster 4 (such as LdbZIP21, 13, and 60), displayed high expression in CK but were significantly suppressed following PIC treatment (Supplementary Table S12). This rapid transcriptional repression suggests that these LdbZIPs may participate to maintain root identity. Conversely, genes grouped in Cluster 1 (13) and Cluster 3 (8) displayed strong auxin responsiveness, characterized by significant upregulation post-induction. LdbZIPs in Cluster 1 (such as LdbZIP62 and 53) transcripts accumulated rapidly during ECI1 and remained abundant through ECI2. As the culture progressed to ECI3, extensive embryogenic callus formation was observed morphologically. Concurrently, genes within Clusters 2 and 3 underwent a pronounced late-stage transcriptional burst, such as LdbZIP26 and 43.

Collectively, these distinct expression patterns revealed stage-dependent and gene-specific regulation of the LdbZIP family during embryogenic callus induction. Different family members are sequentially recruited or repressed to orchestrate the complex progression of lily cell dedifferentiation. Given their tight correlation with critical developmental milestones, candidate genes from Clusters 1–3, including LdbZIP62, 53, 23, 33 and 11, hold potential for future functional validation regarding plant regeneration capacity.

Figure 9: Expression dynamics of LdbZIP genes during embryogenic callus (EC) induction. Heatmap and temporal expression trajectories of LdbZIP genes across four distinct dedifferentiation stages. Roots, which treated by 1 mg L−1 PIC, were sampled at 0 d (CK), 20 d (ECI1, initiation of epidermal cell division), 40 d (ECI2, emergence of small ECs), and 55 d (ECI3, massive proliferation of ECs). Each stage contains two biological replicates. Genes were grouped into five clusters based on their expression patterns. The line charts on the right illustrate the scaled expression trends for each cluster, with the number of genes indicated. Color gradients in the heatmap represent row-scaled Z-scores.

4 Discussion

The bZIP transcription factors constitute one of the most important regulatory protein families in plants, participating in diverse biological processes, including growth and development, abiotic stress responses, and secondary metabolite biosynthesis [40]. Genome-wide identification of bZIP genes has been reported in many plant species beyond model species, such as papaya with 52 [21], Chinese cabbage with 150 [19], and E. salsugineum with 85 [23]. However, the extremely large and highly repetitive genome of lily has long limited systematic genomic studies in this species. The recent release of the first L. davidii var. unicolor genome therefore provides an important opportunity to investigate the bZIP gene family at the whole-genome level.

In this study, 64 LdbZIP genes were identified in the L. davidii var. unicolor genome, and no member of subfamily IV or VIII was detected (Fig. 1). Similarly, the absence of bZIP genes in a specific subfamily has also been reported in flax and papaya [11,21]. Segmental duplication was the predominant mechanisms underlying the expansion of the LdbZIP gene family (Fig. 5 and Supplementary Table S7). This pattern was consistent with that reported in flax and Liriodendron chinense, highlighting the central role of large-scale duplication in the evolution of plant bZIP gene families [11,41]. Lilies possess large and complex genomes whose evolutionary history includes two rounds of Whole Genome Duplication (WGD)-an ancient event shared with other monocot lineages and a subsequent duplication associated with the evolution of Lilium [1,2,42]. Together with extensive transposable element (TE) proliferation, these events have contributed to pronounced genome expansion and structural dynamism. Despite this highly dynamic genomic background, all duplicated LdbZIP gene pairs had Ka/Ks ratios below 1, indicating that they have been subject to purifying selection. Such selective constraints likely facilitated the removal of deleterious nonsynonymous mutations after duplication, thereby preserving the structural integrity and essential regulatory functions of LdbZIP proteins. Purifying selection may therefore have played an important role in maintaining the functional conservation and evolutionary stability of this transcription factor family despite extensive genome duplication, TE accumulation, and genomic reorganization in lilies.

The physicochemical properties of LdbZIPs were highly diverse, reflecting broad functional specialization. Subcellular localization prediction placed all LdbZIP proteins in the nucleus, consistent with their roles as transcription factors and with previous observations in other plant species [21,43]. Interestingly, several members, including LdbZIP25, 32, 50, 54, and 59, were also predicted to localize to the endoplasmic reticulum (ER) or cytoplasm (Supplementary Table S4), suggesting potential dual-localization patterns. Such non-canonical localization may reflect functional specialization. Particularly, nuclear-ER dual localization is closely associated with the unfolded protein response (UPR), a conserved signaling pathway activated by ER stress. When adverse conditions triggered ER stress, bZIPs were activated and mobilized to the nucleus to restore cellular homeostasis [44]. Similar patterns have emerged in pak choi, that 33 BcbZIP proteins were predicted to localize to both nucleus and ER, with co-expression and GO analyses supporting roles in UPR [45]. In rice, the RICE SEED bZIP1 (RISBZ1)/bZIP58 protein interacts with bZIP50 and 60 to antagonistically suppress downstream UPR genes, while simultaneously modulating seed storage protein and starch biosynthesis [46]. In Arabidopsis, AtbZIP17, 28, and 60 operated as core UPR regulators [47,48]. Recent work has further linked UPR-associated bZIPs to thermomorphogenesis, that bZIP17, 28, and 60 cooperated with PHYTOCHROME-INTERACTING FACTOR 4 (PIF4) to promote hypocotyl elongation under warm temperatures [49]. Within the LdbZIP genes, LdbZIP25 and 59, which clustered into the same subfamily with AtbZIP28, exhibited the nuclear-ER dual localization pattern. These observations identify them as candidate components of ER stress-associated signaling. Although heat and biotic stress-induced ER stress and UPR activation are well characterized, the involvement of the classical UPR-bZIP signaling cascade in cold stress response and plant regeneration remains elusive. Notably, LdbZIP32 and 50, which both exhibited nuclear-cytoplasmic dual localization, clustered with the same subfamily as VirE2-INTERACTING PROTEIN 1 (AtVIP1), raising the possibility of conserved functions. In Arabidopsis, VIP1 translocated from the cytoplasm to the nucleus upon flg22 or Agrobacterium treatment and regulates defense-related genes such as PATHOGENESIS-RELATED 1 (PR1) [50]. In tobacco, REPRESSION OF SHOOT GROWTH (RSG), the ortholog of VIP1, and its phosphorylation-related mutants displayed nuclear-cytoplasmic shuttling after Agrobacterium infiltration, with MITOGEN ACTIVATED PROTEIN KINASES (MAPK)-associated phosphorylation sites controlling this process [51]. Together, these dual localizations highlight the sophisticated post-translational mechanisms, which LdbZIPs likely employ to integrate environmental stress signals.

Cis-regulatory elements in promoter play essential roles in transcriptional regulation by mediating interactions between transcription factors and their target genes [52]. The promoters of LdbZIP genes contained abundant cis-regulatory elements associated with responsive to light, plant growth and development, plant hormone and stress (Fig. 4), implying broad involvement of these genes in growth regulation and environmental adaptation. The bZIP domain comprises a basic DNA-binding region and a leucine zipper that mediates homo- or heterodimerization [53]. This architecture enables bZIP dimers to recognize palindromic or pseudopalindromic sequences containing an ACGT core, including the G-box (CACGTG), C-box (GACGTC), and A-box (TACGTA) [54]. However, plant bZIP proteins do not exhibit a uniform affinity hierarchy for these motifs. Binding specificity therefore reflects the combined effects of the basic and hinge regions, nucleotides flanking the ACGT core, and dimer composition. DAP-seq and ChIP-seq analyses have further shown that sequence recognition is modified by genomic and cellular contexts. Using double DAP-seq, it demonstrated that heterodimerization substantially altered the binding profiles of Arabidopsis C/S1 bZIP proteins. The canonical G-box was more frequent in bZIP53 homodimer-specific regions, whereas TGAC half-sites and GCN4-like TGACTCA motifs were commonly associated with heterodimer-specific binding. These altered preferences were linked to ABA-responsive gene regulation by bZIP9-S1 heterodimers and seed-maturation genes targeted by bZIP53-group C heterodimers [55]. Thus, heterodimerization can expand the regulatory repertoire of bZIP proteins by redirecting them toward distinct sequence classes. Accordingly, the presence of ACGT-containing motifs in promoters indicates potential bZIP regulatory sites but does not alone establish direct binding or motif preference.

Expression profiling revealed distinct spatiotemporal specificities among LdbZIPs. 10 of the 64 LdbZIP genes initially lacked detectable expression in 13 tested tissues (Fig. 6). Comparable patterns have been reported for other gene families in L. davidii var. unicolor. 29 of 41 HEAT SHOCK TRANSCRIPTION FACTOR (LdHSF) genes and 23 of 115 LdWRKY genes were undetected or expressed at very low levels in the surveyed tissues [56,57]. However, under specific inductive conditions, several of LdbZIP genes (such as LdbZIP29, 50, 54, and 55) were detected expression, though at low levels (Fig. 7, Fig. 8 and Fig. 9). The low or condition-dependent expression of certain LdbZIP members may be intrinsically linked to the unusual architecture of the Lilium genome. Recent genome assemblies of L. davidii, L. sargentiae, and L. regale consistently revealed massive genome expansion driven by TE proliferation, particularly long terminal repeat (LTR) retrotransposons [1,2,42]. In lilies, TEs frequently accumulated within intronic and intergenic regions, leading to ultra-long genes. While extreme gene length alone could constrain transcriptional elongation, extensive TE insertions also introduce complex epigenetic regulation [2,42]. Epigenomic data from L. regale further showed associations of LTR retrotransposons with DNA methylation and H3K9me2, as well as a relationship between intronic TE density and alternative splicing [42]. These associations suggest that local TE and chromatin features could affect gene expression or transcript processing.

Analysis of the LTCT datasets indicated that LdbZIPs orchestrate dormancy release through tissue- and stage-specific modules, rather than a uniform cold-response program. Although their promoters harbor hormone- and stress-responsive motifs, the lack of significant enrichment relative to the genomic background suggests these elements represent candidate regulatory inputs rather than definitive evidence of direct hormonal control (Fig. 4 and Supplementary Table S6). This interpretation aligns with the established paradigm that bud dormancy was coordinated by a complex interplay of endogenous hormones and environmental cues [58]. Furthermore, the divergent expression patterns observed between SAM of ‘Siberia’ and the scales of ‘Chengse Yangguang’ were biologically highly plausible. Meristematic cells must resume division and drive organ growth, whereas scales primarily function as nutrient sinks with distinct hormonal and metabolic states.

Within scales, the sustained induction of LdbZIP25 and LdbZIP11 was consistent with their possible involvement in reserve mobilization or growth resumption. In contrast, the progressive decline in LdbZIP42, LdbZIP41, and LdbZIP17 transcripts may reflect attenuation of a dormancy-associated transcriptional state. The transient expression of LdbZIP62 could mark the transition from cold perception to developmental reactivation. These interpretations remain tentative, as transcript abundance alone cannot determine whether a gene promotes dormancy release, restrains it, or responds secondarily to the transition. Subsequent evolutionary analysis revealed that LdbZIP62 was a homolog of ABA INSENSITIVE 5 (AtABI5). ABI5, a well-characterized bZIP transcription factor, plays a central role in ABA-mediated seed dormancy and germination. ABI5 directly bound to the promoters of ABA receptor genes PYRABACTIN RESISTANCE 1-LIKE (PYL11/12) [59], while MYB30-INTERACTING E3 LIGASE 1 (MIEL1) mediated proteasomal degradation of ABI5 to fine-tune ABA signaling [60]. In wheat, SIMILAR TO RCD 1 (TaSRO1) also recruited HISTONE DEACETYLASE INTERACTING WITH SRO1 (TaHIS1) to induce histone deacetylation at ABI5 target promoters, thereby epigenetically repressing their expression and facilitating the transition from dormancy to germination [61]. In addition, in rice, OsbZIP58 had been reported to connect temperature-dependent sugar depletion with ABA-responsive transcription, starch mobilization, and premature dormancy release [62]. LdbZIP42 and LdbZIP17 were assigned to the same subfamily as OsbZIP58, raising the possibility that these lily genes operate at an interface between ABA signaling and carbon use. Based on the combined phylogenetic and expression evidence, we hypothesize that LTCT recruits distinct LdbZIP modules in meristematic and storage tissues.

Plant regeneration, underpinned by cellular totipotency, involves dedifferentiation during callus induction and subsequent redifferentiation through organogenesis or somatic embryogenesis. This process forms the foundation of genetic improvement in many crops. Increasing evidence indicate that bZIP transcription factors play important roles in regeneration by modulating hormone signaling and cell fate transitions. In this study, LdbZIP genes exhibited pronounced transcriptional dynamics during bulbil formation and embryogenic callus induction (Fig. 8 and Fig. 9), suggesting their participation in these developmental transitions. During the bulblet formation in L. davidii var. unicolor, the expression of LdbZIP58 and 60, belonging to the same subfamily, was significantly upregulated during the induction process and sharing identical expression patterns, they likely be associated with bulblet formation. LdbZIP34 was upregulated during the early stages of induction, suggesting a potential positive role in the initiation of bulblet formation. In contrast, LdbZIP17, 45, and 22 may act as negative regulators of bulblet formation. Consistent with this finding, several bZIPs have been shown to regulate regeneration across species. In carrot, CAREB1/2, the bZIP transcription factors, showed distinct expression patterns in response to high sucrose and ABA treatments. CAREB2 was rapidly induced, whereas CAREB1 was highly expressed at late developmental stages during somatic embryogenesis. Overexpression of CAREB1 arrests embryo growth at the torpedo stage [14]. In Arabidopsis, AtbZIP59 formed a heterodimeric complex with LBD transcription factors to co-regulate auxin-induced callus formation. Notably, overexpression of AtbZIP59 induced spontaneous callus formation even in the absence of exogenous auxin, and AtbZIP59 and AtLBD16 shared common downstream target genes, establishing the bZIP59-LBD complex as a key regulator of somatic cell reprogramming [25]. In Tartary buckwheat, time-series transcriptome analysis of hypocotyl-derived callus identified three FtbZIP genes strongly associated with callus induction, one of which showed the highest connectivity within its co-expression module [63]. In lily, the WRKY28-ABI5 module promotes embryogenic callus formation by elevating endogenous ABA levels following PIC treatment, providing direct evidence that bZIP factors participate in lily regeneration [64]. Taken together, the dynamic expression profiles of LdbZIP genes observed here support their involvement in hormone-mediated cell fate reprogramming. These candidates represent promising targets for improving regeneration efficiency and establishing more reliable genetic transformation systems in lily.

The findings of this study are based primarily on descriptive and correlation analyses. Although expression profiling provides valuable clues for inferring the potential functions of LdbZIP genes, expression data alone are insufficient to establish their specific biological roles. By examining the expression patterns of LdbZIP genes during lily dormancy and regeneration, we identified several candidate genes and preliminarily confirmed their expression trends using qRT-PCR. These results allowed us to propose hypotheses regarding their potential involvement in the regulation of dormancy and regeneration. Future studies will employ gain- and loss-of-function approaches, including overexpression and genome editing, to genetically validate these candidate genes and further clarify their biological functions and underlying regulatory mechanisms during lily dormancy and regeneration.

5 Conclusion

In conclusion, this study provided the first genome-wide characterization of the bZIP gene family in L. davidii var. unicolor. A total of 64 LdbZIP genes were identified and classified into 11 subfamilies. The LdbZIP proteins showed considerable variation in physicochemical properties and predicted subcellular localization. A wide variety of light-, hormone-, and stress-responsive cis-regulatory elements in their promoters suggests that this gene family may participate extensively in developmental regulation and environmental adaptation. Expression analyses further demonstrated that LdbZIP genes displayed marked tissue specificity and dynamic transcriptional changes during key biological processes. Particularly, the pronounced expression reprogramming observed during dormancy release and regeneration indicates that specific LdbZIP members may play important regulatory roles in bulb developmental transitions and asexual propagation. These findings expand our understanding of the evolutionary organization and functional potential of the bZIP gene family in lily and provide a valuable set of candidate genes for future studies. This study is limited, and the results presented here are mainly descriptive and correlational. Although further functional validation is required, this study establishes a valuable foundation for future studies on the regulatory mechanisms of LdbZIPs underlying dormancy release, regeneration capacity, and adaptive growth in L. davidii var. unicolor.