De novo transcriptome assembly and comparative transcriptomic analysis provide molecular insights into low temperature stress response of Canarium album

A de novo transcriptome analysis was performed in C. album, a temperature sensitive fruit tree in China, after treatment with varied temperatures. A total number of 168,385 transcripts were assembled, comprising of 109,439 unigenes, of which 70,530 were successfully annotated. Compared with control check group (CK), which was treated under 25 °C, the chilling stress (4 °C) treated group (CT), showed about 2810 up-regulated and 2567 down-regulated genes. Whereas, group treated under freezing (− 3 °C) stress (FT) showed an up-regulation and a down-regulation of 1748 and 1459 genes, respectively. GO classification analysis revealed that DEGs related to metabolic processes, single-organism metabolic process, and catalytic activity are significantly enriched in both CT and FT conditions. KEGG pathway enrichment analysis for both CT and FT treatments showed an enrichment of genes encoding or related to glycine/serine and threonine metabolism, alpha-linolenic acid metabolism, carotenoid biosynthesis, photosynthesis-antenna proteins, and circadian rhythm. However, genes related to photosynthesis, carbon fixation in photosynthetic organisms, glutathione metabolism, pyruvate metabolism, nicotinate and nicotinamide metabolism were specifically enriched in CT condition. Nevertheless, FT treatment induced genes related to plant-pathogen interaction, linoleic acid metabolism, plant hormone signal transduction and pentose phosphate pathway. Many of the genes involved in plant hormone signal transduction showed significantly different expression in both FT and CT conditions. However, the change was more evident in FT. Here we present the first of the reports for a de novo transcriptomic analysis in C. album, suggesting that the plant shows differential responses in chilling and freezing temperatures, where the hormone signaling and transduction contribute greatly to FT responses. Our study thus paves way for future research regarding functions of these potentially identified genes.

