Transcriptomic Profiling Analyses Revealed Candidate Genes Under Freezing Stress in Siberian Apricot (Prunus sibirica)

Abstract

Siberian apricot (Prunus sibirica) is an important ecological and commercial woody plant that is negatively affected by spring frosts. However, the mechanisms that control gene expression in adaptation to freezing remain largely unknown. In this work, we investigated the physiological, molecular, and phenotypic characteristics of pistils of two P. sibirica clones that differ in their ability to withstand freezing stress. A total of 14,717 unigenes categorized into 38 functional groups were identified. Additionally, the two P. sibirica clones included 3931 up-regulated and 2070 down-regulated differentially expressed genes (DEGs). Many DEGs are related to Ca²⁺ and MAPK signaling, carbohydrate biosynthesis and metabolism, plant hormone signal transduction, biosynthesis of amino acids, and photosynthesis. The metabolism of carbohydrates, amino acids, lipids, secondary metabolites, plant hormone signal transduction, and terpenoid metabolism were the transcriptome modifications most significantly altered by freezing stress. Real-time quantitative PCR (RT-qPCR) was used to verify the precision of the RNA-seq data. PsbHLH18, PsMYB4, PsMYB44, PsPOD1, and PsCDPK5 may play important roles in the freezing tolerance of the P. sibirica floral organ. This study provides a foundation for further studies on the complex mechanisms of freezing stress response in P. sibirica.

1. Introduction

Siberian apricot (Prunus sibirica L.), a deciduous fruit tree belonging to the Rosaceae family, is found in many different parts of the world, including northern and northeastern China, eastern and southeastern Mongolia, eastern Siberia, and Russian marine regions [1]. P. sibirica can thrive in unfavorable environments with limited rainfall, barren soil, chilly temperatures, and sandy breezes. In semi-arid and dry locations, P. sibirica is the most common and important pioneer forest species [2]. In addition, the seed kernel of P. sibirica, with an oil content of 40%~56%, is an important woody oil plant recognized as a new and promising feedstock for biodiesel [3]. P. sibirica is a typical early-flowering tree species with important ornamental, decorative, and cultural value as a traditional flower species in China [4]. However, this can also make the floral organ more susceptible to spring frosts [5]. Low temperatures have a major negative impact on the flowering and fruiting of P. sibirica, leading to huge yield losses. Nowadays, spring frosts are a serious obstacle to the industrialization of P. sibirica [2,5].

Low temperature is one of the major abiotic stresses that pose a threat to plant growth and development and affect the sustainable development of agroforestry [6,7]. The two types of low-temperature or cold stress are chilling (0~15 °C) and freezing stress (<0 °C) [8]. To withstand this unfavorable factor, plants possess an extensive range of regulatory systems, such as photosynthetic responses, antioxidant systems, cell membrane responses, alterations in metabolism and protein products, and so on [9]. Plants respond to adverse external conditions by converting external stimuli into internal signals and triggering defense mechanisms through transcriptional cascades [6].

There are many signaling pathways in plants to deal with low-temperature stress, such as Ca²⁺, phytohormones, and secondary metabolic pathways. These signaling pathways will stimulate the expression of relevant responsive transcription factors or functional proteins to enhance the cold tolerance [10]. These signaling pathways will induce the expression of relevant genes, which will cause alterations in physiological and biochemical indexes and thus improve cold resistance [8]. Stabilization of osmotic pressure in plants after low-temperature stress is essential to protect cell membranes. Low-molecular-weight solutes, such as soluble sugars, soluble proteins, flavonoids, and proline, can act as better permeability buffers to reduce the degree of damage [11]. Cold stress causes large amounts of reactive oxygen species (ROS) to be produced in the plant. Low concentrations of ROS stimulate the opening of the defense system, while high concentrations of ROS can cause irreversible damage to the plant. At this time, the levels of relevant protective enzymes such as superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) increase rapidly [8,12].

