Transcriptome sequencing and gene expression profiling of Pinus sibirica under different cold stresses.

Cold stress is a major abiotic factor that affects plant growth and geographical distribution. Pinus sibirica is extremely frigostable tree species. To understand the molecular mechanisms of cold tolerance by P. sibirica, physiological responses were analyzed and transcriptome profiling was conducted to the plants treated by cold stress. The physiological data showed that membrane permeability relative conductivity (REC), reactive oxygen species (ROS), malonaldehyde (MDA) content, peroxidase (POD) and catalase (CAT) activity, soluble sugar, soluble protein and proline contents were increased significantly (p < 0.05) in response to cold stress. Transcriptome analysis identified a total of 871, 1397 and 872 differentially expressed genes (DEGs) after cold treatment for 6 h, 24 h and 48 h at -20°C, respectively. The signaling pathway mediated by Ca2+ as a signaling molecule and abscisic acid pathways were the main cold signal transduction pathways in P. sibirica. The APETALA2/Ethylene-Responsive Factor (AP2/ERF) and MYB transcription factor families also play an important role in the transcriptional regulation of P. sibirica. In addition, many genes related to photosynthesis were differentially expressed under cold stress. We also validated the reliability of transcriptome data with quantitative real-time PCR. This study lays the foundation for understanding the molecular mechanisms related to cold responses in P. sibirica.

Introduction

Cold stress is a major abiotic factor that adversely affects plant growth and geographical distribution (Chinnusamy ). Thus, plants have evolved a series of regulatory mechanisms that trigger complex gene expression processes, which in turn trigger a cascade of physiological and biochemical responses to cold stress (Thomashow 1999). Physiological changes mainly include a sharp increase in Ca2+ and abscisic acid (ABA) concentrations (Almadanim ), changes in membrane lipid composition and permeability (Wang ), increases in antioxidant enzymes and antioxidants (Mittler 2002), alterations in plant stomata and photosynthesis (Paredes and Quiles 2015), and the accumulation of osmotic regulatory substances (Wang ).

The molecular mechanism of plant adaption to cold stress is very complex, in which changes in gene expression play a critical role (Hirayama and Shinozaki 2010). Several studies have revealed gene expression changes in response to cold stress using genome-wide transcriptome analysis (Abeynayake , Beike , Chen , Matsui , Wang ), and these changes highlight the importance of transcriptional regulation. For example, in Pinus koraensis transcriptome profiling, 96 differentially expressed genes (41 upregulated and 55 downregulated) were shown to be involved in transcriptional regulation (Wang ). Among these transcription factors, C-repeat binding factors (CBFs) and dehydration-responsive element–binding factors 1 (DREB1) transcription factors (TFs) play a key molecular switching role in cold regulation (Ito ). At present, there are two main cold-induced regulatory mechanisms in plants: ICE1 (inducer of cbf expression 1)-CBF (c-repeat-binding factor) and CBF-independent mechanisms (Hu , Wang ). CBF can activate cold response (COR) genes by binding to cis-acting elements in their promoters (Chinnusamy ). Most COR genes, such as blue copper protein, heat shock protein, osmotic regulators and antioxidant enzymes, protect plant cells from cold damage (Chinnusamy ). In rice seedlings, OsSODB (superoxide dismutase), OsTrx23 (thioredoxin), and OsLti6 (encoding membrane proteins) accumulate in response to cold (Tian ). In recent years, RNA-seq analysis has been widely used to reveal the molecular mechanism of biological processes (Gang , Wang , 2019b).

Although many studies have been conducted in several plant species, there is currently no research on the gene expression profile in response to cold stress in Pinus sibirica. The species is an evergreen conifer belonging to the family Pinaceae and is mainly distributed in Siberia, Russia (49°40ʹ~127°20ʹE, 46°40ʹ~68°30ʹN) (Alantsev 1981). In China, it is only distributed in the northwest of the Altai Mountains, Xinjiang. Due to its excellent wood and nutritious seeds, P. sibirica is widely used in furniture, construction, dried fruit and oil production, and it is extremely cold tolerant (Wang 2011). In some localities, the lowest winter temperature can fall to below –65°C, and –20°C is quite common in most areas where P. sibirica grows, but the species survives and can still grow normally (Wang 2011). To understand the molecular regulatory mechanisms of P. sibirica in tolerating low temperature limitations, exploring the gene expression patterns under cold stress and their relationship with physiological mechanisms is imperative.

