De novo transcriptome sequencing and gene expression profiling of Pinus sibirica under different cold stresses

Abstract

Background: Pinus sibirica is an evergreen conifer tree species with strong cold stress. However, the transcriptional regulation patterns in response to cold stress have not been reported for P. Sibirica. To gain deeper insights into its regulation process of cold tolerance, transcriptome profiling analyses and 12 physiological indices measurement were performed under cold stress (-20 ℃) over time.

Results: More than 54.1 million clean reads were produced, which were assembled into 97,376 unigenes. Among them, 56,994 unigenes had homology with known genes, 36,836 were assigned to 51 GO (gene ontology) categories and 46,972 were assigned to 24 COG (clusters of orthologous group) categories. P. sibirica showed the highest similarity with sequences from Picea sitchensis. In total, 871, 1397 and 872 DEGs (differentially expressed genes) were identified upon exposure to cold for 6 h, 24 h and 48 h at -20 ℃, respectively. Nine physiological indices increased significantly (P<0.05) under cold stress, including membrane permeability, relative conductivity, reactive oxygen species, malonaldehyde, peroxidase activity, catalase activity, soluble sugar, soluble protein and proline content. With extension of the cold stress time, 9 physiological indices generally showed a trend toward first an increase and then a decrease. The net photosynthetic rate, stomatal conductance and transpiration rate in P. sibirica dropped sharply (P<0.05) in response to cold stress, and they were also decreased significantly (P<0.05) with extension of the stress time at -20 ℃.

Conclusions: There were two cold signal transduction pathways in P. sibirica, the Ca²⁺ and ABA (abscisic acid) pathways. The AP2 (ethylene-responsive transcription factor) family and some other transcription factors played an important role in transcriptional regulation. P. sibirica underwent antioxidant and osmotic regulation with changes in the expression of genes related to cold resistance. Photosynthesis was inhibited, and more DEGs associated with photosynthesis were downregulated under cold stress. The DEGs identified in cold signal sensing and transduction and transcriptional, antioxidant and osmotic regulation can provide genetic resources for the improvement of cold-tolerant characters in other conifer tree species and facilitate understanding of molecular control mechanism related to cold responses.

Background

