The plant material utilized for the study was whole floral buds, collected from a 50-year-old tree of Rhododendron delavayi Franch. in the Kunming Golden Temple area under normal/natural conditions. The Temple is situated on the Mingfeng Mountain in Kunming city (102°43′ 5.99″E, 25°02′ 20.00″N) of Yunnan province, Southwest China. The samples were collected in three biological replicates from different branches of the tree at five stages of bud dormancy. The diameter and length of eight sampled floral buds were measured to check their growth rates of progression. The floral buds slowed down their growth starting form late August and entered into the phase of pre-dormancy at the end of September. Then they started to cease the growth process, entering the para-dormancy phase (Dormancy Initiation-DI: September 30th to October 29th), which continues till the start of December with a complete cessation in growth and development, so they achieved the endo-dormancy phase (dormancy maintenance-DM: October 30th to December 2nd). Then, a slow and steady increase in bud size and growth could be observed, as per temperature changes, falling into the eco-dormancy stage (dormancy release initiation-RI: December 3rd to December 27th), which completely diminished in mid-January, giving way to the release phase (full release from dormancy-FR: December 28th to January 14th). The sampled floral buds were immediately kept in liquid nitrogen after being wrapped in aluminum foil and then stored at −80°C for further analyses, including physiological and hormonal content determination and RNA extraction for transcriptome sequencing. The histological assays were conducted by fixing and embedding the sampled dormant bud tissues in wax. This was done to harden the tissues for the easy cutting of buds into different longitudinal and transverse sections. After staining of these tissues, they were examined under a light microscope to distinguish different tissue components.

Estimating and quantifying hormones and physiological indices were carried out to mine for the regulatory mechanisms involved in the onset and release of the dormancy phenomenon. In this regard, four endogenous hormones, namely, ABA, zeatin, IAA, and GAs, were quantified from the sampled floral buds of different dormancy stages using the triple quadrupole LC-MS8040 platform following standard protocols (Xiao et al., 2020). Besides this, seven physiological indices were estimated, comprising soluble starch content (mg/g), soluble sugar content (%), soluble protein content (%), catalase (CAT) activity (U/g min), superoxide dismutase (SOD) activity (U/g.h), proline content (%), and malondialdehyde (MDA) content (mmol/g). The soluble sugar and soluble starch contents were determined by the anthrone colorimetry method (Dubois et al., 1956). The Coomassie Brilliant Blue (CBB) G-250 was utilized for the quantitation of soluble protein contents (Bradford, 1976). For the determination of SOD activity, nitro-blue tetrazolium (NBT) chemical was used via the light reduction method (Kono, 1978). CAT activity was estimated by using the Aebi protocol (Aebi, 1984). The Bates method (Bates et al., 1973) was used to determine the physiological status of plant entering dormancy by detecting the amount of free proline content. The MDA contents were estimated via the protocol mentioned earlier in studies (Wang et al., 2015; Zhang et al., 2021).

Fifteen libraries were constructed from triplicated floral buds randomly selected at the five dormancy stages for transcriptome sequencing. The total RNA extraction was performed by following the protocol of TRIzol (TaKaRa, Dalian, China), and the extracted RNA was treated with DNase-I, Oligo (dT) to obtain mRNA. The quality assessment of extracted RNA regarding RNA integrity and contamination was estimated with the help of the Agilent 2100 Bioanalyzer system (Agilent Technologies, CA, United States) and 1% agarose gel electrophoresis, respectively. For the construction of paired-end libraries, 3 μg of RNA extracted from each sample was utilized and put in fragmentation buffer, followed by the preparation of cDNA from mRNA. The quantification and qualification of the sample library were assessed with the Agilent 2100 Bioanalyzer and the ABI StepOnePlus Real-Time System.

Moreover, according to the instructions recommended by the manufacturer (Illumina^®), libraries were created by utilizing Kit. For RNA-Seq, the Illumina HiSeq technique was used. Furthermore, the paired-end libraries (read length: 150 bp) were sequenced using the Illumina HiSeq 2500 platform.

