A comparison of the low temperature transcriptomes of two tomato genotypes that differ in freezing tolerance: Solanum lycopersicum and Solanum habrochaites

Solanum lycopersicum and Solanum habrochaites are closely related plant species; however, their cold tolerance capacities are different. The wild species S. habrochaites is more cold tolerant than the cultivated species S. lycopersicum. The transcriptomes of S. lycopersicum and S. habrochaites leaf tissues under cold stress were studied using Illumina high-throughput RNA sequencing. The results showed that more than 200 million reads could be mapped to identify genes, microRNAs (miRNAs), and alternative splicing (AS) events to confirm the transcript abundance under cold stress. The results indicated that 21 % and 23 % of genes were differentially expressed in the cultivated and wild tomato species, respectively, and a series of changes in S. lycopersicum and S. habrochaites transcriptomes occur when plants are moved from warm to cold conditions. Moreover, the gene expression patterns for S. lycopersicum and S. habrochaites were dissimilar; however, there were some overlapping genes that were regulated by low temperature in both tomato species. An AS analysis identified 75,885 novel splice junctions among 172,910 total splice junctions, which suggested that the relative abundance of alternative intron isoforms in S. lycopersicum and S. habrochaites shifted significantly under cold stress. In addition, we identified 89 miRNA sequences that may regulate relevant target genes. Our data indicated that some miRNAs (e.g., miR159, miR319, and miR6022) play roles in the response to cold stress. Differences in gene expression, AS events, and miRNAs under cold stress may contribute to the observed differences in cold tolerance of these two tomato species.


