Comparative transcriptome profiling reveals the cold stress responsiveness in two contrasting Chinese jujube cultivars

Low temperature is ( Ziziphus jujuba Mill.) in cold winter and spring. Little is known about the molecular mechanisms for coping with different freezing stress conditions in jujube. To gain insight into the freezing-related molecular changes, we conducted comparative transcriptome analyses from the cold-sensitive cultivar ‘Dongzao’ and cold-tolerant cultivar ‘Jinsixiaozao’ using RNA-Seq. Results In our study, more than 20,000 genes were detected at chilling (4°C) and freezing (-10°C, -20°C, -30°C and -40°C) stress between two cultivars. The numbers of differentially expressed genes (DEGs) between the two cultivars was 1831, 2030, 1993, 1845 and 2137 under the five treatments, respectively. Functional enrichment analysis suggested that metabolic pathway, response to stimulus and catalytic activity were significantly enriched at stronger freezing stress. Among them, nine DEGs participated in the pathway of Ca 2+ signal, thirty-two DEGs were identified to take part in sucrose metabolism, and other DEGs were identified to participate in the regulation of ROS, plant hormones and antifreeze protein. In addition, important transcription factors ( WRKY , AP2 / ERF , NAC and bZIP ) participating in the freezing stress were activated under different degrees of freezing stress. Our research provides a more comprehensive understanding of DEGs involved in freezing stress at the transcriptome levels in Z. jujuba , especially two cultivars with different cold tolerance. These results expanding our understanding on the complex molecular mechanism of jujube, which also provide new insights and candidate genes for genetic improvement of jujube tolerance to freezing stress. Our present study presented a comprehensively analysis of cold responsive genes via comparative transcriptome analysis using two contrasting Chinese jujube cultivars (‘Dongzao’ and ‘Jinsixiaozao’) under different freezing stress. Further analysis of data identified some important DEGs, which involved in galactose metabolism, signal transduction, transcriptional and post-translational regulation, were differentially expressed between two cultivars to resist freezing stress. And the involvement of TFs leads to a series of physiological and cell response, and improved tolerance to freezing stress. The transcriptional level of ‘Jinsixiaozao’ changed greatly during freezing stress may explain they were better cold tolerance compared to ‘Dongzao’. These results provide a reference for our understanding of Chinese jujube adaptation to different freezing stress and provide new insights into the molecular mechanisms of different freezing tolerance in jujube cultivars. and 72℃ for sec. The baseline and threshold cycles (Ct) were automatically determined by the system’s All the experiments were performed in biological triplicate. Transcript levels were normalized against the average expression of the Zjactin gene (GenBank: EU916201).

huge impact on the survival and geographical distribution of plants [3]. More than half of the 350,000 plant species on Earth are grown in the tropics and subtropical regions, such as rice, maize, and tomato which are chilling sensitive. And plants from temperate region are not very tolerant to freezing, they can increase their freezing tolerance by being exposed to chilling, a process known as cold acclimation [4,5].
In previous studies, researchers have revealed that cold exposure results in a variety of alterations in physiology and gene expression patterns, and hundreds to thousands of genes even were identified that are up-and down-regulated in response to cold, including the metabolism of carbohydrates, molecular chaperones, antifreeze proteins, signal transduction (receptor kinase, protein kinase/phosphatase, Ca 2+ binding protein) and regulatory proteins [6][7][8][9]. In addition, transcription factors could regulate genes expression [10][11][12]. Some of these genes were successfully transferred to cold or freezing-sensitive crops, such as Arabidopsis thaliana, rice, to improve yields under adverse growing condition [13,14]. In rice, OsCDPK7 expression is induced by cold and salt stress, and overexpression OsCDPK7 increased tolerance to cold stress [15]. In other species, overexpression of some CBF genes in Arabidopsis and other plants conferred enhanced tolerance to cold stresses [16][17][18][19][20][21]. Also, in barley, TaCBF14 and TaCBF15 could enhanced freezing tolerance [22].
Chinese jujube, native to China, is one of the country's most important fruit crops. In recent years, although global temperature is increasing, abnormal weathers, such as cold extremes, have been observed to increase all over the world, especially in main gain crops and industrial crops production areas [32,33]. These disasters seriously hinder the development of jujube industry [34,35]. To understand the jujubes respond to freezing stress responses is crucial to provide valuable information for breeding cold tolerance jujube cultivars. While many cold stress studies have been conducted in model experimental plants, such as Arabidopsis, rice, potato, grape and apple [36][37][38], little is known about the molecular mechanisms of freezing tolerance in Z. jujuba. Thus, in the present study, we used RNA-seq technology to compared the difference between the two Chinese jujube, the coldsensitive cultivar 'Dongzao' and cold-tolerant cultivar 'Jinsixiaozao', under different degrees of freezing stress. We identified the freezing stress related candidate genes and metabolic pathways involved in chilling and freezing tolerance by comparing two cultivars at different freezing stresses.
Many of the genes are involved in photosynthesis, hormone signaling and plant immunity, and the key genes expression patterns are associated with greater cold tolerance. The objective of this work was to better understand the molecular regulation that occur under cold tolerance using different degrees of low temperature and specific cultivar to cultivar comparisons, which providing a rich genetic resource for further study of cold response in jujube.

