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

Background Low temperature is a major factor influencing the growth and development of Chinese jujube (Ziziphus jujuba Mill.) in cold winter and spring. Little is known about the molecular mechanisms enabling jujube to cope with different freezing stress conditions. To elucidate the freezing-related molecular mechanism, we conducted comparative transcriptome analysis between ‘Dongzao’ (low freezing tolerance cultivar) and ‘Jinsixiaozao’ (high freezing tolerance cultivar) using RNA-Seq. Results More than 20,000 genes were detected at chilling (4 °C) and freezing (− 10 °C, − 20 °C, − 30 °C and − 40 °C) stress between the two cultivars. The numbers of differentially expressed genes (DEGs) between the two cultivars were 1831, 2030, 1993, 1845 and 2137 under the five treatments. Functional enrichment analysis suggested that the metabolic pathway, response to stimulus and catalytic activity were significantly enriched under stronger freezing stress. Among the DEGs, nine participated in the Ca2+ signal pathway, thirty-two were identified to participate in sucrose metabolism, and others were identified to participate in the regulation of ROS, plant hormones and antifreeze proteins. In addition, important transcription factors (WRKY, AP2/ERF, NAC and bZIP) participating in freezing stress were activated under different degrees of freezing stress. Conclusions Our research first provides a more comprehensive understanding of DEGs involved in freezing stress at the transcriptome level in two Z. jujuba cultivars with different freezing tolerances. These results may help to elucidate the molecular mechanism of freezing tolerance in jujube and also provides new insights and candidate genes for genetically enhancing freezing stress tolerance.


Background
Low temperature, as an abiotic stress, not only influences the geographical distribution of many important crops but also negatively impacts the productivity and quality of crops every year [1,2]. Cold stresses are divided into chilling (0-15°C) and freezing (< 0°C) damage, which have a large impact on the survival and geographical distribution of plants [3]. Large amounts of plant species are grown in the tropics and subtropical regions, which are mostly chilling sensitive. Nevertheless, plants can increase their freezing tolerance via cold acclimation [4,5].
Chinese jujube (Ziziphus jujuba Mill.) is one of the most important fruit crops in China and has also been introduced to South Korea, the United States and many other countries [33]. In recent years, possibly because global temperature is increasing, abnormal extreme weather, such as cold extremes, has been observed to increase all over the world [34,35]. Cold stress also seriously hinders the development of the jujube industry; for example, the quality and productivity of fruits are negatively affected when exposed to low temperature. In terms of woody plants, most studies are focused on physiological and biochemical factors, and little is known about the molecular mechanisms of freezing tolerance [36,37]. Saadati et al. [37] analyzed the freezing tolerance of seven olive cultivars by measuring some biochemical and physiological parameters. Nui et al. [38] found some candidate genes potentially involved in freezing tolerance in Prunus persica through highthroughput RNA sequencing. High-throughput RNA sequencing (RNA-seq) technology is useful for the comprehensive analysis of gene functions and metabolic pathways at the transcriptional level. However, it remains unclear how the transcriptome of jujubes responds to freezing stress, and important genes involved in different freezing tolerances are rarely reported. We found that 'Jinsixiaozao' showed a much higher freezing tolerance than 'Dongzao' by LT50 values, which makes them models for investigating the mechanisms of jujube freezing tolerance (Fig. S1). Thus, we used RNA-seq technology to investigate the transcriptome difference between the two cultivars under different degrees of freezing stress. This work provides useful information for understanding the molecular regulation that occurs under freezing tolerance using two cultivars with contrasting freezing tolerance, which provides a rich genetic resource for further research investigating freezing response in jujube and other woody plants.

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 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°C, − 10°C, − 20°C, − 30°C and − 40°C, 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°C and − 40°C. 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, 107,430,353, annotated to the WRKY family, was significantly changed at lower temperatures. In the ABC family, some members (such as 107,429,561, 107,408,533, 107, 404,223, 107,414,655, and 107,416,457) 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,  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 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 Ca 2+ signaling pathway, including upregulated gene (107413690)  . 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°C, 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°C, 39 genes were upregulated and 32 were downregulated. In the case of − 20°C, 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 Ca 2+ signal pathway. Interestingly, among these DEGs, 107,422,102 was significantly downregulated at − 10°C in 'Jinsixiaozao', while it was not changed until the temperature dropped to − 40°C 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 (107,418,844, 107,432,807 and 107,422,868) 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, 107,411,641 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 (107,425,264, 107, 430,302, 107,414,011 and 107,428,511) decreased, which reduced the decomposition of stachyose (Fig. S5).