Unigene annotation. The unigenes of C. album were annotated by using Nr, Nt, Pfam, KOG, Swiss-Prot, KEGG and GO databases ( Table 2). Results showed that 70,530 unigenes, that accounted for almost 64.44% of the total number of unigenes, were annotated in at least one database. The Nr database gave the largest number of unigene annotations (65,274; 59.64%), while those in KEGG database got the least annotations (25,196; 23.02%). In addition, 11,284 (10.31%) unigenes were annotated in all the databases. The largest distribution of annotation similarity ranged from 60 to 80% and 80-95%, accounting for 44.60% and 42.40%, respectively (Fig. 1a). According to RNA-Seq alignment results, Citrus sinensis, C. clementina and Theobroma cacao were the ones having the highest similarity to C. album (Fig. 1b). KEGG classification of the C. album unigenes. The unigenes identified in C. album were annotated in 130 KEGG metabolism pathways that belonged to 5 main categories, including those for cellular processes, environmental information processing, genetic information processing, metabolism, and organismal systems (Fig. 4). The greatest number of annotated unigenes were related to metabolism, followed by those involved in genetic information processing. The least number of unigenes belonged to environmental information processing. In the second level of pathway enrichment, the highest number of unigenes were those that are involved in translation and carbohydrate metabolism and the ones involved in membrane transport were the lowest. The top 5 metabolic pathways, based on the number of unigenes involved, were for carbon metabolism, plant-pathogen interaction, ribosome function, amino acid biosynthesis and RNA transport.    further annotated and classified into 20 sub-categories, where DEGs involved in metabolic processes were the most abundant ones followed by the ones in single-organism process and single-organism metabolic process (Fig. 7a). 2005 DEGs identified in FT treatment were mainly classified into 20 sub-categories, mostly related to metabolic processes in general, and the metabolic processes related to single organism and organic substances. Processes related to protein phosphorylation, general metabolism, metabolism of phosphorus and phosphatecontaining compounds, oxidation-reduction reactions, single-organism metabolic processes, hormonal control, and metabolism, were all found to be enriched in CT and FT treatments (Fig. 7b). Genes involved in cell recognition, pollen-pistil interaction, recognition of pollen, pollination, cell redox homeostasis, single and multicellular organism processes, lipid metabolism, cellular homeostasis and steroid metabolism were specifically enriched in CT. Cellular proteins and other macromolecule modification processes, carbohydrate metabolism and biosyn-     www.nature.com/scientificreports/ host intracellular part, intracellular region of host and extracellular region were found to be strictly CT enriched, while synapse and endoplasmic reticulum lumen were only enriched in FT treatment. A total of 1108 DEGs after CT treatment fell into 20 sub-categories in MF with the most abundant pathway related to catalytic metabolism. 965 of the observed genes that are involved in catalytic, ion binding, transferase, small molecule binding and ATP binding activities were significantly enriched in FT condition. Protein kinases, phosphotransferases with carboxyl group as acceptor, kinases, catalytic proteins, oxidoreductases, calmodulin binding proteins, antioxidant molecules, and chlorophyllase were all enriched in both, CT and FT conditions. CT specifically enriched the metabolic pathways involving various coenzymes, cofactors, nucleic acid binding transcription factors, transcription factor activity-sequence-specific DNA binding factors, heme binding factors, tetrapyrrole binding factors, carbon-carbon lyases, protein disulfide oxidoreductase, disulfide oxidoreductase and peroxidases. While only FT caused enrichment of various transferases, xyloglucan: xyloglucosyl transferase, small molecule binding factors, chlorophyllide A oxygenase [overall activity], ion binding factors, 3-beta-hydroxy-delta 5-steroid dehydrogenase, steroid dehydrogenase, ATP binding factors and cellulose synthase.
KEGG pathway analysis of DEGs. Compared with CK, the 1912 DEGs in CT treatment were found to be involved in 117 pathways whereas, the 1075 DEGs identified in FT were involved in 109 pathways. Out of the top 20 enriched pathways enriched both in CT and FT 15 pathways (Fig. 8), including photosynthesisantenna proteins, glycine/serine and threonine metabolism, alpha-linolenic acid metabolism, photosynthesis, glutathione metabolism, carotenoid biosynthesis, pyruvate metabolism, circadian rhythm-plant, plant hormone signal transduction, pentose phosphate pathway, starch and sucrose metabolism, plant-pathogen interaction, stilbenoid/diarylheptanoid and gingerol biosynthesis, phenylalanine metabolism, glycolysis /gluconeogenesis. Besides, five pathways were specifically enriched in CT, including carbon fixation in photosynthetic organisms, nicotinate and nicotinamide metabolism, citrate cycle (TCA cycle), nitrogen metabolism, and galactose metabolism. Pathways related to linoleic acid metabolism, porphyrin, and chlorophyll metabolism, phenylpropanoid biosynthesis, amino and nucleotide sugar metabolism, and flavonoid biosynthesis were found to be enriched in FT only.
Compared with CK, there were 10 metabolic pathways (involving 232 genes) significantly enriched in CT treatment (p < 0.005). These metabolic pathways include proteins related to photosynthesis-antenna, glycine/serine, threonine and alpha-linolenic acid metabolism, photosynthesis, carbon fixation in photosynthetic organisms, glutathione metabolism, carotenoid biosynthesis, pyruvate metabolism, circadian rhythm-plant, nicotinate and nicotinamide metabolism (Table 3). In these metabolic pathways, 18 DEGs in photosynthesis-antenna proteins were downregulated. Glycine/serine and threonine metabolism contained 17 down and 18 upregulated genes. Alpha linolenic acid metabolism enrichment showed 11 downregulated and 17 upregulated genes. Among the 32 and 41 DEGs involved in photosynthesis and carbon fixation in photosynthetic organisms, 30 and 29 DEGs, respectively, were inhibited. 22 and 12 genes were down and upregulated respectively, that are involved in glutathione metabolism whereas, 13 and 7 down and upregulated genes are involved in carotenoid biosynthesis. 18 genes were found to be inhibited and 25 gens were induced in Pyruvate metabolism. Among the 28 DEGs of circadian rhythm in plant, 25 DEGs were found to be downregulated. Nicotinate and nicotinamide metabolism contained 7 downregulated genes and 9 upregulated genes. Most of DEGs related to photosynthesis or biological rhythm management were downregulated (accounting for 85.71% of the total DEGs in these metabolic pathways) www.nature.com/scientificreports/ whereas, other metabolic pathways regulating photosynthesis-antenna proteins, photosynthesis, carbon fixation in photosynthetic organisms, plant circadian rhythm, were found to be inhibited. The photosynthetic capacity of C. album was found to be decreased after chilling injury, and the biological rhythm adjusted to adapt to the environmental changes. Glutathione metabolism, pyruvate metabolism, nicotinate and nicotinamide metabolism are all closely related to plant antioxidant or oxidative responses. Their significant enrichment after chilling stress indicated that they play important roles in the resistive responses of C. album towards chilling stress. Compared with CK, there were 9 significantly enriched metabolic pathways including 219 DEGs in FT treatment. These metabolic pathways include those related to plant-pathogen interaction, plant circadian rhythm, alpha-linolenic acid and linoleic acid metabolism, glycine/serine and threonine metabolism, photosynthesisantenna proteins, carotenoid biosynthesis, plant hormone signal transduction, and pentose phosphate pathway (Table 4). In these pathways, 52 of the 58 DEGs in plant-pathogen interaction were up-regulated. 20 downregulated genes and 4 up-regulated genes were involved in plant circadian rhythm. Alpha-linolenic acid metabolism pathway showed 11 down-regulated and 7 up-regulated genes, whereas 7 out of 8 DEGs in linoleic acid metabolism were up-regulated. 10 genes showed up-regulation and 10 were down-regulated in glycine/serine and threonine metabolism. 10 genes involved in photosynthesis-antenna proteins were down-regulated. 6 genes involved in carotenoid biosynthesis were down-regulated and 7 were up-regulated. Plant hormone signal transduction pathway had 23 down-regulated and 22 up-regulated genes. 11 genes were down-regulated whereas, 8 were up-regulated in pentose phosphate pathway.
In these significantly enriched metabolic pathways, those of glycine/serine and threonine metabolism, alphalinolenic acid metabolism, carotenoid biosynthesis, photosynthesis-antenna proteins, and plant circadian rhythm  www.nature.com/scientificreports/ were significantly enriched in both CT and FT treatments. Photosynthesis, carbon fixation in photosynthetic organisms, glutathione metabolism, pyruvate metabolism, nicotinate and nicotinamide metabolism were specifically enriched in CT condition, while pathways involved in plant-pathogen interaction, linoleic acid metabolism, plant hormone signal transduction, pentose phosphate pathway were specifically enriched in FT treatment.
Notably, the DEGs involved in circadian rhythms of plant and photosynthesis-antenna proteins were mostly down-regulated, indicating that both chilling and freezing temperatures inhibited these pathways in C. album. Although the number of DEGs in alpha linolenic acid metabolism, glycine/serine and threonine metabolism, carotenoid biosynthesis in chilling and freezing stress was large, the overall differences were not significant. However, compared with the chilling stress, the metabolic pathways for glutathione, pyruvate, nicotinate and nicotinamide were not significantly enriched in FT, suggesting that the changes in these metabolic pathways play more important roles in chilling responses of C. album. Alternatively, freezing treatments caused significant enrichment of genes involved in linoleic acid metabolism, pentose phosphate pathway, and plant hormone signal transduction, indicating that these pathways are involved in resistance to freezing temperature by C. album. There were especially 45 DEGs serving in pathways related to plant hormone signal transduction, suggesting that the hormones are indeed essential in developing resistance to freezing stress in C. album.
Quantitative real time PCR (qRT-PCR) analysis of DEGs. The expression levels of 39 genes, selected from enriched metabolic pathways from KEGG analysis were further studied. Compared with CK, the expression levels of 19 genes in CT and FT were significantly increased, while 20 genes were significantly downregulated both in CT and FT treatments (Fig. 9b). The relative expression levels of most of the genes verified by qRT-PCR were consistent with the RNA-Seq data (Fig. 9a) further validating the results obtained in RNA-Seq.
Verifying the DEGs involved in hormone signal transduction pathways. Under low temperature stress, plant hormone signal transduction pathway was significantly enriched in response to freezing injury but not to the chilling injury. The differential expression of 13 key genes from auxin, cytokinin, gibberellin, abscisic  Table 5. www.nature.com/scientificreports/ acid, ethylene, and salicylic acid signaling pathways were also analyzed by qRT-PCR (Fig. 10). Compared with CK, CT did not significantly change the expression level of Aux/IAA in auxin signal transduction pathway. But its expression was increased more than eightfold in FT. The expression of GH3 increased 16 and 18 times respectively, in CT and FT, while the relative expression level of SAUR was decreased to about 1/6 and 1/2 of that in CK. Although ARR -B from cytokinin signal transduction pathway increased a little in CT, the difference was not significant. However, its relative expression level in FT was increased more than 10 folds. A key gene of gibberellin signal transduction pathway, DELLA was found to be downregulated in CT, but upregulated, nearly 2 folds, in FT condition. PP2C and SnRK2, two key genes of abscisic acid signal transduction pathway, had the similar response pattern to low temperature stress, both being down-regulated in CT and up-regulated 3 and 1.5 times, respectively, in CT and FT. ETR, EBF1/2 and EIN3, three key genes of ethylene signal transduction pathway, also changed expression after exposure to low temperature. The relative expression of ETR in CT increased while that in FT did not change significantly. The expression levels of EBF1/2 and EIN3 consistently decreased in CT but increased 4 and 2 times, respectively, upon FT exposure. TGA and PR-1, the key genes of salicylic acid signal transduction pathway, responded to low temperature stress differently, in which TGA was significantly up-regulated (3 times) in CT, while PR-1 increased 11 and 27 times in CT and FT both, respectively. JAR1 is the key gene of jasmonate signal transduction pathway and the expression level of JAR1 also changed during low temperature stress in C. album. These results further substantiate our findings that auxin, cytokinin, gibberellin, abscisic acid, ethylene, jasmonate and salicylic acid signal transduction pathways are involved in C. album responses to low temperature stress, especially to freezing injury.

