Multiple organ plant materials were obtained from a 3-year-old borneol-type C. burmanni, which was cultivated at the Guangdong Academy of Forestry Sciences (113°37′37.36″N, 23°19′42.08″E, Guangzhou, China). The mother plant originated from a C. burmanni individual with a highly borneol-expressing chemotype, which was selectively bred by our research team in Tianlu Lake Forest Park (113°25′54″N, 23°15′33.4″E, Guangzhou, China). Among the propagated cuttings, this individual was one of the few successfully rooted branches and exhibited the most vigorous growth when sampled. The cuttings were cultivated in a potting mix that was composed of 90% peat soil and 10% perlite (v/v). A rooting powder (including naphthaleneacetic acid, indolebutyric acid, and talcum powder) that was developed by the research team was applied. The climatic conditions included the following: the annual mean temperature was 29.3°C, the mean relative humidity was 77%, extreme highs were 39°C, extreme lows were 3°C, annual precipitation was 1,353.5 mm, and the frost-free period was ≥340 days. The Resource Management Department of Longdong Forest Farm granted permission for the collection of plant materials. The researchers performed hormone metabolomic and transcriptomic analyses on the roots, stems, and leaves, with three biological replicates prepared per tissue. Additionally, nine samples (roots, stems, and leaves) were pooled for full-length transcriptome sequencing. Upon collection, all specimens were swiftly frozen with liquid nitrogen and kept at −80°C.