Discussion
Over the past years, extreme weather has occurred frequently, and the variations in climate have affected the growth and development of jujubes and even caused their death. Multiple responses to cold stress have been reported in many plant species, reflected by variations in transcription or metabolism [39]. Recently, significant advances have been made regarding the understanding of plant cold-signaling perception and transduction; however, we are still far from fully understanding the molecular mechanism under cold signal perception and transduction in plants [40][41][42]. In our study, the comparison of transcriptome profiles between 'Jinsixiaozao' and 'Dongzao' was reported to identify the genes responsible for freezing resistance. We first provided a more comprehensive analysis of changes in important metabolic pathways and transcription factors in jujube. Our research could also provide a reference for other woody plants and facilitate breeding for the 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 in 30 libraries of jujube branches (approximately 87 million reads in 'Dongzao' and 'Jinsixiaozao'). A total of 362 significantly DEGs (323 DEGs in 'Jinsixiaozao') were found to be differentially expressed in 'Dongzao' at different degrees of freezing stress. Compared with the control, most significant DEGs were downregulated. This result is similar to Kou et al. [43] in potato, which suggested that different cultivars have similar responses to cold resistance. However, this finding is different from that of Niu et al.'s [38] study in Prunus persica suffering from freezing tolerance, in which the number of upregulated genes were greater than that of downregulated genes. Moreover, GO classifications of 'Dongzao' and 'Jinsixiaozao' were compared and enriched in functional categories, mainly 'catalytic activity', 'metabolic process', 'cell part' and 'binding'. KEGG pathway analysis revealed that the principal enriched pathways comprised plant hormone signal transduction, metabolic pathways, ribosome and secondary metabolic pathways. With the data provided by the comprehensive analysis two contrasting Chinese jujube cultivars, more important DEGs involved in freezing stress were characterized, which can further contribute to understanding Z. jujuba.

Physiological response-related genes under freezing stress
Some research shows that plants could increase their tolerance by synthesizing numerous protective substances and proteins under these conditions [44][45][46][47]. Raffinose family oligosaccharides (RFOs) are a class of functional oligosaccharides unique to plants, mainly raffinose, stachyose and verbascose. Some studies have found that RFO can improve the ability of plants to cope with adversity [46,48]. The more RFOs accumulate, the stronger the ability of plants to be exposed to low temperatures becomes. Galactose and related enzymes are involved in the RFO metabolic pathway [48]. In our study, galactose metabolism was significantly enriched, and relative genes in this pathway were up-or downregulated to cope with cold stress. Among these genes, 107,426, 117 (stachyose synthetase), 107,415,484 (raffinose synthase) and 107,411,641 (inositol 3-alpha-galactosyltransferase) were significantly upregulated under the stronger freezing treatment in 'Jinsixiaozao', but some of them did not significantly change under the weak freezing stress. This finding may suggest these genes play a 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 the LEA gene (107409127) indicates that LEA genes are involved in freezing tolerance in 'Jinsixiaozao' (− 40°C). 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
The signal transduction pathway is known to play crucial roles in the response to freezing stress; among these, Ca 2+ is an important second messenger in plants [2,50]. The enriched Ca 2+ signaling pathway genes play critical roles in signaling pathways related to the response to freezing stress [38]. In our research, 34 DEGs were identified as related to calcium-binding proteins, calmodulin and Ca 2+ signals. Among them, the gene (107413690) showed a significant upregulation in 'Dongzao', indicating that these genes play important roles 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 [51]. On the one hand, they induce gene expression and protein synthesis to protect cells from stress; on the other hand, they induce oxidative stress [52,53]. ROS can also directly stimulate the expression of genes in activating stressrelated pathways, such as the MAPK cascade pathway, hormone signaling pathway and calcium signaling pathway. In our study, some genes with antioxidant activity were involved in freezing stress. Among these DEGs, six genes (one upregulated gene and five downregulated genes) were involved in redox signaling, such as gene 107, 403,552, which is differentially expressed under four degrees of freezing stress in 'Jinsixiaozao', indicating that they played different roles in cold signal transduction under freezing stress.

