Comparative transcriptomics identifies LlCSLE1 as a key regulator of cellulose accumulation during bulbil formation in lily

Lily (Lilium spp.) is a perennial bulbous flower belonging to the family Liliaceae. It holds considerable importance in the ornamental horticulture industry and agricultural economy because of its multifunctional value in decoration, medicinal applications, and edible uses^[1,2]. Commercial lily propagation primarily relies on scale cuttings. In recent years, bulbil propagation has attracted increasing attention compared to seed-based reproduction, owing to its higher efficiency, shorter growth cycles, and superior preservation of maternal genetic traits^[3−5].

Bulbils, as specialized vegetative reproductive organs, occur in only a limited number of plant species, including Dioscorea batatas, Allium sativum, Titanotrichum oldhamii, Pinellia ternata, Agave tequilana, and selected Lilium cultivars^[6]. The molecular mechanisms underlying bulbil formation exhibit significant parallels with general organogenesis processes, as they are primarily derived from the proliferation and differentiation of axillary meristem cells^[7,8]. The bulbils of Lilium species typically develop in leaf axils, exhibiting spherical or ellipsoidal morphology with smooth surfaces. Their coloration varies according to developmental stages, predominantly appearing green or reddish-brown. Under favorable conditions, these bulbils detach from the mother plant and subsequently undergo root initiation and shoot development to establish new plants^[9]. Certain species within the genus Lilium (e.g., Lilium lancifolium, Lilium sulphureum, and Lilium sargentiae) exhibit a stable bulbil-forming capacity, while others completely lack this trait. Bulbil propagation is recognized as a novel asexual reproduction method with substantial reproductive potential in lilies. However, under natural conditions, the limited bulbil production cannot meet the demands of large-scale lily cultivation. Bulbil formation is co-regulated by phytohormones, environmental factors, and genetic determinants^[10]. Recent studies have demonstrated that increased auxin concentration in leaf axils reduces bulbil formation^[11], while cytokinins play a positive regulatory role^[10]. Furthermore, carbohydrate metabolism, particularly sucrose and starch metabolism, is essential in bulbil initiation and development^[12]. In Lilium lancifolium, both sucrose synthase (LlSusy1) and cell wall invertase (LlCWIN2) have been demonstrated to promote bulbil initiation, while auxin exerts regulatory control over bulbil development through the modulation of carbohydrate metabolism-related genes^[11]. In Aladdin lily, 6-BA treatment during bulbil formation induces significant alterations in starch and sucrose metabolism, with marked upregulation of SUS1 and TPS6 expression. This suggests that sucrose may serve as a critical signaling molecule for bulbil sprouting^[13]. Similarly, sucrose signaling has been demonstrated to promote bulbil growth during bulbil sprouting in yam (Dioscorea spp.), revealing a conserved regulatory mechanism across monocot and dicot species^[14].

The plant cell wall is a polysaccharide-rich extracellular matrix primarily composed of cellulose, hemicellulose, and pectin^[15,16]. This sophisticated structure provides mechanical support and cellular protection and actively participates in numerous critical physiological processes. Substantial evidence demonstrates its essential roles in cell fate determination^[17], organogenesis^[18,19], and plant-environment interactions^[20]. Carbohydrate metabolism provides carbon skeletons and energy for the synthesis of cell wall polysaccharides, and newly fixed photosynthetic carbon is channeled into cell wall biosynthetic pathways mainly in the form of nucleotide sugars such as UDP-glucose^[16,21]. Alterations in carbohydrate metabolic status directly affect the partitioning of carbon flux toward the cell wall, thereby regulating the rate and total amount of cell wall synthesis^[21]. In Arabidopsis, a regulatory network that couples carbohydrate metabolism and cell wall remodeling through signaling nodes such as T6P-TOR drives dynamic cell wall reorganization, providing the necessary mechanical support and structural foundation for organogenesis^[22,23]. Bulbils serve as asexual reproductive organs in lily (Lilium spp.). Their formation depends on active cell division and dynamic cell wall remodeling; however, it remains unclear how carbohydrate metabolism drives bulbil development by regulating cell wall synthesis. The biosynthesis of cell wall polysaccharides is controlled by glycosyltransferase enzymes, with cellulose synthase (CESA) and cellulose synthase-like (CSL) enzymes serving as the key catalytic components that synthesize the β-1,4-linked glycan backbones of cellulose and hemicellulose, respectively^[24]. The CSL superfamily mediates the biosynthesis of various plant cell wall polysaccharides, including cellulose and hemicellulose^[25]. Its subfamilies (CSLA to CSLH) each contain multiple genes, and these genes have different but evolutionarily related functions^[26]. In rice (Oryza sativa), OsCSLD4 plays a critical role in cell wall formation by modulating xylan and cellulose content, and further regulates grain organ size through the control of cell division and expansion^[27]. In barley (Hordeum vulgare), HvCslF3 mediates the synthesis of cell wall polysaccharides and markedly influences root growth and differentiation^[28]. Although numerous CSL genes have been functionally characterized across diverse plant species, whether these genes participate in regulating the initiation and development of axillary meristems (including bulbil formation) remains unclear.