Among the molecular mechanisms of plant response to cold hardiness, the ICE (Inducer of CBF expression)-CBF (C-repeat-binding factor)-COR (cold-regulated genes) module is the most well recognized, important, and reported-on signaling pathway [13]. Various gene family members, like WRKY, bHLH, NAC, and MYB, are involved in this pathway to regulate the expression of these three genes after cold stress [14]. In addition, the mitogen-activated protein kinase (MAPK) cascade is also a major signaling transduction pathway in response to cold stress [15]. MAPK can induce a series of phosphorylation cascade reactions to enhance immunity, increase tolerance, and thus regulate low-temperature resistance. In contrast, cold-responsive protein kinase 1 (CRPK1) can inhibit CBF expression under freezing conditions by phosphorylating the associated proteins for nuclear translocation, thereby subjecting plants to injury [16]. At present, although the genetic and molecular mechanisms behind plants’ reactions to cold stress have been identified and revealed, there are many new signaling pathways to be further explored and excavated. In particular, woody plants have longer growth cycles and more complex cold resistance mechanisms than annual plants.

Transcriptional regulation is an important biological process in plants [17]. Transcriptome sequencing (RNA sequencing) technology allows for the rapid mining of key transcription factors and functional genes, as well as the construction of gene regulatory networks [18]. The transcriptome has been extensively used in woody plants to investigate low-temperature-stress response genes [19]. In Camellia sinensis, a large number of cold-stress differentially expressed genes (DEGs) have been unearthed by comparing cold-tolerant and cold-susceptible cultivars, which can be involved in the CBF-COR pathway to improve cold resistance [20]. Transcriptome sequencing of three cold-tolerant and -susceptible cultivars of mango (Mangifera indica) under low temperature stress identified 22,536 DEGs. These DEGs were mainly involved in phytohormone regulation, MAPK, and ICE-CBF-COR signaling cascades. In addition, they were significantly enriched in the pathways of photosynthesis, flavonoid, and alkaloid biosynthesis [21]. A total of 1544 DEGs were identified in pistils by transcriptome comparison of two apricot (Prunus armeniaca) cultivars differing in cold resistance under spring frost [10]. These genes belonged to the gene families of NAC, WRKY, AP2/ERF, and zinc finger proteins, which were mainly involved in the Ca²⁺ and MAPK signaling pathways. It is noteworthy that plants gradually improve freezing tolerance mainly by inducing cold-responsive genes, which are both influenced by the timing of cold stress induction and are closely related to the range of temperature changes [22].

To gain new insights into freezing-related pathways and genes, the two clones (cold-tolerant ‘CR338’ and cold-sensitive ‘NC28’) demonstrating varying cold sensitivity were investigated via morphological and physiological analyses during the flowering stage of Siberian apricot, and their traits of freezing injury rates were compared following treatment at −4 °C. The transcriptome alterations under low-temperature treatment were analyzed using RNA-seq. The study reveals a potential molecular regulatory pathway involved in freezing stress and predicts the freezing-related candidate genes of P. sibirica, providing a framework for subsequent functional validation. Meanwhile, the findings offer current information and genetic resources for improving the freezing tolerance of P. sibirica.

2. Materials and Methods

2.1. Plant Materials and Treatment

The two seven-year-old P. sibirica clones (‘CR338’ and ‘NC28’) examined in this study were cultivated in the National Long-term Research Base for Apricot Germplasm Conservation and Breeding of Liaoning (longitude 119°32′49″ E, latitude 41°17′27″ N, Kazuo, Liaoning, China). The provenance of the clones ‘CR338’ and ‘NC28’ were from Zalantun, Hulunbeier, Inner Mongolia (longitude 122°12′31″ E, latitude 47°17′7″ N), and Aohan, Chifeng, Inner Mongolia (longitude 119°51′45″ E, latitude 42°17′15″ N), China, respectively. The artificial frost chamber was placed on whole trees of two P. sibirica clones in full bloom. The clone pistil of P. sibirica with a less than 5% freezing-injury rate are cold-tolerant, and those with more than 20% are cold-sensitive. The frost box was set at four temperature gradients of 6 °C (control), −2 °C, −4 °C, and −6 °C for 1 h, and then pistils were collected. Some of the pistils were examined regarding their morphological observations and physiological parameters. Other collected pistils (−4 °C for 1 h) were prepared for an RNA-seq analysis. Three biological replicates were performed using pistils obtained from each clone of the tree for the subsequent analysis. After returning to the lab, the samples were kept in an ultra-low-temperature freezer (−80 °C) to extract RNA. To represent the experimental group, CR338_1, CR338_2, and CR338_3 (cold-tolerant) were used, and to represent the control group, NC28_1, NC28_2, and NC28_3 (cold-sensitive) were used.

