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.
Table 1
Summary of the RNA-seq data collected from ‘Dongzao’ and ‘Jinsixiaozao’.
Sample | Raw data | Clean reads | Q30(%) | GC(%) | Genome map Rate | Gene map Rate | Expressed Gene | Expressed Transcripts | Novel Transcripts | Alternative Splicing | Indel |
DZCK | 89647985 | 87382161 | 86.91 | 44.19 | 72.64% | 74.96% | 20411 | 25818 | 988 | 90011 | 4342 |
DZA | 89648769 | 87565991 | 88.23 | 44.12 | 73.14% | 75.11% | 20362 | 25710 | 954 | 87733 | 8628 |
DZB | 89648625 | 87947768 | 87.87 | 44.33 | 73.84% | 76.32% | 20431 | 25881 | 980 | 88983 | 5811 |
DZC | 89648887 | 87910197 | 88.73 | 44.19 | 73.52% | 75.50% | 20192 | 25458 | 924 | 85154 | 3927 |
DZD | 89649519 | 87899251 | 88.22 | 44.02 | 73.98% | 73.91% | 21546 | 26853 | 1066 | 87281 | 4230 |
JSCK | 89647603 | 87753183 | 87.40 | 43.64 | 73.12% | 74.68% | 20416 | 25724 | 983 | 89086 | 5555 |
JSA | 89648885 | 87822119 | 89.13 | 43.69 | 72.43% | 74.23% | 20264 | 25542 | 983 | 91659 | 4705 |
JSB | 89648594 | 88010731 | 88.92 | 43.91 | 72.85% | 74.63% | 20154 | 25425 | 972 | 86901 | 6013 |
JSC | 89648221 | 87657475 | 87.74 | 43.70 | 72.58% | 73.31% | 19871 | 25011 | 963 | 87767 | 5261 |
JSD | 89647797 | 87561609 | 87.25 | 43.69 | 73.28% | 72.93% | 20648 | 26003 | 1051 | 99095 | 6220 |
DZCK, ‘Dongzao’ treated at 4℃; DZA: ‘Dongzao’ treated at -10℃; DZB: ‘Dongzao’ treated at -20℃; DZC: ‘Dongzao’ treated at -30℃; DZD: ‘Dongzao’ treated at -40℃; JSCK: ‘Jinsixiaozao’ treated at 4℃; JSA: ‘Jinsixiaozao’ treated at -10℃; JSB: ‘Jinsixiaozao’ treated at -20℃; JSC: ‘Jinsixiaozao’ treated at -30℃; JSD: ‘Jinsixiaozao’ treated at -40℃. |
DEGs obtained by comparing two cultivars
In order to understand the differences between different cultivars in response to different low temperature stress, gene expression profiles in ‘Dongzao’ and ‘Jinsixiaozao’ at the same freezing stress were further analyzed though FPKM values. The detailed statistic is shown in Additional file 2. First, we investigated DEGs specific to each freezing stress by comparing ‘Dongzao’ to ‘Jinsixiaozao’. There were 1831 (885 up-regulated genes and 946 down-regulated genes), 2030 (979 up-regulated genes and 1051 down-regulated genes), 1993 (949 up-regulated and 1044 down-regulated genes), 1845(923 up-regulated and 922 down-regulated genes) and 2137(1158 up-regulated and 979 down-regulated genes) DEGs were identified at 4℃, -10℃, -20℃, -30℃ and − 40℃, respectively, suggesting responsive 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 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 stress. Peroxisome pathway were enriched in DZC vs. JSC, and ABC transporters pathway enriched in DZD vs. JSD. In addition, 107430353 annotated to the WRKY family was significantly changed at lower temperature, and in ABC family Genes 107429561, 107408533, 107404223, 107414655, and 107416457 simultaneously exhibited significant up-regulation. The number of DEGs related to kinase activity and carbohydrate metabolism were increased at the stronger freezing stress than the weak.
KEGG pathway classifcation of DEGs under cold stress in ‘Dongzao’ and ‘Jinsixiaozao’, respectively. The Y axis corresponds to KEGG pathway, the X axis shows the enrichment ratio between the number of DEGs enriched in a particular pathway. The color of the dot represents q-value, and the size of the dot represents the number of DEGs mapped to the referent pathway.
Comparison analysis at different degree freezing stress in ‘Dongzao’ or ‘Jinsixiaozao’
To further understand the effect of different degrees freezing stress on jujube, we also compared analysis five stresses in ‘Dongzao’ and ‘Jinsixiaozao’, respectively. In ‘Dongzao’, DZCK as a control, 1291 DEGs were identified at the different degree of freezing stress, though only 55 common DEGs were identified in all freezing stress (Fig. 5, Fig. S2 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 related to Ca2+ signal pathway, including up-regulated gene (107413690) and down-regulated genes (107404917, 107423107, 107425924, 107417689, 107421107, 107422102, 107430853 and 107424939). some related to sucrose metabolism were classified into biological process among the four groups, which including 7 up-regulated and 25 down-regulated. Furthermore, about transcription factor, 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 down changed. 107424282 homologous with grape MYB transcription factor is significantly down-regulated in cold stress, suggesting that it may negatively regulate cold resistance of jujube trees. In addition, NAC transcription factor (107435293) had a higher fold changes at -30℃, indicating that this transcription factor may be regulating the cold resistance of jujube trees in response to low temperature.
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-regulated genes. With the decrease of temperature, 31 genes were up-regulated and 70 genes were down-regulated at -30℃. At -40℃, only 43 genes had bigger changed. Similar to the results of ‘Dongzao’, many genes associated with abiotic stress have been annotated in ‘Jinsixiaozao’, such as sucrose metabolism, Ca2+ signal pathway and so on. Interestingly, among these DEGs, 107422102 was significantly down-regulated at -10℃ in ‘Jinsixiaozao’, however, while it was not changed until the temperature dropped to -40℃ in ‘Dongzao’. That may the different regulatory mechanisms of the two cultivars in response to freezing stress, thus reflecting different cold resistance.
Key transcription factors associated with freezing stress
A number of genes are induced by these stresses, and their products are thought to function in stress tolerance and response. 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 response to freezing stress. All samples FPKM values of these genes are list in Additional files 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 these TFs expressions between the two cultivars revealed that all identified WRKY members, some are highly expressed in ‘Dongzao’ and others are highly expressed in ‘Jinsixiaozao’, indicating different regulation in different pathway. Among them, the expression level of WRKY04 (Gene ID: 107417519), WRKY11 (Gene ID: 107414569), WRKY20 (Gene ID: 107419756) in ‘Jinsixiaozao’ have totally different from ‘Dongzao’. In AP2/ERF gene family, three of 52 AP2/ERF members (Gene ID: 107418844, 107432807 and 107422868) were up-regulated under the light cold stress in ‘Dongzao’, but began to change under the lower temperature ‘Jinsixiaozao’. In terms of NAC and bZIP, NAC19 (Gene ID: 107415997) has higher expression level in DZD than 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 response to different freezing stress in jujube.
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.
Key DEGs in galactose metabolism pathway
The accumulation of stachyose plays an important role in plant resistance to low temperature stress, while galactionl is a key enzyme in the synthesis of stachyose, and α-galactosidase is important in the decomposition process of stachyose. Thus, we further characterized the identified DEGs involved in the metabolism of galactose. In the comparison of five groups, KEGG pathway analysis that the galactose metabolic pathway was significantly enriched under each temperature condition. During the metabolism of galactose, some genes associated with galactose, stachyose, raffinose, and sucrose 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).