In this study, transcriptome sequencing was performed on both natural bulbil-forming lily cultivars ('Hotel California' and 'Flore Pleno') and non-bulbil-forming cultivars ('Red Velvet' and 'Brasilia'). We revealed that differentially expressed genes (DEGs) associated with bulbil formation are predominantly enriched in carbohydrate metabolism and cellular development. Several candidate genes involved in sucrose and starch metabolism, the TCA cycle (Tricarboxylic acid cycle), pyruvate metabolism, cell cycle, and proliferation were identified and validated via qRT-PCR. Further functional validation confirmed that carbohydrate metabolism-related genes, particularly LlCSLE1, play a critical role in promoting cellulose accumulation during bulbil formation. These findings provide new molecular insights into lily bulbil formation, underscoring the importance of dynamic changes in cell wall composition. This study establishes a valuable theoretical foundation for understanding the development of plant vegetative reproductive organs and supports molecular breeding efforts in lilies.

The lily cultivars 'Flore Pleno', 'Hotel California', 'Red Velvet', and 'Brasilia', along with bulbs of Lilium lancifolium, were selected as the experimental materials. All specimens were obtained from the National Lily Germplasm Resource Bank at the Beijing Academy of Agricultural and Forestry Sciences and were cultivated at the Academy's glass greenhouse facility. The plants were maintained in a greenhouse under the following controlled conditions: diurnal temperatures of 18–25 °C, relative humidity of 70%–80%, and a 16 h/d photoperiod. Continuous monitoring of bulbil development was conducted post-planting. A standardized sampling procedure was implemented uniformly upon initial observation of bulbil initiation in varieties characterized by natural bulbil-forming capability. Samples were collected from the upper leaf axils of the 'Flore Pleno' and 'Hotel California' varieties, which naturally form bulbils, and the 'Red Velvet' and 'Brasilia' varieties, which do not.

RNA-seq and detection of relative expression levels of DEGs were performed on axillary tissues of four lily cultivars: 'Hotel California', 'Flore Pleno', 'Red Velvet', and 'Brasilia'. For expression pattern analysis and VIGS experiments, Lilium lancifolium plants were used.

Sampling was conducted using three biological replicates, with each replicate consisting of a pooled sample randomly selected from nine plants. For each plant, three adjacent leaf axils from the same location were collected, immediately placed in liquid nitrogen, and stored at −80 °C for transcriptome sequencing and gene expression analysis.

The formation of bulbils in Lilium lancifolium was divided into three stages: S0 (early vegetative growth stage), S1 (pre-bulbil formation stage), and S2 (bulbil formation stage). S0: leaf axils flat, no swelling; S1: axils enlarged but no visible bulbil primordia; S2: distinct bulbil primordia visible. Tissue samples were collected from the upper leaf axils at each stage, with three biological replicates per stage (each containing leaf axils from nine plants). These samples were employed for the determination of cellulose content, hemicellulose content, and gene expression levels.

For all lily samples, leaf axils measuring 4.5–6.5 mm in length were collected from the upper stem segments, regardless of cultivar or experimental condition.

RNA-seq data processing and de novo assembly

Four varieties used for transcriptome sequencing were sampled at the same developmental stage. After total RNA was qualified, mRNA was enriched with Oligo magnetic beads and fragmented to create cDNA libraries by adding fragmentation buffer. The insert size of the library was assessed using the Agilent 2100 Bioanalyzer (Agilent Technologies, USA) to ensure library quality^[29]. All libraries were sequenced using the paralog-free transcriptome approach on the Illumina NovaSeq platform (Novogene, Beijing, China). The sequenced reads (raw reads) were filtered by removing low-quality reads (NCBI SRA accession number: PRJNA1320984) to obtain clean reads. The Trinity program was used to assemble the clean reads^[30]. After clustering de novo assembled transcripts and counting overlapping reads, unigenes were obtained for further analysis.

In the original analysis, differentially expressed genes were screened using the thresholds of |log₂ (fold change) | ≥ 1 and p < 0.05, as calculated by the OmicShare platform. All data applied FDR correction.

Gene annotation, enrichment, sequencing process, and differential expression analysis

The unigenes were annotated by BLAST against the GenBank non-redundant, Protein family (Pfam), Swiss-Prot, euKaryotic Ortholog Groups (KOG), and Gene Ontology (GO) databases, respectively. Differential expression analysis of genes was further conducted by analyzing GO, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) functional enrichment^[31]. Heatmaps of gene expression, Venn diagram, Principal Component Analysis (PCA), GO, KEGG, and GSEA were performed using the Omicshare Tools, a free online platform for data analysis (www.omicshare.com/tools).

Total RNA extracted from leaf axil tissues of the four lily varieties was subjected to mRNA enrichment using Oligo magnetic beads, followed by fragmentation to construct cDNA libraries. Overlapping reads were clustered to obtain unigenes, which were subsequently used for functional annotation and differential expression analysis.

RNA extraction and quantitative real-time PCR analysis

Following the manufacturer's instructions, total RNA was extracted using the FastPure Plant Total RNA Isolation Kit (Vazyme, China). First-strand cDNA was synthesized according to the user manual of the HiScript III 1^st Strand cDNA Synthesis Kit with gDNA Eraser (Vazyme, China). Quantitative real-time PCR (qRT-PCR) analysis was employed to verify the expression levels of selected genes from the RNA-seq data. The specific primers for qRT-PCR were generated using the online GenScript Real-Time PCR (TaqMan) Primer Design Tool (www.genscript.com/tools/real-time-pcr-taqman-primer-design-tool#). The primer sequences are listed in Supplementary Table S1. 'Flore Pleno' (FP) rRNA was used as the internal reference gene^[32,33]. All qRT-PCR assays were performed with at least three biological replicates. qRT-PCR was conducted using TB Green® Premix Ex Taq™ II (TaKaRa, Japan) in the Bio-Rad CFX96TM Real-Time System (CA, USA).

Determination of cellulose content

Cellulose content was quantified using a Cellulose Assay Kit (BC4280, Solarbio, Beijing, China) according to the manufacturer's instructions. Briefly, cellulose was hydrolyzed into glucose by heating under acidic conditions. The monosaccharides were then dehydrated in concentrated sulfuric acid to generate furfural derivatives, which were then quantified by measuring the blue-green color reaction with anthrone reagent.

Determination of hemicellulose content

Hemicellulose content was determined using a Hemicellulose Assay Kit (BC4445, Solarbio, Beijing, China) following the manufacturer's protocol. Briefly, hemicellulose was hydrolyzed into reducing sugars under the specified conditions. The reducing sugars were then oxidized to sugar acids and other byproducts in an alkaline environment by 3,5-dinitrosalicylic acid (DNS) upon heating, while DNS was reduced to reddish-brown 3-amino-5-nitrosalicylic acid.

Virus-induced gene silencing

The expression of the LlCSLE1 gene was silenced through virus-induced gene silencing (VIGS)^[34]. A 200-nt fragment of LlCSLE1 was amplified using PCR and inserted into the pTRV2 vector to generate the pTRV2-LlCSLE1 construct. Lily bulbs were immersed in the infiltration buffer containing Agrobacterium cells transformed with equal amounts of pTRV1 and either pTRV2 or pTRV2 containing target gene fragments. Lily bulb scales immersed in the bacterial suspension were vacuum-infiltrated at 0.8 MPa for 5 min. This process was repeated twice to improve the infection efficiency. After infiltration, the materials were planted in a greenhouse at 22 °C with a relative humidity of 60%–70%.

Statistical analysis

Three biological replicates were analyzed for each condition. GraphPad Prism (version 7) was used for statistical analysis. Comparisons among three developmental stages were performed using one-way ANOVA followed by Duncan's post-hoc test, with significant differences indicated by letters (a, b, c) in the figures. For comparisons involving only two groups, Student's t-test was applied.

The leaf axils of the lily stem segments were classified into two categories: varieties naturally capable of bulbil formation ('Hotel California' and 'Flore Pleno') and those naturally incapable ('Brasilia' and 'Red Velvet'). Among the bulbil-forming varieties, 'Hotel California' produced an average of 30.18 bulbils, whereas 'Flore Pleno' produced an average of 6.09 (Fig. 1a). Counting the locations of bulbil emergence revealed that bulbils from both varieties primarily appeared in the axils of the upper-middle leaves. In addition, 'Hotel California' (HC) had a greater number of bulbils compared to FP (Fig. 1b).

Figure 1. 

Formation process of different lily bulbil varieties. (a) The number of bulbils was counted when the bulbs of four varieties of lily were planted for 50 d. Scale bar indicates 1 cm. The data are expressed as mean ± SD (n = 11). (b) Statistics of the occurrence position of bulbils for lily varieties that naturally occur in bulbils. L16 represents the 16^th leaf counted from the bottom of the stem segment, and so on, and L128 represents the 128^th leaf counted from the bottom of the stem segment upwards. (c) Phenotypic comparison of leaf axils in bulbil-forming and non-bulbil-forming lily varieties. HC, 'Hotel California'; FP, 'Flore Pleno'; BRA, 'Brasilia'; RV, 'Red Velvet'.

Photographic observations were conducted to further investigate the differences in traits between the varieties. Close-up images of the leaf axils clearly show the presence of roots, scales, and anthocyanin accumulations in the bulbils (Fig. 1c).

Identification of differentially expressed genes related to the formation of lily bulbil

To identify genes potentially involved in bulbil formation, we examined the expression differences of genes potentially associated with bulbil formation in lily. We employed the fragments per kilobase of transcript per million mapped reads (FPKM) method in the RNA-seq analysis to estimate gene expression levels in varieties naturally capable (CBF) and incapable (IBF) of bulbil formation.

The number of upregulated genes was observed to be similar to the number of downregulated genes in the RV-up vs FP-up and RV-up vs HC-up groups, with 23,209 upregulated and 21,670 downregulated genes for the former, and 22,613 upregulated and 25,067 downregulated genes for the latter. In contrast, the upregulated gene counts for the BRA-up vs FP-up and BRA-up vs HC-up groups were 22,926 and 11,968, respectively (Fig. 2a).

