RNA-seq based transcriptomic analysis uncovers α-linolenic acid and jasmonic acid biosynthesis pathways respond to cold acclimation in Camellia japonica

Camellia is a well-known ornamental flower native to Southeast of Asia, including regions such as Japan, Korea and South China. However, most species in the genus Camellia are cold sensitive. To elucidate the cold stress responses in camellia plants, we carried out deep transcriptome sequencing of ‘Jiangxue’, a cold-tolerant cultivar of Camellia japonica, and approximately 1,006 million clean reads were generated using Illumina sequencing technology. The assembly of the clean reads produced 367,620 transcripts, including 207,592 unigenes. Overall, 28,038 differentially expressed genes were identified during cold acclimation. Detailed elucidation of responses of transcription factors, protein kinases and plant hormone signalling-related genes described the interplay of signal that allowed the plant to fine-tune cold stress responses. On the basis of global gene regulation of unsaturated fatty acid biosynthesis- and jasmonic acid biosynthesis-related genes, unsaturated fatty acid biosynthesis and jasmonic acid biosynthesis pathways were deduced to be involved in the low temperature responses in C. japonica. These results were supported by the determination of the fatty acid composition and jasmonic acid content. Our results provide insights into the genetic and molecular basis of the responses to cold acclimation in camellia plants.