Thus, the objective of this study was to elucidate transcriptomic adaptive mechanisms in P. sibirica to cold stress and explore the relationship between gene expression and physiological responses. For this purpose, 5-year-old seedlings showing healthy growth at room temperature (20°C) were exposed to –20°C for different time courses (6 h, 24 h and 48 h). Thereafter, RNA-seq was performed. In addition, the following physiological responses were measured: membrane permeability, relative conductivity (REC), reactive oxygen species (ROS), malonaldehyde (MDA) content, peroxidase (POD) and catalase (CAT) activities, soluble sugar, soluble protein, proline contents, net photosynthetic rate, stomatal conductance and transpiration rate. This study provides insights into the molecular mechanisms of tolerance 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 for future conifer breeding. It could also 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.

Materials and Methods

Plant material and cold treatments

Seeds of P. sibirica were purchased from Novosibirisk (82°55ʹE, 55°2ʹN), Russia, in 2010. The seeds were sown in 2011. After five years, 300 healthy and morphologically uniform seedlings were transferred to the laboratory for growth in a greenhouse (20°C, air humidity from 50% to 65%, 16 h light/8 h dark photoperiod, light intensity 150 μmol·m–2·s–1) for one month until they were fully adapted to the environment. Thereafter, 120 seedlings with good performance were immediately transferred from room temperature to –20°C, and the needles were harvested after 0 (CK), 6, 24 and 48 h, which were achieved by a Haier BC/BD-318HD freezer (10°C~–26°C). For each time point, the harvestable needles included three biological replicates, and each replicate consisted of a pool of needles from 10 independent seedlings. All harvestable needles were frozen immediately in liquid nitrogen and stored in a –80°C freezer for further studies.

Measurement of physiological indices

Membrane permeability and relative conductivity (REC) were measured as described by Daneshmand and Meng , respectively. Reactive oxygen species were determined by the method of Wang and Jiao (2000). Malonaldehyde (MDA) content was measured using a thiobarbituric acid test (Heidarvand and Maali-Amiri 2013). Peroxidase (POD) activity was determined according to the absorbance at 470 nm (Omran 1980). Catalase (CAT) activity was determined by measuring the decomposition of H2O2 at 240 nm as reported by Nazari . Soluble sugar and soluble protein contents were measured as described by Moriyama and Bradford (1976), respectively. The proline content was determined by the method described by Bates .

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 (including control samples) were selected to evaluate various photosynthetic parameters using a photosynthetic apparatus (LI-6400XT, Lico, USA). During measurement, the light intensity and CO2 were set as 1600 μmol·m–2·s–1 and 400 μmol·mol–1, respectively, 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). Changes in the physiological characteristics of P. sibirica under cold stress were examined by one-way ANOVA and Duncan’s multiple comparison if the ANOVA result was significant (p < 0.05). SPSS-17 statistical software was used (SPSS Inc., Chicago, IL, USA).

