Transcriptome Analysis of Apricot Kernel Pistils Reveals the Mechanisms Underlying ROS-Mediated Freezing Resistance

Abstract

Spring frost is a major limiting factor in the production and cultivation of apricot kernels, an ecological and economic dry-fruit tree in China. The frequent occurrence of spring frost often coincides with the blooming period of apricot kernels, resulting in significant damage to floral organs and reductions in yield. We investigated the molecular signature of pistils from two apricot kernel cultivars with different frost-resistance levels using transcriptome data. A total of 3223 differently expressed genes (DEGs) were found between two apricot kernel cultivars under freezing stress, including the bHLH and AP2/ERF-ERF transcription factors. Based on KEGG analysis, DEGs were mostly enriched in the biosynthesis of the secondary metabolites, in the metabolic pathways, and in plant-hormone signal transduction. The co-expression network, which included 81 hub genes, revealed that transcription factors, protein kinases, ubiquitin ligases, hormone components, and Ca²⁺-related proteins coregulated the ROS-mediated freezing response. Moreover, gene interaction relationships, such as ERF109-HMGCR1, ERF109-GRXC9, and bHLH13-JAZ8, were predicted. These findings revealed the regulatory factors for differences in frost resistance between the two tested apricot kernel cultivars and contributed to a deeper understanding of the comprehensive regulatory program during freezing stress. Some of the hub genes identified in this work provide new choices and directions for breeding apricot kernels with a high frost resistance.

1. Introduction

Frost is a common meteorological disaster which refers to a sudden drop in the air and surface temperature to below 0 °C. It has been recognized as a major threat to plant growth, development, and agricultural and forestry productivity [1,2]. Since climate warming has increased temperatures in early spring, perennial plants have become increasingly vulnerable to lower temperatures due to phenological shifts, such as advanced flowering time [3]. Spring frost has been shown to cause irreparable losses to vegetables, fruit trees, and crops [4].

An apricot kernel, an apricot (Prunus armeniaca L.) plant with almonds as its main use, is also a fresh fruit with a unique taste and mainly includes a big flat apricot (Armeniaca vulgaris × sibirica) and Siberian apricot (Armeniaca sibirica L.). It is an important raw material for the food and pharmaceutical industries and is mainly distributed in northern China. Spring frost frequently occurs during their flowering time between late March and mid-April [5]. Among the apricot flower organs, the freezing resistance of the pistils is the weakest, followed by the stamen and petals [6]. Spring frost causes severe damage to apricot kernels’ reproductive organs, resulting in significant yield loss [7]. Frost injury is the main limiting factor in apricot kernel production.

Cold stress is an environmental stress that can be divided into chilling stress (0–15 °C) and freezing stress (<0 °C) (e.g., spring frost) [8]. Plants that suffer from freezing stress have developed sophisticated cold-acclimation mechanisms that improve their freezing tolerance upon exposure to nonlethal lower temperatures [9,10]. A series of cellular responses and molecular strategies are initiated when a plant perceives freezing stress, such as the production of ROS and osmolytes, and changes in the cytosolic Ca²⁺ concentration, hormone content, and gene expression [10,11,12]. CBF genes, which are rapidly induced by low temperatures, play central roles in cold acclimation. Many transcription factors (TFs) (e.g., ICEs, CAMTAs, and MYB15) act as upstream regulators that regulate the expression of CBFs [8,10]. Among them, ICE1, a bHLH TF, is the best-characterized transcriptional activator of CBF genes [13]. Post-translational modifications, such as ubiquitination and phosphorylation, are important for the function of ICE1 in cold tolerance. For example, the protein kinases OST1 and MPK3/6 can phosphorylate ICE1, affecting its transcriptional activity to regulate CBF expression and cold tolerance [14].

In plants, ROS act as a double-edged sword. Excessive ROS accumulation due to stress induces oxidative stress, which can damage plant cells; at this time, ROS-scavenging systems consist of an endogenous defensive mechanism that comprises different enzymatic (e.g., superoxide dismutase, SOD, and catalase, CAT) and nonenzymatic (e.g., glutathione, GSH) antioxidants that are activated to maintain ROS levels [15,16,17]. Moreover, ROS, as messenger molecules, participate in acclimation responses to freezing stress. ROS can interact with different hormones (e.g., ET, JA, and BRs) to control gene expression and induce physiological changes in response to cold stress [18]. In addition, emerging evidence indicates that some key signaling components participate in ROS-mediated stress response processes, such as messenger molecules (e.g., Ca²⁺ and NO), protein kinases (e.g., CIPKs and MAPKs), and TFs (e.g., MYC and MYB) [18,19].

In apricot kernels, research on freezing resistance primarily focuses on the physiological level and differences between different varieties [20,21]. It has been reported that the activities of antioxidant enzymes, such as SOD, in apricots show a change pattern in which they first increase and then decrease under freezing stress, with their activities higher in variety with a strong freezing resistance [7,22]. Our previous study found that many regulators (such as TFs and protein kinases) and some of the genes involved in the oxidation reduction process were regulated in apricot under natural spring frost conditions via transcriptome analysis [7]. However, the underlying relationship and functional mechanism of differentially expressed genes (DEGs) are still largely unknown. We investigated the comparative transcriptome of pistils in this work based on two apricot kernel varieties (‘Weixuan 1’ and ‘Longwangmao’) (Armeniaca vulgaris × sibirica) with different frost-resistance levels under simulated spring frost conditions. Our study aims to analyze the DEGs and different biological processes between ‘Weixuan 1’ and ‘Longwangmao’ and elucidate ROS-mediated molecular mechanisms in response to freezing stress. This may provide new insights into the response mechanisms underlying freezing stress in apricot kernels.

2. Materials and Methods

2.1. Plant Materials and Treatment

For the analysis of apricot kernel pistils’ freezing resistance, two main cultivated varieties were selected, namely ‘Weixuan 1’ and ‘Longwangmao’. ‘Weixuan 1’ is a frost-resistant variety selected from ‘Longwangmao’ through bud mutation [23].