Results
Quality analysis and sequence assembly of RNA-seq data Total RNA extracted from the xylem of jujube branches was employed to build expression libraries for sequencing by using the Illumina High-Seq sequencing system. We obtained 30 libraries reads that each library has about 89 million in total. After filtering, clean reads varied from 87 million to 88 million, which occupied on average 97.47% in total reads of all libraries, were collected (Table 1).
In all samples, the Q30 value was all over 86%, and the GC content was 43%. Genome map rate which accounted for more than 73%, and gene map rate could reach from 72.93-76.32%. A total about 20,000 genes were detected to be expressed in each library. In addition, more novel transcripts were both detected in DZD and JSD. Then, about alternative splicing, each sample has a similar number of event. Indel, JSD has the maximum mark, significantly more than in DZC. Three biological replicates were collected for each cultivar sample under each temperature. Pearson's correlation coefficients among the three technical replicates replication samples were high (γ = 0.91-0.98), indicating that the sequencing quality allowed for further analysis.  In five different degrees of cold stress, as shown in Fig. 2, 854 common DEGs were detected in both cultivars in response to five treatments, that may be suggests these DEGs are involved in the process of dealing with low temperature under different low temperature, and both participate to the same pathways in response to freezing stress. In addition, at -40℃, the number of unique DEGs (600) are significantly higher than other group, indicating that these genes responded differently and thus are candidate contributors to the difference in cold tolerant between 'Dongzao' and 'Jinsixiaozao'.
In this research, RNA sampled from branch for qRT-PCR validation of the sequence based transcription profiles for six DEGs. As shown in Fig.S1, all the genes showed a similar expression pattern between the RNA-seq data and the qRT-PCR results, which indicated that these genes are really expressed and the changes at different temperature are real. Also the results confirmed that the RNA-seq data are highly reliable for further analysis.

Function annotation and enrichment analysis of DEGs
To understand the functions of discovered DEGs between two cultivars, all DEGs were searched against the protein database, followed by Gene Ontology (GO) analysis to assess the genes' putative functions. The top GO terms in each treatment, ranked by the number of DEGs, the stronger freezing stress, varied with freezing treatment stages. The more genes were annotated, and the metabolic processes were modified and changed (Fig. 3). In the comparative analysis of these five groups, GO terms of shared DEGs between the two cultivars could reflect common responses to freezing stress.
Of these GO terms, the majority of DEGs were classified into 'metabolic process', 'cellular process' and 'response to stimulus'; and 'cell', 'cell part' and 'membrane' in cell components; 'catalytic activity', 'binding', 'transporter activity' in molecular function.
In addition, to gain insights into the metabolic or signal pathways involved in freezing stress, 20 common metabolism processes are most enriched at five freezing treatments in two cultivars (Fig. 4).
We also found different numbers of DEGs in the same pathway between two cultivars. Among these top pathways, carbohydrate metabolism pathway, galactose metabolism, fructose and mannose metabolism were associated with many up/down regulated genes. Furthermore, in the comparison of the five groups, there are some pathway enrichment in the stronger freezing stress in 'Jinsixiaozao', such as peroxide of metabolic pathways, which is only at -30℃ significant enrichment, ABC transporter process is at -40℃. In brief, in the functional annotation of DEGs, many genes of different metabolic pathways were involved in the regulation of freezing stress in jujube, and further study on the differential expression pattern in these pathways has a great significance for revealing the cold resistance of jujube.
Whether jujube trees can safely pass the winter is an important factor affecting the yield and quality of the next year. Therefore, in order to further analyze the difference in jujube at extreme temperatures, we conducted a comparative analysis of the two cultivars with temperature of -30℃ and − 40℃. In comparison, we found that more significantly DEGs were found in the stronger freezing  In the comparison of four groups of 'Jinsixiaozao', JSCK as a control, a total of 1668 DEGs were identified (Fig. 5, Fig. S3 and Additional file 5), and 86 common DEGs were present in all comparisons.
Then we focused analysis on the significant DEGs (PPEE < 0.05 and |Log2|≥1), a total of 71 significantly DEGs were identified at -10℃, 39 genes were up-regulated and 32 were down-regulated.
In the case of -20℃, 108 DEGs were identified, including 31 up-regulated genes and 77 down- pathway showed significant differences to cope with the adverse effects of low temperatures.
Through annotation of the target genes and comparison with KEGG database, some genes were shown involving in metabolic pathway and biosynthesis of galactose. As listed in Table S2, among them, 107411641 is a key gene in the process of UDP-galactose to form galactionl, which is significantly up-regulated in stronger freezing stress, that will enhance the synthesis pathway of stachyose. However, in the pathway of stachyose decomposes into galactose and raffinose, FPKM value for 107425264, 107430302, 107414011 and 107428511 showed decrease, which reduces the decomposition of stachyose (Fig. S4).

