XHJ was collected from Chetou Village, Lufeng Town, Xupu County, Hunan Province, China. DHJ was collected from Qinglong Town, Jianshui County, Yunnan Province, China. All these plants were transplanted into the medicinal plant garden of Hunan University of Chinese Medicine in January 2022, and the rhizome, roots, and buds were collected in September 2022, when XHJ began to germinate. These tissues were frozen in liquid nitrogen immediately and stored at -80 °C. All of these plants were identified as Polygonatum kingianum var. grandifolium and Polygonatum kingianu, respectively by Associate Professor Yeman Liu (Hunan University of Chinese Medicine). The authenticated specimens of Polygonatum kingianum var. grandifolium (Code: 20220902-1) and Polygonatum kingianum (Code: 20220902-2) were kept in the TCM and Ethnomedicine Innovation & Development International Laboratory, Innovative Materia Medica Research Institute, Hunan University of Chinese Medicine.

The buds from XHJ and DHJ were vacuum-infiltrated and fixed with 50% FAA (v/v) and were also collected in September 2022. The procedure for paraffin section observation was followed as described in our previous study (Ning et al., 2019).

Total polysaccharide was extracted and tested from the dried XHJ and DHJ rhizomes as in the method described in the Chinese Pharmacopoeia. Three repetitions were conducted in each sample and statistical analysis was performed by SPSS 22.0 software.