Transcriptional regulation-related genes under freezing stress
Transcription factors play various important functions in plant abiotic stress tolerance. Some studies have shown that most cold-regulated transcription factor genes are impaired in response to cold stress in Arabidopsis [54][55][56]. Many studies have shown that some transcription factors play important roles in plant freezing tolerance [57][58][59]. For example, the AP2/ERF factors RAP2.1 and RAP2.6 and the C 2 H 2 -type zinc finger STZ/ZAT10 belong to the CBF regulon [60,61]. Furthermore, NAC was also shown to enhance expression during loquat fruit chilling; however, it could not bind directly to promoters of lignin biosynthesis genes [62]. Many MYBtype transcription factors, such as AtMYB85 [63] and AtMYB88 [64,65], increase cold hardiness by CBFdependent and CBF-independent pathways in apple. In other species, OsMYB4 transgenic lines of Arabidopsis exhibited improved chilling and freezing tolerance with a dwarf phenotype [66]. In this study, the expression patterns of the AP2/ERF, WRKY, NAC and bZIP gene families were analyzed, and we found that different members play different roles in the two cultivars. For example, in the AP2/ERF family, different members of the family are involved in the response to low temperatures. Among these genes, 107,414,104 expression showed significant changes at − 10°C; genes 107,432,697 and 107, 433,390 were at − 20°C; 107,426,263 and 107,433,390 were at − 30°C; however, 107,414,104, 107,432,697 and 107,433,390 changed significantly at the lowest temperatures. Similarly, other TF members are generally divided into two expression patterns. These differences suggest that the two cultivars have different responses to cold stress, even for different degrees of freezing stress. 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 regulators 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 were exposed to different stresses, especially low temperatures [67]. To analyze the variation in AS under cold stress, in this study, we detected that all genes containing introns were alternately spliced under cold stress (Fig. S6). On the one hand, compared with the control, the number of AS events increases slightly as the temperature drops, and 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, post-translational modifications are also important for plant cold responses. Dong et al. [68] found that HOS1 can degrade ICE1 and negatively regulate CBFs and downstream cold response genes by ubiquitination under low temperature stress. In this study, the homologous gene (107429115) was found to be significantly upregulated. In addition, with the enhancement of freezing stress, the pathways of phenylpropanoid biosynthesis and flavonoid biosynthesis were also strengthened. Among the significantly DEGs, genes (107,429,115 and 107,415,362) are glycosyl-transferase-related genes, and the former is significantly upregulated at − 40°C, and the latter is downregulated, which may mean that they are positively and negatively regulated factors in the process of regulating the cold resistance of jujube, which are involved in different regulatory signaling pathways. Moreover, MLO proteins are a unique protein family in plants, and the expression of many MLO genes is affected by various biotic and abiotic stresses [69]. There is a calmodulin binding region at the C-terminus of the MLO protein, which allows the MLO protein to bind to calmodulin as a second messenger to transmit signals to the downstream pathway. The CaM protein binds to the CaMBD of the MLO protein to increase the activity of the MLO protein. In our data, we found that the gene 107,409,553 encoding the MLO protein showed a significant upregulation when 'Jinsixiaozao' was at − 40°C but showed downregulation in 'Dongzao'. This finding may indicate that the genes regulated by calcium are varietally different between the two cultivars in response to freezing stress.
Some studies have shown that cold stress can restrict plant growth and development at different levels, for example the cell and tissue levels. Roots and branches are both the tissues that predominantly perceive and transduce cold signal in cold stress responses [38,40], we comprehensively analyzed the jujube branches at the tissue level which provided valuable information and genetic resources for improving cold stress in jujube. In addition, we found some significant DEGs under higher degrees of freezing stress, which may suggest that the two cultivars have different responses to different degrees of freezing stress in the metabolic pathway and biosynthesis of galactose. Jujubes sense freezing stress through membrane rigidification and cellular changes, leading to an influx of Ca 2+ , reactive oxygen species (ROS) and phytohormone changed. The signaling pathways activate protein kinases and changed various downstream transcription factors family proteins, such as WRKY, AP2/ERF, NAC and bZIP. The activation of TFs triggers the expression of downstream responsive genes to the stress treatment, as well as the induction of biochemical responses leading to freezing tolerance. These results may provide evidence to support that 'Jinsixiaozao' has much higher freezing tolerance than 'Dongzao'. Our understanding of the regulatory mechanisms may provide a reference for further research investigating the role of different tissues in freezing stress and other abiotic stresses.