Discussion
De novo assembly of the C. album transcriptome. High-throughput transcriptome sequencing, a modern biology research tool, is a common and popular method for screening DEGs, constructing metabolism networks and studying molecular regulation in plants 20,21 . In our study, we firstly reported the transcriptomic analysis of C. album, where a total of 18.97 G clean data was obtained and the base error rate, Q20, Q30 and, GC content all are in accordance with the established requirements. Moreover, the N 50 and N 90 values for assembled and spliced genes were 2,543 bp and 648 bp, respectively. 64.44% of unigenes were annotated in at least one of the databases-Nr, Nt, Pfam, KOG, Swiss-prot, KEGG or GO. Therefore, with this clean transcriptomic data in C. album we have tried to provide accurate references for future research in this area.
Comparative transcriptome analysis reveals the responses to low temperature in C. album. Three transcriptome libraries of C. album were constructed and the transcript changes in C. album under chilling and freezing treatments were studied. Compared with control (CK), 2810 genes were upregulated and 2567 were downregulated after chilling stress treatment (CT). Whereas 1748 and 1459 genes were up and down-regulated, respectively, after freezing stress treatment (FT). The numbers of up-regulated genes in CT and FT were both larger than that of down-regulated genes, indicating that C. album might improves cold resistance by positive regulation of these genes. Moreover, commonly up-regulated 1092 genes and 1148 down-regulated genes in both in LT and FT treatments were observed. There was a specific up-regulation and down-regulation of 1718 and 1419 genes in CT, whereas in FT the numbers were 656 and 311, respectively, indicating a mutually exclusive response in these cases. In addition, the number exclusive and non-exclusive DEGs in FT were significantly less than that in CT, which could be the reason behind the intolerance of C. album to freezing injury because of compromised life activities under freezing stress.

