Physiological responses and screening of key genes involved in the cold tolerance of herbaceous peony cultivated over-winter in Harbin, the coldest metropolis of China

Cold stress, encompassing chilling stress (0–15 °C) and freezing stress (< 0 °C), represents a significant abiotic stress causing substantial economic losses in horticultural crops^[1,2]. Typical damage induced by cold temperatures includes impaired photosynthesis, decreased cell membrane fluidity leading to rigidity, destabilization of proteins and protein complexes, and even cell death due to dehydration^[3]. These effects restrict the geographical distribution of plant species and reduce global crop yields, creating worldwide negative impacts on plant systems. The compounded effects of membrane dysfunction, protein denaturation, and structural collapse under cold stress collectively threaten sustainability across diverse ecosystems^[4−6].

Plants have evolved sophisticated regulatory mechanisms to withstand cold stress, such as synthesizing various protective compounds and proteins, including soluble sugars, proline (PRO), cold-resistant proteins, phytohormones, and epigenetics. These substances participate in regulating osmotic potential and cell membrane stability in plants under cold stress^[7]. Additionally, plants activate antioxidant systems by enhancing the activities of enzymes like superoxide dismutase (SOD) and peroxidase (POD), which contribute to reactive oxygen species (ROS) scavenging^[4]. For instance, wheat (Triticum aestivum L.) enhances its antioxidant defense to improve cold resistance, while bamboo (Bambusoideae) and aloe vera (Aloe vera [L.] Burm. f.) boost activities of SOD and POD to elevate amino acid and soluble sugar levels, mitigating cold-induced damage^[8,9]. Furthermore, plants adapt to cold environments through dynamic adjustments in phytohormone levels and balances that modify growth patterns^[10]. Phytohormones such as abscisic acid (ABA) and jasmonic acid (JA) play critical roles in cold adaptation across species like wheat and apple (Malus domestica Borkh.)^[11,12]. The synergistic integration of physiological regulation and metabolic reprogramming constitutes a fundamental adaptive strategy that enables plants to withstand freezing environments^[13]. Plants are also able to modify translation efficiency and methylation status through epigenetics, thereby enhancing cold tolerance via non-nucleotide sequence changes such as structural alterations.

Research on the plant cold tolerance involves morphological observation, physiological measurement, omic sequencing, and bioinformatics-related analysis, and the subsequent identification of crucial genes^[14]. Screening of cold-resistant genes is a pivotal component in such research, bridging a series of studies at morphological, physiological, proteomic, and omics levels, and thus it paves the way for subsequent investigations into gene function and interaction^[15]. Multiple methodologies exist for screening cold-resistant genes. Current mainstream approaches focus on correlating morphological and physiological observations with gene expression data acquired by transcriptome through the weighted gene co-expression network analysis (WGCNA). Adverse environmental conditions activate stress-responsive gene expression through complex signaling pathways. This activation leads to downstream physiological and biochemical adaptations, which makes transcriptomics a key tool for studying stress-related gene expression. This method has been frequently applied in crops such as rice and upland cotton (Gossypium hirsutum L.)^[16,17].

Herbaceous peony (Paeonia lactiflora Pall.) is a world-famous ornamental valued for beautiful flowers, elegant leaves, and medicinal and edible oil uses. It exhibits relatively weak heat tolerance. Therefore, most researches focus on the southward introduction, cultivation adaptability, and heat resistance mechanisms of herbaceous peony^[18−21]. However, the cultivation of this ornamental in high-latitude regions like Northeast China also remains limited. The primary constraint lies in the extreme winter conditions of areas such as Heilongjiang Province, where temperatures routinely plunge −30 °C below, creating deep frozen soil layers. The mechanical stress and cellular dehydration of herbaceous peony's buds and crowns under frozen ground disrupt meristematic tissues, and thus lead to compromised spring sprouting, stunted growth, floral abortion, and low flowering rate^[22]. Despite some preliminary research exploring cold resistance in herbaceous peony, a critical knowledge gap persists due to the insufficient systematic investigations integrating physiological, biochemical, and transcriptomic analyses^[23]. Additionally, previous studies on herbaceous peony cold tolerance also exhibit several other limitations. These include insufficient comparative analyses between low-latitude and high-latitude germplasms and inadequate investigation of tender underground bud organs beneath frozen soil. There is also an incomplete identification of key regulatory genes with limited genetic resources. Therefore, understanding how herbaceous peony crown buds perceive low temperatures and coordinate genetic regulation is crucial^[24].

In this study, two representative herbaceous peony cultivars with contrasting cold-resistances from distinct southern and northern Chinese production regions were selected to explore the physiological regulation of cold-tolerance and mine the crucial associated genes. The study features three novel aspects: comparison of cultivars from totally different climatic regions, conducting experiments in the northernmost Chinese metropolis, and focusing on bud survival under frozen soil instead of traditional research on the aboveground part. This study offers novel insights into plant cold-adaptation mechanisms in extreme high-latitude cold regions. It lays important groundwork for future functional and interaction studies of cold-tolerance genes in herbaceous peonies, and also provides theoretical support for breeding strong cold-resistant herbaceous peony cultivars. These new cultivars could be planted and used in the high-latitude regions.