To attain clean reads, the reads containing short sequence length, i.e., ≤50 bp, adapters, poly-N, and low-quality (≥20% of bases with quality value ≤15) were eliminated via operating the FASTQ program (Cock et al., 2010). To expand the RNA-sequence data utility, contigs of triplicated samples from each dormancy phase were pooled together to assemble them to shape the non-redundant transcripts. The HISAT program v0.1.6-beta (Kim et al., 2015), a fast and sensitive spliced aligner, was utilized to align the transcripts parallel to the reference genome (Zhang L. et al., 2017).

The transcript containing absent features in the reference annotation is termed as “novel transcript” (Zhang L. et al., 2017). It can be a new transcript with totally unknown features or a new isoform of a previously known gene. For the reconstruction of transcripts, StringTie v1. 0. 4 software (Pertea et al., 2015) was utilized, which was further compared against reference annotations using CuffCompare v2. 2. 1. For prediction related to coding potential hidden in the novel transcripts, merging these coding sequences with reference transcripts was accomplished to attain a complete and comprehensive reference. This practice was carried out with the help of CPC v0. 9-r² software (Kong et al., 2007). The downstream analyses were based on this reference.

The number of clean reads obtained after integrating reference transcripts with those of novel coding transcripts mapped onto the genome was estimated using Bowtie2 (v2. 2. 5) (Langmead and Salzberg, 2012). The number of fragments per kilobase of exon model per million mapped reads (FPKM) expected for each unigene was calculated based on the gene’s read count and gene length. The gene expression levels were calculated with the RSEM package (v1. 2. 12) (Li and Dewey, 2011) from RNA-Seq data. The differentially expressed genes (DEGs) between samples from different growth/dormancy stages were identified using the NOIseq package (Tarazona et al., 2011). The significant threshold p-value was estimated based on the false discovery rate (FDR) method to determine the reliability of gene expression. Observing the two-fold difference between FPKM values and FDR ≤0.001, the identified genes were taken as significant DEGs. The Pearson’s correlation among different samples was calculated using the R “cor” function.

The functional annotation related to DEGs was attained via maximum BLAST hit against source databases, viz, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and NCBI Non-redundant (nr) protein sequences. These pathway enrichment analyses were performed with the “Phyper” function of R. To extract the open reading frame (ORF) of each DEG, “get off v. EMBOSS 6. 5.7.0” software (Rice et al., 2000) was used. The obtained ORFs were then aligned against transcription factor domains with the help of “hmmsearch v3.0” (Mistry et al., 2013). The hierarchical clustering of DEGs was performed using R’s “pheatmap” function.

To identify particular modules harboring highly correlated genes, the WGCNA package of R was implemented (Langfelder and Horvath, 2008). The adjacency matrix was generated by using normalized FPKM values. The associations between gene modules and phenotypic data were calculated by utilizing the default settings of the WGCNA package and importing phenotypic data. Further, the Topological Overlap Matrix (TOM) was constructed using the conversion option within the WGCNA package from the already generated adjacency matrix. Next to the expression network building step was the classification of transcripts with more similar patterns regarding their expression into a single module, followed by calculating Eigen/hub/key genes for constructed modules. The calculated eigengenes of modules were exported out via Cytoscape v 3. 8. 2 export (Shannon et al., 2003).

The RNA-seq data from Illumina sequencing was confirmed and validated by quantitative real-time PCR. The total RNA was taken out of floral bud tissues by following the total RNA Kit (BioTake, Beijing, China) manufacturer’s protocol, followed by its treatment with DNase-I. The cDNAs were synthesized from the reverse transcribed RNA with a cDNA Synthesis Kit (TaKaRa). The quantitative real-time PCR machine was used to accomplish the polymerase chain reaction (CFX96 Touch Real-Time PCR system, Biorad, Hercules, CA, United States). The qRT-PCR was accomplished by repeating each reaction three times and from three biological samples (Singh et al., 2021). The PCR was performed under the following conditions: 94°C (4 min) at the initial step, 94°C (15 s) for 40 cycles, 57°C (30 s), 70°C (30 s), 72°C (1 min) for the dissociation step, and a rise of 1°C at a regular interval of 5 s to shift from 75 to 90°C for acquiring the melting curve (Singh et al., 2020). Based on the gene annotation results, 19 candidate genes were selected from the identified DEGs and hub genes for quantitative RT-PCR analysis regarding their key roles in hormone signaling and regulatory pathways of floral growth, development, and bud dormancy (Supplementary Table S1). Specific primers were then designed and synthesized by Oligo v. 7. 0 software for the Rhododendron genes, and the elongation factor “SlEF” was utilized to perform as an internal control (Rotenberg et al., 2006). The specificity of primers was determined by a dissociative curve related to the individual gene. The 2^−ΔΔCt method was deployed to estimate the relative fold change of gene expressions (Livak and Schmittgen, 2001).