Pinus sibirica is an evergreen conifer tree species belonging to the Pinaceae, Pinus, mainly distributed in Siberia, Russia (49°40’ ~ 127°20’ E, 46°40' ~ 68°30’ N) [1]. In China, it is only distributed in the northwest of the Altai mountain, Xinjiang. Due to the excellent wood and nutritious seeds, P. sibirica is widely used in furniture, construction, dried fruit and oil, and it is extremely frigostable [2]. In some localities, the winter minimum temperature can fall to below − 65 °C, and the frost-free period is sometimes less than a month. P. sibirica sometimes encounters late frost, and buds suffer freezing injury, but most can repair themselves [2]. To date, studies on the molecular mechanism of low temperature stress in conifer species are scarce. To understand the molecular regulatory process of P. sibirica in resisting low temperature limitations, exploring the gene expression patterns of P. sibirica under cold stress is imperative. It could provide an ideal model and guidelines for research on the cold tolerance molecular regulation of other cold-sensitive conifer species, and ultimately increase the composition of species in high latitude cold regions and improve ecological benefits.

Cold stress is a major abiotic factor that adversely affects plant growth, development, productivity, quality and geographical distribution [3]. Over time, plants have evolved a series of regulatory mechanisms that trigger complex gene expression processes that have subsequently led to physiological and biochemical modifications to resist cold stress [4]. Physiological changes mainly include a sharp increase in Ca²⁺ and ABA (abscisic acid) concentrations [5], changes in membrane lipid composition and permeability [6], increases in antioxidant enzymes and antioxidants [7], alterations in plant stomatas and photosynthesis [8], and the accumulation of osmotic regulatory materials [6].

The molecular mechanism of plant adaption to the low temperature environment is very complex, in which gene expression changes play a critical role [9]. Numerous studies have revealed the gene expression changes in response to cold stress using genome-wide transcriptome analysis [10, 11, 12, 13] and these changes to highlight the importance of transcriptional regulation. For example, in research on Pinus koraensis transcriptome profiling, a total of 96 differentially expressed genes (41 upregulated, 55 downregulated) encoding 8 transcription factors were shown to be involved in transcriptional regulation [14]. Among these, the C-repeat binding factors (CBF)/dehydration-responsive element–binding factors 1 (DREB1) transcription factor (TF) provide a key molecular switching role in cold regulation [15]. At present, there are manly two kinds of cold-induced regulation mechanisms in plants, one is ICE1-CBFs and the other is not dependent on CBF [16, 17]. CBF can active downstream of cold response (COR) genes by binding to cis-acting elements in their promoters [18]. Most COR genes, such as blue copper protein, heat shock protein, osmotic regulators and antioxidant enzymes protect plant cells from cold damage [4]. In rice seedlings, OsSODB (superoxide dismutase), OsTrx23 (thioredoxin), and OsLti6 (encoding membrane proteins) accumulate in response to cold [19]. Although many such studies have been conducted, there is currently no research on the gene expression profile in response to cold stress in P. sibirica.

To elucidate physiological and transcriptomic adaptive mechanisms in P. sibirica, 5-year-old seedlings showing healthy growth at room temperature (20 ℃) were suddenly exposed to -20 ℃ for different time courses (6 h, 24 h and 48 h). Twelve physiological indices were measured, including membrane permeability, relative conductivity (REC), reactive oxygen species (ROS), malonaldehyde (MDA), peroxidase (POD), catalase (CAT), soluble sugar, soluble protein, proline, net photosynthetic rate, stomatal conductance and transpiration rate. Simultaneously, RNA-seq was performed. This study provides insights into the molecular mechanisms of resistance to cold stress in P. sibirica. The identification of cold response genes such as those encoding signal transduction proteins, TFs, relevant antioxidant enzymes and osmotic substances would provide a valuable genetic resource in future conifer breeding.

Results

RNA sequencing and de novo assembly

In this study, 4 treatments were applied to P. sibirica, each with 3 biological duplicates, resulting in the construction of a total of 12 cDNA libraries using corresponding RNA samples. The cDNA libraries were subjected to paired-end (PE) sequencing using the Illumina HiSeq2000 platform. After filtering out low-quality reads, 541,809,548 obtained clean reads were integrated and assembled into 97,376 unigenes with a mean length of 771 bp. The detailed information is listed in Additional file 1. For the size distribution of unigenes in P. sibirica, the majority of unigenes were between 300 bp and 600 bp in length, only 0.8% of the unigenes were between 2700 bp and 3000 bp, and 2% of the unigenes were > 3000 bp. The number of unigenes decreased as the gene length increased (Fig. 1a). RPKM (reads per kilobase per million mapped reads) values was used to evaluate the expression levels of the unigenes. Seventy percent of the unigenes RPKM values were less than 1.0, and only 2% were greater than 100.0 (Fig. 1b), which showed that most unigenes had low expression.

Functional Annotation Of Unigenes

The unigene sequences were mapped to several protein databases, including Nr, Swiss-Prot, Pfam, eggNOG, COG, GO and KEGG, using BLASTx algorithm-based software with an E-value cut-off of 10^− 5. Approximately 56,994 (58.53%) unigenes were matched to a sequence in at least one database. In this process, the number of unigenes identified in the COG database was the largest, accounting for 48.24% (46,972) of the total unigenes, followed by the Nr (43,213, 44.38%), Swiss-Prot (39,006, 40.06%), GO (36,836, 37.83%), KEGG (32,958, 33.85%) and Pfam (28,699, 29.47%) databases. Only 23,321 (23.95%) unigenes displayed significant homologies in the eggNOG database (Table 1). Based on the BLAST hits in the Nr database, the E-values of the unigenes (50%) ranged from 1e^− 45 to 1e^− 5, and the other half was less than 1e^− 45 (Fig. 2a). In comparison to other species, P. sibirica showed the highest similarity to sequences from Picea sitchensis (11787) followed by Quercus suber (5967), but similarities to sequences from other species were also observed (Fig. 2b).

GO classification can determine the functional categories of unigenes, predict the biological processes they participate in or identify their cellular location. Using Blast2GO, 36,836 unigenes were assigned to 51 functional terms belonging to three categories, biological process, cellular component and molecular function. For biological process, dominant terms were cellular process (28,281, 76.78%), metabolic process (25,261, 68.58%), biological regulation (12,235, 33.21%), response to stimulus (10,597, 28.77%), signaling (3,975, 10.79%); for molecular function, the dominant terms were binding (25,573, 69.42%), catalytic activity (20,572, 55.85%), transporter activity (3,134, 8.51%) and transcription regular activity (1,693, 4.60%); for cellular component, most unigenes were located in cell (31,173, 84.63%) and cell part (31,118, 84.48%), 23,737 (64.45%) unigenes were located in organelle and a small number of unigenes were located in membrane (Fig. 3; Additional file 2).

The COG database is used for homologous protein annotation. We grouped 46,972 unigenes into 24 functional classifications, and the top 10 were (S) “Function unknown” (11,048, 23.52%), (T) “Signal transduction mechanisms” (4,343, 9.25%), (O) “Posttranslational modification, protein turnover, chaperones” (4,175, 8.89%), (G) “Carbohydrate transport and metabolism” (3007, 6.40%), (K) “Transcription” (2,967, 6.32%), (J) “Translation, ribosomal structure and biogenesis” (2,890, 6.15%), (C) “Energy production and conversion” (2,646, 5.63%), (U) “Intracellular trafficking, secretion, and vesicular transport” (2,226, 4.74%), (E) “Amino acid transport and metabolism” (1,753, 3.73%), and (P) “Inorganic ion transport and metabolism” (1,742, 3.71%). However, the smallest functional classifications were (Y) “Nuclear structure” (118, 0.25%) and (N) “Cell motility” (49, 0.10%) (Fig. 4; Additional file 3).

Screening And Classification Of Differentially Expressed Genes (degs)

DEGseq [20] was used to identify DEGs between the cold treatment groups and control group, which conformed to |log2FC|>1.0 and FDR < 0.05. In total, 871, 1,397 and 872 genes were significantly differentially expressed following exposure to cold for 6 h, 24 h and 48 h, respectively. Of the DEGs, 418, 756 and 460 genes were upregulated, while 453, 641 and 412 genes were downregulated at each time point. In the pairwise comparisons of the 6 h, 24 h and 48 h samples, few DEGs were identified. Among them, “24 h vs 6 h” had the smallest number, with only 17 DEGs (9 upregulated, 8 downregulated). There were 401 (224 upregulated, 177 downregulated) and 505 (248 upregulated, 257 downregulated) in the ‘48 h vs 6 h’ and “48 h vs 24 h” comparisons, respectively. Approximately 70% of the DEGs were upregulated or downregulated less than 10-fold.

GO enrichment were performed to investigate the function of the DEGs. GO terms with corrected P-values < 0.05 were identified as significantly enriched. The top 30 significantly enriched GO terms under different stress time are shown in Additional file 4. For biological process, the important GO terms were “ethylene-activated signaling pathway”, “phosphorelay signal transduction system”, “hormone-mediated signaling pathway”, “secondary metabolic process”, “signal transduction”, “signaling”, “response to stimulus”, “intracellular signal transduction” and “calcium-mediated signaling”; for cellular component, the important GO terms were “extracellular region”, “apoplast”, “cell wall” and “external encapsulating structure”; and for molecular function, “terpene synthase activity”, “ADP binding”, “DNA binding transcription factor activity”, “transferase activity, transferring hexosyl groups”, “calcium: cation antiporter activity” and “ glutamate receptor activity” were the important GO terms.

Pathway analysis of DEGs was applied to predict intracellular signaling and metabolic pathways. The pathways with a corrected P-value < 0.05 were identified as significantly enriched. All the significantly enriched pathways of DEGs in P. sibirica under different cold stresses are shown in Additional file 5. The important pathways of DEGs included “spliceosome”, “RNA degradation”, “mRNA surveillance pathway”, “RNA transport”, “biosynthesis of amino acids”, “carbon metabolism”, “pyruvate metabolism” and “basal transcription factors”.

Quantitative Real-time Pcr (qrt-pcr) Analysis

To verify the accuracy of RNA-Seq for P. sibirica, 12 DEGs were randomly selected (6 upregulated, 6 downregulated) for qRT-PCR analysis with specific primers (Table 2). The results showed that expression patterns of 12 unigenes detected via qRT-PCR were highly consistent with the RNA-Seq results (Fig. 5), which suggested that the high-throughput RNA-Seq data were reliable and demonstrated that the DEGs identified based on transcriptome sequencing were available.

Changes In Physiological Indices Under Cold Stress

Physiological indices related to membrane system, including membrane permeability, REC, ROS and MDA, increased significantly (P < 0.05) in response to cold stress. Similarly, physiological indices involved in antioxidant and osmotic regulation, including POD activity, CAT activity, soluble sugar, soluble protein and proline, also increased significantly (P < 0.05) under cold stress. With extension of the cold stress duration, the 9 physiological indices generally showed a trend toward increasing first and then decreasing, with a maximum value at 24 h and decreasing at 48 h (Fig. 6).

Physiological changes in photosynthetic characteristics of P. sibirica under cold stress are shown in Fig. 7. The net photosynthetic rate, stomatal conductance and transpiration rate in P. sibirica seedlings dropped sharply (P < 0.05) in response to cold stress, which were also decreased significantly (P < 0.05) with extension of the stress time at -20 ℃.

Identification of DEGs involved in cold signal sensing and transduction

The DEGs involved in cold signal sensing and transduction are shown in Table 3. A total of 16 genes associated with ABA were notably induced by cold stress, 11 of which were upregulated and mainly encoded ABA 8’-hydroxylase, ABA receptor in response to ABA or participating in ABA biosynthetic process, and 5 were downregulated, mainly in response to ABA. However, it is worth noting that gene 1601162 negatively regulated the abscisic acid-activated signaling pathway. Six genes related to Ca²⁺ were mainly involved in calcium ion binding and transduction. Two calcineurin genes (1 upregulated and 1 downregulated), 2 calmodulin genes (1 upregulated and 1 downregulated) and 1 CBL-interacting protein kinase genes (upregulated) were also induced under cold stress.

Identification of DEGs encoding TFs in response to cold stress

In total, 36 TF genes were identified as DEGs and encoded 5 TFs, including AP2 (ethylene responsive factor), MYB (myeloblastosis), NAC (NAM, ATAF1, ATAF2 and CUC2), ZFP (Zinc finger protein) and bHLH (basic Helix-loop-helix) (Table 4). Of the 16 DEGs encoding AP2 or AP2 receptors, 11 were upregulated and 5 were downregulated. Ten genes encoding MYB changed prominently under cold stress, 6 of which were upregulated, including MYB106, MYBS3, and MYB73. There were also 6 DEGs (3 upregulated, 3 downregulated) encoding ZFP, 2 DEGs (upregulated) encoding NAC, and 2 DEGs (1 upregulated, 1 downregulated) encoding bHLH. With extension of the cold stress time, expression of the genes 1545648 (AP2), 1557638 (AP2), 1574963 (AP2) and 1601328 (ZFP) increased significantly, while expression of the genes 1575658 (MYB), 1583531 (MYB) and 1595496 (NAC) first increased and then decreased, with the highest expression level at 24 h followed by a decrease at 48 h.

Identification of DEGs involved in antioxidative regulation, osmotic regulation and photosynthesis

In this study, 11 DEGs were involved in antioxidative regulation, and 9 DEGs were involved in osmotic regulation to resist cold stress; there were also 7 DEGs involved in photosynthesis (Table 5). Two POD genes (1 upregulated, 1 downregulated), 1 APX (L-ascorbate peroxidase) gene (downregulated), 1 GPX (glutathione peroxidase) gene (upregulated) and 1 CAT gene (upregulated) were identified. GLR (glutamate receptor) and glutaredoxin acted as antioxidants, and the corresponding DEGs were identified. Small molecules, such as proline, soluble sugar, and soluble protein, help mitigate the damage of cold stress by increasing the osmotic pressure of plant cells. One gene encoding proline was significantly upregulated. Eight genes involved in carbohydrate binding and carbohydrate metabolism showed significant changes. With extension of the cold stress time, expression of the genes 1507086 (POD), 1546868 (CAT), 1594284 (proline dehydrogenase 1), 1544490 (carbohydrate metabolism), 1581549 (carbohydrate metabolism), 1585015 (carbohydrate metabolism), and 1591782 (carbohydrate metabolism) first increased and then decreased, demonstrating the highest expression level at 24 h followed by a decrease at 48 h. There were 7 DEGs (3 upregulated, 4 downregulated) associated with photosynthesis in response to cold stress. Their function was mainly to regulate the chlorophyll biosynthetic process, or to encode chlorophyllase and some components of Photosystems I and II.

Discussion

Pinus sibirica is a conifer species with strong cold tolerance. Without a public genome database, detailed information on the genes involved in the cold response are unavailable for P. sibirica. DEG analysis is a good tool to study gene temporal expression regulation [21]. In this study, 12 P. sibirica cDNA libraries of four different cold stress times (three biological replicates per stress time) were constructed for Illumina RNA-sEq. Whole-transcriptome analysis of P. sibirica in response to cold stress was performed. For the 12 samples, more than 52.2 Gb of reads were obtained, which were filtered and de novo assembled into 97,376 unigenes with a N50 of 1,362 bp. Of the unigenes, 58.53% had homologs in at least one of the public databases that we searched. According to the results of the Nr database match, P. sibirica had the strongest homology to Picea sitchensis. Taken together, these results suggested that our P. sibirica EST dataset provided a more adequate transcriptome resource for cold-responsive gene discovery and functional analysis in conifer species. It has also been demonstrated that cold resistance can be enhanced in plants by gene induction and repression in response to cold stress [22, 23]. Some of the P. sibirica DEGs had no annotated homologs in the public database, which indicated that these DEGs might be specific to P. sibirica or that the cold-responsive genes with homologs have not been identified in previous studies with other plant species.

Cold Stress Sensing And Signaling

Thus far, the identification of cold signal sensor in plant cells is unknown. However, cold stress can quickly induce membrane rigidity at microdomains, subsequently influencing protein folding and changing the physical state of the membrane, which are expected to alter metabolic reactions. Therefore, plant cells can sense cold stress through membrane rigidification, protein/nucleic acid conformational changes, and/or metabolite concentration changes [3]. According to a study of physiological indices in P. sibirica, cold stress improved the membrane permeability, leading to cell electrolyte leakage of salt and organic acids into the environmental media, so the REC (relative conductivity) increased [6]. Moreover, the breakdown of membrane permeability led to the accumulation of ROS (reactive oxygen species), which destroyed carbohydrates, some proteins and nucleic acids, causing membrane peroxidation [24, 25]. MDA (malonaldehyde) was the final product of membrane peroxidation, which was not only toxic but also acted as an indicator of membrane damage [26]. These results indicated that cold stress led to damage to the membrane of P. sibirica, changing its physiological state and the conformation of its associated proteins and nucleic acids, thereby sensing the cold signal.

Cytosolic Ca²⁺ and ABA levels can act as second messengers of the cold signal [27, 28]. Calcium accumulation and release in plant cells are mainly achieved through calcium ion channels and maintaining a stable calcium ion concentration is a prerequisite for normal plant development. Cold stress-induced second messenger signatures can be decoded by different pathways. Calcium signatures can be bound by calcium receptor family proteins, such as calmodulins (CaMs), calcineurin, calcineurin B-like (CBL) proteins, CBL-interacting protein kinase (CIPK), or calcium ion binding proteins [29, 30]. The receptor proteins transmit the cold signal and then activate the expression of related genes to induce physiological and biochemical changes in the cell to relieve the stress [31]. In this research, the expression of genes associated with calcium signaling-encoding calcium receptor family proteins changed significantly under cold stress, most of which were upregulated.

ABA regulates multiple processes of plant growth and development, such as seed germination and dormancy, leaf senescence and loss, and seedling growth. In addition, ABA plays vital roles in the adaption of plants to abiotic stress such as cold stress [32]. ABA receptors perceive ABA as a second messenger signature for signal transduction, which triggers an intracellular downstream signaling cascade to finally activate biochemical and physiological responses [33]. In this study, 16 DEGs associated with an ABA pathway were identified in P. sibirica in response to cold stress, 11 of which were prominently upregulated. Their main functions were to encode ABA 8’-hydroxylase and ABA receptors, respond to ABA or participate in the ABA biosynthetic process. Of the 5 downregulated DEGs in response to cold stress, gene 1601162 negatively regulated ABA, which promoted the activation of the ABA signaling pathway in P. sibirica. Other studies have shown that ABA can induce the production of cyclic ADP ribose (cADPR), which controls the release of intracellular calcium storage in plants [34]. Thus, ABA can control the calcium signaling pathway. Both the ABA and calcium signaling pathways were involved in the regulating the cold signal in P. sibirica.

Transcriptional Regulation

Transcriptional regulation plays an important role in the process of resisting cold stress. Transcription factors (TFs) can identify the functional elements on downstream gene promoters to activate or inhibit the expression of target genes, generate the corresponding biochemical and physiological reactions, and then enhance the ability to adapt to cold environment in plants [35, 36]. Low temperature can induce the expression of ethylene-responsive element binding factor/APETALA2 (ERF/AP2)-type transcription factor family genes [37]. DREB/CBF is a member of the AP2/ERF family [37]. Inducer of CBF expression 1 (ICE1) encodes a MYC-type basic-loop-helix (bHLH) transcription factor, can bind to MYC recognition elements in the CBF3 promoter, and induces the expression of CBF3 during cold acclimation [3]. ICE-CBF is a crucial regulatory pathway in response to cold stress. In this study, the expression of genes encoding AP2 and bHLH TFs changed significantly under cold stress, so it was speculated that the ICE-CBF pathway was present in P. sibirica in response to cold stress.

In addition to CBFs, several classes of TFs also play an important role in cold acclimation. Soybean C2H2-type zinc finger protein SCOF-1 were overexpressed in tobacco, and the expression of COR genes and cold tolerance were enhanced in transgenic tobacco plants [38]. The plant-specific factor NAC family plays a key role in environmental stress responses. The cold tolerance of Arabidopsis thaliana was enhanced with TaNAC gene overexpression, and several stress-regulated genes were induced in TaNAC-overexpressing Arabidopsis plants [39]. MYB is the largest family of transcription factors, and its role in cold regulation is a complex process. Overexpression of the cold-regulated rice TF OsMYB4 (an R2R3-type MYB) enhanced cold tolerance in Arabidopsis [40]. Nevertheless, overexpression of TF MYB15 reduced cold tolerance in Arabidopsis through reduced expression of CBF genes [41]. All the DEGs encoding TFs in P. sibirica played a crucial role in the process of resisting cold stress, which might be the main reason for the strong cold tolerance of P. sibirica.

Antioxidant And Osmotic Regulation

Under cold stress, plants perceive the cold signal and transfer it downstream to activate transcriptional regulation, induce target gene expression, and synthesize various antioxidant enzymes and antioxidants to neutralize the excess ROS [42]. Simultaneously, some osmotic regulatory substances are also synthesized in plant cells to maintain the osmotic pressure of the cells, thus protecting the plant from the harm caused by cold stress [42]. In the present study, POD (peroxidase) and CAT (catalase) activity in P. sibirica increased significantly under cold stress, in addition to some small molecular substances such as soluble sugar, soluble protein and proline.

The upregulated POD, CAT and GPX (glutathione peroxidase) genes can improve the cold tolerance of plants [43], which suggests that the upregulated genes encoding the three enzymes relieved cold stress in P. sibirica. Glutamate is an important damage signal in plants, and glutamate receptor-like protein family members act as plant damage receptors, which can induce an increase in the calcium ion concentration in plant cells and inhibit stress by regulating calcium ion inflow [44]. The upregulated glutamate receptor genes indicated that cold stress caused damage to P. sibirica, and P. sibirica alleviated and inhibited the damage by regulating the high expression of glutamate receptor genes. Glutaredoxin is a glutathione-dependent redox enzyme, which is widely found in plants and plays an important role in plant growth and development, as well as in resisting adversity [45]. In this study, 2 significantly upregulated glutaredoxin DEGs played an important role in resisting the oxidative stress induced by cold stress in P. sibirica.

Under cold stress, the accumulation of some small molecular substances in plant cells, such as proline, soluble sugar, soluble protein and other carbohydrates, can enhance the osmotic pressure of cells, thus alleviating the damage caused by cold stress [46]. In this study, 8 DEGs encoding proline and substances involved in carbohydrate binding and metabolism were significantly upregulated, which indicated that P. sibirica could resist cold stress damage by increasing the osmotic regulatory substances in the cells. With extension of the cold stress time, expression of the genes 1507086 (POD), 1546868 (CAT), 1594284 (proline dehydrogenase 1), 1544490 (carbohydrate metabolism), 1581549 (carbohydrate metabolism), 1585015 (carbohydrate metabolism), and 1591782 (carbohydrate metabolism) first increased and then decreased, with the highest expression level observed at 24 h followed by a decrease at 48 h. These results supported the measurements of physiological indicators.

Photosynthesis Under Cold Stress

Photosynthesis is an essential metabolic process for plant growth and development, and it is very sensitive to cold stress [47]. In this study, the net photosynthetic rate and stomatal conductance were inhibited, and more DEGs associated with photosynthesis were downregulated. Cold stress inhibited some genes related to the regulation of the chlorophyll biosynthetic process [48]. The decline in photosynthesis is affected by the suppression of photosynthesis-related enzyme activity [49]. The gene related to the activity of chlorophyllase activity was repressed in P. sibirica. The photochemical reaction of photosynthesis comprises two optical systems, PS I and PS II, which play a decisive role in the photosynthetic capacity of plants [23]. The genes encoding photosystem II reaction center X protein (PsbX) and Photosystem I reaction center subunit VI-2 were downregulated in P. sibirica in response to cold stress. Even plants with strong cold tolerance demonstrate an inhibition of photosynthesis under cold stress [50].

Conclusions

This study provided a comprehensive description of transcriptomic responses to cold stress in the leaves of P. sibirica with different stress times. Based on physiological data and transcriptomic DEGs data, a cold response model in P. sibirica was roughly summarized. Cold signal was first perceived by the plant cell membrane, with membrane rigidity at microdomains, and then the cold signal was transmitted downstream through Ca²⁺ and ABA signaling pathways to activate relevant TFs. The ICE-CBF pathway played an important role in the transcriptional regulation of cold stress, as well as the CBF-independent pathway. The action of the TFs triggered a cascade of downstream COR gene, synthesizing various antioxidant enzymes and osmotic protective substances to modulate cellular metabolism homeostasis and improve tolerance to cold stress. Photosynthesis is an essential metabolic process in plants, and it was affected by cold stress even in P. sibirica with strong cold tolerance. The rough cold response model and DEGs involved in the model may facilitate studies of molecular breeding improvement of conifer species.

Methods

Cdna Library Preparation And Illumina Sequencing For Transcriptome Analysis

The total RNA was extracted using a total RNA extraction kit (TaKaRa, Beijing, China) according to the instruction manual. RNA quality was verified by RNase-free agarose gel electrophoresis, and the concentration was measured using a biological analyzer (2100, Agilent, USA) at 260 nm and 280 nm. The RNA with 260 nm/280 nm ratios between 1.8 and 2.0 was used for subsequent experiments. A total of 12 cDNA libraries were constructed, which were sequenced using a DNA sequencer (HiSeq™ 2000, Illumina, USA), and the sequencing strategy was PE150.

De novo assembly and functional annotation of Illumina sequencing

Low-quality reads (more than 5% ambiguous nucleotides, or more than 15% bases with a Q-value ≤ 19) and adaptor-containing reads were filtered and removed. Clean reads were obtained for de novo assembly. Trinity assembly software was used, which overlapped the clean reads to generate contigs (contiguous sequences) [52]. Next, the reads were compared back to the contigs. The distance between different contigs from the same transcript was determined according to paired-end reads. The contigs were linked together by Trinity to obtain the unigene sequences that could not be extended at both ends. All the assembled P. sibirica unigenes were aligned to several protein databases using the BLASTX algorithm with an E-value < 1e^− 5, such as the Swiss-Prot protein, Protein families (Pfam), evolutionary genealogy of genes: Nonsupervised Orthologous Groups (eggNOG) protein, Clusters of Orthologous Groups of proteins (COG), Gene Ontology (GO) protein and Kyoto Encyclopedia of Genes and Genomes (KEGG) protein databases. The Blast2GO program was used to obtain the gene ontology (GO) annotation, and the KEGG database was used to analyze related gene functions in cellular processes and some protein products of metabolic processes [53].

Identification And Annotation Of Differentially Expressed Genes (degs)

Reads associated with each unigene were mapped to the transcriptome using the alignment software Bowtie 0.128. Unigene expression was calculated by counting and normalizing the number of mapped reads of each unigene into a reads per kb per million reads (RPKM) value. The thresholds were set with a false discovery rate < 0.001 and an absolute value of log₂ ratio > 1 to identify significant differences [54]. GO function analysis and KEGG pathway analysis were performed for the DEGs [53].

Quantitative Real-time Pcr (qrt-pcr) Validation

The 12 total RNA were extracted from the samples, and 1 µg RNA from each sample was reverse-transcribed to synthesize cDNA using the ReverTre Ace®qPCR RT Kit (Toyobo, Osaka, Japan). Each of the generated cDNAs was diluted 10 times as the qRT-PCR template. qRT-PCR was performed with a DNA Engine OpticonTM 2 Real-Time System (Bio-Rad, USA), and the reaction mixture was composed of 10 µl 2 × SYBR Green Realtime PCR Master mix (Toyobo, Osaka, Japan), 2.5 µl cDNA, 0.5 µl upstream primer, 0.5 µl downstream primer and deionized water in a final volume of 20 µl. The PCR was conducted at 94 ℃ for 30 s, followed by 45 cycles of 94 ℃ for 12 s, 54 ℃ for 30 s and 72 ℃ for 30 s. The expression levels of the selected genes were determined using the 2^−ΔΔCt algorithm with the P. sibirica Tubulin alpha (TUBA) gene as an internal control [55]. Deionized water was used as the no-template control. Each sample comprised three biological replicates, and the data are presented as means ± standard errors (SE) (n = 3). The primer sequences of the selected genes are listed in Table 2.

Measurement Of Physiological Indices

Membrane permeability and relative conductivity (REC) were measured as described by Daneshmand et al. [56] and Meng et al. [57], respectively. Reactive oxygen species were determined by the method of Wang and Jiao [58]. Malonaldehyde (MDA) was measured using a thiobarbituric acid test [59]. Peroxidase (POD) activity was determined according to the absorbance at 470 nm [60]. Catalase (CAT) activity was determined by measuring the decomposition of H₂O₂ at 240 nm as reported by Nazari et al. [61]. Soluble sugar and soluble protein were measured as described by Moriyama et al. [62] and Bradford [63], respectively. Proline was determined by the method according to Bates et al. [64].

After cold stress treatment, the seedlings were restored for 30 minutes in the greenhouse, and then the needles on the second round of branches for each seedling were selected to evaluate various photosynthetic parameters using a photosynthetic apparatus (LI-6400XT, Lico, USA). During measurement, the light intensity and CO₂ were set as 1600 µmol·m^− 2·s^− 1 and 400 µmol·mol^− 1, while other environmental factors were not controlled. The output data of the instrument mainly included the instantaneous net photosynthetic rate (mol·m^− 2·s^− 1), transpiration rate (mol·m^− 2·s^− 1), and stomatal conductance (mol·m^− 2·s^− 1).

Abbreviations

ABA: abscisic acid; CBF: C-repeat binding factors; DREB1: dehydration-responsive element-binding factors 1; TF: transcription factor; ICE: Inducer of CBF Expression; COR: cold-regulated; REC: relative conductivity; ROS: reactive oxygen species; MDA: malonaldehyde; POD: peroxidase; CAT: catalase; PE: paired-end; RPKM: reads per kilobase per million mapped reads; DEGs: differentially expressed genes; qRT-PCR: quantitative real-time PCR; AP2: ethylene responsive factor; MYB: myeloblastosis; NAC: NAM, ATAF1, ATAF2 and CUC2; ZFP: Zinc finger protein; bHLH: basic Helix-loop-helix; APX: L-ascorbate peroxidase; GPX: glutathione peroxidase; GLR: glutamate receptor.

Declarations

Ethic approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All transcriptome data sets used and/or analysed during the current study are available from the in the US National Library of Medicine, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA613137.

Competing interests

The authors declare that the research was conducted in absence of any competing interests.

Funding 

Study design，the labor service, travel expense, data analysis and interpretation were supported by the Opening Project of State Key Laboratory of Tree Genetics and Breeding in China (No. K2019103). The funding body did not have any role in writing of the manuscript.

Authors ' contributions

FW analyzed the datasets and wrote the manuscript. SC, ZL, YY and MT performed the experiments and helped interpret the datasets. XZ conceived and planned the experiments and reviewed and commented on the original manuscript. All authors read and approved the final manuscript.

Acknowledgments

We are grateful for the experimental platform provided by the state key laboratory of tree genetics and breeding, and for the experimental help provided by the laboratory students.

References

1. Alantsev NK. The stone pine. Moscow: Forestry Industry Press. 1981;p. 60-61.
2. Wang C. Study on the introduction and seed origin experiment of Pinus sibirica. Thesis for M.S., Northeast Forestry University, Supervisor: Jiang J. 2011;p. 1-7.
3. Chinnusamy V, Zhu JK, Sunkar R. Gene regulation during cold acclimation in plants. Methods in Molecular Biology. 2010;639:39-55.
4. Thomashow, Michael F. Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu. Rev. Plant Physiol. Plant Mol Biol. 1999;50:571–599.
5. Almadanim MC, Gonçalves NM, Rosa MTG, Alexandre BM, Cordeiro AM, Rodrigues M, et al. The rice cold-responsive calcium-dependent protein kinase OsCPK17 is regulated by alternative splicing and post-translational modifications. Biochimica Et Biophysica Acta. 2018;1865(2):231.
6. Wang F, Liang D, Pei X, Zhang Q, Zhang P, Zhang J, et al. Study on the physiological indices of Pinus sibirica and Pinus koraiensis seedlings under cold stress. J Forestry Res. 2018;30:1-11.
7. Mittler R. Oxidative stress, antioxidants and stress tolerance. Trends in Plant Science. 2002;7(9):405-10.
8. Paredes M, Quiles MJ. The effects of cold stress on photosynthesis in hibiscus plants. PloS One, 2015;10(9):e0137472.
9. Hirayama T, Shinozaki K. Research on plant abiotic stress responses in the post-genome era: past, present and future. Plant J. 2010;61(6):1041-52.
10. Matsui A, Ishida J, Morosawa T, Mochizuki Y, Kaminuma E, Endo TA et al. Arabidopsis transcriptome analysis under drought, cold, high-salinity and aba treatment conditions using a tiling array. Plant Cell Physiol. 2008;49(8):1135-1149.
11. Chen SY, Huang X, Yan XQ, Liang Y, Wang YZ, Li XF, et al. Transcriptome analysis in sheepgrass (Leymus chinensis): a dominant perennial grass of the Eurasian Steppe. PLoS One. 2013;8(7), e67974.
12. Abeynayake SW, Byrne S, Nagy I, Jonavičienė K, Etzerodt TP, Boelt B, et al. Changes in Lolium perenne transcriptome during cold acclimation in two genotypes adapted to different climatic conditions. BMC Plant Biol. 2015;15:250-264.
13. Beike AK, Lang D, Zimmer AD, Wust F, Trautmann D, Wiedemann G, et al. Insights from the cold transcriptome of Physcomitrella patens: global specialization pattern of conserved transcriptional regulators and identification of orphan genes involved in cold acclimation. New Phytol. 2015;205(2):869–81.
14. Wang F, Chen S, Liang DY, Qu GZ, Chen S, Zhao XY. Transcriptomic analyses of Pinus koraiensis under cold stresses. BMC Genomics. 2020;21:10.
15. Ito Y, Katsura K, Maruyama K, Taji T, Kobayashi M, Seki M, et al. Functional analysis of rice DREB1/CBF-type transcription factors involved in cold-responsive gene expression in transgenic rice. Plant Cell Physiol. 2006;47(1):141-53.
16. Wang Y, Jiang CJ, Li YY, Wei CL Deng W. CsICE1 and CsCBF1: two transcription factors involved in cold responses in Camellia sinensis. Plant Cell Rep. 2012;31:27-34.
17. Hu Y, Jiang L, Wang F, Yu D. Jasmonate regulates the inducer of cbf expression-C-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. The Plant Cell. 2013;25 (8): 2907-24.
18. Chinnusamy V, Zhu J, Zhu JK. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12(10):444–51.
19. Tian Y, Zhang H, Pan X, Chen X, Zhang Z, Lu X, et al. Overexpression of ethylene response factor TERF2 confers cold tolerance in rice seedlings. Transgenic Res. 2011;20(4):857-866.
20. Wang L, Feng Z, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq Data. Bioinformatics. 2010;26: 136-8.
21. Hao J, Guo H, Shi X, Wang Y, Wan Q, Song Y, et al. Comparative proteomic analyses of two Taxus species (Taxus × media and Taxus mairei ) reveals variations in the metabolisms associated with paclitaxel and other metabolites. Plant Cell Physiol. 2017;58:1878-90.
22. Zhan XQ, Yang L, Wang D, Zhu JK, Lang ZB. De novo assembly and analysis of the transcriptome of Ocimum americanum var. pilosum under cold stress. BMC Genomics. 2016;17:209-21.
23. Fu JJ, Miao YJ, Shao LH, Hu TM, Yang PZ. De novo transcriptome sequencing and gene expression profiling of Elymus nutans under cold stress. BMC Genomics 2016;17:870-89.
24. Foyer CH, Noctor G. Redox homeostasis and antioxidant signaling: a metabolic interface between stress perception and physiological responses. Plant Cell. 2005;17:1866-75.
25. Senthil-Kumar M, Kumar G, Srikanthbabu V, Udayakumar M. Assessment of variability in acquired thermo tolerance: potential option to study genotypic responseand the relevance of stress genes. J Plant Physiol. 2207;164:111-25.
26. Rakei A, Maali-Amiri R, Zeinali H, Ranjbar M. DNA methylation and physio-biochemical analysis of chickpea in response to cold stress. Protoplasma. 2016;253, 61-76.
27. Knight MR. Signal transduction leading to low-temperature tolerance in Arabidopsis thaliana. Philosophical Transactions of The Royal Society B Biological Sciences. 2002;357, 871-5.
28. Shi Y, Yang S. ABA regulation of the cold stress response in plants. https://www.researchgate.net/publication/285919571_ABA_Regulation_of_the_Cold_Stress_Response_in_Plants. 2014;p. 89-108.
29. Klimecka M, Muszyńska G. Structure and functions of plantcalcium-dependent protein kinases. Acta Biochimica Polonica. 2007;54, 219-33.
30. Reddy VS, Reddy AS. Proteomics of calcium-signaling components in plants. Phytochemistry 2004;65:1745-76.
31. Sheen J. Ca²⁺-dependent protein kinase and stress signal transduction in plants. Science. 1996;274:1900-02.
32. Himmelbach A, Yang Y, Grill E. Relay and control of abscisic acid signaling. Curr Opin Plant Biol. 2003;6:470–479.
33. Venis M. Methods in receptor research. In: Venis M, editor. Hormone binding sites in plants. New York: Longman. 1985;p. 24-40.
34. Wu Y, Kuzma J, Marechal E, Graeff R, Lee HC, Foster R, et al. Abscisic acid signaling through cyclic ADP ribose in plants. Science. 1997;278, 2126-2130.
35. Ashraf M. Inducing drought tolerance in plants: recent advances. Biotechnol Adv. 2010;28:169-183.
36. Lee BH, Henderson DA, Zhu JK. The Arabidopsis cold-responsive transcriptome and its regulation by ICE1. Plant Cell. 2005;17:3155-75.
37. Bai B, Wu J, Sheng WT, Zhou B, Zhou LJ, Zhuang W, et al. Comparative analysis of anther transcriptome profiles of two different rice male sterile lines genotypes under cold stress. Int J Mol Sci. 2015;16:11398-416.
38. Kim JC, Lee SH, Cheong YH, Cheolcmin Y, Cho MJ. A novel cold-inducible zinc finger protein from soybean, scof-1, enhances cold tolerance in transgenic plants. Plant J. 2001;25, 247-59.
39. Mao X, Zhang H, Qian X, Li A, Zhao G, Jing R. TaNAC2, a NAC-type wheat transcription factor conferring enhanced multiple abiotic stress tolerances in Arabidopsis. J Exp Bot. 2012;63:2933-46.
40. Pasquali G, Biricolti S, Locatelli F, Baldoni E, Mattana M. Osmyb4 expression improves adaptive responses to drought and cold stress in transgenic apples. Plant Cell Rep. 2008;27:1677-86.
41. Agarwal M, Hao Y, Kapoor A, Dong CH, Fujii H, Zheng X. 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: 37636-45.
42. Gill SS, Tuteja N. Polyamines and abiotic stress tolerance in plants. Plant Signal Behav. 2010;5:26-33.
43. Baek KH, Skinner DZ. Alteration of antioxidant enzyme gene expression during cold acclimation of near-isogenic wheat lines. Plant Sci. 2003;165, 0-1227.
44. Sivaguru M, Ezaki B, He ZH, Tong HY, Osawa H, Baluska F, et al. Aluminum-induced gene expression and protein localization of a cell wall-associated receptor kinase in Arabidopsis. Plant Physiol. 2003;132: 2256-66.
45. Daniel L, Ema O, Paula S, Marcela S, Xavier J, Loreto H. (2012). Glutaredoxin GRXS13 plays a key role in protection against photooxidative stress in Arabidopsis. J Exp Bot. 2012;63:503-15.
46. Farhangi-Abriz S, Torabian S. Antioxidant enzyme and osmotic adjustment changes in bean seedlings as affected by biochar under salt stress. Ecotox Environ Safe. 2017;137:64-70.
47. Zheng YL, Feng YL, Lei YB, Yang CY. Different photosynthetic responses to night chilling among twelve populations of Jatropha curcas. Photosynthetica 2009;47:559-66.
48. Tewari AK, Tripathy BC. Temperature-stress-induced impairment of chlorophyll biosynthetic reactions in cucumber and wheat. Plant Physiol. 1998;117:851-8.
49. Allen DJ, Ort DR. Impacts of chilling temperature on photosynthesis in warm-climate plants. Trends Plant Sci. 2001;6:36-42.
50. Oquist G, Huner NPA. Photosynthesis of overwintering evergreen plants. Annu Rev Plant Biol. 2003;54:329-55.
51. Zhao XY, Wang C, Li SC, Hou W, Zhang SQ, Han GJ, et al. Genetic variation and selection of introduced provenances of Siberian Pine (Pinus sibirica) in frigid regions of the Greater Xing’an Range, northeast China. J Forestry Res. 2014;25(3):549–56.
52. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Trinity: reconstructing a full-length transcriptome without a genome from RNA-Seq data. Nat Biotechnol. 2011;29(7):644–52.
53. Chunmei H. Identification of genes involved in biosynthesis of mannan polysaccharides in dendrobium officinale by rna-seq analysis. Plant Mol Biol. 2015;3:219-231.
54. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 2010;26:139–140.
55. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using Real-Time quantitative PCR and the 2^-ΔΔCt method. Methods. 2001;25(4):402-408.
56. Daneshmand F, Arvin MJ, Kalantari KM. Physiological responses to NaCl stress in three wild species of potato in vitro. Acta Physiol Plant. 2010;32:91-101.
57. Meng P, Bai X, Li H, Song X Zhang X. Cold hardiness estimation of Pinus densiflora var. zhangwuensis based on changes in ionic leakage, chlorophyll fluorescence and other physiological activities under cold stress. J Forestry Res. 2015;26,641-649.
58. Wang SY, Jiao HJ. Scavenging capacity of berry crops on superoxide radicals, hydrogen peroxide, hydroxyl radicals, and singlet oxygen. J Agr Food Chem. 2000;48:5677-84.
59. Heidarvand L, Maali-Amiri R. Physio-biochemical and proteome analysis of chickpea in early phases of cold stress. J Plant Physiol. 2013;170:459-469.
60. Omran RG. Peroxide levels and the activities of catalase, peroxidase, and indoleacetic acid oxidase during and after chilling cucumber seedlings. Plant Physiol. 1980;65:407-408.
61. Nazari M, Amiri RM, Mehraban FH, Khaneghah HZ. Change in antioxidant responses against oxidative damage in black chickpea following cold acclimation. Russ J Plant Physl. 2012;59:183-189.
62. Moriyama U, Tomioka R, Kojima M, Sakakibara H, Takenaka C. Aluminum effect on starch, soluble sugar, and phytohormone in roots of Quercus serrata Thunb seedlings. Trees. 2016;30:405-413.
63. Bradford MM. A rapid and sensitive method for the quantitation of microgram quantities of protein utilizing the principle of protein-dye binding. Anal Biochem. 1976;72:248–54.
64. Bates LS, Waldren RP, Teare ID. Rapid determination of free proline for water-stress studies. Plant Soil. 1973;39:205-207.

Tables

Table 1 List of P. sibirica transcriptome annotations

Table 2 Primer sequences

Table 3 The DEGs involved in cold signal sensing and transduction

Table 4 The DEGs encoding TFs in response to cold stress in P. sibirica

Table 5 The DEGs involved in antioxidation mechanisms and photosynthesis under cold stress in P. sibirica

Additional Files

Additional file 1: Summary statistics of P. sibirica RNA-seq data (XLSX).

Additional file 2: GO classifications of P. sibirica transcriptome (XLSX).

Additional file 3: COG functional classification of P. sibirica transcriptome (XLSX).

Additional file 4: GO classification of DEGs under different stress time (XLSX).

Additional file 5: KEGG enrichment of DEGs under different stress time (XLSX).