Two representative cultivars with contrasting provenance and cold tolerance were selected as experimental materials: 'Hang Baishao' (hereafter abbreviated as HB), an herbaceous peony cultivar originating from the lowest latitude region in China^[24]. and 'Hongling Chijin' (hereafter abbreviated as HL), a cold-resistant cultivar of northern origin, identified through introduction and screening in Harbin (44°04′~ 46°40′ N, 125°42′~ 130°10′ E) (Fig. 1a, b, and Table 1). These two were introduced to the open experimental nursery of the Harbin Academy of Agricultural Sciences in Songbei District, Harbin, Heilongjiang Province, China. Daily temperatures during the period from November 2022 to January 2023 are all shown in Fig. 1c–e, and the temperature data referenced in this analysis are accessible at https://tianqi.2345.com/. The buds of two cultivars in the underground frozen soil layer were sampled on November 10, 2022, December 10, 2022, and January 10, 2023. Three biological replicates per cultivar were collected at each sampling stage. All samples were immediately frozen in liquid nitrogen and stored at −80 °C for subsequent experiments.

Figure 1. 

Experimental materials: (a) 'Hang Baishao' (HB). (b) 'Hongling Chijin' (HL). Temperature statistics for Songbei district, Harbin, Heilongjiang Province, China. (c) November 2022. (d) December 2022. (e) January 2023.

Morphological observations for budbreak and growth after an extremely cold winter in Harbin

Plant height, stem width, bud height, and bud diameter were collected using standardized measuring tapes. Plant height was measured vertically from the substrate to the tallest leaf apex. Stem width was assessed at the main stem base, and bud dimensions were recorded at their median points to ensure representative sampling. Based on previous studies on the comprehensive adaptive evaluation of herbaceous peony, the flower number per plant, i.e., the number of fully bloomed flowers, was confirmed during the peak flowering stage^[25,26].

Physiological and biochemical assays

Soluble sugar contents were determined using a Soluble Sugar Content Assay Kit (Bixbio, Beijing, China). Soluble protein concentrations were measured with a Bradford Protein Assay Kit (Beyotime, Shanghai, China). Stress-related enzyme activities and proline (PRO) contents were quantified using specific detection kits: peroxidase (POD) Assay Kit, superoxide dismutase (SOD) Assay Kit, and Proline Assay Kit (all from Jiancheng, Nanjing, China). Phytohormone levels in buds were analyzed via the enzyme-linked immunosorbent assay (ELISA), targeting nine phytohormones across five categories: abscisic acid (ABA); auxins: indoleacetic acid (IAA) and indolepropionic acid (IPA); gibberellins (GA₃ and GA₄); cytokinins: zeatin riboside (ZR) and dihydrozeatin riboside (DHZR); brassinolide (BR); and methyl jasmonate (MeJA)^[27].

RNA extraction, cDNA synthesis, and quality assessment

Total RNA was extracted from both cultivars using the EASYspin Plus Complex Plant RNA Kit (RN53, Aidlab, Beijing, China). Reverse transcription to synthesize cDNA was performed with the PrimeScript™ RT Reagent Kit with gDNA Eraser (Perfect real-time) (RR047A, Takara, Japan). RNA purity and concentration were evaluated using a NanoDrop 2000 spectrophotometer (Thermo SCIENTIFIC, America), while RNA integrity was precisely assessed with a LabChip GX Touch HT system.

Transcriptome sequencing and annotation

The experimental samples, consisting of three biological replicates each of HB and HL collected during three winter periods (November 2022 to January 2023), were processed for transcriptome analysis. All samples were shipped to the Biomarker Technologies Corporation (Beijing, China) for the transcriptome library construction and sequencing using the Illumina NovaSeq 6000 platform with PE150 sequencing mode. Sequence assembly was conducted via Trinity (v 2013), yielding 55,498 assembled transcripts (Unigenes or genes). Sequences of these transcripts were aligned against the NCBI non-redundant protein database (NR), Swiss-Prot, Clusters of Orthologous Groups (COG), Eukaryotic Orthologous Groups (KOG), evolutionary genealogy of genes: Non-supervised Orthologous Groups (eggNOG 4.5), and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases by using DIAMOND software. KEGG Orthology (KO) annotations were obtained through KOBAS. Gene Ontology (GO) analysis was performed using InterProScan with integrated InterPro databases. Amino acid sequences of Unigenes were predicted and aligned against the Pfam database via HMMER software to acquire functional annotations.

qRT-PCR and detecting the accuracy of transcriptome Unigene expression levels

To verify the expression trend of transcripts in the transcriptome and investigate expression patterns of selected genes, this study designed primers targeting transcript sequences of candidate genes (Table 2) and performed quantitative expression assays following omics data acquisition. Quantitative real-time PCR (qRT-PCR) was conducted using the NovoStart® SYBR qPCR SuperMix Plus (GenStar BioSolutions, Shanghai, China) on a Bio-Rad CFX96 system (Bio-Rad Laboratories, Shanghai, China). The thermal cycling protocol included: 95 °C for 1 min (initial denaturation); 40 cycles of 95 °C for 20 s (denaturation), 52 °C for 20 s (annealing), and 72 °C for 30 s (extension); followed by a melt curve analysis (95 °C for 5 s, 65 °C for 5 s, and 95 °C for 5 s). The α-Tubulin gene (ATUBA) was used as the reference gene, and relative gene expression levels were calculated using the 2^−ΔΔCᴛ method^[28,29].

