Comparative transcriptome analysis unravels distinct regulatory mechanisms governing fiber development in red silk cotton tree (Bombax ceiba)

Abstract

Bombax ceiba a multipurpose tree that has recently garnered increasing interest due to the unique characteristics of its fibers, which hold promise for commercial textile applications. These fibers, derived from the inner epidermis of the fruit, are fine, lightweight, highly hollow, and exhibit excellent thermal insulation. Notably, they lack natural twist—an uncommon feature among plant-derived fibers. While this feature contributes to their distinctive structure, it also results in low tensile strength, rendering them unsuitable for conventional spinning processes and limiting their utility in textile manufacturing. In this study, we performed a comprehensive comparative transcriptome analysis to elucidate the regulatory mechanisms underlying fiber development in B. ceiba. Our findings revealed a distinct transcriptional landscape in the inner fruit epidermis, characterized by the enrichment of genes involved in hemicellulose and lignin biosynthesis. Furthermore, we identified a markedly reduced expression of genes associated with microtubule (MT) organization and orientation, which may influence the patterning of cellulose microfibril deposition and contribute to the suppression of fiber twist. These molecular signatures likely underpin the unique morphology and mechanical properties of B. ceiba fibers. Altogether, our study provides novel insights into the genetic and molecular basis of fiber development in B. ceiba, offering a foundation for future strategies aimed at improving fiber quality and expanding the industrial application of this and other natural fiber-producing species.

Background

Bombax ceiba L. (Malvaceae), commonly known as the red silk cotton tree or kapok tree, is a fast-growing deciduous species widely distributed across tropical and subtropical regions, including Southeast Asia, Africa, and parts of Central America. It plays an important ecological role in maintaining biodiversity and supporting ecosystem stability [1]. In addition to its ecological importance, B. ceiba possesses significant economic value. Its timber is used in construction and furniture, while various parts of the plant serve medicinal and industrial purposes [1,2,3]. Notably, the tree produces unique fruit fibers that have attracted growing attention for their promising potential in sustainable material production. B. ceiba fibers are derived from the inner pericarp (epidermis) of mature fruit and exhibit remarkable physical characteristics, including high hollowness, low density, hydrophobicity, excellent thermal insulation, and buoyancy [4]. These properties have enabled their traditional use as stuffing in mattresses, pillows, upholstery, and life vests [5]. Compared to upland cotton (Gossypium hirsutum L.), B. ceiba fibers are contain lower cellulose (~ 64%) and higher proportions of lignin (~ 13%) and xylan, which may contribute to their relatively low tensile strength [6, 7]. This mechanical weakness renders B. ceiba fibers unsuitable for spinning, thereby limiting their large-scale application in textile industries. Although fiber strength is known to be influenced by secondary cell wall thickness and cellulose deposition [8], the molecular basis of these traits in B. ceiba remains poorly understood.

Over the past two decades, extensive research on upland cotton has uncovered numerous transcription factors and regulatory pathways involved in fiber development, particularly in the seed outer epidermis [9, 10]. Key regulators such as GhMYB25, GhHD1, and GhWRKY16 are essential for fiber initiation, while GhMYB2 and GhMYB3 are involved in fiber elongation [11, 12]. Moreover, plant hormones like auxin, gibberellins, and brassinosteroids have been shown to coordinate cell expansion and elongation during fiber growth [9]. Genes such as GhCESA4/7/8, GhKNL1 (KNOTTED1-LIKE), and GhFSN1 (fiber secondary cell wall-related NAC1) have been identified as central players in cellulose biosynthesis and secondary wall thickening, processes critical for achieving high fiber strength and quality [13,14,15]. Despite these advances in cotton, the fiber of B. ceiba develops in a distinct anatomical origin—the inner fruit epidermis—and the genetic and transcriptional regulation underlying its formation remains largely unexplored.

A better understanding of B. ceiba fiber development is essential for elucidating the diversity of plant fiber formation and for improving the industrial utilization of its unique fibers. In this study, we conducted a time-series transcriptome analysis of developing B. ceiba fibers to investigate the molecular regulatory networks involved in their formation. Our findings reveal a distinct transcriptional program in the inner fruit epidermis, highlighting differences in cell wall biosynthesis, microtubule orientation, and fiber architecture when compared with cotton. These insights expand our understanding of fiber diversity in plants and lay the groundwork for potential genetic improvement and sustainable application of B. ceiba fibers in eco-friendly industries. Given their natural properties and biodegradability, B. ceiba fibers also hold promise for future development as eco-friendly alternatives to synthetic fibers in packaging, insulation, and absorbent materials. With advances in fiber engineering and composite technology [16], there is growing potential to modify or enhance B. ceiba fiber properties for use in biodegradable textiles, non-woven fabrics, and green construction materials. Exploring their integration into circular bioeconomy frameworks may further promote the sustainable utilization of this underexploited natural resource.

Materials and methods

Plant materials

Wild B. ceiba trees were sampled from Yuanmou Township, Chuxiong Prefecture, Yunnan Province, China (coordinates: E101°53′27.76″, N25°40′50.0″). The species was taxonomically identified by Professor Bin Tian of Southwest Forestry University. A voucher specimen has been deposited in the Herbarium of Southwest Forestry University under the accession number TianB-2021133-SWFC. As B. ceiba is not an endangered or protected species in China, no specific permissions were required for the collection of plant materials, and the field sampling did not involve privately owned or protected land. For transcriptome analysis, tissues were collected at five key developmental stages: ovary walls at 0 days post anthesis (DPA), developing fibers at 10, 20, 35, and 45 DPA, and ovary tissue lacking fiber at 45 DPA (Hereafter referred to as 45 ODPA). Sampling was conducted during the flowering and early fruit development season, from February to April 2018. Each developmental stage included three biological replicates, with each replicate representing a separate, healthy tree. All harvested samples were immediately frozen in liquid nitrogen on-site and subsequently stored at −80℃ until RNA extraction.