. MDA and relative electrical conductivity change in different C. japonica cultivars that subjected to natural low temperature in 2013 winter. LT0 and LT40 represent plant that subjected to low temperature for 0 day and 40 days, respectively. The bars represent the standard error (n = 3). * and ** indicate a significant difference between LT0 and LT40 at P < 0.05 and P < 0.01 levels, respectively (two tailed T-test).
Scientific RepoRts | 6:36463 | DOI: 10.1038/srep36463 were lower than those in the other cultivars (Fig. 1). These results suggested that the C. japonica cultivar 'Jiangxue' is more tolerant to cold stress than other cultivars. To elucidate how 'Jiangxue' addresses cold stress, we used 'Jiangxue' leaves during cold acclimation to obtain the transcriptomic response data. RNA sequencing and de novo transcriptome assembly. Plants in genus Camellia usually have large genome sizes, and thus far, the draft genome of C. japonica has not been sequenced and assembled. Therefore, to gain insight into the molecular mechanisms of cold tolerance in C. japonica, we used RNA sequencing in this study.
RNA was isolated from the leaves of 'Jiangxue' at 2, 8, 24, 72 and 168 h of 4 °C cold acclimation and control treatments. Three biological replicates were sampled. In total, 18 RNA samples were subjected to paired-end RNA sequencing. A total of 1,031 million raw reads were generated. After the low-quality reads were removed, and the adapter sequences were trimmed, 1,006 million (125.84 Gb) clean reads were obtained with an average of 55.9 million reads (6.99 Gb) for each sample (Supplementary Table 2).
The Trinity sequence assembler was used, and 367,620 transcript sequences and 207,592 unigenes were generated with transcript lengths ranging from 201 bp to 14,487 bp. The average length of an assembled transcript was 860 bases, and the N50 length was 1,415 bases (Supplementary Table 3). The length statistics of the assembled transcripts and unigenes are shown in Supplementary Fig. 1. Functional annotation of the assembled transcriptome. To predict and analyse the function of the assembled unigenes, we assessed the non-redundant sequences using a BLASTX search against the following databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (protein family database), Swiss-Prot (a manually annotated and reviewed protein sequence database), GO (Gene Ontology), COG (cluster of orthologous groups) and KEGG (Kyoto Encyclopaedia of Genes and Genomes). After the analysis, 65,150 (31.38%), 44,101 (21.24%), 47,011 (22.64%), 48,079 (23.16%) and 48,916 (23.56%) unigenes returned BLAST results and showed identity with sequences in the Nr, Nt, Pfam, SwissProt and GO databases, respectively. Overall, 84,910 (40.9%) unigenes were significantly matched to known genes in the public databases mentioned above (Supplementary Table 4).
To determine the potential functions of unigenes, we used GO assignments to classify the predicted C. japonica genes, and 48,916 unigenes were assigned to three major functional categories (Biological Process, Cellular Component and Molecular Function) and 56 subcategories (Fig. 2a). In terms of Biological Processes, 'cellular processes' and 'metabolic processes' were the top two GO terms, which indicated that the plants were undergoing rapid cell growth and were metabolically active. In the Molecular Function category, the unigenes were predominantly assigned to the 'binding' and 'catalytic activities' groups. In the binding subset, 'organic cyclic compound binding' and 'heterocyclic compound binding' were the most common groups. In the 'catalytic activity' subset, the two major groups were 'transferase activity' and 'hydrolase activity' . In the Cellular Component category, the unigenes were frequently assigned to 'cell' , 'cell junction' , 'macromolecular complex' and 'organelle' .
To identify the biological pathways in the annotated C. japonica sequences, we annotated the unigenes to the reference pathways in the KEGG using KeggArray software, and 23,476 unigenes were assigned to five specific pathways, including 'Cellular Processes' , 'Environmental Information Processing' , 'Genetic Information Processing' , 'Metabolism' , and 'Organism Systems' (Fig. 2b). Among these pathways, the 'translation' cluster represented the largest group, followed by 'carbohydrate metabolism' and 'signal transduction' .
To classify the orthologous gene products, 23,514 unigenes were subdivided into 26 COG classifications (Fig. 2c). Among these classifications, the cluster of 'general function prediction only' represented the largest group, followed by 'posttranslational modification, protein turnover, chaperones' and 'translation, ribosomal structure and biogenesis' . The two categories involving 'cell motility' and 'unnamed protein' represented the smallest COG classifications.
Identification of genes involved in the cold response. The clean data from each sample were mapped onto the assembled transcriptome, and gene expression levels were estimated using RSEM (RNA-Seq by Expectation Maximization) for each sample. Genes with FPKM (fragments per kilobase of transcript per million fragments mapped) values equal to or larger than 0.3 were defined as expressed.
To investigate genes involved in the cold response of C. japonica during cold acclimation, we identified the differentially expressed genes (DEGs) between cold-acclimated and non-acclimated samples using the DESeq R package, with an adjusted q value < 0. 005 (Fig. 3a).
The number of genes up-or down-regulated at different sampling time points during cold acclimation is shown in Fig. 3b. The number of DEGs increased rapidly at 24 h of cold acclimation and reached to 19,908 at 72 h. In addition, there were more up-regulated genes than down-regulated genes during this period. The number of DEGs decreased at 168 h of cold acclimation, and the number of up-regulated genes was nearly equal to the number of down-regulated (Fig. 3b).
DEGs were divided into three groups based on their expression patterns by Genesis based on the K-means method ( Supplementary Fig. 2). Type I genes, which included 13,663 genes, were positively modulated during cold acclimation and were divided into ten sub-clusters. Type II genes, which included 11,335 genes, were negatively affected during this process and were divided into nine sub-clusters. Type III represents genes showing complex expression patterns and was divided into four sub-clusters, containing 2710 genes.
To validate the expression profiles obtained by RNA-seq, we used qRT-PCR to confirm the expression of 35 genes that showed different levels during cold acclimation. The relative expression levels measured by qRT-PCR were converted to fold changes (cold-acclimated/non-acclimated) to enable a direct comparison with the RNA-seq data. The Pearson correlation coefficient was calculated by SPSS to assess the correlation between different platforms. Overall, the qRT-PCR measurements were highly correlated with the RNA-seq results  Global gene regulation during cold acclimation. To survive in low temperatures, plants require an effective response to adapt to the low temperature. A rapid reaction may predominantly rely on a gene network involving TFs, protein kinases, plant hormone signal transduction pathways, unsaturated fatty acid biosynthesis, and JA biosynthesis pathways.