Based on bud morphology, five different growth cessation phases were identified, i.e., pre-dormancy, para-dormancy, endo-dormancy, eco-dormancy, and release from dormancy. The observed differences in floral buds’ size and growth rate before, during, and after the dormancy period were diagrammatically and graphically illustrated. The pictorial summary exhibited paraffin and freehand longitudinal sections of floral buds from September to the end of October. The diagrams showed that the floral buds had completed the differentiation among the inflorescence parts before entering the dormancy phase (Figure 2).

The meteorological details, including average temperature (°C), ranged daily and fortnightly during the experimental period, were presented in Figures 3A,B. The graph depicted the end of September as an initiation point for the size and growth rate decline. The cessation in growth started at the end of October, marked as dormancy initiation (DI). The dormancy maintenance (DM) phase continued till the start of December; then, a change in the size of the floral bud was observed from mid-to-end of December, shifting the bud growth rate towards the dormancy release initiation-RI-Phase. Finally, dormancy broke during the mid-to-end of January, bringing back the floral bud in its normal growth period (Figures 3C,D).

We estimated the content of four types of endogenous hormones, viz, ABA, IAA, cytokinins and/or zeatin, Gas, etc., in floral bud tissues during the mentioned bud dormancy phases. Moreover, seven physiological indices, i.e., soluble protein content, soluble starch content, soluble sugar content, CAT, SOD, proline content, and MDA, were also quantified during the bud dormancy phases. The ABA and zeatin hormones increased during dormancy phases, while IAA and GAs content lowered during dormancy phases from October to December (Figure 4). Similarly, an increase in soluble starch, soluble sugar content, soluble protein content, CAT, proline content, and MDA content was observed. However, SOD decreased during this period of dormancy (Figure 4).

The bud samples were collected in triplicate from the five dormancy stages. The number of raw reads obtained after RNA sequencing was about 51.55 million reads, filtered via removal of adaptors and ambiguous or low-quality reads. Consequently, approximately 49 million (95%) clean reads were obtained (Supplementary Table S2). On average, 4.43 Gb of clean data was obtained from each bud sample, with a Q20% greater than 99.07% and a Q30% greater than 96.85% (Supplementary Table S2). The clean base data contained GC content ranging between 46.8 and 48.37% (Supplementary Table S3). After assembling, almost 82.8% of the clean reads were aligned against the reference genome (Zhang L. et al., 2017), using the HISAT program (Kim et al., 2015). The analysis provided 65,531 transcripts, and 39,166 novel transcripts were detected (Supplementary Table S4).

The novel transcripts were then merged with reference transcripts to attain a complete reference from which the clean reads were mapped to calculate the expression of genes. The Venn diagram depicts the expressed transcripts identified uniquely at each dormancy phase (504, 514, 307, 240, and 253 uniquely expressed transcripts at pre-dormancy, para-dormancy, endo-dormancy, eco-dormancy, and dormancy release phases, respectively). Except for the dormancy release , a total of 23,537 expressed transcripts were found common at all the four dormancy phases (Figure 5). Before moving towards comparative transcriptome analyses, the data of each transcriptome was evaluated regarding quantitation and quality control. A correlation analysis was conducted for this purpose. Almost all the bud samples at different dormancy phases revealed a highly significant correlation among themselves (Supplementary Figure S1). Such quality of results paved the way further towards transcriptome profiling of genes related to the regulation of floral bud dormancy.

