Next Article in Journal
Positive and Negative Work Practices of Forest Machine Operators: Interviews and Literature Analysis
Next Article in Special Issue
Transcriptome Analysis Reveals the Molecular Mechanisms Associated with Flower Color Formation in Camellia japonica ‘Joy Kendrick’
Previous Article in Journal
Integrating Social Forestry and Biodiversity Conservation in Indonesia
Previous Article in Special Issue
Transcriptome Analysis Provides Insights into the Mechanisms of Starch Biosynthesis in the Kernels of Three Chestnut Cultivars
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptome and Expression Analysis of Genes Related to Regulatory Mechanisms in Holly (Ilex dabieshanensis) under Cold Stress

1
Jiangsu Key Laboratory for the Research and Utilization of Plant Resources, Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, Nanjing 210014, China
2
Fuyang Academy of Agricultural Sciences, Fuyang 236065, China
3
Zhejiang Provincial Key Laboratory of Forest Aromatic Plants-Based Healthcare Functions, Zhejiang A & F University, Hangzhou 311300, China
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Forests 2022, 13(12), 2150; https://doi.org/10.3390/f13122150
Submission received: 16 October 2022 / Revised: 11 December 2022 / Accepted: 13 December 2022 / Published: 15 December 2022
(This article belongs to the Special Issue Genetic Regulation of Growth and Development of Woody Plants)

Abstract

:
Ilex dabieshanensis (K. Yao and M. B. Deng) is not only an important economic tree species, but also has the characteristics of evergreens in all seasons, as well as strong cold resistance. In order to understand the molecular mechanism of holly’s response to cold stress, we used transcriptome analysis to identify the main signaling pathways and key genes involved in cold stress. The result showed that 5750 differentially expressed genes (DEGs) were identified under different cold treatment times compared with the control (cold—0 h). The Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs showed that seven phytohormone signal transduction were the most highly enriched, including abscisic acid (ABA), ethylene (ET), cytokinin (CK), auxin (IAA), gibberellin (GA), jasmonate (JA), and brassinosteroids (BR). In addition, proline metabolism, arginine metabolism, flavonoid biosynthesis, and anthocyanin biosynthesis were also implicated in response to cold stress. The weighted gene co-expression network analysis (WGCNA) showed that the genes in two modules were significantly up-regulated after 12 h and 24 h treatments, suggesting these two module genes may participate in the cold stress. The gene ontology (GO) results of the two module genes showed that calcium, scavenging reactive oxygen species, and nitric oxide might act as signaling molecules to regulate cold tolerance in holly. By calculating the connectivity and function prediction of genes in the two modules, five genes (evm.TU.CHR2.244, evm.TU.CHR1.1507, evm.TU.CHR1.1821, evm.TU.CHR2.89, and evm.TU.CHR2.210) were identified as the key hub genes of I. dabieshanensis response to cold stress. These results provided candidate genes and clues for further studies on the molecular genetic mechanism of cold stress in holly.

1. Introduction

Ilex L. (holly) is the largest woody dioecious genus of angiosperms, with approximately 600 species and wide geographical distribution. Many holly species are economically important trees used in various ways, including manufacturing beverages, medicines, and wood [1,2,3]. Ilex dabieshanensis (K. Yao and M. B. Deng) (I. dabieshanensis) is a species of Ilex, which was first discovered and named by the Institute of Botany, Chinese Academy of Sciences in Jiangsu Province. This plant is a precious native tree species with multiple uses in medicine, ornamental landscaping, and tea. I. dabieshanensis also has unique ecological characteristics, such as pruning resistance, wide adaptability, cold resistance (Figure 1), pollution resistance, and fewer pests and diseases [4]. Most of the research conducted on holly involves processed products, extract functions, origin and distribution, physiological characteristics, cultivation methods, and ornamental applications [5,6,7,8,9]. A recent article has reported the whole genome of I. polyneura, which provides a necessary basis for advancing research on Ilex to the molecular level [7]. As a perennial woody plant, I. dabieshanensis undergoes seasonal changes and is evergreen throughout the year. Due to the cyclical changes in seasons, plants are inevitably affected by the adversity of low temperatures. Low temperatures are one of the most harmful environmental stressors encountered by vascular plants [10,11]. Throughout long-term evolution, plants have gradually developed unique biological characteristics, physiological functions, and other adaptive mechanisms, which have given them the ability to adapt to and resist low-temperature conditions. These adaptations vary by species, geographic distribution, and even developmental stage. Each plant has different enrichment pathways during different cold stress periods, and every gene has a unique mechanism in response to cold stress [12]. Currently, there are no reports on the cold tolerance mechanism of holly.
When plants are subjected to low-temperature stress, they undergo a series of physiological and biochemical reactions and gene expression regulation from sensing low-temperature signals, which produce cold resistance. At the physiological level, many substances or protective proteins are synthesized in the plant, such as soluble sugars, proline, and cold resistance proteins [13,14]. These substances or proteins regulate osmotic potential, ice crystal formation, cell membrane stability, and scavenging reactive oxygen species (ROS) [15,16]. At the molecular level, according to the function of low temperature-induced expression genes, genes can be divided into two categories: regulatory genes, which play a role in regulating gene expression and signal transduction, including transcription factors (TFs) such as C-repeat binding factors (CBF) and protein kinases that sense and transmit signals such as mitogen-activated protein kinase (MAPKs) [17,18,19,20,21]; and cold-regulated (COR) genes, which are early responsive to dehydration (ERD) genes, responsive to dehydration (ED) genes, low-temperature induced (LTI) genes, and cold inducible (KIN) genes [22,23]. Previous studies found that the inducer of CBF expression (ICE)-CBF-COR pathway plays a crucial role in plant cold tolerance and is one of the most widely reported pathways [24]. In most plant species, the ICE-CBF-COR pathway is induced by cold stress, activating the expression of downstream genes encoding osmotic regulators [25]. In detail, ICE1 can be released from DELLA-binding Jasmonate ZIM-domain (JAZ) to induce the expression of CBF3 under low-temperature stress. CBFs, also known as dehydration-responsive element-binding proteins (DREBs), are members of the APETALA2/ethylene-responsive factor (AP2/ERF) family that activates cold-responsive (COR) gene expression by binding to cis-elements on the COR gene promoter [26]. Winter barley has a vernalization gene (VRN1) that can synergize with the CBF gene, inducing a stronger cold resistance than spring barley [27]. ICE acts upstream to induce and regulate the expression of CBF [28,29,30]. Two homologs of the ICE gene (ICE1 and ICE2) have been characterized in multiple plant species, and their role in cold tolerance has been elucidated [31]. In Arabidopsis, many COR genes are activated or repressed by the action of CBF1/2/3 and have been shown to directly or indirectly increase cold stress tolerance in plants [32]. Many TFs that regulate cold signaling and stress have been identified, such as CBF1-3, ICE1-2, CAMTA3, MYB15, COR15a, COR15b, etc. [33]. In addition, some protein kinases (MEKK1-MKK1/2-MPK4) were found to induce the expression of CBF, especially CBF2, which is involved in the mitogen-activated protein kinase (MAPK) cascade [34,35]. In detail, cold stress receptor kinase-like 1 (CRKL1) is regulated by calcium/calmodulin binding, and CRLK1/2 interacts with MEKK1, a MAPK module that responds to lower temperatures [19,20]. MEKK1 then phosphorylates MKK2, which activates MPK4/6, forming an upstream pathway, CRLK1-MEKK1-MKK2-MPK4-MPK3/6, that enhances CBF gene expression [21,36].
Plant hormones play an important role in plant growth and development and protect plants from abiotic stresses, such as cold stress. It has been reported that many hormones regulate cold stress by regulating CBF. One of the most important hormones in cold stress is abscisic acid (ABA) [37]. Previous studies have shown that ABA can induce the expression of CBF and COR genes [38]. Gibberellin (GA) plays a powerful role in the ICE-CBF-COR pathway [39]. It is also understood that anthocyanin content in plants is higher under cold stress conditions. Brassica napus transformed with the AtGA2ox8 gene has a reduced amount of bioactive GA and can produce more anthocyanins in winter [40]. Ethylene (ET), ABA, and jasmonate (JA) induce the expression of ethylene-responsive factor (ERF) genes. ERFs can bind to the GCC box and DRE elements at low temperatures to enhance plant tolerance to cold stress [41]. In addition, methyl jasmonate can induce the expression of guaiacol peroxidase (POD), catalase (CAT), and ascorbate peroxidase (APX)-related genes that reduce reactive oxygen species, enhancing bell pepper tolerance to low-temperature stress [33]. Two brassinosteroids (BR)-responsive transcription factors, Brassinazole-resistant 1 (BZR1) and CES (CESTA), are also direct regulators of CBF [39,42].
Due to these differences in signal transduction pathways and metabolisms, the plant mechanisms associated with cold stress are complex [43]. Studying transcriptional changes in plants during cold stress is crucial for understanding and identifying plants’ molecular mechanisms under cold stress. In this study, we designed cold treatments at different time markers to investigate the transcriptome changes of I. dabieshanensis. A comprehensive comparative analysis identified the major signaling pathways and key genes involved or shared by cold stress. These results broadened our understanding of the molecular mechanisms by which holly responds to low temperatures and provided information for further studies on key gene functions.

2. Materials and Methods

2.1. Plant Materials and Experimental Design

The research material is the annual cutting seedlings of I. dabieshanensis, and the cuttings are all derived from the same mother plant. This ensures the consistency of the genetic background and physiological state of the experimental material. The seedlings with consistent growth (about 10 leaves) were selected and placed in an artificial climate room at 25 °C for 4 weeks. After 4 weeks, they were treated at a low temperature of 4 °C, and samples were taken at 0 h, 3 h, 6 h, 9 h, 12 h, and 24 h, respectively [44], and the leaves of 6 plants were mixed as a sample in each time period, and three biological replicates were set for each treatment. Samples were immediately stored in liquid nitrogen. A part was directly used for transcriptome sequencing, and the rest was stored in a −80 °C ultra-low temperature freezer for later use.

2.2. Total RNA Extraction and mRNA Libraries Construction

Total RNA was extracted from leaves using a plant RNA purification reagent. We used the Illumina TruseqTM RNA sample prep Kit method to construct the library. Detailed experimental steps are as follows. The starting amount of 1 μg total RNA was used to build the library; after the mRNA was isolated by magnetic beads, the mRNA was interrupted by ions (TruseqTM RNA sample prep Kit); double-stranded cDNA synthesis, filling, adding A at the 3′ end, and connecting the index adapter (TruseqTM RNA sample prep Kit) were performed; library enrichment and PCR amplification for 15 cycles were performed; 2% agarose gel recovery of the target band (Certified Low Range Ultra Agarose); TBS380 (Picogreen) quantification was performed, mixed according to the data ratio; cBot Bridge PCR amplification was performed on the top to generate clusters; and the Illumina Novaseq sequencing platform was used for 2×150 bp sequencing.