Figure 2. 

DEGs between the axils of the upper leaves of different lily varieties. (a) Number of DEGs compared in pairs. (b) Venn diagram showing the number of DEGs among different comparisons. (c) PCA between samples. (d–g) Volcano plots showing DEGs with common trends in different comparisons. (h) Heatmap showing the relative expression levels of DEGs with 5,121 annotations. The blue-to-red color scale on the right indicates the FPKM values from low to high.

Venn diagrams were generated, and principal component analysis (PCA) was conducted on the DEGs on 5,121 individual genes with functional annotations involved in the bulbil-forming process across the four varieties (Fig. 2b, c). A total of 111,978 DEGs were identified in the pairwise comparisons, with 44,878, 51,679, 68,775, and 43,696 DEGs found in the comparisons of RV-up vs FP-up, RV-up vs HC-up, BRA-up vs FP-up, and BRA-up vs HC-up, respectively (Fig. 2b). PCA revealed similarities in DEG expression between HC and FP and more distinct differences between CBF and IBF. Hierarchical clustering separated the bulbil-forming varieties (HC and FP) from the non-bulbil-forming ones (BRA and RV), which was further confirmed by the PCA plot showing that the two groups formed distinct clusters (Fig. 2c). This indicates that the identified DEGs are predominantly associated with bulbil-forming capacity rather than varietal background.

We compared the expression of these 5,121 genes across the groups and generated volcano plots, revealing a similar number of significantly differentially expressed genes between groups. The volcano plots show that a comparable number of significantly upregulated and downregulated genes existed between the comparisons, suggesting that both activation and repression of specific genes contribute to bulbil formation (Fig. 2d–g). Figure 2h presents the heatmap of the 5,121 DEGs, indicating their expression profiles and the hierarchical clustering results. Most DEGs exhibited opposite expression patterns between CBF and IBF. Overall, 71.76% of DEGs showed a high number of transcripts in CBF and a low number in IBF.

GO annotation and KEGG pathway analysis

A total of 5,121 unigenes were aligned with the GO database and classified into the categories of Biological Process, Molecular Function, and Cellular Component (Supplementary Fig. S1). The results indicate notable enrichment of genes in Biological Process, whereas Molecular Function and Cellular Component contained fewer enriched genes (2,375 and 2,788, respectively). Within the Biological Process category, cellular process and metabolic process were notably enriched, comprising 1,037 and 900 DEGs, respectively. Binding was the most enriched category in Molecular Function, with 1,029 DEGs identified. These findings suggest that during lily bulbil formation, significant differences in genes related to cell division and metabolic processes may play a critical role in initiating bulbil development (Supplementary Fig. S2a).

KEGG pathway analysis was conducted to describe the metabolic pathways associated with the DEGs^[31]. Supplementary Fig. S2b presents the top 19 KEGG pathway enrichment results. Metabolism had the highest enrichment of DEGs, followed by genetic information processing, with 878 DEGs enriched in metabolic pathways and 137 in genetic information processing pathways. A total of 102 DEGs were relatively enriched in carbohydrate metabolism. In addition, three pathways related to biological metabolism—amino acid metabolism, lipid metabolism, and biosynthesis of other secondary metabolites—were also relatively enriched, indicating that these genes may play a role in energy transformation and the synthesis of cell components during bulbil formation.

To further illustrate the expression trends of DEGs, we analyzed the FPKM values for the 5,121 DEGs across the four varieties. Supplementary Fig. S3 depicts the corresponding trend analysis plot. DEGs with the highest enrichment showed moderate down-regulation, suggesting that a significant number of DEGs were upregulated during bulbil formation in CBF varieties (Supplementary Fig. S3).

We subsequently performed KEGG pathway analysis on the DEGs identified in the trend analysis. The results aligned well with the previous KEGG results (Supplementary Fig. S2b). Metabolism exhibited the highest enrichment of DEGs, with carbohydrate metabolism being the most enriched among the three groups, comprising 24, 29, and 12 DEGs, respectively.

Expression profiling of DEGs associated with carbohydrate metabolism

Based on the results of the GO and KEGG analyses, we further investigated the DEGs enriched in the carbohydrate metabolism pathway. A total of 102 DEGs were identified from the RNA-seq data. After examining the 102 DEGs enriched in carbohydrate metabolism, we selected those with the highest enrichment in key pathways (propanoate metabolism, TCA cycle, pentose and glucuronate interconversions, starch and sucrose metabolism, and pyruvate metabolism). Among these, 30 genes with significant differential expression, clear functional annotations, and reliable expression levels were chosen for heatmap visualization in Fig. 3, representing the core carbohydrate metabolic branches most relevant to bulbil formation. Supplementary Fig. S4 presents the heatmap of the expression patterns of the carbohydrate-related DEGs. The enrichment of DEGs in carbohydrate metabolism pathways (Fig. 3) underscores the significance of carbohydrate metabolism during bulbil formation.

Figure 3. 