cDNA library preparation and Illumina sequencing for transcriptome analysis

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. 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 and sequenced using a DNA sequencer (HiSeqTM 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 (https://github.com/trinityrnaseq/trinityrnaseq/wiki) to overlap the clean reads to generate contiguous sequences (contigs) (Grabherr ). The reads were then 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 RefSeq Nonredundant (Nr) database (https://www.ncbi.nlm.nih.gov/refseq/), the Swiss-Prot protein (https://www.expasy.org/resources/uniprotkb-swiss-prot), Protein families (http://pfam.xfam.org/), Nonsupervised Orthologous Groups (eggNOG) protein (eggnog.embl.de/version_3.0/), Clusters of Orthologous Groups of proteins (COG) (https://www.ncbi.nlm.nih.gov/research/cog-project/), Gene Ontology (GO) protein (http://www.geneontology.org), and Kyoto Encyclopedia of Genes and Genomes (KEGG) protein (https://www.genome.jp/kegg/) 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 (He ). All transcriptome data sets were deposited in NCBI (National Center for Biotechnology Information), https://www.ncbi.nlm.nih.gov/bioproject/PRJNA613137.

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 reads per kb per million reads (RPKM) values. The thresholds were set with a false discovery rate <0.05 and an absolute value of log2 ratio >1 to identify significant differences (Robinson ). GO function analysis and KEGG pathway analysis were performed for the DEGs (He ).

Quantitative Real-Time PCR (qRT-PCR) validation

Twelve total RNAs 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 Real-time 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. PCR was conducted at 94°C for 30 s, followed by 45 cycles of 94°C for 12 s, 54°C for 30 s and 72°C 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 (Livak and Schmittgen 2001). Deionized water was used as the no-template control. Each sample comprised three biological replicates, and the data are presented as the means ± standard errors (SE) (n = 3). The primer sequences of the selected genes are listed in Table 1.

Table 1.

Primer sequences used in this study

Unigenes	Forward primer	Reverse primer
TUBA	CCAGTTTGTTGATTGGTGTCC	ACGGCTCTCTGAACCTTGG
1585055	CGCAGCAACAACAACACAGAGTTC	GTCTTTGTCTCTGCTTGGCTTAGC
1586991	TGCAGTAAGCATGGAAGCAACC	CTTTGTGGCAGAGCAAATGCCT
1606092	AGTGCCGGACCTTGGAATAATG	GGTGATGGCACTGTTGAGAATGG
1503936	GGCGATACGCTGGAGTACAATAAG	CCTTGTATTAACTGTTGCCGTCCC
1544490	GGATTAGGCTACACATCTGCCTG	GTTGAGGGTGGTTACGTGCAT
1573575	AGGGAATGCTTCGAGTGGGAG	GGAAAGGCACCGTCTGATGGA
1580119	GCCATCTCCTACACTACCTTGTTTC	CAATGAAGGCGGTGATGTCTCT
1563093	CCAATCACAAGTTGGCTCCCATC	CAATCTGGTCTCTCTGTTCGTCCT
1487288	GCCAAAGGCACTAGAGGAAACAAG	GAGACGAGCAAGCGAGAAGCAA
1591073	AGAACGTCCTGGATCACTTCTG	GGAAGCAGCAATAGCTCAACGATC
1558637	TTAGGGCAGGTGGAAGAGTAGAAC	GCTTATGGGTCTTCTTGGGCTTTCT
1601162	ATTCGAGCATGGCATCACTTCC	GAGCCTAACCATTGACCAACCATG

Results

Physiological responses to cold stress

Physiological responses related to the membrane system, including membrane permeability, REC, ROS and MDA content, increased significantly (p < 0.05) in response to cold stress. Similarly, physiological mechanisms involved in antioxidant and osmotic regulation, including POD activity, CAT activity, soluble sugar, soluble protein and proline contents, increased significantly (p < 0.05) under cold stress. With extension of the cold stress duration, the physiological responses generally showed an increasing trend first and then decreasing, with a maximum value at 24 h and decreasing at 48 h (Fig. 1). The net photosynthetic rate, stomatal conductance and transpiration rate in P. sibirica seedlings dropped sharply (p < 0.05) in response to cold stress and decreased significantly (p < 0.05) with extension of the stress time at –20°C (Fig. 2).

Fig. 1.

Changes in physiological indices in P. sibirica under cold stress. (a) Membrane permeability; (b) Relative electric conductivity (REC); (c) Reactive oxygen species (ROS); (d) Malonaldehyde content (MDA); (e) Peroxidase activity (POD); (f) Catalase activity (CAT); (g) Soluble sugar content; (h) Soluble protein content; and (i) Proline content. Each value represents the mean ± SE, and different letters indicate significant differences at the 0.05 level according to Duncan’s multiple range test. The X-axis shows the time (hours) of the –20°C treatment.

Fig. 2.

Photosynthetic characteristics of P. sibirica under cold stress. (a) Net photosynthetic rate; (b) Stomatal conductance; and (c) Transpiration rate. Each value represents the mean ± SE, and different letters indicate significant differences at the 0.05 level according to Duncan’s multiple range test. The X-axis shows the time (hours) of the –20°C treatment.

RNA sequencing and de novo assembly

In this study, after filtering out low-quality reads, 541,809,548 clean reads were obtained, which were integrated and assembled into 97,376 unigenes with a mean length of 771 bp. Detailed information is listed in Supplemental Table 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. 3a). It was found that 70% of the unigene RPKM values were less than 1.0, and only 2% were greater than 100.0 (Fig. 3b), which showed that most unigenes had low abundance.

Fig. 3.

Illumina sequencing of P. sibirica. (a) The length distribution of unigenes in P. sibirica. (b) The expression distribution of unigenes in P. sibirica. Purple, blue, red, and green represent the percentages of unigenes with RPKM values greater than 100, less than 1.0, between 1 and 10, and between 10 and 100 for all unigenes, respectively.

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 cutoff of 10–5. A total of 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 2). 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 (Supplemental Fig. 1a). In comparison to other species, P. sibirica showed the highest similarity to sequences from Picea sitchensis followed by Quercus suber, but similarities to sequences from other species were also observed (Supplemental Fig. 1b, Supplemental Table 2).

