- Research article
- Open Access
Comprehensive profiling of alternative splicing landscape during cold acclimation in tea plant
BMC Genomics volume 21, Article number: 65 (2020)
Alternative splicing (AS) may generate multiple mRNA splicing isoforms from a single mRNA precursor using different splicing sites, leading to enhanced diversity of transcripts and proteins. AS has been implicated in cold acclimation by affecting gene expression in various ways, yet little information is known about how AS influences cold responses in tea plant (Camellia sinensis).
In this study, the AS transcriptional landscape was characterized in the tea plant genome using high-throughput RNA-seq during cold acclimation. We found that more than 41% (14,103) of genes underwent AS events. We summarize the possible existence of 11 types of AS events, including the four common types of intron retention (IR), exon skipping (ES), alternative 5′ splice site (A5SS), and alternative 3′ splice site (A3SS); of these, IR was the major type in all samples. The number of AS events increased rapidly during cold treatment, but decreased significantly following de-acclimation (DA). It is notable that the number of differential AS genes gradually increased during cold acclimation, and these genes were enriched in pathways relating to oxidoreductase activity and sugar metabolism during acclimation and de-acclimation. Remarkably, the AS isoforms of bHLH transcription factors showed higher expression levels than their full-length ones during cold acclimation. Interestingly, the expression pattern of some AS transcripts of raffinose and sucrose synthase genes were significantly correlated with sugar contents.
Our findings demonstrated that changes in AS numbers and transcript expression may contribute to rapid changes in gene expression and metabolite profile during cold acclimation, suggesting that AS events play an important regulatory role in response to cold acclimation in tea plant.
Low temperature is one of the most important environmental factors affecting plant growth, development and geographical distribution. Cold acclimation (CA) is an important mechanism that has been widely reported to improve cold resistance of plants by modulating numerous physiological and biochemical processes [1,2,3]. For instance, cold acclimation improves the tolerance of North American Rhododendron from − 7 °C to − 53 °C . During CA, the contents of soluble sugars including sucrose, fructose, glucose and raffinose increase significantly, and these sugars are thought to osmotically stabilize membranes [5, 6]. Unsaturated fatty acids also play an important role in cold resistance, which has previously been demonstrated in tobacco . In addition to increases of these substances, many transcription factors and oxidoreductase regulatory mechanisms are also present in plants during cold acclimation. For example, The CBF (C-repeat-binding factor) transcription factor regulates COR (cold-regulated) genes which function in the response to low temperature . Plants can also defend against active oxygen through the protection system of oxidoreductases during cold acclimation . Cold acclimation grants the ability to withstand low temperature and plays an important role in the growth and development of plants.
Alternative splicing (AS) is the process by which different mRNA splicing isoforms are produced from a single mRNA precursor by different splicing methods.AS can generate multiple transcript and protein isoforms from the same gene [10, 11]. Evidence is accumulating that AS plays a crucial role in a variety of plant development processes and stress responses, including regulation of flowering , defense response to pathogenic bacteria , and abiotic stress response [14, 15]. The targets of AS in rice were found to include three Serine/arginine-rich (SR) protein-encoding genes regulating phosphorus (P) uptake and remobilization in a highly nutrient-specific manner, thus demonstrating that AS plays a critical role in maintaining mineral nutrient homeostasis . Previous studies have reported that temperature-associated alternative splicing is an important mechanism involved in the regulation of flowering and the plant circadian clock [12, 17]. Under cold stress, hundreds of genes in A.thaliana showed transcriptional changes due to rapidly occurring AS and AS can affect several cold-responsive transcription factors and RNA binding proteins [18, 19]. In durum wheat, low temperature promotes intron retention in two early cold-regulated genes . In temperature signaling cascades in plants, AS is considered as a way of perceiving temperature fluctuations and modulating transcription factor activity, perhaps by linking regulation of gene expression with PEPi (Peptide interference) and/or NMD (Nonsense-mediated decay) mechanisms . However, most of these studies focus on the low temperature, and few address on plant cold acclimation, especially the AS changes after de-acclimation.
The tea plant (Camellia sinensis) is a perennial evergreen woody crop and its leaves are used for making tea beverage, which is one of the three popular nonalcoholic beverages consumed worldwide [22, 23]. It is popular among most consumers due to its good taste and health-promoting effects [24, 25]. Tea plants are vulnerable to low temperatures in winter, especially in northern China . In tea plant, the contents of sucrose, glucose and fructose were found to be constantly elevated during cold acclimation . In addition, significantly increased expression of CBF (C-repeat-binding factor) and DHN (dehydrin) occurs during cold acclimation . Both increases in sugar contents and expression of related genes can improve the cold resistance in tea plant. But these studies have not focused on the effect of AS on cold stress in tea plants. In our previous study, we found that some AS events were tissue specific in stem and root; we also found that some AS isoforms were the major transcripts involved in the flavonoid synthesis pathway, which suggested AS is positively correlated with the contents of catechins . We also reported that six CsLOXs (lipoxygenase) varied significantly in relative abundances under the different stresses . However, these studies on AS emphasize largely on secondary metabolism.
In this study, we investigated AS events during cold acclimation with genome-wide analysis in tea plant, which detected a large number of AS occurrences. Meanwhile, variations of AS and the related biological functions were analyzed during cold acclimation. Our results indicated that AS may regulate gene expression and contents of metabolites during cold acclimation. A number of different AS patterns were found being involved in transcriptional regulation during the process of cold acclimation in tea plants, especially the change and function of AS at de-acclimation. This provides a better understanding of the functions of AS in tea plants responding to cold acclimation.
Global identification and classification of AS events
The RNA-seq data was used to investigate AS events at different periods of cold acclimation (Additional file 4: Table S4). Briefly, clean reads were initially mapped to the tea plant reference genome, and AS events were then identified using the AStalavista tool. From this analysis, a total of 63,329 AS events were identified from 14,103 genes in all samples (Fig. 1). Among them, IR was the most abundant type (18,231, 28.79%), followed by A3SS (9561, 15.10%), A5SS (6465, 10.21%) and ES (6187, 9.77%) (Fig. 1). In addition, a large number of other types of AS were detected in all samples because of multiple splicing modes occurring on a single transcript. These results were consistent with those of previous reports in other plant species . During the sharp cooling treatment of tea plants, the numbers of AS events under low-temperatures [cold acclimation of 6 h at 10 °C, day/night temperature (CS), cold acclimation of 7 days at 10/4 °C, day/night temperature (CA1) and cold acclimation of 7 days at 4/0 °C, day/night temperature (CA2)] were significantly increased relative to AS events in those growing under normal temperature (NA, 25/20 °C, day/night temperature); thus, the amount of CS (30,212), CA1 (30,325) and CA2 (30,552) was about 3000 more than NA (27,234), respectively. It is noteworthy that the number of AS events (24,616) tremendously decreased under temperature recovery condition compared with both CA treatment and NA groups, implying that many AS genes had adapted to low-temperature environment and resulted in reduced occurrence of AS events during DA (de-acclimation) of tea plant.
Considering the isoforms in the annotated loci, most of the splicing junctions (SJs) were found to reside in the coding regions (or coding sequence) (CDS) (137,112), while fewer SJs were observed in 5′ untranslated region (5′UTR) (3143) and 3′ untranslated region (3′UTR) (3030), highlighting that coding proteins were intensively influenced by AS (Fig. 2a). At the splice event level, the number of known splice annotation was the highest followed by partial novel. Whereas, the novel splice annotation accounted for the greatest number at the splice junction level (Fig. 2b). With regard to the splicing donor-acceptor sites, a majority (96.50%) were found to be the canonical GU-AG, followed by GC-AG with only 1.68% (Fig. 2c and Additional file 1: Table S1). The noncanonical AU-AC splice site pair, specifically (U12-type introns) are thought to have important regulatory roles , and accounted for only 0.53% (Fig. 2c). The 5′ss and 3′ss sequences are shown in Fig. 2d. We located the branch site and branch point A for almost all of the SJs (Fig. 2e). The average distance from the branch point A to 3′ss was 54.5 ± 86.6 bp in length. Most introns were found to be 51 bp - 150 bp in length, the average intron length was 2527.04 ± 12,234.5 bp, and the median was 500 bp. The branch point A offset from the 5’ss was positively correlated with intron length (Fig. 2f).
Among 14,103 identified AS genes at five different time points, 3779 AS genes were conserved during cold acclimation (Fig. 3a). Some genes were also found to exhibit AS specificity; for example, 1204 DA-specific AS genes were found in DA (de-acclimation). Interestingly, the amount of DA-specific AS genes were significantly more than the number of NA- (588), CS- (642), CA1- (484) and CA2-specific genes (694), suggesting the stimulation of abundant novel AS genes by cold acclimation. To explore the biological functions influenced by AS, KEGG enrichment analyses of 14,103 AS genes were performed (Fig. 3b). The result indicated that the process of spliceosome, RNA transport and RNA degradation were highly enriched, which were in accordance with a previous study .
Dynamic characterization and analysis of AS events during cold acclimation
To investigate variation in types of AS events during cold acclimation, we calculated the proportion of the four main AS events (IR, A5SS, A3SS and ES) among the five periods. However, the differences observed among these AS events were insignificant (Additional file 10: Figure S2). We then counted the number of differentially AS genes (DAGs) at different time points (Fig. 4) based on TPM value. It was demonstrated that the number of DAGs gradually increased during cold acclimation (Fig. 4a), with a total of 940, 1552, 2264 and 2025 DAGs observed in the NA vs CS, NA vs CA1, NA vs CA2 and NA vs DA groups, respectively. Notably, the number of unique DAGs (1142) in NA vs DA group was far greater than in other groups, which implied the occurrences of novel specific AS genes in response to DA.
To further determine the biological functions of tea plant AS genes involved in the two transition periods (normal temperature to low temperature and de-acclimation to non-acclimation), we performed GO enrichment analyses of DAGs in these two time points. The highest enrichment of 15 molecular functions is shown according to the p-value (p < 0.01) (Fig. 4b and c). When the tea plants were initially subjected to cold stress, the DAGs were mainly enriched on the activity of oxidoreductase and the binding of substances. Interestingly, pathways relating to sugar, serine and oxidoreductase activity were significantly enriched after cold acclimation compared with non-acclimation.
Characterization and expression pattern of AS genes during cold acclimation
To better understand the putative impact of AS during cold acclimation, AS genes associated with cold stress were investigated. We explored and identified genes involved in regulatory and metabolic pathways, including the ICE-CBF-COR pathway, proline metabolism, sugar metabolism pathway, oxidoreductase, abscisic acid pathway these are all involved in the defense response to cold stress. Expression levels of many transcription factors including CsMYBs, CsbHLHs and CsWRKYs were changed at low temperature. (Additional file 2: Table S2). Overall, about 35.9% (42/117) AS events were dramatically enriched in genes within the sugar metabolism pathway, while only 18.6% (24/129) AS events were observed in oxidoreductase clusters (Additional file 2: Table S2). The expression patterns of these AS transcripts were then analyzed (Additional file 11: Figure S3).
To verify the accuracy of the high throughput RNA-seq data, some AS transcripts were selected for validation using RT-PCR and with primers flanking the AS site (Fig. 5 Additional file 3: Table S3). All amplified PCR products were sequenced and verified (Additional file 3 Table S3) which were in consistent with the results of gels. These AS events mainly presented as IR, A3SS and complex types. Among them, IR was the dominant type, and sequence analysis revealed that most IR type AS transcripts could result in introduction of a premature stop codon (PTC); these transcripts may be translated into truncated proteins, thereby imparting structural and functional diversity between AS and non-AS transcripts. For instance, four different AS transcripts were found in SUS (sucrose synthase) gene, its full-length transcript (genome annotated transcripts) encodes a protein possessing 838aa and a complete sucrose synthase domain which potentially exhibits catalytic activity. Nevertheless, a PTC was observed in one SUS AS transcript (CsSUS-3) which encodes only 206 aa of the protein, implying that this remarkable variation may significantly affect the structure and function of this protein (Fig. 5). Additionally, similar truncated proteins were found in CsCOR (cold-regulated gene), CsRS (raffinose synthase gene) and CsPOD (peroxidase gene).
To further explore the biological function of these AS transcripts during cold acclimation, their expression patterns were determined using qRT-PCR. The expression levels of AS transcripts vary significantly among different time periods (Fig. 6). For example, the expression of full-length CsPOD2 increased under low temperature while its AS isoform showed an opposite expression pattern; expression of the CsRS-2 isoform was significantly higher during CS and CA2 than the other time periods. Meanwhile, CsbHLH-1 expressed in CS and CA1 but not the other time periods. The AS isoform CsbHLH-2 was the predominant transcript and expression was higher under low temperature (Fig. 6).
The relationship between AS transcripts and metabolites during cold acclimation
To explore the relationship between AS transcripts and metabolites, the metabolome data was initially obtained by LC–MS using the same treated leaf materials with six independent biological replicates. A total of 19,305 substances were identified, including 8344 annotated metabolites another 590 metabolites which were further identified by LC-MS/MS (Fig. 7a). Furthermore, 5062 metabolites with > 2-fold change and p-value < 0.05 were defined as differential metabolites; of these, we identified 2911, 3309 and 2478 differential metabolites in the NA vs CA2, CA2 vs DA and NA vs DA groups, respectively (Fig. 7b). In addition, a total of 137 differential metabolites were finally identified and classified by LC-MS/MS, of which the majorities were flavonoids, sugars and fatty acids. We analyzed the changes of these three metabolite categories together with glucose, fructose, sucrose and raffinose measured in the previous study  during cold acclimation (Fig. 7c).
We found that a large number of AS variants were involved in sugar metabolism (Fig. 8a), and sugars were also among the metabolites with differential abundances. Therefore, correlation analysis was performed between these AS transcripts and the contents of sugars. The results showed differences in the correlation between sugar content and AS transcripts (Fig. 8 and Additional file 12: Figure S4). For example, four transcripts of CsSUS (TEA025243) were positively correlated with sucrose content; among them, three transcripts (TEA025243.2, TEA025243.3 and TEA025243.4) were positively correlated with glucose and fructose content (Fig. 8c). Moreover, the expression levels of two SPS transcripts (TEA030791) showed opposite trends with sucrose content (Fig. 8c). Interestingly, the full-length CsRS-1 (TEA013682.1) is strongly correlated with raffinose content, while its AS isoform CsRS-2 (TEA013682.2) is correlated with levels of glucose, sucrose and fructose rather than raffinose (Fig. 8b). We also analyzed the domains of CsSUS and CsRS transcripts and found that the domains of TEA025243.3, TEA025243.4 and TEA013682.2 were incomplete due to the premature termination codon. What is more interesting is that the TEA025243.2 transcript of CsSUS gene possessed a transmembrane region (Additional file 13: Figure S5).
Alternative splicing (AS) plays an important role in plant responses to biotic and abiotic stresses such as cold stress and pathogen defense [13, 33]. In tea plant, previous studies indicate AS is involved in important secondary metabolism pathways such as flavonoid and linalool biosynthesis [29, 34]. Cold stress is one of the most severe abiotic stresses in tea plant, and a series of studies have been focusing on tea plant gene regulation during cold stress and cold acclimation [27, 35, 36]. Through cold acclimation, tea plants can alleviate the harm of low temperature and enhance their tolerance to cold stress . However, little information is known about AS in tea plant during cold acclimation. In this study, we investigated the total AS events during cold acclimation in tea plant through high-throughput RNA-seq. The numbers and expression levels of AS events varied greatly during cold acclimation, and some AS transcripts were strongly correlated with sugar content (correlation coefficient > 0.7). This study provides novel understanding of the roles of AS for cold acclimation in tea plant.
Previous studies have reported that AS changes are influenced by developmental and environmental factors [37, 38]. For instance, during the development of seed, endosperm and embryo, the AS isoforms change significantly in maize . In wheat, 251, 6618 and 7541 AS isoform changes were detected in response to drought stress, heat stress and combined stress in 9312 stress-responsive events, respectively . But little is known about the number of AS changes in response to cold acclimation. During cold acclimation, we found that the number of AS transcripts increased about 3000 in response to low temperature, but decreased by about 6000 during DA (Fig. 2); this implies that low temperature can induce new AS event and AS genes may be adapted to low-temperature environment, which would explain the reduced occurrence of AS events during DA. In addition, we made a comparison with two acclimatedplants of Longjing 43 (SRA database number: SRA061043) and Yinghong 9 (SRA database number: SRP108833), and found that the number of AS was also increased during CA2 in Yinghong 9 and CA3 in Longjing43, but with different extents, probably because the varieties and treatment conditions were different (Additional file 6: Table S6). Low temperature can affect not only the variation of AS number but also the expression level of AS transcripts to generate a large number of DAGs. In our study, the number of DAGs increased gradually with increasing cold acclimation (Fig. 4). This indicates that plants need to produce more functional AS to resist the gradually decreasing temperatures. Similarly in A.thaliana, DAGs also change with time at low temperatures, and these DAGs also gradually increased over the course of 4 °C treatment . The variation in the number of DAGs was also reported in drought-treated maize, and the number of DAGs increased significantly after 24 h of drought treatment compared with a 6 h of drought treatment . In addition, the number of DAGs was greater at de-acclimation compared with non-acclimation, which also suggests that the improvement of plant tolerance to low temperature after acclimation may be regulated by a large number of DAGs. To mitigate cold damage, the plant can eliminate reactive oxygen and maintain a redox balance by antioxidant metabolism [9, 42]. Our results showed that the oxidoreductase activity was enriched both at NA vs CS and NA vs DA (Fig. 4). Furthermore, many of the DAGs in NA vs DA are related to sugar metabolism; In A.thaliana, sugar metabolism also significantly increased during cold acclimation . Therefore, the ability of tea plant to tolerate low temperatures may be enhanced by changing the number and expression of AS, which thereby regulates oxidoreductase activity and sugar metabolism. Most unique DAGs are obtained from NA vs DA group, while the majority of differences on the metabolic level are in CA2 vs DA group. It is consistent with the reality that change of metabolic levels was lagged in comparison with the change of gene expression during transcription. Moreover, metabolite levels are often affected by various regulators at the same time in a complex regulate network, and the similar inconsistency between metabolites and transcripts has been reported .
The occurrence of AS is influenced by low temperature. For example, cold-induction resulted in alternative splicing of MtJMJC5 (JmjC domain-containing proteins) in Medicago truncatula, and the AS of MtJMJC5 could be reversed after the temperature returns to normal . In A.thaliana, plant biological clocks are also regulated by temperature associated AS events . Two intron-retained COR genes were previously found to be induced at low temperatures in durum wheat . In this study, we also found that the COR genes underwent AS events. In addition, a large number of AS events were found in the sugar metabolism and antioxidant pathways such as CsSUS, CsRS and CsPOD. All of these pathways have been associated with plant cold resistance. The ICE-CBF-COR regulatory pathway is one of the major ways how plants acquire resistance to cold [44, 45]. Plants can also protect their membrane structures and reduce cold damage through the synergistic action of superoxide dismutase (SOD), peroxidase (POD) protective enzymes and sugar metabolism [46, 47]. Our results suggest that AS may play crucial regulatory roles in these pathways in order to enhance the cold resistance of plants. Moreover, the tobacco transcription factor NtbHLH123 positively regulates the CBF Pathway and reactive oxygen species homeostasis to impart resistance to cold stress . However, the effect of AS on bHLH transcription factors during cold acclimation has not been previously explored. Most significantly, the CsbHLH isoforms were the predominant transcripts under cold treatments, and had higher expression levels than the full-length transcripts; this implies that the two AS isoforms may play a more important regulatory role and function for cold stress. As a precedent for this hypothesis, the Athaliana AS isoform FLM-β (β isoform of FLOWERING LOCUS M, a MADS-box transcription factor gene) is predominately formed at low temperatures and prevents precocious flowering by interacting with SVP (the Short Vegetative Phase gene, an MADS-box transcription factor) . The bHLH104 transcription factor has been reported to be involved in the sugar and ABA signaling pathways in A.thaliana . We found that the bHLH gene of TEA018946.1 has a high sequence similarity with that of A.thaliana, and thus, it was speculated that this bHLH AS isoform may be involved in some signaling pathways and in the downstream gene expression regulation during cold acclimation. The same variation trends of expression levels of this bHLH were also found in the acclimation tea plant of longjing43 and yinghong 9 under cold condition (Additional file 7: Table S7). In addition, The comparison of the AS events in the sugar metabolism and antioxidant pathways between Shuchazao and Longjing43 or Yinghong 9 was made respectively, and it was found that many AS events were the same (Additional file 8: Table S8).
It is well established that metabolites contents changed dramatically during cold acclimation [50, 51]. For example, most flavonols and anthocyanins accumulated during cold acclimation in A.thaliana , and the accumulation of trienoic fatty acid during acclimation is one of the necessary conditions for normal development of leaves at low temperature in tobacco . We found that sugar, flavonoids and fatty acids were three major classes of metabolites, and almost all of them were increased at low temperature, while quinic acid, tridecadienoic acid and octadecadienoic acid decreased. However, these reports only focused on the changes of metabolites and gene expression level during cold acclimation. In comparison, we found that these metabolites may be regulated by AS, and some metabolites remained at higher levels during de-acclimation; therefore, the presence of these metabolites may increase the resistance of plants to low temperature after acclimation. In addition, sugar has been widely studied as a class of metabolites playing a role in cold acclimation [5, 53]. In tea plant, the effects of cold acclimation on sugar metabolism and sugar-related gene expression has been reported , and in their study, the expression profiles of a large number of sugar-related genes were verified. However, a large number of AS events in sugar-related genes have not been reported. In our correlation analysis of AS and sugar content, levels of the full-length transcript of CsRS-1 (TEA013682.1) had a significant positive correlation with raffinose levels, while levels of its AS isoform CsRS-2 (TEA013682.2) were positively correlated with glucose, sucrose and fructose (Fig. 8). It is possible that the truncated proteins formed by the retention of introns did not have the catalytic activity of raffinose synthase, and the absence of raffinose promoted the synthesis of other soluble sugars. Changes in the functions of truncated proteins have also been reported in A.thaliana. For example, the JAZ gene can generate truncated protein due to the retention of introns, which weakens the interaction with JA receptor CORONATINE INSENSITIVE1 in the presence of the active JA-Ile conjugate; this imparts resistance to proteasomal degradation [54, 55]. We found that AS transcripts and full-length transcript of SUS have different correlations with metabolites and their domain were changed. This suggests that changes in the domain cause changes in gene function. For example, AS transcripts HAB1.1 and HAB1.2 play opposing roles in ABA-mediated seed germination and ABA-mediated post-germination developmental arrest in A. thaliana .
In conclusion, we analyzed the AS events during cold acclimation of tea plant based on the transcriptome. The results of this study revealed the dynamic changes in AS and its possible regulatory function during cold acclimation. Based on these results, we propose that changes in abundances and expression levels of AS transcripts may influence sugar metabolism and the antioxidant enzyme system to impart resistance to cold stress. Meanwhile, some transcription factors may regulate downstream genes or signal pathways through AS (Fig. 9). This study provides an essential insight into AS events involved in tea plant cold acclimation, and highlights the importance of AS regulatory function during low temperature stress.
Plant materials and growth conditions
One-year-old clone plants propagated with cuttings of (Camellia sinensis var. sinensis ‘Shuchazao’) grown in soil were obtained from the Dechang Tea Plantation in in Shucheng County, Anhui Province (116° 56′ 24″ E, 31° 27′N), China. The tea plant ‘Shuchazao’ (accession no. GS2002008; Agricultural Plant Variety Name Retrieval System, Ministry of Agriculture and Rural Affairs, China) is a national tea cultivar of China, which belongs to the subspecies Camellia sinensis var. sinensis of Camellia sinensis (L) O. Kuntze. Then the tea plants (30-40 cm) were transferred to a growth chamber at temperatures of 25 °C/20 °C for daytime/night-time. The culture conditions and treatments were described previously . The tea plants were subjected to cold acclimation and de-acclimation with the treatments as follows: Ten well-grown plants were collected for treatments of non-acclimation (NA), CS (sampling after 10 °C for 6 h), CA1 (sampling after 10/4 °C, day/night for 7 days), and CA2 (sampling after 4/0 °C, day/night for another 7 days). Finally, the treated plants of Shuchazao were placed at normal temperature (25/20 °C, day/night temperature) for 7 days of de-acclimation (DA). The leaf samples (CA1, CA2, DA) were collected at the same time (9:00 am). Three biological leaf samples were obtained at each time point. All the samples were immediately frozen in liquid nitrogen and then maintained at − 80 °C. However, the leaf samples were used in this experiment only had two biological replicates due to the consumption during previous experiments.
RNA extraction and sequencing
Total RNA was extracted from tea leaves using the Total RNA Purification Kit (NorgenBiotek Corporation, Canada) according to the manufacturer’s protocol. The quality and quantity of each RNA extract was determined using agarose gel electrophoresis and the Nanodrop 2500 (Thermo Fisher Scientific, US). RNA-Seq libraries were constructed as described previously  and Paired-end 150-bp sequencing was performed for the qualified library on a HiSeq X Ten machine. During sequencing, mRNA was enriched from total RNA using magnetic beads with Oligo (dT). Reads containing low-quality (low quality base, linker contamination and with unknown base N) were filtered by SOAPnuke (version 1.5.6) software , and the remaining reads (clean sequencing data) were used for subsequent analyses. The specific parameters are as follows: (1) Removing the reads containing the adapters (−f, −r); (2) Removing reads with unknown base N greater than 5% (−n); (3) Removing low-quality reads (defined as those containing > 20% bases with a quality value < 10) (−l, −q). The RNA-Seq data are publicly available at the SRA database of NCBI (https://www.ncbi.nlm.nih.gov/) under project accession number PRJNA387105.
Alignment and analysis of the sequencing data
The clean sequencing data were mapped to the tea plant reference genome  using HISAT2 (v.2.0.5) with default parameter. Transcripts were assembled by StringTie software (v.1.3.3b) with default settings. The transcripts per million reads (TPM) were calculated by StringTie and applied to quantify the transcript expression level. Ultimately, the assembled transcripts were aligned to the annotation by using Cuffcompare , and only transcripts with class_code ‘=’ or ‘j’ were considered as isoforms of known genes.
For the genic splice junctions, we took those coming from protein-coding genes and examined the predicted splice junctions within the coding sequence, 5′UTR or 3′UTR. If one extremity of the junction was in the coding sequence and the other in the UTR, they were classified as 5′UTR-CDS or CDS-3′UTR. The branch site and polypyrimidine tract for each intron were identified using a motif searching method implemented in a self-written PERL script, as described in previous study . Sequence logos were created using WebLogo 3 (http://weblogo.threeplusone.com/) [62, 63].
Differential expression analysis was performed using the DEGseq R package (1.34.1) . The analysis results using q-values < 0.01 and fold-change > 2 were considered to be differentially expressed. We defined these differentially expressed genes that undergo AS as differentially expressed AS genes (DAGs). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were conducted using OmicShare tool (http://omicshare.com/). GO terms with corrected P < 0.01 were considered significantly enriched by differential expression of AS genes and KEGG analysis with the Fisher’s exact test and false discovery rate (FDR) correction, visualizing the enrichment result with the R (v3.3.2) software. Correlation analysis of sugar levels and AS transcripts was also performed by R with ‘corrplot’ package (0.84).
Characterization of alternative splicing
To identify the type of AS events, the AStalavista tool  (http://astalavista.sammeth.net/) was employed using the gtf files assembled from the sequencing data by StringTie. Major types of AS events including intron retention (IR, AS code: 1^2-,0), exon skipping (ES, AS code: 1–2^, 0), alternative 3′ splice site (A3SS, AS code: 1-,2-), alternative 5′ splice site (A5SS, AS code: 1^,2^), IR1 + IR2 (AS code: 1^2–3^4-,0), A5SS or A3SS (AS code: 1^3-,2^4-), IR1 or IR2 (AS code: 1^2-,3^4-), A5SS + A3SS (AS code: 1^4-,2^3-), ES1 + ES2 (AS code: 1–2^3–4^,0), mutually exclusive exons (MXE, AS code: 1–2^,3–4^), A5SS + ES + A3SS (AS code: 1^6-,2^3–4^5-) were extracted from the output files and counted, respectively.
Validation of alternative splicing
Total RNA was isolated from tea plants at different time points of cold acclimation as described above. First-strand cDNA was synthesized using a PrimeScript RT Reagent Kit (cat 6110A, Takara, Japan) according to the manufacturer’s instructions. The alternative splicing isoforms associated with cold stress were viewed using IGV software . Gene structures were analyzed with online website Gene Structure Display Server (GSDS 2.0, http://gsds.cbi.pku.edu.cn/index.php). The AS transcripts were validated by RT-PCR according to a previous study . Details of all the relevant primers are listed in Additional file 5: Table S5. Conserved domain search services (CD Search, https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART (http://smart.embl-heidelberg.de/) was used to predict the domain of AS transcripts with the protein sequences.
Quantitative real-time PCR (qRT-PCR)
For qRT-PCR, first-strand cDNA was synthesized from total RNA with the PrimeScript RT Reagent Kit (cat RR036A, Takara, Japan) using the manufacturer’s protocols. PCR amplification and thermal cycling conditions were performed according to our previous study . The CsGAPDH gene was selected as the internal control. The relative gene expression values were analyzed using the 2-△Ct method . All reactions were run in technical triplicates for each sample. The relevant primers are listed in Additional file 5: Table S5.
Availability of data and materials
The data sets supporting the results of this article are available at the SRA database of NCBI (https://www.ncbi.nlm.nih.gov/) under project accession number PRJNA387105.
3′ untranslated region
5′ untranslated region
Alternative 3′ splice site
Alternative 5′ splice site
Cold acclimation of 7 days at 10/4 °C, day/night temperature
Cold acclimation of 7 days at 4/0 °C, day/night temperature
Cold stress of 6 h at 10 °C, day/night temperature
De-acclimation of 7 days at 25/20 °C, day/night temperature
Differentially AS genes
Kyoto Encyclopedia of Genes and Genomes
Premature stop codon
Quantitative real time polymerase chain reaction
Reverse transcription polymerase chain reaction
Thomashow MF. Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Biol. 1999;50(1):571–99.
Fowler S, Thomashow MF. Arabidopsis transcriptome profiling indicates that multiple regulatory pathways are activated during cold acclimation in addition to the CBF cold response pathway. Plant Cell. 2002;14(8):1675–90.
Hannah MA, Heyer AG, Hincha DK. A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Genet. 2005;1(2):e26.
Wei H, Dhanaraj AL, Rowland LJ, Fu Y, Krebs SL, Arora R. Comparative analysis of expressed sequence tags from cold-acclimated and non-acclimated leaves of Rhododendron catawbiense Michx. Planta. 2005;221(3):406–16.
Kaplan F, Kopka J, Sung DY, Zhao W, Popp M, Porat R, Guy CL. Transcript and metabolite profiling during cold acclimation of Arabidopsis reveals an intricate relationship of cold-regulated gene expression with modifications in metabolite content. Plant J. 2007;50(6):967–81.
Bolouri-Moghaddam MR, Le Roy K, Xiang L, Rolland F, Van den Ende W. Sugar signalling and antioxidant network connections in plant cells. FEBS J. 2010;277(9):2022–37.
Kodama H, Horiguchi G, Nishiuchi T, Nishimura M, Iba K. Fatty acid desaturation during chilling acclimation is one of the factors involved in conferring low-temperature tolerance to young tobacco leaves. Plant Physiol. 1995;107(4):1177–85.
Shi Y, Ding Y, Yang S. Molecular regulation of CBF signaling in cold acclimation. Trends Plant Sci. 2018;23(7):623–37.
Schöner S, Krause GH. Protective systems against active oxygen species in spinach: response to cold acclimation in excess light. Planta. 1990;180(3):383–9.
Zhang R, Calixto CP, Marquez Y, Venhuizen P, Tzioutziou NA, Guo W, Spensley M, Entizne JC, Lewandowska D, Ten Have S. A high quality Arabidopsis transcriptome for accurate transcript-level analysis of alternative splicing. Nucleic Acids Res. 2017;45(9):5061–73.
Laloum T, Martín G, Duque P. Alternative splicing control of abiotic stress responses. Trends Plant Sci. 2018;23(2):140–50.
Posé D, Verhage L, Ott F, Yant L, Mathieu J, Angenent GC, Immink RG, Schmid M. Temperature-dependent regulation of flowering by antagonistic FLM variants. Nature. 2013;503(7476):414.
Liu J, Chen X, Liang X, Zhou X, Yang F, Liu J, He SY, Guo Z. Alternative splicing of rice WRKY62 and WRKY76 transcription factor genes in pathogen defense. Plant Physiol. 2016;171(2):1427–42.
Ali G, Reddy A. Regulation of alternative splicing of pre-mRNAs by stresses. In: Nuclear pre-mRNA Processing in Plants: Springer Berlin Heidelberg; 2008. p. 257–75.
Gassmann W. Alternative splicing in plant defense. In: Nuclear pre-mRNA processing in plants: Springer Berlin Heidelberg; 2008. p. 219–33.
Dong C, He F, Berkowitz O, Liu J, Cao P, Tang M, Shi H, Wang W, Li Q, Shen Z. Alternative splicing plays a critical role in maintaining mineral nutrient homeostasis in Rice (Oryza sativa). Plant Cell. 2018;30(10):2267–85.
James AB, Syed NH, Bordage S, Marshall J, Nimmo GA, Jenkins GI, Herzyk P, Brown JW, Nimmo HG. Alternative splicing mediates responses of the Arabidopsis circadian clock to temperature changes. Plant Cell. 2012;24(3):961–81.
Calixto CP, Guo W, James AB, Tzioutziou NA, Entizne JC, Panter PE, Knight H, Nimmo HG, Zhang R, Brown JW. Rapid and dynamic alternative splicing impacts the Arabidopsis cold response transcriptome. Plant Cell. 2018;30(7):1424–44.
Gallegos JE. Alternative splicing plays a major role in plant response to cold temperatures. Plant Cell. 2018;30(7):1378–9 tpc. 00430.02018.
Mastrangelo AM, Belloni S, Barilli S, Ruperti B, Di Fonzo N, Stanca AM, Cattivelli L. Low temperature promotes intron retention in two e-cor genes of durum wheat. Planta. 2005;221(5):705–15.
Seo PJ, Park M-J, Park C-M. Alternative splicing of transcription factors in plant responses to low temperature stress: mechanisms and functions. Planta. 2013;237(6):1415–24.
Gramza-Michałowska A. Functional aspects of tea Camellia sinensis as traditional beverage. In: Functional Properties of Traditional Foods: Springer US; 2016. p. 353–69.
Liu S, An Y, Li F, Li S, Liu L, Zhou Q, Zhao S, Wei C. Genome-wide identification of simple sequence repeats and development of polymorphic SSR markers for genetic studies in tea plant (Camellia sinensis). Mol Breed. 2018;38:1–13.
Sueoka N, Suganuma M, Sueoka E, Okabe S, Matsuyama S, Imai K, Nakachi K, Fujiki H. A new function of green tea: prevention of lifestyle-related diseases. Ann N Y Acad Sci. 2001;928(1):274–80.
Zaveri NT. Green tea and its polyphenolic catechins: medicinal uses in cancer and noncancer applications. Life Sci. 2006;78(18):2073–80.
Vyas D, Kumar S. Tea (Camellia sinensis (L.) O. Kuntze) clone with lower period of winter dormancy exhibits lesser cellular damage in response to low temperature. Plant Physiol Biochem. 2005;43(4):383–8.
Yue C, Cao H-L, Wang L, Zhou Y-H, Huang Y-T, Hao X-Y, Wang Y-C, Wang B, Yang Y-J, Wang X-C. Effects of cold acclimation on sugar metabolism and sugar-related gene expression in tea plant during the winter season. Plant Mol Biol. 2015;88(6):591–608.
Ban Q, Wang X, Pan C, Wang Y, Kong L, Jiang H, Xu Y, Wang W, Pan Y, Li Y. Comparative analysis of the response and gene regulation in cold resistant and susceptible tea plants. PLoS One. 2017;12(12):e0188514.
Zhu J, Wang X, Xu Q, Zhao S, Tai Y, Wei C. Global dissection of alternative splicing uncovers transcriptional diversity in tissues and associates with the flavonoid pathway in tea plant (Camellia sinensis). BMC Plant Biol. 2018;18(1):266.
Zhu J, Wang X, Guo L, Xu Q, Zhao S, Li F, Yan X, Liu S, Wei C. Characterization and alternative splicing profiles of the Lipoxygenase gene family in tea plant (Camellia sinensis). Plant Cell Physiol. 2018;59(9):1765–81.
Filichkin SA, Priest HD, Givan SA, Shen R, Bryant DW, Fox SE, Wong W-K, Mockler TC. Genome-wide mapping of alternative splicing in Arabidopsis thaliana. Genome Res. 2010;20(1):45–58.
Patel AA, McCarthy M, Steitz JA. The splicing of U12-type introns can be a rate-limiting step in gene expression. EMBO J. 2002;21(14):3804–15.
Leviatan N, Alkan N, Leshkowitz D, Fluhr R. Genome-wide survey of cold stress regulated alternative splicing in Arabidopsis thaliana with tiling microarray. PLoS One. 2013;8(6):e66511.
Liu GF, Liu JJ, He ZR, Wang FM, Yang H, Yan YF, Gao MJ, Gruber MY, Wan XC, Wei S. Implementation of CsLIS/NES in linalool biosynthesis involves transcript splicing regulation in Camellia sinensis. Plant Cell Environ. 2018;41(1):176–86.
Wang X-C, Zhao Q-Y, Ma C-L, Zhang Z-H, Cao H-L, Kong Y-M, Yue C, Hao X-Y, Chen L, Ma J-Q. Global transcriptome profiles of Camellia sinensis during cold acclimation. BMC Genomics. 2013;14(1):415.
Zheng C, Wang Y, Ding Z, Zhao L. Global transcriptional analysis reveals the complex relationship between tea quality, leaf senescence and the responses to cold-drought combined stress in Camellia sinensis. Front Plant Sci. 2016;7:1858.
Szakonyi D, Duque P. Alternative splicing as a regulator of early plant development. Front Plant Sci. 2018;9:1174.
Calixto CP, Tzioutziou NA, James AB, Hornyik C, Guo W, Zhang R, Nimmo HG, Brown JW. Cold-dependent expression and alternative splicing of Arabidopsis long non-coding RNAs. Front Plant Sci. 2019;10:235.
Thatcher SR, Danilevskaya ON, Meng X, Beatty M, Zastrow-Hayes G, Harris C, Van Allen B, Habben J, Li B. Genome-wide analysis of alternative splicing during development and drought stress in maize. Plant Physiol. 2016;170(1):586–99.
Liu Z, Qin J, Tian X, Xu S, Wang Y, Li H, Wang X, Peng H, Yao Y, Hu Z. Global profiling of alternative splicing landscape responsive to drought, heat and their combination in wheat (Triticum aestivum L.). Plant Biotechnol J. 2018;16(3):714–26.
Mei W, Liu S, Schnable JC, Yeh C-T, Springer NM, Schnable PS, Barbazuk WB. A comprehensive analysis of alternative splicing in paleopolyploid maize. Front Plant Sci. 2017;8:694.
Song L, Guo S. Study on the relations of esterase, peroxidase and cold hardiness of winter rape (Brassica napus L.). J Anhui Agric Sci (China). 1995;23:4.
Shen Y, Wu X, Liu D, Song S, Liu D, Wang H. Cold-dependent alternative splicing of a Jumonji C domain-containing gene MtJMJC5 in Medicago truncatula. Biochem Biophys Res Commun. 2016;474(2):271–6.
Wang Y, Jiang C-J, Li Y-Y, Wei C-L, Deng W-W. CsICE1 and CsCBF1. Two transcription factors involved in cold responses in Camellia sinensis. Plant Cell Rep. 2012;31(1):27–34.
Yn J, Zhai S, Wang W, Ding X, Guo Z, Bai L, Wang S. Identification of genes from the ICE–CBF–COR pathway under cold stress in Aegilops–Triticum composite group and the evolution analysis with those from Triticeae. Physiol Mol Biol Plants. 2018;24(2):211–29.
Cao Y-Y, Yang M-T, Li X, Zhou Z-Q, Wang X-J, Bai J-G. Exogenous sucrose increases chilling tolerance in cucumber seedlings by modulating antioxidant enzyme activity and regulating proline and soluble sugar contents. Sci Hortic. 2014;179:67–77.
Zhou M, Li W, Zheng Y, Lin P, Yao X, Lin J. CbRCI35, a cold responsive peroxidase from Capsella bursa-pastoris regulates reactive oxygen species homeostasis and enhances cold tolerance in tobacco. Front Plant Sci. 2016;7:1599.
Zhao Q, Xiang X, Liu D, Yang A, Wang Y. Tobacco transcription factor NtbHLH123 confers tolerance to cold stress by regulating the NtCBF pathway and reactive oxygen species homeostasis. Front Plant Sci. 2018;9:381.
Min J-H, Park C-R, Jang Y-H, Ju H-W, Lee K-H, Lee S, Kim CS. A basic helix-loop-helix 104 (bHLH104) protein functions as a transcriptional repressor for glucose and abscisic acid signaling in Arabidopsis. Plant Physiol Biochem. 2019;136:34–42.
Vandenbussche B, Leurdian S, Verdoodt V, Gysemberg M, De Proft M. Changes in sugar content and fatty acid composition of in vitro sugar beet shoots after cold acclimation: influence on survival after cryopreservation. Plant Growth Regul. 1999;28(3):157–63.
Schulz E, Tohge T, Zuther E, Fernie AR, Hincha DK. Flavonoids are determinants of freezing tolerance and cold acclimation in Arabidopsis thaliana. Sci Rep. 2016;6:34027.
Schulz E, Tohge T, Zuther E, Fernie AR, Hincha DK. Natural variation in flavonol and anthocyanin metabolism during cold acclimation in Arabidopsis thaliana accessions. Plant Cell Environ. 2015;38(8):1658–72.
Burbulis N, Jonytiene V, Kupriene R, Blinstrubiene A. Changes in proline and soluble sugars content during cold acclimation of winter rapeseed shoots in vitro. J Food Agric Environ. 2011;9:371–4.
Chung HS, Cooke TF, DePew CL, Patel LC, Ogawa N, Kobayashi Y, Howe GA. Alternative splicing expands the repertoire of dominant JAZ repressors of jasmonate signaling. Plant J. 2010;63(4):613–22.
Moreno JE, Shyu C, Campos ML, Patel LC, Chung HS, Yao J, He SY, Howe GA. Negative feedback control of jasmonate signaling by an alternative splice variant of JAZ10. Plant Physiol. 2013;162(2):1006–17.
Wang Z, Ji H, Yuan B, Wang S, Su C, Yao B, Zhao H, Li X. ABA signalling is fine-tuned by antagonistic HAB1 variants. Nat Commun. 2015;6:8138.
Li Y, Wang X, Ban Q, Zhu X, Jiang C, Wei C, Bennetzen JL. Comparative transcriptomic analysis reveals gene expression associated with cold adaptation in the tea plant Camellia sinensis. BMC Genomics. 2019;20(1):624.
Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, Li Y, Ye J, Yu C, Li Z. SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience. 2017;7(1):gix120.
Wei C, Yang H, Wang S, Zhao J, Liu C, Gao L, Xia E, Lu Y, Tai Y, She G. Draft genome sequence of Camellia sinensis var. sinensis provides insights into the evolution of the tea genome and tea quality. Proc Natl Acad Sci. 2018;115(18):E4151–8.
Li Y, Dai C, Hu C, Liu Z, Kang C. Global identification of alternative splicing via comparative analysis of SMRT-and Illumina-based RNA-seq in strawberry. Plant J. 2017;90(1):164–76.
Xie B-B, Li D, Shi W-L, Qin Q-L, Wang X-W, Rong J-C, Sun C-Y, Huang F, Zhang X-Y, Dong X-W. Deep RNA sequencing reveals a high frequency of alternative splicing events in the fungus Trichoderma longibrachiatum. BMC Genomics. 2015;16(1):54.
Schneider TD, Stephens RM. Sequence logos: a new way to display consensus sequences. Nucleic Acids Res. 1990;18(20):6097–100.
Crooks GE, Hon G, Chandonia J-M, Brenner SE. WebLogo: a sequence logo generator. Genome Res. 2004;14(6):1188–90.
Wang L, Feng Z, Wang X, Wang X, Zhang X. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics. 2009;26(1):136–8.
Foissac S, Sammeth M. ASTALAVISTA. dynamic and flexible analysis of alternative splicing events in custom gene datasets. Nucleic Acids Res. 2007;35(suppl_2):W297–9.
Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14(2):178–92.
Wang Y-N, Tang L, Hou Y, Wang P, Yang H, Wei C-L. Differential transcriptome analysis of leaves of tea plant (Camellia sinensis) provides comprehensive insights into the defense responses to Ectropis oblique attack using RNA-Seq. Funct Integr Genomics. 2016;16(4):383–98.
Zhao S, Wang X, Yan X, Guo L, Mi X, Xu Q, Zhu J, Wu A, Liu L, Wei C. Revealing of MicroRNA involved regulatory gene networks on Terpenoid biosynthesis in Camellia sinensis in different growing time points. J Agric Food Chem. 2018;66(47):12604–16.
We thank Dechang Tea Plantation in Shucheng County, Anhui Province, China. And Tea Plant Cultivar and Germplasm Resource Garden in Guohe Town, Anhui Agricultural University for providing samples of tea plants.
This work was supported by the National Key Research and Development Program of China (2019YFD1001601, 2018YFD1000601), the National Natural Science Foundation of China (31171608), the Special Innovative Province Construction in Anhui Province (15czs08032), the Special Project for Central Guiding Science and Technology Innovation of Region in Anhui Province (2016080503B024).
Ethics approval and consent to participate
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Table S1.
The number and percentages of splicing donor-acceptor di-nucleotide usages among all transcripts.
Additional file 2: Table S2.
Characteristics of cold stress-responsive AS genes in tea plant.
Additional file 3: Table S3.
The genome and sequencing sequence of the AS transcript.
Additional file 4: Table S4.
The data for each sequencing sample.
Additional file 5: Table S5.
Sequences of primer pairs used for verifying AS transcripts and qRT-PCR analysis.
Additional file 6: Table S6.
Alternative splicing analysis of Longjing 43 and Yinghong 9.
Additional file 7: Table S7.
Expression of bHLH in Shuchazao, Longjing 43 and Yinghong 9.
Additional file 8: Table S8.
Comparison of AS events in Shuchazao, Longjing43 and Yinghong 9.
Additional file 9: Figure S1.
The picture of tea plantation and tea plants.
Additional file 10: Figure S2.
Percentages of four main AS types during cold acclimation.
Additional file 11: Figure S3.
Expression analysis of AS transcripts during cold acclimation.
Additional file 12: Figure S4.
Correlation analysis of AS transcriptions and sugar content.
Additional file 13: Figure S5.
Domain analysis of SUS and RS in AS transcripts.
Additional file 14: Figure S6.
Expression analysis of qRT-PCR and sequencing results.
Additional file 15: Figure S7.
Verification of alternatively spliced isoforms.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Li, Y., Mi, X., Zhao, S. et al. Comprehensive profiling of alternative splicing landscape during cold acclimation in tea plant. BMC Genomics 21, 65 (2020). https://doi.org/10.1186/s12864-020-6491-6
- AS isoforms
- Camellia sinensis
- Cold adaptation
- Cold stress