Heatmap of the relative expression levels of DEGs potentially associated with bulbil formation. Carbohydrate metabolism pathways related to bulbil formation. A total of 30 significantly differentially expressed genes were found to be involved in carbohydrate metabolism during bulbil formation. Red indicates an upward adjustment, and blue indicates a downward adjustment. Glucose-6-P; fructose-6-P; FBP; glyceraldehyde-3-P; 6-P-gluconolactone; 6-P-gluconate; ribose-5-phosphate; PEP; oxaloacetate; acetyl CoA; citrate; α-Ketoglutarate; succinate; propanoate; malate.

Notably, the heatmap reveals significant differences in the expression of most DEGs in 'Hotel California', suggesting that they may directly or indirectly influence bulbil number (Fig. 3). Among the carbohydrate metabolism-related DEGs shown in Fig. 3, genes involved in starch and sucrose metabolism, the TCA cycle genes including phosphoenolpyruvate carboxykinase PCK, and pyruvate metabolism genes exhibited consistently higher expression levels in the bulbil-forming varieties, particularly HC, compared to the non-bulbil-forming varieties. In contrast, genes associated with pentose and glucuronate interconversions showed relatively lower expression in the bulbil-forming group. This distinct expression pattern highlights that the activation of specific carbohydrate metabolic pathways, especially those supplying energy and carbon skeletons, is closely linked to bulbil formation. Recent studies indicate that the formation and development of lily bulbils are closely related to carbohydrate metabolism, particularly starch and sucrose metabolism^{[12,24,35,36]}. Furthermore, O-glycosidic biomolecules constitute the cell wall, an important cellular exoskeleton, primarily composed of cellulose, pectin, hemicellulose, and lignin. Bulbil formation involves cell proliferation and differentiation, accompanied by dynamic cell wall modifications that facilitate cell division and growth.

Expression profiling of DEGs associated with cell proliferation and division

To comprehensively analyze the expression data of DEGs in the transcriptome, we identified DEGs with gene name annotations and ultimately screened 2,485 DEGs for GSEA. We compared annotated DEGs using the same paired grouping described previously.

GSEA was performed to assess differences in signaling pathways between the CBF and IBF groups. Our results indicate that cell division and cell proliferation were significantly enriched in the CBF group (Supplementary Fig. S5). The gene set associated with cell proliferation was particularly enriched in GSEA, emphasizing changes in the expression of genes related to cell proliferation during bulbil formation. Specifically, 77 genes were enriched in the cell division pathway, with 27 showing significant enrichment. In addition, 43 genes were enriched in the cell cycle regulatory pathway, with 22 of these showing significant enrichment.

GSEA revealed cyclins to be notably enriched in the cell division and proliferation pathways, exhibiting significant upregulation in the CBF group. Cyclins accounted for 17% and 27% of the DEGs in the two pathways, respectively (Supplementary Fig. S6b, d).

Validation of DEGs using qRT-PCR analysis

To verify the authenticity of the RNA-seq data, we selected 20 genes related to carbohydrate metabolism, cell proliferation, and cell division for qRT-PCR analysis. These genes included the transcription factors MYB3R4, SCL28, and WOX4, as well as beta-1,4-glucanase (At4g02290), sucrose synthase (SUS4, SUS5), cytosolic phosphorylase (PHS2), phosphoenolpyruvate carboxykinase (PCK), cellulose synthase-like (CSLE1), sucrose-phosphate phosphatase (SPP2), β-xylosidase (Xyl3A), cyclins (CYCB, CYCD), Kip-related protein (KRP4), and histone acetyltransferase (HAM1).

Most selected genes showed consistent expression trends between qRT-PCR and RNA-seq data, except for SUS4 and SUS5. The opposite trends observed for these two genes may be explained by functional divergence of sucrose synthase isoforms: SUS catalyzes the reversible conversion of sucrose and a nucleoside diphosphate to the corresponding nucleoside diphosphate-glucose and fructose, and different isoforms (SUS4 and SUS5) may play distinct roles during bulbil formation (Fig. 4).

Figure 4. 

Expression analysis of selected unigenes in leaf axils of different lily cultivars. (a) Unigenes of carbohydrate metabolism pathways. (b) Unigenes of cell proliferation and division pathways. Blue and purple columns represent the CBF and IBF groups, respectively. The qRT-PCR values were detected using the 2^−ΔΔCᴛ method. Error bars show the standard deviation (SD) of three independent biological replicates. Lowercase letters indicate significant differences (p < 0.05).

Changes in cell wall composition and associated gene expression in the leaf axil during bulbil formation

Among the candidate genes determined to be associated with lily bulbil formation, we identified the rarely reported LlCSLE1, which is the focus of this study. We conducted an in-depth investigation into its specific molecular mechanisms and biological functions during bulbil formation in Lilium lancifolium.