Table 2.

List of P. sibirica transcriptome annotations

Public databases	Number of unigene hits	Percentage (%)
Nr	43,213	44.38
Swiss-Prot	39,006	40.06
Pfam	28,699	29.47
eggNOG	23,321	23.95
COG	46,972	48.24
GO	36,836	37.83
KEGG	32,958	33.85
ALL	56,994	58.53

GO classification can be used to describe biological processes, functional categories, and cellular locations of unigenes. 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%), and 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. 4, Supplemental Table 3).

Fig. 4.

GO classifications of the P. sibirica transcriptome.

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

Fig. 5.

COG functional classification of the P. sibirica transcriptome.

Screening and classification of differentially expressed genes (DEGs)

DEseq was used to identify DEGs between the cold treatment groups and the control group, which conformed to |log2FC| > 1.0 and FDR < 0.05 (Wang ). After cold treatment for 6 h, 24 h and 48 h, 871, 1397 and 872 genes, respectively, showed significant differences in expression compared with the control group (Supplemental Table 5). 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 expressed (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 was performed to investigate the function of the DEGs. GO terms with corrected p < 0.05 were identified as significantly enriched. The top 30 significantly enriched GO terms under different stress times are shown in Supplemental Table 6. For biological process, the important GO terms were ethylene-activated signaling pathway, phosphorlay 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.

KEGG pathway analysis of DEGs was applied to predict intracellular signaling and metabolic pathways. The significantly enriched KEGG pathway for DEGs was identified with a corrected p-value < 0.05. All the significantly enriched pathways of DEGs in P. sibirica under different cold stresses are shown in Supplemental Table 7. 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. The results showed that the expression patterns of 12 unigenes detected via qRT-PCR were highly consistent with the RNA-Seq results (Fig. 6), which suggested that the high-throughput RNA-Seq data were reliable and demonstrated that the DEGs identified based on transcriptome sequencing were available.

Fig. 6.

qRT-PCR analysis of DEGs in P. sibirica under cold stress. Transcript levels of 12 randomly selected DEGs, including 6 upregulated (a) and 6 downregulated (b) DEGs. The Y-axis on the left shows the relative gene expression levels (2–ΔΔCt) analyzed by qRT-PCR (blue lines), while the Y-axis on the right shows the corresponding RNA-seq expression data (red dotted lines). The X-axis shows the time (hours) of the –20°C treatment.

Identification of DEGs involved in cold signal sensing and transduction

The DEGs involved in cold signal sensing and transduction are shown in Table 3. Studies have shown that ABA (abscisic acid) is the key plant hormone regulating plant responses to abiotic stresses such as drought, cold, heat and salt, so we identified the genes involved in ABA synthesis (Ding ). 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 of which were downregulated, mainly in response to ABA. However, it is worth noting that gene 1601162 negatively regulated the abscisic acid-activated signaling pathway. In addition, calcium is a key ion for signal transduction in response to various environmental stresses (Liu ), and we also identified genes associated with calcium ions. Six genes related to the Ca2+ signaling pathway were mainly involved in the binding and transduction of calcium ions. Two calcineurin genes (1 upregulated and 1 downregulated), 2 calmodulin genes (1 upregulated and 1 downregulated) and 1 CBL-interacting protein kinase gene (upregulated) were also induced during cold stress.

Table 3.

DEGs involved in cold signal sensing and transduction

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), as shown in 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, the expression of the genes 1545648 (AP2), 1557638 (AP2), 1574963 (AP2) and 1601328 (ZFP) increased significantly, while the expression of the genes 1575658 (MYB), 1583531 (MYB) and 1595496 (NAC) first increased and then decreased, with the highest expression level after 24 h followed by a decrease after 48 h of cold stress.