During hibernation, flower branches from cold-tolerant ‘Weixuan 1’ (CtW) and cold-sensitive ‘Longwangmao’ (CsL) were collected from the Apricot Germplasm Resource nursery in Shanxi, China. Flower branches were brought to full bloom via hydroponic cultures in an incubator (20 °C) and were treated with temperatures of −2 °C, −3 °C, and −4 °C for 1 h. The cooling method involved reducing the temperature from 20 °C to 2 °C at a rate of 10 °C/0.5 h, and then reducing the temperature from 2 °C to the treatment temperature at a rate of 3 °C/h in a low-temperature incubator. The flower branches were then placed in the incubator (20 °C) to recover for 3 h; that is, when the browning of the pistils showed no significant change. The pistils undergoing and not undergoing freezing treatments were collected for transcriptome sequencing and quantitative real-time PCR (qRT-PCR) analysis in liquid nitrogen and were stored at −80 °C.

2.2. RNA-seq

The pistils of CsL and CtW undergoing and not undergoing freezing treatments were used for RNA exaction. According to the differences in the temperature treatments between CsL and CtW, the samples were named CsL1 (20 °C), CsL2 (−2 °C), CsL3 (−3 °C), CsL4 (−4 °C), CtW1 (20 °C), CtW2 (−2 °C), CtW3 (−3 °C), and CtW4 (−4 °C), and every sample contained three biological replicates (every five flower branches represented a biological repeat, containing 50–70 pistils). Total RNA was extracted with the RNAprep Pure Plant Plus Kit (polysaccharide- and polyphenolic-rich; Tiangen, Beijing, China). RNA quality estimation, including purity, integrity, and concentration, was checked by 1% agarose gel, a Nano Photometer spectrophotometer (Implen, Westlake Village, CA, USA), a Qubit 2.0 Fluorometry (Life Technologies, Carlsbad, CA, USA), and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). High-quality RNAs were used to construct cDNA libraries sequenced by Illumina paired-end sequencing technology on an Illumina HiSeq platform (Illumina, San Diego, CA, USA) at Metwell Biotechnology Co., Ltd. (Wuhan, China). The raw data were cleaned by removing adapter sequences, the reads with more than 10% N bases, and low-quality reads with a percentage of base quality value ≤ 5 exceeding 50%. The clean reads were mapped onto the apricot (Prunus armeniaca L.) reference genome using the HISAT2 software using the default parameters [24].

2.3. Identification of DEGs and Enrichment Analyses

The number of reads mapped to each gene were counted using HTSeq v0.6.1, and the fragments per kilobase of transcript per million fragments mapped (FPKM) of each gene were calculated based on the gene length and the read count mapped to the gene. The DEG analysis between two groups was performed using the DESeq2 R package. The false discovery rate (FDR) was obtained from p-values adjusted using the Benjamini–Hochberg method. The DEGs were screened with |log2Fold Change| >= 1 and FDR < 0.01. KEGG enrichment analyses of the DEGs were conducted using KOBAS software.

2.4. Weighted Gene Co-Expression Network and Hub Genes Analysis

The antioxidant enzymes (POD, SOD, and CAT) were determined using water-soluble tetrazolium salt-(WST-1), guaiacol-, and hydrogen-peroxide-based methods, respectively [25]. The mixed pistils undergoing and not undergoing freezing treatments were used to analyze the activities of the antioxidant enzymes, and each experiment contained four biological replicates. SPSS 25 was used to analyze the data of the enzyme activities by one-way ANOVA and Duncan’s multiple comparison analysis.

The activity of the antioxidant enzymes and DEGs was used for weighted gene co-expression network analysis (WGCNA). WGCNA was conducted in R using the default parameters. The FPKM values of the DEGs were normalized, and Pearson’s correlation coefficient was calculated for each pair of genes to construct an adjacency matrix. Gene modules were identified based on TOM and converted by the adjacency matrix, using the WGCNA package. The correlation between modules and traits was estimated. The hub genes within a selected module were screened by kME (intra-module connectivity) > 0.9 and GS (gene significance) > 0.2. Cytoscape (version 2.7.2) was used to visualize the relationships between hub genes.

2.5. Protein–Protein Interaction (PPI) Network Prediction

For the PPI analysis, hub genes associated with the activity of antioxidant enzymes were used to retrieve genes interacting with STRING 11 in Prunus armeniaca var. bungo. A required confidence score (combined score) greater than 0.4 was used as the threshold for the interaction. The disconnected genes were hidden in the network. The PPI network was constructed using STRING and was further analyzed using Cytoscape (version 2.7.2).

2.6. Verification of qRT-PCR Analysis

Nine genes were selected from the hub genes related to the antioxidant enzyme for qRT-PCR analysis. Total RNA extraction was performed with a TIANGEN kit (Beijing, China), and cDNA was synthesized using a reverse transcription kit (Takara Dalian, Japan). A qRT-PCR test containing three biological replicates was conducted using the KAPA SYBR FAST qPCR Master Mix (Kapa Biosystems, Boston, MA, USA) on a LightCycler 480 II Real-time PCR Instrument (Roche) according to the manufacturer’s protocols. In that test, 18S was used as a reference gene. The primers are listed in Table S1.

3. Results

3.1. Transcriptome Analysis of Apricot Kernel Pistils under Freezing Stress

RNA-seq data were generated for eight different freezing-treated and untreated pistils to explore the molecular mechanism of the freezing resistance between CsL and CtW pistils. In total, 120.78 Gb of clean reads were obtained from 24 libraries, ranging from 4.16 to 6.02 Gb per library, with an average GC content of 45.81% (Table 1). The high Q30 (>91.73%) represents bases with error rates < 0.1% for libraries showing high-quality RNA-seq. The rate of clean reads mapped to the apricot genome ranged from 87.55% to 94.74%, in which uniquely mapped reads exceeded 84.85%.

3.2. Identification and Functional Analysis of DEGs of Apricot Kernel Pistils under Freezing Stress

In total, 5206, 2032, and 3223 DEGs were identified under freezing treatment in CsL, CtW, and CsL vs. CtW. Among all of the DEGs, 509 shared DEGs that were found in the three comparison groups (CsL, CtW, and CsL vs. CtW). In the three CsL freezing-treatment groups, 438 genes were up-regulated, and 68 genes were down-regulated. In the three CtW freezing-treatment groups, 367 genes were up-regulated, and 142 genes were down-regulated; in CsL vs. CtW, 355 genes were up-regulated, and 157 genes were down-regulated (Figure 1a).