Transcriptome profiling

Total RNA was extracted using the RNAprep Pure Tissue Kit (Tiangen, Beijing, China) from all samples collected at six developmental stages (0, 10, 20, 35, and 45 DPA, and 45 ODPA), each with three biological replicates. RNA quality and integrity were assessed using the Agilent 2100 Bioanalyzer System, and only samples with high RNA Integrity Number (RIN) scores were selected for transcriptome sequencing. Library construction and sequencing were performed by NextOmics Biosciences Co., Ltd (Wuhan, China) on the Illumina HiSeq 4000 platform. Paired-end sequencing (2 × 150 bp) was conducted, generating an average of approximately 20–30 million clean reads per sample. Raw reads were subjected to quality control using fastp (v0.20.0) [17] to remove adapters and low-quality reads, yielding high-quality clean reads for downstream analyses. These clean reads were aligned to the newly assembled B. ceiba v2.0 reference genome (accession number GWHFSMB00000000.1), which has been deposited to the Genome Warehouse (GWH), using HISAT2 (v2.2.1) [18]. Gene expression quantification were performed using StringTie (v2.2.1) [19], and expression levels were subsequently analyzed with Ballgown (v2.36.0) [20].

To assess data quality and sample relationships, gene expression matrices were used to conduct correlation analysis and hierarchical clustering via the tidyverse (v2.0.0) package in R. Temporal expression dynamics were explored using time-series clustering with ClusterGVis (https://github.com/junjunlab/ClusterGVis), enabling the identification and visualization of stage-specific gene expression patterns. Functional enrichment analysis of Gene Ontology (GO; release 2024.11.03) terms was conducted using the clusterProfiler R package [21]. Gene co-expression networks were constructed for genes within each expression cluster based on the UMI (unique molecular identifier) count matrix, using GENIE3 (v1.26.0) [22], and the networks were visualized with Cytoscape (v3.9.1) [23].

Identification of fiber development-related genes and comparison expression levels analysis with upland cotton

To identify candidate genes involved in fiber development, protein sequences from the upland cotton and Arabidopsis were used as queries to perform homology searches against the B. ceiba protein database using BLASTP [24], with an E-value cutoff of < 1e-5. Expression data for homologous genes during fibers development in upland cotton were obtained from the COTTONGEN database (https://www.cottongen.org). Comparative expression profiles between B. ceiba and upland cotton were visualized using the heatmap package in R.

RNA isolation and qRT-PCR

Total RNA was isolated from developing fibers of B. ceiba at 20, 35, and 45 DAP, as well as from ovary tissue at 45 DPA (45 ODPA), using the RNAprep Pure Tissue Kit (Tiangen, Beijing, China), following the manufacturer's protocol. First-strand cDNA was synthesized from 1μg of total RNA using the TransScript All-in-One First-Strand cDNA Synthesis SuperMix for qRT-PCR (TransGen Biotech, Beijing, China). Quantitative real-time PCR (qRT-PCR) was conducted using the PerfectStart Green qPCR SuperMix Kit with the following thermal cycling conditions: an initial denaturation 95℃ for 5 min, followed by 35 cycles of denaturation at 95℃ for 5 s, annealing at 58℃ for 30 s, and extension at 72℃ for 30 s. A melting curve analysis was subsequently performed to confirm the specificity of the amplified products. The BcActin2 gene was used as an internal reference for normalization. Relative transcript levels were quantified using the 2−ΔΔCt method [25]. To ensure reliability and reproducibility, each sample was analyzed with three independent biological replicates, each comprising three technical replicates.

Results

Definition of B. ceiba fiber developmental stages

As previous reported [26], the development of B. ceiba fibers is closely associated with capsule growth, beginning with flower bud maturation (Fig. 1a), followed by pollination (Fig. 1b), ovary expansion (Fig. 1c), and subsequent capsule development and ripening (Fig. 1d-g). Based on our observations in combination with earlier studies [26], B. ceiba fiber development can be categorized into five successive stages: initiation (0–10 DPA; Fig. 1c-d), elongation and expansion (11–20 DPA; Fig. 1e), transition (21–35 DPA; Fig. 1f), secondary cell wall (SCW) thickening (36–45 DPA; Fig. 1g), and maturation (46–55 DPA; Fig. 1h). During the initiation phase, the inner surface of the ovary remains smooth, and the epidermal trichome cells lining this surface begin to differentiate into fiber cells [26]. In contrast to the fiber development observed in upland cotton, where each fiber arises from a single epidermal cell on the outer ovule surface, B. ceiba exhibits a distinct clustered initiation mechanism [26]. As demonstrated in earlier studies, trichome cells in B. ceiba undergo several rounds of cell division prior to fiber cell differentiation, giving rise to clusters of fiber cells rather than solitary fibers [26]. Beginning around 11 DPA, short fuzz fibers emerge on the inner ovary wall, Marking the onset of the elongation and expansion stage. This is followed by a period of rapid fiber cell elongation and enlargement that peaks at approximately 20 DPA. During the transition stage (21–35 DPA), fiber cells begin to increase in width, and cell wall thickening is initiated [27]. From 36 to 45 DPA, although fiber length and cell wall thickness show minimal further increase, the dry mass of the fiber rises steadily, suggesting active deposition of secondary wall materials (Fig. 1g). After 45 DPA, the fibers enter the Maturation phase, characterized by progressive dehydration and structural stabilization. By 55 DPA (Fig. 1h), mature B. ceiba fibers exhibit a smooth, fluffy texture, and a silk-like sheen, commonly earning them the nickname “soft gold”. Compared to upland cotton, B. ceiba fibers possess distinct biochemical and morphological features. Several reports have shown that B. ceiba fibers are enriched in non-cellulosic polysaccharides, particularly xylan (a major component of hemicellulose) as well as lignin [4, 28]. Morphologically, these fibers are typically shorter, more fragile, and contain a hollow lumen with relatively smooth surfaces [29]. In contrast, upland cotton fibers are longer and mechanically more robust, with well-developed SCWs that contribute to their superior textile quality [30]. To investigate gene expression dynamics across these B. ceiba fiber developmental phases, fiber tissues were harvested at 0, 10, 20, 35, and 45 DPA for RNA-seq analysis, representing key transitional points in B. ceiba fiber development (Fig. 1).

[IMAGE OMITTED: SEE PDF]

Transcriptomic dynamics during B. ceiba fiber development

To investigate gene expression dynamics during fiber development in B. ceiba, RNA-sequencing (RNA-seq) was performed on fiber tissues collected at five developmental stages (0, 10, 20, 35, 45 DPA), as well as on ovary tissue without fibers at 45 DPA (designated as 45 ODPA) as a control. Three biological replicates were prepared for each of the six conditions, resulting in the construction of 18 RNA-seq libraries. In total, 1,642.62 million paired-end reads were generated, of which approximately 94.61% were successfully aligned to the new reference genome B. ceiba v2.0 (accession number GWHFSMB00000000.1) (Table S1), indicating that the sequencing data were of high quality and suitable for downstream transcriptomic analyses.

To enhance the accuracy of gene expression profiling, genes with low expression (FPKM < 0.5 in at least one biological replicate) were filtered out, resulting in the removal of between 12,177 and 18,659 genes per sample. Among the retained genes, 34.81–40.41% were moderately expressed (1 ≤ FPKM < 10), 13.19–20.70% were highly expressed (10 ≤ FPKM < 100), and only a small fraction (1.60–1.93%) exhibited very high expression levels (FPKM ≥ 100) during fiber development (Fig. 2c). Notably, differential expression analysis was conducted through pairwise comparisons among the six tissue types using the criteria |log2 (fold change)|≥ 2 and false discovery rate (FDR) ≤ 0.05. A total of 11,103 non-redundant differentially expressed genes (DEGs) were identified across the fiber development stages. The largest number of DEGs were observed in the comparisons between 10 DPA vs. 45 ODPA (10,119 DEGs) and 0 DPA vs. 45 ODPA (9,104 DEGs), suggesting significant transcriptional reprogramming during early fiber initiation. In contrast, fewer DEGs were detected between 10 DPA vs. 20 DPA (1,270 DEGs) and 20 DPA vs. 35 DPA (1,919 DEGs), reflecting relatively stable transcriptomic profiles during mid-developmental transitions (Fig. 2a and Table S2). Analysis of DEG overlaps revealed shared transcriptional signatures between adjacent stages: 118 DEGs were common to both 0 DPA vs. 10 DPA and 10 DPA vs. 20 DPA; 88 between 10 DPA vs. 20 DPA and 20 DPA vs. 35 DPA; and 99 between 20 DPA vs. 35 DPA and 35 DPA vs. 45 DPA. Notably, 1,321 genes were co-regulated between 35 DPA vs. 45 DPA and 45 DPA vs. 45 ODPA, while 1,031 genes were shared between 45 DPA vs. 45 ODPA and 0 DPA vs. 45 ODPA. These patterns support the presence of five distinct yet overlapping phases of fiber development in B. ceiba (Fig. 1c).

[IMAGE OMITTED: SEE PDF]

Identification of spatial expression trends across fiber transcriptomes

To elucidate the molecular mechanisms governing fiber development in B. ceiba, from initiation to maturation, we performed K-means clustering on all expressed genes (14,239 genes with FPKM > 1 in at least one of the six sampled time points). This analysis grouped the genes into six distinct clusters (C1-C6), each exhibiting a stage-specific expression profile (Fig. S1, Table S3). The top three Gene Ontology (GO) terms enriched in each cluster are shown in Fig. 3.

[IMAGE OMITTED: SEE PDF]

Cluster 1 (C1) and 2 (C2) contained genes predominantly expressed at the early fiber developmental stages—0 and 10 DPA. C1, comprising 1,960 genes, is associated with the fiber initiation stage and enriched in GO terms related to early cellular responses, such as response to chitin (GO: 0010200, P value = 6.44E-36), CCR4-NOT complex (GO: 0030014, P value = 1.46E-05), and calmodulin binding (GO: 0005516, P value = 1.40E-05) (Fig. 3, Fig. S2). In contrast, C2 includes 3,367 genes expressed at a slightly later initiation phase and is enriched in GO categories related to ribosome assembly (GO: 0042255, P value = 2.35E-22), ATP binding (GO: 0005524, P value = 1.33E-09), and protein folding chaperone (GO: 0044183, P value = 7.49E-10), indicating a preparatory shift toward active fiber elongation (Fig. 3, Fig. S3). The C3, containing 2,527 genes, is predominantly expressed at 20 DPA, marking the fiber elongation and expansion phase. Genes in this cluster are enriched in functions related to chromatin remodeling and plastid activity, including nucleosome organization (GO: 0034728, P value = 1.24E-14), plastid thylakoid membrane (GO: 0055035, P value = 1.39E-23), and electron transfer activity (GO: 0009055, P value = 2.3 × 10−4) (Fig. 3, Fig. S4). C4, consisting of 1,814 genes, is characterized by preferential expression at 35 DPA and represents the transition phase from primary to SCW formation (Fig. 3). Enriched GO terms include cell wall organization or biogenesis (GO: 0071554, P value = 1.16E-12), proton-transporting two-sector ATPase complex (GO: 0016759, P value = 1.21 × 10−3), and hexosyltransferase activity (GO: 0016758, P value = 1.128 × 10−4), all were crucial for secondary wall initiation (Fig. 3, Fig. S5). C5 comprises 2,037 genes that are highly expressed at 45 DPA, a stage associated with SCW thickening. These genes are primarily involved in xylan biosynthesis (GO: 0045492, P value = 6.58E-18), cell wall biogenesis (GO: 0042546, P value = 1.28E-12), and hemicellulose metabolic processes (GO: 0010410, P value = 3.22E-14), marking the peak of SCW deposition (Fig. 3, Fig. S6). Finally, C6 contains 2,534 genes specifically expressed in ovary tissues at 45 DPA, a developmental stage characterized by the absence of fiber formation. This cluster is enriched in GO terms of mRNA processing (GO: 0006397, P value = 2.12E-09), nuclear body (GO: 0016604, P value = 6.90E-07), and protein homodimerization activity (GO: 0042803, P value = 5.27E-05) (Fig. 3, Fig. S7), suggesting transcriptional regulation and developmental readiness of the ovary. Collectively, the functional enrichment and clustering analyses delineate the fiber developmental process in B. ceiba into five distinct stages: initiation, elongation and expansion, transition, secondary wall formation, and maturation. These findings underscore dynamic temporal changes in gene expression and provide novel insights into the genetic regulation of fiber cell development in B. ceiba.

Prediction of hub genes through time-ordered gene co-expression networks

To identify key genes and transcription factors (TFs) involved in fiber development from initiation to maturation, we constructed time-ordered gene co-expression networks using GENIE3 [22] based on above six expression clusters (C1-C6). During the initiation stage, C1 genes were mainly associated with major biological processes such as plant hormones signaling transduction, cell wall biosynthesis, Ca2+ signal transduction, and trichome branch initiation (Fig. 4a). These genes likely play pivotal roles in early cell division and trichome clusters formation. Notable examples include calmodulin-related calcium sensor CML42 (Calmodulin-like protein 42, BceiT011556.1), calmodulin (CaM, BceiT037931.1), KCBP-interacting Ca2+ binding protein (KIC, BceiT035984.1), IRREGULAR TRICHOME BRANCH 2 (ITB2, BceiT017465.1), SQUAMOSA promoter binding protein-like 13 (SPL13, BceiT005924.1), SCARECROW-like 11 (SCL11, BceiT002336.1), and Repressor of GA 2 (RGA2, BceiT036494.1). In C2, genes appeared to regulate the transition from inner pericarp protrusions to fiber cells, likely contributing to fiber initiation based on morphological observations. These genes were involved in auxin signaling, fiber cell proliferation, and microtubule (MT) organization, including auxin transport proteins PIN-FORMED1 (PIN1, BceiT011552.1/BceiT008741.1), Auxin response factor 6 (ARF6, BceiT007278.1), CYCLIN D3;2 (CYCD3;2, BceiT006923.1), 65-kDa microtubule-associated protein 1 (MAP65-1, BceiT030580.1), and IQ67-domain protein 26 (IQD26, BceiT005107.1) (Fig. S8). C3 genes were implicated in fiber elongation and expansion—key processes in fiber development—and included Gibberellic acid -stimulated Arabidopsis 7 (GASA7, BceiT007753.1), EXPANSIN5 (EXPA5, BceiT007210.1), EXORDIUM-LIKE 1 (EXL1, BceiT010690.1), Fiber protein 15 (Fb15, BceiT033427.1), and Xyloglucan endotransglucosylase 4 (XET4, BceiT028551.1), as illustrated in Fig. S9. In C4, genes were primarily associated with fiber cell transition, including roles in cell wall remodeling, lignin biosynthesis, and lipid metabolism. Representative genes include Xyloglucan endotransglucosylase/hydrolase 9 (XTH9, BceiT003667.1), PROTODERMAL FACTOR 1 (PDF1, BceiT017805.1), Peroxidase 52 (PRX52, BceiT024840.1), and GPI-anchored lipid transfer protein 1 (LTPG1, BceiT011022.1) (Fig. S10). During the final stages—fiber maturation and SCW thickening—genes involved in xylan and cellulose biosynthesis and assembly were upregulated. These included IRREGULAR XYLEM 15 LIKE (IRX15L, BceiT007783.1), IRREGULAR XYLEM 10/Glucuronoxylan glucuronosyltransferase 2 (IRX10/GUT2, BceiT006210.1), Cellulose synthase 4 (CESA4, BceiT020256.1), WLIM1 (BceiT000278.1), Chitinase-like 2 (CTL2, BceiT029770.1), and Fasciclin-Like Arabinogalactan-proteins 12 (FLA12, BceiT033294.1) (Fig. 4b). These results suggest that the unique physical properties of B. ceiba fibers are largely shaped by their specific lignin and hemicellulose composition.

[IMAGE OMITTED: SEE PDF]

Conservation and divergence of fiber development between B. ceiba and upland cotton

To uncover the molecular basis of the unique fiber characteristics of B. ceiba, we compared its RNA-seq profiles with RNA-seq data from developing fibers of upland cotton (TM-1). Although numerous key genes involved in fiber development have been characterized in upland cotton, we observed pronounced differences in expression dynamics between the two species across developmental stages, including initiation, elongation, and SCW formation (Fig. 5). The MYB-bHLH-WD40 (MBW) regulatory complex, which governs trichome initiation in Arabidopsis and upland cotton [31], exhibited stronger expression in B. ceiba relative to TM-1. Specifically, MYB77 (BceiT000249.1), MYC2 (BceiT012396.1), and TRANSPARENT TESTA GLABRA1 (TTG1, BceiT009830.1) reached peak expression during the initiation stage in both B. ceiba and TM-1 (Fig. 5), supporting a conserved role in specifying epidermal cell fate in the fruit inner epidermis (B. ceiba) and the seed outer epidermis (upland cotton) [32, 33]. Genes associated with cellulose and hemicellulose biosynthesis—such as CESA3 (BceiT029697.1), xyloglucan xylosyltransferase 1 (XXT1, BceiT022988.1), and UDP-glucose dehydrogenase 2 (UGD2, BceiT027249.1)—were also exhibited higher expression levels in B. ceiba than in TM-1 during the initiation stage (Table S4), which may contribute to the higher hemicellulose content of B. ceiba fiber. Furthermore, genes implicated in trichome branching and cell division, including GENERAL CONTROL NON-REPRESSED PROTEIN 5 (GCN5, BceiT010663.1), ORESARA 15 (ORE15, BceiT012202.1), CML42 (BceiT011556.1), ITB2 (BceiT017465.1), and SPL14 (BceiT018111.1), showed high expression in B. ceiba fibers during initiation (Fig. 4a, Table S4). Similarly, central regulators of the gibberellin (GA) pathway, such as SCL11 (BceiT002336.1), DELLA protein RGA1 (BceiT028859.1), and RGA-LIKE1 (RGL1, BceiT013627.1), were strongly expressed in this stage (Fig. 4a). In contrast, during the elongation and expansion phases, several genes crucial for cell wall loosening and actin cytoskeleton organization—such as EXPA5 (BceiT007210.1), XET4 (BceiT028551.1), PDF1 (BceiT017805.1), profilin 3 (PFN3, BceiT002731.1), and TEOSINTE BRANCHED1/CYCLOIDEA/PROLIFERATING CELL FACTOR 14 (TCP14, BceiT013201.1)—were expressed at relative lower levels in B. ceiba compared with TM-1 (Table S4). However, genes regulating cortical microtubules (MTs) orientation and cell wall remodeling, including Fragile Fiber 1a (FRA1a, BceiT007785.1), XTH7 (BceiT015982.1), and α-xylosidase 7 (XYL7, BceiT014204.1), displayed conserved expression patterns in both species. Notably, genes involved in lignin biosynthesis—such as Phenylalanine Ammonia Lyase 2 (PAL2, BceiT009127.1), 4-coumarate-CoA ligase 1 (4CL1, BceiT031608.1), Cinnamoyl CoA Reductase 1 (CCR1, BceiT025580.1), Flavonoid 3-Hydroxylase (F3H, BceiT016224.1), Cinnamyl Alcohol Dehydrogenase 5 (CAD5, BceiT014744.1), and PRX52 (BceiT024840.1)—were consistently upregulated in B. ceiba compared to TM-1 (Fig. S10 and Table S4). During SCW formation, genes related to lignin and hemicelluloses (e.g., xylans) biosynthesis also showed elevated expression in B. ceiba, including IRX15 (BceiT006510.1), IRX15L (BceiT007783.1), PARVUS (BceiT038158.1), and FRA8 (BceiT033111.1) (Fig. 5a). These transcriptional shifts may alter fiber strength and crystallinity by limiting cellulose microfibril aggregation (Fig. 4b and Table S4). Conversely, several genes central to cellulose biosynthesis and microtubule organization were expressed at lower levels in B. ceiba fibers compared to TM-1. These included CESA4 (BceiT020256.1), KORRIGAN (KOR, BceiT028094.1), WLIM1a (BceiT000067.1), and FLA12 (BceiT033294.1) [34,35,36,37], as well as MT-associated genes such as FRA1b (BceiT004374.1), IQD10 (BceiT030843.1), and WAVE DAMPENED2-LIKE 1 (WDL1, BceiT023894.1/BceiT029987.1) (Fig. 5b, Table S4).

[IMAGE OMITTED: SEE PDF]

To validate our results, we reanalyzed raw transcriptome data from a previous study (PRJNA908266), which included endocarp, ovule, and fiber tissues collected from 8.05 to 16.00 cm in length [30]. The results were consistent with our findings (Table S5), confirming the reliability of candidate gene identification and their expression patterns. In summary, we constructed a comprehensive regulatory network spanning four key stages of B. ceiba fiber development, integrating hub TFs and key enzyme-encoding genes (Fig. 4b). This network highlights both conserved and divergent regulatory mechanisms between B. ceiba and upland cotton, providing a molecular explanation for their distinct fiber properties (Fig. 6).

[IMAGE OMITTED: SEE PDF]

Validation of key gene expression profilesvia qRT-PCR

To verify the accuracy of the transcriptomic data, the expression profiles of 12 candidate genes were examined across three developmental stages of B. ceiba fibers (20, 35, and 45 DPA), as well as in ovary tissues at 45 DPA (45 ODPA), using qRT-PCR analysis (Fig. 7, Table S6). Genes associated with cell division and trichome development—such as ARF6 (BceiT007278.1), XXT1 (BceiT022988.1), and SCL11 (BceiT002336.1)—displayed high expression at 20 or 35 DPA. This observed temporal pattern indicates that the genes function not only in fiber initiation but also in elongation and in facilitating the initiation-to-elongation transition (Fig. 6). In contrast, EXPA5 (BceiT007210.1), a gene involved in cell wall loosening, showed elevated transcript levels from 35 to 45 DPA, consistent with the period of active fiber elongation and expansion (Fig. 6). Additionally, genes associated with MTs orientation and cell wall remodeling—such as IQD13 (BceiT030380.1) and IQD21 (BceiT011789.1)—showed peak expression at 35 DPA, coinciding with the fiber transition phase, during which dynamic reorganization of cortical MTs and cellulose microfibril reorientation are believed to be key drivers of fiber shaping and cell wall architecture. Moreover, numerous genes implicated in SCW biosynthesis demonstrated markedly increased expression from 35 to 45 DPA, consistent with the increased accumulation of cellulose, hemicellulose, and lignin during this period. Overall, the qRT-PCR results closely mirrored the RNA-seq data, confirming the reliability and robustness of our transcriptomic analysis.

[IMAGE OMITTED: SEE PDF]

Discussion

Natural plant fibers have long served as essential raw materials in the textile and materials processing industries. B. ceiba produces a typical cellulosic fiber, yet its distinct features—such as a thin cell wall, reduced mechanical strength, and poor spinnability—set it apart from conventional commercial fibers like those of upland cotton [27]. Despite its industrial relevance, the molecular mechanisms governing fiber development and cell morphogenesis in B. ceiba remain largely unexplored.

To elucidate the morphological and functional divergence between B. ceiba and upland cotton fibers, we propose a hypothetical model outlining the key regulatory factors involved in the formation of B. ceiba fibers. Unlike the well-characterized fiber initiation in cotton, where single trichome cells expand directly into fibers, B. ceiba fibers arise from clusters of trichome cells generated via successive cell divisions during the initiation phase [26, 27]. This developmental distinction is mirrored at the transcriptomic level, where genes associated with cell division or trichome cell proliferation—including CML42 (BceiT011556.1), histone acetyltransferase GCN5 (BceiT010663.1), SPL14 (BceiT018111.1), and PLATZ transcription factor ORE15 (BceiT012202.1)—exhibited elevated expression in B. ceiba than in upland cotton during the initiation stage, such as the the calmodulin-related calcium sensor [38, 39]. Nevertheless, the precise regulatory pathways and gene interactions involved in these atypical cell division events remain to be elucidated. The hollow morphology and thin cell walls of B. ceiba fibers may reflect enhanced activity of genes related to cell wall loosening and expansion. Notably, genes such as EXPA5 (BceiT007210.1), XET4 (BceiT028551.1), and Fb15 (BceiT033427.1) were more highly expressed in B. ceiba than in upland cotton. Conversely, genes integral to cellulose biosynthesis and secondary wall assembly—such as CESA8 (BceiT031769.1), KOR (BceiT028094.1), IRX9 (BceiT001194.1), and FLA12 (BceiT033294.1)—were expressed at lower levels in B. ceiba than in upland cotton. These findings indicate a trade-off between fiber extensibility and cell wall rigidity, which likely contributes to the mechanical characteristics of B. ceiba fibers.

Another defining feature of B. ceiba fibers is the absence of natural twist, which contributes to their lower strength and spinnability relative to upland cotton fibers, making it difficult to spin into yarns [40]. Previous studies have implicated cellulose synthase complexes (CSC, particularly CESA8) in fiber twist formation [13]. Moreover, xylan—a hemicellulose component enriched in B. ceiba fibers—may hinder microfibril aggregation by forming hydrogen-bonded coatings on cellulose, thereby reducing crystallinity and potentially preventing helical twisting [41, 42]. Supporting this hypothesis, xylan biosynthetic genes such as IRX15 (BceiT006510.1, BceiT035498.1), IRX15L (BceiT007783.1, BceiT016061.1), and PARVUS (BceiT038158.1), showed higher expression levels in B. ceiba than in upland cotton fibers (Fig. 5a). Moreover, fiber twist is also influenced by cortical MTs organization, which directs the deposition of cellulose microfibrils [43]. Our analysis revealed lower expression of MT-associated genes—FRA1b (BceiT004374.1), IQD10 (BceiT030843.1), and WDL1 (BceiT023894.1)—in B. ceiba compared to upland cotton. This reduced expression may impair helical cell growth, contributing to the observed fiber morphology. Collectively, these results support our proposed model, in which cortical MT dynamics, cell wall composition, and specific transcriptional regulators coordinate the unique developmental trajectory and structural features of B. ceiba fibers (Fig. 6).

Conclusion

This study offers a comprehensive transcriptomic analysis of fiber development in B. ceiba, revealing key molecular and structural differences from upland cotton. The unique fiber characteristics of B. ceiba—such as enhanced flexibility, reduced strength, and lack of twist—appear to be linked to differential expression of genes associated with cell division, wall modification, and microtubule dynamics. In particular, the composition of secondary cell walls—marked by a distinct balance of hemicellulose and reduced cellulose crystallinity—plays a central role in defining these traits. The identification of regulatory genes and pathways unique to B. ceiba fiber formation provides a foundation for future genetic and biotechnological approaches aimed at improving fiber quality in other species. Furthermore, B. ceiba represents a promising alternative source of natural fiber with distinctive physicochemical properties, potentially suitable for niche textile and industrial applications. Future research should prioritize optimizing the cultivation, harvesting, and post-processing of B. ceiba fibers to enhance their economic and industrial viability.

Data availability

Raw sequencing reads of transcriptome data have been deposited at the Genome Sequence Archive (GSA) database in the National Genomics Data Center of China National Center for Bioinformation/Beijing Institute of Genomics under the project number PRJCA034554 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA034554).

Abbreviations

ADF1:

Actin depolymerizing factor 1

ARF6:

Auxin response factor 6

B. ceiba :

Bombax ceiba

bHLH:

Basic helix-loop-helix

BLASTP:

Basic local alignment search tool for proteins

CAD:

Cinnamyl alcohol dehydrogenase

CaM:

Calmodulin

CCR2:

Cinnamoyl-CoA reductase 2

CESA:

Cellulose synthase A

CSC:

Cellulose synthase complex

CTL2:

Chitinase-like 2

CYCD3;2:

Cyclin D3;2

DEG:

Differentially expressed gene

DPA:

Days post anthesis

EXL:

Exordium-like

EXO:

Exordium

EXPA5:

Expansin A5

F3H:

Flavanone 3-hydroxylase

F-actin:

Actin filaments

FDR:

False discovery rate

FLA12:

Fasciclin-like arabinogalactan protein 12

FPKM:

Fragments per kilobase of exon model per million mapped fragments

FRA:

Fragile fiber

FSN1:

Fiber secondary cell wall-related NAC1

GASA7:

Gibberellic acid-stimulated Arabidopsis 7

GCN5:

General control non-derepressible 5

GO:

Gene ontology

GSA:

Genome sequence archive

GUT2:

Glucuronoxylan glucuronosyltransferase 2

IQD:

IQ67 domain protein

IRX :

Irregular xylem

IRX15L :

Irregular xylem 15-like

ITB2:

Irregular trichome branch2

KIC :

KCBP-interacting Ca2+ binding protein

KNL1:

Knotted1-like

KOR:

Korrigan

LTPG1:

GPI-anchored lipid transfer protein 1

MAP65-1:

microtubule-associated protein 65-1

MBW:

MYB-basic helix-loop-helix-WD40

MT:

Microtubule

NAC:

NAM, ATAF1/2 and CUC2 domain-containing transcription factor

ORE15:

Oresara 15

PAL2:

Phenylalanine ammonia lyase 2

PDF1:

Protodermal factor 1

PFN3:

Profilin 3

PIN1:

PIN-formed protein 1

PM:

Plasma membrane

PRX52:

Peroxidase 52

RGA2:

Repressor of gibberellin acid 2

RGL1:

RGA-like1

SCL11:

Scarecrow-like protein 11

SCW:

Secondary cell wall

SPL:

SQUAMOSA promoter binding protein-like

TCP14:

Teosinte branched1/cycloidea/proliferating cell factor 14

TF:

Transcription factor

TM-1:

Texas marker-1

TTG1:

Transparent testa glabra 1

UGD2:

UDP-Glucose dehydrogenase 2

WDL:

Wave dampened2-like

WLIM :

LIN-11, Isl1, and MEC-3 domain protein

XET4:

Xyloglucan endotransglucosylase 4

XTH:

Xyloglucan endotransglucosylase/hydrolase

XXT:

Xyloglucan xylosyltransferase

XYL7:

α-Xylosidase 7

References

Komati A, Anand A, Nagendla NK, Madhusudana K, Mudiam MKR, Babu KS, et al. Bombax ceiba calyx displays antihyperglycemic activity via improving insulin secretion and sensitivity: identification of bioactive phytometabolomes by UPLC-Qtof-MS/MS. J Food Sci. 2022;87(4):1865–81.

Namsa ND, Tag H, Mandal M, Kalita P, Das AK. An ethnobotanical study of traditional anti-inflammatory plants used by the Lohit community of Arunachal Pradesh, India. J Ethnopharmacol. 2009;125(2):234–45.

Wang LP, Xie SY, Jiang X, Xu CN, Wang YQ, Feng JF, et al. Therapeutic effects of Bombax ceiba flower aqueous extracts against loperamide-induced constipation in mice. Pharm Biol. 2023;61(1):125–34.

Zheng Y, Wang JT, Zhu YF, Wang AQ. Research and application of kapok fiber as an absorbing material: a mini review. J Environ Sci (China). 2015;27:21–32.

Baraniak J, Kania-Dobrowolska M. Multi-purpose utilization of Kapok fiber and properties of Ceiba pentandra tree in various branches of industry. J Nat Fibers. 2023;20(1):2192542. https://doi.org/10.1080/15440478.2023.2192542

Yaqoob A, Bashir S, Rao AQ, Shahid AA, Jenkins J. Transformation of α-EXPA1 gene leads to an improved fibre quality in Gossypium hirsutum. Plant Breed. 2020;139(6):1213–20.

Sartika D, Syamsu K, Warsiki E, Fahma F, Arnata IW. Nanocrystalline cellulose from Kapok fiber (Ceiba pentandra) and its reinforcement effect on alginate hydrogel bead. Starch – Stärke. 2021;73:2100033.

Hinchliffe DJ, Meredith WR, Delhorn CD, Thibodeaux DP, Fang DD. Elevated growing degree days influence transition stage timing during cotton fiber development resulting in increased fiber-bundle strength. Crop Sci. 2011;51:1683–92.

Huang G, Huang JQ, Chen XY, Zhu YX. Recent advances and future perspectives in cotton research. Annu Rev Plant Biol. 2021;72:437–62.

Wen XP, Chen ZW, Yang ZR, Wang MJ, Jin SX, Wang GD, et al. A comprehensive overview of cotton genomics, biotechnology and molecular biological studies. Sci China Life Sci. 2023;66(10):2214–56.

Wang NN, Li Y, Chen YH, Lu R, Zhou L, Wang Y, et al. Phosphorylation of WRKY16 by MPK3-1 is essential for its transcriptional activity during fiber initiation and elongation in cotton (Gossypium hirsutum). Plant Cell. 2021;33(8):2736–52.

Tian Y, Zhang TZ. MIXTAs and phytohormones orchestrate cotton fiber development. Curr Opin Plant Biol. 2021;59:101975.

Wen XP, Zhai YF, Zhang L, Chen YJ, Zhu ZY, Chen G, et al. Molecular studies of cellulose synthase supercomplex from cotton fiber reveal its unique biochemical properties. Sci China Life Sci. 2022;65(9):1776–93.

Wang Y, Li Y, Gong SY, Qin LX, Nie XY, Liu D, et al. Ghknl1 controls fiber elongation and secondary cell wall synthesis by repressing its downstream genes in cotton (Gossypium hirsutum). J Integr Plant Biol. 2022;64(1):39–55.

Zhang J, Huang GQ, Zou D, Yan JQ, Li Y, Hu S, et al. The cotton (Gossypium hirsutum) NAC transcription factor (FSN1) as a positive regulator participates in controlling secondary cell wall biosynthesis and modification of fibers. New Phytol. 2018;217(2):625–40.

Farooq T, Zain-ul-Hudda, Qadir F, Fatima R, Muneer R, Ali A, Shehzad S, Imran A, Sohail S, Iftikhar A, Mazhar M, Ghaffor A, Zaib AA, Anwar Z, Ijaz A, Misbah F, Rasool G, Firdous H, Razzaq A. Overexpression of gene ΔGH_A07G1537 associated with fiber quality in upland cotton (Gossypium hirsutum L.) Through pollen tube transformation method. Agrobiological Records. 2024;17:129–37.

Chen SF, Zhou YQ, Chen YR, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884-90.

Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15.

Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. Stringtie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290–5.

Frazee AC, Pertea G, Jaffe AE, Langmead B, Salzberg SL, Leek JT. Ballgown bridges the gap between transcriptome assembly and expression analysis. Nat Biotechnol. 2015;33(3):243–6.

Wu TZ, Hu EQ, Xu SB, Chen MJ, Guo PF, Dai ZH, et al. Clusterprofiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Camb). 2021;2(3):100141.