Table 4.

DEGs encoding TFs in response to cold stress in P. sibirica

Identification of DEGs involved in antioxidative and osmotic regulation and photosynthesis

In this study, 11 DEGs were involved in antioxidative regulation, 9 DEGs were involved in osmotic regulation to resist cold stress and 7 DEGs were 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, the expression of the genes 1507086 (POD), 1546868 (CAT), 1594284 (proline dehydrogenase 1), 1544490, 1581549, 1585015, and 1591782 (carbohydrate metabolism) first increased and then decreased, demonstrating the highest expression level after 24 h followed by a decrease after 48 h of cold stress. 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.

Table 5.

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

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 is unavailable for P. sibirica. DEG analysis is a tool to study temporal gene expression regulation (Hao ). 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. Transcriptomics studies related to cold stress have been widely carried out in broadleaf tree species such as Taxus chinensis, Populus euphratica, and Hevea brasiliensis, but there are few transcriptomics studies on cold stress related to coniferous species such as P. sibirica. Therefore, the P. sibirica 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 (Fu , Zhan ). 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 cold-responsive genes with homologs have not been identified in previous studies with other plant species.

Cold stress sensing and signaling

Cold stress can quickly induce membrane rigidity at microdomains, subsequently influencing protein folding and changing the physical state of the membrane. Therefore, plant cells can sense cold stress through membrane rigidification, protein/nucleic acid conformational changes, and/or metabolite concentration changes (Chinnusamy ). A study on the physiological responses of P. sibirica to cold stress revealed increased membrane permeability, leading to cell electrolyte leakage of salt and organic acids into the environmental media and thus an increase in relative conductivity (REC) (Wang ). Moreover, the breakdown of membrane permeability led to the accumulation of reactive oxygen species (ROS), which destroyed carbohydrates, some proteins and nucleic acids, causing membrane peroxidation (Foyer and Noctor 2005, Senthil-Kumar ). In addition, MDA is the final product of membrane peroxidation, which is an indicator of membrane damage (Rakei ). Thus, the cell membrane of P. sibirica is first affected under cold stress, a large number of osmotic substances in the cell membrane flow out, and the metabolism of P. sibirica is imbalanced. In addition, the accumulation of reactive oxygen species and malondialdehyde also changes the membrane system of cell membranes.

Cold stress-induced second messenger signatures can be decoded by different pathways, such as cytosolic Ca2+ and ABA levels (Knight 2002, Shi and Yang 2014). 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. 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 (Klimecka and Muszyńska 2007, Reddy and Reddy 2004). Receptor proteins transmit cold signals and then activate the expression of related genes to induce physiological and biochemical changes in the cell to relieve stress (Sheen 1996). 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, and it plays vital roles in the adaption of plants to abiotic stresses such as cold stress (Himmelbach ). 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 (Venis 1985). In this study, 16 DEGs associated with the 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, the 1601162 gene 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 (Wu ). Thus, ABA can control the calcium signaling pathway. Both the ABA and calcium signaling pathways were involved in 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, thereby generating the corresponding biochemical and physiological reactions and enhancing the ability to adapt to cold environments in plants (Ashraf 2010, Lee ). Low temperature can induce the expression of ethylene-responsive element binding factor/APETALA2 (ERF/AP2)-type transcription factor family genes (Bai ), and DREB/CBF is a member of the AP2/ERF family (Bai ). 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 (Tian ). 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. For instance, soybean C2H2-type zinc finger protein SCOF-1 was overexpressed in tobacco, and the expression of COR genes and cold tolerance were enhanced in transgenic tobacco plants (Kim ). 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 (Mao ). 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 (Pasquali ). Nevertheless, overexpression of TF MYB15 reduced cold tolerance in Arabidopsis through reduced expression of CBF genes (Agarwal ). 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.