DEGs were characterized using KEGG databases to understand their biological roles in CsL and CtW. There were 15, 21, and 15 significantly enriched pathways (p < 0.05) in CsL, CtW, and CsL vs. CtW, respectively. In CsL, most of the DEGs were enriched in plant-hormone signal transduction (7.17%), the MAPK signaling pathway (4.89%), the biosynthesis of secondary metabolites (25.40%), and plant–pathogen interaction (9.83%) (Figure 1b). In CtW, most DEGs were enriched in the biosynthesis of secondary metabolites (26.96%), plant–pathogen interaction (11.65%), plant-hormone signal transduction (6.28%), and phenylpropanoid biosynthesis (5.76%) (Figure 1c). In CsL vs. CtW, more DEGs were enriched in the biosynthesis pathways of the secondary metabolites (27.90%), in the metabolic pathways (45.42%), in plant-hormone signal transduction (6.24%), and in phenylpropanoid biosynthesis (6.00%) (Figure 1d).

3.3. Differentially Expressed Transcription Factors of Apricot Kernel Pistils under Freezing Stress

The transcriptional regulation of cold stress has been widely studied in model plants. To identify the TFs involved in apricot kernels’ response to freezing stress, we analyzed differentially expressed TFs in CsL and CtW under freezing stress. As shown in Figure S1, 423 DEGs were assigned to 50 TF families in CsL; 151 DEGs were assigned to 33 TF families in CtW; and 242 DEGs were assigned to 48 TF families in CsL vs. CtW. Genes belonging to the NAC, AP2/ERF-ERF, MYB, WRKY, and bHLH TF families in the CsL, CtW, and CsL vs. CtW groups accounted for more than 40%. Furthermore, there were 44 shared TFs, including 10 AP2/ERF-ERF TFs, 6 bHLH TFs, 5 MYB TFs, 3 NAC TFs, and 2 WRKY TFs in the three groups (Figure 2). Out of these TFs, the expression level of most genes under freezing stress was up-regulated compared to the control, except for the down-regulated PARG01786 (RAP2.4), PARG12349 (MYB6), PARG06699 (MYB21), and PARG30216 (NAC25) TFs in CtW and the PARG29164 (WRKY70) TF in CsL (Figure 2b). These findings demonstrate that these TFs, especially NAC, AP2/ERF-ERF, MYB, WRKY, and bHLH, may govern the transcriptional changes through both transcriptional activation and repression in response to freezing stress.

3.4. The Co-Expression Network Analysis of DEGs Related to the Antioxidant Enzyme Activity

Freezing stress produces excessive ROS scavenged by antioxidant mechanisms, such as enzymatic and nonenzymatic systems, to regulate cold resistance in plants [26,27]. We employed WGCNA to detect the co-expressed genes associated with antioxidant enzyme activities, including POD, SOD, and CAT, that may be involved in regulating freezing resistance in CsL and CtW (Table S2). The network was constructed with 3223 DEGs, and nine co-expressed gene modules labeled with different colors were determined (Figure 3a). Then, the correlations between the module eigengenes (MEs) and antioxidant enzyme parameters were analyzed. Only the magenta module with 173 genes was significantly associated with SOD activity (r = 0.73, p = 0.04) (Figure 3b). In addition, the module membership and GS of the magenta module were highly correlated (Table 2), further demonstrating that genes in the magenta module were significantly associated with SOD activity.

Subsequently, the DEGs in the magenta module that exhibited the strongest correlation with the SOD parameter were further analyzed. These genes were significantly up-regulated in CsL3 (treated with −3 °C) and in CtW2 (treated with −2 °C) (Figure S2), implying that the genes in CtW were modulated earlier than those in CsL in order to adapt to and defend against freezing stress when suffering from lower temperatures. Hub genes within the magenta module, referring to the most highly connected nodes within the module, were used to construct the gene network. A total of 81 hub genes containing 77 up-regulated and 4 down-regulated genes in CsL vs. CtW were identified based on kME > 0.9 and GS > 0.2 (Figure 4; Table S3). The co-expression network had ten hormone-related genes and five Ca²⁺-related genes: ERF025, ERF109, JAZ8, JAZ10B, CAMBP25, and PBP1-like genes were highly connected within the module. Moreover, twelve TFs, including one CBF/DREB subfamily member, three WRKY members, three bHLH members, and two NAC members, were related to the SOD parameter, and NAC090, bHLH35, ZAT11, and GTE12 had higher connections within the module. In addition, twelve post-translational modified proteins with five protein kinases, one phospholipase PP2C25, five E3 ubiquitin-protein ligases, and the F-box protein SKIP27, were also connected with the SOD parameter; the kME of the PUB21-like, SKIP27, EFR, and PARG02353 (LRR receptor-like serine/threonine-protein kinase At3g47570) were high. Some structural functional genes, such as the redox genes HMCGR1 and GRXC9, and the jasmonic acid (JA) synthesis genes, namely AOS1 and OPR3, were correlated with the SOD parameter (Figure 4; Table S3).

In addition, the co-expression relationship of TF–TF, TF–post-translational modified protein, and TF–structural functional genes were identified in the gene network (Figure 4). The TFs ERF025, WRK18, and ERF109, as well as the E3 ubiquitin-protein ligases ATL7 and PUB23, showed a higher correlation with CBF3. The repressor proteins JAZ8 and JAZ10B in JA signaling were co-expressed with NAC090, WRKY53, ERF109, PUB21-like, and SKIP27. A co-expression relationship between HMCGR1 and the TFs (ERF109, WRKY40, NAC090, bHLH35, and ZAT11) was found in the gene network. These relationships can act as a reference for the research on the freezing response.

3.5. Interaction Network Analysis of Hub Genes Related to Antioxidant Enzyme Activity