2.2. The Extent of the Pistil Damage Under Freezing Stress

Due to the release of phenolic compounds, the tissues injured by frostbite turn brown-black in color [23]. The pistils obtained at each sampling were incubated at 18 °C following a one-hour exposure to cold stress (−4 °C). To determine the freezing-injury rate, three duplicates of 300 pistils were examined for every clone.

2.3. Physiological Analysis

The primary enzymatic scavengers of hydrogen peroxide (H₂O₂) are SOD and POD, both of which play an important role in this system that scavenges ROS in plants [8]. The activities of SOD and POD were determined by nitroblue tetrazolium (NBT) and guaiacol, respectively [2]. POD may oxidize guaiacol in the presence of H₂O₂ to create a brown material that has a maximum absorbance of 470 nm. The enzyme extract was made by mixing 2 mL of extract with around 0.2 g of frozen pistils in each duplicate, which was then centrifuged at 8000× g for 8 min at 4 °C.

Soluble sugar (SS) concentration was ascertained using the anthrone colorimetric technique [24]. Soluble protein (SP) concentration was determined by the Coomassie Brilliant Blue G-250 staining approach [25].

Malondialdehyde (MDA) content was determined using the thiobarbituric acid reactive substances method [2]. Absorbances were measured at 450, 532, and 600 nm to reduce the possibility of interference from other compounds. A conductivity meter was used to measure the cell membrane permeability (CMP) [2]. The CMP was initially measured with a conductivity meter after 0.3 g of pistils was thoroughly mixed with distilled water for four hours. After that, it was immersed in boiling water for 30 min. After resting for four hours, the data were collected again using a conductivity meter.

2.4. RNA Extraction, cDNA Library Construction, and Sequencing

Transcriptome sequencing was performed using three biological replicates of pistils from two P. sibirica clones. Approximately 1.0 g of pistils were thoroughly ground in liquid nitrogen, followed by the addition of Trizol reagent, and total plant RNA was extracted using the RNA Extraction Kit (Tiangen Biotech Co., Ltd., Beijing, China). Briefly, fresh pistils were thoroughly ground and Buffer SL was added; the mixture was centrifuged, and the supernatant was transferred to filter column CS set. After centrifugation, anhydrous ethanol was added to the supernatant and the mixture was transferred together to column CR3 set. After centrifugation, Buffer RW1, DNase I, 1, and Buffer RW were added, respectively. Finally, RNA was obtained by dissolving it with 70 μL of RNase-Free ddH₂O. The integrity of RNA was assessed using the RNA Nano 6000 Assay Kit for the Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA). The TruSeq PE Cluster Kit v3-cBot-HS (Illumina Inc., California, CA, USA) was used to sequence the cDNA libraries, producing 150 bp paired-end reads through the use of sequencing-by-synthesis technology. Using the P. armeniaca genome as a reference genome, sequence alignment and follow-up analyses were carried out [26]. Hisat2 v2.0.5 tools for mapping RNA-Seq reads to the reference genome. The transcriptome sequencing and assembly were performed by Novogene Technology Co., Ltd. (Beijing, China).

2.5. Transcriptome Assembly

Illumina CASAVA 1.8 CASAVA was used to convert the raw picture data from sequencing into sequence reads, and R language was used to filter the raw data. After removing adaptor sequences, reads with more than 10% ambiguous bases, reads with unknown nucleotides, low-quality reads (quality score < 30), and clean reads were produced. The Trinity software v2.15.2 reassembles the reads into complete transcripts, generating a collection of transcripts for further analysis and annotation. At the same time, the expression level of each transcript is calculated based on the reads coverage information in the RNA-Seq data. By using the Corset hierarchy for clustering, redundancy was eliminated. The constructed unigenes functioned as the pistil transcriptome’s reference sequences. The readcount of genes was normalized using fragments per kilobase of exon per million mapped fragments (FPKM).

2.6. Identification of Differentially Expressed Genes (DEGs)

Identification of DEGs were carried out using the DESeq program (1.12.0). It was used to construct a DataSet object to normalize the input matrix data. The resulting p-values were modified using Benjamini and Hochberg’s method for regulating the false discovery rate (FDR). Differences in expression were defined as significant DEGs with FDR < 0.05 and log2 (fold change) > 1.