The gene expression profiles of sampled buds for growth and dormancy were estimated by the DEGseq algorithm to discover the differentially expressed genes (DEGs). A total of 13,542 DEGs were detected based on all possible pairwise comparisons of bud dormancy phases (Supplementary Table S5). Details of identified DEGs associated with following sets of comparison were, viz, pre-vs-para: 1878 (up: 24%, down: 76%), pre-vs-endo: 3667 (up: 23%, down: 77%), pre-vs-eco: 3356 (up: 33%, down: 67%), pre-vs-rele: 1408 (up: 41%, down: 59%), para-vs-endo: 391 (up: 40%, down: 60%), para-vs-eco: 641 (up: 63%, down: 37%), para-vs-rele: 680 (up: 70%, down: 30%), endo-vs-eco: 216 (up: 53%, down: 47%), endo-vs-rele: 805 (up: 83%, down: 17%):, and eco-vs-rele: 490 (up: 67%, down: 33%) (Figure 6 and Supplementary Table S6).

Comparative transcriptomics identified 405 DEGs conserved between different comparison groups of floral bud dormancy phases (Supplementary Table S6). The conserved genes were selected based on the expression conserved between at least two comparisons. Among the 405 genes, 4, 4, 29, 30, 14, 21, 2, and 300 were identified for floral bud dormancy, flowering initiation, IAA biosynthesis, embryogenesis, ABA biosynthesis, GA synthesis, zeatin, cytokinesis, and transcription regulation, respectively (Supplementary Table S7). The four floral bud dormancy transcript genes were PORCN, LEA29, TRMT2A, and ATRX.

We further explored the conserved genes based on KEGG annotation and identified four genes encoding FLOWERING LOCUS T (FLT) and elongation factor 1 alpha-like protein. All four genes related to floral initiation were down-regulated during all dormancy phases. A total of 29 genes were associated with auxins. Most genes like the SAUR family, RPL9e, and GH3 were related to IAA, depicting their role in sustaining the dormancy period (Supplementary Table S8). Similarly, comparative transcriptomics suggested genes related to embryogenesis, viz, GST, ATRX, SLC39A1, SERK1, PGM, LEA29. They were down-regulated during all floral bud dormancy stages while up-regulated in the release phase, suggesting their role in the flower developmental process except for PGM and LEA29, which were up-regulated in early dormancy phases, inhibiting the floral development process. The KEGG annotation corresponding to these genes emphasized their roles in embryogenesis during dormancy and dormancy release. The remaining 20 genes were down-regulated during the initial dormancy stages, including para-dormancy, endo-dormancy, and eco-dormancy (Supplementary Table S8).

We further studied the genes associated with ABA synthesis, GA synthesis, cytokinesis, and transcription regulation. Embryonic ABA synthesis played a critical role during dormancy. The identified DEGs concerning ABA synthesis, i.e., ASR2, ABFs, which were down-regulated during the deep dormancy period, including endo-dormancy, para-dormancy, and eco-dormancy, induced the high ABA production to promote dormancy. They were up-regulated at the release of dormancy to lower the production of ABA for the breaking of dormancy. However, the transcripts of the SAUR family and RPL9e showed up-regulated expression patterns during dormancy. Late embryogenesis abundant (LEA) protein was identified with an up-regulated expression pattern throughout the dormancy period. Few gibberellin-induced proteins (related to GA synthesis) such as GID1 were identified with down-regulated expression patterns during the deep dormancy period but up-regulated at the dormancy release phase (Supplementary Table S7).

Among 301 transcription factors (TF), the most prominent TFs of identified DEGs were zinc finger (38), CAD (26), MADS-box (24), MYB (17), WRKY (9), HSP (7), MYC2 (6), and SPT5 (6). Most TFs depicted a trend of down-regulation during dormancy while up-regulated at the dormancy release phase. However, temperature-sensitive TFs were up-regulated during the dormancy period but did not express themselves at later stages of dormancy release (Supplementary Table S7).