A PPI regulatory network of 81 hub genes in the magenta module was constructed to explain the potential regulatory mechanism in apricot kernels responding to freezing stress. As shown in Figure 5, JAZ8 has a direct relationship with bHLH13, WRKY18/40, ERF109, and GATA25, and especially with bHLH13, in addition to the genes in the JA signal pathways (JAZ1/10/12, NINJA, and COI1) and the JA biosynthesis genes (OPR3, AOC3/4, and AOS). bHLH13 plausibly interacts with the JA-related genes (JAZ1/8/10, NINJA, COL1, OPR3, and AOS1), GATA25, and bHLH92, whereas bHLH35 may only interact with JAZ1 and bHLH92. There may also be interactions between ERF109 and WRKY18/40, HMGCR1, GRXC9, and PP2C25. Moreover, WRKY18/40 and bHLH92, GRXC9 and JAZ12, and RGI3 and RGL1 were found to have direct relationships. PP2C25, NPR4, CCR4, ZAT11, and PBP1 were all related to PP2Cc. The JA-related genes (JAZ8, AOS1, and OPR3) and bHLH13 were considered to be key genes in the network, suggesting that they may play an important role in the ROS-mediated freezing response.

3.6. qRT-PCR Validation of Key DEGs Involved in the Response of Apricot Kernels to Freezing Stress

We selected nine key genes in the magenta module for qRT-PCR analysis to verify the expression pattern’s accuracy for these hub genes in the transcriptome data. With the exception of NAC090 in CtW3 and CtW4, the expression levels of ERF109, ZAT11, PBP1, NAC090, PP2C25, and PUB21 were higher under lower temperatures, higher in CsL3 or CsL4 in ‘Longwangmao’, and higher in CtW2 in ‘Weixuan 1’ (Figure 6). JAZ8 was induced by freezing stress in CsL, whereas it had no significant change or decrease in CtW. In addition, the expression levels of bHLH35 and OPR3 obviously increased in CsL4 in ‘Longwangmao’ and decreased in CtW3 in ‘Weixuan 1’. The qRT-PCR results were consistent with the RNA-Seq results, indicating the reliability of the RNA-Seq data.

4. Discussion

4.1. Transcriptional Regulation in Freezing Response of Apricot Kernels

TFs are important regulators for controlling gene expression to modulate the stress response. Many TFs regulate the expression of cold-stress-responsive genes (CORs) and adapt the tolerance of plants to cold stress, such as the AP2/ERF, NAC, WRKY, MYB, and bZIP TFs [28]. In the pistils of CsL and CtW, these families were also the main differentially expressed TFs under freezing stress (Figure S1). Recent studies show that CBFs elevate antioxidant enzymes to regulate cold tolerance [29,30]. Consistent with these observations, one CBF/DREB TF, CBF3, was regulated by freezing stress and was found to be associated with antioxidant enzyme activity (SOD) in apricot kernels (Figure 4). No ICEs, which are key inducers of CBFs expression, were differentially expressed in CtW and CsL under freezing stress, whereas other bHLH TFs (such as bHLH13 and bHLH35) were induced by freezing stress and were correlated with antioxidant enzyme activity (Figure 2 and Figure 4). Moreover, bHLH13 and bHLH35 were predicted to interact with JAZ8 or JAZ1, suggesting that they link with JA signaling to regulate the SOD-mediated freezing response.

In addition to the ICE1-CBF pathway, many other TFs regulate plant resistance in a CBF-dependent or CBF-independent way [10]. For example, MYB15 inhibits the expression of CBFs and negatively regulates the freezing tolerance of Arabidopsis thaliana [31]. In soybean, the overexpression of GmNAC20 increases the activity of antioxidant enzymes and enhances cold tolerance via the CBF-COR pathway [32,33]. On the contrary, some CBF-independent TFs, such as MYB73, WRKY33, and ZAT10, function parallel to CBFs and modulate COR expression [10,28,30]. In this work, some TFs, including ERF109, bHLH35, and WRKY18, were found to be related to antioxidant enzyme activity and to be co-expressed with CBF3 (Figure 4). Furthermore, ERF109 may directly regulate the expression of redox proteins HMGCR1 and GRXC9 in response to freezing stress (Figure 5). Our results suggest that these TFs regulate the antioxidant-enzyme-mediated freezing response by modulating CBFs.

4.2. Post-Translational Regulation in the Freezing Response of Apricot Kernels

Besides transcriptional regulation, post-translational modifications, such as phosphorylation and ubiquitination, also play important roles in the responses of plants to cold [10]. Many protein kinases and phosphatases have been confirmed to regulate the cold tolerance of plants, such as OST1/SnRK2.6, MPK3, MPK6, and BIN2 [34]. In Arabidopsis, a type 2C phosphatase, ABI1, dephosphorylates protein kinase OST1 to repress kinase activity and to negatively regulate freezing tolerance [35]. In the pistils of apricot kernels, five protein kinases (RGI3, EFR, CCR4, GSO2, and At3g47570-like) and one phosphatase (PP2C25) were identified as hub genes in the magenta module: they were highly correlated with antioxidant enzyme activity under freezing stress (Figure 3 and Figure 4). Previous studies have shown that OsCPK24 functions as a positive regulator of cold tolerance by phosphorylating and inhibiting OsGrx10 to improve glutathione (an antioxidant) levels [36]. The PPI network predicted the direct relationship between PP2C25-ERF109, PP2Cc-ZAT11, and RGI3-RGL1 (Figure 5), indicating that these protein kinases and phosphatase may control the tolerance of apricot kernels to freezing stress by regulating the phosphorylation level of the targets related to the antioxidant process.

Ubiquitination and protein-degradation mediated by the E3 ubiquitin ligases are also important for cold signaling. E3 ubiquitin ligases such as HOS1, PUB25, ATL78, and ATL80 have been extensively studied to be involved in cold stress in plants [34,37,38]. Two U-box type E3 ubiquitin ligases, PUB25 and PUB26, improve their E3 ligases activity by the phosphorylation of OST1 and function as negative regulators in response to cold stress by targeting MYB15 for degradation [39]. Four PUBs E3 ubiquitin ligases (PUB21, PUB21-like, PUB23, and PUB35) were differentially expressed in apricot kernel pistils during freezing stress and were highly related to antioxidant enzyme activity (Figure 4). OsSRFP1, a RING finger E3 ligase, negatively regulates the activity of antioxidant enzymes and cold stress tolerance in Oryza sativa [40]. However, there are few studies on the mechanism in ubiquitin ligase that regulates antioxidant enzyme activity. In addition to PUB35, other E3 ubiquitin ligases have co-expressive relationships with ERF109 (Figure 4), suggesting that these E3 ubiquitin ligases regulate ROS homeostasis and the freezing resistance of apricot kernels through ERF109.

