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 library reads, and each library had approximately 89 million reads. After filtering, clean reads varied from 87 million to 88 million, which occupied on average 97.47% of the total reads of all libraries (Table 1). In all samples, the Q30 value was over 86%, and the GC content was 43%. The genome map rate accounted for more than 73%, and the gene map rate reached 72.93% to 76.32%. A total of approximately 20,000 genes were detected to be expressed in each library. In addition, more novel transcripts were detected in both DZD and JSD. Then, regarding alternative splicing, each sample had a similar number of events. Indel, JSD has the maximum mark, significantly more than in DZC. Pearson’s correlation coefficients among the three biological replicates were high (γ=0.91-0.98).
DEGs obtained by comparing two cultivars
To understand the differences between different cultivars in response to freezing stress, gene expression profiles in ‘Dongzao’ and ‘Jinsixiaozao’ at the same level of freezing stress were further analyzed though FPKM values. The detailed statistics are shown in Additional file 2. We investigated DEGs specific to each freezing stress by comparing ‘Dongzao’ to ‘Jinsixiaozao’. There were 1831 (885 upregulated genes and 946 downregulated genes), 2030 (979 upregulated genes and 1051 downregulated genes), 1993 (949 upregulated and 1044 downregulated genes), 1845 (923 upregulated and 922 downregulated genes) and 2137 (1158 upregulated and 979 downregulated genes) DEGs were identified at 4℃, -10℃, -20℃, -30℃ and -40℃, respectively, suggesting the responsiveness of these genes to freezing and cold resistance (Fig.1 and Additional file 3).
In five different degrees of cold stress, as shown in Fig. 2, 854 common DEGs were detected in both cultivars, which may suggest that these DEGs are involved in the process of dealing with low temperature under different low temperatures, and both participate in the same pathways in response to freezing stress. In addition, at -40°C, the number of unique DEGs (600) was significantly higher than that in the other groups, indicating that these genes may be related to the difference in cold tolerance between ‘Dongzao’ and ‘Jinsixiaozao’.
In this study, total RNA from branches was used for qRT-PCR validation. As shown in Fig. S2, six DEGs showed similar expression pattern between the qRT-PCR data and the RNA-seq results, demonstrating that the RNA-seq data is highly reliable for further analysis.
Function annotation and enrichment analysis of DEGs
To understand the functions of the discovered DEGs between the two cultivars, all DEGs were searched against the protein database, followed by Gene Ontology (GO) analysis to assess the putative functions of the genes. Among the GO terms of each treatment, the stronger freezing stress had a larger proportion of DEGs, and the more metabolic processes were modified and changed (Fig. 3). In the comparative analysis of five groups, the DEGs of these GO terms 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 elucidate the metabolic pathways involved in freezing stress, 20 common metabolic processes were most enriched at the five freezing treatments in the two cultivars (Fig. 4). Among these top pathways, the carbohydrate metabolism pathway, galactose metabolism, fructose and mannose metabolism were associated with many up/downregulated genes. Furthermore, in the comparison of the five groups, there was some pathway enrichment in the stronger freezing stress in ‘Jinsixiaozao’, such as peroxide of metabolic pathways, which was only significantly enriched at -30°C, and the ABC transporter process was enriched at -40°C. In brief, 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 great significance for revealing the cold resistance of jujube.
To further analyze the difference in jujube at extreme temperatures, we conducted a comparative analysis of the two cultivars with temperatures of -30℃ and -40℃. In comparison, we found that more significantly DEGs were found in the stronger freezing stress. The peroxisome pathway was enriched in DZC vs. JSC, and the ABC transporters pathway was enriched in DZD vs. JSD. In addition, 107430353, annotated to the WRKY family, was significantly changed at lower temperatures. In the ABC family, some members (such as 107429561, 107408533, 107404223, 107414655, and 107416457) simultaneously exhibited significant upregulation. The ratio of DEGs related to kinase activity and carbohydrate metabolism also increased under the stronger freezing stress than under weak freezing stress.
Comparison analysis at different degrees of freezing stress in ‘Dongzao’ and ‘Jinsixiaozao’
To understand the different effects on jujube at different degrees of freezing stress, we also compared four freezing stresses in ‘Dongzao’ and ‘Jinsixiaozao’. In ‘Dongzao’, DZCK was used as a control, and 1291 DEGs were identified at different degrees of freezing stress, although only 55 common DEGs were identified in all freezing stresses (Fig. 5, Fig. S3 and Additional file 4). Then, these DEGs were further screened and 327 significant DEGs (PPEE<0.05 and |Log2|≥1) were identified. Of these significant DEGs, some genes were related to the Ca2+ signaling pathway, including upregulated gene (107413690) and downregulated genes (107404917, 107423107, 107425924, 107417689, 107421107, 107422102, 107430853 and 107424939). Some DEGs related to sucrose metabolism were classified into biological processes among the four groups, including 7 upregulated and 25 downregulated. Furthermore, regarding transcription factors, through the analysis of DEGs, we found that the genes belonging to the NAC, WRKY, MYB, AP2/ERF and bHLH families all showed a significant up- or downregulation. A gene (107424282) homologous to the grape MYB transcription factor is significantly downregulated under cold stress, suggesting that it may negatively regulate the cold resistance of jujube trees. In addition, the NAC transcription factor (107435293) had a higher fold change at -30℃, indicating that this transcription factor may regulate the cold resistance of jujube trees in response to freezing stress.
In the comparison of four groups of ‘Jinsixiaozao’, JSCK as a control, a total of 1668 DEGs were identified (Fig. 5, Fig. S4 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 upregulated and 32 were downregulated. In the case of -20℃, 108 DEGs were identified, including 31 upregulated genes and 77 downregulated genes. With decreasing temperature, 31 genes were upregulated and 70 genes were downregulated at, -30°C. At -40°C, only 43 genes had larger changes. Similar to the results of ‘Dongzao’, many genes associated with abiotic stress have been annotated in ‘Jinsixiaozao’, such as sucrose metabolism and the Ca2+ signal pathway. Interestingly, among these DEGs, 107422102 was significantly downregulated at -10℃ in ‘Jinsixiaozao’, while it was not changed until the temperature dropped to -40℃ in ‘Dongzao’. This finding may be due to the different regulatory mechanisms of the two cultivars in response to freezing stress, thereby reflecting different cold resistance.
Key transcription factors associated with freezing stress
In this study, we further analyzed the expression of key transcription factors (AP2/ERF, WRKY, NAC and bZIP), which have been reported to be involved in freezing stress. All FPKM values of these genes are listed in Additional file 6. In all libraries, we identified a total of 28 WRKY members, 52 AP2/ERF members, 22 NAC members and 15 bZIP members in two cultivars (Fig. 6). Comparison of the expression of these TFs between the two cultivars identified WRKY members; some members were highly expressed in ‘Dongzao’, and others were highly expressed in ‘Jinsixiaozao’. Among the TFs, the expression levels of WRKY04 (107417519), WRKY11 (107414569), and WRKY20 (107419756) in ‘Jinsixiaozao’ were completely different from those in ‘Dongzao’. In the AP2/ERF gene family, three of 52 AP2/ERF members (107418844, 107432807 and 107422868) were upregulated under light cold stress in ‘Dongzao’ but began to change under the lower temperature ‘Jinsixiaozao’. In terms of NAC and bZIP, NAC19 had a higher expression level in DZD than in other samples, while in ‘Jinsixiaozao’, the expression level was higher at all five temperatures; bZIP members are generally divided into two expression patterns. The observed expression differences of these transcription factors between the two cultivars suggest that TFs may play important roles in the response to different freezing stresses in jujube.
Key DEGs in the galactose metabolism pathway
The accumulation of stachyose plays an important role in plant resistance to low temperature stress, and α-galactosidase is important in the decomposition process of stachyose. In the comparison of the five groups, the galactose metabolic pathway was significantly enriched under each temperature condition. Thus, we further characterized the identified DEGs involved in the metabolism of galactose. During the metabolism of galactose, some genes associated with the galactose, stachyose, raffinose, and sucrose pathways showed significant differences in coping with freezing stress. As listed in Table S2, among these genes, 107411641 is a key gene in the process of UDP-galactose, which is significantly upregulated under stronger freezing stress, which will enhance the synthesis pathway of stachyose. However, in the pathway of stachyose decomposition into galactose and raffinose, the FPKM value of genes (107425264, 107430302, 107414011 and 107428511) decreased, which reduced the decomposition of stachyose (Fig. S5).