Physiological responses and expression of associated genes

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 excess ROS (Gill and Tuteja 2010). 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 (Gill and Tuteja 2010). In the present study, POD and CAT activities in P. sibirica increased significantly under cold stress, in addition to some small molecular substances such as soluble sugar, soluble protein and proline. We found that the expression of the genes 1507086 (POD), 1546868 (CAT), 1594284 (proline dehydrogenase 1), 1544490, 1581549, 1585015 and 1591782 (carbohydrate metabolism) first increased and then decreased, with the highest expression level observed after 24 h followed by a decrease after 48 h of cold stress. Upregulated POD, CAT and GPX (glutathione peroxidase) genes can improve the cold tolerance of plants (Baek and Skinner 2003). 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 (Sivaguru ). 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 that is widely found in plants and plays an important role in plant growth and development and in resisting adversity (Laporte ). In this study, two glutaredoxin DEGs were significantly upregulated and possibly 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 (Farhangi-Abriz and Torabian 2017). In this study, 8 DEGs encoding proline and substances involved in carbohydrate binding and metabolism were significantly upregulated. This indicates that P. sibirica could resist cold stress damage by increasing the osmotic regulatory substances in the cells.

In this study, the net photosynthetic rate and stomatal conductance were inhibited, and more DEGs associated with photosynthesis were downregulated. Photosynthesis is an essential metabolic process for plant growth and development, and it is very sensitive to cold stress (Zheng ). Cold stress inhibited some genes related to the regulation of the chlorophyll biosynthetic process (Tewari and Tripathy 1998). The decline in photosynthesis is attributed to the suppression of photosynthesis-related enzyme activity (Allen and Ort 2001). Our study showed that the gene related to chlorophyllase activity was repressed in P. sibirica. We also observed that 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. It should be noted that the photochemical reaction of photosynthesis comprises two optical systems, PS I and PS II, which play decisive roles in the photosynthetic capacity of plants (Fu ). Plants with strong cold tolerance even demonstrate an inhibition of photosynthesis under cold stress (Oquist and Huner 2003).

Conclusion

In this study, high-quality transcriptome data of P. sibirica were obtained, which reflected the gene expression changes of P. sibirica under cold stress. This database has a total of 52.2 Gb and contains 97376 assembled unigenes. The reliability and accuracy of P. sibirica transcriptome data were confirmed by qRT-PCR with 12 random primers. Under cold stress, the contents of REC, ROS, MDA and POD and other products increase in response to cold stress. In transcriptome studies, upregulated and downregulated DEGs work together to resist cold stress. A large number of DEGs involved in biological processes and molecular and cellular metabolic pathways have been identified. The main DEGs are involved in the recognition and transduction of low-temperature signals, transcriptional regulation and antioxidant processes. In addition, transcription factors such as AP2, MYB and NAC also regulate the expression of low temperature-related genes. This study provides valuable genetic resources for the discovery of cold tolerance genes and gene functions of P. sibirica and lays a foundation for genetic improvement of coniferous species.

Author Contribution Statement

FW was a major contributor in writing the manuscript. SC analyzed the data and make figures. KWC drafted the manuscript and substantially revised it. ZML participated in RNA extraction and performed RT-qPCR assay. YCY and MT participated in the design of the study and analyzed data. XYZ conceived of the study, participated in its design and data interpretation, and revised the manuscript critically.

Supplemental Figure

Supplemental Tables

Supplemental Table 5