To investigate the cytological basis of bulbil formation, we determined the content of major cell wall components (cellulose and hemicellulose) in leaf axil tissues at different developmental stages and analyzed the expression pattern of LlCSLE1 associated with this process. The cellulose content in the leaf axil exhibited a progressive and significant increase with bulbil development (Fig. 5a, b). At the S0 stage, the cellulose content was relatively low and increased markedly as development progressed, reaching its peak at the S2 stage. This demonstrates that the biosynthesis and accumulation of cellulose represent a key event in the morphogenesis of bulbils. In contrast to the change trend of cellulose, the hemicellulose content remained stable with no significant changes during the S0 and S1 stages (Fig. 5c). However, at the S2 stage, the hemicellulose content rapidly peaked and was significantly higher than that in the stages prior to bulbil formation. The expression level of LlCSLE1 in the upper leaf axils of Lilium lancifolium during bulbil formation was determined by qRT-PCR. The results indicate that the expression pattern of LlCSLE1 from stage S0 to S2 was highly synchronized with the changes in cellulose content, exhibiting a parallel increasing trend and reaching a maximum at the bulbil formation stage (Fig. 5d).

Figure 5. 

Characteristics of the cell wall components and LlCSLE1 expression in L. lancifolium at different bulbil development stages. (a) Morphological observations of leaf axils at different stages of bulbil formation. S0, S1, and S2 represent three distinct growth phases. Scale bar indicates 1 cm. (b) Cellulose content during different developmental stages of bulbil formation in Lilium lancifolium. (c) Hemicellulose content during different developmental stages of bulbil formation in Lilium lancifolium. (d) Relative expression levels of LlCSLE1 during different developmental stages of bulbil formation in Lilium lancifolium. Different letters (a, b, c) indicate significant differences among groups (p < 0.05).

These results reveal that during bulbil formation, the contents of both cellulose and hemicellulose increase significantly in the leaf axil. Furthermore, the expression of the key cellulose synthase-like gene, LlCSLE1, is synchronously upregulated with cellulose accumulation. This suggests that LlCSLE1 plays an important role in regulating the cell wall construction and reinforcement required for bulbil initiation.

Loss-of-function of LlCSLE1 inhibits bulbil formation and reduces cellulose content in lily

To examine whether the observed correlation implies causality, we specifically silenced LlCSLE1 in the bulb of Lilium lancifolium using VIGS to investigate the impact of its loss-of-function on bulbil formation and cell wall composition. Compared with the control group, the expression level of LlCSLE1 in the leaf axils of the experimental group was significantly suppressed (Fig. 6a, b). Phenotypic observations revealed that gene silencing significantly reduced the number of bulbils formed (Fig. 6c). Furthermore, the determination of the cell wall components in the leaf axils of silenced plants revealed a significant decrease in cellulose content compared with the control group (Fig. 6d). This is directly correlated with the suppression of LlCSLE1 expression. No significant differences were observed in hemicellulose content between the experimental and control groups (Fig. 6e).

Figure 6. 

LlCSLE1 affects the formation of bulbils in L. lancifolium. (a) Middle and upper stem segments of silenced and TRV2 plants form a bulbil phenotype after 60 d. Leaves were removed, leaving only the petiole base to provide a clearer view of the phenotype. This manipulation was applied equally to both groups and performed only at the endpoint, without affecting normal growth or bulbil development during the experimental period. At least 13 independent plants in the TRV2 or silenced groups were used. White asterisks denote the bulbils. The scale bar indicates 1.0 cm. (b) Relative expression of LlCSLE1; ** indicates significant differences at p < 0.01. (c) Bulbil formation amount, * indicates significant differences at p < 0.05. (d) Cellulose dry weight content in the leaf axil, * indicates significant differences at p < 0.05. (e) Hemicellulose dry weight content in the leaf axil, 'ns' indicates no significant differences at p > 0.05. The data are expressed as mean ± SD of three biological replicates. Student's t-test was used for statistical analysis.

These results demonstrate that LlCSLE1 is crucial for cellulose biosynthesis during lily bulbil formation. Its functional loss specifically impedes cellulose accumulation, thereby inhibiting bulbil development without affecting hemicellulose biosynthesis, providing genetic evidence that LlCSLE1 positively regulates this process by modulating cellulose biosynthesis.

The formation process of lily bulbils involves complex metabolic pathways and regulatory mechanisms^[24]. Among these, carbohydrate metabolism is crucial for the formation, growth, development, and yield of bulbils. In our functional analysis, genes related to carbohydrate metabolism were significantly enriched (Supplementary Fig. S2), exhibiting differential expression across various pathways and general upregulation in the CBF group (Fig. 3). This enhancement directly supports the increased demand for carbon skeletons and energy required for cell wall synthesis, particularly cellulose accumulation, during bulbil initiation.