TFs.
A total of 1455 DEGs encoding putative TFs were identified in C. japonica during cold acclimation. These TFs were classified into 53 different common families based on the classification of their Arabidopsis homologues (Supplementary Table 5). Among these families, the WRKY family was the largest group (294, 20%), followed by the NAC family (239, 16%), the bHLH family (96, 7%), the ethylene-responsive TF (ERF) family (72, 5%) and C2H2 (65, 4%) ( Supplementary Fig. 3). In the WRKY and bHLH TF families, the number of up-regulated genes was nearly equal to that of the down-regulated genes. In the ERF and C2H2 families, the number of up-regulated genes was substantially greater than that of the down-regulated genes; while there were more down-regulated genes than up-regulated genes in the NAC family ( Supplementary Fig. 4).
Protein kinases. Protein phosphorylation plays a pivotal role in plant responses to cold stress 25 . In our study, 74 MAP kinase kinase kinases (MPKKK), 5 MAP kinase kinases (MPKK), 11 MAP kinases (MPK), 17 calcium dependent protein kinases (CPK) and 31 CBL-interacting protein kinases (CIPK) were identified in C. japonica during cold acclimation based on MPK, MPKK, MPKKK, CPK and CIPK gene family information in Arabidopsis 26,27 using the BLASTP programme (Supplementary Table. 6). In the MPKKK family, the number of up-regulated genes was slightly more than that of the down-regulated genes. Interestingly, most of DEGs in the CIPK family, the CPK family, the MPK family and the MPKK family were induced by low temperature in C. japonica ( Supplementary Fig. 5). These results indicated that protein phosphorylation may be induced by low temperature during cold acclimation in C. japonica.
To confirm the results of the transcriptomic analysis, we verified the expression of several C18 unsaturated fatty acid biosynthesis genes using qRT-PCR. The expressions of CjFabG1, CjSAD1, CjSAD2 and CjFAD8 in 'Jiangxue' were induced during cold acclimation, especially after 24 h, while the CjFAD2 transcripts were transiently inhibited in the early period (2 and 8 h) of cold acclimation and were then induced after 24 h (Fig. 6a). In addition, the expression of these genes in two other C. japonica cultivars 'Desire' and 'Nuccio's Bella Rossa' , which were more sensitive to cold stress compared with 'Jiangxue' , were also determined. The expression of CjFabG1, CjSAD1, CjSAD2 and CjFAD2 were slightly inhibited at the early period (2 and 8 h) and then induced after 24 h, while the CjFAD8 transcripts were almost continuous induced. However, the induction degrees of these genes by low temperature in 'Desire' and 'Nuccio's Bella Rossa' were lower than that of in 'Jiangxue' (Fig. 6a). These results indicated that the unsaturated fatty acid biosynthesis pathway was activated by cold acclimation in C. japonica.
To further validate the above expression data, we verified the transcript levels of six JA biosynthesis genes using qRT-PCR. As shown in Fig. 6b, the expression of CjLOX3, CjAOC2, CjAOC3, CjAOS1, CjOPR3 and CjOPCL1 was strongly up-regulated during cold acclimation in 'Jiangxue' leaves, consistent with the results of expression profiling. However, the expression of these genes in 'Desire' and 'Nuccio's Bella Rossa' during cold acclimation just slightly up-regulated (Fig. 6b). Together with the activated JA signalling transduction pathway ( Supplementary Fig. 6a), the whole JA signalling pathway may play key roles in cold acclimation in C. japonica. Both α-linolenic acid and JA content were enhanced in C. japonica leaves during cold acclimation.
To confirm that the unsaturated fatty acid biosynthesis pathway is important for the cold response in C. japonica, we measured the leaf fatty acid compositions in different C. japonica cultivars during cold acclimation using gas chromatography. As shown in Table 1, in 'Jiangxue' , the content of palmitic acid (16:0), stearic acid (18:0) and oleic acid (18:1) gradually decreased during cold acclimation, while the content of α -linolenic acid (18:3) increased significantly during the same process, consistent with the up-regulated expression of unsaturated fatty acid biosynthesis genes. Similar tendencies were found in 'Desire' and 'Nuccio's Bella Rossa' , but the changes were not as significantly as that of in 'Jiangxue' . Moreover, the degree of fatty acid unsaturation in all C. japonica cultivars gradually increased during cold acclimation, while both the degree of fatty acid unsaturation and its increment in 'Jiangxue' were remarkable higher than that of 'Desire' and 'Nuccio's Bella Rossa' (Fig. 7a), which was consistent with the results from the gene expression analysis and the determination of fatty acid composition. These results indicated that the unsaturated fatty acid biosynthesis pathway, particularly the biosynthesis of α -linolenic acid (18:3), was involved in the cold response and was activated by cold acclimation to enhance cold tolerance in C. japonicas. c (18:3) is the major substrate of JA biosynthesis 10 , and JA-related pathways, including the JA biosynthesis and JA signal transduction pathways, were activated during cold acclimation in this study. Therefore, we determined whether JA content was also increased during the cold acclimation process. The free JA levels in different C. japonica cultivars during cold acclimation were measured. In 'Jiangxue' , the free JA content increased sharply after 2 h of 4 °C acclimation and was maintained at a high level from 2 h to 24 h compared with non-acclimated leaves. Then, at 72 h cold acclimation, the free JA contents decreased to the same level as that of the non-acclimated  (Fig. 7b). Similar changes in JA content were found in 'Desire' and 'Nuccio's Bella Rossa' , while the increments of JA content in 'Desire' and 'Nuccio's Bella Rossa' were lower than that of in 'Jiangxue' . These results suggested that the activation of JA biosynthesis pathway at the early stage of cold acclimation may contribute to the cold tolerant in C. japonica.