2.7. GO and KEGG Enrichment Analysis

Implementation of Gene Ontology (GO) enrichment analysis was performed using the clusterProfiler (3.4.4) R package. DEGs were thought to significantly enrich GO terms if the corrected p-value was less than 0.05. The Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database tool for comprehending the higher-level operations and features of the biological system (http://www.genome.jp/kegg/, accessed on 21 March 2024). To test the enrichment analysis of DEGs in KEGG pathways, the clusterProfiler (3.4.4) R package was utilized.

2.8. Reverse Transcriptase–Quantitative PCR (RT-qPCR)

Twenty genes were chosen for an RT-qPCR investigation to assess the validity of the transcriptome sequencing results. The 18SrRNA gene was used as an internal reference gene [5]. The primer sequences are shown in Table S1. Gene-specific primers were designed using Primer Premier 5.0. To synthesize first-strand cDNA, approximately 1.0 μg of total RNA was reverse-transcribed using the EasyScript One-Step gDNA Removal and cDNA Synthesis SuperMix (TransGen, Beijing, China). The RT-qPCR evaluation was performed with the StepOne Real-Time PCR System using SuperReal PreMix Plus SYBR Green (TransGen, Beijing, China). The PCR reaction program consisted of the following steps: an initial denaturation at 95 °C for 30 s, followed by 40 cycles of 95 °C for 10 s, 60 °C for 15 s, and 60 °C for 1 min, finishing with 95 °C for 15 s. The relative levels of gene expression were determined using the 2^−ΔΔCT method [27]. Three independent biological replicates were used in this analysis. The SPSS v26.0 program was used to conduct a statistical analysis. The GraphPad Prism v8.4.0 tools were used to make bar graphs.

3. Results

3.1. Morphological of P. sibirica Floral Organ to Freezing Stress

Damaged apricot tissues leak and generate phenolic compounds under freezing conditions, giving them a brown-black appearance [10]. In the absence of adversity stress, there were no apparent changes between the two P. sibirica clones in this investigation, based on an examination of the damaged floral organs (Figure 1a). However, there were some significant morphological differences between ‘CR338’ and ‘NC28’ for the floral organ and pistils collected after freezing stress. ‘NC28’ browned in the petals and pistils more severely than ‘CR338’ (Figure 1a). The freezing-injury rates of ‘CR338’ and ‘NC28’ pistils were 3.79% and 35.33%, respectively (Figure 1b). The freezing-injury rates of the ‘CR338’ pistils were significantly lower than those of the ‘NC28’ pistils, implying that the clone ‘CR338’ was more freezing-tolerant than the clone ‘NC28’ (Figure 1b).

3.2. Physiological Responses of P. sibirica Pistils Under Cold Stress

To investigate the freezing tolerance of P. sibirica pistils under cold stress, six physiological indices of two clones were determined. As presented in Figure 2a,b, compared with the control (6 °C), the SOD and POD activities markedly increased at −4 °C but subsequently decreased at −6 °C in both P. sibirica clones. The SOD and POD content of the ‘CR338’ pistils was consistently higher than that of the ‘NC28’ pistils. Similarly, the soluble sugar and soluble protein contents of the ‘CR338’ pistils were higher than those of the clone ‘NC28’ pistils after freezing stress (Figure 2c,d). Furthermore, in response to freezing stress (−4 °C~−2 °C), the MDA content and electrical conductivity were higher in the ‘NC28’ pistils than in the ‘CR338’ pistils (Figure 2e,f). When the temperature dropped to −6 °C, the MDA content and electrical conductivity were not significantly different in the two P. sibirica clones (Figure 2e,f).

3.3. RNA Isolation, Library Construction, and Sequencing

In order to clarify the molecular mechanisms underlying the differences in the freezing tolerance of the two analyzed P. sibirica clones, transcriptome sequencing was performed with the Illumina PE system. Each sample had total RNA concentrations ranging from 436 ng/uL to 948 ng/μL: OD260/280 ranged from 2.15 to 2.19, OD260/230 ranged from 1.91 to 2.26, and the RIN value ranged from 5.9 to 6.7 (Table S2). The results demonstrated that the total RNA gathered was of sufficient quality and met the experimental conditions. Following the removal of low-quality reads, 44,131,376~47,158,396 clean reads were identified. The Q30 percentages for the six samples ranged from 87.64% to 93.48%, while the clean reading rate was between 97.59% and 98.86%. The percentage of the mapped rate varied from 89.68% to 94.13%, while the GC content ranged from 45.79% to 46.26% (Table S3). The findings indicated that the volume of data satisfied the quality requirements needed for further analysis.

3.4. Unigene Annotation and Functional Classification

With a GO analysis, 14,717 unigenes were categorized into 38 functional groups that belonged to the 3 primary GO categories (cellular component, molecular function, and biological process) (Figure 3, Table S4). The top three subcategories within the biological process (BP) category were metabolic process (6025 unigenes; GO:0008152), cellular process (5021 unigenes; GO:0009987), and single-organism process (4577 unigenes; GO:0044699). In addition, there were four subcategories within the BP category with more than 500 unigenes, which were in biological regulation, localization, regulation of biological process, and response to stimulus. The main subcategories within the cellular component (CC) category were cell part (1359 unigenes; GO:0044464) and cell (1359 unigenes; GO:0005623). The two largest subcategories in the molecular function (MF) category were binding (8467 unigenes; GO:0005488) and catalytic activity (6297 unigenes; GO:0003824).

3.5. Analysis of the Differentially Expressed Genes (DEGs)

DEGs are defined as those that are expressed and have an adjusted p-value < 0.05 and log2 (fold change) > 1. A total of 6001 genes related to freezing stress in pistil were differentially expressed with 3931 up-regulated genes and 2070 down-regulated genes. Many of these DEGs belong to the WRKY, NAC, bHLH, MADS, MYB, Homeobox, F-box, ZEP, SBT, ERF, CAMTA, and MDH gene families related to freezing resistance (Figure 4, Table S5).

The DEGs were categorized into four clusters, and genes in the same cluster had similar trends of expression level changes under different treatment conditions (Figure 5, Table S6). The trends of DEGs in different resistant clones can be categorized as slowly down-regulated (sub_cluster_1), sharply up-regulated (sub_cluster_2), sharply down-regulated (sub_cluster_3), and slowly up-regulated (sub_cluster_4).

3.6. GO Enrichment Analysis of DEGs

A GO functional classification analysis of the DEGs between the cold-tolerant clone ‘CR338’ and the cold-sensitive clone ‘NC28’ was conducted. The results are shown in Figure 6. Within the BP category, DEGs were mainly enriched in transmembrane transport (143 genes; GO:0055085), cell communication (116 genes; GO:0007154), and cellular response to stimulus (104 genes; GO:0051716). In addition, there are many classifications related to stimulus and stress response, e.g., cellular response to DNA damage stimulus (9 genes; GO:0006974), response to endogenous stimulus (21 genes; GO:0009719), response to stress (70 genes; GO:0006950), response to oxidative stress (24 genes; GO:0006979), response to hormone (21 genes; GO:0009725), and photosynthesis (16 genes; GO:0015979). In the MF category, the DEGs were mainly enriched in transmembrane transporter activity (135 genes; GO:0022857), transferase activity, transferring glycosyl groups (127 genes; GO:0016757), and hydrolase activity, acting on glycosyl bonds (117 genes; GO:0016798). Furthermore, this category includes a large number of functional annotations that are associated with low-temperature stress, e.g., calcium ion binding (49 genes; GO:0005509), antioxidant activity (28 genes; GO:0016209), calcium-dependent phospholipid binding (5 genes; GO:0005544), and protein serine/threonine kinase activity (4 genes; GO:0004674). In the CC category, the DEGs were mainly distributed in nucleus (86 genes; GO:0005634), cytoplasmic part (57 genes; GO:0044444), and cell periphery (45 genes; GO:0071944) (Table S7). The findings suggested that these functional classifications served a core role in the response to cold stress in P. sibirica.

3.7. KEGG Pathway Enrichment Analysis of DEGs

The KEGG pathway cluster analysis was performed on DEGs between the cold-tolerant clone ‘CR338’ and the cold-sensitive clone ‘NC28’ of P. sibirica to investigate the primary metabolic pathways. The DEGs were mapped to 114 biological pathways, and 20 of them were significantly enriched (Figure 7). They mainly included plant–pathogen interactions (77 genes; pmum04626), plant hormone signal transduction (75 genes; pmum04075), carbon metabolism (62 genes; pmum0120), starch and sucrose metabolism (53 genes; pmum00500), biosynthesis of amino acids (45 genes; pmum01230), glycolysis/gluconeogenesis (40 genes; pmum00010), MAPK signaling pathway-plant (39 genes; pmum04016), oxidative phosphorylation (39 genes; pmum00190), photosynthesis (19 genes; pmum00195), flavonoid biosynthesis (17 genes; pmum00941), circadian rhythm-plant (16 genes; pmum04712), photosynthesis–antenna proteins (10 genes; pmum00196), and so on (Table S8). These pathways might play a crucial role under cold stress in P. sibirica.

3.8. Validation of DEGs by RT-qPCR

In order to verify the accuracy of the DEG expression levels, 20 genes (10 up-regulated and 10 down-regulated genes) were selected for an RT-qPCR analysis (Figure 8). In the cold-tolerant clone ‘CR338’ and the cold-sensitive clone ‘NC28’, PsCDPK5, PsbHLH18, PsHB1, PsbZIP1, PsVRN1, PsPOD1, PsEARLI1, PsPPO, PsPODA2, and PsWD40 expression levels were markedly up-regulated, whereas PsMYB4, PsERF13, PsMYB30, PsSIZ1, PsMYB44, PsSABP2, PsHHP1, PsEIN3, PsERA8, and PsEDS5 expression levels were markedly down-regulated. RT-qPCR results for all 20 genes agreed with the RNA-seq results, confirming the accuracy of the RNA-seq investigation (Table S5).

4. Discussion

In recent decades, extreme weather events have been occurring more frequently and with greater intensity, and agroforestry disasters have become more severe as a result of global climate change [28]. Spring frost is an event in which temperatures plummet below 0 °C due to a sudden drop in temperature, severely harming plant growth and development [29,30]. Development of apricot reproductive organs is often subject to spring frosts [10]. This will have a huge impact on apricot flowering and fruiting. ‘CR338’ is one of the main P. sibirica clones and is widely recognized as a cold-tolerant clone, while the clone ‘NC28’ is extremely vulnerable under cold stress (Figure 1). At present, the cold tolerance of clone ‘CR338’ is superior to that of clone ‘NC28’, mainly in phenotypic characteristics and physiological and biochemical indexes measured (Figure 1 and Figure 2). In this study, we examined the cold resistance (artificial frost simulation experiment) of two P. sibirica clones from the morphological, physiological, and molecular aspects. The findings laid a theoretical foundation for further research on the mechanisms behind P. sibirica responses to natural frost in the field.

The main damage caused by spring frost to the reproductive organs of plants is the browning of the tissues, which leads to the inhibition of their growth [31]. Under freezing conditions, water loss occurs in plant organs, where water will rapidly escape the cell and then form ice nuclei, which diffuse into the extracellular plasma, leading to cell dehydration [32]. If the ice crystals diffuse into the cell, irreversible damage can occur. In this study, the floral organs and pistils of Siberian apricot clones ‘CR338’ and ‘NC28’ subjected to low-temperature stress turned brown, with the clone ‘NC28’ showing more pronounced browning and a higher rate of freezing of the pistils (Figure 1a). These results are consistent with those of previous studies [10].

Plants subjected to freezing stress produce ROS [33]. Excess ROS can harm plant tissues or even kill them. Plants produce ROS-scavenging antioxidant enzymes (SOD and POD) to protect cells from damage [34]. Freezing stress causes changes in cellular osmotic potential. At this time, plants increase the content of SS and SP to maintain the osmotic pressure [35,36]. On the other hand, the formation of ROS disrupts the cell membrane system and increases the MDA level and electrical conductivity due to lipid peroxidation [2]. It was discovered in this study that SOD and POD activities of clones ‘CR338’ and ‘NC28’ were significantly increased under freezing stress (Figure 2a,b). Clone ‘CR338’ had stronger antioxidant enzyme activity than ‘NC28’, which could promptly and effectively scavenge excessive ROS and effectively inhibit lipid peroxidation in cell membranes (Figure 2a,b). Meanwhile, ‘CR338’ synthesized more soluble sugars and soluble proteins than ‘NC28’, suggesting that ‘CR338’ may have promoted the expression of related genes and induced more osmoregulatory factors to resist freezing stress (Figure 2c,d). These findings are in line with those of previous studies [2]. With the temperature decreasing, both MDA content and electrical conductivity kept increasing in clones ‘CR338’ and ‘NC28’, and the content of ‘NC28’ was higher than that of ‘CR338’, indicating that it was more severely damaged (Figure 2e,f). When the temperature reached −6 °C, the MDA content and electrical conductivity between the two clones did not significantly change (Figure 2e,f), suggesting that this temperature had a serious impact for floral organ development of P. sibirica. The −4 °C can be used as a critical temperature to measure the strength of frost resistance of the floral organs of P. sibirica, and similar results were found in P. armeniaca [10].

Transcription factors play an important role in plants’ responses to low-temperature stress [37]. In this study, many of the highly expressed DEGs belonged to the WRKY and bHLH gene families (Figure 4 and Table S5). This is consistent with previous findings that many members of the WRKY and bHLH families are involved in the floral organ development and pistil frost resistance of P. sibirica [2,5]. The expression of CsbHLH18 in sweet orange (Citrus sinensis) was consistently up-regulated after cold stress, and it could reduce the ROS content by regulating the expression of CsPOD, which could maintain the intracellular balance of ROS and improve the cold tolerance [38]. In this study, we also found that PsbHLH18 was significantly expressed in clone ‘CR338’ (Table S5 and Figure 8), suggesting that it may play a role in the freezing resistance of P. sibirica. VRN1 belongs to the MADS-box transcription factor, and the expression level of HvVRN1 in Hordeum vulgare increased with the prolongation of the cold treatment time during low-temperature-induced seed germination [39]. In this study, the expression of PsVRN1 in clone ‘CR338’ was significantly higher than that of ‘NC28’ (Table S5 and Figure 8), suggesting that this gene may positively regulate plant freezing resistance. In breeding practice, PsVRN1 can be genetically transformed into P. sibirica to improve the freezing resistance of seedlings, which will facilitate its introduction. The MYB gene family is also closely related to plant growth and development and adversity stress [40]; in this study, the expression of some members was significantly decreased in ‘CR338’, such as PsMYB4 and PsMYB44 (Table S5 and Figure 8), suggesting that they may negatively regulate frost resistance in P. sibirica. These results are consistent with those of previous studies [41]. During fruit development in banana (Musa acuminate), low-temperature stress inhibits MaMYB4 expression, which affects its interaction with MaHDA2, promotes the expression of MaFADs, and enhances fatty acid biosynthesis [42]. Therefore, genetic transformation and functional validation of PsMYB4 appear to be necessary because it may be able to improve both flower and fruit cold hardiness. In addition to transcription factors, an increase in related protease synthesis genes also enhances freezing tolerance in plants [43]. PsPOD1 was a key regulatory gene for peroxidase synthesis and was significantly expressed in freezing-tolerant clones, which is also consistent with the results shown by physiological indicators (Figure 2 and Figure 8). In sugarcane (Saccharum spp. hybrid), overexpression of the cytoskeletal microtubule synthesis gene SoTUA induced the expression of SoPOD1, which was able to significantly increase the cold resistance of seedlings, and this characteristic was stably inherited to the next generation [44]. Considering this feature, PsPOD1 can be used for early identification of pistil freezing resistance in different clones of P. sibirica, which will further save breeding time. AtSIZ1 and AtHHP1 overexpressed transgenic plants in Arabidopsis thaliana enhance plant freezing tolerance [45,46]. MdSIZ1 in apples (Malus domestica) interacts with MdMYB1 via sumoylation, which in turn induces anthocyanin synthesis to improve plant cold resistance [47]. The expression of CpSIZ1 in wintersweet (Chimonanthus praecox) is distinctly tissue-specific, and it is significantly expressed mainly at the late stage of leaf and flower development to improve the low-temperature tolerance [48]. Genetic transformation of the PuICE1 from pear (Pyrus ussuriensis) into tomato enhances cold tolerance by binding to the PuHHP1 protein, which increases the chlorophyll content, protective enzyme activity, and proline content, and decreases electrolyte leakage [49]. In contrast, this study revealed that PsSIZ1 and PsHHP1 negatively regulated freezing resistance in P. sibirica (Figure 8). This result suggests that gene functions may diverge in woody and herbaceous plants or different stages of growth and development. In conclusion, the discovery of these genes provides new perspectives and ideas for freezing-resistant breeding and germplasm selection in P. sibirica.

Ca²⁺ plays a key role as a second messenger in response to cold stress in plants [50]. A sharp drop in temperature alters the homeostasis of the cell membrane and causes it to rigidify, which leads to a rapid increase in the concentration of Ca²⁺-receiving signals in the cytoplasm through several calcium sensors (CaM, CML, and CDPKs). The CAMTA transcription factor is induced to express itself under low-temperature stress, where it exists as a CaM-binding protein and interacts with other related genes to improve cold tolerance [51]. In this study, most genes related to Ca²⁺ signaling (calcium ion binding and calcium-dependent phospholipid binding) were regulated by freezing stress (Table S5 and S7), implying that Ca²⁺-mediated signaling pathways may be important for P. sibirica responses to spring frost. In wild grapes (Vitis amurensis), overexpression of VaCDPK20 significantly increased cold tolerance in the plant [52]. In this study, PsCDPK5 was more highly expressed in ‘CR338’ than in ‘NC28’ (Figure 8). Accordingly, PsCDPK5 may have a key role in the early stages of spring frost.

Plants rapidly synthesize fatty acids and amino acids under freezing stress [53,54]. The formation of these substances contributes to the stability of the cell membrane and the osmotic balance within cells. Lipids are not only a major component of the cell membrane but also play a role as signaling substances in response to cold stress [55]. At this time, lipids adsorb proteins associated with antifreeze to the cell membrane, enabling them to exercise their functions by changing the conformation and activity of proteins and metabolites. Numerous DEGs in the current study were mostly primarily involved in lipid metabolic process, lipid biosynthetic process, and biosynthesis of amino acids (Figure 7 and Table S8). These results suggest that the antifreeze clones may enhance their protection by regulating the synthesis and metabolism of fatty acids and amino acids. Carbohydrates are important osmoprotectors against cellular dehydration under adverse conditions, including low-temperature stress [56]. In this study, many DEGs were mainly distributed in cellular carbohydrate biosynthetic and metabolic processes and in starch and sucrose metabolism (Figure 7 and Table S8). Physiological indexes also found that the SS content of the clone ‘CR338’ was significantly higher than that of ‘NC28’ (Figure 2c). These findings provide strong evidence that carbohydrates may have an important protective role for P. sibirica. Cold stress decreases the rate of photosynthesis in forest trees by inhibiting enzyme activities and slowing down stomatal opening and CO₂ uptake rates, which will lead to slower plant growth and lower yields [57,58]. In this study, we found that the number of up-regulated expression related to photosynthesis was significantly higher than the number of down-regulated expression (Figure 6 and Figure 7). This suggests that the freezing-resistant clones may have induced an increased number of genes related to photosynthesis to ensure normal plant development.

5. Conclusions

In this study, two P. sibirica clones that were subjected to freezing stress had their phenotypic and physiological traits examined. Compared with the clone ‘NC28’, the clone ‘CR338’ had a lower pistil freezing rate, MDA content, and electrical conductivity, and higher SOD, POD, soluble sugar, and protein content, implying that ‘CR338’ was more tolerant to freezing stress than ‘NC28’. Transcriptome sequencing was used for the regulation and expression profiling of genes in P. sibirica pistil under freezing stress. Many DEGs are closely related to response to stress, oxidative stress, and hormones. They are mainly involved in Ca²⁺ and MAPK signaling, carbohydrate biosynthesis and metabolism, biosynthesis of amino acids, and photosynthesis. PsbHLH18, PsMYB4, PsMYB44, PsPOD1, and PsCDPK5 may play important roles in the freezing resistance of P. sibirica flower organs. These important candidate genes can subsequently be genetically transformed to validate gene function and observe their adaptability to environmental changes, or they can be used for introgression or early screening and characterization of frost-resistant clones. Taken together, these findings provide theoretical support for further studies on the complex mechanisms of spring frost response in P. sibirica and other economic forest species. They are also important for the early identification of excellent floral organ freezing-resistant germplasm resources and the breeding of new freezing-resistant P. sibirica varieties.