Discussion
Over the past years, extreme weather occurs frequently, the uncontrollable climate affected the growth and development of jujube, even caused their death. Multiple responses to cold stress have been reported in many plant species, reflected by the variations in transcription or metabolism [39,30]. The cold tolerance of plants depends on cellular signal transduction pathways. Recently, significant advances have been made regarding the understanding of plant cold-signaling perception and transduction [40][41][42]. In our study, the systematic transcriptome profile between 'Jinsixiaozao' and 'Dongzao' was reported to identify the genes responsible for freezing resistance, which could help the breeding for future development of improved cold-tolerance in jujube cultivars.

Identification of DEGs in two cultivars
In two cultivars, some DEGs were investigated and hundreds of freezing-related genes were identified

Physiological responses related genes under freezing stress
Freezing stress may cause plant death. Some research shows that, under this conditions, plants could increase their tolerance by synthesizing numerous protective substances and proteins [44][45][46][47].
Raffinose family oligosaccharides (RFOs) are a class of functional oligosaccharides unique to plants, mainly including raffinose, stachyose and verbascose. Some study found that RFO can improve the ability of plants to cope with adversity. The more RFOs accumulate, the stronger the ability of plants to be exposed to low temperatures. Galactose and related enzymes are involved in the entire RFO metabolic pathway [48]. In our study, galactose metabolism was significantly enriched, and relative genes in this pathway have up or down regulation to cope with cold stress. Among these genes, 107426117 (stachyose synthetase), 107415484 (raffinose synthase) and 107411641 (inositol 3-alphagalactosyltransferase) were significantly up-regulated under stronger freezing treatment in 'Jinsixiaozao', however, some of them did not significant change under the weak freezing stress. That may be suggested these genes play more important role in galactose metabolism in response to cold, especially under stronger freezing stresses. In addition, late embryogenesis abundant protein (LEA) could respond to abiotic stresses such as low temperature and drought [49]. In our study, the upregulation of LEA (Gene ID:107409127) indicates that LEA genes are involved in freezing tolerance in 'Jinsixiaozao' (-40℃). The function of LEA protein in this study should be further studied as a regulatory mechanism of freezing response.

Cold signal transduction related genes under freezing stress
Signal transduction pathway is known to play crucial roles in response to the stress of freezing, among these, Ca 2+ is an important second messenger in plant [50,51]. In our research, 34 DEGs were identified in response genes related to calcium-binding proteins, calmodulin and Ca 2+ signal.
Among them, gene (Gene ID: 107413690) showed a significant up-regulation in 'Dongzao', which indicating these genes play important role in jujube resistance to freezing stress and the function of these genes should be further studied.
In addition to Ca 2+ , ROS also play important roles [52]. On the one hand, they induce genes expression and protein synthesis to protect cells from stress; on the other hand, they induce oxidative stress [53,54]. Meanwhile, ROS also can directly stimulate the expression of genes in activating stress-related pathways, such as MAPK cascade pathway, hormone signaling pathway and the calcium signaling pathway. In our study, some genes with antioxidant activity are involved in freezing stress.
Among these DEGs, six genes (one up-regulated gene and five down-regulated genes) were involved in the redox signaling, such as gene 107403552, which are differentially expressed in four degrees of freezing stress in 'Jinsixiaozao', indicating they played different roles in cold signal transduction under freezing stress.