Huynh-Thu VA, Geurts P. dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data. Sci Rep. 2018;8(1):3384.

Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

McGinnis S, Madden TL. BLAST: at the core of a powerful and diverse set of sequence analysis tools. Nucleic Acids Res. 2004;32:20–5.

Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25(4):402–8.

Ding YH, Gao W, Qin Y, Li XP, Zhang ZN, Lai WJ, et al. Single-cell RNA landscape of the special fiber initiation process in Bombax ceiba. Plant Commun. 2023;4(5):100554.

Zhao GJ, Ge L, Ma HC, Yang JJ, Huang D, Ping P. Kapok capsule formation and fiber development process in Yuanjiang dry-hot valley. Chin J Appl Ecol. 2014;25(12):3443–50.

Gao J, Liu YJ, Zhao T, Chen JB. Composition, Structure and Property Analysis of Calotropis Gigantea, Kapok and Cotton Fibers. J Donghua Univ (Natural Science). 2012;38(2):151-5.

Gómez-Maqueo X, Gamboa-deBuen A. The biology of the genus Ceiba, a potential source for sustainable production of natural fiber. Plants. 2022. https://doi.org/10.3390/plants11040521.

Shao L, Jin SK, Chen JW, Yang GS, Fan R, Zhang ZY, et al. High-quality genomes of Bombax ceiba and Ceiba pentandra provide insights into the evolution of Malvaceae species and differences in their natural fiber development. Plant Commun. 2024;5(5):100832.