Pathways involved in both chilling and freezing responses in C. album. GO classification and
KEGG enrichment analysis are widely used tools for transcriptomic research. The chilling and freezing responses of C. album, as seen in these analyses separately, showed some similarities but bore differences as well. The GO results showed that the shared responses mostly include DEGs related to various metabolic processes (BP category) including single-organism metabolic process, as well as different molecules catalyzing these processes. In contrast, DEGs solely related to FT, were specifically enriched in organic substance metabolic process (BP category), ion binding, transferase activity, small molecule binding and ATP binding, etc., that particularly regulate freezing acclimation in C. album. The GO classification results indicated that the responses of C. album to chilling injury and freezing injury are quite different. Pathways that were significantly enriched after CT and FT treatments, as seen in KEGG pathway analysis were related to or involved in lysine/serine and threonine metabolism, alpha-linolenic acid metabolism, carotenoid biosynthesis, photosynthesis-antenna proteins, and plant circadian rhythm. However, those specifically enriched after CT are photosynthesis, carbon fixation in photosynthetic organisms, glutathione metabolism, pyruvate metabolism, nicotinate, and nicotinamide metabolism. FT treatment particularly enriched pathways regulating plant-pathogen interaction, linoleic acid metabolism, plant hormone signal transduction and pentose phosphate pathway. These results indicated that some of metabolic pathways were enriched in both the stresses (CT and FT) whereas some were quite exclusive to only one kind of stress exposure.
Previous reports have shown that the changes in pathways controlling plant circadian rhythm could help the plant in adjusting to a change in temperature such as cold stress 22 . This involves increasing expression levels of key genes involved in these pathways 8,22,23 . Serine and threonine protein kinases (STPK) like SpkB, SpkD, SpkE and SPkG, that function as transcriptional regulators have also been shown to be involved in responses towards low temperature stresses in different organisms 24 . Plant-pathogen interaction is closely related to cold acclimation of grape plant that is put under cold stress (4 °C). Such a treatment increases expression levels of CNGC, CMLs, JAZ1, Rboh, FLS2, BAK, MEKK1, MKKs, RPM1, RPS5, RIN4, PBS1, WRKYs, PR1 and MIN7 as a response to cold resistance 25 . The accumulation of linolenic acid in peach fruits improves its tolerance to low temperature stress 26 . Low temperature leads to decrease in CO 2 absorption in many cold-sensitive plants, resulting in the degradation of photochemical efficiency during photosynthesis 27,28 . Presence of carotenoids affects the cold resistance of plants and may also be related to the changes in auxin content, oxidative damage, and membrane permeability 29 . Plant hormone signal transduction pathways, involving ABA, SA, JA, ETH, GA, IAA and CTK are widely recognized to play important part in cold and disease resistance, heat tolerance and other biological or abiotic stress processes in plants 30 . In addition, linoleic acid, pentose phosphate, glutathione, pyruvate, and nicotinamide [31][32][33][34][35] have all been shown to be closely related to low temperature stress responses or cold resistance in plants. All these past studies suggest that enrichment of such metabolic or signal transduction pathways are closely related to low temperature stress responses in C. album. Aux/IAA, GH3 and SAUR are the key genes for Aux/IAA-TIR1 nuclear signaling pathway. Aux/IAA dimers with Auxin Response Factor thus inhibiting its transcriptional regulation functions. GH3 encodes enzymes catalyzing auxin and amino acid coupling, and hence maintain auxin homeostasis through feedback regulation and SAUR affects functioning in some highly unstable transcriptional lines. The cold acclimation process in Arabidopsis thaliana and low temperature stress in soybean involves changes in the expression levels of auxin signal Cytokinin is an important plant hormone that regulates plant growth and development as well as responses to abiotic stress though key genes like ARR -A and ARR -B. The expression of ARR18 gene, a positive regulator of low temperature responses, is significantly increased in rubber plant after low temperature treatment 5 . We also found an increased expression level of ARR-B gene in C. album after chilling and freezing stress, which positively correlated with the extent of temperature downregulation.

The plant hormone metabolism and signal transduction mediate FT tolerance in
A protein named DELLA is the core component of gibberellin signal transduction pathway and plays a negative role in regulating GA signal transduction pathway 48 . Gibberellin plays an important role in plant responses to low temperature stress 46,47 . DELLA enables plant survival through adverse environmental conditions by regulating plant hormones 49 . The tolerance of plants to stress depends on the amount of DELLA in vivo. When the active gibberellin is limiting, DELLA level increases and the tolerance to stress improves. In tomato, the crosstalk of PIF4 (Phytochrome interaction 4) and DELLA regulates CBF transcription and hormone homeostasis during cold stress response by integrating light and temperature signaling and sharing of hormone pathways and transcription regulators thus enabling better adaptation of plants to cold stress 50 . We also found an increased expression of DELLA in qRT-PCR after freezing injury, suggesting its role in negatively regulating gibberellin content in vivo in C. album to adapt to low temperature stress.
SnRK2 is a plant specific protein kinase, which plays an important role in ABA signal transduction and osmotic stress regulation 51 . Low temperature could induce the SnRK2 gene in tobacco and wheat 52,53 and overexpression of TaSnRK2.4 in Arabidopsis improves its cold resistance through a series of physiological changes such as reduction in rate of water loss, enhanced stability of cell membrane and improvement in the photosynthetic and osmotic potential 54 . PP2C is a serine/threonine residue protein phosphatase and a core negative regulator of ABA signal transduction. In Arabidopsis, PP2C has also been shown to be closely related to the low temperature acclimation 55 . An enrichment of ABA signaling pathway in our study further confirms its involvement in low temperature stress responses and subsequently validates our results too.
Being the key enzyme for jasmonate signal transduction pathway, JAR1 catalyzes the last step forming molecularly active JA-Ile from the combination of jasmonate and amino acids, especially isoleucine 56 . LOX, AOS and JMT are other key enzymes in jasmonate biosynthesis or signal transduction pathway. Previous studies have shown that low temperature could induce the expression of LOX, AOC and JAR1 genes in Artemisia annua, thus increasing the level of endogenous jasmonate, which consequently induces the expression of ETR1, ETR2 and AP2/ERF 57 . We also found an increased expression of these genes in C. album, as response to chilling or freezing injury.
ETR, EBF1/2 and EIN3 are genes serving in ethylene signal transduction pathway. After cold stress in kiwifruit and papaya, expression levels of some of the ETR and EIN3 members are increased, with a slight change in the expression level of EBF as well 58,59 . As a key member of salicylic acid signal transduction pathway, transcription factor TGA is involved in temperature stress responses in some plants 60 , while PR-1 gets expressed in others 61 . Our results corroborate with these earlier findings further strengthening the role of ethylene and salicylic acid signal transduction pathways in cold stress responses in C. album.

Materials and methods
Plant materials and low temperature treatments. Healthy, strong and uniform C. album. cv. 'Fulan 1' seedlings of height of ~ 20 cm were selected and kept in growth chambers (GXZ-280C, Ningbo, China) maintained at 25 °C, 60% to 80% relative humidity, and a photoperiod of 12 h (2000 ± 200 lx) for one week. Then the seedlings were subjected to three different temperature treatments, i.e. freezing stress treatment (FT, − 3 °C) in freezer, chilling stress treatment (CT, 4 °C) in freezer and control check treatment (CK, 25 °C) in growth chambers for 24 h. Ten seedlings were collected for each treatment and each experiment was repeated 3 times. During www.nature.com/scientificreports/ the whole treatment process, all seedlings were grown in the same conditions as formal, and the samples were simultaneously collected in dark after treatment. Compared with the seedlings in CK (Fig. 11a), the leaves from seedlings in FT experienced dehydration and gradually wilted (Fig. 11c), while there was no obvious phenotypic change in leaves from seedlings grown in CT (Fig. 11b). After the 24 h temperature treatment of the seedlings, the 3rd and 4th leaves from top to base were collected (from all the replicative experiments) and instantly frozen in liquid nitrogen. The samples were stored in the freezer at -80 °C till further use. The collection of plant material in this research, complied with relevant institutional, national, and international guidelines and legislation.
Total RNA extraction and detection. Total RNA was extracted by E.Z.N.A. Plant RNA Kit (OMEGA bio-tek, USA) as per the manufacturer's instructions. The RNA was checked for quality and quantity by first running it on a 1.5% agarose gel, followed by a test for purity, concentration and integrity using the TU-1810 spectrophotometer (Puxi, CHN), NanoDrop OneC Flurometer (Thermo Fisher Scientific, USA) and the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA). Only high-quality total RNA was used for further experiments including the library construction.
Library preparation for RNA-Seq. 1