Transcriptional regulation related genes under freezing stress
Transcription factors (TFs) play important functions in plant development and stress tolerance. Many transcriptomic data depicted some TFs that govern the massive and highly coordinated transcriptional changes under different freezing stress. Some studies show that most cold-regulated transcription factor genes are impaired in response to cold stress in Arabidopsis [55][56][57]. A lot of research shows that some transcription factors play important role in plant freezing tolerance [58][59][60]. Such as AP2/ERF factors RAP2.1 and RAP2.6, and the C 2 H 2 -type zinc finger STZ/ZAT10 belong to the CBF regulon [61][62]. Furthermore, NAC was also shown to enhance expression during loquat fruit chilling-induced, however, it could not bind directly to promoters of lignin biosynthesis genes [63]. Many MYB-type transcription factors, for example AtMYB85 [64], AtMYB88 [65,66], increases cold hardiness by CBF-dependent and CBF-independent pathways in apple. In other species, OsMYB4 transgenic lines of Arabidopsis exhibited improved chilling and freezing tolerance with a dwarf phenotype [67]. In this study, we validated AP2/ERF, WRKY, NAC and bZIP gene families and analyzed the expression patterns, found that different members play different roles in two cultivars. other TF members are generally divided into two expression patterns. These differences suggest that the two cultivars have different responses for cold stress, even for different degree of low temperature. We did not go further to verify the function of these genes, but listed the significant differences in this work. Thus, the emphasis of future research could focus on the role of these genes in metabolic pathways, and understand the regulatory mechanism of genes.
In addition, some studies have found that alternative splicing (AS) events are important post-transcriptional regulation and are required for reprogramming gene expression under environmental stress. Previous studies showed that the number of AS events was significantly higher compared to control plants when plants exposed to different stresses, especially low temperatures [68]. To analyze the variation of AS under cold stress, in this study, we detected that all genes containing intron were alternately spliced under cold stress (Fig. S5). On the one hand, compared with the control, the number of AS events increases slightly with the temperature drops, we hypothesize that plants can fine-tune gene expression by AS events in response to temperature changes. On the other hand, most AS genes were not differentially regulated by stress, suggesting that AS may represent an independent gene regulatory layer that responds to cold stress.

Post-translational regulation related genes under freezing stress
In addition to the above, post-translational modifications also are important for plant cold responses.
Dong et al. [69] found that HOS1 can degrade ICE1 and negatively regulating CBFs and downstream cold response genes by ubiquitination under low temperature stress. In this study, the homologous gene (Gene ID: 107429115) was found to be significantly up-regulated. In addition, with the enhancement of freezing stress, the pathways of phenylpropanoid biosynthesis and flavonoid biosynthesis were also strengthened. Among the significantly DEGs, genes (Gene ID: 107429115 and 107415362) are glycosyl-transferase related genes, and the former is significantly up-regulated at -40℃, and the latter is down-regulated, which may mean that they are positively and negatively regulated factor in the process of regulating the cold resistance of jujube, which involved in different regulatory signaling pathways. Moreover, MLO proteins are a unique proteins family to plants, and the expression of many MLO genes is affected by various biotic and abiotic stresses [70]. There is a calmodulin binding region at the C-terminus of the MLO protein, which allows the MLO protein to bind were firstly processed through in-house perl scripts. In this step, clean data were obtained by removing reads containing adapter, reads containing ploy-N and low quality reads from raw data. At the same time, Q20, Q30 and GC content the clean data were calculated.  Figure S3, GO analysis of DEGs under different degree cold stress in 'Jinsixiaozao'. The X and Y axes correspond to GO terms and the number of DEGs. Figure S4, DEGs in galactose metabolism.
Different boxes represent different genes in galactose metabolism. The red boxes represent upregulated genes, and the green boxes represents down-regulated genes. Figure S5, Alternative Splicing Event and gene number at the same degree cold stress between two cultivars. Table S1, Primers for qRT-PCR. Table S1, DEGs involved in galactose metabolism pathways.      Expression profiles of the jujube AP2/ERF, WRKY, bZIP and NAC genes in two cultivars. Data were normalized based on the mean expression value of each gene in all samples analysed.
Red and green boxes indicate high and low expression levels, respectively, for each gene.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.