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 Ca2+ 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.