Metabolite profiling was performed using liquid chromatography–tandem mass spectrometry (LC-MS/MS), according to methods in Chen et al. (2013) with modifications (China). Nine samples were frozen, dried using a vacuum freeze-dryer (Scientz-100F, Ningbo Scientz Biotechnology Co., Ltd., Ningbo, China), and pulverized into a fine powder using a grinding machine with zirconium oxide beads (30Hz, 1.5min). Powder (100 mg) was weighed and sequentially added with 1 mL methanol/water/formic acid (15:4:1, v/v/v) and 10 μL internal standard (100 ng/mL). After vortexing (10min) and centrifugation (12,000 rpm/min, 5min, 4°C), supernatants were concentrated, redissolved in 80% methanol (100 μL), and filtered using a 0.22-μm membrane filter (SCAA-104, 0.22-μm pore size; ANPEL, http://www.anpel.com.cn/) for analysis.

Chromatographic separation was achieved using a Waters ACQUITY UPLC HSS T3 C18 column (1.8 µm, 100 × 2.1mm) with mobile phases A (0.04% acetic acid in water) and B (0.04% acetic acid in acetonitrile) (Niu et al., 2014). The elution protocol was as follows: 0min, A:B = 95:5 (V:V); 1min, A:B = 95:5 (V:V); 8min, A:B = 5:95 (V:V); 9min, A:B = 5:95 (V:V); 9.1min, A:B = 95:5 (V:V); and 12min, A:B = 95:5 (V:V) (flow rate, 0.35 mL/min; injection volume, 2 μL; column temperature, 40°C). MS detection employed a QTRAP^®6500+LC-MS/MS system (https://sciex.com.cn/) in multiple reaction monitoring (MRM) mode with an ESI Turbo ion spray interface (voltage, ± 5500/4500V; temperature, 550°C; curtain gas, 35psi) (Niu et al., 2014). Data were acquired using Analyst 1.6.3 (AB Sciex, Framingham, MA, USA) and analyzed using MultiQuant 3.0.3. According to the Metware Database, a triple quadrupole mass spectrometer was employed for qualitative analysis of the data, while quantitative analysis of the data was performed using the MRM mode in a triple quadrupole mass spectrometer (MS2 tolerance, 20 ppm; Retention Time (RT) error, 0.2min).

Hormone standard solutions (0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50, 100, 200, and 500 ng/mL) were prepared, with l-tryptophan (TRP) and salicylic acid 2-O-β-glucoside (SAG) standard curves ranging from 0.2 to 10,000 ng/mL (20 times dilution). The mass spectral peak intensities of the quantitative signals of hormone standards were recorded at each concentration, and standard curves were generated ( Supplementary Table 1 ).

To enhance the normal distribution, metabolite data were log2-transformed for normalization. The prcomp statistical function from the R software (version 4.2.2) was employed for principal component analysis (PCA) (base package, version 3.5.1) of metabolites across nine samples. Subsequently, hierarchical cluster analysis (HCA) was conducted, and the ComplexHeatmap package (Gu et al., 2016) from the R software was utilized to construct heatmaps. In addition, the R software package MetaboAnalystR (https://gitcode.com/gh_mirrors/me/MetaboAnalystR) was employed to construct an orthogonal partial least squares discriminant analysis (OPLS-DA) model to select significantly differential hormones among distinct samples, with Fold Change ≥ 2 and Fold Change ≤ 0.5 as the selection thresholds. The Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg) was employed for the enrichment analysis of the significantly differential hormones. The absolute hormone levels were calculated based on the peak area ratios with a standard curve using the linear equation. Data visualization was conducted using GraphPad Prism v6.01 (https://www.graphpad.com).

The extraction of the total RNA of each sample was performed using the FastPure^® Universal Plant Total RNA Isolation Kit. A 1% agarose gel electrophoresis (180V) was employed to assess RNA integrity and potential DNA contamination. The OD260/230 and OD260/280 values were measured using the NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA) to determine the purity of total RNA. Next, precise determination of the RNA concentration was conducted using the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and evaluation of RNA integrity was conducted employing the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

The Oligo(dT) magnetic bead approach was applied to construct a full-length transcriptome library in this research. The mRNAs that contained polyA were enriched with Oligo(dT) magnetic beads. Subsequently, a SMARTer PCR cDNA synthesis kit was utilized to transform mRNA to cDNA, which was further enriched for synthesis. The optimal conditions for downstream PCR were determined by optimizing experimental procedures. Large-scale PCR was performed using Oligo(dT) bead-selected fragments to construct a SMRTbell library. The cDNA of full-length was linked to sequencing adapters. Subsequent to the annealing process, the SMRTbell template was utilized as sequencing primers bound to the polymerase. The PacBio Sequel platform was subsequently employed to conduct a 10-hour sequencing run, which utilized P6-C4 chemistry. The raw data of the full-length sequencing have been deposited in the NCBI Sequence Read Archive (SRA) with project number PRJNA1287254.

The initial sequencing data were filtered using the PacBio official software package SMRTlink v8.0 (https://www.pacb.com/support/software-downloads/) with the following parameters: the minimum length was configured as 50 bp, the maximum length was adjusted to 15,000 bp, and the minimum number of full passes was set to one. After correcting the subreads, Circular Consensus Sequences (CCSs) were obtained. CCSs were categorized as Full-Length (FL) sequences and non-Full-Length (nFL) sequences based on the presence of 5′, 3′ primers and polyA tail signals. Within the FL sequences, the Full-Length Non-Chimeric (FLNC) sequences corresponding to the same transcript were analyzed using a hierarchical clustering algorithm [hierarchical nlog(n)] to facilitate the removal of redundant sequences and the derivation of cluster consensus sequences. Subsequently, cluster consensus sequences were polished and yielded consensus sequences with high fidelity. Consensus sequences with high accuracy were utilized via the Arrow (Chin et al., 2013) for the generation of polished consensus sequences for subsequent analyses.

The polished consensus sequences were employed using the GMAP v2017-06-20 (Genomic Mapping and Alignment Program) software to align to the reference genome. The structurally annotated transcripts were subsequently employed for functional annotation using five primary databases. Fictional annotations of these transcripts were performed against KEGG, NR (Non-Redundant Protein Database), and Swiss-Prot (a manually annotated and reviewed protein sequence database) using Diamond blastx v0.8.36 (https://github.com/bbuchfink/diamond). An annotation of Pfam (Protein Family Analysis and Modeling) was performed utilizing Hmmscan v3.1b2 (http://hmmer.org/download.html). Functional annotation of Gene Ontology (GO) was performed by applying NR annotation data and the Blast2GO software (Wu and Watanabe, 2005). Finally, the number of successfully annotated transcripts was statistically analyzed.

The TAPIS (Conesa et al., 2005) software was applied to classify and characterize the full-length transcripts of C. burmanni. After aligning to the polished consensus reference genome, the high-quality isoforms were obtained by further correction, clustering, and redundancy removal. Subsequently, transcript characteristics, novel gene and novel transcript identification, novel gene database annotation, and transcription factor analysis and prediction were conducted. The prediction of plant transcription factors was achieved using iTAK (https://github.com/kentnf/iTAK/) (Abdel-Ghany et al., 2016). Coding potential prediction of PacBio sequencing data was carried out using CNCI (Sun et al., 2013), PLEK (Li et al., 2014), the CPC software (Kong et al., 2007), and the Pfam database (Finn et al., 2016).

The integrity and potential DNA contamination of the RNA sample were assessed employing 1% agarose gel electrophoresis (180V). Subsequently, the NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, USA) was employed to ascertain the OD260/230 and OD260/280 ratios to assess the purity of the RNA. Accurate quantification of RNA concentration was performed using the Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and precise quantification of RNA integrity was performed using the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA).

In this study, 1 μg of total RNA was utilized as the primary sample to construct the sequencing library of transcriptome utilizing the NEBNext^® Ultra™ RNA Library Prep Kit (Illumina). Before library construction, the Qubit 2.0 Fluorometer was used to carry out initial quantification and then dilute the library to 1.5 ng/μL. The effective concentration of the library was precisely quantified via qRT-PCR. Upon passing quality control, the Illumina sequencing pipeline was initiated. The raw data of the RNA sequencing have been deposited in the NCBI SRA with project number PRJNA1286930.

To generate clean reads for further analysis, raw sequencing data were treated to filter out sequences of adapters and reads of poor quality using Fastp v0.19.3 (Chen et al., 2018). After downloading the reference genome and its annotation files, the index was generated using the software HISAT v2.1.0 (Kim et al., 2015). Subsequently, the clean reads were aligned against the reference genome. Fragments Per Kilobase per Million mapped reads (FPKM) was used to evaluate the expression level of genes. DESeq2 v1.22.1 and edgeR v1.22.1 were utilized to assess differentially expressed genes among the sample groups. The Benjamini and Hochberg procedure was executed to correct for multiple hypothesis testing by adjusting the probability values (p-values) to yield the false discovery rate (FDR). The threshold for DEGs was determined in accordance with the following criteria: |log2(fold change)| ≥ 1 and q-value < 0.05.

The numbers of DEGs were counted for pairwise comparisons among roots, stems, and leaves as the three groups. The NR, GO, KEGG, Swiss-Prot, and Pfam databases were utilized to annotate all DEGs. Next, the DAVID database (https://david.ncifcrf.gov/home.jsp) was employed for the functional enrichment analysis of GO and KEGG pathways of all DEGs.

To forecast putative TFs that connected with DEGs, which were relevant to hormone signal transduction within three organs in C. burmanni, the protein sequences of these DEGs were loaded into PlantTFDB v5.0 (https://planttfdb.gao-lab.org/prediction.php). This database identifies TFs that were based on input sequences, family classification rules, and predefined thresholds (Tian et al., 2020).

Based on transcriptomic data, the cor() function within the stats package of the R software (version 4.2.2) was utilized to predict Pearson’s product-moment correlation coefficients (PCCs) between TFs and DEGs (Bai et al., 2022). The correlations with statistical significance (p < 0.05) and absolute values > 0.8 were retained. Co-expression networks were visualized using Cytoscape (version 3.9.1) (Kohl et al., 2011). This study details the comprehensive bioinformatics workflow, as depicted in Supplementary Figure 1 , which encompasses analyses of full-length transcriptomics, transcriptomics, and metabolomics.

The 10 most abundant hormones, including GA1, GA3, gibberellin A20 (GA20), indole-3-acetyl-l-valine methyl ester (IAA-Val-Me), N ⁶-isopentenyl-adenine-9-glucoside (iP9G), meta-topolin-9-glucoside (mT9G), ortho-topolin-9-glucoside (oT9G), 2-oxindole-3-acetic acid (OxIAA), 4-[[(9-beta-d-glucopyranosyl-9H-purin-6-yl)amino]methyl]phenol (pT9G), and TRP, of different organs were selected for weighted gene co-expression network analysis (WGCNA) using the WGCNA 1.71 package from the R software (version 4.2.2), with FPKM values (R² > 0.85) set at an appropriate threshold (Wang et al., 2022). A hierarchical clustering tree was generated based on the correlation of intergenic expression levels between gene expression and the abundance of hormones. Dynamic tree cutting analysis was utilized to segregate genes exhibiting significant association with the 10 selected hormones into discrete modules, followed by the amalgamation of analogous modules.

The investigation utilized the FastPure^® Universal Plant Total RNA Isolation Kit for the extraction of RNA and employed the HiScript III RT SuperMix for qPCR (+gDNA wiper) for the reverse transcription process. A total of 12 differentially expressed genes associated with auxin, gibberellin, and cytokinin were randomly chosen and confirmed through qRT-PCR experiments. Root, stem, and leaf samples were collected, and triplicate technical and biological replicates were conducted. The qRT-PCR assay used Actin7 as the internal standard gene. The primer sequences were designed based on the study by Shi et al. (2024). A specific primer was designed using the PrimerQuest™ Tool (https://sg.idtdna.com/pages/products/qpcr-and-pcr). The detailed sequences of the primer are presented in Supplementary Table 2 . The reaction condition was set as follows: 10 µL 2× ChamQ Universal SYBR qPCR Master Mix, 1.0 µL cDNA, 1.0 µL Primer1 (10 µM), 1.0 µL Primer2 (10 µM), and 8.2 µL ddH₂O. The procedures of qRT-PCR cycling were configured in the following manner: pre-denaturation at 95°C for 30 seconds, denaturation at 95°C for 10 seconds, and annealing at 60°C for 30 seconds, with 40 cycles. The comparative level of DEG expression was quantified through the 2^−ΔΔCT approach (Livak and Schmittgen, 2001).

Metabolomic analysis was conducted to comprehensively identify the significantly differential hormones present in the root, stem, and leaf tissues of C. burmanni. Differential hormones were observed in the roots, stems, and leaves based on the PCA ( Figure 2B ). In our investigation, 72 differential hormones were identified ( Supplementary Table 3 ), including 32 CTK hormones (44.44%), 19 auxin hormones (26.38%), seven GA hormones (9.7%), nine jasmonic acid (JA) hormones (12.5%), two salicylic acid (SA) hormones (2.8%), two ABA hormones (2.8%), and one ethylene (ETH) hormone (1.4%). Hormones exhibiting a Fold Change of ≥2 or ≤0.5 were deemed to be significantly differentially expressed, resulting in the identification of a total of 70 such hormones. In this study, auxin hormones exhibited significant enrichment within the LEAF category, with a total of 11 instances, while cytokinin hormones were predominantly found in the STEM category, with 14 instances. Additionally, gibberellin hormones were significantly enriched in STEM, with a count of three (refer to Figure 2D ). Notably, differential metabolites tryptamine (TRA) and trans-zeatin riboside (tZR) were not categorized as significantly differential hormones, as indicated in Supplementary Table 4 . Venn diagram analysis of significantly differential hormones enriched in the three differential organs showed that 21 significantly differential hormones were commonly enriched across all three tissues ( Figure 2A ), including 10 CTK hormones [6-benzyladenine (BAP), 6-benzyladenosine (BAPR), cis-zeatin (cZ), cis-zeatin-O-glucoside riboside (cZROG), N ⁶-isopentenyladenine (IP), kinetin (K), mT9G, oT9G, pT9G, and trans-zeatin-O-glucoside (tZOG)], five IAA hormones [IAA, indole-3-acetyl-l-tryptophan (IAA-Trp), 3-indoleacetamide (IAM), 3-indolebutyric acid (IBA), and indole-3-carboxaldehyde (ICAld)], two GA hormones [GA3 and gibberellin A24 (GA24)], two JA hormones [JA and 3-oxo-2-(2-(Z)-pentenyl) (OPC-4)], and two SA hormones (SA and SAG). Notably, 3-indoleacetonitrile (IAN) was a unique, significantly differential hormone in the LEAFvsROOT group.

The metabolic pathways enriched by significantly differential metabolites in the three organs were revealed through KEGG enrichment analysis. All groups exhibited enrichment in the following hormone-related pathways: “alpha-Linolenic acid metabolism” (ko00592), “Diterpenoid biosynthesis” (ko00904), “Plant hormone signal transduction” (ko04075), “Tryptophan metabolism” (ko00380), and “Zeatin biosynthesis” (ko00908). However, the pathways “Carotenoid biosynthesis” (ko00906) and “Phenylalanine, tyrosine and tryptophan biosynthesis” (ko00400) were exclusively accumulated in ROOTvsLEAF and ROOTvsSTEM ( Supplementary Figures 2A–C ). To further investigate the variations in hormone content across differential organs, 10 hormones that were relatively highly expressed in three distinct organs were selected for comparative analysis, including GA1, GA3, GA20, IAA-Val-Me, iP9G, mT9G, oT9G, OxIAA, pT9G, and TRP. The results demonstrated that LEAF contained higher levels of GA20, iP9G, mT9G, and TRP compared to other organs. STEM exhibited higher levels of GA1, IAA-Val-Me, OxIAA, and pT9G, while ROOT showed higher levels of GA3 and oT9G ( Figures 2C, D ).

This investigation performed full-length transcriptome analysis to obtain high-quality and as complete as possible transcript sequences in C. burmanni. After sequencing, 43.79-GB raw data and 384,075 polymerase reads were obtained, with a polymerase read N50 of 2,931 bp. The distribution of polymerase read lengths is demonstrated in Figure 3A . After filtering raw data, 18,738,519 subread sequences were obtained. Upon self-alignment correction of each subread sequence, 354,555 CCSs were acquired. The length distribution of CCSs matched the expected pattern, as illustrated in Supplementary Figure S3A . According to the 5′, 3′ primers and tail signals of polyA, the CCSs were categorized into FL and nFL sequences. A total of 316,567 FL sequences and 37,988 nFL sequences were obtained.

Among the FL sequences, 315,605 FLNC and 962 Full-Length Chimeric (FLC) sequences were separated. The distribution of FLNC sequence lengths is depicted in Supplementary Figure S3B . First, after clustering and deduplicating, FLNC sequences were utilized to obtain consensus sequences. Subsequently, consensus sequences were adjusted using the Arrow software; 34,245 polished consensus sequences were obtained, and the distribution of polished consensus sequence lengths is displayed in Figure 3B , which demonstrates that the majority of these sequences were approximately 2,000 bp in length. The Genomic Mapping and Alignment Program (GMAP) was employed to align between polished consensus sequences and the reference genome, resulting in 30,027 (87.68%) successfully mapped consensus sequences (total mapped).

Among the 15,483 structurally annotated transcripts, 15,043 (97.16%) had annotations no fewer than one within the NR, GO, KEGG, Pfam, or Swiss-Prot database ( Supplementary Table 5 ). As illustrated in Figure 3C , the annotated transcripts of the highest quantity are in the NR database, and those of the lowest quantity are in the GO database. Additionally, 10,122 (65.37%) transcripts were annotated in all five databases (NR, GO, KEGG, Swiss-Prot, and Pfam) ( Figure 3D ).

In the NR database, 14,934 transcripts (96.45%) were annotated. The alignment results of these transcripts with C. burmanni revealed that the four species with the highest similarity were Nelumbo nucifera (7,157, 20.90%), Vitis vinifera (1,272, 8.22%), Elaeis guineensis (1,151, 7.43%), and Phoenix (817, 5.28%) ( Supplementary Figure 4A ). To uniformly describe gene functions, 10,785 transcripts (69.66%) were subdivided into three GO categories: biological process (BP), cellular component (CC), and molecular function (MF). In the BP category, the top three subgroups were “cellular process” (GO:0009987), “metabolic process” (GO:0008152), and “single-organism process” (GO:0044699). In the CC category, the top three subgroups were “cell” (GO:0005623), “cell part” (GO:0044464), and “organelle” (GO:0043226). In the MF category, the top three subgroups were “binding” (GO:0005488), “catalytic activity” (GO:0003824), and “transporter activity” (GO:0005215) ( Supplementary Figure 4B ).

A total of 14,868 (96.02%) transcripts were successfully annotated in the KEGG database, which were classified into six branches: Cellular Processes, Environmental Information Processing, Genetic Information Processing, Human Diseases, Metabolism, and Organismal Systems, which elucidated the roles and metabolism routes of gene expression products in cells. The five most abundant subgroups were “Signal transduction” (ko02010, ko04016, ko04070, and ko04075) under Environmental Information Processing, “Carbohydrate metabolism” (ko01200) under Metabolism, “Transport and catabolism” (ko04136, ko04144, ko04145, and ko04146) under Cellular Processes, and “Folding, sorting and degradation” (ko03018, ko03050, ko03060, ko04120, ko04122, and ko04141) and “Translation” (ko00970, ko03008, ko03010, ko03013, ko03015, and ko03040) under Genetic Information Processing ( Supplementary Figure 4 C).

TFs represent a class of protein molecules that specifically occupy the promoter of target genes to quantitatively modulate the expression of genes with temporal precision and spatial specificity. They currently serve as potent tools widely applied in genetic engineering. We utilized the iTAK software (Zheng et al., 2016) to perform predictive analysis of transcription factors, which identified 922 TFs from 29 gene families ( Figure 4 ). The bHLH family exhibited the highest membership (63 members, 6.83%), followed by WRKY (53 members, 5.75%), C3H (48 members, 5.21%), MYB-related (47 members, 5.09%), and NAC families (44 members, 4.77%). In contrast, the AUX/IAA, HSF, and C2C2-Dof families all displayed the lowest representation, which contained 14 members (1.51%). These findings will provide significant implications for the investigation of transcriptional regulation in related biological processes and facilitate the construction of an accurate, high-quality, and comprehensive transcription factor library.

Transcriptome sequencing analysis of RNA samples from three organs was performed to discover the DEGs that are associated with hormone signal transduction in the root, stem, and leaf organs of C. burmanni, with three biological repeats, which were performed to confirm the validity of the data. The sequencing yielded raw reads ranging from 40.13 to 47.94 Mb. After quality filtering, error frequency examination, and Guanine-Cytosine (GC) composition allocation assessment, clean reads from the nine cDNA libraries ranged from 38.72 to 47.77 Mb, totaling 61.04 Gb of clean bases, with an average of 6.78 G per sample. The Q20 base percentage for all samples consistently exceeded 97%, and the Q30 base percentage exceeded 93% ( Supplementary Table 6 ). For each library, the count of reads that demonstrated successful alignment to the reference genome varied between 30,110,309 (70.76%) and 43,870,063 (93.60%), while uniquely mapped reads ranged from 29,131,403 (68.46%) to 41,971,195 (89.55%) ( Supplementary Table 7 ). PCA ( Figure 5B ) demonstrated statistically significant separation among samples from distinct organs, with specimens derived from identical tissues clustering together, which indicated substantial divergence between the three differential groupings. Correlation analysis of C. burmanni organs revealed excellent reproducibility across all sample groups ( Supplementary Figure S5A ). Additionally, 812 DEGs were commonly enriched in all three comparative groups ( Figure 5A ). In the ROOTvsLEAF group, 7,479 DEGs exhibited significantly differential expression, including 4,422 upregulated (59.13%) and 3,057 downregulated (40.87%) ( Supplementary Figure S5B ). The STEMvsLEAF group contained 4,264 DEGs, with 1,661 upregulated (38.95%) and 2,603 downregulated (61.05%) ( Supplementary Figure S5C ). The ROOTvsSTEM group displayed 5,061 DEGs, comprising 2,792 upregulated (55.17%) and 2,269 downregulated (44.83%) ( Supplementary Figure S5D ).

The enrichment analyses of GO and KEGG were implemented to elucidate the biological implications of DEGs in this study. The enrichment analysis of GO categorized the functions of DEGs, which demonstrated that 20 GO terms with the greatest enrichment classified into three major GO categories: BP, CC, and MF. Specifically, the BP category contained 14 substantially enriched GO terms, the CC category comprised four, and the MF category included two. Within the three differential groups, the highest quantity of DEGs was enriched in the CC term “cellular anatomical entity” (GO:0110165), followed by the BP term “cellular process” (GO:0009987), the MF term “binding” (GO:0005488), the BP term “metabolic process” (GO:0008152), and the MF term “catalytic activity” (GO:0003824) ( Supplementary Figures 6A–C ). The KEGG pathway analysis provided additional functional insights into the pathways linked to DEGs. The KEGG results revealed that these functional pathways were primarily categorized into five major categories: Cellular Processes, Environmental Information Processing, Genetic Information Processing, Metabolism, and Organismal Systems. Among three distinct groups, the most abundant subpopulations were “Metabolic pathways” (ko01100) and “Biosynthesis of secondary metabolites” (ko01110) under the Metabolism category, as well as “Plant−pathogen interaction” (ko04626) under the Organismal Systems category ( Supplementary Figures 7A–C ). KEGG enrichment analysis revealed that the DEGs were predominantly enriched in the following pathways: “Biosynthesis of secondary metabolites” (ko01110), “Photosynthesis” (ko00195), “Photosynthesis-antenna proteins” (ko00196), and “Porphyrin metabolism” (ko00860) ( Figures 5C–E ).

To investigate the DEGs that are associated with hormone signal transduction across differential organs, this study selected 10 highly expressed hormones from multiple organs, including three IAA hormones (IAA-Val-Me, OxIAA, and TRP), three GA hormones (GA1, GA3, and GA20), and four CTK hormones (iP9G, mT9G, oT9G, and pT9G). Based on transcriptome profiling data, WGCNA was applied with a filtering threshold set at 0.85 (85% of genes were excluded, retaining 15% of the data for analysis). By applying hierarchical clustering analysis, all genes were partitioned into 10 modules, each containing genes with similar expression patterns and labeled with distinct colors ( Figure 6A ). The five largest modules were turquoise (845 genes), blue (576 genes), brown (340 genes), yellow (329 genes), and green (261 genes) ( Figure 6B ).

Seven modules exhibiting significant correlations with the hormone signaling pathway were identified from the WGCNA results. Specifically, the green module demonstrated a significantly positive correlation with TRP (p < 0.01) and displayed strongly negative correlations with GA1, IAA-Val-Me, OxIAA, and pT9G (p < 0.01). The turquoise module demonstrated a positive correlation with TRP (p < 0.05); exhibited significantly positive correlations with GA20, iP9G, and mT9G (p < 0.01); and displayed negative correlations with GA1, OxIAA, and pT9G (p < 0.05) and a strongly negative correlation with IAA-Val-Me (p < 0.01). The yellow module demonstrated a significantly positive correlation with GA3 and oT9G (p < 0.01). The blue module displayed positive correlations with OxIAA and pT9G (p < 0.05), exhibited strongly positive correlations with GA1 and IAA-Val-Me (p < 0.01), and demonstrated a negative correlation with TRP (p < 0.05). The black module exhibited negative correlations with GA20, iP9G, and mT9G (p < 0.05). The brown module displayed significantly positive correlations with GA1 and pT9G (p < 0.01) and exhibited positive correlations with IAA-Val-Me and OxIAA (p < 0.05). The pink module demonstrated strongly positive correlations with OxIAA and pT9G (p < 0.01), exhibited positive correlations with GA1 and IAA-Val-Me (p < 0.05), and displayed a negative correlation with TRP (p < 0.05) ( Figure 6B ). Furthermore, seven significantly correlated modules were subjected to the KEGG pathway analysis using the graphics package of R (version 3.5.1). The results demonstrated that the turquoise module (13 genes) contained the highest count of DEGs connected to hormone signal transduction, whereas the yellow module (three genes) exhibited the lowest quantity of DEGs associated with hormone signaling pathways ( Figure 6C ).

Transcriptomic and metabolomic analyses were integrated to construct a heatmap, which was linked to phytohormone signal transduction and metabolic regulation pathways to explore DEGs across three organs of C. burmanni. A total of 50 DEGs related to hormone signal transduction were determined ( Supplementary Table 8 ). Specifically, 19 DEGs were linked to auxin, with eight upregulated in leaves, seven in stems, and four in roots; 12 DEGs were associated with cytokinin, including four with highly expressed levels in leaves, five in stems, and three in roots; 19 DEGs were related to gibberellin, among which six exhibited elevated expression in leaves, two in stems, and 11 in roots ( Figure 7A ). Twelve candidate DEGs were randomly chosen to perform qRT-PCR assays for validating the patterns of gene expression. The experimental results confirmed consistency between the observed expression profiles and bioinformatics predictions ( Figure 7B ).

The Mantel test correlation heatmap analysis revealed that the association patterns of three phytohormones (IAA, CTK, and GA) with significantly DEGs in the root, stem, and leaf organs exhibited remarkable tissue specificity. In leaves (LEAF), GA demonstrated the most prominent correlation effects, showing highly significant associations with nine DEGs (p < 0.01) and significant correlations with eight DEGs (p < 0.05), with the number of significantly correlated genes substantially exceeding that of IAA and CTK. IAA exhibited significant correlations with only two genes (p < 0.05) ( Figure 8A ). In roots (ROOT), both GA and CTK displayed extensively regulatory capacities. GA correlated with 10 DEGs at highly significant levels (p < 0.01) and five DEGs at significant levels (p < 0.05), while CTK showed one DEG at highly significant levels (p < 0.01) and nine DEGs with significant correlations (p < 0.05). The number of IAA-associated genes was fewer, with only six exhibiting significant/highly significant levels ( Figure 8B ). In stems (STEM), GA dominated the regulatory network, which exhibited highly significant correlations with 10 DEGs (p < 0.01) and significant associations with two DEGs (p < 0.05). CTK displayed significant correlations with only two genes (p < 0.05), but IAA showed no significant DEG associations ( Figure 8C ).

The DEGs that were associated with the three hormones (IAA, CTK, and GA) were employed to predict and screen TFs using PlantTFDB v5.0 (Tian et al., 2020) ( Supplementary Table 9 ). Among these, TRI1, AUX/IAA, GH3, SAUR, CRE1, AHP, A-ARR, and GID1 were excluded during the prediction process since they function as auxiliary regulatory proteins instead of transcription factors. The prediction results encompassed five transcription factor families, including ARF, bHLH (PIF3/4), G2-like (DELLA), GRAS (GLK/KAN1/2/HH2O/APL/FT), and ARR-B.

To investigate the relationships among DEGs, phytohormones, and TFs in the IAA, GA, and CTK signaling pathways, correlation analysis was conducted. The screening criteria for analysis were established at p < 0.05, and the absolute correlation coefficient surpassed 0.8. The results revealed 39 differential phytohormones and 46 DEGs. Specifically, 15 phytohormones and 18 DEGs were associated with IAA, four phytohormones and 18 DEGs were linked to GA, and 20 phytohormones and 10 DEGs were correlated with CTK. Among the six DEGs participating in the IAA signaling pathway, the expression of Cbur03G008810, Cbur03G017910, Cbur04G003980, Cbur10G012060, Cbur10G026020, and Cbur11G010450 exhibited positive correlations with the ARF transcription factor. For the nine DEGs engaged in the GA signaling pathway, Cbur05G008030 and Cbur10G025060 showed positive correlations with the bHLH transcription factor, whereas Cbur03G007900 and Cbur10G008890 displayed negative correlations with bHLH. Additionally, Cbur02G000190, Cbur04G005370, and Cbur05G003960 demonstrated positive correlations with the GRAS transcription factor, while Cbur01G019820 and Cbur07G013060 exhibited negative correlations with GRAS. Among the eight DEGs that participated in the CTK signaling pathway, Cbur12G017510, Cbur11G010330, Cbur01G030160, Cbur03G017710, Cbur07G000830, and Cbur04G006160 showed positive correlations with G2-like, whereas Cbur04G007700 displayed a negative correlation with the ARR-B transcription factor. Conversely, Cbur01G029630 exhibited a positive correlation with ARR-B ( Figure 9 ).