The weighted gene co-expression network analysis (WGCNA) is the analytical tool to define correlation networks among variables for connecting them and clustering the highly co-expressed ones for shaping modules. The FPKM values of expressed genes related to the four stages were utilized in WGCNA to expose highly correlated genes. An adjacency matrix was constructed based on DEGs, and 16 discrete modules of genes related to phenotypic, physiological, and biochemical traits expressed for the regulation of floral bud dormancy were generated. Each module represented a cluster of genes depicted with arbitrary colors in the form of a cluster gram coupled with a network heatmap (Figure 7).

Out of the 16 modules, six significantly correlated with phenotypic, physiological, and biochemical traits. It suggests a potential role in the regulation of floral bud dormancy. The MEblue, MEpink4 and MEthistle4 modules displayed the maximum positive correlations with the traits investigated. The details of transcript-based modules were graphically elaborated in Figure 7. Selection of genes from these three modules was accompanied by estimation of their nodes and edges with the help of the R package WGCNA to visualize gene networks. Visualization of gene networks was performed by exporting the generated information into Cytoscape. We further moved towards recruiting the top 40 putative genes from the three significantly correlated gene modules, thus contributing to primary functions in the gene networks (Figure 7). The MEblue module comprised six top eigengenes (hub genes), viz, DUH013357.1, DUH022017.1, DUH011830.1, DUH010324.1, DUH013714.1, and DUH028183.1. The MEpink4 module revealed seven hub genes, i.e., BGI_novel_G003073, DUH016333.2, DUH003125.1, DUH025079.1, DUH003517.1, DUH002558.1, and BGI_novel_G004463. The METhistle4 module also contained seven hub genes, i.e., DUH019196.1, DUH008230.1, DUH0238828.1, DUH001404.1, DUH011979.1, DUH012339.1 and BGI_novel_G002949 (Figure 7). A hub gene from the blue module DUH022017.1 encoded Glycine-rich RNA binding protein (GR-RBP), which was produced as a cold shock response. DUH011830.1 encoded heavy metal-transporting ATPases (HMAs) involved in flower development by absorbing minerals or heavy metals from intracellular membrane-bound organelles. DUH010324.1 hub gene from the blue module encodes pectin methylesterase inhibitor (PMEI), with functions in the metabolism of sucrose, starch, and carbohydrate catabolism by regulating the activity of enzymes.

A hub gene from the pink4 module, viz, DUH003125.1, encodes a late embryogenesis abundant (LEA_2) protein with a function in transport and catabolism of protein from other plant parts to floral buds to regulate the dormancy phenomenon. DUH016333.2 and DUH003517.1 from pink4 module encodes a SERPIN protein and transmembrane protein, playing functions of RNA transport in defense mechanisms and metabolism of terpenoids and polyketides. DUH011979.1 hub gene from the thistle4 module encoded Thioredoxin (TRX) protein with oxidoreductase activity during the floral bud dormancy period. Another hub gene, DUH023828.1 from the thistle4 module, encodes Photosystem II 5 kDa protein (PSBT) with a significant role in dark reactions of photosynthesis, but does not use sunlight energy directly from the sun. One hub gene from the thistle4 module, DUH008230.1, encoded Methanol O-anthraniloyltransferase (AMAT), has a function of linolenic acid metabolism and related transferase activity. There are five hub genes, i.e., BGI_novel_G004463, DUH019196.1, DUH002558.1, BGI_novel_G003073, and DUH001404.1, with either uncharacterized proteins or with unknown functions.

The accuracy of the transcriptome sequencing data is the prerequisite for identifying differentially expressed genes and subsequent enrichment analysis of GO and KEGG functions. For transcriptome data validation, 19 genes with diverse expression patterns at the five floral bud dormancy stages were chosen for the quantitative PCR analysis. qRT-PCR analysis revealed that the trend of gene expression by transcription measured (FPKM) and relative expression by the qRT-PCR method was consistent, and a commendable correlation, i.e., R ² = 0.8129, has been detected between the pair of qRT-PCR data and RNA-seq data (Figure 8). This amount of correlation suggests the reliability of RNA-seq data gathered in the current bundle of studies.