Ramsay NA, Glover BJ. MYB-bHLH-WD40 protein complex and the evolution of cellular diversity. Trends Plant Sci. 2005;10(2):63–70.

Wen XP, Chen ZW, Yang ZR, Wang MJ, Jin SX, Wang GD, et al. A comprehensive overview of cotton genomics, biotechnology and molecular biological studies. Sci China Life Sci. 2023 ;66(10): 2214-56. https://doi.org/10.1007/s11427-022-2278-0

Timilsina R, Kim Y, Park S, Park H, Park SJ, Kim JH, et al. Oresara 15, a PLATZ transcription factor, controls root meristem size through auxin and cytokinin signalling-related pathways. J Exp Bot. 2022;73(8):2511–24.

Han LB, Li YB, Wang HY, Wu XM, Li CL, Luo M, et al. The dual functions of WLIM1a in cell elongation and secondary wall formation in developing cotton fibers. Plant Cell. 2013;25(11):4421–38.

Maloney VJ, Samuels AL, Mansfield SD. The endo-1,4-β-glucanase Korrigan exhibits functional conservation between gymnosperms and angiosperms and is required for proper cell wall formation in gymnosperms. New Phytol. 2012;193(4):1076–87.

Ma Y, Shafee T, Mudiyanselage AM, Ratcliffe J, MacMillan CP, Mansfield SD, et al. Distinct functions of FASCILIN-LIKE ARABINOGALACTAN PROTEINS relate to domain structure. Plant Physiol. 2023;192(1):119–32.