Previous studies have shown that sucrose synthase (SuSy) activity influences plant growth, shoot apical meristem function, and leaf morphology^[37]. In the phloem, SuSy helps to maintain the balance between sucrose and its breakdown products, providing hexoses for energy production in companion cells and serving as substrates for complex carbohydrates^[38]. The overexpression of ARK1 in poplar stems inhibits leaf differentiation and internode elongation, with 41% of genes in the stem, including At4g02290 and SUS5, showing marked upregulation. These genes encode proteins involved in the synthesis or modification of the extracellular matrix^[39]. In our study, the upregulation of At4g02290 and SUS5 in bulbiferous varieties (Fig. 4) suggests that sucrose metabolism is specifically activated to provide UDP-glucose, the direct substrate for cellulose synthesis, thereby linking carbohydrate flux to the subsequent action of cellulose synthase-like proteins such as LlCSLE1. Thus, our results indicate that, in pathways involving sucrose and starch metabolism, the three homologous genes At4g02290, a carbonyl hydrolase, and SUS5, a sucrose synthase, were significantly upregulated in 'Hotel California' and 'Flore Pleno' (Fig. 4). These genes may play key roles in carbohydrate metabolism and energy conversion processes related to bulbil formation, potentially channeling carbon toward cell wall polysaccharide biosynthesis. Carbohydrate metabolism, particularly sucrose and starch, is essential for the growth of lotus rhizomes, potatoes, and cassava tubers^[40−42]. This study emphasizes the significance of sucrose and starch metabolism in the formation of bulbils, supporting previous studies that confirmed the roles of sucrose in osmoregulation, carbon sourcing, and signal transduction during bulbil formation, ultimately promoting bulbil development^[24,36]. GSEA also revealed significant enrichment of cell cycle and proliferation pathways in bulbil-forming varieties (Supplementary Fig. S5, S6), indicating that active cell division accompanies carbohydrate metabolism. The TCA cycle genes, including PCK, were upregulated in these varieties (Fig. 3), suggesting enhanced energy supply is required for bulbil formation, qRT-PCR validation confirmed the transcriptomic trends (Fig. 4).

In addition, our study revealed the significant expression of the TCA cycle in the transcriptome (Fig. 3). The TCA cycle is a major component of mitochondrial metabolism in eukaryotes, including higher plants. Our results align with previous findings that suggest that the TCA cycle plays a crucial role in the transition from the vegetative to reproductive stage, with its specific interactions potentially related to the metabolic and energy requirements of the buds^[43−45]. The TCA cycle generates ATP and carbon intermediates that are essential for the energy-intensive process of cellulose biosynthesis and cell wall remodeling during bulbil formation, providing the necessary metabolic support for LlCSLE1 function. Most studies on the initiation of lateral tissues focus on plant hormones, carbon sources, sugar signals, or exogenous stimuli^[13]. However, research on the TCA cycle and mitochondrial formation remains limited. Nonetheless, our findings suggest that TCA cycle activity, together with sucrose metabolism, contributes to the metabolic environment that enables LlCSLE1-mediated cellulose accumulation and bulbil development.

In the context of sucrose metabolism, we noted a discrepancy between qRT-PCR and RNA-seq data for two sucrose synthase isoforms, SUS4 and SUS5. While the transcriptome suggested elevated expression of both genes in bulbil-forming varieties, qRT-PCR did not fully corroborate this trend. Such inconsistency may arise from differences in primer sensitivity, transcript isoform detection, or sequence polymorphisms among the four lily varieties used for RNA-seq, given that the reference transcriptome was assembled de novo from the same samples. Post-transcriptional regulations such as alternative splicing or differential RNA stability could also contribute to platform-dependent quantification biases.

LlCSLE1 positively regulates the formation of bulbils in lilies, a process mediated through its promotion of cellulose accumulation

Stage-specific analysis of bulbil formation in Lilium lancifolium and functional validation via gene silencing revealed the important role of LlCSLE1 in bulbil development. Our findings provide novel insights into the molecular mechanism underlying bulbil initiation and underscore the close association between dynamic changes in cell wall composition and the formation and development of bulbils.

First, we established a strong link between the development of bulbils-an asexual reproductive organ-and dynamic changes in cell wall composition. Cellulose and hemicellulose, as major structural polysaccharides of the plant cell wall, play a crucial role in determining cell wall architecture and function^[46]. Changes in their content directly influence the structural integrity and physiological properties of the cell wall. We observed a continuous increase in cellulose content in the leaf axils as bulbil development progressed, accompanied by a significant accumulation of hemicellulose during the bulbil formation stage (Fig. 5b, c). These results strongly suggest that leaf axillary cells actively synthesize and deposit cell wall materials, providing the structural foundation and mechanical support for bulbil morphogenesis. Bulbil formation, as a specialized mechanism of asexual propagation in plants, requires rapid cell division and differentiation. The dynamic remodeling of the cell wall, including polysaccharide accumulation, serves as a critical foundation supporting this process^[11]. In addition, other plants, such as potato tuberization^[47] and onion bulb development^[48], exhibit this process with high similarity to de novo organogenesis. In these developmental processes, the accumulation of cellulose and hemicellulose is also accompanied by an increased rate of cell proliferation, thereby providing essential structural support and material reserves for organ formation. However, cellulose and hemicellulose exhibited distinct accumulation patterns, indicating that although both are cell wall polysaccharides, they may play different physiological roles during bulbil formation.