Identification of differentially expressed genes

Differentially expressed genes (DEGs) were identified using DESeq2 with stringent thresholds: p_adjust value < 0.001 and |log₂ fold change| > 1. DEGs from different comparison groups were pooled via Venn analysis, followed by time-series clustering using the MFuzz algorithm implemented on the Bioinformatics Online Platform (www.bioinformatics.com.cn). KEGG functional enrichment analysis was performed using TBtools (v2.225)^[30].

Adopting the WGCNA to screen the key genes involved in cold tolerance

The WGCNA package in R was employed to construct co-expression networks by integrating morphological and physiological indices measured previously and transcriptomic data from bud samples. Gene modules were correlated with cold resistance-related traits through Pearson correlation coefficients. Key modules demonstrating significant correlations were subjected to KEGG functional enrichment analysis using TBtools (v2.225). Hub genes within these modules were identified and ranked by degree values (a connectivity metric) using Cytoscape (v3.10.1) for the co-expression network visualization^[31].

Data analysis

Data processing, statistical analysis, and graphical visualization were performed using Excel (2016), SPSS (v25.0), R (v4.4.1), TBtools (v2.225), Cytoscape (v3.10.1), GraphPad Prism (v10.0), the Bioinformatics Online Platform (www.bioinformatics.com.cn), and the Biomarker Cloud Platform (https://international.biocloud.net). One-way analysis of variance (ANOVA) was used for statistical comparisons, with results expressed as mean ± standard deviation (SD) from at least three biological replicates and statistical significance set at p < 0.05.

As herbaceous perennials, two cultivars exhibited distinct morphological differences. HL demonstrated greater morphological traits compared to those of HB (p < 0.001). Specifically, HL exhibited a plant height of 89.70 cm vs 73.24 cm for HB, a crown width of 81.40 cm vs 72.07 cm, 4.30 vs 3.47 stems per plant, and 12.67 vs 5.87 flowers per plant (Fig. 2a–d). HL displayed a more robust growth habit with larger overall plant dimensions and more abundant flowering than HB. Bud length differed significantly between cultivars throughout overwintering stages, with HL exhibiting consistently longer buds than HB (Fig. 2e). Bud diameter showed no significant differences, though HL buds were marginally thicker than HB at each stage (Fig. 2f).

Figure 2. 

Morphological and physiological data in 'Hang Baishao' (HB) and 'Hongling Chijin' (HL) during the winter of Harbin. (a) Plant height. (b) Crown width. (c) Number of stems per plant. (d) Number of flowers per plant. (e) Bud length. (f) Bud diameter. (g) Soluble sugar content. (h) Soluble protein content. (i) PRO content. (j) SOD activity. (k) POD activity. All data are the means of three replicates with standard deviations, and different letters indicate significant differences among the data according to Duncan's multiple range test. * indicates significant difference at the 0.05 probability level, ** indicates extremely significant difference at the 0.01 probability level, *** indicates extremely significant difference at the 0.001 probability level, **** indicates extremely significant difference at the 0.0001 probability level, and ns indicates non-significant. Different lowercase letters above the bars of Fig. 2g–k show significant differences at the 0.05 probability level; n = 3.

Soluble sugar and protein quantifications

Plants synthesize protective compounds such as soluble sugars and proteins to enhance cold-stress tolerance. Soluble sugar content increased significantly in both cultivars as cold stress intensified from November to January, with HL maintaining significantly higher levels than HB at all time points. For instance, in November, HL contained about 140 mg/g FW compared to 110 mg/g FW in HB; by January, levels rose to approximately 280 mg/g FW in HL vs 210 mg/g FW in HB (Fig. 2g). Soluble protein concentrations remained elevated between 6–7 mg/g FW throughout the sampling periods, with HB showing slightly higher values than HL at each stage (Fig. 2h).

Peroxidase activity assays

Physiological assessments of ROS accumulation and enzymatic scavenging systems revealed distinct patterns between HB and HL. Across all sampling stages, HB exhibited significantly higher PRO, SOD, and POD activities compared to HL (Fig. 2i–k). Throughout the sampling period from November to January, HB consistently showed higher levels of PRO, SOD, and POD activities compared to HL. PRO in HB peaked at approximately 870 μg/g FW in January, significantly above HL's 370 μg/g FW (Fig. 2i). Similarly, SOD activity in HB rose to about 4,200 U/g FW, surpassing the 1,880 U/g FW in HL (Fig. 2j). POD activity in HB exhibited a biphasic pattern, increasing to nearly 1.65 U/g FW in December before declining, while HL maintained stable levels around 0.1 U/g FW (Fig. 2k). Both cultivars exhibited rising PRO accumulation and SOD activity with prolonged cold exposure. Notably, POD activity in HB displayed a biphasic response, initially increasing, followed by a decline, contrasting with the sustained enzymatic profile in HL.

Phytohormone profiling

Dynamic phytohormone fluctuations can reflect plant adaptation to stress. During cold exposure, the eight phytohormones measured in the HB cultivar exhibited different response patterns. GA₃, IAA, IPA, and DHZR levels exhibited sustained increases throughout the experimental period (Fig. 3b, f, g, i). ABA, GA₄, MeJA, and ZR concentrations displayed biphasic patterns, with initial declines followed by subsequent increases (Fig. 3a, c, d, h). BR levels showed an initial rise followed by a reduction (Fig. 3e). Under identical experimental conditions, the same phytohormones in the HL cultivar exhibited different response patterns. ABA, GA₃, BR, IAA, IPA, ZR, and DHZR levels peaked mid-experiment before declining (Fig. 3a, b, e–i). GA₄ and MeJA concentrations progressively decreased over time (Fig. 3c, d).

Figure 3. 

Hormone content in three periods in 'Hang Baishao' (HB) and 'Hongling Chijin' (HL) under cold stress. (a) ABA content. (b) GA₃ content. (c) GA₄ content. (d) JA content. (e) BR content. (f) IAA content. (g) IPA content. (h) ZR content. (i) DHZR content. Different lowercase letters above the bars of Fig. 3a–i show significant differences at the 0.05 probability level; n = 3.

Comparative analysis shows that at all sampling stages, the concentrations of ABA, GA₃, GA₄, ZR, IPA, and DHZR were consistently higher in the HB cultivar than in the HL cultivar. Specifically, the ABA levels in HB were higher than in HL across all three stages, with concentrations in HB ranging between 150–200 ng/g FW, while those in HL remained between 100–150 ng/g FW. A notable exception was observed in November 2022, when the HL cultivar exhibited elevated levels of MeJA, BR, and IAA compared to the HB cultivar. Specifically, MeJA levels in HB showed an initial decrease followed by an increase as cold stress intensified, reaching approximately 36 ng/g FW by January, whereas in HL, they exhibited a declining trend. Both HB and HL showed an initial increase followed by a decrease in BR levels as cold stress progressed; the peak level in HB reached over 7 ng/g FW, while in HL it was only about 5 ng/g FW. However, the phytohormone concentrations in the other two stages were still higher in the HB cultivar (Fig. 3a–i).

Transcriptome sequencing overview

The length distributions of Unigene, CDS, and protein sequences collectively reflect key structural and functional features of the transcriptome. Most Unigenes are short, with the highest abundance in the 300–400 nt range, decreasing gradually as length increases, though a slight rebound occurs beyond 3,000 nt (Fig. 4a). Similarly, the CDS length distribution follows a declining trend as sequence length increases, with a steep drop in frequency observed among shorter sequences before leveling off into a slower descent. Toward the tail end, there is a slight uptick in CDS counts (Fig. 4b). Meanwhile, the protein sequence length distribution displays a distinct 'high left, low right' pattern, where short protein sequences vastly outnumber longer ones, and their frequency drops sharply as length increases (Fig. 4c).

Figure 4. 

Statistics and functional analysis of DEGs in HB and HL transcriptome during the winter of Harbin. (a) Length distribution of Unigene. (b) Length distribution of complete CDS. (c) Length distribution of complete protein. (d) Number of DEGs in three periods of HB and HL under cold stress. (e) Venn diagram of DEGs in three periods of HB and HL under cold stress. (f) KEGG enrichment of DEGs.

Statistical analysis of Unigene annotation

COG organizes prokaryotic orthologs into 26 categories, while KOG performs a similar classification for eukaryotes (4,852 groups) (Table 3). The expanded eggNOG database integrates these frameworks with additional viral proteins. Pfam employs HMM modeling to characterize protein domains, and GO standardizes gene function classification through three ontologies (Molecular Function, Cellular Component, Biological Process). Finally, KEGG enables systems-level analysis by mapping gene products to biological pathways and networks (Table 3).

Identification and KEGG enrichment of DEGs

DEGs between two cultivars at various stages were screened from the transcriptome and subjected to the Venn analysis. The November sampling period showed the highest downregulation number among all stages, with 2,296 upregulated vs 3,798 downregulated genes identified from the total 6,094 DEGs detected between HB and HL cultivars. While in January, 6,112 DEGs were identified between HB and HL, including 3,071 upregulated and 3,041 downregulated genes, representing the highest upregulated number among all stages (Fig. 4d). Cross-comparison across three stages revealed 1,558 overlapping DEGs (Fig. 4e). A total of 10,700 DEGs (differentially expressed genes) were further analyzed, while stage-specific DEGs proved more numerous in January (2,805) and November (2,655), significantly exceeding December's count (1,067). KEGG enrichment analysis demonstrates significant enrichment in pathways including the starch and sucrose metabolism, Plant-pathogen interaction, Plant hormone signal transduction, MAPK signaling pathway-plant, and Zeatin biosynthesis (Fig. 4f). These findings would appear to suggest that herbaceous peony responds to cold stress through coordinated regulation of multiple biological processes and metabolic pathways.

Validation of several key DEGs' expression

Based on the functional annotation, nine candidate genes associated with sugar metabolism (Fig. 5a–c), ROS regulation (Fig. 5d–f), and phytohormone signaling (Fig. 5g–i), were selected for the validation of expression by qRT-PCR. The concordance between qPCR-derived expression patterns and transcriptomic expression profiles (Fragments Per Kilobase of exon per Million mapped reads, FPKM) confirmed the reliability of sequencing results. These nine candidate genes are closely associated with some important physiological regulation under extremely cold stress, so they will be subjected to further in-depth analysis in subsequent discussions.

Figure 5. 

The qRT-PCR analysis and transcriptomic expression of cold-tolerance related genes in 'Hang Baishao' (HB) and 'Hongling Chijin' (HL). (a) APL3. (b) BMY3. (c) GBSS1. (d) GSH2. (e) GRXS17. (f) MSD1. (g) SPY. (h) ABF3. (i) IAA9. The bar chart represents the gene relative expression level, and the line chart represents the transcriptome expression level. FPKM: Fragments per kilobase of exon model per million mapped reads. Different lowercase letters above the bars of Fig. 5a–i show significant differences at the 0.05 probability level; n = 3.

As cold stress intensified, the expression levels of APL3, BMY3, GBSS1, and GRXS17 consistently increased in both HB and HL, with significantly higher expression observed in HL compared to HB, indicating their positive roles in responding to cold stress (Fig. 5a–c, e). In contrast, although MSD1 expression rose with decreasing temperature, its levels were significantly higher in HB than in HL. While IAA9 expression declined with temperature, it remained significantly higher in HL than in HB; meanwhile, the expression profiles of GSH2, SPY, and ABF3 displayed irregular and non-directional shifts across temperature regimes, showing no clear correlation with stress intensity or cultivar-specific resilience (Fig. 5d, f–i). These divergent and often contradictory patterns suggest that the functional contributions of MSD1, IAA9, GSH2, SPY, and ABF3 to cold stress adaptation are likely complex and context-dependent, warranting further mechanistic investigation (Table 4).

Clustering of DEGs and analysis of enriched KEGG pathways

To further explore the expression patterns of DEGs in HB and HL during the three periods, all DEGs were divided into four clusters, presenting four trends (Fig. 6): with the temperature drop, Cluster 1 was highly expressed in HB, where it increased, but showed no significant change in HL. Cluster 2 increased and was highly expressed in HL, but showed no significant change in HB. Cluster 3 decreased in both HB and HL, but was highly expressed in HB. Cluster 4 decreased in both HB and HL, but was highly expressed in HL.

Figure 6. 

Clustering and KEGG enrichment analysis of DEGs of HB and HL during the winter of Harbin; The x-axis labels represent cultivar names combined with sampling months, using the following nomenclature: HBN denotes the HB cultivar sampled in November (N), HBD indicates the HB cultivar in December (D), and HBJ corresponds to the HB cultivar in January (J). Similarly, HLN, HLD, and HLJ refer to the HL cultivar sampled in November, December, and January, respectively.

The KEGG analysis of each cluster shows Cluster 1 is significantly enriched in plant-pathogen interaction, environmental adaptation, protein families: signaling and cellular processes; Cluster 2 is significantly enriched in the environmental information process pathway; Cluster 3 is significantly enriched in the carbohydrate metabolism, metabolism, and phenylpropanoid biosynthesis pathways; Cluster 4 is significantly enriched in the DNA replication proteins and biosynthesis of other secondary metabolism pathways (Fig. 6).

Screening the key modules involved in herbaceous peony's response to cold stress based on the WGCNA

To further determine the specific modules and genes related to cold resistance-related traits between HB and HL, a WGCNA analysis was conducted. This analysis combined transcriptome sequencing data with 14 physiological indicators. These indicators, which included levels of soluble sugar, antioxidant enzyme activity, and phytohormone contents, were obtained from previous observations (Figs 2g–k, 3). Results show that 21 different modules were obtained. Analysis of the relationship between the module characteristic genes and physiological parameters related to cold resistance indicate that ME black is significantly positively correlated with POD, SOD, PRO, ABA, GA3, ZR, GA4, and IPA; while ME greenyellow is significantly negatively correlated with solution protein, POD, SOD, PRO, ABA, IAA, GA₃, BR, DHZR, and IPA (Fig. 7).

Figure 7. 

WGCNA analysis. (a) Cluster dendrogram. (b) Module-trait relationships.

Expression of key genes in two modules and genes related to the aforementioned physiological and biochemical indicators

Figure 8 highlights some key candidate genes selected from two modules for further investigation. Ten genes with the highest average expression levels were selected from the Me black and Me green-yellow modules, respectively (Fig. 8). These genes were annotated as XERO1, GSTU9, LEA14, NADP-ME4, TRA2, etc., which may play key roles in regulating the cold resistance of herbaceous peony cultivars (Table 5).

Figure 8. 

Transcriptome expression levels of genes potentially involved in cold stress response are shown for 'Hang Baishao' (HB) and 'Hongling Chijin' (HL). (a)–(j) Genes from Me black modules, and (k)–(t) Me greenyellow modules are displayed in the left and right panels, respectively. Different lowercase letters above the bars of (a)–(t) show significant differences at the 0.05 probability level; n = 3.

Notably, highly expressed genes, including XERO1 (a dehydrin gene), MT3 (limiting oxidative damage), and LTI65 (a low-temperature-responsive protein gene), are all functionally associated with plant stress responses, particularly cold acclimation, though exhibiting distinct expression dynamics. The sustained high expression of XERO1 is consistent with a role in basal protection, likely contributing to cellular stability under cold stress. Whereas the induction of LTI65 likely represents a targeted response to cold-induced dehydration. While the expression of XERO1 and LTI65 was elevated in HL compared to HB, implicating them as potential contributors to its strong cold tolerance, their expression patterns did not demonstrate a perfectly synchronized upward trend correlating with progressive decreases in temperature.

NADP-ME4 enhances the scavenging of active oxygen species by upregulating key antioxidant genes. In this study, NADP-ME4 expression was consistently higher in HB than HL across all three sampling stages, suggesting that HB experienced greater oxidative stress under low temperatures, thereby triggering a stronger compensatory antioxidant response^[32].

A similar expression pattern was observed for MT3, which was upregulated in both cultivars in response to increasingly low temperatures, with particularly high transcript levels observed in the buds of HB in January, indicating its close involvement in the peony's adaptation to cold stress. Contrary to expectations, however, its expression was significantly higher in the less cold-tolerant cultivar HB than in the more resilient HL. It may suggest its intense activation in the sensitive cultivar is a response to heightened oxidative damage. Given these intriguing expression profiles, further functional studies are warranted to elucidate the specific mechanisms underlying their roles in cold response.

Identifying crucial genes in the cold stress-related modules

To further investigate the cold tolerance-related genes within the two modules under both HB and HL conditions, a broadly targeted screening of the transcripts was performed. This screening focused on those associated with sugar and phytohormone signaling pathways to refine the analysis. These selected transcripts were then imported into Cytoscape to construct co-expression networks, where the degree value served as the primary metric for mapping interactions (Fig. 9). The top 15 genes by degree in each module were labeled to facilitate subsequent analysis and discussion (Table 6 and 7).

Figure 9. 

Co-expression networks of important genes in the (a) ME black, and (b) ME yellowgreen modules. In the network visualization, both the size and color intensity of each gene node are determined by its degree value. Larger node sizes and darker colors represent higher degree values.

Notably, in the ME black module, ANTHRANILATE SYNTHASE 2 (ASA2), a gene strongly associated with phytohormone signaling and involved in indole-3-acetic acid synthesis, stands out. Its central position suggests it may coordinate developmental adjustments under cold stress. Another significant gene in this module is CHLOROPLAST-LOCALIZED SEC14-LIKE PROTEIN 1 (CPSFL1), which also shows strong connections to hormone signaling and functions as a negative regulator in the ABA signaling pathway, directly participating in hormonal response mechanisms. Two other genes with high degree values in this module were AT5G55530 and AT1G66860, both of which lack annotation in The Arabidopsis Information Resource; however, their prominent network positions highlight them as promising candidates for future functional characterization. Additionally, DOF2.1, a DOF transcription factor containing a conserved zinc finger DNA-binding domain, was identified as another high-degree hub, suggesting its potential role as a key transcriptional regulator in the cold stress response.

Within the ME yellowgreen module, the gene with the highest degree value is WRKY DNA-BINDING PROTEIN 39 (WRKY39). This gene is a member of the WRKY transcription factor family and is a well-studied transcription factor gene known for its role in plant stress resistance. The CHLOROPLAST GRPE 2 (CGE2) gene also ranked high in degree, suggesting its importance in chloroplast protein import and plastid-to-nucleus signaling during cold adaptation. The emergence of these highly connected hub genes suggests that cold tolerance in herbaceous peony may be coordinated through the integration of multiple synergistic mechanisms.

The process of plant resistance to low-temperature stress is dually regulated by internal signals and external environmental changes. The generation and transmission of internal signals, as well as their regulation of cold resistance, are influenced by various elements such as sugars, proline, calcium ions, and protective proteins^[6]. Soluble sugars are the primary substances for plant cold resistance. It functions to raise cellular osmotic potential and improve water retention, while also acting as a vital winter nutrient reserve that supports other physiological and biochemical processes essential for coping with low-temperature stress. The findings in iris (Iris japonica) demonstrate a dual role for sugar molecules as both signals and energy reserves in regulating growth transition through IjFLCL1-mediated pathways^[33]. In addition, the content and conversion rate of starch during overwintering are closely associated with cold hardiness^[34].

Meanwhile, cold treatment significantly increased soluble sugar content in both sugarcane (Saccharum officinarum L.) varieties (ROC22 and GT42), with GT42 showing greater accumulation and higher cold resistance than ROC22 at the bud stage^[35]. Furthermore, salicylic acid treatment enhanced cold hardiness in Magnolia wufengensis by increasing soluble sugar accumulation, which improves freezing tolerance during natural cold acclimation^[36]. In this study, HL demonstrated significantly superior morphological characteristics. These included longer bud length, greater plant height, larger crown width, increased stem number per plant, and higher flower count per individual. HL also has significantly higher soluble sugar content compared to HB, showing its stronger cold resistance, demonstrating the critical role of soluble sugars in mediating plant adaptive responses to cold stress (Fig. 2a–e, g). These results are consistent with previous studies conducted in sugarcane and Magnolia wufengensis^[35,36]. The enriched starch and sucrose metabolism pathways likely facilitated increased soluble sugar levels through enhanced starch hydrolysis, thereby enhancing cold stress resistance in herbaceous peony.

Some studies suggest a positive correlation between increased protein content and plant cold resistance. A study demonstrated a positive correlation between increased soluble protein content and enhanced cold resistance in annual pitaya branches during 7 d cold stress^[37]. Similarly, HB and HL exhibited an upward trend in soluble protein content under escalating cold stress, mirroring the response observed in annual pitaya shoots. However, this positive correlation does not hold consistently across all plant species and developmental stages. For example, the wheat of cold-sustainable genotypes maintained consistently high soluble protein levels both under optimal conditions and after cold exposure^[38]. Also in this study, soluble protein content remained consistently high in both cultivars across stages, with HB slightly higher than HL (Fig. 2h). The results suggest that soluble proteins likely contributed to enhanced cold tolerance in both cultivars, but since no significant difference was observed between them, this parameter may not serve as a reliable indicator for evaluating and determining herbaceous peony cold resistance. The soluble protein content remained relatively stable in herbaceous peony, potentially due to increased synthesis of stress-resistant proteins and reduced production of other protein types through the MAPK signaling pathway.

Metabolism of antioxidant enzymes closely regulates the cold resistance of herbaceous peony

Plants accumulate excessive free radicals and ROS that damage membrane structures under various stresses^[4]. SOD and POD are primary antioxidant enzymes, while PRO acts both as a ROS scavenger and an osmotic regulator to prevent cellular dehydration. In this study, both HB and HL showed increased SOD/POD activities and PRO content under cold stress. However, HB demonstrates significantly higher SOD/POD activities and greater PRO content than HL across all stages (Fig. 2i–k). Interestingly, studies on banana (Musa spp.) revealed that while PRO content increased in all cultivars during chilling, cold-tolerant banana varieties accumulated less PRO than cold-sensitive ones^[39]. This observation shows a certain similarity with the findings of the present study. This similarity may suggest that HB, with lower cold tolerance compared to HL, experiences more pronounced oxidative damage.

Regulation of traditional phytohormones in cold tolerance of herbaceous peony

The five major traditional phytohormones, namely ABA, auxin, gibberellin, cytokinin, and ethylene, are critically involved in plant resistance to cold stress and other environmental challenges. ABA, a core phytohormone in cold stress, activates cold acclimation and enhances key stress-mitigating components, including proline, soluble proteins, and SOD activity^[40,41]. Exogenous application of ABA has been shown to significantly enhance cold tolerance in cucumber (Cucumis sativus L.)^[41]. Similarly, studies in Bermudagrass (Cynodon dactylon (L). Pers.) demonstrate that ABA treatment effectively strengthens cold stress tolerance, particularly in cold-resistant cultivars^[42]. Moreover, ABA has been found to interact with other phytohormones, such as strigolactones and brassinosteroids, to positively regulate cold tolerance in tomato^[43]. ABA content remained consistently high in both HB and HL cultivars, suggesting its critical role in the low-temperature stress response of herbaceous peony (Fig. 3a).

Auxins play a key role in cold stress responses. For example, overexpression of microRNA393 improved cold tolerance and tillering of switchgrass (Panicum virgatum L.) through regulation of auxin signaling transduction^[44]. Building on prior evidence of auxin's role in cold adaptation, the analysis of HB and HL reveals that their hormonal dynamics also operate through auxin pathways. Specifically, this study revealed that IAA, IPA, and DHZR showed sustained increases in HB, but followed an initial rise, followed by a decline in HL, consistent with cold adaptation traits (Fig. 3f, g, i). This implies that both cultivars responded to cold stress through substantial increases in cytokinins and auxins, with HL exhibiting stronger cold resistance by ultimately reducing these phytohormones after initial elevation.

Meanwhile, in HL, GA₃ first increased then decreased, whereas in HB, GA₃ consistently rose and was significantly higher than in HL (Fig. 3b). This confirms that HL enhances cold resistance by suppressing GA signaling, representing a strategic balance between growth and stress adaptation.

Regulation of novel phytohormones in cold tolerance of herbaceous peony

In addition to the traditional five major classes of phytohormones, the scientific community has identified a variety of novel phytohormones over recent decades, such as JA, BR, salicylic acid, and strigolactones, all of which have been demonstrated to play significant roles in enhancing plant resistance to cold and other environmental stresses. Furthermore, exogenous MeJA significantly promoted cold tolerance of tomato^[45]. Experimental studies demonstrate that exogenous jasmonates (JAs) enhance cold tolerance in blood orange (Citrus sinensis L. Osbeck)^[46]. This data showed JA content continuously decreased in HL but exhibited an initial decline followed by an increase in HB, suggesting JA plays a prominent role in the early cold stress phase before diminishing over time in HL (Fig. 3d).

BR, a natural stress mitigator, exhibits dual functions in growth regulation and abiotic stress protection. For instance, BR activates the ascorbate-glutathione cycle through an ABA-independent regulatory system, scavenging excess ROS under cold stress in grapevine (Vitis vinifera L.), thereby enhancing their cold resistance^[47]. Sodium nitrophenate plays a role in mediating the BR signal to regulate the cold tolerance of cucumber (Cucumis sativus L.)^[48]. BR content in both HB and HL initially increased, then declined, peaking during early cold stress before decreasing temporally (Fig. 3e). This response pattern is not unique; under cold stress, BR levels in grapevine, HB, and HL all exhibited a similar trend: an initial increase to a peak, followed by a subsequent decline.

Regulation of cold adaptation by key genes involved in sugar metabolism, ROS regulation, and phytohormone signaling

Plants employ diverse molecular strategies to enhance cold tolerance, involving key genes that regulate sugar metabolism, redox homeostasis, and phytohormone signaling. Specifically, APL3 in Spirodela polyrhiza promotes starch accumulation to support turion formation and energy storage^[49]. While BMY3 in rhubarb (Rheum rhabarbarum L. Polygonacea family) facilitates starch-to-sugar conversion to enable dormancy release and cold adaptation^[50]. GSH2 enhances chilling tolerance in species like tomato by synthesizing glutathione to mitigate oxidative stress, dependent on energy availability and regulatory cascades. Similarly, GRXs in tea plant (Camellia sinensis (L.) Kuntze) and banana (Musa acuminata) maintain redox balance under low-temperature stress, with members like MaGRX17 specializing in antioxidative defense^[51,52]. APL3, BMY3, and GRXS17 all exhibited increased expression levels in both HL and HB under progressively severe cold stress, with even higher expression observed in the more cold-tolerant cultivar HL (Fig. 5a, b, e). This indicates their active roles in responding to cold stress, consistent with the functions of their orthologs in Spirodela polyrhiza, rhubarb, banana, tomato, rice, and tea plant^[49−52].

In contrast, MSD1 is down-regulated across Arabidopsis ecotypes under cold conditions, suggesting context-specific roles^[53]. Meanwhile, ABF3 in pear (Pyrus pyrifolia White Pear Group) integrates ABA and GA signaling to ensure winter survival and synchronized bud-break^[54]. In Moso bamboo (Phyllostachys edulis [Carrière] J. Houz.), IAA9 promotes stress adaptation via hypomethylation, which facilitates the integration of auxin signaling^[55]. These mechanisms collectively underscore the complexity and species-specificity of plant cold adaptation; however, the coordinated regulatory functions of genes in response to cold stress require further elucidation.

Regulatory functions of key cold-responsive genes in herbaceous peony

The expression levels of XERO1 and LTI65 were higher in HL compared to HB, suggesting their potential functional importance in cold tolerance, though their precise mechanisms remain to be fully elucidated. Similarly, MT3 enhances cold tolerance in banana (Musa acuminata Colla cv. Cavendish) through sustained antioxidant defense that mitigates oxidative damage under low temperature^[56]. NADP-ME4 enhances the expression of genes encoding superoxide dismutase, peroxidase, and pyrroline-5-carboxylate synthase, thereby boosting the plant's ability to scavenge ROS. This activity reduces the damage caused by cold stress^[32]. The expression patterns of genes such as XERO1, LTI65, and MT3 in HL and HB did not fully align with initial expectations, suggesting that their roles in responding to cold stress require further investigation (Fig. 8a, b, i).

By regulating the expression of target genes, WRKY transcription factors are increasingly recognized by emerging research as critical players in multiple facets of plant biology, from growth and development to adaptation to environmental stresses^[57]. For example, a comparative genomic analysis of three Pterocarya species (specifically between P. stenoptera and P. macroptera), WRKY39 was identified as a key cold-tolerance gene under environmental selection during the speciation of Pterocarya trees^[58]. The specific functions of WRKY transcription factors in cold tolerance of herbaceous peony remain unclear; however, given the likely conservation of their regulatory pathways, a key objective is to identify how peony WRKYs integrate upstream cold-stress and hormone signaling to enhance cold tolerance. PCA and network analyses demonstrated ASA2's effectiveness in enhancing cold tolerance by strengthening photosynthesis-antioxidant coordination. This treatment upregulated key genes, reducing oxidative damage and preserving chloroplast function in common bean (Phaseolus vulgaris L.)^[59]. The function of these genes in the cold stress response of herbaceous plants is yet to be explored, although their roles in cold tolerance have been preliminarily established in other species. In Arabidopsis thaliana, CPSFL1 potentially facilitates vesicle formation by transferring lipids such as PIP and PtdOH, which are critical for thylakoid biogenesis and carotenoid transport. It is hypothesized that by maintaining the integrity of these photosynthetic components, CPSFL1 may indirectly contribute to the plant's cold tolerance^[60]. Future investigations could therefore focus on the synergistic relationships among these transcription factors and key genes to identify critical regulatory nodes within the cold tolerance network of herbaceous peony.

Physiological differences are consistent with HL's superior overwintering phenotypes compared to HB. Transcriptomic analysis revealed that differential gene expression in starch and sucrose metabolism, plant-pathogen interactions, and plant hormone signal transduction pathways in cold adaptation. Notably, genes such as XERO1 (a dehydrin gene), LTI65 (a low-temperature-responsive protein gene), NADP-ME4 (involved in energy metabolism), PLA2A (a lipid acyl hydrolase), WRKY39 (a member of the WRKY transcription factor), and ASA2 (involved in plant hormone biosynthesis) were identified as potential key regulators for further functional investigation. The innovation lies in evaluating herbaceous peonies in the extreme cold of a major metropolitan area. Unlike controlled laboratory settings, this approach provides a more authentic assessment of frost damage. It also screened southern against local cultivars to identify cold-resilient varieties. This study establishes a foundation for elucidating the physiological mechanisms of cold stress adaptation in herbaceous peony and provides essential groundwork for future functional validation and interaction studies of cold tolerance-related genes.