Kumar M, Mishra L, Carr P, Pilling M, Gardner P, Mansfield SD, et al. Exploiting cellulose synthase (CESA) class specificity to probe cellulose microfibril biosynthesis. Plant Physiol. 2018;177(1):151–67.

Yang BJ, Sun YB, Minne M, Ge YH, Yue QR, Goossens V, Mor E, Callebaut B, Bevernaege K, Winne JM, Audenaert D, Rybel BD. SPL13 controls a root apical meristem phase change by triggering oriented cell divisions. Science. 2024;386(6723): eado4298. https://doi.org/10.1126/science.ado4298.

Wang TY, Jia QM, Wang W, Hussain S, Ahmed S, Adnan, Zhou DX, Ni ZF, Wang SC. GCN5 modulates trichome initiation in Arabidopsis by manipulating histone acetylation of core trichome initiation regulator genes. Plant Cell Rep. 2019;38(6):755–65.

Kaur J, Sawhney RS, Singh H, Sandhu M, Sohal RS, Singh A, et al. Comparative analysis of triboelectric behavior of natural waste biomaterials and green energy harvesting. Energy Environ. 2023 ;36(3):1–12. https://doi.org/10.1177/0958305X231193873.

Bringmann M, Landrein B, Schudoma C, Hamant O, Hauser MT, Persson S. Cracking the elusive alignment hypothesis: the microtubule-cellulose synthase nexus unraveled. Trends Plant Sci. 2012;17(11):666–74.

Berglund J, Mikkelsen D, Flanagan BM, Dhital S, Gaunitz S, Henriksson G, et al. Wood hemicelluloses exert distinct biomechanical contributions to cellulose fibrillar networks. Nat Commun. 2020;11(1):4692.

He JB, Zhao XH, Du PZ, Zeng W, Beahan CT, Wang YQ, et al. KNAT7 positively regulates xylan biosynthesis by directly activating IRX9 expression in Arabidopsis. J Integr Plant Biol. 2018;60(6):514–28.