4.3. Ca²⁺ and Hormone Signaling in Freezing Response of Apricot Kernels

Ca²⁺ is an important secondary messenger in plant response to cold stress. Previous research has shown that calmodulin (CAM) activity is essential for the expression of CORs and CAMTAs, which harbor conserved CAM-binding sites and activate CBF expression [10,41]. In this study, two CAM-binding proteins, CAMBP60B and CAMBP25, were induced by freezing stress and were found to be associated with antioxidant enzyme activity in apricot kernels (Figure 4). However, our results did not observe a co-expression relationship between CAMBPs and CBFs. There is the possibility that CAMBP60B and CAMBP25 are involved in the freezing response in apricot kernel pistils via a CBF-independent pathway.

Plant hormones that play key roles in cold stress tolerance by regulating ROS balance have been found [18]. Cold-activated BZR1, a positive transcriptional factor in BRs signaling, directly promotes RBOH1 expression and H₂O₂ production [42]. In Arabidopsis and peas, SLs positively regulate chilling tolerance via increasing glutathione and ascorbate accumulation [43]. Two JAZ proteins (JAZ8 and JAZ10B), two ethylene-responsive TFs (ERF109 and CRF4), one DELLA protein RGL1, and the SA receptor NPR4 were found to be highly associated with antioxidant enzyme activity in apricot kernel pistils under freezing stress (Figure 4). Consistent with these findings, JA, SA, ET, and GA were also involved in the response to cold stress; PtrERF109 directly promotes the expression of PtrPrx1 to improve peroxidase activity [12,18,44]. Moreover, the direct relationships between JAZ8 and bHLH13, ERF109, WRKY18, and WRKY40 were predicted, suggesting that JA signaling cooperates with other TFs and hormones to affect the ROS-mediated freezing response. These potential functional genes involved in freezing stress can provide a choice for genetic-engineering-assisted breeding through gene-editing technology and a direction reference for molecular-marker-assisted breeding.

5. Conclusions

In the present study, we investigated gene co-expression networks in the SOD-mediated response to freezing stress in the pistils of two apricot kernel cultivars with a different level of frost resistance. One gene network was identified to correlate with the antioxidant enzyme SOD activity under freezing stress. The direct relationship of regulatory and functional hub genes within this network were predicted. Our study confirmed some novel hub genes and potential mechanisms underlying the variation in the freezing resistance of apricot kernels.

Author Contributions

X.L. and L.W. designed the research. X.L. performed the experiments and analyzed the data. X.L. and Y.Y. wrote the manuscript. H.X., D.Y. and Q.B. provided helpful comments on the work and manuscript. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data supporting the results are concluded in the article and supplementary information files. The RNA-seq data were deposited in NCBI SRA under the accession number PRJNA832066.

Acknowledgments

We thank the forestry and seedling workstation in Yuyang District, Yulin City, for the ‘Weixuan 1’ and ‘Longwangmao’ material.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Venn diagram and significantly enriched KEGG pathways of differently expressed genes (DEGs) under freezing treatment. (a) Venn diagram showing the shared and specific number of up-regulated and down-regulated DEGs identified in CtW, CsL, and CtW vs. CsL. CtW: cold-tolerant ‘Weixuan 1’; CsL: cold-sensitive ‘Longwangmao’. (b,c) Significantly enriched KEGG pathways of DEGs in CtW (b), CsL (c), and CtW vs. CsL (d).

Figure 2. Analysis of differently expressed transcription factors (TFs) under freezing treatment. (a) Venn diagram showing the shared and specific number (and ratio) of differently expressed TFs identified in CtW, CsL, and CtW vs. CsL. (b) Heatmap of overlapping TFs between CtW, CsL, and CtW vs. CsL. Each column represents the gene expression at different temperatures (20 °C, −2 °C, −3 °C, and −4 °C) in CtW and CsL.

Figure 3. WGCNA of DEGs and antioxidant enzyme activity. (a) Cluster dendrogram presents nine co-expression modules labeled with different colors. (b) Correlation heatmap between DEGs and antioxidant enzyme activity. Rows correspond to modules. The color (from blue to red) indicates the correlation value (from −1 to 1).

Figure 4. The co-expression network of hub genes in the magenta module. Red nodes refer to Ca2+-related genes; green nodes refer to genes involved in hormone signaling; purple nodes refer to genes involved in post-translational modification; blue nodes refer to TFs. Black words indicate up-regulated genes; green words indicate down-regulated genes. The size of the nodes is based on intra-modular connectivity (kME). The width of the edges is based on weight.

Figure 5. PPI regulatory network of hub genes in the magenta module. Red nodes refer to hub genes identified in the magenta module; blue nodes refer to genes in the genome that interacted with hub genes. The size of the nodes is based on degree. The width of the edges is based on the confidence score.

Figure 6. The relative expression levels of nine hub genes by qRT-PCR. 18S was used as the internal control. Error bars indicate SD.

Table 1

Summary of mapping transcriptome reads to reference sequence.

Variety	Treatment	Sample	Clean Reads	GC (%)	Q20 (%)	Q30 (%)	Read Mapped (%)	Unique Mapped (%)
‘Longwangmao’ (CsL)	CsL1	CsL11	60,245,518	45.92	97.21	92.15	94.14	91.47
		CsL12	48,251,936	45.93	97.01	91.73	93.41	90.81
		CsL13	48,950,240	46.25	97.06	91.86	93.94	91.31
	CsL2	CsL21	47,482,406	45.86	97.66	93.19	94.74	92.12
		CsL22	45,957,672	45.72	97.67	93.26	93.90	90.94
		CsL23	54,500,124	45.50	97.50	92.88	93.69	90.63
	CsL3	CsL31	56,350,610	45.30	97.64	93.21	93.29	90.70
		CsL32	56,506,172	45.78	97.65	93.25	92.27	89.72
		CsL33	55,646,854	45.85	97.67	93.31	92.45	89.86
	CsL4	CsL41	54,666,004	45.78	97.66	93.34	87.55	84.85
		CsL42	44,453,966	45.72	97.60	93.14	93.55	91.11
		CsL43	49,831,556	45.90	97.58	93.12	94.34	91.61
