Five-year-old healthy H. rhamnoides cv. ‘Shenqiuhong’ plants grown in Xinjiang, China, were used in this study. Five independent plants were selected as biological replicates.

Axillary buds were collected at three morphologically distinct developmental stages: T0 (dormant stage), T1 (budburst stage, 1–2 mm in length), and T2 (branch elongation stage, 5–10 mm in length), following established protocols (Ni et al., 2017; Yao et al., 2022; Tang et al., 2023). Axillary buds at different developmental stages were sampled from the same set of five plants, constituting a paired sampling design across stages.

For each plant at each developmental stage, axillary buds from multiple branches were pooled to generate one biological replicate (n = 5). Bud sampling and tissue separation were conducted based on external morphological characteristics, including bud size, shape, and visible developmental status. During dissection, only a minimal amount of immediately adjacent nodal tissue was retained when unavoidable, particularly at early developmental stages when axillary buds were small and closely associated with the stem, in order to minimize the inclusion of elongated internode or mature stem tissues.

All samples were harvested between 09:00 and 11:00 to minimize diurnal effects, immediately flash-frozen in liquid nitrogen, and stored at –80 °C until RNA extraction.

To evaluate transcriptional responses to auxin, a decapitation and auxin replacement assay was performed. Plants were assigned to three treatments: (i) Control: apical bud retained; (ii) Decapitation: apical bud removed; (iii) Decapitation + IAA: 10 μM IAA applied immediately to the cut surface. The first pair of axillary buds below the apex, along with adjacent stem tissues, was collected at 0, 3, 6, 12, 24, and 48 hours after treatment. Three biological replicates were included for each treatment–time combination, yielding a total of 54 samples. Sample handling followed the same procedure as Section 2.1.1 (Prusinkiewicz et al., 2009; Balla et al., 2016; Beveridge et al., 2000; Shimizu-Sato et al., 2009).

Total RNA was extracted using the FastPure Plant Total RNA Isolation Kit (Vazyme, China). RNA purity and concentration were determined with a NanoDrop 2000 spectrophotometer, and integrity was assessed using an Agilent 2100 Bioanalyzer. Samples with RIN ≥ 7.0 and A260/A280 ratios of 1.8–2.2 were used for library preparation.

Strand-specific libraries were constructed from poly(A)+ mRNA enriched using oligo(dT) magnetic beads, including fragmentation, cDNA synthesis, end repair, A-tailing, adapter ligation, and PCR amplification. Library quality was verified using Qubit and qPCR. Qualified libraries were sequenced on an Illumina platform to generate 150 bp paired-end reads.