Second, we demonstrated that the cellulose synthase-like gene LlCSLE1 is required for cellulose accumulation and positively regulates bulbil formation, suggesting it plays an important role in the cellulose biosynthetic pathway during lily bulbil development. The CSL gene family constitutes a large family and serves as a core regulator of plant cell wall polysaccharide biosynthesis. Among them, CSL genes have been frequently reported to be involved in the synthesis of hemicelluloses such as xyloglucan and mannan. For example, AtCSLD in Arabidopsis thaliana regulates the synthesis of xylan, a major hemicellulosic component, thereby influencing the elongation and differentiation of root tip cells^[49]. In this study, the expression pattern of LlCSLE1 was found to be closely synchronized with both cellulose accumulation and bulbil formation. More importantly, silencing LlCSLE1 via VIGS not only directly led to a reduction in cellulose content but also resulted in a significant decrease in bulbil number, while hemicellulose levels remained unaffected (Fig. 6). These results clearly demonstrate that LlCSLE1 specifically regulates cellulose biosynthesis rather than hemicellulose formation, and that the cellulose synthesis mediated by this gene serves as a critical limiting factor in bulbil formation. Insufficient cellulose production may lead to structurally unstable cell walls in the leaf axil, which are unable to support the initiation or sustained growth of the bulbil primordium, ultimately inhibiting bulbil development. This finding broadens our understanding of the functional scope of the CSL gene family, suggesting that certain CSL genes may be involved in regulating cellulose biosynthesis during specific organogenesis processes. This functional specificity raises interesting questions about the underlying biochemical mechanism. Although the present study does not resolve the precise molecular mechanism, it clearly establishes a necessary role for LlCSLE1 in cellulose deposition and bulbil formation, thereby providing a foundation for future investigations into the mechanistic divergence between CSL and CESA functions. Several hypothetical models could be considered: LlCSLE1 might directly catalyze β-1,4-glucan chain elongation, similar to canonical cellulose synthases (CESAs), or it might interact with other cell wall biosynthetic complexes to channel substrate flux toward cellulose. Another possibility is that LlCSLE1 acts as a regulator of CESA activity in the leaf axil. Further biochemical and genetic studies, for example, in vitro enzyme activity assays or protein-protein interaction analyses, will be required to test these hypothetical mechanisms and determine which of them operates during bulbil formation. Although the present study does not resolve the precise molecular mechanism, it clearly establishes a necessary role for LlCSLE1 in cellulose deposition and bulbil formation, thereby providing a foundation for future investigations into the mechanistic divergence between CSL and CESA functions. Similar functions have been reported for CSL genes in other species. For example, transient overexpression of ZjCSLA1 in jujube fruit resulted in increased cellulose content, indicating that ZjCSLA1 plays a crucial role in cellulose biosynthesis during fruit development^[50]. The OsCSLF6 gene in rice (Oryza sativa) is directly involved in cellulose biosynthesis, and the upregulation of its expression significantly increases cellulose content in the cell walls of stem tissues^[51]. However, these examples mostly involve CSLA or CSLF subfamilies, whereas our study adds a new dimension by implicating an E-subgroup CSL gene in cellulose accumulation during an asexual reproductive process. The unaltered hemicellulose content following gene silencing, combined with the reduced bulbil number, suggests that the accumulation of hemicellulose during the bulbil formation stage is likely a consequence rather than a cause of bulbil development. The temporal dynamics of cell wall components further support this interpretation. Cellulose content progressively increased throughout bulbil development, whereas hemicellulose content remained relatively stable at S0 and S1 but rose sharply only at S2 (Fig. 5b, c). This sequential pattern implies that cellulose accumulation provides early structural support for bulbil primordium initiation, while the late surge of hemicellulose may contribute to subsequent bulbil expansion and maturation. Notably, the sharp increase in hemicellulose at S2 aligns with the observation that silencing LlCSLE1 did not affect hemicellulose levels (Fig. 6e), further confirming that hemicellulose accumulation is not directly regulated by LlCSLE1 and likely occurs secondary to bulbil formation.

Based on these findings, we propose a working model: during bulbil formation, the LlCSLE1 gene is specifically activated and upregulated, promoting the rapid synthesis and accumulation of cellulose in the leaf axil. The accumulation of cellulose creates a necessary structural microenvironment for the initiation of the bulbil primordium and serves as an important driving force for bulbil development. In contrast, hemicellulose accumulation is likely regulated independently by other genes and supports subsequent bulbil expansion and maturation. However, this study still leaves some questions unanswered-how cellulose accumulation is precisely confined to the leaf axil and which upstream signaling pathways regulate this process remain to be elucidated.

Our study reveals key carbohydrate metabolic pathways involved in lily bulbil formation, highlighting the sucrose and starch metabolism pathways as well as the TCA cycle. Transcriptomic analysis indicated the significant enrichment of processes related to glucose metabolism, and cell proliferation and division. Quantitative real-time PCR results from these varieties validated the transcriptomic data and confirmed the functional role of the LlCSLE1 gene. We found that LlCSLE1 plays an indispensable role in bulbil formation by specifically regulating cellulose biosynthesis. This discovery broadens our understanding of the cell wall biological mechanisms underlying the development of plant asexual reproductive organs (Fig. 7).

Figure 7. 

Schematic diagram of the mechanism by which LlCSLE1-mediated cellulose synthesis promotes bulbil formation.