‘Weixuan 1’(CtW)	CtW1	CtW11	41,850,108	45.93	97.30	92.29	92.11	89.62
		CtW12	41,637,772	46.20	97.32	92.34	94.16	91.59
		CtW13	50,142,906	46.10	97.29	92.29	93.79	91.23
	CtW2	CtW21	45,960,022	45.49	97.62	93.05	94.17	91.66
		CtW22	55,251,738	45.38	97.75	93.40	94.46	91.98
		CtW23	51,158,772	45.88	97.61	93.14	94.37	91.71
	CtW3	CtW31	43,518,372	45.87	97.72	93.45	93.96	91.45
		CtW32	60,004,540	45.87	97.44	92.79	93.50	91.03
		CtW33	60,219,916	45.69	97.55	93.03	94.03	91.44
	CtW4	CtW41	45,356,754	45.84	97.65	93.25	90.25	87.83
		CtW42	46,029,912	45.81	97.83	93.67	94.45	91.78
		CtW43	43,832,166	45.84	97.70	93.37	93.16	90.58

Table 2

Correlations between module membership and gene significance of modules.

Module	No. of Genes	Module Membership vs. Gene Significance
		Correlation (r2)	p-Value
Orange	997	−0.0094	0.77
Dark orange	526	−0.2	3.8 × 10−6
Blue	734	−0.65	2.4 × 10−89
Tan	270	−0.51	2.8 × 10−19
Violet	51	−0.36	0.0095
Dark olive green	121	0.027	0.77
Turquoise	230	0.04	0.55
Magenta	173	0.58	6.2 × 10−17
Grey	121	0.036	0.7

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f13101655/s1, Figure S1: The distribution of differently expressed transcription factor families in CsL (a), CtW (b), and CsL vs. CtW (c). CtW: cold-tolerant ‘Weixuan 1’; CsL: cold-sensitive ‘Longwangmao’. Figure S2: The expression analysis of DEGs in the magenta module. (A) Heatmap of DEGs in the magenta module. (B) The overall expression level of eigengenes identified in the magenta module for each sample. Table S1: List of primers used for qRT-PCR. Table S2: The antioxidant enzyme activities in apricot kernel pistils under freezing stress. Table S3: The annotation of hub genes in the magenta module.

References