Raw reads were processed using DataFilterQC (v0.0.2; parameters: −x 5, −l 30, −n 3) to remove adapters and low-quality sequences. Clean reads were aligned to the Hippophae rhamnoides L. reference genome, together with the corresponding gene structure annotation files, downloaded from the China National GeneBank DataBase (CNGBdb; https://db.cngb.org/), using HISAT2 (v2.2.1; parameters: –dta-cufflinks –no-unal –un-conc-gz) (Kim et al., 2019). BAM files were sorted, merged, and indexed using SAMtools (v0.1.19) (Li et al., 2009).

Gene-level read counts were obtained using featureCounts (v2.0.2) (Liao et al., 2014). Differential expression analysis was performed with edgeR (v3.3.3) using TMM normalization and a negative binomial GLM (Robinson et al., 2010). Genes with |log₂FC| ≥ 1 and FDR < 0.05 were considered differentially expressed; for comparisons with fewer than 100 differentially expressed genes (DEGs), thresholds were relaxed to |log₂FC| ≥ 1 and p < 0.05. Transcriptome assembly and quantification were conducted using StringTie (Pertea et al., 2015) and GFF Utilities (Pertea and Pertea, 2020).

Gene models derived from the H. rhamnoides genome annotation were standardized and assigned unique identifiers using AGAT (Another Gff Analysis Toolkit), resulting in AGAT-based gene IDs (e.g., agat-gene-22040) used throughout this study.

Functional annotation was performed against NR, Swiss-Prot, and KEGG databases. GO and KEGG enrichment analyses were conducted using topGO (Alexa and Rahnenfuhrer, 2023) and KOBAS (Bu et al., 2021), with hypergeometric tests and Benjamini–Hochberg correction (FDR < 0.05). Transcription factors were identified using iTAK (Zheng et al., 2016), PlantTFDB, and HMMER (Potter et al., 2018).

qRT-PCR was performed on RNA from axillary buds at T0, T1, and T2 (n = 5) (Tang et al., 2023; Wang et al., 2023; Xie et al., 2024). cDNA synthesis used HiScript II Q RT SuperMix (+gDNA wiper). Reactions were run on a QuantStudio system with ChamQ Universal SYBR qPCR Master Mix. Primers annealed at 54°C, and 18S rRNA was the reference. Relative expression was calculated via 2^-ΔΔCt. One-way ANOVA with Tukey’s post-hoc test (p < 0.05) assessed significance; three technical replicates were performed.

FPKM values for the six candidate genes were extracted from the IAA treatment time-course (0–48 h), log₂(FPKM + 1)-transformed, and visualized in R using ggplot2. Temporal expression profiles revealed rapid, stage-specific responses, highlighting their roles as key mediators in auxin signaling during axillary bud development (Wójcikowska and Gaj, 2017; Shen et al., 2014; Markakis et al., 2013; Ren et al., 2015).

To systematically characterize transcriptional changes during axillary bud development in H. rhamnoides, we performed comparative transcriptome analyses across three key developmental stages (T0, T1, and T2). Sampling was conducted as illustrated in Figure 1, with five biological replicates per stage, each comprising pooled axillary buds from three branches of the same plant. DEGs were identified using thresholds of |log₂FC| ≥ 1 and FDR < 0.05.

A total of 5,452 DEGs were identified between T0 and T1, accounting for approximately 49% of all expressed genes, indicating extensive transcriptional reprogramming during initial bud activation. Notably, the comparison between T2 and T0 also revealed the same total number of DEGs (5,452); however, the composition of upregulated and downregulated genes differed, indicating continued but reorganized transcriptional activity during branch elongation. In contrast, only 213 DEGs (67 upregulated, 146 downregulated) were detected between T1 and T2, suggesting transcriptional stabilization during subsequent branch elongation.

Stage-specific transcriptional shifts and high reproducibility among replicates were visualized through hierarchical clustering heatmaps and volcano plots (Figures 4A, C–E). The heatmap illustrated global expression patterns for the T1 vs. T0 comparison (Figure 4A), while volcano plots highlighted statistically significant DEGs across comparisons (Figures 4C–E). Detailed GO and KEGG enrichment results for pairwise comparisons are provided in the Supplementary Figures, with Supplementary Figure S1 showing T0 vs. T1, Supplementary Figure S2 showing T2 vs. T0, and Supplementary Figure S3 showing T1 vs. T2. Venn diagram analysis further revealed 4,409 DEGs shared between T1 vs. T0 and T2 vs. T0, with 1,043 unique to each comparison, underscoring the coexistence of common and stage-specific transcriptional programs (Figure 4B). A complete list of differentially expressed genes from all pairwise comparisons is provided in Supplementary Table S2. These DEGs established a foundation for subsequent time-series clustering using Mfuzz and hub gene identification via WGCNA, enabling the capture of dynamic expression patterns and the identification of core regulators governing stage-specific developmental transitions.

To interpret the biological significance of these transcriptional changes, we performed GO and KEGG enrichment analyses for all identified DEGs.

In the T1 vs. T0 comparison, upregulated genes were significantly enriched in biological processes such as xylem development, organ morphogenesis, auxin response, and microtubule motor activity. The enrichment of xylem development-related pathways at the T0 and T1 stages is consistent with early axillary bud activation, during which vascular reconnection between the bud and the main stem is progressively established. This enrichment may reflect vascular-related developmental reprogramming during early bud activation, rather than necessarily indicating stem tissue contamination. KEGG analysis further indicated their involvement in plant hormone signal transduction and zeatin biosynthesis, suggesting a pronounced activation of hormonal signaling and cellular restructuring during the transition from dormancy to bud burst. In contrast, downregulated genes were associated with defense response, secondary metabolism, and oxidation–reduction processes, showing enrichment in phenylpropanoid biosynthesis, terpenoid biosynthesis, and the MAPK signaling pathway. These patterns reflect a metabolic shift from stress defense toward growth-oriented processes, highlighting the roles of hormonal reactivation and energy metabolism in bud release from dormancy.

For the T2 vs. T0 comparison, upregulated genes were primarily related to microtubule–kinetochore organization, chromatin assembly, and auxin transport. Corresponding KEGG pathways included plant hormone signal transduction, phenylpropanoid biosynthesis, and flavonoid biosynthesis, indicating active cell proliferation and differentiation during branch formation. Downregulated genes were mainly involved in plasma membrane composition, cell wall modification, and hormone responses, with significant enrichment in plant–pathogen interaction, MAPK signaling, and zeatin biosynthesis. It is noteworthy that the total number of DEGs in the T2 vs. T0 comparison (5,452) was identical to that in the T1 vs. T0 comparison, underscoring extensive transcriptional reprogramming across early and late developmental stages.

In the T1 vs. T2 comparison, where only a limited number of DEGs were identified (213), upregulated genes were enriched in abscisic acid (ABA) metabolism, osmotic stress response, and transcription factor activity. Downregulated genes were associated with auxin signaling, phospholipase C activity, and plasma membrane function. Key KEGG terms included plant hormone signal transduction, ABC transporters, and TCA cycle, implying a shift toward energy supply and environmental adaptation as buds progress from activation to elongation, with overall gene expression stabilizing during this phase.

Collectively, these enrichment results reveal stage-specific regulation of hormone signaling, metabolic activity, and cell proliferation during axillary bud development in H. rhamnoides. The initial bud activation stage (T0 → T1) is characterized by enhanced auxin- and cytokinin-related signaling and changes in energy- and growth-associated pathways, which are consistent with transcriptional reprogramming during bud outgrowth.

Detailed enrichment results for each comparison are provided in Supplementary Figures S1. These visualizations support stage-specific transcriptional dynamics and functional interpretations.

To investigate the temporal expression dynamics during axillary bud development in H. rhamnoides, we performed stage-resolved temporal clustering using the Mfuzz algorithm. This approach identified genes sharing similar expression trajectories, capturing coordinated regulatory programs underlying developmental transitions. Based on standardized expression profiles, genes were grouped into eight distinct clusters (Clusters 1–8) as shown in Figure 5 and detailed in Supplementary Table S3.

The gene counts, expression patterns, and cluster classifications are summarized in Supplementary Table S3 and visually represented in Figures 5I, J. Three major temporal expression modes were prominently represented: Cluster 8 (2,003 genes, ~18% of DEGs) exhibited an early-activation pattern, with expression rapidly increasing from T0 to T1, suggesting a pivotal role in bud initiation. Cluster 6 (3,294 genes, ~30%) displayed sustained activation, with gradual upregulation from T0 to T2, consistent with branch elongation and active cell proliferation. By contrast, Cluster 4 (2,677 genes, ~24%) showed a suppression–release pattern, characterized by transient downregulation from T0 to T1 followed by partial recovery, indicative of stress adaptation and energy metabolism in later developmental stages.

In summary, Mfuzz clustering effectively delineated temporal expression trajectories and revealed that distinct gene groups dominate at different developmental phases. These clusters capture dynamic transcriptional changes across developmental stages, providing a foundation for subsequent identification of hub genes and regulatory modules.

To interpret the biological significance of these transcriptional changes, GO and KEGG enrichment analyses for all identified temporal clusters are shown in Supplementary Figure S4. The complete gene lists for each temporal expression cluster are provided in Supplementary Table S3.

The early-activation module (Cluster 8) was significantly enriched in genes involved in auxin signaling, transcriptional regulation, and cell fate determination. Expression peaked at T1 after a sharp increase from T0, followed by a slight decline at T2 (T0 < T1 > T2). Most candidate branching-related genes were concentrated in this cluster, underscoring its central role in initiating axillary bud outgrowth. Notably, the BRC1 homolog in H. rhamnoides was included in Cluster 8, exhibiting a rapid increase in expression from T0 to T1 during budburst, followed by a slight decrease at T2. This temporal pattern mirrors other early-activation genes, including key auxin-responsive regulators such as ARF, IAA16, and SAUR36, suggesting that BRC1 may act as an integrator of hormonal and transcriptional signals during early bud activation (Aguilar-Martínez et al., 2007; González-Grandío et al., 2013).

The sustained-activation module (Cluster 6) comprised genes related to cell division, microtubule organization, mitotic spindle formation, and ribosome biogenesis. Expression levels increased progressively from T0 to T2 (T0 < T1 < T2), supporting continuous branch growth and meeting the biosynthetic demands of developing buds.

In contrast, the suppression–release module (Cluster 4) was enriched in genes associated with reactive oxygen species (ROS) detoxification, ABA signaling, and mitochondrial energy metabolism. Expression decreased sharply from T0 to T1 and remained low through T2 (T0 > T1 ≈ T2), reflecting a transient suppression phase followed by metabolic adaptation during later stages of bud development.

Collectively, these results reveal a temporally coordinated transcriptional program in which hormonal signaling, transcriptional regulation, energy metabolism, and stress responses are orchestrated to regulate axillary bud outgrowth and branch formation. The inclusion of BRC1 and other early auxin-responsive genes in Cluster 8 highlights a hub network potentially governing dormancy release and the initiation of branching. The correspondence between Mfuzz cluster patterns and functional annotations provides a robust framework for subsequent identification of hub genes and regulatory modules.

Although Cluster 8 represented the most prominent early-activation pattern, a small subset of genes in Cluster 2 also showed transient upregulation at T1, but did not contribute substantially to the final set of core branching-related candidates.

To complement the time-series clustering analysis and further characterize coordinated gene expression patterns during axillary bud development, weighted gene co-expression network analysis (WGCNA) was performed using normalized FPKM values of all expressed genes. This approach clusters genes based on expression similarity across samples and facilitates the identification of co-expression modules with shared expression dynamics.

Genes with low expression variance (variance < 0.1) were excluded prior to network construction. Hierarchical clustering of samples revealed no apparent outliers, indicating that all samples were suitable for WGCNA (Figure 6A). A soft-thresholding power of 6 was selected to approximate a scale-free network topology (scale-free R² ≈ 0.9) while maintaining sufficient mean connectivity (Figure 6B). Based on the topological overlap matrix (TOM), hierarchical clustering combined with dynamic tree cutting identified multiple co-expression modules. After merging highly similar modules, four major modules—blue, brown, green, and turquoise—were retained for subsequent analyses (Figure 6C). Genes within each module exhibited highly correlated expression profiles, suggesting coordinated transcriptional behavior. The gene composition of each merged WGCNA module is provided in Supplementary Table S4.

To examine the relationship between co-expression modules and axillary bud developmental stages, module eigengene expression patterns were compared across T0, T1, and T2 (Figure 6E). The blue module showed relatively higher eigengene expression at T1 compared with T0 and T2, indicating an association with transcriptional changes occurring during bud activation. The brown module displayed a gradual increase in eigengene expression from T0 to T2, consistent with sustained transcriptional activity during bud outgrowth and branch development. In contrast, the turquoise module exhibited elevated eigengene expression at T2, suggesting potential involvement in later developmental processes. The green module maintained relatively stable expression across all stages, implying roles in basic cellular functions.

Correlation analysis among module eigengenes revealed distinct relationships among the identified modules (Figure 6D). Some 5modules showed positive correlations, whereas others exhibited negative correlations, reflecting differences in their overall expression trends across developmental stages. Collectively, these results indicate that WGCNA effectively captured structured co-expression modules with clear stage-associated expression patterns.

Functional enrichment analysis was conducted to explore the biological relevance of the four major modules (Figure 6F). The blue module was enriched in transcription factors (e.g., MYB77-like and GATA18-like) and ribosomal protein-related genes, suggesting a role in transcriptional regulation and translational capacity during axillary bud activation. The brown module was primarily enriched in translation-related processes and protein biosynthesis, consistent with its progressive upregulation during sustained bud growth. The green module was associated with redox regulation and amino acid metabolism, indicating potential involvement in maintaining cellular metabolic homeostasis. The turquoise module was enriched in protein phosphorylation and kinase activity, and included highly connected hub genes such as serine/threonine-protein kinase BSK3 and Rac-like GTP-binding protein 5, highlighting its potential role in signaling processes during later stages of bud development.

Overall, WGCNA identified biologically coherent co-expression modules and candidate hub genes associated with different stages of axillary bud development, providing a network-level perspective that complements the temporal expression patterns revealed by Mfuzz clustering.

To validate the reliability of the RNA-seq data, we examined the expression patterns of six key candidate genes—ARF, IAA16, SAUR36, PP2C, MYC, and NAC41—across three developmental stages (T0, T1, T2) using qRT-PCR. The qRT-PCR results were highly consistent with the RNA-seq data (Figure 7), confirming the accuracy of the transcriptomic analysis. Both methods revealed significant upregulation of these genes during the transition from bud dormancy (T0) to activation (T1), followed by stabilization or a slight decline during branch elongation (T2). Although minor differences in absolute expression levels were observed, these variations fell within acceptable ranges and did not affect the biological interpretation.

This strong concordance verifies the robustness of the RNA-seq dataset and provides experimental support for the potential regulatory roles of these genes in axillary bud development. Primer sequences, amplification conditions, and detailed qRT-PCR procedures are provided in the Supplementary Materials, Primer sequences and raw qRT-PCR data are provided in Supplementary Table S6.

To further investigate the roles of these genes in auxin-mediated bud outgrowth, we analyzed their expression dynamics following apical decapitation and exogenous IAA application (Figure 8). Decapitation significantly induced the expression of all six hub genes (ARF, IAA16, SAUR36, PP2C, MYC, and NAC41) within 3–48 hours, effectively releasing their transcriptional repression. This response pattern closely mirrored their activation during the natural T0-to-T1 transition, indicating that decapitation transcriptionally mimics the natural early activation process of axillary buds.

Importantly, exogenous IAA treatment partially or completely reversed the decapitation-induced upregulation, restoring gene expression to near-control levels. This “induction–rescue” response is consistent with auxin-mediated regulation of apical dominance and supports an association between apical auxin status and the transcriptional behavior of these candidate genes. Furthermore, distinct temporal response profiles were observed among genes: ARF, IAA16, and SAUR36 showed rapid and transient induction, characteristic of early auxin-responsive elements, whereas PP2C and NAC41 displayed delayed and sustained upregulation, suggesting roles in downstream signaling or hormonal crosstalk.

In summary, expression analyses under both natural development and hormone treatment conditions not only validate the transcriptomic findings but also identify a core set of auxin-responsive regulators crucial for axillary bud release in H. rhamnoides. These results establish a mechanistic framework for understanding the molecular regulation of branching architecture in woody plants.