Clustering, sequencing, and data analysis. Clustering of the indexed-coded samples was done via
cBot Cluster Generation System that uses TruSeq PE Cluster Kit v3-cBot-HS from Illumina. The library preparations were sequenced on Illumina Hiseq-2000 platform and paired-end reads were generated after the cluster generation. Followed processing of fastq formatted raw data (raw reads) by internal perl scripts, and clean data (clean reads) were obtained by removing reads with adapter, ploy-N and low-quality raw data. Subsequently, Q 20 , Q 30 , GC-content and sequence duplication level for the clean data were calculated. The high-quality clean data were then used for downstream analysis. Transcriptome was assembled using Trinity with min_kmer_cov default value of 2 and all other parameters were set to their respective default values 63 . Gene function was annotated using Nr, Nt, Pfam, KOG, Swiss-Prot, KEGG and GO databases. The read counts for each sequenced library was adjusted by edgeR program package through one scaling normalized factor. Differential expression analysis of two samples was performed using the DEGseq R package as shown previously 64 . P value was adjusted by using the q value 65 . q-value < 0.005 and |log 2 (fold change) |> 1 was set as the threshold for calculating significantly differential expression of a gene.
GO enrichment analysis of the differentially expressed genes (DEGs) was performed according to the GOseq R packages based Wallenius non-central hyper-geometric distribution 22 . The statistical enrichment of DEGs in KEGG pathways was tested by KOBAS software as shown previously 66,67 . qRT-PCR analysis of DEGs. TransScript All-in-One First-Strand cDNA Synthesis SuperMix kit (TransGen Biotech, CHN) was used to synthsize high quality cDNA for qRT-PCR analysis. cDNA samples of CK, CT and FT were mixed equally, and then gradiently diluted 4 times to obtain a standard curve (dilutions: 10 -1 , 40 -1 , 160 -1 and 640 -1 ) using a suitable annealing temperature, that was screened according to the amplification efficiency. Eppendorf Realplex 4 and TransStart Top Green qPCR SuperMix were used to detect the relative expression of the selected 39 genes (including 13 plant hormone signaling and transduction pathway related genes, i.e. AUX/ IAA, GH3 and SAUR involving in auxin signaling pathway, ARR -B involving in cytokinin signaling pathway, DELLA involving in gibberellin signaling pathway, PP2C and SnRK2 involving in abscisic acid signaling pathway, ETR, EBF1/2 and EIN3 involving in ethylene signaling pathway, TGA and PR-1 involving in salicylic acid signaling pathway and JAR1 involving in jasmonate signaling pathway) in different treatments by using primers as listed in www.nature.com/scientificreports/ dissolution curve. The reaction was set as with 10 μL 2 × TransStart Top Green qPCR SuperMix, 0.4 μL Passive Reference Dye (50 ×), 0.2 μM each upstream and downstream primer, and 0.5 μM cDNA template. ACTB7 and TUB5 were used as reference genes 68 . SPSS19.0 was used for analyzing the results statistically. TBtools was used to draw heatmaps according to the relative expression and FPKM value of the genes.

Gene ID Gene name Primer sequence (5'-3') Description
Cluster-7340.50923 FBA  CCC ATC TGC CCT TAC TGT TC  Fructose-bisphosphate aldolase  GGT TAG CCT CGG TGT TCT Table 5. List of qRT-PCR primers used in the study.