Conclusions
Our present study first provides a comprehensive analysis via transcriptome sequencing between two Chinese jujube cultivars ('Dongzao' and 'Jinsixiaozao') under freezing stress. The important genes involved in signal transduction and transcriptional and posttranslational regulation were differentially expressed between these two cultivars. Some different members of the AP2/ERF, WRKY, NAC and bZIP gene families were significantly up/downregulated in the two cultivars. In addition, some DEGs, such as 107,426,117 and 107,411,641, involved in galactose metabolism were significantly upregulated under higher degree freezing treatment in 'Jinsixiaozao' than in 'Dongzao', which further highlights the importance of sugar metabolism involved in freezing tolerance in plants. The present results may explain why 'Jinsixiaozao' has considerably higher freezing tolerance than 'Dongzao' and may elucidate the molecular mechanisms of freezing tolerance. The important DEGs discovered in this study may facilitate further study on freezing tolerance mechanisms in jujube and other woody plants.

Plant materials and cold treatment
The current year's branches of cold-sensitive cultivar 'Dongzao' (DZ) and cold-tolerant cultivar 'Jinsixiaozao' (JS) were collected from the National Key Base for Improved Chinese Jujube Cultivar (Cangzhou, China). The collected branches were placed at 4°C for 10 h and used as controls (DZCK and JSCK), the others were treated and maintain 10 h at the freezing temperature: -10°C (DZA and JSA), − 20°C (DZB and JSB), − 30°C (DZC and JSC) and − 40°C (DZD and JSD). Then the xylem of the branches was collected and snap-frozen in liquid nitrogen. Three replicates were used for each treatment.

RNA-seq library construction and sequencing
Total RNA were extracted from the 30 samples using the Plant RNA Kit (Omega) according to the manufacturer's instructions. Three biologic replicates for each treatment were used for analyzing transcriptome. The library products were ready for sequencing analysis using an Illumina Genome Analyzer at Beijing Genomics Institute (Shenzhen, China) in 2018. Raw data of fastq format were firstly processed through in-house perl scripts to obtain clean data by removing the reads containing adapter, ploy-N and low-quality reads. All the subsequent analyses were based on the clean data. Reference genome and gene model annotation files were based on Ziziphus jujuba (assembly ZizJuj_1.1) on NCBI. The Zscore method using the p-value as a statistical significance index [70] was applied to identify differentially expressed genes.

Identification of DEGs, GO and KEGG enrichment analysis
A cluster analysis was performed according to Eisen et al. [71]: the log2 of TPM for each gene was used for the hierarchical clustering analysis. Gene Ontology (GO) categorization was classified into molecular function, biological process, and cellular component. We used KOBAS software to test the statistical enrichment of differential expression genes in KEGG pathways (http:// www.genome.jp/kegg/). Four types of AS events including exon skipping, intron retention, alternative 5′ splice site and alternative 3′ splice site were identified using TopHat [72].

Real-time quantitative PCR
Six randomly selected DEGs were chosen for validation by real-time quantitative PCR (qRT-PCR). The sequences of the primer pairs (designed using Beacon Designer 7.2) were listed in Table S1. Realtime PCR was performed in the presence of TB Green Premix Ex Taq II (Tli RNaseH Plus) (TaKaRa, Dalian, China) and read on an IQ5 real-time PCR instrument (Bio-Rad, Hercules, CA, USA). The amplification was performed as previously reported [73], and all the experiments were performed in biological triplicate. Transcript levels were normalized against the average expression of the Zjactin gene (GenBank: EU916201).