2.3. mRNA Sequence Data Processing

The original image data obtained by Illumina sequencing was converted into sequence data through base calling, and the results are stored in the FASTQ file format. The FASTQ file is the most original data file, and the file contains the sequence information of the sequencing read and the sequencing quality information. The raw sequencing data were first filtered by using trimmomatic (ILLUMINACLIP:adapters.fa:2:30:10 SLIDINGWINDOW:4:15 MINLEN:75) (http://www.usadellab.org/cms/index.php?page= trimmomatic) (accessed on 12 August 2021), so as to obtain high-quality sequencing data (clean data) to ensure the smooth progress of subsequent analysis (Bioproject accession number: PRJNA897070). Next, the high-quality sequencing sequences obtained after quality control were compared with the designated I. polyneura reference genome (https://ngdc.cncb.ac.cn/search/?dbId=gwh&q=GWHBDNW00000000&page=1) (accessed on 7 April 2022), using Hisat2 (http://://ccb.jhu.edu/software/hisat2/faq.shtml) (accessed on 7 April 2022) with default parameters. The FPKM value (fragments per kilobase of exon model per million mapped reads) was calculated using featureCounts (http://subread.sourceforge.net/) (accessed on 7 April 2022).

2.4. Screening of Differentially Expressed Genes

Use edgeR software was used to analyze differentially expressed genes, and it was used to compare different cold treatment groups in pairs to obtain Cold-3 h vs. Cold-ck (0 h), Cold-6 h vs. Cold-ck, Cold-9 h vs. Cold-ck, Cold-12 h vs. Cold-ck, Cold-24 h vs.Cold-ck, Cold-6 h vs. Cold-3 h, Cold-9 h vs. Cold-3 h, Cold-12 h vs.Cold-3 h, Cold-24 h vs. Cold-3 h, Cold-9 h vs. Cold-6 h, Cold-12 h vs. Cold-6 h, Cold-24 h vs. Cold-6 h, Cold-12 h vs. Cold-9 h, Cold-24 h vs. Cold-9 h, and Cold-24 h vs. Cold-12 h as differentially expressed genes. The fold change (FC) and false discovery rate (FDR) were used as the screening conditions for differentially expressed genes, and the screening threshold was FDR ≤ 0.05 and FC ≥ 2. The differential expression was analyzed visually.

2.5. Differentially Expressed Gene Enrichment Analysis

The software Goatools (https://github.com/tanghaibao/GOatools) (accessed on 7 April 2022) was used for GO enrichment analysis, and the method was Fisher’s exact test. The p-values were corrected using the FDR (false discovery rate) multiple testing method to control for the calculated false positive rate. Typically, this GO function was considered significantly enriched when the p-value ≤ 0.05. Generally, the top 30 results of the GO enrichment analysis were selected as the main nodes of the directed acyclic graph, and the associated GO term was displayed together through the inclusion relationship. The depth of the color represents the degree of enrichment. In our project, we drew DAG diagrams of biological process, molecular function, and cellular component, respectively.
The software KOBAS (http://kobas.cbi.pku.edu.cn/kobas3/?t=1) (accessed on 7 April 2022) was used for KEGGPATHWAY enrichment analysis. The calculation principle was the same as that of GO functional enrichment analysis, and Fisher’s exact test was used for calculation. In order to control the calculated false positive rate, the BH (FDR) method was used for multiple testing. The calculation formula was the same as the previous section, and the p value was 0.05 as the threshold. The KEGG pathway that met this condition was defined as the KEGG pathway that was significantly enriched in differentially expressed genes.

2.6. Construction of Weighted Gene Co-Expression Network

The WGCNA algorithm is a common algorithm for constructing gene co-expression networks, and the R language package [45] was used for analysis. The WGCNA algorithm first assumes that the gene network obeys a scale-free distribution and defines the gene co-expression correlation matrix and the adjacency function formed by the gene network, and then it calculates the dissimilarity coefficients of different nodes and builds a hierarchical clustering tree accordingly. Different branches of the clustering tree represent different gene modules. To identify biologically meaningful modules, module features are used to calculate correlation coefficients. Principal component analysis (PCA) was performed on all genes in each co-expression module, and the principal component 1 (PC1) was called the eigengene of this module (Module eigengene, ME) in order to screen the cold stress-related specificity module, as well as to calculate the correlation coefficient r and p value of ME values of each module and different traits (here, different cold stress time). In this study, a module was considered to be a specific module if the correlation coefficient between its ME value and the trait |r| > 0.70 and p < 0.05. The top 30 genes connected in the module were selected as hub genes (highly connected genes). In this study, the connection relationship with edge weight > 0.3 was selected, and the hub gene interaction network was visualized using Cytoscape_3.3.0.

2.7. Validation of Gene Expression by qRT-PCR

To verify the RNA-seq results, we selected five DEGs for quantitative reverse transcription PCR (RT-qPCR) analysis. Total RNA of 18 samples was isolated by Trizol reagent following the protocol. The cDNA was generated by Primerscript RT reagent Kit with gDNA Erase (Takara), according to the manufacturer’s protocol. Quantitative PCR (qPCR) was performed using SYBR Premix Ex Tag (Takara), with actin as an internal reference gene. The relative expression levels of five DEGs normalized to the expression level of the internal reference control were calculated using the 2−∆∆Ct method [46]. The primers are listed in Table S1.

3. Results

3.1. Differential Analysis of Gene Expression in I. dabieshanensis under Cold Stress

Taking the genome annotation file as a reference, the FPKM value of each transcript in each sample was counted, and this value was used as the transcript’s expression level. The significant differences in expression were analyzed for different cold treatments [0 h (h), 3 h, 6 h, 9 h, 12 h, and 24 h] to study the gene expression of I. dabieshanensis under cold stress. The significantly different gene (DEG) screening conditions were FDR ≤ 0.05 and FC ≥ 2. A total of 5750 DEGs were obtained, and differential expression was visualized (Table S2). Compared to no cold stress (CK), gene expression was significantly up-regulated or down-regulated during different cold stress time markers. As the period of stress increased, the number of up-regulated DEGs and down-regulated DEGs increased (Figure 2a–e). Compared to 3 h, the number of DEGs under stress for 6 h, 9 h, 12 h, and 24 h were 241, 555, 1859, and 2228, respectively (Figure 2f–i). Compared to 6 h, the number of DEGs under stress for 9 h, 12 h, and 24 h were 238, 1737, and 1937, respectively (Figure 2j–l). Compared to 9 h, the number of DEGs under stress for 12 h and 24 h was 1330 and 1636, respectively (Figure 2m,n). Compared to 12 h, the number of DEGs under stress for 24 h was 1566 (Figure 2o). These results indicated that the expression of some genes was significantly different between different stress treatments, and the number of DEGs showed an upward trend as the treatment time gap increased.

3.2. GO Enrichment Analysis of DEGs

To fully understand the gene function classification of holly in response to cold stress, GO annotation was first performed on DEGs. GO annotation includes three aspects: biological process (BP), cellular component (CC), and molecular function (MF). The results showed that, compared to 0 h, DEGs under different cold treatments were annotated in BP, CC, and MF, with differences. From the analysis of level 2, the DEGs of 3 h vs. 0 h were mainly annotated to 8, 3, and 17 GO terms for MF, CC, and BP, respectively (Figure S1). The DEGs of 6 h vs. 0 h were mainly annotated to 8, 3, and 19 GO terms for MF, CC, and BP, respectively (Figure S1). The DEGs of 9 h vs. 0 h were mainly annotated to 8, 3, and 19 GO terms for MF, CC, and BP, respectively (Figure S1). The DEGs of 12 h vs. 0 h were mainly annotated to 9, 3, and 18 GO terms for MF, CC, and BP, respectively (Figure S1). The DEGs of 24 h vs. 0 h were mainly annotated to 8, 3, and 20 GO terms for MF, CC, and BP, respectively (Figure S1). Although the GO terms were different, they were mainly annotated to catalytic activity, intracellular entities, cellular anatomical entities, cellular processes, biological regulation, response to stimuli, and metabolic processes (Figure S1).
Canonically, GO functions are considered to be significantly enriched when the p-value ≤ 0.05. According to the GO enrichment analysis, compared to 0 h, the number of GO-enriched terms of DEGs under stress for 3 h, 6 h, 9 h, 12 h, and 24 h were 160 (BP: 124, CC: 17, MF: 19), 240 (BP: 184, CC: 17, MF: 39), 239 (BP: 197, CC: 10, MF: 32), 329 (BP: 271, CC: 6, MF: 52), and 319 (BP: 224, CC: 10, MF:76), respectively (Table 1). With prolonged cold-treatment time, the terms of DEGs related to BP and MF were increasingly abundant, while the terms related to CC decreased.
In terms of MF, DEGs under cold stress at different times were significantly enriched in “DNA-binding transcription factor activity” (Figure 3) and “calcium ion binding” (Tables S2–S7). According to statistics, there are 425 transcription factors in DEGs, which are distributed in 42 transcription factor families. The larger transcription factor families are ERF (58), MYB (33), NAC (32), bHLH (28), WRKY (28), C2 H2 (24), GRAS (21), MYB_ Related (21), HD-ZIP (14), G2-like (13), Trihelix (11), bZIP (10), Dof (10), HSF (10), etc. (Figure S2). In addition, the activities of some enzymes may be affected differently in different stress stages (Tables S3–S7). For example, “galactosidase activity” and “trehalose-phosphatase activity” were mainly affected in the early cold treatment stage, while “ubiquitin-protein transferase activity” and “UDP-glucosyltransferase activity” were mainly affected in the later stages. In terms of CC, GO terms are few, and DEGs were mainly enriched in “plant-type cell wall”, “plasmodesma”, “amyloplast”, “intrinsic component of plasma membrane”, and “anchored component of membrane” (Figure 3a–c).
The most abundant GO terms for BP were cold-treatment-induced processes related to hormone metabolism and signal transduction pathways. Under cold stress at different time markers, DEGs were significantly enriched in “response to hormone” (Figure 3). Some DEGs under cold stress were also significantly enriched in “hormone transport”, “hormone metabolic process”, “Regulation of hormone levels”, “regulation of hormone metabolic process”, and “regulation of hormone biosynthetic process” (Tables S3–S7). In more detail, cold-induced DEGs were significantly enriched in a variety of hormone-related terms. For example, cold-induced DEGs were involved in ten IAA-related GO terms, including “response to auxin”, “auxin-activated signaling pathway”, “auxin transport”, “auxin homeostasis”, and “auxin efflux”, etc. (Tables S3–S7). Six ABA-related GO terms were significantly enriched, including “response to abscisic acid”, “abscisic acid-activated signaling pathway”, and “abscisic acid metabolic process”, etc. (Tables S3–S7). Four JA-related GO terms were significantly enriched, including “response to jasmonic acid”, “jasmonic acid mediated signaling pathway”, and “regulation of jasmonic acid mediated signaling pathway”, etc. (Tables S3–S7). Two ET-related GO terms were significantly enriched, including “response to ethylene” and “ethylene-activated signaling pathway” (Tables S3–S7). Two GA-related GO terms were significantly enriched, including “response to gibberellin” and “regulation of gibberellic acid mediated signaling pathway”) (Tables S3–S7). Two CK-related GO terms were significantly enriched, including “response to cytokine” and “cytokinein-activated signaling pathway” (Tables S3–S7). Three salicylic acid-related GO terms were significantly enriched, including “response to salicylic acid”, “regulation of salicylic acid mediated signaling pathway”, and “regulation of salicylic acid metabolic process” (Tables S3–S7). Additionally, three BR-related GO terms were significantly enriched, including “response to brassinosteroid”, “brassinosteroid biosynthetic process”, and “brassinosteroid homeostasis” (Tables S3–S7). In addition, “response to larrikin”, “response to water”, “response to temperature stimulus”, “response to cold”, and “response to chitin” were also highly enriched (Figure 3b–e).

3.3. KEGG Enrichment Analysis of DEGs

The results of the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis showed that the DEGs induced by cold treatment were significantly enriched in various metabolic pathways, including carbohydrate metabolism, biosynthesis of other secondary metabolites, energy metabolism, amino acid metabolism, xenobiotics biodegradation and metabolism, lipid metabolism, glycan biosynthesis and metabolism, and metabolism of terpenoids and polyketides (Tables S8–S12). Among them, “Pentose and glucuronate interconversions” and “Starch and sucrose metabolism” were significantly enriched in the early stages of cold stress, but not during later stages. In contrast, some metabolic pathways were significantly enriched in the middle and late stages of cold stress and not significantly enriched in earlier periods, such as “Arginine and proline metabolism”, “Carotenoid biosynthesis”, “Cysteine and methionine metabolism”, “Galactose metabolism”, “Inositol phosphate metabolism”, “Photosynthesis—antenna proteins”, “Fatty acid elongation”, and “Flavone and flavonol biosynthesis” (Tables S8–S12).
In terms of cellular processes, DEGs were significantly enriched in “Cellular community—eukaryotes”, “involving Gap junction” and “Tight junction” (Figure 4a–d). In organismal systems, DEGs were significantly enriched in the “Toll-like receptor signaling pathway”, “NOD-like receptor signaling pathway”, “Toll and Imd signaling pathway”, “Plant-pathogen interaction”, “Mineral absorption”, and “Neurotrophin signaling pathway” (Figure 4). In organismal systems, cold stress only induced significant enrichment of “Mineral absorption” (Figure 4e) and “Plant-pathogen interaction” (Figure 4d,e) in the middle and late stages.
Notably, signal transduction pathways in environmental information processing were significantly enriched, including “Plant hormone signal transduction”, “MAPK signaling pathway”, “NF-kappa B signaling pathway”, and “Two-component system” (Figure 4). In particular, “Plant hormone signal transduction” was a significantly enriched pathway of DEGs in each cold stress stage and was also the most significantly enriched pathway, except for during the 3 h time marker (Figure 4). The “MAPK signaling pathway” and the “NF-kappa B signaling pathway” were significantly enriched in four of the cold stress stages (Figure 4a–e).

3.4. The Expression Patterns of DEGs Response to Cold Tolerance in I. dabieshanensis

To understand the gene expression profile of I. dabieshanensis in response to cold stress, we used Delight Hwarari’s figure [24] as a reference to show the expression pattern of DEGs under cold stress. The results showed that various hormone signaling pathways, MAPK signaling pathways, and some “cold-responsive” genes played essential roles during cold stress in holly. We found that seven hormone (IAA, ABA, CK, ET, GA, BR, and JA) signal transduction-related genes were up-regulated or down-regulated to varying degrees by cold stress (Figure 5, Table S13). Some MAPK signaling-related genes were up-regulated, such as four serine/threonine protein kinase genes (evm.TU.CHR18.1145, evm.TU.CHR2.571, evm.TU.CHR8.589, and evm.TU.CHR9.656), one HSP gene (evm.TU.CHR2.705), and one mitogen-activated protein kinase gene (evm.TU.CHR13.25) (Figure 5, Tables S14 and S15). The ICE-CBF-COR pathway-related gene ICE1 (evm.TU.CHR20.178) gene and two CBF genes (evm.TU.CHR4.1892 and evm.TU.CHR8.33) were shown to be up-regulated by cold stress (Figure 5, Table S16). Among some “response to cold” genes, one ERF gene (evm.TU.CHR11.908) and one peroxidase gene (evm.TU.CHR17.727) were up-regulated (Figure 5, Table S15). In contrast, more “response to cold” genes were significantly down-regulated, such as C-repeat/DRE binding factor 1 (evm.TU.CHR19.317), dehydrin, and dehydration-responsive element-binding protein 1D (Figure 5, Table S15).

3.5. Construction and Analysis of Weighted Gene Co-Expression Network

We performed a weighted gene co-expression network analysis (WGCNA) to identify genes associated with cold stress. All the DEGs obtained by cold treatment of holly during 0 h (CK), 2 h, 6 h, 12 h, and 24 h were selected, and the WGCNA package in the R software environment was used to construct the weighted gene co-expression network. According to Figure S3, a soft threshold β = 6 was chosen to construct the co-expression network. The modules were divided by the dynamic shearing method and merged into 11 co-expression modules by choosing a cutting height of 0.25. Different colors represented different modules, and the grey modules represented genes that could not be classified into modules (Figure 6a). Two modules were significantly associated with cold traits (showing regular changes with increasing cold treatment time). The black module (r = 0.86, p = 5 × 10−6) was significantly positively correlated with cold stress for 12 h, and the blue module (r = 0.91, p = 1 × 10−7) was significantly positively correlated with cold stress for 24 h (Figure 6b). We found that black and blue modules’ gene expression profiles significantly differed between cold-stressed and unstressed samples. The genes of the black module were significantly up-regulated by cold stress for 12 h, while the genes of the blue module were significantly up-regulated by cold stress for 24 h (Figure S4). These results suggested that genes in these two modules may trigger and participate in the cold mediation process within the plant.
To find the key genes involved in low-temperature stress, GO and KEGG enrichment analyses were performed on the genes identified in the black and blue modules. GO enrichment analysis showed that some GO terms were common to the black and blue modules and were arranged in the front, such as “response to acid chemical”, “response to oxygen-containing compound”, “response to organic substance”, and “response to acid chemical”, etc. (Figure 7a,c). In addition, some GO terms were also significantly enriched in the black and blue modules, such as “response to cold”, “calcium ion binding”, and some hormone-related terms (“response to abscisic acid”, “cellular response to jasmonic acid stimulus”, and “response to ethylene”) (Tables S17 and S18). Some GO terms, such as “UDP-glycosyltransferase activity”, “UDP-glucosyltransferase activity”, “glucosyltransferase activity”, “quercetin 3-O-glucosyltransferase activity”, and “quercetin 7-O-glucosyltransferase activity”, were unique to the blue module and had a high degree of enrichment (Figure 7a,c). KEGG pathway enrichment analysis showed that “Plant hormone signal transduction” and “Plant-pathogen interaction” were significantly enriched pathways for both black and blue modules (Figure 7b,d). “Toll-like receptor signaling pathway”, “NF-kappa B signaling pathway”, “Toll and Imd signaling pathway”, “MAPK signaling pathway”, and “NOD-like receptor signaling pathway” were significantly enriched in the black module, but were not significantly enriched in the blue module (Tables S19 and S20). Conversely, some pathways were significantly enriched in the blue module but not significantly or not enriched in the black module, such as “Fatty acid elongation”, “Flavone and flavonol biosynthesis”, and “Flavonoid biosynthesis” (Tables S19 and S20).
Next, the highly connected genes of the black and blue modules were used as hub genes to construct gene co-expression networks for holly cold stress, in which each node represented a gene and the lines between genes represented the co-expression relationship. The internal connectivity of a gene in a module reflected the network location and importance of the gene to the module. In the hub gene co-expression network of the black module, the connectivity and edge weight of nine genes were high, such as evm.TU.CHR2.70, evm.TU.CHR1.2647, evm.TU.CHR1.2184, evm.TU.CHR1.1707, evm.TU.CHR2.566, evm.TU.CHR2.244, evm.TU.CHR2.177, evm.TU.CHR1.1507, and evm.TU.CHR2.210 (Figure 8a). Through functional annotation of the genes, we found that evm.TU.CHR2.244 and evm.TU.CHR1.1507 were transcription factors (Table S21). In contrast, evm.TU.CHR2.227, a transcription factor, had low connectivity and edge weight values in this network. Additionally, three genes (evm.TU.CHR1.1351, evm.TU.CHR1.1079, and evm.TU.CHR2.696) had higher connectivity but lower edge weight values (Figure 8a). Eleven genes (evm.TU.CHR1.1761, evm.TU.CHR1.2160, evm.TU.CHR2.63, evm.TU.CHR1.1762, evm.TU.CHR2.89, evm.TU.CHR1.1567, evm.TU.CHR1.2325, evm.TU.CHR2.974, evm.TU.CHR1.30, evm.TU.CHR1.2240, and evm.TU.CHR1.2003) in the blue module had high connectivity and edge weights (Figure 8b), and evm.TU.CHR1.2003 was functionally annotated as a transcription factor (Table S22). Another transcription factor, evm.TU.CHR1.1821, was found to have the highest number of connections in this network, but the edge weights were all low (Figure 8b, Table S22). The identified hub genes had close interactions or indirect cross-regulation, suggesting that these genes may play crucial roles in regulating cold stress.

3.6. Validation of RNA-Seq Data by qRT-PCR Analysis

In order to confirm our calculation and analysis results, we selected five important hub genes from the black and blue modules and performed a quantitative reverse-transcription PCR (qRT-PCR) analysis of their expression levels under cold treatment (0 h (ck), 3 h, 6 h, 9 h, 12 h, and 24 h). The results showed that the change trend of transcription level based on qRT-PCR was basically consistent with RNA-Seq (Figure 9). Compared with the control, the expression levels of five genes were significantly up-regulated by cold stress for 12 h or 24 h (Figure 9). This indicates that the RNA-Seq data obtained in this study are reliable and can support the above transcriptome analysis.

4. Discussion

Many previous studies have investigated and described the molecular mechanisms of plant cold tolerance [47,48,49,50,51,52]. The current understanding of the complex mechanisms of cold resistance in holly is limited. This study used transcriptome sequencing to elucidate the cold tolerance mechanisms of holly undergoing cold stress by comparing and analyzing the DEGs of holly under cold treatment (3 h, 6 h, 9 h, 12 h, 24 h, and untreated control conditions).

4.1. Expression Patterns of Plant Hormone Signal Transduction-Related Genes in the Regulation of Cold Stress in I. dabieshanensis

Cold signaling and other signaling pathways that coordinate stress adaptation responds to plant growth and development. Hormonal signals, such as BR, JA, ET, IAA, GA, and ABA, have been reported to regulate cold tolerance [24]. For example, low temperature induces the expression of CBF to inhibit growth, reducing the GA content by stimulating the expression of the GA 2-oxidase gene, thereby enhancing the accumulation of DELLA proteins [53]. The GRETCHEN HAGEN 3 (GH3) family gene OsGH3-2 encodes an enzyme that catalyzes the binding of IAA to amino acids. OsGH3-2 overexpression can enhance cold resistance and reduce ABA levels in rice. F-box protein transport inhibitor response 1 (TIR1) is a receptor for auxin [54], and as putative orthologs of TIR1, the transcript levels of OsAFB2 and OsTIR1 decreased under cold stress [55]. Small auxin-up RNAs (SAURs) are negative regulators of auxin synthesis and transport [56], and nine OsSAUR genes were significantly down-regulated under cold stress in Oryza sativa, and two OsSAUR genes were induced by cold stress [55]. In Arabidopsis thaliana, the cytokinin receptor histidine kinases AHK2-AHK4 and Arabidopsis histidine phosphor-transferase (AHP) may play negative roles under adverse environmental conditions [57]. Most type A response regulators (A-ARRs) are partially redundant negative regulators of cytokinin signaling, induced by cytokinins and cold stress [58]. Ethylene negatively regulates cold signaling through direct transcriptional control of the COR, CBF, and A-ARR genes [59]. Overexpression of ethylene-responsive factor 13 (ERF13) can significantly increase the cold tolerance of birch [60]. Drought stress leads to higher expression levels of SIMKK (salt stress-induced MAPK kinase) genes in Medicago sativa L. and Medicago arborea L. When plants encounter cold stress, they trigger a series of intracellular responses that help them tolerate stress by producing BRs. Activated BZR 1 can induce cold tolerance in plants through a CBF-dependent or CBF-independent pathway [61]. The Arabidopsis Tetracycline hydrochloride (TCH) gene encoding a calmodulin-related protein and xyloglucan endotransglycosylase is known to be up-regulated after cold shock [62]. The basic helix-loop-helix (bHLH) transcription factor MYC2, a master regulator of the core JA signaling machinery, is involved in methyl jasmonate (MeJA)-induced cold tolerance in banana fruit in synergy with MaICE1 [63].
Plant hormone signal transduction was this study’s most highly enriched pathway, excluding cold stress at 3 h (Figure 3). Expression results of differentially expressed genes indicated that seven hormones (ABA, ET, CK, IAA, GA, JA, and BR) were involved in response to cold stress (Figure 5). Under cold stress, the expression of a large number of DELLA genes in holly increased significantly (Figure 5, Table S13). In particular, evm.TU.CHR17.303 was up-regulated 5.9-fold under cold stress at 24 h, and evm.TU.CHR4.1916 was up-regulated 7.2-fold under cold stress at 9 h (Figure 5, Table S13). The expression levels of two ABA receptor genes, pyrabactin resistance (PYR)/PYR1-like (PYL) genes (evm.TU.CHR2.528 and evm.TU.CHR12.725) were up-regulated by cold stress. The evm.TU.CHR2.528 gene was not expressed at 0 h, but its FPKM values were 14.18 and 14.01 at 12 h and 24 h, respectively (Figure 5, Table S13). The above results indicate that holly can also respond to cold stress under low temperatures by regulating the accumulation of DELLA proteins. The expression levels of the two GH3 genes (evm.TU.CHR4.280 and evm.TU.CHR5.340) in holly were significantly changed at the later stage of cold stress (12 h and 24 h), but the expressions were opposite. Two putative TIR1s (evm.TU.CHR1.1393 and evm.TU.CHR4.2549_evm.TU.CHR4.2550) genes were down-regulated by more than 1-fold at 24 h cold treatment (Figure 5, Table S13). Two SAUR genes (evm.TU.CHR5.1764 and evm.TU.CHR19.132) were up-regulated, and one SAUR gene (evm.TU.CHR11.1410) was down-regulated by cold in I. dabieshanensis (Figure 5, Table S13). Under cold stress, I. dabieshanensis adapted to unfavorable environmental factors by regulating auxin homeostasis, signal transduction, and transport. One putative AHP gene (evm.TU.CHR8.1461) was consistent with previous findings, and its expression decreased gradually as cold stress time increased (Figure 5, Table S13). The putative B-ARR genes evm.TU.CHR10.112 and evm.TU.CHR5.1852 were slightly up-regulated by cold stress (Figure 5, Table S13), which was consistent with previous studies. However, the putative A-ARR genes evm.TU.CHR11.1113 and evm.TU.CHR15.84 in I. dabieshanensis were up-regulated and down-regulated (Figure 5, Table S13). Therefore, further exploration is needed to determine the exact role of the cytokinin signaling pathway in holly freezing tolerance. The SIMKK gene (evm.TU.CHR13.25) was up-regulated by cold stress at different times, especially at 24 (up-regulated by about 4 times) (Figure 5, Table S13). The expression levels of two EBF1/2 genes (evm.TU.CHR16.373 and evm.TU.CHR 3.2027) were higher than the control at 12 h and 24 h, especially evm.TU.CHR3.2027, which increased by 1.28 times and 8.43 times, respectively (Figure 5, Table S13). One ERF gene, evm.TU.CHR14.26, was slightly down-regulated by cold stress at 3 h, 6 h, and 9 h, but not significantly. Then evm.TU.CHR14.26 was up-regulated by cold stress at 12 h and 24 h (4.36 times at 24 h). The other two ERF genes were not expressed at 0 h, 3 h, 6 h, and 9 h, but were expressed at 12 h or 24 h (Figure 5, Table S13). In conclusion, I. dabieshanensis’s response to cold stress may involve the participation of ethylene signals, mainly by up-regulating the expression of SIMKK, EBF1/2, and ERF genes. One BZR 1/2 gene (evm.TU.CHR4.816) was induced by cold stress and was up-regulated by 1.23 times after 24 h of cold stress (Figure 5, Table S13). Consistent with previous findings, the expression of the four TCH 4 genes (evm.TU.CHR10.943, evm.TU.CHR3.2328, evm.TU.CHR3.1162, and evm.TU.CHR3.1174) were all up-regulated to varying degrees (Figure 5, Table S13). These results suggest that BR signaling was also involved in mediating the process of low-temperature tolerance in I. dabieshanensis. Under cold stress, the expression levels of four holly MYC2 genes (evm.TU.CHR16.819, evm.TU.CHR5.2033, evm.TU.CHR1.2022, and evm.TU.CHR6.1409) increased to varying degrees at different stress times, especially evm.TU.CHR16.819, evm.TU.CHR5.2033, evm.TU.CHR16.819, and evm.TU.CHR5.2033, which increased by more than two times at 12 h or 24 h (Figure 5, Table S13). This result suggests that JA signaling may positively regulate cold tolerance in I. dabieshanensis through the up-regulation of MYC2 expression. Overall, seven hormonal signaling pathways were involved in holly’s response and tolerance to cold stress. Crosstalk may exist between these hormonal signals, so further research is needed to clarify these relationships.

4.2. Metabolic Pathways during Cold Stress Response in I. dabieshanensis

During cold stress, the oxidation reaction in chloroplasts becomes more intense, resulting in an increase in glycolate content and conversion to glyoxylic acid, which is accompanied by the accumulation of a reactive oxygen species (ROS), hydrogen peroxide (H2O2) [12]. In this study, “Glyoxylate and dicarboxylate metabolism” (ko00630) were highly enriched for DEGs at 9 h of cold stress (Figure 4c) and were also enriched (but not significantly) at other stages of cold stress (Tables S8–S12). This observation indicated that glyoxylic acid metabolism changed the most in I. dabieshanensis under cold stress for 9 h. Proline protects and stabilizes ROS scavenging enzymes and activates alternative detoxification pathways when plants are exposed to various abiotic stresses [64]. Arginine also significantly reduces chilling symptoms and H2O2 accumulation [65]. The DEGs at 12 h and 24 h of I. dabieshanensis cold stress were significantly enriched in “Arginine and proline metabolism” (ko00330) (Figure 4d,e). In particular, 11 genes were enriched at 24 h, and all were up-regulated, of which two genes (evm.TU.CHR2.2709 and evm.TU.CHR5.126) were up-regulated by more than four times (Tables S12 and S16). The red peel is associated with the accumulation of anthocyanins and flavonoids, exhibiting increased resistance to cold and pathogens [66]. In our study, the DEGs of each stage of cold stress were enriched in “Flavonoid biosynthesis” (ko00941) (Tables S8–S12) and were significantly enriched at 6 h, 9 h, and 24 h of cold stress (Figure 4b,c,e). “Anthocyanin biosynthesis” (ko00942) was enriched at 6 h, 9 h, and 24 h (Tables S9, S10 and S12) and was significantly enriched at 6 h and 9 h (Figure 4b,c). This result suggests that the biosynthesis of anthocyanins and flavonoids may play an important role in fruit reddening and low-temperature resistance in I. dabieshanensis. The metabolism of glycolate, arginine, proline, flavonoid, and anthocyanin may regulate the process by which holly responds and tolerates cold stress.

4.3. Peroxisome Pathway during Cold Stress Response in I. dabieshanensis

Usually, low temperatures induce oxidative stress in plants. Mitochondria (in the Ubichinon (UQ)/Ubichinol (UQH2) cycle) activate the plant’s internal antioxidant system consisting of scavenging enzymes and antioxidants to scavenge excess ROS [67,68]. In this study, “Peroxisome” (ko04146) was enriched in each stage of I. dabieshanensis cold stress, but not significantly (Tables S8–S12). The 24 h DEGs under cold stress were most enriched in “Peroxisome” (ko04146) (Tables S8–S12), including six up-regulated genes (two SOD genes, two MPV17 genes, one Epoxide Hydrolase 2 (EPHX2) gene, and one long-chain acyl-coenzyme A synthase (ACSL) gene), and one was down-regulated Nudix Hydrolase 7 (NUDT7) (Figure S5). MPV17 directly scavenges excess ROS [69]. Two MPV17 analogs, evm.TU.CHR18.801 and evm.TU.CHR5.351, were up-regulated by 1.8 and 2.1 times, respectively, at 24 h of cold stress (Table S16). In Brassica rapa L., many cold-induced genes were enriched in the “Peroxisome”, and MPV17 and ACSL were up-regulated [69]. Overall, the peroxisomal pathway is an important cellular process for holly to respond to and defend against cold rapidly.

4.4. Gene Co-Expression Network in the Regulation of Cold Stress in I. dabieshanensis

In this study, 11 co-expression network modules were obtained. Two black and blue modules were significantly associated with cold stress in I. dabieshanensis. Cold-induced film hardening is considered a primary cold-induced event [70]. Plasma membrane fluidity correlates with the proportion of desaturated fatty acids [71]. In this study, “Fatty acid elongation” (ko00062) and “Glycerolipid metabolism” (ko00561), as lipid metabolism-related pathways, were significantly enriched in the blue module. This result suggests lipid metabolism changes in I. dabieshanensis under cold stress, which is consistent with previous studies. Sensing and intracellular transduction of stress signals are critical for the adaptation and survival of organisms to environmental stressors [72]. The GO enrichment analysis results of the black and blue module genes revealed that there were “cellular response to stimulus” (GO:0051716), “cell communication” (GO:0007154), “cell surface receptor signaling pathway” (GO:0007166), and “intracellular signal transduction” (GO:0035556) terms (Tables S17 and S18). This result indicated that stress signaling was also an important component of cold stress tolerance in I. dabieshanensis. Messenger molecules, such as calcium (Ca2+), ROS, and nitric oxide (NO), are known to be involved in plant responses to cold stress [73]. In this study, “response to oxygen-containing compound” (GO:1901700), “response to nitrogen compound” (GO:1901698), and “calcium ion binding” (GO:0005509) were enriched in two cold-related modules (Tables S17 and S18). This result shows that Ca2+, ROS, and nitric oxide (NO) may regulate cold tolerance as signaling molecules in I. dabieshanensis.
Gene regulatory networks mediated by transcription factors play an essential role in plant resistance to various stresses. ERFs serve as a key regulatory center, integrating ET, ABA, JA, and redox signaling in plant responses to many abiotic stressors [41]. HvWRKY1 was reported to play an important role during abiotic stress, possibly regulating cold and drought stress in barley [74]. MYBS3 represses the well-known DREB1/CBF-dependent cold signaling pathway in rice, and this repression appears to operate at the transcriptional level [75]. MpMYBS3, as a key transcription factor in cold signaling, confers cold tolerance in bananas [76]. This study found three genes with higher connectivity in the two cold-related modules (evm.TU.CHR2.244, evm.TU.CHR1.1507, and evm.TU.CHR1.1821) were annotated as “ERF025-like”, “WRKY1b”, and “MYBS3-like”, respectively (Figure 9, Tables S21 and S22). The prediction results of promoter elements of genes co-expressed with these transcription factors show that many genes contain one or more transcription factor action sites in their promoters, while some genes also have some hormone and cold response sites. This result further proves the credibility of the co expression network and the importance of these three key transcription factors (Figure S6). In addition to transcription factors, many receptor-like protein kinases (RLKs) are involved in abiotic stress responses, including ABA response, calcium signaling, and antioxidant defense, implicated in drought, salt, and cold stress [77]. Previous studies have found that the A-type ARR serves as a critical node that integrates ethylene and cytokinin signaling to regulate plant responses to environmental stressors [78]. In addition, it has been reported that cold-induced A-type ARRs play a negative regulatory role in cold stress signaling by inhibiting ABA response, independent of the CBF/DREB pathway [79]. Two highly connected hub genes (evm.TU.CHR2.210 and evm.TU.CHR2.89) in two cold-related modules were functionally described as “a receptor-like serine/threonine- protein kinase”, and “Two-component response regulator-like”, respectively (Figure 9, Tables S21 and S22). In conclusion, five I. dabieshanensis genes (evm.TU.CHR2.244, evm.TU.CHR1.1507, evm.TU. CHR1.1821, evm.TU.CHR2.89, and evm.TU.CHR2.210) have high connectivity and edge weight value, and their analogs have also been shown to play important roles in cold stress. These genes may be the key to holly responding to and tolerating cold stress.

5. Conclusions

Our current study shows that cold stress leads to dramatic changes in gene expression profiles. There were differences in the reaction process of holly under different cold treatment time. Proline metabolism, arginine metabolism, flavonoid biosynthesis, anthocyanin biosynthesis, and various hormone signaling pathways may play an important role in the response of holly to cold stress. The results of WGCNA analysis indicated that some holly genes may be the key genes regulating cold tolerance. These results provide candidate genes and insights for studying the regulatory mechanisms of holly in response to low temperature conditions.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f13122150/s1, Figure S1: GO annotations of cold-induced DEGs; Figure S2: Analysis of transcription factors; Figure S3: Soft ower; Figure S4: Expression of black and blue modules eigengene; Figure S5: Expression of genes associated with the peroxisomal pathway; Figure S6: Analysis of promoter elements of genes co expressed with TF in black and blue module co expression network; Table S1 Primers used for qRT-PCR analysis; Table S2: all.DEG.fpkm; Table S3: Cold-3_vs_Cold-ck.DEG.Go.enrich; Table S4: Cold-6_vs_Cold-ck.DEG.Go.enrich; Table S5: Cold-9_vs_Cold-ck.DEG.Go.enrich; Table S6: Cold-12_vs_Cold-ck.DEG.Go.enrich; Table S7: Cold-24_vs_Cold-ck.DEG.Go.enrich; Table S8: Cold-3_vs_Cold-ck.DEG.Ko.enrich; Table S9: Cold-6_vs_Cold-ck.DEG.Ko.enrich; Table S10: Cold-9_vs_Cold-ck.DEG.Ko.enrich; Table S11: Cold-12_vs_Cold-ck.DEG.Ko.enrich; Table S12: Cold-24_vs_Cold-ck.DEG.Ko.enrich; Table S13: FPKM values and annotations of genes related to Plant hormone signal transduction; Table S14: FPKM values and annotations of genes associated with the MAPK signaling pathway; Table S15: FPKM values and annotations of genes associated with response to cold; Table S16: FPKM values and annotations of other genes associated with cold stress; Table S17: black.Go.enrich; Table S18: blue.Go.enrich; Table S19: black.Ko.enrich; Table S20: blue.Ko.enrich; Table S21: CytoscapeInput-edges-black; Table S22: CytoscapeInput-edges-blue.

Author Contributions

Methodology, H.L. and T.Z.; software, H.L. and X.L.; investigation, H.L., T.Z., X.C., X.L., Y.L., B.Z. and H.C.; resources, H.C.; data curation, H.L., X.W. and T.Z.; writing—original draft preparation, H.L. and H.C.; writing—review and editing, H.C.; visualization, H.L. and X.C.; supervision, H.C. and B.Z.; project administration, H.C. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Jiangsu Agricultural Science and Technology Innovation Fund (grant number: CX(21)3020), the Opening Project of Zhejiang Provincial Key Laboratory of Forest Aromatic Plants-based Healthcare Functions (SLFX202205), the Open Fund of Jiangsu Key Laboratory for the Research and Utilization of Plant Resources (grant number: JSPKLB202209), and the Independent Research Project of Jiangsu Province (grant number: JSPKLB202036).

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Yuan, L.; Wu, H.; Zhang, C.; Wang, Y.; Huang, Q.; Fan, S.; Su, T. The complete plastid genome sequence of Ilex suaveolens (H. Lév.) Loes, the most abundant medicinal holly in Mount Huangshan. Mitochondrial DNA Part B 2021, 6, 468–469. [Google Scholar] [CrossRef] [PubMed]
  2. Loizeau, P.-A.; Barriera, G.; Manen, J.-F.; Broennimann, O. Towards an understanding of the distribution of Ilex L.(Aquifoliaceae) on a world-wide scale. Biol. Skr 2005, 55, 501–520. [Google Scholar]
  3. Yao, X.; Song, Y.; Yang, J.B.; Tan, Y.H.; Corlett, R.T. Phylogeny and biogeography of the hollies (Ilex L., Aquifoliaceae). J. Syst. Evol. 2021, 59, 73–82. [Google Scholar] [CrossRef]
  4. Zhou, T.; Ning, K.; Mo, Z.; Zhang, F.; Zhou, Y.; Chong, X.; Zhang, D.; El-Kassaby, Y.A.; Bian, J.; Chen, H. Complete chloroplast genome of Ilex dabieshanensis: Genome structure, comparative analyses with three traditional Ilex tea species, and its phylogenetic relationships within the family Aquifoliaceae. PLoS ONE 2022, 17, e0268679. [Google Scholar] [CrossRef] [PubMed]
  5. Bracesco, N.; Sanchez, A.; Contreras, V.; Menini, T.; Gugliucci, A. Recent advances on Ilex paraguariensis research: Minireview. J. Ethnopharmacol. 2011, 136, 378–384. [Google Scholar] [CrossRef]
  6. Kungel, P.T.; Correa, V.G.; Corrêa, R.C.; Peralta, R.A.; Soković, M.; Calhelha, R.C.; Bracht, A.; Ferreira, I.C.; Peralta, R.M. Antioxidant and antimicrobial activities of a purified polysaccharide from yerba mate (Ilex paraguariensis). Int. J. Biol. Macromol. 2018, 114, 1161–1167. [Google Scholar] [CrossRef]
  7. Yao, X.; Lu, Z.; Song, Y.; Hu, X.; Corlett, R.T. A chromosome-scale genome assembly for the holly (Ilex polyneura) provides insights into genomic adaptations to elevation in Southwest China. Hortic. Res. 2022, 9, uhab049. [Google Scholar] [CrossRef]
  8. Rakocevic, M.; Medrado, M.; Martim, S.; Assad, E. Sexual dimorphism and seasonal changes of leaf gas exchange in the dioecious tree Ilex paraguariensis grown in two contrasted cultivation types. Ann. Appl. Biol. 2009, 154, 291–301. [Google Scholar] [CrossRef]
  9. Jacques, R.A.; Arruda, E.J.; de Oliveira, L.C.; de Oliveira, A.P.; Dariva, C.; De Oliveira, J.V.; Caramão, E.B. Influence of agronomic variables on the macronutrient and micronutrient contents and thermal behavior of mate tea leaves (Ilex paraguariensis). J. Agric. Food Chem. 2007, 55, 7510–7516. [Google Scholar] [CrossRef]
  10. Theocharis, A.; Clément, C.; Barka, E.A. Physiological and molecular changes in plants grown at low temperatures. Planta 2012, 235, 1091–1105. [Google Scholar] [CrossRef]
  11. Zhou, M.; Chen, H.; Wei, D.; Ma, H.; Lin, J. Arabidopsis CBF3 and DELLAs positively regulate each other in response to low temperature. Sci. Rep. 2017, 7, 1–13. [Google Scholar] [CrossRef]
  12. Ritonga, F.N.; Chen, S. Physiological and molecular mechanism involved in cold stress tolerance in plants. Plants 2020, 9, 560. [Google Scholar] [CrossRef]
  13. Kaplan, F.; Guy, C.L. β-Amylase induction and the protective role of maltose during temperature shock. Plant Physiol. 2004, 135, 1674–1684. [Google Scholar] [CrossRef]
  14. Kaplan, F.; Kopka, J.; Sung, D.Y.; Zhao, W.; Popp, M.; Porat, R.; Guy, C.L. Transcript and metabolite profiling during cold acclimation of Arabidopsis reveals an intricate relationship of cold-regulated gene expression with modifications in metabolite content. Plant J. 2007, 50, 967–981. [Google Scholar] [CrossRef]
  15. Lee, B.-h.; Lee, H.; Xiong, L.; Zhu, J.-K. A mitochondrial complex I defect impairs cold-regulated nuclear gene expression. Plant Cell 2002, 14, 1235–1251. [Google Scholar] [CrossRef]
  16. Dong, C.-H.; Zolman, B.K.; Bartel, B.; Lee, B.-h.; Stevenson, B.; Agarwal, M.; Zhu, J.-K. Disruption of Arabidopsis CHY1 reveals an important role of metabolic status in plant cold stress signaling. Mol. Plant 2009, 2, 59–72. [Google Scholar] [CrossRef]
  17. Yin, J.; Yi, H.; Chen, X.; Wang, J. Post-translational modifications of proteins have versatile roles in regulating plant immune responses. Int. J. Mol. Sci. 2019, 20, 2807. [Google Scholar] [CrossRef]
  18. Jofuku, K.D.; Den Boer, B.; Van Montagu, M.; Okamuro, J.K. Control of Arabidopsis flower and seed development by the homeotic gene APETALA2. Plant Cell 1994, 6, 1211–1225. [Google Scholar] [CrossRef]
  19. Wei, X.; Liu, S.; Sun, C.; Xie, G.; Wang, L. Convergence and divergence: Signal perception and transduction mechanisms of cold stress in Arabidopsis and Rice. Plants 2021, 10, 1864. [Google Scholar] [CrossRef]
  20. Lin, L.; Wu, J.; Jiang, M.; Wang, Y. Plant mitogen-activated protein kinase cascades in environmental stresses. Int. J. Mol. Sci. 2021, 22, 1543. [Google Scholar] [CrossRef]
  21. Yang, T.; Chaudhuri, S.; Yang, L.; Du, L.; Poovaiah, B. A calcium/calmodulin-regulated member of the receptor-like kinase family confers cold tolerance in plants. J. Biol. Chem. 2010, 285, 7119–7126. [Google Scholar] [CrossRef] [PubMed]
  22. Shi, Y.; Ding, Y.; Yang, S. Molecular regulation of CBF signaling in cold acclimation. Trends Plant Sci. 2018, 23, 623–637. [Google Scholar] [CrossRef] [PubMed]
  23. Alves, M.S.; Fontes, E.P.; Fietto, L.G. Early responsive to dehydration 15, a new transcription factor that integrates stress signaling pathways. Plant Signal. Behav. 2011, 6, 1993–1996. [Google Scholar] [CrossRef] [PubMed]
  24. Hwarari, D.; Guan, Y.; Ahmad, B.; Movahedi, A.; Min, T.; Hao, Z.; Lu, Y.; Chen, J.; Yang, L. ICE-CBF-COR signaling cascade and its regulation in plants responding to cold stress. Int. J. Mol. Sci. 2022, 23, 1549. [Google Scholar] [CrossRef] [PubMed]
  25. Zhang, F.; Jiang, Y.; Bai, L.; Zhang, L.; Chen, L.; Li, H.; Yin, Y.; Yan, W.; Yi, Y.; Guo, Z. The ICE-CBF-COR pathway in cold acclimation and AFPs in plants. Middle East J. Sci. Res. 2011, 8, 493–498. [Google Scholar]
  26. Borba, A.R.; Serra, T.S.; Górska, A.; Gouveia, P.; Cordeiro, A.M.; Reyna-Llorens, I.; Kneřová, J.; Barros, P.M.; Abreu, I.A.; Oliveira, M.M. Synergistic binding of bHLH transcription factors to the promoter of the maize NADP-ME gene used in C4 photosynthesis is based on an ancient code found in the ancestral C3 state. Mol. Biol. Evol. 2018, 35, 1690–1705. [Google Scholar] [CrossRef]
  27. Mickelbart, M.V.; Hasegawa, P.M.; Bailey-Serres, J. Genetic mechanisms of abiotic stress tolerance that translate to crop yield stability. Nat. Rev. Genet. 2015, 16, 237–251. [Google Scholar] [CrossRef]
  28. Shu, Y.; Li, W.; Zhao, J.; Zhang, S.; Xu, H.; Liu, Y.; Guo, C. Transcriptome sequencing analysis of alfalfa reveals CBF genes potentially playing important roles in response to freezing stress. Genet. Mol. Biol. 2017, 40, 824–833. [Google Scholar] [CrossRef]
  29. Thomashow, M.F. Plant cold acclimation: Freezing tolerance genes and regulatory mechanisms. Annu. Rev. Plant Biol. 1999, 50, 571–599. [Google Scholar] [CrossRef]
  30. Guo, J.; Ren, Y.; Tang, Z.; Shi, W.; Zhou, M. Characterization and expression profiling of the ICE-CBF-COR genes in wheat. PeerJ 2019, 7, e8190. [Google Scholar] [CrossRef]
  31. Verma, R.K.; Kumar, V.V.S.; Yadav, S.K.; Kumar, T.S.; Rao, M.V.; Chinnusamy, V. Overexpression of Arabidopsis ICE1 enhances yield and multiple abiotic stress tolerance in indica rice. Plant Signal. Behav. 2020, 15, 1814547. [Google Scholar] [CrossRef]
  32. Li, X.; Liu, C.; Zhao, Z.; Ma, D.; Zhang, J.; Yang, Y.; Liu, Y.; Liu, H. COR27 and COR28 are novel regulators of the COP1–HY5 regulatory hub and photomorphogenesis in Arabidopsis. Plant Cell 2020, 32, 3139–3154. [Google Scholar] [CrossRef]
  33. Wang, S.; Yang, C.; Zhao, X.; Chen, S.; Qu, G.-Z. Complete chloroplast genome sequence of Betula platyphylla: Gene organization, RNA editing, and comparative and phylogenetic analyses. BMC Genom. 2018, 19, 950. [Google Scholar] [CrossRef]
  34. Liu, Y.; Zhou, J. MAPping kinase regulation of ICE1 in freezing tolerance. Trends Plant Sci. 2018, 23, 91–93. [Google Scholar] [CrossRef]
  35. Matsukura, S.; Mizoi, J.; Yoshida, T.; Todaka, D.; Ito, Y.; Maruyama, K.; Shinozaki, K.; Yamaguchi-Shinozaki, K. Comprehensive analysis of rice DREB2-type genes that encode transcription factors involved in the expression of abiotic stress-responsive genes. Mol. Genet. Genom. 2010, 283, 185–196. [Google Scholar] [CrossRef]
  36. Teige, M.; Scheikl, E.; Eulgem, T.; Dóczi, R.; Ichimura, K.; Shinozaki, K.; Dangl, J.L.; Hirt, H. The MKK2 pathway mediates cold and salt stress signaling in Arabidopsis. Mol. Cell 2004, 15, 141–152. [Google Scholar] [CrossRef]
  37. Shinozaki, K.; Yamaguchi-Shinozaki, K. Molecular responses to drought and cold stress. Curr. Opin. Biotechnol. 1996, 7, 161–167. [Google Scholar] [CrossRef]
  38. Rihan, H.Z.; Al-Issawi, M.; Fuller, M.P. Upregulation of CBF/DREB1 and cold tolerance in artificial seeds of cauliflower (Brassica oleracea var. botrytis). Sci. Hortic. 2017, 225, 299–309. [Google Scholar] [CrossRef]
  39. Wang, D.-Z.; Jin, Y.-N.; Ding, X.-H.; Wang, W.-J.; Zhai, S.-S.; Bai, L.-P.; Guo, Z.-F. Gene regulation and signal transduction in the ICE–CBF–COR signaling pathway during cold stress in plants. Biochemistry 2017, 82, 1103–1117. [Google Scholar] [CrossRef]
  40. Zhou, B.; Lin, J.; Peng, W.; Peng, D.; Zhuo, Y.; Zhu, D.; Huang, X.; Tang, D.; Guo, M.; He, R. Dwarfism in Brassica napus L. induced by the over-expression of a gibberellin 2-oxidase gene from Arabidopsis thaliana. Mol. Breed. 2012, 29, 115–127. [Google Scholar] [CrossRef]
  41. Müller, M.; Munné-Bosch, S. Ethylene response factors: A key regulatory hub in hormone and stress signaling. Plant Physiol. 2015, 169, 32–41. [Google Scholar] [CrossRef] [PubMed]
  42. Barrero-Gil, J.; Salinas, J. CBFs at the crossroads of plant hormone signaling in cold stress response. Mol. Plant 2017, 10, 542–544. [Google Scholar] [CrossRef] [PubMed]
  43. Tolosa, L.N.; Zhang, Z. The role of major transcription factors in Solanaceous food crops under different stress conditions: Current and future perspectives. Plants 2020, 9, 56. [Google Scholar] [CrossRef] [PubMed]
  44. Liu, Y.; Jiang, Y.; Lan, J.; Zou, Y.; Gao, J. Comparative transcriptomic analysis of the response to cold acclimation in Eucalyptus dunnii. PLoS ONE 2014, 9, e113091. [Google Scholar] [CrossRef]
  45. Langfelder, P.; Horvath, S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinform. 2008, 9, 559. [Google Scholar] [CrossRef]
  46. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  47. Pagter, M.; Alpers, J.; Erban, A.; Kopka, J.; Zuther, E.; Hincha, D.K. Rapid transcriptional and metabolic regulation of the deacclimation process in cold acclimated Arabidopsis thaliana. BMC Genom. 2017, 18, 731. [Google Scholar] [CrossRef]
  48. Song, Y.; Liu, L.; Wei, Y.; Li, G.; Yue, X.; An, L. Metabolite profiling of adh1 mutant response to cold stress in Arabidopsis. Front. Plant Sci. 2017, 7, 2072. [Google Scholar] [CrossRef]
  49. Roux, B.; Bolot, S.; Guy, E.; Denancé, N.; Lautier, M.; Jardinaud, M.-F.; Fischer-Le Saux, M.; Portier, P.; Jacques, M.-A.; Gagnevin, L. Genomics and transcriptomics of Xanthomonas campestris species challenge the concept of core type III effectome. BMC Genom. 2015, 16, 975. [Google Scholar] [CrossRef]
  50. Kim, S.; Ahn, S.Y.; Yun, H.K. Transcriptome analysis of grapevine shoots exposed to chilling temperature for four weeks. Hortic. Environ. Biotechnol. 2016, 57, 161–172. [Google Scholar] [CrossRef]
  51. Zeng, X.; Xu, Y.; Jiang, J.; Zhang, F.; Ma, L.; Wu, D.; Wang, Y.; Sun, W. Identification of cold stress responsive microRNAs in two winter turnip rape (Brassica rapa L.) by high throughput sequencing. BMC Plant Biol. 2018, 18, 52. [Google Scholar] [CrossRef]
  52. Wang, J.; Yao, L.; Li, B.; Meng, Y.; Ma, X.; Lai, Y.; Si, E.; Ren, P.; Yang, K.; Shang, X. Comparative proteomic analysis of cultured suspension cells of the halophyte Halogeton glomeratus by iTRAQ provides insights into response mechanisms to salt stress. Front. Plant Sci. 2016, 7, 110. [Google Scholar] [CrossRef]
  53. Achard, P.; Gong, F.; Cheminant, S.; Alioua, M.; Hedden, P.; Genschik, P. The cold-inducible CBF1 factor–dependent signaling pathway modulates the accumulation of the growth-repressing DELLA proteins via its effect on gibberellin metabolism. Plant Cell 2008, 20, 2117–2129. [Google Scholar] [CrossRef]
  54. Tan, X.; Calderon-Villalobos, L.I.A.; Sharon, M.; Zheng, C.; Robinson, C.V.; Estelle, M.; Zheng, N. Mechanism of auxin perception by the TIR1 ubiquitin ligase. Nature 2007, 446, 640–645. [Google Scholar] [CrossRef]
  55. Du, H.; Liu, H.; Xiong, L. Endogenous auxin and jasmonic acid levels are differentially modulated by abiotic stresses in rice. Front. Plant Sci. 2013, 4, 397. [Google Scholar] [CrossRef]
  56. Kant, S.; Bi, Y.-M.; Zhu, T.; Rothstein, S.J. SAUR39, a small auxin-up RNA gene, acts as a negative regulator of auxin synthesis and transport in rice. Plant Physiol. 2009, 151, 691–701. [Google Scholar] [CrossRef]
  57. Shi, Y.; Yang, S. ABA regulation of the cold stress response in plants. In Abscisic Acid: Metabolism, Transport and Signaling; Springer: Berlin/Heidelberg, Germany, 2014; pp. 337–363. [Google Scholar]
  58. To, J.P.; Haberer, G.; Ferreira, F.J.; Deruere, J.; Mason, M.G.; Schaller, G.E.; Alonso, J.M.; Ecker, J.R.; Kieber, J.J. Type-A Arabidopsis response regulators are partially redundant negative regulators of cytokinin signaling. Plant Cell 2004, 16, 658–671. [Google Scholar] [CrossRef]
  59. Shi, B.; Ni, L.; Zhang, A.; Cao, J.; Zhang, H.; Qin, T.; Tan, M.; Zhang, J.; Jiang, M. OsDMI3 is a novel component of abscisic acid signaling in the induction of antioxidant defense in leaves of rice. Mol. Plant 2012, 5, 1359–1374. [Google Scholar] [CrossRef]
  60. Lv, K.; Li, J.; Zhao, K.; Chen, S.; Nie, J.; Zhang, W.; Liu, G.; Wei, H. Overexpression of an AP2/ERF family gene, BpERF13, in birch enhances cold tolerance through upregulating CBF genes and mitigating reactive oxygen species. Plant Sci. 2020, 292, 110375. [Google Scholar] [CrossRef]
  61. Ali, M.; Ihsan, M.; Zafar, H.; Rauf, Q.; Akhtar, M. Brassinosteroid biosynthesis, stress resistance in plants, and application of brassinosteroids in plant biotechnology. J. Biol. Regul. Homeost. Agents 2018, 32, 1457–1459. [Google Scholar]
  62. Purugganan, M.M.; Braam, J.; Fry, S.C. The Arabidopsis TCH4 xyloglucan endotransglycosylase (substrate specificity, pH optimum, and cold tolerance). Plant Physiol. 1997, 115, 181–190. [Google Scholar] [CrossRef] [PubMed]
  63. Wang, J.; Song, L.; Gong, X.; Xu, J.; Li, M. Functions of jasmonic acid in plant regulation and response to abiotic stress. Int. J. Mol. Sci. 2020, 21, 1446. [Google Scholar] [CrossRef] [PubMed]
  64. Khedr, A.H.A.; Abbas, M.A.; Wahid, A.A.A.; Quick, W.P.; Abogadallah, G.M. Proline induces the expression of salt-stress-responsive proteins and may improve the adaptation of Pancratium maritimum L. to salt-stress. J. Exp. Bot. 2003, 54, 2553–2562. [Google Scholar] [CrossRef] [PubMed]
  65. Babalar, M.; Pirzad, F.; Sarcheshmeh, M.A.A.; Talaei, A.; Lessani, H. Arginine treatment attenuates chilling injury of pomegranate fruit during cold storage by enhancing antioxidant system activity. Postharvest Biol. Technol. 2018, 137, 31–37. [Google Scholar] [CrossRef]
  66. Sivankalyani, V.; Feygenberg, O.; Diskin, S.; Wright, B.; Alkan, N. Increased anthocyanin and flavonoids in mango fruit peel are associated with cold and pathogen resistance. Postharvest Biol. Technol. 2016, 111, 132–139. [Google Scholar] [CrossRef]
  67. Guo, Z.; Ou, W.; Lu, S.-y.; Zhong, Q. Differential responses of antioxidative system to chilling and drought in four rice cultivars differing in sensitivity. Plant Physiol. Biochem. 2006, 44, 828–836. [Google Scholar] [CrossRef]
  68. Goulas, E.; Schubert, M.; Kieselbach, T.; Kleczkowski, L.A.; Gardeström, P.; Schröder, W.; Hurry, V. The chloroplast lumen and stromal proteomes of Arabidopsis thaliana show differential sensitivity to short- and long-term exposure to low temperature. Plant J. 2006, 47, 720–734. [Google Scholar] [CrossRef]
  69. Ma, L.; Coulter, J.A.; Liu, L.; Zhao, Y.; Chang, Y.; Pu, Y.; Zeng, X.; Xu, Y.; Wu, J.; Fang, Y. Transcriptome analysis reveals key cold-stress-responsive genes in winter rapeseed (Brassica rapa L.). Int. J. Mol. Sci. 2019, 20, 1071. [Google Scholar] [CrossRef]
  70. Örvar, B.L.; Sangwan, V.; Omann, F.; Dhindsa, R.S. Early steps in cold sensing by plant cells: The role of actin cytoskeleton and membrane fluidity. Plant J. 2000, 23, 785–794. [Google Scholar] [CrossRef]
  71. Martinière, A.; Shvedunova, M.; Thomson, A.J.; Evans, N.H.; Penfield, S.; Runions, J.; McWatters, H.G. Homeostasis of plasma membrane viscosity in fluctuating temperatures. New Phytol. 2011, 192, 328–337. [Google Scholar] [CrossRef]
  72. Long, Y.; Li, L.; Li, Q.; He, X.; Cui, Z. Transcriptomic characterization of temperature stress responses in larval zebrafish. PLoS ONE 2012, 7, e37209. [Google Scholar] [CrossRef]
  73. Ding, Y.; Shi, Y.; Yang, S. Advances and challenges in uncovering cold tolerance regulatory mechanisms in plants. New Phytol. 2019, 222, 1690–1704. [Google Scholar] [CrossRef]
  74. Mangelsen, E.; Kilian, J.; Berendzen, K.W.; Kolukisaoglu, Ü.H.; Harter, K.; Jansson, C.; Wanke, D. Phylogenetic and comparative gene expression analysis of barley (Hordeum vulgare) WRKY transcription factor family reveals putatively retained functions between monocots and dicots. BMC Genom. 2008, 9, 194. [Google Scholar] [CrossRef]
  75. Su, C.-F.; Wang, Y.-C.; Hsieh, T.-H.; Lu, C.-A.; Tseng, T.-H.; Yu, S.-M. A novel MYBS3-dependent pathway confers cold tolerance in rice. Plant Physiol. 2010, 153, 145–158. [Google Scholar] [CrossRef]
  76. Dou, T.-X.; Hu, C.-H.; Sun, X.-X.; Shao, X.-H.; Wu, J.-H.; Ding, L.-J.; Gao, J.; He, W.-D.; Biswas, M.-K.; Yang, Q.-S. MpMYBS3 as a crucial transcription factor of cold signaling confers the cold tolerance of banana. Plant Cell Tissue Organ Cult. (PCTOC) 2016, 125, 93–106. [Google Scholar] [CrossRef]
  77. Ye, Y.; Ding, Y.; Jiang, Q.; Wang, F.; Sun, J.; Zhu, C. The role of receptor-like protein kinases (RLKs) in abiotic stress response in plants. Plant Cell Rep. 2017, 36, 235–242. [Google Scholar] [CrossRef]
  78. Shi, Y.; Tian, S.; Hou, L.; Huang, X.; Zhang, X.; Guo, H.; Yang, S. Ethylene signaling negatively regulates freezing tolerance by repressing expression of CBF and type-A ARR genes in Arabidopsis. Plant Cell 2012, 24, 2578–2595. [Google Scholar] [CrossRef]
  79. Jeon, J.; Kim, N.Y.; Kim, S.; Kang, N.Y.; Novák, O.; Ku, S.-J.; Cho, C.; Lee, D.J.; Lee, E.-J.; Strnad, M. A subset of cytokinin two-component signaling system plays a role in cold temperature stress response in Arabidopsis. J. Biol. Chem. 2010, 285, 23371–23386. [Google Scholar] [CrossRef]
Figure 1. I. dabieshanensis in the snow.
Figure 1. I. dabieshanensis in the snow.
Forests 13 02150 g001
Figure 2. Volcano-plots of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; (e) 24 h vs. CK; (f) 6 h vs. 3 h; (g) 9 h vs. 3 h; (h) 12 h vs. 3 h; (i) 24 h vs. 3 h; (j) 9 h vs. 6 h; (k) 12 h vs. 6 h; (l) 24 h vs. 6 h; (m) 12 h vs. 9 h; (n) 24 h vs. 9 h; (o) 24 h vs. 12 h. The abscissa is the fold change value of the gene or transcript expression difference between two samples (sample 1/sample 2), and the ordinate is the statistical test value of the difference in gene or transcript expression change; that is, the higher the p-value, the more significant the difference in expression. The values of the abscissa and ordinate are all logarithmic. Each dot in the figure represents a specific gene or transcript; the red dots represent significantly up-regulated genes, the blue dots represent significantly down-regulated genes, and the black dots represent non-significant differential genes. The dots on the left represent differentially down-regulated genes, and the dots on the right are genes whose expression is differentially up-regulated.
Figure 2. Volcano-plots of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; (e) 24 h vs. CK; (f) 6 h vs. 3 h; (g) 9 h vs. 3 h; (h) 12 h vs. 3 h; (i) 24 h vs. 3 h; (j) 9 h vs. 6 h; (k) 12 h vs. 6 h; (l) 24 h vs. 6 h; (m) 12 h vs. 9 h; (n) 24 h vs. 9 h; (o) 24 h vs. 12 h. The abscissa is the fold change value of the gene or transcript expression difference between two samples (sample 1/sample 2), and the ordinate is the statistical test value of the difference in gene or transcript expression change; that is, the higher the p-value, the more significant the difference in expression. The values of the abscissa and ordinate are all logarithmic. Each dot in the figure represents a specific gene or transcript; the red dots represent significantly up-regulated genes, the blue dots represent significantly down-regulated genes, and the black dots represent non-significant differential genes. The dots on the left represent differentially down-regulated genes, and the dots on the right are genes whose expression is differentially up-regulated.
Forests 13 02150 g002
Figure 3. GO enrichment map of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; and (e) 24 h vs. CK. The horizontal axis represents the enrichment factor, and the ratio of differential genes enriched to a GO term to the number of background genes obtained by sequencing is presented. The vertical axis represents the enriched function of the GO term: the larger the circle, the more enriched. The number of differential genes for this function is relatively greater. The color spectrum from blue to red represents uncorrected p-values.
Figure 3. GO enrichment map of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; and (e) 24 h vs. CK. The horizontal axis represents the enrichment factor, and the ratio of differential genes enriched to a GO term to the number of background genes obtained by sequencing is presented. The vertical axis represents the enriched function of the GO term: the larger the circle, the more enriched. The number of differential genes for this function is relatively greater. The color spectrum from blue to red represents uncorrected p-values.
Forests 13 02150 g003
Figure 4. KEGG enrichment map of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; (e) 24 h vs. CK. The horizontal axis represents the enrichment factor, that is, the ratio of the number of differential genes enriched to a certain KEGG pathway to the number of background genes obtained by sequencing; the vertical axis represents the enriched function of the KEGG pathway: the larger the circle, the more enriched. The number of differential genes for this function is relatively more. The color spectrum from blue to red represents uncorrected p-values.
Figure 4. KEGG enrichment map of cold-induced DEGs. (a) 3 h vs. CK (0 h); (b) 6 h vs. CK; (c) 9 h vs. CK; (d) 12 h vs. CK; (e) 24 h vs. CK. The horizontal axis represents the enrichment factor, that is, the ratio of the number of differential genes enriched to a certain KEGG pathway to the number of background genes obtained by sequencing; the vertical axis represents the enriched function of the KEGG pathway: the larger the circle, the more enriched. The number of differential genes for this function is relatively more. The color spectrum from blue to red represents uncorrected p-values.
Forests 13 02150 g004
Figure 5. Expression of genes related to response to cold stress in I. dabieshanensis. The original figure is from Hwarari’s figure, and it was modified according to our results. The squares in each column of the heatmap represent cold treatment at 0 h (CK), 3 h, 6 h, 9 h, 12 h, and 24 h, respectively, from left to right.
Figure 5. Expression of genes related to response to cold stress in I. dabieshanensis. The original figure is from Hwarari’s figure, and it was modified according to our results. The squares in each column of the heatmap represent cold treatment at 0 h (CK), 3 h, 6 h, 9 h, 12 h, and 24 h, respectively, from left to right.
Forests 13 02150 g005
Figure 6. WGCNA analysis of DEGs. (a) The clustering dendrogram of DEGs identifying the WGCNA modules. (b) The correlation of the identified modules with the different cold treatments.
Figure 6. WGCNA analysis of DEGs. (a) The clustering dendrogram of DEGs identifying the WGCNA modules. (b) The correlation of the identified modules with the different cold treatments.
Forests 13 02150 g006
Figure 7. Function enrichment analysis of black and blue module genes. (a) GO enrichment map of black module genes. (b) KEGG enrichment map of black module genes. (c) GO enrichment map of blue module genes. (d) KEGG enrichment map of blue module genes.
Figure 7. Function enrichment analysis of black and blue module genes. (a) GO enrichment map of black module genes. (b) KEGG enrichment map of black module genes. (c) GO enrichment map of blue module genes. (d) KEGG enrichment map of blue module genes.
Forests 13 02150 g007
Figure 8. Co-expression regulatory network analysis of the black and blue module. Co-expression regulatory network analysis of the black and blue module. (a) Co-expression network constructed by 30 hub genes in black modules; (b) Co-expression network constructed by 30 hub genes in blue modules. Each circle represents a gene. The red circles represent transcription factors. The larger the circle, the higher the number of connections. The darker the black line connecting two nodes, the higher the edge weight value.
Figure 8. Co-expression regulatory network analysis of the black and blue module. Co-expression regulatory network analysis of the black and blue module. (a) Co-expression network constructed by 30 hub genes in black modules; (b) Co-expression network constructed by 30 hub genes in blue modules. Each circle represents a gene. The red circles represent transcription factors. The larger the circle, the higher the number of connections. The darker the black line connecting two nodes, the higher the edge weight value.
Forests 13 02150 g008
Figure 9. The expression analysis of five genes from the I. dabieshanensis core cold stress responsive genes. The relative expression of five genes was obtained by qRT-PCR based on 2−∆∆Ct method [46], while FPKM values of five genes were obtained by RNA-Seq. Data are from three biological replicates and three technical replicates. Results are the means±standard error of the three.
Figure 9. The expression analysis of five genes from the I. dabieshanensis core cold stress responsive genes. The relative expression of five genes was obtained by qRT-PCR based on 2−∆∆Ct method [46], while FPKM values of five genes were obtained by RNA-Seq. Data are from three biological replicates and three technical replicates. Results are the means±standard error of the three.
Forests 13 02150 g009
Table 1. GO enrichment analysis statistics of DEGs (p-value ≤ 0.05).
Table 1. GO enrichment analysis statistics of DEGs (p-value ≤ 0.05).
SampleDEG-TypeGO Term_
Enrich_Num.
BPCCMF
Cold-3_vs_Cold-ckDEG-all1601241719
Cold-3_vs_Cold-ckDEG-up1621271817
Cold-3_vs_Cold-ckDEG-down373250
Cold-6_vs_Cold-ckDEG-all2401841739
Cold-6_vs_Cold-ckDEG-up2151651832
Cold-6_vs_Cold-ckDEG-down10790116
Cold-9_vs_Cold-ckDEG-all2391971032
Cold-9_vs_Cold-ckDEG-up2251731339
Cold-9_vs_Cold-ckDEG-down706514
Cold-12_vs_Cold-ckDEG-all329271652
Cold-12_vs_Cold-ckDEG-up3392741055
Cold-12_vs_Cold-ckDEG-down605604
Cold-24_vs_Cold-ckDEG-all3102241076
Cold-24_vs_Cold-ckDEG-up3953031973
Cold-24_vs_Cold-ckDEG-down7944530
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Li, H.; Zhou, T.; Chong, X.; Lu, X.; Li, Y.; Zheng, B.; Wang, X.; Chen, H. Transcriptome and Expression Analysis of Genes Related to Regulatory Mechanisms in Holly (Ilex dabieshanensis) under Cold Stress. Forests 2022, 13, 2150. https://doi.org/10.3390/f13122150

AMA Style

Li H, Zhou T, Chong X, Lu X, Li Y, Zheng B, Wang X, Chen H. Transcriptome and Expression Analysis of Genes Related to Regulatory Mechanisms in Holly (Ilex dabieshanensis) under Cold Stress. Forests. 2022; 13(12):2150. https://doi.org/10.3390/f13122150

Chicago/Turabian Style

Li, Huihui, Ting Zhou, Xinran Chong, Xiaoqing Lu, Yunlong Li, Bingsong Zheng, Xiaolong Wang, and Hong Chen. 2022. "Transcriptome and Expression Analysis of Genes Related to Regulatory Mechanisms in Holly (Ilex dabieshanensis) under Cold Stress" Forests 13, no. 12: 2150. https://doi.org/10.3390/f13122150

APA Style

Li, H., Zhou, T., Chong, X., Lu, X., Li, Y., Zheng, B., Wang, X., & Chen, H. (2022). Transcriptome and Expression Analysis of Genes Related to Regulatory Mechanisms in Holly (Ilex dabieshanensis) under Cold Stress. Forests, 13(12), 2150. https://doi.org/10.3390/f13122150

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop