Transcriptome differences between two sister desert poplar species under salt stress

Background Populus euphratica Oliv and P. pruinosa Schrenk (Salicaceae) both grow in dry desert areas with high summer temperatures. However, P. euphratica is distributed in dry deserts with deep underground water whereas P. pruinosa occurs in deserts in which there is underground water close to the surface. We therefore hypothesized that these two sister species may have evolved divergent regulatory and metabolic pathways during their interaction with different salt habitats and other stresses. To test this hypothesis, we compared transcriptomes from callus exposed to 24 h of salt stress and control callus samples from both species and identified differentially expressed genes (DEGs) and alternative splicing (AS) events that had occurred under salt stress. Results A total of 36,144 transcripts were identified and 1430 genes were found to be differentially expressed in at least one species in response to salt stress. Of these DEGs, 884 and 860 were identified in P. euphratica and P. pruinosa, respectively, while 314 DEGs were common to both species. On the basis of parametric analysis of gene set enrichment, GO enrichment in P. euphratica was found to be significantly different from that in P. pruinosa. Numerous genes involved in hormone biosynthesis, transporters and transcription factors showed clear differences between the two species in response to salt stress. We also identified 26,560 AS events which were mapped to 8380 poplar genomic loci from four libraries. GO enrichments for genes undergoing AS events in P. euphratica differed significantly from those in P. pruinosa. Conclusions A number of salt-responsive genes in both P. euphratica and P. pruinosa were identified and candidate genes with potential roles in the salinity adaptation were proposed. Transcriptome comparisons of two sister desert poplar species under salt stress suggest that these two species may have developed different genetic pathways in order to adapt to different desert salt habitats. The DEGs that were found to be common to both species under salt stress may be especially important for future genetic improvement of cultivated poplars or other crops through transgenic approaches in order to increase tolerance of saline soil conditions. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-337) contains supplementary material, which is available to authorized users.


Background
Salinity and drought stresses are the two most important environmental factors limiting plant growth and development in semiarid and arid areas [1]. Over 100 countries in the world have been identified as being affected by salinity [2], and the scale of the problem seems to be increasing at an alarming rate [3]. Salinity, together with drought, has far-reaching implications for food security, economic sustainability and the irreplaceable biodiversity of any affected area, and it is anticipated that these challenges will be exacerbated by the projected impact of climate change. The effects of water-insufficiency stresses have been studied extensively; they limit water and micronutrient uptake and lead to closure of stomata, decline in carbon metabolism, stunted growth, ion/salt toxicity and reduced yield [3,4].
For plants to survive under such conditions, they must sense and respond to these abiotic stresses rapidly and in a complex manner [5], through signalling and regulatory pathways [3,4,6] mediated by abscisic acid [7] or ethylene [8], generally resulting in altered expression of transcription factors [9], and in many cases in increased expression of genes encoding products required for osmoregulation, cell protection and/or acclimation [10][11][12][13][14][15]. These modifications may lead to changes in signal transduction, ionic homeostasis, scavenging of reactive oxygen species, accumulation of compatible solutes and growth regulation [3,6,[16][17][18]. A common strategy for the identification of overall changes in gene expression under salt stress is to compare the transcriptomes of the targeted species or cultivars using microarrays and/or RNA-Seq technologies [19]. A plethora of comparisons between salt-sensitive and salt-tolerant cultivars of model and non-model plant species, including Arabidopsis [20][21][22], rice [23], poplar [24][25][26][27], tomato [28], potato [29], Medicago truncatula [30], sugarcane [31] and olive [32], have been reported to date. These studies have identified more than 30 families of transcription factors and numerous enzyme-encoding genes involved in responses to salt stress [33,34]. However, overall changes in gene expression and physiological responses to salt stress vary greatly between different species, particularly between sensitive and non-sensitive pairs of related species [35][36][37][38][39]. It is often difficult to ascertain whether these differences were caused by divergence during the course of evolution or were brought about through adaptive differentiation. It is therefore of interest to compare the overall changes in gene expression that occur in sister species under salt stress, as this will minimise phylogenetic effects.
Here we examine differences in the transcriptomes of two sister desert poplar species under salt stress. Populus serves as a model for elucidating physiological and molecular mechanisms of stress tolerance in tree species [40][41][42]. Both P. euphratica and P. pruinosa grow in dry deserts with high summer temperatures [43][44][45][46]. Both species can tolerate high salinity and survive NaCl concentrations of more than 300 mM [47] in nutrient solution, and P. euphratica has been used as a model species for studying abiotic responses to salt or drought stress [27,[48][49][50]. In addition to differences in leaf and hair morphology between the two species, they also occur in different types of habitat. P. euphratica is found in dry deserts with deep underground water while P. pruinosa is distributed in deserts where the underground water is closer to the surface, and therefore more accessible, but also saltier near ancient or extant rivers. It is likely that these two species have diverged due to ecological differentiation, in spite of ongoing gene flow [46].
In order to test whether regulatory and metabolic pathways in these two species have diverged during their adaptive interactions with salt and other stresses, the transcriptomes of callus subjected to 24 h of salt stress, and control callus samples, from P. euphratica and P. pruinosa were compared in order to identify differentially expressed genes (DEGs) and alternative splicing (AS) events that occurred in response to salt stress. Our results revealed that these two poplar species have both common and species-specific patterns of gene expression under salt stress. The dynamic transcriptome expression profiles of these sister species under salt stress obtained in this study may provide useful insights to inform further analyses of the mechanism of high salinity tolerance in plants. In addition, the genes found to be differentially expressed under salt stress in both species may facilitate the identification of key genes as potentially suitable targets for biotechnological manipulation with the aim of improving poplar salt tolerance.

Results and discussion
Analysis and mapping of Illumina-Solexa sequencing tags We used the Illumina-Solexa sequencing platform to sequence the P. euphratica [27] and P. pruinosa [51] transcriptomes obtained from the four treatments, including two unstressed callus samples as controls (P. euphratica control callus, PeuC; P. pruinosa control callus, PprC) and two salt-stressed callus samples as treatments (P. euphratica salt-stressed callus, PeuS; P. pruinosa salt-stressed callus, PprS). After removing low-quality sequences and trimming adapter sequences,~28 million 75-bp paired-end clean reads were generated from each of the cDNA libraries in the Illumina Genome Analyzer runs (Table 1). These tags from the four digital gene expression (DGE) libraries were mapped to the available P. trichocarpa transcript sequences. Approximately 80% of the tags had matches. Most (79.2-82.4%) of the tags with matches were unique tags (matching only one poplar locus), while the remainder (~17.6-20.8%) were non-unique (matching more than one poplar locus) or unaligned. For more detailed investigation of gene expression in the different treatments, only unique tags were used in the analysis. In total, 36,144 transcripts were identified from the four conditions. The transcripts identified accounted for 80.3% of the 45,033 annotated genes in poplar. In both control and salt stress treatments, the numbers of mapped genes in P. euphratica (33,528 and 32,508 genes) were found to be similar to those in P. pruinosa (32,996 and 33,055 genes, respectively) ( Table 1). We further compared the mapped genes among the four treatments (PeuC, PeuS, PprC and PprS), and found that~89.1% of them were present in at least two treatments ( Figure 1).

DEGs in the two species under salt stress
To identify global transcriptional changes occurring under salt stress, we applied four independent metrics to identify genes that were differentially expressed between the 24-h salt-stressed callus and control callus samples in P. euphratica and P. pruinosa. For each metric, we selected those DEGs whose expression profiles met three criteria: (i) the FPKM value was ≥1 in either of the libraries, (ii) log 2 (FPKM salt /FPKM control ) was > 1 or < −1, and (iii) the adjusted p-value (FDR) was < 0.05. In this study, DEGs with higher expression levels in salt-stressed callus when compared with control callus samples were termed 'up-regulated' , while those with lower expression levels in saltstressed callus were termed 'down-regulated'. There were 471 and 593 genes identified by all metrics as being upregulated in P. euphratica and P. pruinosa, respectively, and 413 and 267 genes identified by all metrics as downregulated in P. euphratica and P. pruinosa, respectively ( Figure 2). There were more up-regulated DEGs in P. pruinosa than in P. euphratica, while there were more downregulated DEGs in P. euphratica than in P. pruinosa.
The DEGs identified were classified into eight clusters according to their expression patterns ( Figure 3, Additional file 1). Of these eight clusters, four were up-regulated or down-regulated exclusively in a single species, as follows: up-regulated exclusively in P. euphratica (272 DEGs) or in P. pruinosa (394); down-regulated exclusively in P. euphratica (298) or in P. pruinosa (152). The remaining four clusters consisted of genes that were up-or down-regulated in the two species; two of these clusters showed similar co-regulation patterns whereas the other two showed opposing regulation patterns. In the two clusters with similar co-regulation patterns, 198 DEGs were co-up-regulated and 114 DEGs were co-downregulated in the two species. Within the co-up-regulated clusters, only one transcript (POPTR_0013s12880.1) was undetectable in the calli of the two species under unstressed conditions (Additional file 1), suggesting that this gene is expressed specifically under salt stress in both species. In the two clusters with opposing patterns of regulation, only 1 DEG was up-regulated in P. euphratica but down-regulated in P. pruinosa, and only 1 DEG was downregulated in P. euphratica but up-regulated in P. pruinosa. This result suggested that our integrated DEG identification was sensitive and reliable.

Confirmation of differentially expressed candidate genes by qRT-PCR analysis
To confirm the gene expression inferred from RNA-seq, a total of 21 candidate DEGs with salt-related process were selected for the qRT-PCR analyses, comprising 7 DEGs exclusively regulated in a single species, 8 co-upregulated and 6 co-down-regulated in the two species ( Figure 4). Although the exact change did not exactly match each other, the expression trends of all 21 genes from qRT-PCR and Illumina-Solexa RNA sequencing analyses were largely consistent (Pearson's correlation coefficient r = 0.8), demonstrating the reliability of the RNA-seq results ( Figure 4).

Gene functional categories of two species under salt stress
Firstly, an overview of the main results was obtained by WEGO and the DEGs were assigned to GO terms in the three component ontologies ( Figure 5). Then, groups of genes with functions involved in salt responses were identified using parametric analysis of gene set enrichment (PAGE) ( Table 2). GO enrichment in P. euphratica was significantly different from that in P. pruinosa. In the Cellular Component ontology, 'apoplast' (GO:0044464) appeared to respond to salt stress in both species; while 'cell part' (GO:0044464) and 'cell' (GO:0005623) were enriched only in P. euphratica; whereas 'extracellular region' (GO:0005576), 'external encapsulating structure' (GO:0030312) and 'cell wall' (GO:0005618) were enriched only in P. pruinosa. In the Molecular Function ontology, 'cofactor binding' (GO:0048037), 'coenzyme binding' (GO:0050662), 'peptidase inhibitor activity' (GO:0030414) and 'endopeptidase inhibitor activity' (GO:0004866) were enriched in both species, while another nine terms from  the Molecular Function ontology were enriched exclusively in P. euphratica. Six terms from the Biological Processes ontology were enriched exclusively in P. euphratica and three terms were enriched in both species. The GO terms enriched in P. euphratica were related to responses to stress and metabolic processes, and the most highly enriched term was 'response to stress' (GO:0006950). We also used singular enrichment analysis (SEA) to identify functional groups of genes differentially expressed in the two species under salinity (Additional file 2). GO enrichment for genes  up-regulated or down-regulated exclusively in P. euphratica was significantly different from that in P. pruinosa. The detected differences suggested that these two desert poplars might have developed different genetic pathways for adaptation to differentiated salty desert habitats.

Differences in expression of hormone-related genes in the two species under salt stress
Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database as our source of annotations, 583 out of 803 Populus genes annotated as being involved in  Gene ontology (GO) annotation of salt-responsive genes compared between P. euphratica and P. pruinosa. WEGO was used to produce the graph. We divided the sets into the three major GO domains: biological process, cellular component and molecular function, and the number (right y-axis) and percentage (left y-axis) of genes were calculated.
hormone biosynthesis [52] were detected in the four libraries and 59 of these genes were differentially expressed in either P. euphratica or P. pruinosa under salt stress (Additional file 3). Among these hormone biosynthesisrelated DEGs, 37 were identified in P. euphratica, of which 36 were up-regulated and one was down-regulated during salt stress, while 39 were identified in P. pruinosa, including 36 up-regulated and 3 down-regulated DEGs; only 17 were co-regulated in the two species. Hierarchical clustering of the 59 hormone-related DEGs showed overall differences between P. euphratica and P. pruinosa in response to salt stress ( Figure 6). Under salt stress, most hormone-related DEGs were co-up-regulated or co-down-regulated in both species. Interestingly, three ABA metabolism-related genes (ABA1, POPTR_0007s10980.1; and two genes encoding 9-cis-epoxycarotenoid dioxygenase (NCEDs), POPTR_ 0011s11370.1 and POPTR_0001s40420.1) were up-regulated exclusively in P. pruinosa. The 1430 DEGs were analyzed by parametric analysis of gene set enrichment (PAGE). 2 Onto, Gene Ontology domain; P, biological process; F, molecular function; C, cellular component. 3 Mean, mean log 2 expression ratio, >0 represents up-regulation, <0 represents down-regulation.  CHX18), potassium transporter 6 (AT1G70300, KUP6) and ABC transporter (AT1G66950, ABCG39), respectively, were co-up-regulated in both species. However, POPTR_ 0005s04660.1 and POPTR_0014s12700.1, which are homologous to Arabidopsis thaliana sodium/hydrogen exchanger 2 (AT3G05030, NHX2) and potassium transporter 5 (AT4G13420, HAK5) genes, were up- regulated only in P. euphratica. In contrast, POPTR_ 0004s23680.1 and POPTR_0013s08110.1, which are homologous to Arabidopsis thaliana chloride channel protein CLC-c (AT5G49890, CLC-C) and potassium transporter 2 (AT2G40540, KT2) genes respectively, were up-regulated exclusively in P. pruinosa. These results corroborate previous findings [27,[53][54][55], and confirm that genes encoding proteins such as sodium and potassium ion transmembrane transporters, and chloride channel and ABC transporters, which are important for maintaining and re-establishing homeostasis in the cytoplasm, are induced to high levels in response to salinity stress [16].

Differences in expression of transcription factor genes in the two species under salt stress
We identified 4016 transcription factors in Populus trichocarpa and classified them into 92 families (Additional file 5) based on published annotations. A total of 115 genes that were differentially regulated in either P. euphratica or P. pruinosa during salt stress were categorized as transcription factors (Additional file 6). Of these, 59 DEGs were identified in P. euphratica, including 24 that were upregulated and 35 down-regulated during salt stress, while 73 DEGs were identified in P. pruinosa, of which 52 were up-regulated and 21 down-regulated. Only 17 DEGs were co-regulated in both species (Table 3, Additional file 6). Several of the transcription factors, such as AP2/ERF and bZIP, which are known to be induced by stress in model plant species (Arabidopsis thaliana and rice) [56,57], were highly expressed in response to salinity stress in P. euphratica or P. pruinosa.

The co-up-regulated DEGs in the two species under salt stress and allele mining
A total of 198 co-up-regulated DEGs in the two species were identified in salt stress (Additional file 1) and the important ones were selected and listed in Table 4. The candidate genes identified in the present study contained both the previously reported salt-responsive genes and some species-specific ones. Of them, most genes were involved in and highly enriched in functional categories such as response to stress, signal transduction, transmembrane transport, transcriptional regulation and basic metabolic processes (Additional files 1 and 2). These findings are beneficial to allele mining of two poplar species related to their common or differentiated response to stressed habitats in the future. Allele mining based on the candidate genes were found to be important in dissecting naturally selected allelic variations that controlled differentiated traits [58,59]. In addition, promoters are found to play a key role in gene regulation, and any change in these regions will change gene expression and the controlled traits. Therefore, the identified variations through such an approach may be mainly located in the promoter regions [60]. Overall, the co-up-regulated DEGs identified in the present study provide critical genetic bases for further allele mining, functional analyses and transgenic practices for developing the salt-tolerant poplars and crops.

A comparison of DEGs identified by our results and other transcriptome studies of the salt-stressed poplars
In order to test the consistency of DEGs across different treatments and approaches, we compared DEGs between our results and other available transcriptome studies of the salt stressed poplars. Ottow et al. [48] examined changes in transcript levels of various genes known to be involved in salt or general stress signaling or adaptation in P. euphratica leaves by dot-blot expression. They identified nine genes with significant changes in response to salt stress. Some of them were be confirmed in the present study, for example, galactinol synthase 2 (GolS2, POPTR_ 0013s00730.1), calcineurin B-like protein 4 (CBL 4, POPTR_0015s01550.1), alternative oxidase 1A (POPTR_ 0012s01630.1) and 1-aminocyclopropane-1-carboxylate oxidase (POPTR_0011s00970.1) (Additional file 1). Galactinol synthase (GolS) catalyzes the first step in the biosynthetic pathway of raffinose oligosaccharides using galactose and myo-inositol as substrates and this gene was also up-regulated in plants under cold, heat, drought, and salt stress [21,61,62]. Significant increases in galactinol synthase and alternative oxidase after salt stress point to shifts in carbohydrate metabolism and suppression of reactive oxygen species in mitochondria under salt stress [48]. In addition, Gu et al. [63] identified 54 genes with altered transcript accumulation in the salt-stressed P. euphratica by microarray hybridization. The genes of them, responsible for hydroxyproline-rich glycoprotein, carbonic anhydrase 2, cytochrome P450, aquaporin, sucrose synthase and aspartate aminotransferase were confirmed in these present study. These genes were also revealed to be salt-responsive in other studies [26,27,64]. The drought responses of plants are similar to those in response to salinity because both stresses lead to physiological water deficit [65]. Bogeat-Triboulot et al. [66] provided a comprehensive analysis of P. euphratica subjected to gradual soil water depletion, and observed 110 regulatory and protective genes involved in long-term response to drought. Similar results were also found by Cohen et al. [25] and Tang et al. [67]. Among them, those genes involved in metabolites of proline, raffinose, galactose, inositol and sucrose under drought stress were found to have changed their expressions in response to salt stress in the presnet study. An increase in galactinol, raffinose and stachyose content may have improved osmoprotection and ROS scavenging when poplars were stressed by drought or salt. However, in the present study, we identified numerous more transcripts with significant up-regulations in both poplars when stressed, including UDP-glycosyltransferase-like protein, FADbinding and BBE domain-containing protein, putative nucleoredoxin 1, and glyceraldehyde-3-phosphate dehydrogenase. All these newly identified genes should have also played an important role during salt adaptation of two species. Their functions and molecular mechanisms need further clarifications in the future.

Alternative splicing of transcripts in the two species under salt stress
Finally, to investigate the role of alternative splicing (AS) in response to salt stress, we conducted a survey of transcript isoforms across the four libraries and examined six common types of 'alternative splicing events'. For each of these event types, reads deriving from specific regions can be used to identify the expression of one alternative isoform or the other (Table 5). We identified 26,560 AS events which were mapped to the 8380 poplar genomic loci from the four libraries, suggesting that 20% of 40,668 loci in poplar are potentially subject to AS. This observed AS percentage is comparable to the percentage of the genes shown to undergo AS in A. thaliana (21.8%) and rice (21.2%) [68]. In both control and salt stress treatments, the number of genes exhibiting AS events in P. euphratica (6662 and 5850 genes) is similar to that in P. pruinosa (6192 and 5765 genes), respectively. In addition to the AS loci (4115) common to both species, around 346 and 243 of the loci show AS events only in either P. euphratica or P. pruinosa in response to salt stress (Figure 7). We further classified those genes that underwent AS events on the basis of functional ontology. GO enrichment for the genes displaying AS events exclusively in P. euphratica was significantly different from that for genes undergoing AS in P. pruinosa (Additional file 7).

Conclusions
Our transcriptional profiling analysis revealed numerous genes that were differentially expressed in both P. euphratica and P. pruinosa under salt stress. The differential expressions of the selected genes inferred from RNA-seq were confirmed by qRT-PCR data. Gene ontology analyses of these DEGs suggested that GO enrichment in P. euphratica was significantly different from that in P. pruinosa. We found that numerous genes involved in hormone biosynthesis, or encoding transporters or transcription factors, showed different expression patterns between these two species under salt stress. These differences suggest that these two desert poplars may have developed species-specific pathways for adaptation to salinity during the course of ecological speciation in their different salty desert habitats. The results of our comparative analyses imply that different species, even sister species, may employ different genetic pathways to cope with salt stress. This suggests that it may be more difficult than previously anticipated to design salt-tolerant plant cultivars [69,70]. In order to develop cultivars with high salt tolerance, particular attention should be paid to those genes that are differentially expressed in two or more different species under salt stress. Such genes can be used to facilitate genetic improvement of crops, including cultivated poplars, for growth on saline soils.

Gene expression data
Paired-end RNA-seq reads for control callus and saltstressed callus of P. euphratica and P. pruinosa, which were obtained by Qiu et al. [27] and Zhang et al. [51], respectively, were downloaded from the NCBI sequence read archive (accession numbers SRX025571, SRX025568, SRX245887 and SRX245885). We cultured P. euphratica and P. pruinosa calli induced from the shoot under the same conditions. We then replaced the growth medium for one set with the fresh medium and the same medium but supplemented with 100 mM NaCl (salt stress) for another set. We harvested both sets of calli 24 h later. The calli from P. euphratica and P. pruinosa had the same subculture generation and time and they were highly comparable in terms of physiological state. After RNA extraction and quality determination, we constructed the paired-end cDNA libraries with insert sizes of 200 base pair (bp), and then sequenced the cDNA using an Illumina (San Diego, CA, USA) Genome Analyzer platform according to the manufacturer's protocols with a read length of 75 bp in two lanes. Image output data from the sequencer was transformed into raw sequence data by base calling.
Raw reads generated by Illumina Genome Analyzer were initially processed to obtain clean reads. We first cleaned raw sequence reads by removing exact duplicates from both sequencing directions. We further cleaned reads by removing adapter sequences as well as reads with too many (>8) unknown base calls (N), low complexity, and low-quality bases (>50% of the bases with a quality score ≤5). Cleaned reads from each library were used for later differential expression analysis in this study.

Initial mapping of reads
To determine the level of gene expression, Bowtie2 [71] was used to align RNA-seq reads from the control and salt-stressed samples to transcript sequences from Populus trichocarpa Torr. & A. Gray [41], using annotation files downloaded from http://www.phytozome.net/poplar (JGI Populus trichocarpa v2.2). No more than a 1 bp mismatch was allowed when taking into account differences between the two species. Reads that mapped to reference sequences from multiple genes were filtered out. The remaining clean reads, which were considered to be distinct, were used for further analysis. Transcript abundances were calculated using eXpress [72], which outputs read counts and the number of fragments per kilobase of exon per million fragments mapped (FPKM) [73]. Transcripts with FPKM values < 1 in both libraries were filtered out and not subjected to further analysis.

Identification of differentially expressed genes
To identify differentially expressed genes (DEGs) in control callus and salt-stressed callus from P. euphratica and P. pruinosa, we applied four independent, widely used tools: Cuffdiff [73], DESeq [74], edgeR [75], and EBSeq [76]. Cuffdiff takes a nonparametric, annotation-guided approach to estimating the means and variances of transcript FPKM values under different conditions, using Student's t-tests to identify differentially expressed transcripts [73]. In contrast, DESeq, edgeR and EBSeq estimate the means and variances of raw read counts under a negative binomial distribution and use exact tests to identify differentially expressed transcripts. The main difference between DESeq, edgeR and EBSeq is that they use different statistical approaches to estimate variance [74][75][76]. After the p-values for each expressed genes were obtained by the four tools, the false discovery rate (FDR) was used to justify the p-value by the function p.adjust in R. Sequences were deemed to be differentially expressed if log 2 (FPKM salt /FPKM control ) > 1 or < −1, and the adjusted p-value (FDR) was < 0.05 as identified by all four metrics.

Functional annotation through BLAST2GO and KEGG
Gene Ontology (GO) terms were assigned to the identified genes by the blast2GO pipeline [77] using NCBI databases, followed by functional classification using the WEGO software package [78]. For the comparative analysis of DEGs between P. euphratica and P. pruinosa in response to salinity, singular enrichment analysis (SEA) and parametric analysis of gene set enrichment (PAGE) were performed using the agriGO program (http://bioinfo.cau. edu.cn/agriGO) [79] with the default parameters, using the P. trichocarpa gene models as background, followed by multiple testing with Bonferroni correction (corrected P-value < 0.05). PermutMatrix (Version 1.9.3; http://www. lirmm.fr/~caraux/PermutMatrix/index.html) was used to cluster genes related to plant hormone biosynthesis according to their mean normalized intensity values [80].

Validation of DEG Expression with Quantitative Real-time PCR (qRT-PCR)
In order to validate the reliability of RNA-Seq experiments, a total of 21 candidate DEGs highly related to salt stress were selected for qRT-PCR test. These genes were chosen for the qRT-PCR analysis based on two criteria: (i) gene's expression patterns between these two species under salt stress should be similar; (ii) it should have only one BLAST hit when searching against genes of Arabidopsis thaliana to exclude paralogs. A total of 0.5 μg of DNase I-treated total RNA was converted into single-stranded cDNA using a Prime-Script 1st Strand cDNA Synthesis Kit (TaKaRa, Dalian, China). The cDNA templates were then diluted 20-fold before use. The quantitative reaction was performed on a CFX96 Real-Time PCR Detection System (Bio-Rad, Singapore) using SYBR Premix Ex Taq™ (TaKaRa, Dalian, China). PCR amplification was performed under the following conditions: 30 s at 95°C, followed by 40 cycles of 95°C for 15 s, 60°C for 30 s and then 72°C for 20 s. All primers were designed using PRIMER3 software and were listed in Additional file 8. Three biological replicates based the calli cultured from different individuals with the same subculture and physiological state were performed in order to exclude sampling errors. The relative expression levels of the selected DEGs normalized to an internal reference gene actin was calculated using 2 -ΔΔCt method [81].

Identification of alternative splicing
We prepared a database of all possible splice junctions between annotated exons in each selected gene and identified new possible junctions using TopHat [82]. We combined these two databases, removing any redundancy between them, and then extracted junction sequences of width 65 bases on each side from all the above junctions. To evaluate which of these junctions were validated by our Illumina reads, we aligned reads from each library separately against the junction sequences, allowing up to one mismatch (in a read of 75 bp). If at least two reads aligned to a splice junction, we considered it to be validated. Six different types of alternative mRNA processing events were analysed [83]. We first considered skipped exons (SE), in which one or more exons are spliced out of the mature message, and retained introns (RI), in which one or more introns are included in the message. Also included were alternative 5' splice site (A5SS) and alternative 3' splice site (A3SS) events, which are particularly difficult to interrogate by microarray analysis because the variably included region is often quite small. Finally, alternative last exons (ALEs) in which alternative use of a pair of polyadenylation sites results in distinct terminal exons, and alternative first exons (AFEs), where alternative promoter use results in mRNA isoforms with distinct 5' UTRs, were considered.