Background
Some plants increase their cold tolerance to deal with low temperatures; this phenomenon is termed cold acclimation. During this process, various biochemical and physiological changes occur in plants, which make plants more cold tolerant. Plants have different abilities to deal with low temperatures. Plants that have adapted to cold environments increase their cold tolerance in response to low but nonfreezing temperatures. By contrast, plants that have adapted to subtropical and tropical climates, such as maize, rice, and tomato, generally have little cold tolerance and are unable to acclimate to cold temperatures [1].
The roles of cold-regulated genes in plant cold acclimation show that differential expression of genes is related to different cold adaption abilities of plants. In A. thaliana and Chorispora bungeana, many alterations in gene expression begin within minutes of transferring plants to a low temperature [15][16][17][18][19]. Moreover, some studies have demonstrated that differential expression of a cold-responsive gene is caused by differences in cold tolerance in plants [20,21]. For example, there are considerable differences in cold-responsive genes in Solanum tuberosum and Solanum commersonii, which are closely related but have different cold tolerances [22].
A large number of cold-related genes have been identified using transcript analysis techniques, and the products of cold-related genes, including regulatory proteins and functional proteins, are thought to play key roles in gene expression and signal transduction [3,15,17,[23][24][25][26][27][28]. The expression and alternative splicing (AS) of some serine/arginine-rich (SR) genes, which encode splicing factor proteins that are vital for splicing and constitutive expression, vary under cold stress [29][30][31][32][33]. Cold stress affects the expression of splicing factors; therefore, the splicing of precursor-mRNAs (pre-mRNAs) of other genes are altered under cold stress. AS of pre-mRNA is an important mechanism for increasing transcriptome and proteome variety in eukaryotes. AS has been confirmed widely at the functional level in A. thaliana, rice, and maize [34][35][36]. AS may be regulated spatially and developmentally under environmental stress [33,[37][38][39]; thus, AS could play an important role under cold stress or other abiotic stress.
MicroRNAs (miRNAs) have been discovered in plants [40][41][42], changing our perception of the mechanisms of gene expression, transcription, and translation. In plants, [21][22][23][24] miRNAs are negative regulators of gene expression [43]. The pool of miRNAs in plants is highly diverse [44][45][46]. Many studies have indicated that the key role of miRNAs is regulating organ development and biological processes [47][48][49]. MiRNAs are also associated with abiotic stress responses [50][51][52]. In accordance with their regulatory roles, many miRNAs target genes that have roles in developmental patterning and show unique development-specific, tissue-specific, and stress-induced expression patterns [47,53,54]. However, to date, only 44 annotated tomato miRNAs have been deposited in the miRBase database v19.0, and only a few miRNA targets have been confirmed experimentally. At present, it is unknown whether important regulators such as miR-NAs play a vital role in tomato's response to cold stress.
The cultivated tomato species (Solanum lycopersicum) suffers from cold stress at all stages of growth and development; by contrast, the wild tomato species (Solanum habrochaites) grows well at low temperatures [55][56][57]. To understand the molecular basis underlying why S. habrochaites can acclimate to cold and survive freezing temperatures, whereas S. lycopersicum cannot, we report the results of an RNA sequencing (RNA-seq) transcriptome and miRNA analysis of RNA populations obtained from cold-treated leaves of the two plants. The results showed that many changes in the S. lycopersicum and S. habrochaites transcriptomes occur in plants transferred from warm to cold conditions. We predicted that at least 21 % and 23 % of genes were cold responsive in S. lycopersicum and S. habrochaites, respectively. A gene ontology (GO) term enrichment analysis of the data indicated that many GO terms ("abiotic stimulus response", "ethylene stimulus response") were significantly enriched in the cold-responsive genes between the two species. Our data also provided an evaluation of AS between S. lycopersicum and S. habrochaites. RNA-seq identified many annotated introns, known AS, and 75,885 novel splice junctions. We identified 89 miRNA sequences and 423 targets of 83 miRNAs. Our data showed that some miR-NAs (e.g., miR159, miR319, and miR6022) play important roles under cold stress in the two species. These results provided a new insight into the roles of miRNAs under cold stress in these two closely related species under cold stress. Thus, the differences in gene expression, AS events, and miRNAs under cold stress may contribute to the differences in the cold tolerance between S. lycopersicum and S. habrochaites.

Results
Phenotypic and physiological responses to cold stress Solanum lycopersicum and S. habrochaites leaf tissue were chosen to study cold responses. The degree of cold stress was identified by malondialdehyde (MDA) content, proline content, peroxidase (POD) activity, and catalase (CAT) activity exchange in the leaves. S. habrochaites exhibited less severe wilting than S. lycopersicum after 10 days of treatment at 4°C (Additional file 20: Figure S1A-D). Cold stress caused significant increased MDA content, proline content, POD activity, and CAT activity in these plants (Additional file 20: Figure S1E-H).

Solanum lycopersicum and S. habrochaites transcriptome analyses
The transcriptomes of S. lycopersicum (C) and S. habrochaites (Tsh) under cold stress were analyzed by RNAseq using the Illumina Genome Analyzer II. After growth at 25°C for 12 weeks, plants were moved to 4°C for 0, 1, and 12 h, and the total RNA from leaves was extracted and analyzed. More than 200 million reads were produced, with approximately 33.3 million reads from each sample. We aligned the reads to the entire reference genome sequence of tomato (version SL2.40; http://solgenomics.net/) using the TopHat tool. The tolerance was set to allow two mismatches at most for reads in each alignment. Using these criteria, 96.32-97.21 % of the reads mapped uniquely to a genomic location, and 2.79-3.68 % of the reads were filtered out as multiple-mapped or low-quality reads. Alignment of the reads to tomato cDNAs demonstrated that 70 % of the tomato genome-annotated cDNAs had a sequence represented by Illumina RNA-seq reads (Table 1). Compared with the annotated tomato genome, the RNA-seq data revealed that 92.5-95 % of the reads that matched to the genome mapped to annotated genic regions, but only 5-7.5 % of the reads mapped to intergenic regions (Additional file 20: Figure S2). The depth of coverage along the length of the transcripts reduced towards the 5′ termini for RNA-seq data derived from the full-length cDNA libraries (Additional file 20: Figure S3). De novo assembly was performed using the Trinity method with default parameters [58]. Overviews of the assembly results are shown in Additional file 18 and Additional file 19. The reads were assembled into 68,051 nonredundant unigenes (>200 bp) in S. habrochaites.
To verify the RNA-seq data, some genes whose expressions increased, some that decreased, and some that showed no change in abundance were chosen for realtime PCR (qPCR) under cold stress. The results of RNA-seq and qPCR were similar (Additional file 20: Figure S8), showing the same general expression trends by qPCR as were revealed by RNA-seq.
To identify S. lycopersicum and S. habrochaites miRNAs that affect gene regulation under cold stress, six miRNA libraries were constructed from the leaves of S. lycopersicum and S. habrochaites that were or were not treated with cold. The six miRNA libraries were sequenced using high-throughput RNA-seq and yielded approximately 5.4 million raw reads in each sample. We excluded the poorquality reads and those whose length was smaller than 14 nucleotides from further analysis. Finally, we obtained approximate 4.2 million non-redundant reads (14-24 nucleotides) in each sample (Table 1).

Differential expression and GO enrichment
To study the impact of cold stress on gene expression in S. lycopersicum (C) and S. habrochaites (Tsh), the transcript abundance of each gene was identified by the reads per kilobase of transcript per million reads mapped method (Additional file 1, Additional file 20: Figure S4). To compare the transcriptomes in S. lycopersicum and S. habrochaites under cold stress, a heat map was generated to present the transcript abundance for all differentially expressed genes (DEGs) under cold stress at 0, 1, and 12 h (Fig. 1, Additional file 20: Figure  S5-7). The results showed that a series of changes in gene expression in S. lycopersicum and S. habrochaites occur when plants are moved from warm to cold conditions. Moreover, the gene expression patterns for S. lycopersicum and S. habrochaites were dissimilar. For example, cluster A genes were little affected at 1 h in S. lycopersicum and returned to low transcript abundance levels at 12 h of cold stress; cluster B genes were unaffected in S. lycopersicum at 1 h of cold stress, but were highly increased at 12 h; cluster C or D genes were little affected after cold stress in S. lycopersicum, but were affected in S. habrochaites ( Fig. 1). We used a threshold of a minimum 2-fold change in abundance between any two time points to define DEGs in the following analysis (Fig. 2, Table 2, Additional file 3). The results showed that~4 % (sample C1 vs. C0), 10 % (sample C12 vs. C0),~5 % (sample Tsh1 vs. Tsh0), and~8 % (sample Tsh12 vs. Tsh0) of the unigenes were cold induced; and~2 % (sample C1 vs. C0),~9 % (sample C12 vs. C0),~6 % (sample Tsh1 vs. Tsh0), and 9 % (sample Tsh12 vs. Tsh0) were cold repressed. In S. lycopersicum, transcripts for 1,256 and 3,350 unigenes increased at 1 and 12 h, respectively, and 804 unigenes increased at both time points tested; transcripts for 856 and 3,022 unigenes decreased at 1 and 12 h, respectively, and 339 unigenes decreased at both time points tested (Fig. 2, Table 2, Additional file 3, Additional file 4). In S. habrochaites, transcripts for 1,725 and 2,940 unigenes increased at 1 and 12 h, respectively, and 722 unigenes increased at both time points tested; transcripts for 1,967 and 3,126 unigenes decreased at 1 and 12 h, respectively, and 1,000 unigenes decreased at both time points tested. Moreover, in S. habrochaites, transcripts for 3,608, 2,813, and 3,549 unigenes increased at 0, 1, and 12 h, respectively, compared with S. lycopersicum at same time points; and transcripts for 3,897, 3,592, and 3,815 unigenes decreased at 0, 1, and 12 h, respectively. In sum, the gene expression profiles in S. lycopersicum and S. habrochaites changed under cold stress to different degrees; however, there were some overlapping genes that were regulated by low temperature in both tomato species.
We analyzed the genes that we determined to be responsive to cold at 1 h. The GO terms enriched in each species were comparable (Additional file 3) and were generally related to "response to abiotic stimulus". From the heat map, it was also obvious that S. lycopersicum was less affected by cold than S. habrochaites at 1 h. The expressions of genes that were enriched in GO categories corresponding to "cell wall metabolism" were increased under cold stress in S. lycopersicum, but the opposite result was observed in S. habrochaites. We observed a similar contrast in the GO category "response to organic substance". In the GO categories "response to chitin", "response to carbohydrate stimulus", and "DNAbinding WRKY", there was a significant enrichment in S. lycopersicum, but S. habrochaites showed no enrichment. For the categories "chloroplast", "transit peptide", "pentatricopeptide repeat", "phenylpropanoid metabolic process", "flavonoid metabolic process", and "amino acid derivative biosynthetic process", no significant enrichment was observed in S. lycopersicum, but enrichment was observed in S. habrochaites (Additional file 3).
We then compared responses to cold at 12 h. The analysis of GO terms for cold-regulated genes suggested that the categories "response to organic substance", "response to endogenous stimulus", "response to hormone stimulus", "response to abscisic acid stimulus", "pentatricopeptide repeat", "response to abiotic stimulus", "response to ethylene stimulus", "serine/threonine-protein kinase", "phenylpropanoid metabolic process", "amino acid derivative biosynthetic process", "lignin biosynthetic process", and "flavonoid metabolic process" were enriched in both S. lycopersicum and S. habrochaites (Additional file 3). In the case of the GO category "UDP-glucuronosyl/UDP-glucosyltransferase", there was significant enrichment for S. lycopersicum, but not for S. habrochaites.
Alternative splicing in S. lycopersicum and S. habrochaites To study the effect of cold stress on AS in S. lycopersicum and S. habrochaites, we compared splicing events   (Fig. 3a) and receptor-like protein (RLK) (Solyc01g007980.2) (Fig. 3b), are shown in Fig. 3. The TopHat-generated S. lycopersicum LOB domain protein 38 mRNA model predicted a distinct AS event yielding a splice isoform that retains full-length intron 1 (Fig. 3a, Additional file 7). An analysis of RLK, a putative resistance protein with an antifungal domain, provides another example of an IntronR event in plants. The depth of coverage of the intron 3 splice junction was confirmed by RNA-seq (Fig. 3b, Additional file 7). Dense micro-read coverage of intron 3 in the RLK transcript contrasted with the low coverage of other introns, indicating that intron 3 may be retained in some mature RLK transcripts (Fig. 3b).
We tried to identify differences in altered AS events between the two species at 0, 1, and 12 h of cold treatment at 4°C. Then, 270 (sample Tsh0 vs.  (Table 4). Table 4 shows that AS occurred more frequently in genes regulated in response to cold at 12 h than in genes at 1 h.
The AS event of diacylglycerol acyltransferase (Solyc02g068240.2) under cold stress is shown in Fig. 3c. The TopHat-generated S. habrochaites diacylglycerol acyltransferase mRNA model predicted a distinct AS event that yielded a splice isoform that retains intron 4 ( Fig. 3c). Accumulation of the no IntronR 4-containing transcripts decreased approximately three-fold under cold treatment. Other examples of cold stress-associated AS genes (SR45a, SR30) are provided in Additional file 20: Figure S9.
We compared the functions of the AS genes that were regulated in response to cold at 1 h and 12 h with the DEGs. Cold-regulated differentially expressed AS genes overlapped with DEGs in S. lycopersicum (C) and S. habrochaites (Tsh), and these genes were in the GO categories "dephosphorylation" and "phosphoprotein phosphatase activity" (Additional file 8), suggesting these activities were present in both plants. In the case of the GO categories "detection of light stimulus", "phenylpropanoid metabolic process", "response to cadmium ion", "phosphoinositide binding", and "heat shock protein binding", there was significant enrichment for S. lycopersicum, but S. habrochaites showed no enrichment. For the categories "carboxylic acid catabolic process", "proteolysis", "cell death", "reproductive developmental process", and "ethylene mediated signaling pathway", no enrichment was observed in S. lycopersicum, but significant enrichment was observed in S. habrochaites (Additional file 8).

Impact of cold stress on miRNAs in tomato
To identify miRNAs in tomato, we analyzed miRNAs by BLAST searches against the tomato genome sequence by The sequences corresponding to the known non-coding RNAs (tRNAs, rRNAs, small nucleolar and small nuclear RNAs) were filtered out using BLASTn to search the Rfam database (http://rfam.xfam.org/) (Additional file 20: Figure S11). The remaining sequences were assigned as either other endogenous small RNAs or miRNA candidates and used for a fold-back structure prediction. We compared the unique miRNAs with the miRBase  Fig. 4 The total number of differentially alternative splicing genes (DASG) in S. lycopersicum (C) and S. habrochaites (Tsh) at 0, 1, and 12 h of cold treatment database (version 19.0). In this analysis, the miRNAs were required to show a perfect or nearly perfect match (mismatch ≤ 1) to known miRNAs. After these analyses, 112 unique miRNA were obtained as novel miRNA candidates (Additional file 12). A large number of miRNA sequences were produced by Illumina sequencing, permitting us to confirm the relative abundance of miRNAs in tomato. To study the expression dynamics of miRNAs and their potential roles in gene expression regulation in S. lycopersicum and S. habrochaites, the transcript abundance of each miRNA was evaluated by transcripts per million (TPM). The TPM of the miRNAs varied from 0 to 27,670 (miR396a, sample C1), suggesting that the expression of miRNAs varied greatly in tomato (Additional file 13). MiR159 and miR396a were the most abundant miRNAs in the six sequencing datasets. According to the TPM, some miRNAs (miR159, miR396a, miR396b, miR482b, and miR6022) were highly expressed in tomato, with a TPM of greater than 100 each. MiR6027, miR171a, miR482, miR319, and miR1919a were moderately expressed and had a TPM between 10 and 100. MiR5303, miR169b, miR1916, miR171c, and miR395a represent miRNAs with low expression and a TPM of less than 10 (Additional file 13). Sequence analysis showed that the relative abundance of some members in one miRNA family changed considerably in tomato. For example, the TPM for miR396a was 9,433, whereas the TPM for miR396b was only 4,347 (Additional file 13).
To detect the effect of cold stress on miRNAs, the expression of miRNAs in S. lycopersicum and S. habrochaites seedlings with and without cold treatment was examined. Fourteen (sample C1 vs. C0), eight (sample C12 vs. C0), two (sample Tsh1 vs. Tsh0), and four (sample Tsh12 vs. Tsh0) of the miRNAs were cold induced; seven (sample C1 vs. C0), six (sample C12 vs. C0), five (sample Tsh1 vs. Tsh0), and eight (sample Tsh12 vs. Tsh0) of the miRNAs were cold repressed (Fig. 5, Additional file 14). In response to cold treatment, the most significant change was observed for miR169c, whose expression level increased approximately 35-fold in sample C1 compared with C0. The expressions of some miR-NAs in S. habrochaites were opposite to those in S. lycopersicum under cold stress. For example, miR1919a-miR1919c, and miR396b were upregulated under cold stress for 1 h in S. lycopersicum, whereas they were downregulated in S. habrochaites (Additional file 14). In contrast, miR172a and miR172b were downregulated by cold stress for 1 h in S. lycopersicum, while they upregulated in S. habrochaites by cold stress for 12 h.
We used psRNATarget (http://plantgrn.noble.org/ psRNATarget/) to predict targets for the miRNAs. For the miRNAs that were annotated as described above, we identified 423 mRNA targets (Fig. 6, Additional file 15). From Fig. 6 it was also evident that S. lycopersicum was more affected by 1 h of cold than S. habrochaites. To further characterize the role of the miRNAs in response to cold, we examined the target list for genes that could be related to the cold response and that were either induced or repressed by cold, based on our Illumina results (Additional file 17). For example, one of the predicted targets was the transcript of the homeodomain leucine zipper class I (HD-Zip I) protein (ATHB13, AT1G69780, Solyc02g087840.2) (Additional file 17). ATHB13 is induced in S. lycopersicum after cold treatment for 12 h based on our sequencing data (Additional file 3). The miRNA predicted to target ATHB13 is miR6022. Our sequencing data showed that miR6022 was downregulated in S. lycopersicum after cold stress for 12 h (Additional file 14). Based on our sequencing data, we did not find differential expression of miR6022 after cold treatment for 1 h in S. lycopersicum. The induction of ATHB13 under cold stress for 12 h correlates with miR6022 repression by cold, suggesting that ATHB13 levels are post-transcriptionally regulated by this miRNA in response to cold. Thus, miR6022/ ATHB13 represents an abiotic stimulus module that could be important for the cold response in S. lycopersicum leaves. Other examples of cold stress-associated miRNAs (miR159, miR319) are provided in Additional file 17.
For comprehensive annotation, all putative target genes in each sample were analyzed by GO terms using the DAVID program. An analysis of GO enrichment for the targets revealed that target functions were enriched in many different biological processes (Additional file 16). Among the mRNA targets that were upregulated in response to cold at 1 h and 12 h, comparable coldregulated mRNA target expression was observed between S. lycopersicum (C) and S. habrochaites (Tsh) in relation to the GO terms (Additional file 16), which included "ATP binding" and "nucleotide binding" in both species. In the case of the GO categories "leaf development", "shoot development", "CCAAT-binding factor", "CBF", "regulation of RNA metabolic process", "cell death", "gene silencing", "immune response", "flower development", "ATPase activity", and "leucine-rich repeat", there was significant enrichment in S. lycopersicum, but not in S. habrochaites. For the categories "response to cold" and "hormone stimulus", no enrichment was observed in S. lycopersicum, but significant enrichment was observed in S. habrochaites (Additional file 16). Among the targets that were determined as downregulated in response to cold at 1 h and 12 h, an analysis of GO enrichment showed "defense response", "ATP binding", and "nucleotide binding" in both S. lycopersicum and S. habrochaites. In the case of the GO categories "leucine-rich repeat", "reproductive structure development", "meristem development", "intracellular signaling cascade", and "response to hormone stimulus", there was no significant enrichment in S. habrochaites, but S. lycopersicum showed enrichment.

Discussion
Plants have different abilities to deal with low temperatures. The cultivated tomato (S. lycopersicum) suffers from cold stress, but the wild species (S. habrochaites) does not [55][56][57]. RNA-seq of cold-stressed S. lycopersicum leaves and a comparison with the transcriptome of S.
habrochaites are presented here. The results revealed the effects of cold stress on transcript abundance in S. lycopersicum and S. habrochaites; 21 % and 23 % of transcripts in S. lycopersicum and S. habrochaites, respectively, are cold regulated. There is a large overlap in the genes that were cold responsive in both plant species, but the results indicated many differences in the cold-responsive genes of the two species (Figs. 1, 2). The diversity of GO categories that were enriched in cold-stressed S. lycopersicum and S. habrochaites (Additional file 3) indicated the complexity of the response.
For cold-regulated DEGs of S. lycopersicum and S. habrochaites, some similar clusters of GO categories "response to abiotic stimulus" was found in both plants (Additional file 3), confirming earlier observations [20]. However, in response to cold stress in S. lycopersicum at 1 h, many genes encoding proteins associated with the abiotic stimulus response showed increased transcript abundance, and a few genes showed decreased transcript abundance (Additional file 3).
The data also suggested that some GO terms overlap in cold-treated S. habrochaites, but not in S. lycopersicum (Additional file 3). Some photosynthesis-related GO terms Fig. 6 The total number of targets of differentially expressed microRNAs (DEmiRNAs) that were either cold-induced or cold-repressed in S. lycopersicum (C) and S. habrochaites (Tsh) Fig. 5 The total number of microRNAs (miRNAs) that were either cold-induced or cold-repressed in S. lycopersicum (C) and S. habrochaites (Tsh) were significantly enriched among the upregulated genes in S. habrochaites under cold stress at 1 h. The results indicated that photosynthesis-related genes in cold regulatory programs might contribute to transient cold tolerance of S. habrochaites. In 2012, Liu et al. studied the transcriptomes of S. lycopersicum and S. habrochaites under cold stress at 3 days using GeneChip [20]. Their experiments indicated that S. lycopersicum showed more severe inhibition of photosynthesis than S. habrochaites during chilling stress.
The data showed that the GO category "cell wall metabolism" was enriched with genes with increased expression in cold-treated S. lycopersicum at 1 h, but the opposite result was observed for cold-treated S. habrochaites at both 1 and 12 h (Additional file 3). A GO term enrichment analysis indicated that "cell wall metabolism" was significantly depressed in the long term by cold stress in S. habrochaites, but was transiently induced by cold stress in S. lycopersicum. The present findings in S. habrochaites are similar to those of Fowler and Thomashow [15], who used a microarray chip to analyze the transcriptome of Arabidopsis under cold stress.
Hormones are signaling molecules that play key roles in regulating gene expression under cold stress [15,20]. RNA-seq analysis showed that many genes related to abscisic acid, ethylene, auxin, jasmonic acid, and gibberellin were regulated by cold stress in S. lycopersicum and S. habrochaites (Additional file 3). Two hormone-related GO terms, "response to abscisic acid stimulus" and "response to ethylene stimulus", were significantly enriched among the differentially expressed genes in S. lycopersicum and S. habrochaites under cold stress (Additional file 3). In S. lycopersicum, our data show that ABArelated GO terms were significantly enriched among the upregulated genes under cold stress at 12 h and were not enriched at 1 h. In contrast, ABA-related GO terms were significantly enriched among the upregulated and downregulated genes under cold stress at multiple time points in S. habrochaites. The production of ethylene also has been associated with cold stress [63][64][65]. In S. lycopersicum, our data show that ethylene-related GO terms were significantly enriched among the upregulated genes under cold stress at multiple time points. In contrast, ethylene-related GO terms were only significantly enriched among the upregulated genes under cold stress at 12 h, but not at 1 h, in S. habrochaites.
Transcription factors (TFs) play a key role in the regulation of gene expression under abiotic and biotic stresses in plants. The RNA-seq results showed that many TFs were regulated under cold stress in S. lycopersicum and S. habrochaites, and some members of the CBF were upregulated under cold stress at 1 h or 12 h ( Fig. 7; Additional file 3). Two additional cold-regulated genes were found that encode transcription factors: a homolog of Cys2/His2-type zinc-finger protein (ZAT10) and a zinc finger protein involved in high light and cold stress (ZAT12). Indeed, the RNA-seq results showed that the transcript levels of a MYC-type bHLH transcription factor (ICE1) increased under cold stress at 1 h, and decrease after 12 h.
AS of pre-mRNA is an important mechanism to increase transcriptome and proteome variation in eukaryotes. Previous examinations of AS events under abiotic stress showed that some AS events are stress related [34,66,67]. Different coverage of mRNA isoforms in RNA-seq was observed under abiotic stress, likely reflecting the regulation of AS events. Here, we used RNA-seq to identify the AS events in S. lycopersicum and S. habrochaites that differed under normal conditions and cold stress treatment. Compared with other methods, RNA-seq supplies a wide and deep sequencing of the transcriptome, providing experimental confirmation of splice junctions and AS events with low false-positive rates. Our data provide an exceptional and impartial evaluation of AS in S. lycopersicum and S. habrochaites. The results were similar to those of Filichkin et al. [34], who used A. thaliana RNA-seq data to compare the specific abiotic stress transcriptomes of A. thaliana. The authors identified 6,000 novel AS events within the introns of 3,120 genes.
Our RNA-seq analysis of the S. lycopersicum and S. habrochaites transcriptomes suggests numerous genes with AS may be associated with cold stress (Additional file 7). The expressions of 121 (sample C1 vs. C0), 522 (sample C12 vs. C0), 112 (sample Tsh1 vs. Tsh0), and 553 (sample Tsh12 vs. Tsh0) AS genes were increased under cold stress; and the expressions of 110 (sample C1 vs. C0), 140 (sample C12 vs. C0), 111 (sample Tsh1 vs. Tsh0), and 122 (sample Tsh12 vs. Tsh0) genes decreased (Fig. 4, Additional file 7). Most of the genes identified as undergoing cold-induced AS were transcripts whose levels remained constant under cold stress. Thus, despite the lack of change in transcript expression level, their coding ability could be very different.
Recently, miRNAs have been identified as new players in plant tolerance to abiotic stress, such as cold, heat, high salinity, drought, oxidative, hypoxia, and UV B [68]. Many studies have attempted to understand the roles of miRNAs in the response to cold in several plants, including A. thaliana [69], Brachypodium distachyon [70], Oryza sativa [71], and Populus trichocarpa [72]. In this study, sequencing was used to confirm the genome-wide miRNA expression patterns of S. lycopersicum and S. habrochaites under cold stress.
Some miRNAs in different plant species present different expression patterns under cold stress. For example, the expression of miR172 was inhibited after cold stress at 1 h in S. lycopersicum, but was induced in S. habrochaites after cold stress at 12 h (Additional file 14). Additionally, miR172 was upregulated in A. thaliana [69] and B.
distachyon [70]. Similarly, miR170/171 expression was upregulated under cold stress at 1 h in S. lycopersicum, while the transcript level of miR170/171 was decreased under cold stress at 12 h in S. habrochaites. The transcript level of miR170/171 was upregulated in A. thaliana [69] and downregulated in Oryza sativa [71] and Populus trichocarpa [72] in response to cold.

Conclusions
S. lycopersicum and S. habrochaites are closely related plant species, but their cold tolerances are different. In recent years, our research group has investigated the cold-tolerance mechanisms of these plants at the physiological and molecular levels. Here, we studied the transcriptomes of cold stressed leaves of S. lycopersicum and S. habrochaites. We obtained 68,051 assembled unigenes, and many cold-regulated genes were detected, representing useful resources for gene cloning to improve cold tolerance of crops. Furthermore, the comparison of the functional networks of cold-regulated genes in S. lycopersicum and S. habrochaites provided information that could help us to identify the differences Fig. 7 Diagram of cold-responsive transcriptional network in plant. Solid arrows indicate activation, whereas lines ending with a bar show negative regulation. Abbreviations: CBF, C-repeat binding factor (an AP2-type transcription factor); ICE1, inducer of CBF expression 1 (a MYC-type bHLH transcription factor); LOS2, low expression of osmotically responsive genes 2 (a bifunctional enolase with transcriptional repression activity); ZAT12, a zinc finger protein involved in high light and cold acclimation; ZAT10, related to Cys2/His2-type zinc-finger proteins found in higher plants; KIN1, encodes protein kinase APK2a; ERD10, encodes a gene induced by low temperature and dehydration in cold-tolerance mechanisms between S. lycopersicum and S. habrochaites. We found that 21 % and 23 % of genes were differentially expressed between the cultivated and wild tomato species, respectively, when plants were transferred from warm to cold temperatures. An AS analysis suggested that the relative abundance of isoforms of S. lycopersicum and S. habrochaites significantly shifted under cold stress. In addition, certain miRNAs (e.g., miR159, miR319, and miR6022) play roles in the response to cold stress. Thus, differences in cold regulatory mechanisms may contribute to the differences in cold tolerance of these two tomato species.

Methods
Plant material and cold stress conditions S. habrochaites LA1777 was supplied by Tomato Genetics Research Center (University of California, Davis, USA). S. lycopersicum 'glamor' and S. habrochaites were grow at 25°C with 16-h light and 8-h dark cycles for 8 weeks before harvesting. To avoid changes caused by the circadian rhythm, all cold stress treatments were started at 4°C at 12 PM under light and continued for 0 (untreated control), 1, and 12 h.

Physiological responses to cold stress
The MDA content was assayed as described by Campos [73]. The free proline content was determined according to the method described by Zhang et al., [70]. POD activities were determined as described by Quiroga [74]. CAT activity was assayed as described by Yao [75].

Total RNA extraction and library preparation
The total RNA from leaves was extracted using the TRIzol reagent (Invitrogen) and digested with RQ1 DNase (Promega) to remove genomic DNA. The quality and integrity of the total RNA were detected using a Smart-Spec plus Spectrophotometer (Bio-Rad) and 1.5 % agarose gel electrophoresis. Polyadenylated mRNAs were purified and concentrated using oligo(dT)-conjugated magnetic beads before being used for directional RNAseq library preparation. Purified mRNAs were fragmented at 95°C and subjected to end repair and 5′ adaptor ligation, followed by reverse transcription using randomized hexamers and an RT primer with a 3′ adaptor. Purified cDNAs were amplified, and 200-500 bp PCR products were quantified and purified. RNA-seq libraries were prepared and applied to an Illumina Genome Analyzer IIx system for 80 nt single-end sequencing by ABlife Inc. (Wuhan, China).
For small RNAs, 3 μg of total RNA was used for small RNA library preparation using the Balancer NGS Library Preparation Kit for small/microRNA (GnomeGen), following the manufacturer's instructions. The purified small RNA libraries were quantified using a Qubit Fluorometer (Invitrogen) and used for cluster generation and 36 nt single-end sequencing analysis using the Illumina Genome Analyzer IIx system.
Processing, mapping of Illumina reads, detection of alternative splicing, and de novo assembly The adapter sequences and low-quality bases at the 3′ ends were removed from the RNA-seq reads generated by the Illumina Genome Analyzer. Reads whose lengths were more than 20 bp were used for further analysis. The reads were mapped to the tomato genome using TopHat [59], which allows the confirmation of AS events. In total, 172,910 junction sites were identified. We categorized the AS events into different types according to the exon structures, using the ABLas-1 package (ABlife Inc., Wuhan, China). These categories included exon skipping, intron retention, alternative 5′ splice site, alternative 3′ splice site, and mutually exclusive exons, as described by Wang et al. [76]. Reads were assembled separately from each S. lycopersicoides library using the Trinity method [58]. There are three software modules in Trinity: Inchworm, Chrysalis, and Butterfly, which were used to process the RNA-seq reads sequentially. First, the Inchworm program assembled the reads to contigs. Second, the Chrysalis program clustered the minimal overlapping contigs. Third, the Butterfly program constructed the transcripts. Finally, the multiple sequence alignment tool BLAST was used to cluster the transcripts by similarity of the right match length [49]. The coding sequences (CDS) of the unigenes were predicted using EMBOSS (http://emboss.sourceforge.net/) and the longest CDS was considered as the complete CDS of a unigene.

Differential expression and GO enrichment
The expression level of genes were evaluated and normalized using the reads per kilobase of transcript per million reads mapped method [77]. Unigene expression was analyzed using the Bioconductor package with the edgeR and Bayseq methods. A very stringent cutoff, normalized fold change ≥2 or ≤0.5 and a P-value (P ≤0.01), was used to identify cold-regulated DEGs. A GO analysis was employed to predict gene function and calculate the functional distribution frequency, using the DAVID package (ABlife Inc., Wuhan, China).

qPCR analysis
To validate the transcript abundance of genes measured by RNA-seq, we performed qPCR using Power SYBR Green Mastermix in an Applied Biosystems 7500 Real-Time PCR System. The RNAs from S. habrochaites and S. lycopersicum used in RNA-seq were reverse transcribed into cDNAs. The 14 primer pairs used are listed in Additional file 2. The ACTIN gene was used as a reference in these experiments. Three technical replicates were used for qPCR. The single amplicons were confirmed by melting curve analysis and gel electrophoresis of the final product. The cycle threshold (CT) value of each gene was normalized to the reference gene to detect the relative fold changes in each sample, which was calculated using the ΔΔCT method, as described previously [78].

Identification of miRNAs in tomato
High-quality small RNA reads ranging from 14 to 24 nucleotides were acquired from the raw data. Adaptor sequences and low-quality tags were removed to detect known and novel miRNAs in tomato. Small RNA reads were used to search the Rfam database and NCBI database to remove non-coding RNAs, such as rRNA, tRNA, snRNA, and snoRNA. The remaining sequences were searched in the miRBase database v19.0, with no mismatch permitted, to identify conserved mature miRNA orthologs. Small RNAs that did not map to any miRNAs in miRBase database were analyzed as novel miRNAs using miRDeep2 (developed by ABlife Inc.).

Availability of supporting data
The raw RNA-seq data supporting the result of this article is available in the Sequence Read Archive (SRA), with accession numbers SRX1013429 and SRX1014317.