Discussion
Cold tolerance is an important trait in plants because it limits the geographical distribution of wild species and the growth performance and yield of crop plants. Unfortunately, most species in the genus Camellia are difficult to overwinter in high altitude cold regions because of their weak freezing tolerance 21 . Exploring the genetic resources and understanding the molecular and physiological mechanisms of chilling-resistant Camellia are needed to develop an effective strategy for genetic breeding of cold tolerance in camellia plants. In this study, we performed an in-depth transcriptomic analysis of the responses to low temperature in C. japonica. Our results revealed that complex gene networks, especially the unsaturated fatty acid and JA biosynthesis pathways, are involved in the cold acclimation process in C. japonica, further expanding our knowledge of how camellia plants increase cold-resistance after cold acclimation.
With the application of next-generation sequencing technology, plant genomic studies have been rapidly performed in recent years 29 . However, for species without a sequenced genome, such as the non-model plant species C. japonica, transcriptome sequencing is a rapid approach to study the regulatory pathways and molecular mechanisms of various complex physiological processes. Although several transcriptomic analysis in tea plants  (Camellia sinensis) have been reported 23,30 , transcriptomic studies in C. japonica are still rare. In this study, we generated approximately 126 Gb of paired-end reads and obtained a reference transcriptome consisting of 367,620 contigs with a mean size of 860 bp and a N50 of 1415 bp (Supplementary Table 2 30 , the present study obtained more complete annotation information. TFs always act as master switches by controlling the expression of a series of genes to regulate different plant developmental processes or responses to biotic or abiotic stresses 31 . The WRKY family is one of the best-characterized classes of plant TFs, and many WRKY TFs have been shown to be involved in the cold response 32 . In this study, the WRKY family was the largest group involved in the cold acclimation process ( Supplementary Figs 2 and 3; Supplementary Table 5), which indicated that increasing the tolerance to cold stress in C. japonica might require the complex mechanisms of signalling and transcriptional reprogramming controlled by WRKY proteins. The ERF family TFs are generally considered to be mediators of ethylene-related responses 33 . Genome-wide expression analyses of the AP2/ERF family genes indicated that many ERF family genes are also induced by low or high temperature 34 . In the present study, 72 DEGs encoding ERF family proteins were identified during cold acclimation and most of these genes were up-regulated (Supplementary Fig. 3; Supplementary Table 5), indicating that the induction of ERF family genes was required during the cold acclimation process in C. japonica.
MPK cascades convey stress signals from receptors to specific effectors to regulate gene expression, cell activities, and protein functions in various developmental and adaptive processes 35,36 . In our study, 75 MPKKK, 11 MPK and 5 MPKK genes were found to be involved in signal transduction upon low temperature stress (Supplementary Fig. 4; Supplementary Table 6), indicating that MAPK cascades play an important role in the cold acclimation process in C. japonica. CPK and CIPK, two protein kinase families known as Ca 2+ sensors, have been shown to act as positive regulators of the plant cold response 27,37,38 . In our study, 31 CIPK and 17 CPK genes were differentially expressed during cold acclimation, and most of these genes were up-regulated (Supplementary Fig. 4; Supplementary Table 6). These results suggest that CPKs and CIPKs promote cold tolerance in C. japonica.
Phytohormones, such as ABA, JA, ET, SA, auxin, GA, CK, and BR, appeared to play critical roles in the complex signalling cascades and were essential for the plant adapation to biotic stresses by mediating a wide range of adaptive responses 39 . The ABA level increased in response to low temperatures in Arabidopsis 40 , and transcriptome analysis revealed that approximately 10% of ABA-responsive genes are also responsive to cold stress 41 . In the ABA signalling pathway, the ABA receptor, PP2C and SnRK2 form a signalling complex referred to as the ' ABA signalosome' , which can turn on ABA signals by phosphorylation of downstream factors, such as the AREB/ ABF TFs or membrane proteins involving ion channels 42 . As shown in this study, 24 genes in the ABA signalling pathway, including ABA receptor genes, SNRK2 genes, PP2C genes and ABF genes (Supplementary Fig. 5a; Supplementary Table 7), were differentially expressed during cold acclimation. These results thus indicate that the ABA-dependent cold signalling pathway may play important roles in C. japonica. JA and its derivatives are essential signalling molecules in plant developmental regulation and environmental stress responses 10 . A previous study showed that MaMYC2s, a key regulator of JA signalling, is involved in MeJA-induced chilling tolerance in banana fruit by interacting and functionally coordinating with MaICE1 13 . In the present study, four MYC genes, including two MYC2 genes and two MYC4 genes, were significantly induced during the cold acclimation process (Supplementary Fig. 5a; Supplementary Table 7), which indicates that the same regulatory pathways may exist in C. japonica. Furthermore, most DEGs in the JA signalling pathway were up-regulated in our study, indicating that JA may play a positive role in C. japonica adaptation to low temperatures.
Thus far, no specific receptor has been identified in plants in response to cold stress. However, the plasma membrane itself could be the primary sensor of temperature fluctuations 36,43 . Studies of the Arabidopsis fad2 mutant have shown that low temperatures can be sensed by plant cells via membrane rigidification 3,8 . Other FAD genes were shown to be important for the membrane stabilization of plants in low temperature environments by maintaining high levels of polyunsaturated fatty acids 6,44,45 . In the present study, the expression of genes involved in several steps of the C18 unsaturated fatty acid biosynthesis pathway was up-regulated during cold acclimation (Figs 5a and 6a); as a result, α -linolenic acid (18:3), the end-product of this pathway, accumulated (Table 1). These results thus indicated that the C18 unsaturated fatty acid biosynthesis pathway was globally actived during cold acclimation in C. japonica. Moreover, a previous study showed that overexpression of the ω -3 FAD genes BnFAD3 and StFAD7, which catalyse the conversion of linoleic acid (18:2) to linolenic acid (18:3) in tomato, increased the linolenic acid (18:3) level, while the level of linoleic acid (18:2) declined 46 . Consistent with the results from that study, our results showed that the precursors of linoleic acid (18:2) and α -linolenic acid (18:3), including palmitic acid (16:0), stearic acid (18:0) and oleic acid (18:1), were decreased in C. japonica leaves during cold acclimation (Table 1), while the degree of fatty acid unsaturation increased (Fig. 7a). Cold stress can change the membrane fatty acid compositions in many plants 5,24 . Therefore, the regulation of C18 fatty acid desaturation likely plays an important role in membrane stabilization under cold stress in C. japonica.
The substrate of JA biosynthesis is α -linolenic acid (18:3) 10,11 . A study of tomato showed that loss of function of a chloroplast FAD gene, LeFAD7, reduced the α -linolenic acid (18:3) content and impaired JA accumulation in defence responses 47 , suggesting that α -linolenic acid (18:3) in the chloroplast pools is required for stress-stimulated JA synthesis. Free α -linolenic acid (18:3) is oxygenated by lipoxygenase enzymes (LOX) and Scientific RepoRts | 6:36463 | DOI: 10.1038/srep36463 then converted to 12-oxo-phytodienoic acid (OPDA) through the combined action of allene oxide synthase (AOS) and allene oxide cyclase (AOC). OPDA is subsequently converted to JA by reduction and three cycles of β -oxidation 10,48 . In our study, the expression of nine JA biosynthesis genes, LOX3, AOS1, AOS2, AOC2, AOC3, OPR3, OPCL1, ACX2, and ACAA1, were up-regulated during cold acclimation (Figs 5b and 6b). In addition, our data further confirmed that the free JA content was increased at the early stage of cold acclimation in C. japonica leaves (Fig. 7b). Previous studies have shown that JA and its derivatives are important regulators of the plant responses to abiotic stresses 10 . The induced expression of JA biosynthesis genes and the increased free JA content in our study indicated that JA and its signalling pathway may play critical roles in the early stage of cold acclimation in C. japonica. Interestingly, a previous study showed that the ratio of unsaturated/saturated fatty acid in MeJA-treated loquat fruit was significantly higher than that in control fruit, along with reduced injury symptoms in MeJA-treated fruits 49 , indicating that JAs may also influence unsaturated fatty acid biosynthesis. In our study, both α -linolenic acid (18:3) and JA biosynthesis were enhanced during the cold acclimation process. However, whether α -linolenic acid (18:3) or JA acts as the upstream regulatory factor during the cold acclimation process in C. japonica will need to be determine in future studies. Determination of MDA and relative electric conductivity. The MDA content analysis was performed using the thiobarbituric acid method. First, 0.5 g fresh leaf samples were homogenized in 5 ml extraction buffer (10% trichloroacetic acid, w/v). After centrifugation at 10,000× g for 15 min, a 2 ml aliquot of the supernatant was heated with 2.0 ml of 0.5% (w/v) thiobarbituric acid at 100 °C for 30 min. Absorbance of the supernatants was measured at 450, 532, and 600 nm. The MDA content was calculated as described by Draper et al. 50 .

Methods
The relative electric conductivity was determined following the methods as described by Wang et al. 23 . Briefly, leaves were washed with deionized water. Then, 0.5 g midvein-excluded leaf samples were placed in closed vials containing 20 ml of deionized water and incubated at 25 °C on a rotary shaker for 24 h. The electrical conductivity of the solution (L1) was determined. Samples were then heated to 100 °C for 20 min and the final electrical conductivity (L2) was determined after equilibration at 25 °C. The relative electrical conductivity was defined as follows: relative electrical conductivity (%) = (L1/L2 × 100%).
Sampling for RNA-seq and RNA preparation. Leaves  Library preparation and RNA sequencing. RNA samples were sent to Beijing Novogene Bioinformatics Technology Co., Ltd., where the libraries were produced and sequenced. Briefly, the RNA-seq libraries were generated using the NEBNext ® Ultra ™ RNA Library Prep Kit for Illumina ® (NEB, USA) following manufacturer's recommendations. mRNA was purified from total RNA using poly-T oligo-linked magnetic beads. Then, purification fragmentation buffer was added to cleave the mRNA molecules into short fragments. First-strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase (RNase H-). Second-strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. NEBNext Adaptors with hairpin loop structures were ligated for hybridization. To select cDNA fragments of preferentially 150-200 bp in length, we purified the library fragments using the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 μ l USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C prior to PCR. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, universal PCR primers and index (X) primer. Finally, the PCR products were purified (AMPure XP system), and library quality was assessed on an Agilent Bioanalyzer 2100 system. The resulting cDNA library was sequenced using an Illumina HiSeq platform. The original data set was deposited in the NCBI Sequence Read Archive (accession no. SRP076436).
De novo transcript assembly and functional annotation of unigenes. Raw reads generated from sequencing were firstly processed through in-house Perl scripts. In this step, clean reads were obtained by removing reads containing adapter sequences, reads containing poly-N sequences and low quality reads from the raw data. First, bases with Phred quality scores lower than 20 were trimmed from the 3′ end of the sequence until reaching a base with a higher quality score (≥ 20). If the read length was shorter than 50 bp, it was discarded. Second, reads were further filtered by the criterion that 70% of the bases in one read must have high quality scores (≥ 20). Third, only paired-end reads were used for further assembly. De novo transcript assembly was conducted using Trinity software 51 , which consisted of three successive software modules: Inchworm, Chrysalis, and Butterfly.
The unigene sequences of the four tea plant cultivars were searched using BLASTX against the Nt, Nr, Pfam, Swiss-Prot, GO, COG, and KEGG databases (E-value ≤ 1E-5) to retrieve protein functional annotatio ns based on sequence similarity.
Gene expression quantification and differential expression analysis. To quantify transcript abundance, the sequenced pair-end reads were mapped onto the assembled transcriptome, and the read count for each gene was obtained from the mapping results. Mapped reads were used for quantification by RSEM software 52 . Gene or isoform abundance was represented by the fragment per kilobase of transcript per million fragment mapped (FPKM) value, and those transcripts with FPKM values equal to or larger than 0.3 were defined as expressed. Prior to differential gene expression analysis, for each sequenced library, the read counts were adjusted by the edgeR programme package through one scaling normalized factor. Differential expression analysis of two treatments was performed using the DEGseq R package 53 . Three independent biological replicates for each treatment were analysed. The P value was adjusted using a q value (Storey et al. 2003). Q value < 0.005 was set as the threshold for significantly different expression.
Validation of RNA-seq analysis by qPCR. qPCR assays were performed to confirm the RNA-seq results.
C. japonica leaves acclimated at 4 °C for 0, 2, 8, 24, 72 and 168 h were sampled and used for RNA extracting. 5 μ g of total RNA was treated with 10 U of DNase I (Thermo Scientific, USA) to remove residual DNA and then used for reverse transcription with TransScript First-Strand cDNA Synthesis Super Mix (TransGen). qRT-PCR was performed as described previously 54 . Three independent biological replicates for each sample and three technical replicates for each biological replicate were analyzed. The C. Japonica 18S RNA gene was used as the internal control. All of the primers that were used are listed in Supplementary Table 8.
Fatty acid analysis. C. japonica leaves acclimated at 4 °C for 0, 2, 8, 24, 72 and 168 h were sampled and used for lipid extraction and profiling. The fatty acid extraction was performed following a described protocol 55 . Briefly, samples of ~70 mg fresh weight were ground and then heated at 90 °C in 2.5% (v/v) H 2 SO 4 in methanol for 90 min in screw-capped tubes. After the addition of 500 μ l of hexane containing 0.01% butylated hydroxytoluene, fatty acids were extracted into the organic phase by shaking and the tubes were centrifuged at low speed. The samples were then analyzed using an Agilent 7890 series gas chromatograph. Samples (5 μ l) of the organic phase were separated by gas chromatography on a 30 m × 0.25 mm, 0.25 μ m film Agilent J&W GC column (USA) and quantified using a flame ionization detector (FID). The gas chromatograph was programmed for an initial temperature of 180 °C for 1 min followed by an increase of 10 °C/min to 220 °C; this final temperature was maintained for a further 4 min. Glyceryl triheptadecanoate (Sigma-Aldrich: T2151) was used as the internal standard. JA determination. Sample preparation and JA content quantitation were performed as described previously 56 . Leaves of the C. japonica acclimated at 4 °C for 0, 2, 8, 24, 72 and 168 h were sampled, and samples of ~100 mg fresh weight were frozen with liquid nitrogen and stored at − 80 °C before analysis. The tissues were ground in 800 μ l of pre-cooled extraction buffer (methanol:ddH 2 O:acetic acid, 80:19:1, v/v/v) supplied with 9,10-dihydro-JA (DHJA; Olchemin: 014 5321) as internal standard with ceramic beads (3 mm) using a Tissue Lyser (JXFSTPRP-192; China) for 90 s; the mixture was incubated at 4 °C for 12 h with shaking in the dark. After centrifugation (10,000× g) at 4 °C for 10 min, the supernatant was transferred to a fresh tube; the pellet was mixed with another 400 μ l of pre-cooled extraction buffer without internal standard, shaken for 2 h at 4 °C in the dark, and centrifuged. After drying the combined supernatants with nitrogen gas, the pellet was resuspended in 300 μ l of pre-cooled 10% methanol and filtered using a syringe-facilitated 13 mm diameter nylon filter with a pore size 0.22 μ m. The prepared samples were detected using a 4000Q-TRAP HPLC-MS system (AB SCIEX) and the parameter settings were followed as described previously 57 . Pure JA (Sigma-Aldrich: 14631) was used as external standard to generate the standard curve to quantify the JA content.

Statistics.
Each graphical plot represents the results of multiple independent experiments (n ≥ 3), and the values are means ± SE. Statistical significance of physiological parameters (MDA and relative conductivity) was determined using two-tailed unpaired Student's t-tests, and P values pf < 0.05 were considered statistically significant. The fatty acid compositions in different sampling points were compared using a one-way ANOVA and shortest significant range (P < 0.05) post hoc analysis, and values not sharing a common letter were considered statistically significant.