1. Vitasse, Y.; Schneider, L.; Rixen, C.; Christen, D.; Rebetez, M. Increase in the risk of exposure of forest and fruit trees to spring frosts at higher elevations in Switzerland over the last four decades. Agric. For. Meteorol.; 2018; 248, pp. 60-69. [DOI: https://dx.doi.org/10.1016/j.agrformet.2017.09.005]

2. Yang, J.; Huo, Z.; Wang, P.; Wu, D.; Ma, Y.; Yao, S.; Dong, H. Process-based indicators for timely identification of apricot frost disaster on the warm temperate Zone, China. Theor. Appl. Climatol.; 2021; 1463, pp. 1143-1155. [DOI: https://dx.doi.org/10.1007/s00704-021-03762-2]

3. Chen, L.J.; Xiang, H.Z.; Miao, Y.; Zhang, L.; Guo, Z.F.; Zhao, X.H.; Lin, J.W.; Li, T.L. An overview of cold resistance in plants. J. Agron. Crop. Sci.; 2014; 200, pp. 237-245. [DOI: https://dx.doi.org/10.1111/jac.12082]

4. Augspurger, C.K. Reconstructing patterns of temperature, phenology, and frost damage over 124 years: Spring damage risk is increasing. Ecology; 2013; 94, pp. 41-50. [DOI: https://dx.doi.org/10.1890/12-0200.1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23600239]

5. Liu, Y.C.; Chen, W.Z.; Liu, W.S.; Liu, N.; Liu, S. Palynological study on the origin and systematic evolution of kernel-using apricots. Acta Hortic. Sin.; 2010; 37, pp. 1377-1387. [DOI: https://dx.doi.org/10.16420/j.issn.0513-353x.2010.09.002]

6. Ning, C.; Meng, Q.; Li, S.; Li, Y.; Yang, J. Comparative study of cold resistance in Kernel-apricot germplasm floral organs. J. Agric. Univ. Hebei; 2010; 33, pp. 37-41. [DOI: https://dx.doi.org/10.1016/S1002-0160(10)60014-8]

7. Yu, D.; Liu, X.; Cui, Y.; Bi, Q.; Zhao, Y.; Li, D.; Yu, H.; Wang, L. Comparative transcriptome combined with morpho-physiological analyses revealed candidate genes potentially for differential cold tolerance in two contrasting apricot Prunus armeniaca L. cultivars. Trees; 2020; 34, pp. 1205-1217. [DOI: https://dx.doi.org/10.1007/s00468-020-01991-y]

8. Shi, Y.; Ding, Y.; Yang, S. Molecular regulation of CBF signaling in cold acclimation. Trends Plant Sci.; 2018; 23, pp. 623-637. [DOI: https://dx.doi.org/10.1016/j.tplants.2018.04.002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29735429]

9. Thomashow, M.F. Plant cold acclimation: Freezing tolerance genes and regulatory mechanisms. Annu. Rev. Plant Biol.; 1999; 50, pp. 571-599. [DOI: https://dx.doi.org/10.1146/annurev.arplant.50.1.571] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15012220]

10. Ding, Y.; Shi, Y.; Yang, S. Advances and challenges in uncovering cold tolerance regulatory mechanisms in plants. New Phytol.; 2019; 222, pp. 1690-1704. [DOI: https://dx.doi.org/10.1111/nph.15696] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30664232]

11. Liu, C.; Cao, M.; Zhou, D.; Peng, W. Drought, salt and temperature stress-induced metabolic changes in plant. Biotechnol. Bull.; 2013; 4, pp. 1-7. [DOI: https://dx.doi.org/10.13560/j.cnki.biotech.bull.1985.2013.04.028]

12. Shi, Y.; Ding, Y.; Yang, S. Cold signal transduction and its interplay with phytohormones during cold acclimation. Plant Cell Physiol.; 2015; 561, pp. 7-15. [DOI: https://dx.doi.org/10.1093/pcp/pcu115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25189343]

13. Chinnusamy, V.; Ohta, M.; Kanrar, S.; Lee, B.H.; Hong, X.; Agarwal, M.; Zhu, J.K. ICE1: A regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev.; 2003; 17, pp. 1043-1054. [DOI: https://dx.doi.org/10.1101/gad.1077503]

14. Ding, Y.; Shi, Y.; Yang, S. Molecular regulation of plant responses to environmental temperatures. Mol Plant.; 2020; 13, pp. 544-564. [DOI: https://dx.doi.org/10.1016/j.molp.2020.02.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32068158]

15. Choudhury, F.K.; Rivero, R.M.; Blumwald, E.; Mittler, R. Reactive oxygen species, abiotic stress and stress combination. Plant J.; 2017; 90, pp. 856-867. [DOI: https://dx.doi.org/10.1111/tpj.13299] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27801967]

16. Hasanuzzaman, M.; Bhuyan, M.H.M.B.; Zulfiqar, F.; Raza, A.; Mohsin, S.M.; Mahmud, J.A.; Fujita, M.; Fotopoulos, V. Reactive oxygen species and antioxidant defense in plants under abiotic stress: Revisiting the crucial role of a universal defense regulator. Antioxidants; 2020; 9, 681. [DOI: https://dx.doi.org/10.3390/antiox9080681]

17. Wani, K.I.; Naeem, M.; Castroverde, C.D.M.; Kalaji, H.M.; Albaqami, M.; Aftab, T. Molecular mechanisms of nitric oxide NO signaling and reactive oxygen species ROS homeostasis during abiotic stresses in plants. Int. J. Mol. Sci.; 2021; 2217, 9656. [DOI: https://dx.doi.org/10.3390/ijms22179656]

18. Devireddy, A.R.; Zandalinas, S.I.; Fichman, Y.; Mittler, R. Integration of reactive oxygen species and hormone signaling during abiotic stress. Plant J.; 2021; 105, pp. 459-476. [DOI: https://dx.doi.org/10.1111/tpj.15010]

19. Zandalinas, S.I.; Fichman, Y.; Devireddy, A.R.; Sengupta, S.; Azad, R.K.; Mittler, R. Systemic signaling during abiotic stress combination in plants. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 13810-13820. [DOI: https://dx.doi.org/10.1073/pnas.2005077117]

20. Tian, J.; Wei, A.; Liu, Y.; Wang, S.H.; Liu, Y. Research progress on cold resistance of kernel apricot. J. Nor. For. Univ.; 2018; 331, pp. 174-178. [DOI: https://dx.doi.org/10.3969/j.issn.1001-7461.2018.01.28]

21. Kaya, O.; Kose, C.; Esitken, A.; Gecim, T.; Donderalp, V.; Taskin, S.; Turan, M. Frost tolerance in apricot Prunus armeniaca L. receptacle and pistil organs: How is the relationship among amino acids, minerals, and cell death points?. Int. J. Biometeorol.; 2021; 65, pp. 2157-2170. [DOI: https://dx.doi.org/10.1007/s00484-021-02178-x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34324064]

22. Yang, Y.; Yang, F.; Song, Y.; Li, J.; Lü, Z. Principal component analysis on cold resistance of different apricot varieties in Shanxi. J. Shanxi Agric. Sci.; 2020; 48, pp. 1913-1915. [DOI: https://dx.doi.org/10.3969/j.issn.1002-2481.2020.12.11]

23. Li, Y.; Gao, L.; Ren, S.; Fu, Y.; Yang, F.; Yang, J. A new cold-resistant kernel apricot variety “Waixuan1”. Acta Hortic. Sin.; 2010; 371, pp. 155-156. [DOI: https://dx.doi.org/10.16420/j.issn.0513-353x.2010.01.009]

24. Jiang, F.; Zhang, J.; Wang, S.; Yang, L.; Luo, Y.; Gao, S.; Zhang, M.; Wu, S.; Hu, S.; Wang, Y. The apricot (Prunus armeniaca L.) genome elucidates Rosaceae evolution and beta-carotenoid synthesis. Hortic. Res.; 2019; 6, 12. [DOI: https://dx.doi.org/10.1038/s41438-019-0215-6]

25. Yu, D.; Liu, X.; Bi, Q.; Zhao, Y.; Ju, J.; Li, D.; Yu, H.; Wang, L. Response of pistils in two contrasting apricot Prunus armeniaca L. cultivars to low temperature stress and recovery. J. Nor. For. Univ.; 2021; 491, pp. 1-5. [DOI: https://dx.doi.org/10.3969/j.issn.1000-5382.2021.01.001]

26. Li, K.; Du, Y.; Miao, Y. Future challenges in understanding ROS in plant responses to abiotic stress. Sci. China Life Sci.; 2016; 5912, pp. 1343-1344. [DOI: https://dx.doi.org/10.1007/s11427-016-0362-7]

27. Dreyer, A.; Dietz, K.J. Reactive oxygen species and the redox-regulatory network in cold stress acclimation. Antioxidants; 2018; 711, 169. [DOI: https://dx.doi.org/10.3390/antiox7110169]

28. Mehrotra, S.; Verma, S.; Kumar, S.; Kumari, S.; Mishra, B.N. Transcriptional regulation and signalling of cold stress response in plants: An overview of current understanding. Environ. Exp. Bot.; 2020; 180, 104243. [DOI: https://dx.doi.org/10.1016/j.envexpbot.2020.104243]

29. Maurya, N.K.; Goswami, A.K.; Singh, S.K.; Prakash, J.; Goswami, S.; Chinnusamy, V.; Talukdar, A.; Pradhan, S.; Kumari, A. Studies on expression of CBF1 and CBF2 genes and anti-oxidant enzyme activities in papaya genotypes exposed to low temperature stress. Sci. Horticult.; 2020; 261, 108914. [DOI: https://dx.doi.org/10.1016/j.scienta.2019.108914]

30. Ritonga, F.N.; Ngatia, J.N.; Wang, Y.R.; Khoso, M.A.; Farooq, U.; Chen, S. AP2/ERF, an important cold stress-related transcription factor family in plants: A review. Physiol. Mol. Biol. Plants; 2021; 279, pp. 1953-1968. [DOI: https://dx.doi.org/10.1007/s12298-021-01061-8]

31. Agarwal, M.; Hao, Y.; Kapoor, A.; Dong, C.H.; Fujii, H.; Zheng, X.; Zhu, J.K. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing tolerance. J. Biol. Chem.; 2006; 281, pp. 37636-37645. [DOI: https://dx.doi.org/10.1074/jbc.M605895200] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17015446]

32. Hao, Y.J.; Wei, W.; Song, Q.X.; Chen, H.W.; Zhang, Y.Q.; Wang, F.; Zou, H.F.; Gang, L.; Tian, A.G.; Zhang, W.K. Soybean NAC transcription factors promote abiotic stress tolerance and lateral root formation in transgenic plants. Plant J.; 2011; 682, pp. 302-313. [DOI: https://dx.doi.org/10.1111/j.1365-313X.2011.04687.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21707801]

33. Yarra, R.; Wei, W. The NAC-type transcription factor GmNAC20 improves cold, salinity tolerance, and lateral root formation in transgenic rice plants. Funct. Integr. Genom.; 2021; 213, pp. 473-487. [DOI: https://dx.doi.org/10.1007/s10142-021-00790-z] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34191184]

34. Chen, X.; Ding, Y.; Yang, Y.; Song, C.; Wang, B.; Yang, S.; Guo, Y.; Gong, Z. Protein kinases in plant responses to drought, salt, and cold stress. J. Integr. Plant Biol.; 2021; 63, pp. 53-78. [DOI: https://dx.doi.org/10.1111/jipb.13061] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33399265]

35. Ding, Y.; Li, H.; Zhang, X.; Xie, Q.; Gong, Z.; Yang, S. OST1 kinase modulates freezing tolerance by enhancing ICE1 stability in Arabidopsis. Dev. Cell; 2015; 32, pp. 278-289. [DOI: https://dx.doi.org/10.1016/j.devcel.2014.12.023]

36. Liu, Y.; Xu, C.; Zhu, Y.; Zhang, L.; Chen, T.; Zhou, F.; Chen, H.; Lin, Y. The calcium-dependent kinase OsCPK24 functions in cold stress responses in rice. J. Integr. Plant Biol.; 2018; 60, pp. 95-110. [DOI: https://dx.doi.org/10.1111/jipb.12614]

37. Kim, S.J.; Kim, W.T. Suppression of Arabidopsis RING E3 ubiquitin ligase AtATL78 increases tolerance to cold stress and decreases tolerance to drought stress. FEBS Lett.; 2013; 587, pp. 2584-2590. [DOI: https://dx.doi.org/10.1016/j.febslet.2013.06.038]

38. Suh, J.Y.; Kim, W.T. Arabidopsis RING E3 ubiquitin ligase AtATL80 is negatively involved in phosphate mobilization and cold stress response in sufficient phosphate growth conditions. Biochem. Biophys. Res. Commun.; 2015; 4634, pp. 793-799. [DOI: https://dx.doi.org/10.1016/j.bbrc.2015.06.015]

39. Wang, X.; Ding, Y.L.; Li, Z.Y.; Shi, Y.T.; Wang, J.L.; Hua, J.; Gong, Z.Z.; Zhou, J.M.; Yang, S.H. PUB25 and PUB26 promote plant freezing tolerance by degrading the cold signaling negative regulator MYB15. Dev. Cell; 2019; 51, pp. 222-235. [DOI: https://dx.doi.org/10.1016/j.devcel.2019.08.008]

40. Fang, H.; Meng, Q.; Xu, J.; Tang, H.; Tang, S.; Zhang, H.; Huang, J. Knock-down of stress inducible OsSRFP1 encoding an E3 ubiquitin ligase with transcriptional activation activity confers abiotic stress tolerance through enhancing antioxidant protection in rice. Plant Mol. Biol.; 2015; 874–875, pp. 441-458. [DOI: https://dx.doi.org/10.1007/s11103-015-0294-1]

41. Doherty, C.J.; Van Buskirk, H.A.; Myers, S.J.; Thomashow, M.F. Roles for Arabidopsis CAMTA transcription factors in cold-regulated gene expression and freezing tolerance. Plant Cell; 2009; 21, pp. 972-984. [DOI: https://dx.doi.org/10.1105/tpc.108.063958] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19270186]

42. Fang, P.; Yan, M.; Chi, C.; Wang, M.; Zhou, Y.; Zhou, J.; Shi, K.; Xia, X.; Foyer, C.H.; Yu, J. Brassinosteroids act as a positive regulator of photoprotection in response to chilling stress. Plant Physiol.; 2019; 180, pp. 2061-2076. [DOI: https://dx.doi.org/10.1104/pp.19.00088] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31189657]

43. Cooper, J.W.; Hu, Y.; Beyyoudh, L.; Yildiz Dasgan, H.; Kunert, K.; Beveridge, C.A.; Foyer, C.H. Strigolactones positively regulate chilling tolerance in pea and in Arabidopsis. Plant Cell Environ.; 2018; 41, pp. 1298-1310. [DOI: https://dx.doi.org/10.1111/pce.13147] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29341173]

44. Wang, M.; Dai, W.; Du, J.; Ming, R.; Dahro, B.; Liu, J.H. ERF109 of trifoliate orange Poncirus trifoliata L. Raf. contributes to cold tolerance by directly regulating expression of Prx1 involved in antioxidative process. Plant Biotechnol. J.; 2019; 17, pp. 1316-1332. [DOI: https://dx.doi.org/10.1111/pbi.13056] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30575255]