The rhizome, roots, and buds of XHJ and DHJ were chosen for targeted metabolomics of phytohormones analysis, which were abbreviated as X-rhizome, X-root, X-bud, D-rhizome, D-root, and D-bud, respectively. Each group had three biological replicates, and each biological replicate contained more than 10 plants. This part of the experiments was carried out by MetWare Biotechnology Co., Ltd. (Wuhan, China). In total, 88 phytohormones were detected, including 2 ABAs, 26 auxins, 36 CTKs, 10 GAs, 9 jasmonic acids (JAs), 2 salicylic acids (SAs), 2 strigolactones (SLs) and 1 ethylene (ETH). All samples were vacuum-freeze dried and ground to powder separately. The sample extracts were analyzed by a UPLC-Q-TOF-MS/MS system (ExionLC™ AD, https://sciex.com.cn/, QTRAP^®6500+, https://sciex.com.cn/). The conditions for the high-performance liquid chromatography (HPLC) were as follows: chromatographic separation was conducted on Waters’ ACQUITY UPLC HSS T3 C18 column (1.8 µm, 100 mm × 2.1 mm), using water with 0.04% acetic acid and acetonitrile with 0.04% acetic acid as mobile phase A and B, respectively. The elution program followed was: 95:5 V(A)/V(B) at 0 min, 95:5 V(A)/V(B) at 1 min, 5:95 V(A)/V(B) at 8 min, 5:95 V(A)/V(B) at 9 min, 95:5 V(A)/V(B) at 9.1 min, and 95:5 V(A)/V(B) at 12 min. The flow rate was 0.35 mL/min, the column temperature was 40 °C, and the injection volume was 2 μL. Mass data acquisition was carried out in ion spray voltage in positive (5500 V)/negative (-4500 V) mode, with an electrospray ionization temperature of 550 °C. Curtain gas was set at 35 psi. The content of phytohormones was determined using an internal standard method.

Principal component analysis (PCA) of different samples for comparison was conducted by ‘prcomp’ in R software (www.r-project.org/). The DAPs were determined by log2 fold change ≥ 2 or ≤0.5 and variable importance in projection ≥ 1. The identified DAPs were further analyzed using the KEGG database.

The X-rhizome, X-root, X-bud, D-rhizome, D-root, and D-bud were also used for transcriptomic analysis. Each group had three biological replicates, and each biological replicate contained more than 10 plants. A TRIzol RNA extracting kit was used to extract the total RNA. RNA concentration and integrity were detected. The mRNAs with polyA tails were enriched by Oligo(dT) beads. The double-stranded cDNA was purified using AMPure XP beads and repaired at the end: adding an A-tail and the sequencing linkers. A cDNA library was constructed, and the Illumina Novaseq6000 system platform was used to sequence. CASAVA base recognition was conducted to convert the image data into raw data. The fastq software was used to obtain clean reads (Chen et al., 2018), and the clean reads were stitched using the Trinity software (Grabherr et al., 2011). The transcriptome’s overall quality was assessed for completeness using the Benchmarking Universal Single-Copy Orthologs tool (BUSCO version 5.4.3) (Seppey et al., 2019). The gene expression levels were measured by FPKM (Chen et al., 2021). DESeq2 software was used to identify the DEGs with p-values < 0.05 and |log2 fold changes| ≥ 1. The KEGG database was used to perform the DEG enrichment analysis.

ReverTra Ace^® qPCR-RT Master Mix with gDNA Remover (TOYOBO) was used for total RNA reverse transcription according to the manufacturer’s protocol. Elongation factor 1-α was used as an internal reference (Li D. D. et al., 2022). Specific primers were designed by Primer 5.0 software. An Applied Biosystems QuantStudio 5 system with the Hieff^® qPCR SYBR Green Master Mix (Low Rox Plus) (YEASEN) was used to perform the qRT-PCR. The relative expression was calculated by the 2^-ΔΔCt method (Livak and Schmittgen, 2001). All primer sequences are listed in Supplementary Table S1 .

The DEGs identified by the transcriptomic analysis were used to obtain the co-expression network modules by the WGCNA package in R. The co-expression modules were generated via the automatic network construction function with the default parameters: mergeCutHeight was 0.25, powerEstimate was 18, and minModuleSize was 50. The eigengene values were calculated in each module, which was applied to search for the association with the main compounds in XHJ and DHJ. The module member-ship value was obtained by using the signedKME function of the WGCNA package to analyze the correlation between gene expression and eigengene modules. These networks were visualized using Cytoscape v.3.7.2 software.

Both XHJ and DHJ belong to the whorled leaf group in Polygonatum Mill. However, the flower color of XHJ is greenish-white, while it is red in DHJ ( Figures 1A, B ). The rhizome morphology of XHJ and DHJ was also different ( Figures 1C, D ) as the rhizome node of DHJ was more rounded. The paraffin sections observation showed that the rhizome bud of XHJ in autumn was completely differentiated, which was easy to see in the floral and leaf primordium ( Figure 1E ), but the rhizome bud of DHJ in autumn was undifferentiated and stayed dormant ( Figure 1F ). The total polysaccharide content of XHJ and DHJ rhizome was 102.25 mg/g and 90.91mg/g, respectively ( Figure 1G ).

To fully clarify the different rhizome bud germination mechanisms of XHJ and DHJ in autumn, a targeted phytohormone metabolomic analysis was conducted in XHJ rhizome buds, XHJ rhizomes, XHJ roots, DHJ rhizome buds, DHJ rhizomes, and DHJ roots. Finally, 63 phytohormones were identified in these tissues ( Figure 2A ; Supplementary Table S2 ). The heatmap showed that the X-rhizomes and D-rhizomes were located in the same cluster, the X-buds and D-buds were on the adjacent branch, and the X-roots and D-roots were on another cluster, which indicated that the phytohormones in the same tissues were relatively similar. However, the content of GA19, GA7, GA24, GA4, GA15, and GA9 was significantly higher in the X-buds, and the content of mT9G, cZR, and ABA was significantly higher in the D-buds. The content of IP, IAA-Ala, SA, and SAG was higher in the D-rhizomes. Moreover, the content of iP7G, Strigol, cZ, OxIAA, IAA-Glu, GA20, GA3, JA, JA-ILE, and JA-Phe was particularly higher in the X-roots ( Figure 2A ; Supplementary Table S2 ). Compared with the D-buds, 14 upregulated and 17 downregulated DAPs were identified in the X-buds. Compared with the D-rhizomes, 6 upregulated and 13 downregulated DAPs were identified in the X-rhizomes. Compared with the D-roots, 18 upregulated and 14 downregulated DAPs were identified in the X-roots ( Figure 2B ; Supplementary Table S3 ).

To reveal the different shooting mechanisms between XHJ and DHJ in autumn, RNA-seq of buds, rhizomes, and roots of them were profiled, which generated a total of 818.6 million raw reads ( Supplementary Table S4 ). The final assemblies were well annotated and made available. The transcript N50 value was 1159 bp, the unigene N50 value was 1364 bp, and the total BUSCO score was 98.1% ( Supplementary Table S4 ). A good repetition among all biological replicates in each group was obtained ( Supplementary Figure S1 ). In total, 298,251 unigenes were annotated using the KEGG, NR, Swiss-Prot, GO, COG/KOG, TrEMBL, and Pfam databases ( Supplementary Table S5 ). The DEGs were analyzed and a total of 43,860, 41,520, and 45,276 upregulated DEGs, and 42,420, 42,741, and 41,758 downregulated DEGs (P < 0.05) were identified in D-buds vs X-buds, D-rhizomes vs X-rhizomes, and D-roots vs X-roots, respectively ( Figures 3A–C ). There were 11,074, 11,570, and 15,603 DEGs specifically differentially expressed in the buds, rhizomes, and roots, respectively ( Figure 3D ). Hierarchical cluster analysis (HCA) of all the DEGs was conducted to reveal the global gene expression pattern in different groups ( Figure 3E ). As a result, all the XHJ samples were significantly separated from the DHJ samples.

According to the KEGG pathway analysis, in D-buds vs X-buds, ‘Plant hormone signal transduction’ (ko04075) and ‘Photosynthesis - antenna proteins’ (ko00196) were significantly enriched ( Supplementary Table S6 ). In D-rhizomes vs X-rhizomes, ‘Ribosome biogenesis in eukaryotes’ (ko03008) and ‘Spliceosome’ (ko03040) were significantly enriched ( Supplementary Table S6 ). In D-roots vs X-roots, ‘Proteasome’ (ko03050) and ‘Protein processing in endoplasmic reticulum’ (ko04141) were significantly enriched ( Supplementary Table S6 ).

Furthermore, 11 DEGs related to plant hormone metabolism were selected for qRT-PCR validation ( Supplementary Figure S2 ), including Cluster-21187.19 (YUCCA8-like gene), Cluster-65218.0 (GA20 oxidase 1-D-like), Cluster-56959.5 (GA3-beta-dioxygenase 1-like), and so on. The expression patterns of the most selected DEGs were consistent with the RNA-seq data, which suggested the RNA-seq data in the present study were reliable.

In transcriptome, PCA results showed that PC1 and PC2 explained 43.28% and 8.95% of the total variance, respectively ( Figure 4A ). X-roots and D-roots could be distinguished from X-rhizomes, D-rhizomes, X-buds, and D-buds by PC1 and PC2. X-rhizomes were difficult to distinguish from X-buds by PC1. D-rhizomes were also difficult to distinguish from D-buds by PC1 and PC2 ( Figure 4A ). In the metabolome, the PCA results showed that PC1 and PC2 explained 38.84% and 26.94% of the total variance, respectively ( Figure 4B ). X-roots and D-roots could also be distinguished from other samples by PC1 and PC2. X-rhizomes were difficult to distinguish from X-buds by PC1. D-rhizomes were also difficult to distinguish from D-buds by PC1 and PC2. These results were consistent with the PCA results in the transcriptome ( Figure 4B ). However, the distributions of the X-rhizome, D-rhizome, X-bud, and D-bud groups were more concentrated than those in transcriptome.

To illustrate the relationships between DAPs and DEGs, a nine-quadrant maps were drawn ( Figures 4C–E ). The results showed that most DAPs (36 in D-buds vs X-buds, 35 in D-rhizomes vs X-rhizomes, and 32 in D-roots vs X-roots) and DEGs (18,888 in D-buds vs X-buds, 31,850 in D-rhizomes vs X-rhizomes, and 37,166 in D-roots vs X-roots) were in quadrants 2, 4, 6, and 8, indicating no correlation between these metabolites and genes. DAPs and DEGs in quadrants 3 and 7 were positively correlated ( Figures 4C–E ; Supplementary Table S7 ). In D-buds vs X-buds, 6 DAPs (including tZR, GA9, GA24, IAA, IAA-Leu, and IAA-Asp) and 15,339 DEGs were up-regulated in quadrant 3, and 13 DAPs (including ABA, 3-Indole acetamide, salicylic acid 2-O-β-glucoside, and so on) and 49,478 DEGs were downregulated in quadrant 7 ( Figure 4C ; Supplementary Table S7 ). In D-rhizomes vs X-rhizomes, 4 DAPs (including GA3, GA9, GA24, and IAA-Leu) and 7,997 DEGs were upregulated in quadrant 3, and 12 DAPs (including SA, JA, JA-ILE, trans-zeatin (tZ), and so on) and 30,857 DEGs were downregulated in quadrant 7 ( Figure 4D ; Supplementary Table S7 ). In D-roots vs X-roots, 7 DAPs (including GA3, OxIAA, JA-Phe and so on) and 9,607 DEGs were upregulated in quadrant 3, and 6 DAPs (including tZ, cZ9G, TRA, and so on) and 44,327 DEGs were downregulated in quadrant 7 ( Figure 4E ; Supplementary Table S7 ). The common KEGG pathways enriched by the DEGs and DAPs are shown in Figures 5A–C ( Supplementary Table S8 ). Interestingly, ‘Diterpenoid biosynthesis’ (ko00904) and ‘Plant hormone signal transduction’ (ko04075) were collectively enriched in all tissues. Among these DAPs, GAs were significantly different between XHJ and DHJ. As is well known, GAs belong to ‘Diterpenoid’, so the DEGs involved in diterpenoid metabolism were further focused on. In ‘Diterpenoid biosynthesis’ (ko00904), the relative expression levels of the related DEGs are shown in Figure 5D , including the ent-copalyl diphosphate synthase gene (ent-CPS), ent-kaurene synthase gene (KS), ent-kaurene oxidase gene (KO), ent-kaurenoic acid oxidase gene (KAO), GA 3β-dioxygenase gene (GA3β-ox), GA 2β-dioxygenase gene (GA2β-ox), GA20 oxidase gene (GA20ox), momilactone-A synthase gene (MAS), geranyllinalool synthase gene (GES), and GA13 oxidase gene (GA13ox). Among these, ent-CPS, KO, KAO, and GES showed relatively higher expression in XHJ, while KS, MAS, GA3β-ox, GA2β-ox, GA20ox, and GA13ox were upregulated in DHJ ( Figure 5D ). In addition, the relative expression of unigenes related to ETH and SL metabolites were upregulated in XHJ, such as ACC oxidase gene, ACC synthase gene, DWARF27 (D27) isomerase gene, Carotenoid cleavage dioxygenase 7 (CCD7), and CCD8 ( Supplementary Table S9 ).

To gain further insight into the gene regulatory network and identify potential key regulatory factors related to auxin, CTK, GA, and JA metabolism in XHJ and DHJ, a WGCNA was conducted. A total of six co-expression modules were identified according to similar expression patterns ( Figure 6A ; Supplementary Table S10 ). The heatmap of module-trait correlations showed the accumulation of transcripts in the brown module was highly correlated with most of the phytohormones, such as IAA-Glu, GA20, and GA3 ( Figure 6B ).

To generate the potential regulatory network associated with phytohormone metabolism, the structural genes in the brown module were further analyzed (absolute value of Pearson’s correlation coefficient > 0.9, p-value < 0.01) ( Figure 6A ; Supplementary Table S10 ). We identified 12 structural genes related to the auxin metabolic process in the brown module, including Cluster-133654.6 (auxin efflux carrier family protein: PIN-LIKES 3-like), Cluster-78230.2 (SAUR family protein: auxin-induced protein 6B-like), Cluster-106117.2 (GO:2000012, regulation of auxin polar transport), and so on ( Supplementary Table S11 ), which had a high correlation with the accumulation of IAA-Asp, IAA-Glu, ICAld, MEIAA, and OxIAA. By correlating the accumulation patterns of transcripts and the potential binding affinity of the structural genes associated with auxin metabolism, we identified 62 TFs (the top three gene families were AP2/ERFs, WRKY, and NAC), whose expression was highly correlated with the 12 auxin-related structural genes in the brown module and constructed a correlation network ( Figure 7A ; Supplementary Table S12 ).

In addition, five CTK metabolism-related structural genes were identified, including two CTK-activated signaling pathway genes (GO:0009736: Cluster-103520.0 and Cluster-117775.6), two gene responses to CTK (GO:0009735: Cluster-124082.2 and Cluster-124082.3), and one CTK dehydrogenase gene (CKX, Cluster-98370.0), which had a high correlation with the accumulation of cZ, iP7G, DZ, and KR ( Supplementary Table S11 ). Additionally, 58 TFs (the top three gene families were AP2/ERFs, WRKY, and NAC) were identified with a high correlation with the five CTK metabolism-related structural genes in the brown module and we constructed a correlation network ( Figure 7B ; Supplementary Table S12 ).

Moreover, nine structural genes related to the JA metabolic process in the brown module were identified, including four JA ZIM domain-containing proteins (Cluster-103094.0, Cluster-67629.11, Cluster-67629.12 and Cluster-67629.4), one JA-induced oxygenase 2 (Cluster-124125.15), one JA methylesterase (Cluster-127170.5), two JA-induced proteins (Cluster-22029.1 and Cluster-26158.3), and one JA O-methyltransferase (Cluster-90404.0), which had a high correlation with the accumulation of OPC-6, JA-ILE, and JA-Phe ( Supplementary Table S11 ). Furthermore, 59 TFs (such as AP2/ERFs, WRKY, and NAC) were identified to be highly correlated with the nine JA metabolism-related structural genes in the brown module and we constructed a correlation network ( Figure 7D ; Supplementary Table S12 ).