- Research article
- Open Access
Transcriptomic analyses of Pinus koraiensis under different cold stresses
BMC Genomics volume 21, Article number: 10 (2020)
Pinus koraiensis is an evergreen tree species with strong cold resistance. However, the transcriptomic patterns in response to cold stress are poorly understood for P. koraiensis. In this study, global transcriptome profiles were generated for P. koraiensis under cold stress (− 20 °C) over time by high-throughput sequencing.
More than 763 million clean reads were produced, which assembled into a nonredundant data set of 123,445 unigenes. Among them, 38,905 unigenes had homology with known genes, 18,239 were assigned to 54 gene ontology (GO) categories and 18,909 were assigned to 25 clusters of orthologous groups (COG) categories. Comparison of transcriptomes of P. koraiensis seedlings grown at room temperature (20 °C) and low temperature (− 20 °C) revealed 9842 differential expressed genes (DEGs) in the 6 h sample, 9250 in the 24 h sample, and 9697 in the 48 h sample. The number of DEGs in the pairwise comparisons of 6 h, 24 h and 48 h was relatively small. The accuracy of the RNA-seq was validated by analyzing the expression patterns of 12 DEGs by quantitative real-time PCR (qRT-PCR). In this study, 34 DEGs (22 upregulated and 12 downregulated) were involved in the perception and transmission of cold signals, 96 DEGs (41 upregulated and 55 downregulated) encoding 8 transcription factors that regulated cold-related genes expression, and 27 DEGs (17 upregulated and 10 downregulated) were involved in antioxidant mechanisms in response to cold stress. Among them, the expression levels of c63631_g1 (annexin D1), c65620_g1 (alpha-amylase isozyme 3C), c61970_g1 (calcium-binding protein KIC), c51736_g1 (ABA), c58408_g1 (DREB3), c66599_g1 (DREB3), c67548_g2 (SOD), c55044_g1 (CAT), c71938_g2 (CAT) and c11358_g1 (GPX) first increased significantly and then decreased significantly with the extension of stress time.
A large number of DEGs were identified in P. koraiensis under cold stress, especially the DEGs involved in the perception and transmission of cold signals, the DEGs encoding TFs related to cold regulation and the DEGs removing ROS in antioxidation mechanisms. The transcriptome and digital expression profiling of P. koraiensis could facilitate the understanding of the molecular control mechanism related to cold responses and provide the basis for the molecular breeding of conifers.
Cold stress is one of the most important abiotic stresses that adversely affects plant growth and development, crop yield and quality, and geographic distribution . Plants have frequently suffered sudden cold stress, such as early or late frost in nature; subsequently, they would enable a diverse set of response mechanisms to protect against damage . The duration of stress is also a test of plant cold tolerance, which involves various cellular response mechanisms . Studies of the mechanisms that improve cold resistance have suggested the importance of a wide range of physiological, biochemical, cellular and molecular processes, and these processes have been associated with the regulation of gene transcription .
At present, low-temperature signal transduction has been widely studied, and the clearest pathway was the C-repeat (CRT)-binding factors (CBF) signal pathway, which is also known as the dehydration-responsive element-binding factors 1 (DREB1) signal pathway and is ABA (abscisic acid)-independent . Many transcription factors are involved in this signal pathway, including CBF1 / DREB1B, CBF2 / DREBlC, CBF3 / DREB1A, CBF4 / DREBlD, DREB1E, DREB1F. ICE1 (Inducer of CBF Expression 1), is located upstream of CBF, and together, they jointly regulate the expression of a spectrum of cold-regulated (COR) genes, through CBF binding to the cis-acting element (CRT/DRE) that contains a core conserved sequence of CCGAC [6, 7]. The CBF-COR pathway constitutes the predominant cold signaling pathway in plants, and the CBF gene is regulated positively by ICE1 (Inducer of CBF Expression 1). However, HOS1 (High Expression of Osmotically Responsive Gene 1) and MYB15 (myeloblastosis 15) negatively regulate the CBF genes, which provides a more complete understanding of the complexity of CBF-mediated cold signaling [8,9,10].
The expression patterns of cold-responsive genes were different for different plant species during exposure to cold [11,12,13]. Through the research on Arabidopsis transcriptome profiling, a total of 306 genes were identified as cold-responsive genes, with 218 genes increasing and 88 genes decreasing, while the studies on Cassava reported that 508 transcripts were identified as early cold-responsive genes, in which 319 sequences had functional descriptions [14, 15]. There have been many similar reports [16, 17]. The number of identified genes involved in cold stress response has been increasing, but the function of most genes have not been revealed. Only 12% of the cold-responsive genes were likely regulated by the CBF transcription factor; therefore, it was predicted that there was a CBF-independent pathway to respond to cold stress in plants . However, to date, there have been few studies on the CBF-independent pathway.
Genome-wide transcriptome analysis is a useful strategy for revealing the molecular mechanism of gene expression, and it can improve the efficiency of identifying the genes of interest. RNA-Seq is a high-throughput DNA sequencing approach, which generates a large amount of transcriptome data for both model and non-model species . This approach has been widely used to analyze the cold stress response of many plants. For example, the gene expression patterns were identified in Arabidopsis under drought, cold, high-salinity and ABA-treatment conditions . Comparative transcriptome analysis on two tobacco cultivars (cold-tolerant NC567 and cold-sensitive Taiyuan8) showed that the important COR genes were specifically induced during cold stress in NC567 . The transcriptome analysis of sunflower identified the candidate genes involved in response to chilling and salt stresses . P. koraiensis is a famous mixed fruit and wood forest with strong cold tolerance [21, 22]. It is the main tree species in the cold temperate zone and contains abundant cold resistance genes. It is an important material for the study of cold hardiness and acquiring cold resistance genes of coniferous tree species. Thus, it is appropriate and valuable to explore the responsive genes under sudden cold stress in Pinus koraiensis through transcriptome sequencing.
P. koraiensis is an evergreen tree belonging to Pinaceae, Pinus, which is mainly distributed in the northeastern part of China, the Korean peninsula, south of the Russian Far East and Honshu, Japan (124°38′ ~ 140°20′ E, 33°50′ ~ 52°40′ N) . P. koraiensis has a strong cold resistance, and it can surmount the extreme low temperature of − 40 °C in its natural growth state. Studies have shown that cold stress resulted in an increase or decrease in the abundance of transcripts associated with several metabolic pathways, and the expression data further suggested the involvement of both the CBF-dependent and independent pathways in the cold responses [5, 14]. In this study, seedlings of P. koraiensis that show healthy growth at room temperature suddenly underwent low temperature stress at − 20 °C with stress times of 6, 24 and 48 h, which did not experience cold acclimation. Exploring the expression pattern of genes under sudden cold stress and obtaining the differential expressed genes using RNA-seq and digital expression profiling would provide a valuable genetic resource for cold resistance genes of interest in future conifer breeding process.
RNA sequencing and de novo assembly
RNA sequencing was used to investigate the transcriptional changes of P. koraiensis under cold stress. In total, twelve cDNA libraries were constructed using RNA extracted from P. koraiensis needles, which were exposed to low temperature (− 20 °C) for 0 h, 6 h, 24 h and 48 h, respectively. The cDNA libraries were subjected to paired-end (PE) sequencing by the Illumina HiSeq2000 platform. After filtering out low-quality reads, a total of 763,995,954 clean reads were obtained. The clean reads were de novo assembled into contigs using the Trinity program . A total of 150,528 contigs consisting of 182,705,782 bp, with N50 length of 1951 bp, were obtained. Based on the paired-end sequence information, 123,445 unigenes consisting of 137,320,368 bp with N50 length of 1778 bp were obtained. The expressional levels of the unigenes were obtained by mapping the sequencing reads to the unigenes and normalized using RPKM (reads per kilobase per million mapped reads) method .
Functional annotation of unigenes
The unigene sequences were mapped to public databases using BLASTn or BLASTx with a cut-off of E-value of 10− 5. The databases used included NCBI Nucleotide sequence database (Nt), NCBI nonredundant database (Nr), the Universal protein (UniProt) database, the Clusters of Orthologous Groups of proteins (COG) database, the Protein families (Pfam) database, the evolutionary genealogy of genes: Nonsupervised Orthologous Groups (eggNOG) database, the Gene Ontology (GO) database, and Kyoto Encyclopedia of Genes and Genomes (KEGG) database. In total, 38,905 unigenes (31.52%) could be matched to a sequence in at least one of the databases mentioned above. The number of unigenes hits (E-value <1e− 5) was 31,997 (25.92%) and 28,163 (22.81%) in the Nr and UniProt database, respectively, followed by 22,644 (18.34%) in the Pfam database, 18,909 (15.32%) in the COG database, 18,239 (14.78%) in the GO database, 17,684 (14.33%) in the Nt database, 16,963 (13.74%) in the eggNOG database and 14,201 (11.50%) in the KEGG database (Table 1). The E-value distribution of annotation based on the Nr database is shown in Fig. 1a, which indicated that the E-values of unigenes (46%) ranged from 1e− 45 to 1e− 5, and a greater number of unigenes (54%) showed an E-value <1e− 45, which revealed strong homology. According to a search on the Nr database, the unigene sequences had the strongest BLASTx matches with the gene sequences from Picea sitchensis (27%), Beauveria bassiana (10%), Amborella trichopoda (6%), Ricinus communis (4%), Nelumbo nucifera (4%), Sordaria macrospora (3%), Pinus taeda (2%), Physcomitrella patens (2%), Vitis vinifera (1%) and Elaeis guineensis (1%) (Fig. 1b).
Go analysis was widely used to predict the uncharacterized unigene sequences, and 18,239 unigenes were classified as 54 functional groups. The functional groups were further divided into three categories (‘biological process’, ‘molecular function’, and ‘cellular component’) (Fig. 2) using BLAST2GO . For biological process, most of the unigenes were again classified as ‘biosynthetic process’ (6769), ‘transport’ (4853) and ‘transcription’ (4468); the percentages of total unigenes were 37.11, 26.61 and 24.50%, respectively, followed by ‘metabolic process’ (3705; 20.31%), ‘catabolic process’ (3610; 19.79%). The percentage of others were less than 10%. In the ‘molecular function’ category, the major subcategories were ATP binding (7260; 39.80%), DNA binding (5779; 31.68%) and metal-ion binding (5081; 27.86%), while ‘integral to membrane’ (7942; 43.54%), ‘plasma membrane’ (6593; 36.15%) and ‘cytoplasm’ (5427; 29.75%) were the most representative subcategories in the ‘cellular component’ category (Additional file 1).
COG database is a protein database. The protein sequences came from the genomes of bacteria, plants and animals. According to the sequence similarity, the 18,909 unigene sequences in P. koraiensis were matched to the COG database, which could be grouped into 25 categories (Additional file 2). The top 10 classes were: (R) ‘General function prediction only’ (3025; 16.00%), (E) ‘Amino acid transport and metabolism’ (2117; 11.20%), (C) ‘Energy production and conversion’ (1795; 9.49%), (S) ‘Function unknown’ (1778; 9.40%), (L) ‘Replication, recombination and repair’ (1739; 9.20%), (J) ‘Translation, ribosomal structure and biogenesis’ (1698; 8.98%), (O) ‘Posttranslational modification, protein turnover, chaperones’ (1558; 8.24%), (G) ‘Carbohydrate transport and metabolism’ (1488; 7.66%), (P) ‘Inorganic ion transport and metabolism’ (1439; 7.61%) and (K) ‘Transcription’ (1429; 7.56%). However, the smallest classes were (W) ‘Extracellular structures’ (4; 0.02%) and (Y) ‘Nuclear structure’ (1; 0.01%) (Fig. 3; Additional file 2).
Differential expressed genes (DEGs)
DEGseq  was used to identify differential expressed genes (DEGs) in cold treatment groups in comparison to control group with specified thresholds. In total, 9842, 9250 and 9697 genes were differentially expressed when exposed to cold for 6 h, 24 h and 48 h, respectively. Of the DEGs, 5444, 5241 and 5286 genes were upregulated, while 4398, 4009 and 4411 genes were downregulated at each time point. In addition to identification of DEGs under cold stress, we also identified genes differentially expressed between different time points of cold stress. The numbers of DEGs in the pairwise comparisons of the 6 h, 24 h and 48 h samples was relatively small. A total of 1516 genes were upregulated and 1286 genes were downregulated in the comparison between the 24 h and 6 h samples, 1784 genes were upregulated and 1589 genes were downregulated in the comparison between the 48 h and 6 h samples, and 1364 genes were upregulated and 1226 genes were downregulated in the comparison between the 48 h and 24 h samples, respectively (Table 2).
Gene ontology analysis of DEGs
GO enrichment was performed to investigate the function of DEGs. GO terms with corrected P-values < 0.05 were identified as significantly enriched. In total, 484, 371 and 543 GO terms were significantly enriched after cold treatment of 6 h, 24 h and 48 h, respectively. The top 30 enriched GO terms were shown in Additional file 3. We investigated the enriched GO terms at each time point and found that the enriched biological events of 6 h were similar to 24 h. However, the enriched biological events of 48 h were different from 6 h and 24 h. At time points 6 h and 24 h, the enriched biological processes include ‘response to acid chemical’, ‘response to biotic stimulus’, ‘response to stress’, ‘defense response’, and ‘single-organism developmental process’. The result indicated that plenty of DEGs were involved in the perception and transmission of cold signals. Meanwhile, the enriched molecular functions include ‘kinase activity’, ‘protein kinase activity’, ‘ADP binding’, ‘phosphotransferase activity, alcohol group as acceptor’ and ‘protein serine/threonine kinase activity’, which suggested that the relevant enzyme activities and products might change under cold stress. These results showed a complex regulatory cold stress response and indicated that the changes in the biological process might be very important in response to cold stress in P. koraiensis.
Pathway enrichment analysis of DEGs
KEGG analysis of the DEGs was showed in Additional file 4. The result showed that a total of 46 pathways were significantly affected under cold stress (Q-value < 0.05) in the 6 h sample, among which the top four abundant pathways were ‘plant-pathogen interaction’, ‘plant hormone signal transduction’, ‘starch and sucrose metabolism’ and ‘phenylpropanoid biosynthesis’, followed by ‘protein processing in endoplasmic reticulum’, ‘hippo signaling pathway’ and ‘spliceosome’. For either the 24 h sample or 48 h sample, the top four pathways were the same as those described above. The DEGs in the same pathways in all three stress periods should be of great concern. However, in the pairwise comparison among the 24 h, 6 h and 48 h samples, ‘plant-pathogen interaction’ and ‘phenylpropanoid biosynthesis’ were both included in the top four abundant pathways. Interestingly, ‘metabolism of xenobiotics by cytochrome P450’ could also be noted, because it was enriched in each comparison group.
Quantitative real-time PCR (qRT-PCR) analysis
To verify the accuracy of the RNA-Seq data under cold stress in P. koraiensis, 12 randomly selected DEGs (6 upregulated and 6 downregulated) were used for qRT-PCR analysis of transcript abundance with specific primers (Table 3). The result showed that the expression patterns of 11 of 12 unigenes detected via qRT-PCR were highly consistent with the RNA-Seq result; only c51762_g2 was specific (Fig. 4), which suggested that the high-throughput RNA-Seq data was reliable and demonstrated that the DEGs identified based on transcriptome sequencing were available.
Gene expression associated with perception and transmission of cold signals in P. koraiensis
The perception and transmission of cold signals is crucial for plants during cold stress. A total of 34 unigenes involved in the perception and transmission of cold signals were significantly differentially expressed (22 up- and 12 downregulated) under cold stress (Fig. 5). Of the 34 cold signal related genes, 31 (21 upregulated and 10 downregulated) acted as calcium ion receptors, 3 (1 upregulated and 2 downregulated) were involved in ABA synthesis and binding. The calcium ion receptor genes mainly encoded annexin, alpha-amylase, calcium-binding protein, calmodulin, calcium-dependent protein kinase, calcineurin B-like protein (CBL), CBL-interacting serine/threonine-protein kinase, CBL-interacting protein kinase and mitogen-activated protein kinase (Additional file 5).
The expression levels of four genes, including annexin D1 (c63631_g1), alpha-amylase isozyme 3C (c65620_g1), calcium-binding protein KIC (c61970_g1) and abscisic acid 8′ -hydroxylase 1 (c51736_g1), first increased significantly and then decreased significantly with the extension of stress time (Fig. 5, Additional file 5).
Transcription factors (TFs) in response to cold stress in P. koraiensis
To explore the TFs in response to cold stress in P. koraiensis, 8 TFs were identified, including 96 DEGs (41 upregulated and 55 downregulated) (Fig. 6). Among them, the number of genes for ethylene responsive factor (AP2) was the highest, with 30 genes (5 upregulated and 25 downregulated), followed by MYB with 24 genes (15 upregulated and 9 downregulated), NAM, ATAF1, ATAF2 and CUC2 (NAC) with 16 genes (4 upregulated and 12 downregulated) and zinc finger protein (ZFP) with 11 genes (6 upregulated and 5 downregulated). The number of other transcription factors was small, with 9 basic Helix-loop-helix (bHLH) genes (6 upregulated and 3 downregulated), 4 WRKY genes (3 upregulated and 1 downregulated), 1 upregulated ethylene-insensitive 3 (EIN3) gene and 1 upregulated vascular plant one any zinc-finger protein (VOZ) gene (Additional file 5).
The expression level of the same gene was different for different cold-stress periods. The expression of c58408_g1 and c66599_g1, two DREB genes first increased significantly and then decreased significantly with the extension of stress time (Fig. 6, Additional file 5).
Antioxidant enzymes and antioxidants in response to cold stress in P. koraiensis
In this study, 27 DEGs (17 upregulated and 10 downregulated) are shown in Fig.7, which were associated with the reactive oxygen species (ROS) family. These components of the ROS family were antioxidant enzymes and antioxidants including ascorbate peroxidase (APX), superoxide dismutase (SOD), catalase (CAT), glutathione peroxidase (GPX), blue copper protein, glutaredoxin and heat shock protein. The expression level of the same gene was different under different cold stress time. It was noteworthy that the expression level of some genes first increased significantly and then decreased significantly with the extension of stress time. These genes included c67548_g2 (SOD), c55044_g1 (CAT), c71938_g2 (CAT) and c11358_g1 (GPX) (Fig. 7, Additional file 5).
The general molecular regulation in response to cold stress for plants
The molecular regulation of plants against cold stress was showed in Fig. 8. When plants are subjected to cold stress, the cell membrane first sensed the cold signal with membrane rigidification and other changes, inducing second messengers (Ca2+, ABA, etc.) accumulation to transmit cold signals downstream . Calcium receptor is the most important cold signal sensor, and the calcium signaling pathway is the most important signal transduction pathway [28, 29]. Some of the downstream cold signals directly activate the expression of target genes and the others induce the transcriptional regulatory network to activate the expression of a series of cold-related genes, and thus synthesize antioxidant enzymes, condensate protective substances and osmotic regulatory substances, etc. . These substances work together to regulate the metabolic balance of substances and energy, and ultimately resist or adampt cold stress.
Perception and transmission of cold signals
Woody plants also first show cold signal recognition and transmission in response to cold stress . In addition to the very important calcium ion signal conduction pathway, there are some other conduction pathways, which are mainly related to ABA .
In this study, 31 DEGs were found to be related to calcium ion binding and were calcium ion receptors, among which 21 were up-regulated under cryogenic treatment, while only 3 DEGs were related to ABA. It was preliminarily speculated that calcium ion signaling pathway was also the main pathway of intracellular cryogenic signal transduction in P. koraiensis. There exist calcium ion receptors in plants, mainly including calmodulins (CaMs), calc-dependent protein kinases (DPKs), calcineurin B-like proteins (CBLs) and CBL-interacting protein kinases (CIPKs) [31, 32], which bind to calcium ions, thus changing the conformation of the proteins and regulating cold-related gene expression or directly transmitting cold signals to the downstream target gene [33, 34]. The expression of the above-mentioned calcium ion receptor genes in P. koraiensis was significantly changed under low temperature stress.
Annexin is a calcium-dependent phospholipid binding protein found in plants and animals. It is involved in many life activities regulated by calcium ions, such as signal transduction and calcium ion channel formation, and participates in cold-resistance reactions . Four kinds of annexin were found in cold-resistant wheat, among which, annexin P39 and P225 might be signal transducers in the cold signal transduction or regulators of cytoplasmic calcium ion concentration . 2 upregulated annexin genes (c56125_g1, c63631_g1) were involved in calcium-ion binding to alleviate cold stress in P. koraiensis, and similar results were also observed in studies on rice at low temperature . Alpha-amylase can participate in the binding of calcium ions to regulate the recognition and transmission of cold signals to alleviate cold stress. It can also catalyze starch into various soluble sugars, including dextrin, oligosaccharides, a small amount of maltose and glucose, etc. to regulate cell osmotic pressure to alleviate cold stress . Low temperature induced the expression of 3 alpha-amylase genes in P. koraiensis, and the expression first increased significantly and then decreased significantly with the extension of low temperature stress time. These results support the previous experiment by the team, in which the content of soluble sugar increased significantly at − 20 °C and a long stress time induced the soluble sugar content to increase first and then decrease .
Transcription factors in response to cold stress
In order to survive, plants form complex and efficient regulatory networks to resist and adapt to cold stress, in which transcriptional regulation plays a key role. Transcription factors regulate the expression of a series of genes and play a key role in plant abiotic adversity response networks by binding to cis-acting elements in the promoter region.
AP2/ERF is a large gene family, which is divided into four subfamilies, namely, AP2, RAV, DREB and ERF . Among them, DREB and ERF are closely related to cold stress, and DREB1/CBF (C-repeat binding factor) plays a role of molecular switch in transcriptional regulation [40, 41]. There were 4 DREB DEGs (c51648_g1 named DREB3A, c58408_g1 named DREB3B, c66599_g1 named DREB3C, c64188_g3 named DREB3D) and 12 ERF DEGS all belong to the AP2 family in this study, which could serve as candidate genes for studying the cold response of conifer species. The expression of DREB3A-C significantly increased with the decrease of the temperature, and their expression were up-regulated by 9.1, 8.1 and 5.6 times at − 20 °C for 6 h, respectively. The expression of DREB3D was down-regulated at − 20 °C. CBF gene family of Arabidopsis also included 4 genes, which were CBF1–4, respectively. When CBF overexpression in response to cold, most genes were positively regulated, but a few genes were inhibited by CBF overexpression . Whether DREB3A-D in P. koraiensis had similar function need further study.
Studies have shown that MYB15 in Arabidopsis was an R2R3 transcription factor and that overexpression of MYB15 resulted in decreased expression of the CBF gene and negative regulation of plant cold resistance . Both c45624_g1 and c65145_g1 in P. koraiensis were R2R3 MYB transcription factors, but their expression was up-regulated under cold stress, so further studies were needed. MYB96 has been reported to positively regulate plant cold resistance . Therefore, the regulation of MYB is a complex process under cold stress. NAC and bHLH TFs can be induced by low temperature, and they play important roles in plant cryogenic regulation network; the overexpression of PbeNAC1 and PtrbHLH could improve the cold resistance of Pyrus betulaefolia and Nicotiana tabacum, respectively [45, 46]. Four NAC and six bHLH DEGs were induced in P. koraiensis under cold stress .
In addition to the AP2/ERF, MYB, NAC and bHLH TFs, the DEGs encoding the remaining 4 TFs (WRKY, EIN3, ZFP and VOZ) were also identified in P. koraiensis in response to cold stress, which suggested that these TFs also might be important regulators that triggered a cascade of downstream gene expression. However, little research had been conducted on the function of these TFs, and further study of function characterization involved in cold tolerance in plants is still be needed.
Antioxidation mechanisms in response to cold stress
Hypothermia injury is mainly caused by oxidative stress caused by ROS accumulation. The antioxidant defense mechanism can remove excess ROS to protect plants from the harm for cold-tolerant plants . The antioxidant defense mechanism consists of antioxidant enzymes and antioxidants.
SOD is the first enzyme in antioxidant action. Its main function is to remove O2− and produce H2O2 at the same time. CAT can degrade excess H2O2 and other ROS by enzymatic action . APX is one of the key enzymes for H2O2 clearance, mainly existing in chloroplasts, where CAT enzyme does not exist. Therefore, APX is the key enzyme for H2O2 clearance in chloroplasts . GPX is a kind of peroxidase containing sulfhydryl, which can remove peroxides such as H2O2 in plants, so as to avoid ROS damage to plants . Therefore, SOD, CAT, APX and GPX are important antioxidant enzymes in plants. Their high expression can improve the antioxidant capacity of plants and enhance their low temperature tolerance. The lower the temperature, the higher the expression [51, 52]. There were some significantly differentially expressed SOD, CAT, APX and GPX genes in P. koraiensis under cold stress, among which the up-regulated genes played an important role. Blue copper protein, glutaredoxin and heat shock protein are all antioxidant substances with antioxidant functions [53, 54]. The expressions of c67548_g2 (SOD), c55044_g1 (CAT), c71938_g2 (CAT) and c11358_g1 (GPX) first increased significantly and then decreased significantly with the extension of cold stress time, which indicated that the longer the low-temperature stress time was, the more seriously the plant was hurt, and the higher the expression of antioxidant enzymes was. However, when the stress time reached a certain level, the plant might adapt to this stress state and the expression level decreased . In the previous report on the physiological indices of P. koraiensis seedlings, low temperature (− 20 °C) enhanced the activity of SOD and CAT, and SOD activity increased first and then decreased with the extension of stress time , which further demonstrated that low temperature stress could be resisted by increasing the expression of antioxidant enzymes genes, thereby increasing the activity of antioxidant enzymes in P. koraiensis.
In this study, a Pinus koraiensis dataset comprising 123,445 unigenes was generated by high-throughput sequencing, and the dynamic changes in gene expression were observed under cold stress for different stress times. A large number of DEGs were identified, especially the DEGs involved in the perception and transmission of cold signals, DEGs encoding TFs related to cold regulation and the DEGs removing ROS in antioxidation mechanisms. The transcriptome and digital expression profiling of P. koraiensis could facilitate the understanding of the molecular control mechanism related to cold responses and provide the basis for the molecular breeding of conifers. The DEGs of P. koraiensis without annotation in the Nr database were what we would study in the future, which might be specific to P. koraiensis.
Plant materials and cold treatments
Seeds of P. koraiensis were obtained from Kaishan village seeds orchard (129°45′ E, 42°40′ N) located on the east side of Changbai Mountain in Longjing city, northeast China. The average annual temperature, temperature in January and temperature in July are 5.3 °C, − 13.4 °C and 21.2 °C, respectively . The lowest temperature in Longjing was − 34.8 °C . The seeds were obtained in 2010 and germinated in 2011 in a greenhouse. After 4 years, 300 healthy seedlings were taken back to the laboratory growing in the greenhouse (20 °C, air humidity from 50% to 65%, 16 h/d light cycle, light intensity 150 μmol·m− 2·s− 1) for 1 month until they fully adapted to the environment. For cold treatment, 120 well-performed and similar seedlings were immediately placed at − 20 °C extreme low temperature from the ambient temperature (20 °C, control (CK)), and the needles were harvested after 0 (CK), 6, 24 or 48 h, which was achieved by the freezer Haier BC/BD-318HD (10 °C ~ − 26 °C). For each time point, the obtained needles included three biological replicates, and each replicate was a mixture from 10 independent seedlings.
Total RNA was extracted using the plant total RNA extraction kit (TaKaRa, Beijing, China) by following the detailed instructions from the manual. The total RNA extracted was detected by a biological analyzer (2100, Agilent, USA). Next, 20 μg of high-quality RNA was selected from each of the twelve RNA samples (each treatment including three biological replicates) for the construction of cDNA libraries, and the remaining high-quality RNA was used for qRT-PCR analysis.
cDNA library construction and RNA-sequencing
mRNA was enriched using magnetic beads with Oligo (dT), and fragmentation buffer was added to break the mRNA into short fragments (200–700 bp). The short mRNA fragments were used to synthesize double-chain cDNAs using random hexamer primers, buffer, dNTPs, RNase H and polymerase I. The cDNAs were purified using a PCR purification kit (QiaQuick, USA) and were washed with EB buffer for terminal repairing and poly(A) addition. Different cDNA fragments with various sizes were separated via agarose gel electrophoresis and were enriched via PCR amplification to construct cDNA libraries. The obtained 12 cDNA libraries were sequenced by DNA sequencer (HiSeq™ 2000, Illumina, USA), and the sequencing strategy was PE150.
Data filtering and de novo assembly
Data filtering was performed through the following steps for the obtained reads by sequencing: first, adaptor-containing reads were removed; second, the reads containing more than 5% ambiguous nucleotides were removed; and finally, the low-quality reads that contained more than 15% bases with Q-value ≤19 were removed, and the clean reads were obtained for de novo assembly. Trinity was the assembly software used, which linked the reads with overlapping sequences into much longer contiguous sequences; these longer sequences are termed as ‘contigs (contiguous sequences)’. 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 assembled by Trinity to obtain the unigene sequences that could not be extended at both ends.
All the nonredundant unigenes were subjected to BLASTn or BLASTx alignment (E-value < 10− 5) against nucleic acid database Nt or various public protein databases that included Nr, UniProt, KEGG and COG. The proteins with the best comparison result were selected to determine the sequence direction of the unigenes. If the comparison results were inconsistent, the sequence direction of the unigenes was determined based on the priority of Nr, UniProt, KEGG and COG. The ESTScan software was used to determine the sequence direction of the unigenes when the above databases were not suitable for alignment .
Sequence annotation and analysis of DEGs
The Blast2GO software was used to generate the GO terms based on the Nr annotation for the nonredundant unigenes , and the KEGG database was used to determine metabolic pathways of the unigenes . Twelve cDNA libraries of P. koraiensis that was treated for 0, 6, 24, and 48 h at − 20 °C were sequenced, and the raw data were analyzed via the above methods. The expression of unigenes was calculated by RPKM, which can eliminate the influence on the expression of calculated genes . The calculated gene expression level was directly used to compare the differences in gene expression between different cDNA libraries. To determine significant differences, a false discovery rate < 0.001 and an absolute value of log2FoldChange > 1 were set as thresholds. Unigenes with different RPKM values and conforming to thresholds were identified as DEGs  . The DEGs mentioned in this paper are all significant. Next, GO function analysis and KEGG Pathway analysis were performed on DEGs [58, 59].
The twelve remaining total RNA samples were reverse-transcribed with the ReverTre Ace®qPCR RT Kit (Toyobo, Osaka, Japan). The reverse transcription reaction contained 1 μg of RNA, 2 μl of 5 × RT Buffer, 0.5 μl of Primer Mix, 0.5 μl of Enzyme Mix and deionized water in a final volume of 10 μl; the reaction was conducted at 37 °C for 15 min and 98 °C for 5 min. Each of the generated cDNAs was diluted 10 times as the qRT-PCR template. qRT-PCR was performed with a DNA Engine Opticon™ 2 Real-Time System (Bio-Rad, USA), and the reaction was composed of 10 μl of 2 × SYBR Green Realtime PCR Master mix (Toyobo, Osaka, Japan), 2.5 μl of cDNA, 0.5 μl of upstream primer, 0.5 μl of downstream primer and deionized water in a final volume of 20 μl. Meanwhile, the 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 level of the selected genes was determined by the 2-ΔΔCt algorithm, and the P. koraiensis Tubulin alpha (TUBA) gene was used as an internal control . Each sample had three biological replicates, and the data were presented as the means ± standard errors (SE) (n = 3). The primer sequences of the selected genes are listed in Table 3.
Availability of data and materials
The datasets 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/PRJNA510863.
Ethylene responsive factor
C-repeat (CRT)-binding factors
Calcineurin B-like proteins
CBL-interacting protein kinases
The Clusters of Orthologous Groups of proteins database
Differential expressed genes
Calc-dependent protein kinases
Dehydration-responsive element-binding factors
The evolutionary genealogy of genes: Nonsupervised Orthologous Groups database
The Gene Ontology database
High Expression of Osmotically Responsive Gene
Inducer of CBF Expression
Kyoto Encyclopedia of Genes and Genomes database
NAM, ATAF1, ATAF2 and CUC2
NCBI nonredundant database
NCBI Nucleotide sequence database
The Protein families database
Quantitative real-time PCR
Reactive oxygen species
Reads per kilobase per million mapped reads
The Universal protein database
Vascular plant one any zinc-finger protein
Zinc finger protein
Chinnusamy V, Zhu JH, Zhu JK. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007;12:444–51.
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–6.
Wang M, Zhang X, Liu JH. Deep sequencing-based characterization of transcriptome of trifoliate orange (Poncirus trifoliata (L.) Raf.) in response to cold stress. BMC Genomics. 2015;16(1):555.
Zhu JH, Dong CH, Zhu JK. Interplay between cold-responsive gene regulation, metabolism and RNA processing during plant cold acclimation. Curr Opin Plant Biol. 2007;10(3):290–5.
Ito Y, Katsura K, Maruyama K, Taji T, Kobayashi M, Seki M, Shinozaki K, Yamaguchi-Shinozaki K. 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.
Medina J, Catalá R, Salinas J. The CBFs: three Arabidopsis transcription factors to cold acclimate. Plant Sci. 2011;180:3–11.
Nakashima K, Ito Y, Yamaguchi-Shinozaki K. Transcriptional regulatory networks in response to abiotic stresses in Arabidopsis and grasses. Plant Physiol. 2009;149:88–95.
Wu L, Chen X, Ren H, Zhang Z, Zhang H, Wang J, Wang XC, Huang R. ERF protein JERF1 that transcriptionally modulates the expression of abscisic acid biosynthesis-related gene enhances the tolerance under salinity and cold in tobacco. Planta. 2007;226(4):815–25.
Zhou MQ, Shen C, Wu LH, Tang KX, Lin J. CBF-dependent signaling pathway: a key responder to low temperature stress in plants. Crit Rev Biotechnol. 2011;31:186–92.
Qin F, Shinozaki K, Yamaguchi-Shinozaki K. Achievements and challenges in understanding plant abiotic stress responses and tolerance. Plant Cell Physiol. 2011;52:1569–82.
Tian DQ, Pan XY, Yu YM, Wang WY, Zhang F, Ge YY. De novo characterization of the Anthurium transcriptome and analysis of its digital gene expression under cold stress. BMC Genomics. 2013;14:827.
Wang XC, Zhao QY, Ma CL, Zhang ZH, Cao HL, Kong YM, et al. Global transcriptome profiles of Camellia sinensis during cold acclimation. BMC Genomics. 2013;14:415.
Wang JM, Yang Y, Liu XH, Huang J, Wang Q, Gu JH, et al. Transcriptome profiling of the cold response and signaling pathways in Lilium lancifolium. BMC Genomics. 2014;15:203.
Fowler S, Thomashow MF. Arabidopsis transcriptome profiling indicates that multiple regulatory pathways are activated during cold acclimation in addition to the CBF cold response pathway. Plant Cell. 2002;14:1675–90.
An D, Yang J, Zhang P. Transcriptome profiling of low temperature-treated cassava apical shoots showed dynamic responses of tropical plant to cold stress. BMC Genomics. 2012;13(1):64.
Matsui A, Ishida J, Morosawa T, Mochizuki Y, Kaminuma E, Endo TA, Okamoto M, Nambara E, Nakajima M, Kawashima M, Satou M. Arabidopsis transcriptome analysis under drought, cold, high-salinity and aba treatment conditions using a tiling array. Plant Cell Physiol. 2008;49(8):1135–49.
Winfield MO, Lu C, Wilson ID, et al. Plant responses to cold: transcriptome analysis of wheat. Plant Biotechnol J. 2010;8(7):749–71.
Martin JA, Wang Z. Next-generation transcriptome assembly. Nat Rev Genet. 2011;12:671–82.
Hu R, Zhu X, Xiang S, Liu Z. Comparative transcriptome analysis revealed the genotype specific cold response mechanism in tobacco. Biochem Biophys Res Commun. 2016;469(3):535–41.
Fernandez P, Rienzo JD, Fernandez L, Hopp HE, Paniego N, Heinz RA. Transcriptomic identification of candidate genes involved in sunflower responses to chilling and salt stresses based on cDNA microarray analysis. BMC Plant Biol. 2008;8(1):11–28.
Wang F, Zhang QH, Tian YG, et al. Comprehensive assessment of growth traits and wood properties in half-sib Pinus koraiensis families. Euphytica. 2018;214(11):201–16.
Wang F, Wang YX, Wang CL, et al. Variation of the growth, fruiting and resistanceto desease and insect of the half-sib families of Pinus koraiensis superior trees. Chin J Appl Ecol. 2019;30(5):1679–86.
Ma JL, Zhuang LW, Chen D, Li JW. The distribution of Pinus Koraiensis. J Northeast Forestry Univ. 1992;5:40–8.
Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.
Mortazavi A, Williams BA, McCue K, et al. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8.
Conesa A, Götz S, García-Gómez JM, et al. Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics. 2005;21(18):3674–6.
Wang L, Feng Z, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2010;26(1):136–8.
Ma Y, Dai XY, Xu YY, Luo W, Zheng XM, Zeng DL, et al. COLD1 confers chilling tolerance in rice. Cell. 2015;160(6):1209–21.
Wu FZ, Wang HX, Xu GH, et al. Research progress on the physiological and molecular mechanisms of woody plants under low temperature stress. Scientia Silvae Sinicae. 2015;51(7):116–25.
Xiong LM, Ishirani M, Lee H, Zhu JK. The Arabidopsis LOS5/ABA 3 locus encodes a molybdenum cofactor sulfurase and modulates cold stress-and osmotic stress-responsive gene expression. Plant Cell. 2001;13(9):2063–83.
Klimecka M, Muszyńska G. Structure and functions of plantcalcium-dependent protein kinases. Acta Biochim Pol. 2007;54(2):219–33.
Reddy VS, Reddy AS. Proteomics of calcium-signaling components in plants. Phytochemistry. 2004;65(12):1745–76.
Yang T, Poovaiah BW. Calcium/calmodulin-mediated signal network in plants. Trends Plant Sci. 2003;8(10):505–12.
Sathyanarayanan PV, Poovaiah BW. Decoding Ca2+ signals in plants. Crit Rev Plant Sci. 2004;23(1):1–11.
Zhang N, Shang ZL. Annexin in plant cells. Plant Physiol Commun. 2010;46(3):277–83.
Gerke V, Moss SE, Gerke V, Moss SE. Annexins: from structure to function. Physiol Rev. 2002;82(2):331–71.
Hashimoto M, Toorchi M, Matsushita K, Iwasaki Y, Komatsu S. Proteome analysis of rice root plasma membrane and detection of cold stress responsive proteins. Protein Pept Lett. 2009;16(6):685–97.
Ling GT. Food additive handbook. 3rd ed. Beijing: China light industry press; 2012. p. 681–713.
Wang F, Liang DY, Pei XN, Zhang QH, Zhang P, Zhang JQ, Lu ZM, Yang YC, Liu GF, Zhao XY. Study on the physiological indices of Pinus sibirica and Pinus koraiensis seedlings under cold stress. J Forestry Res. 2018. https://doi.org/10.1007/s11676-018-0833-0.
Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K. AP2/ERF family transcription factors in plant abiotic stress responses. Biochim Biophys Acta. 2012;1819(2):86–96.
Licausi F, Ohme-Takagi M, Perata P. APETALA2/ethylene responsive factor (AP2/ERF) transcription factors: mediators of stress responses and developmental programs. New Phytol. 2013;199(3):639–49.
Matysik J, Alia Bhalu B, Mohanty P. Molecular mechanisms of quenching of reaction oxygen species by proline under stress in plants. Curr Sci. 2002;82:525.
Agarwal M, Hao Y, Kapoor A, Dong CH, Fujii H, Zheng X, Zhu JK. 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(49):37636–45.
Guo L, Yang H, Zhang X, Yang S. Lipid transfer protein 3 as a target of MYB96 mediates freezing and drought stress in Arabidopsis. J Exp Bot. 2013;64(6):1755–67.
Jin C, Li KQ, Xu XY, Zhang HP, Chen HC, Chen YH, Hao J, Wang Y, Huang XS, Zhang SL. A novel NAC transcription factor, PbeNAC1, of Pyrus betulifolia confers cold and drought tolerance via interacting with PbeDREBs and activating the expression of stress-responsive genes. Front Plant Sci. 2017;8:1049.
Huang XS, Liu JH. A basic helix-loop-helix transcription factor PtrbHLH of Poncirus trifoliata confers cold tolerance and modulates POD-mediated scavenging of H2O2. Plant Physiol. 2013;162(2):1178–94.
Gill SS, Tuteja N. Polyamines and abiotic stress tolerance in plants. Plant Signal Behav. 2010;5(1):26–33.
Garratt LC, Janagoudar BS, Lowe KC, et al. Salinity tolerance and antioxidant status in cotton cultures. Free Radic Biol Med. 2002;33(4):502–11.
Sun WH. Cloning and functional analysis of LetAPX under environment stress. Thesis for Doctor, Shandong Agricultural University, Supervisor: Meng QW. 2008; p.26–40.
Qiao XR, Zhang JY. Research progress on GPX in plants. Biotechnol Bull. 2016;32(9):7–13.
Moses JC. Cu/Zn superoxide dismutase activity and respective gene expression during cold acclimation and freezing stress in barley cultivars. Biol Plant. 2012;56(4):693–8.
Milla MAR, Maurer A, Huete AR, Gustafson JP. Glutathione peroxidase genes in Arabidopsis are ubiquitous and regulated by abiotic stresses through diverse signaling pathways. Plant J. 2010;36(5):602–15.
Holmgren A. Antioxidant function of thioredoxin and glutaredoxin systems. Antioxid Redox Signal. 2000;2(4):811–20.
Ye XY, Ma FW. Advances in the research of heat shock proteins in plants. Acta Agriculturae Boreali-occidentalis Sinica. 2004;13(2):109–14.
Liu YJ, Cao HX, Zhang RL. Effect of low temperature stress on physiological changes in oil palm (Elaesis guineensis Jacg.) under different time. Bull Bot Res. 2014;25(3):549–56.
Liang DY, Jin YZ, Zhao GJ, Dong YH, Leng WW, Chen CL, Wang H, Zhao XY. The study on the variation of the growth and wood characters of 50 Pinus Koraiensis clones. J Beijing Forestry Univ. 2016;38(6):51–9.
Iseli C, Jongeneel CV, Bucher P. ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. ISMB. 1999;99:138–47.
Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.
Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25–9.
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–8.
We acknowledge all the members of the research group.
Study design, sequencing, data collection, data analysis and data interpretation were supported by the Forestry Science and Technology Development Project (KJZXSA2019046). The labour service, travel expense et al. were supported by the Innovation Project of State Key Laboratory of Tree Genetics and Breeding (2018FT6). Both funding bodies of this study did not have any role in writing of the manuscript.
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Go analysis of the unigenes for P. koraiensis in response to cold stress.
COG annotation of the unigenes for P. koraiensis in response to cold stress.
Gene Ontology clustering of DEGs.
KEGG analysis of the DEGs.
All the DEGs in each comparison.
About this article
Cite this article
Wang, F., Chen, S., Liang, D. et al. Transcriptomic analyses of Pinus koraiensis under different cold stresses. BMC Genomics 21, 10 (2020). https://doi.org/10.1186/s12864-019-6401-y
- Pinus koraiensis
- Cold stress
- Differential expressed genes