Transcriptomic profiling of the salt-stress response in the wild recretohalophyte Reaumuria trigyna

1Soil salinity is one of the most significant abiotic stresses limiting plant growth, productivity, and geographical distribution. Approximately 20% of the world’s cultivated lands, and nearly 50% of irrigated lands, are affected by salinity. High salt levels impose two stresses on plants: osmotic stress, resulting from the lowered water availability because of the high osmotic pressure in the soil; and ionic stress arising from solute imbalance, in which there are increased levels of Na+ and Cl- in the cytosol and changes to the intracellular K+/Na+ ratio. Halophytes are the natural inhabitants of highly saline soils, and have evolved resistance strategies including efficient control of the uptake and compartmentalization of salt ions, synthesis of organic ‘compatible’ solutes, and unique morphological structures such as succulent leaves, salt glands, and bladders. Such remarkable adaptations to salty soil make halophytes excellent candidates to investigate salt-tolerance mechanisms and to identify effective salt-response genes in this type of flora.

2Studies on the salt-tolerance mechanisms of halophytes have been carried out recently in several model plants, such as Thellungiella halophila, Mesembryanthemum crystallinum, Suaeda, and Populus. Excellent progress has been made in understanding the physiological and molecular mechanisms underlying stress tolerance. However, knowledge about halophytes and their physiology is rather limited with respect to their wide taxonomic occurrence. It remains unclear whether halophytes have unique tolerance mechanisms, and whether such mechanisms are taxonomically linked, and/or have evolved to respond to interactions between salinity and other environmental variables. Therefore, studies on several species from different salt-stressed environments should be performed to answer these questions.

3Reaumuria trigyna (Reaumuria Linn genus, family Tamaricaceae) is an endangered dicotyledonous shrub with the features of a recretohalophyte. It is regarded as a living fossil owing to its Tethys Ocean origin. This species is endemic to the Eastern Alxa–Western Ordos area (106°27’E–111°28’E, 39°13’N–40°52’N, elevation 1500–2100 m), a salinized desert in Inner Mongolia, China. The area is characterized by its high soil salinity (up to 0.7% salts), hyper-drought conditions (annual average precipitation of 140.9–302.2 mm), and low temperature (annual average temperature of 6.0–9.2°C). To adapt to the semi-desert and salty soil environment, R. trigyna has developed a distinct morphology that is characterized by succulent leaves and sunken stomata. These features are typical of a recretohalophyte—a halophyte with salt excretion glands. According to our previous studies, the concentrations of Na+, Cl-, and SO42- in the soil where R. trigyna grows are higher than those of K+, Mg2+, and Ca2+. The Na+ concentration in the soil is as high as 0.6 mg/g, with a Na+/K+ ratio of approximately 16. R. trigyna excretes Na+ and Cl- via its multicellular salt glands, but little K+ is excreted, with a Na+/K+ ratio of only 0.44. When seedlings were treated with 400 mM NaCl, R. trigyna showed a rapid increase in glutathione (GSH) from nearly 4.77 to 9.62 nmol·g-1 FW, and showed significant increases in the activities of superoxide dismutase (SOD) and peroxidase (POD) (P < 0.05). As a typical recretohalophyte growing in this area, R. trigyna could be considered as a reference species to represent those that survive in the severely adverse environment of the Ordos desert. Knowledge about the unique adaption strategies employed by R. trigyna will provide valuable references to the salt-tolerance mechanisms in these species and possibly in other halophytes growing in different environments. However, these objectives are difficult to achieve without detailed genetic and sequence information for this species. Fortunately, powerful high throughput sequencing technologies are now available to sequence the genome of plants de novo. Therefore, it is now possible to obtain the required sequences of R. trigyna to address important questions about its physiology, such as its mechanisms of salt-tolerance.

4Next-generation sequencing (NGS) technology has developed rapidly in recent years. It not only provides rapid, cost-effective, and comprehensive analyses of complex nucleic acid populations for model plants or species closely related to model plants, but also provides opportunities to analyze non-model plants whose genomes have never been sequenced. This technology has been used widely in comparative transcriptomics to identify differences in transcript abundance among different cultivars, organs, and different treatment conditions. In the present study, we generated a transcriptome dataset to provide genetic information to explore the salt-tolerance mechanisms of R. trigyna using the Illumina HiSeq™ 2000 platform. We compared transcriptomes of salt-stressed and control plants to identify genes showing transcriptional changes, and identified the functions of the transcripts and the KEGG pathways that showed changes. We speculate that the assembled, annotated transcriptome sequences and transcript abundance patterns will provide a valuable genetic resource for further investigations of the molecular mechanisms of salt tolerance in this species, and possibly in other recretohalophytes.

5Two sequencing libraries were prepared from control (C21) and NaCl-treated samples (T43) to investigate the transcriptomic responses to salt-stress in R. trigyna. In total, 26.51 million raw reads were generated from C21 and 28.17 million raw reads were generated from T43 (Sequence Read Archive accession number: SRP008320; Table 1). The raw reads had an average length of 90 bp. Among all the raw reads, more than 92% had Phred-like quality scores at the Q20 level (an error probability of 1%). These were used for de novo assembly. As a result, 459601 (C21) and 476229 (T43) contigs, and 97752 (C21) and 101931 (T43) scaffolds, were generated. After removal of redundancy, 68076 (C21) and 71194 (T43) unigenes were clustered using the Gene Indices Clustering Tools (GICT) program (Table 1, deposited in the Transcriptome Shotgun Assembly Sequence database. Accession numbers are summarized in Additional files 1 and 2). The random distribution of reads in the above-assembled unigenes indicated that the 3’ ends of all assembled unigenes contained relatively fewer reads compared with other positions within these unigenes, which showed greater and more even distributions (Figure 1).

6Using the CAP3 assembler, we obtained 65340 all-unigenes with an average length of 562 bp and an N50 of 770 bp by combining C21 (25.6 million) and T43 (27.3 million) clean reads (Table 1). Of the 65340 unigenes, 34745 were clustered from the assembled C21 (43506) and T43 (44558) unigenes. The remaining 30595 sequences were singletons with 14764 from C21 and 15831 from T43 unigenes. The length distributions of the unigenes from both C21 and T43 were similar with sequences in the 200–300 bp range making up 21.55% and 23.2% of C21 and T43 unigenes, respectively. The length of all-unigenes was dramatically increased following further assembly of the C21 and T43 unigenes, resulting in assembled sequences with lengths of > 300 bp (Figure 2). Ultimately, 54331 gap-free sequences were obtained by assembling all of the C21- and T43-unigenes.

7All assembled high-quality unigenes were first blasted against the NCBI non-redundant (nr) database using BLASTX with a cut-off E-value of 10-5. Of the 65340 all-unigenes, 35271 (53.9%) returned at least one match at the E-value > 10-5. Because of the lack of genome and EST information for R. trigyna or closely related taxa, 46.2% of the unigenes did not match to known genes in the database. Based on sequence homology, 13552 unigenes were categorized into 44 Gene Ontology (GO) terms (Figure 3, Additional file 3). In each of the three main GO classifications, i.e., biological process, cellular component, and molecular function, “cellular process”, “metabolic process”, “cell”, “cell part” and “binding” terms were dominant among the returned terms. Many of the identified unigenes were classified in the categories “response to stimulus”, “organelle” and “catalytic activity”, whereas only a few genes belonged to “biological adhesion”, “cell killing”, “locomotion”, “viral reproduction”, “virion” and “virion part”. Interestingly, 894 unigenes were classified into the category “transporter activity” and 71 unigenes into the category “antioxidant activity”.

8To further evaluate the integrity of our transcriptome library and the effectiveness of our annotation process, unigene sequences were subjected to Clusters of Orthologous Groups (COG) classification. Out of 35271 nr hits, 10968 sequences showed a COG classification (Figure 4, Additional file 3). Among the 25 COG categories, the cluster for “general function prediction only” (3179, 28.98%) was the largest group, followed by “replication, recombination and repair” (1681, 15.32%), “transcription” (1621, 14.78%), “signal transduction mechanisms” (1257, 11.46%), and “carbohydrate transport and metabolism” (1089, 9.93%). The categories “nuclear structure” (8, 0.07%) and “extracellular structures” (7, 0.06%) had the fewest corresponding genes.

9The 35271 annotated sequences were mapped to the reference canonical pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG). Among those, 15995 unigenes were assigned to 119 KEGG pathways (Additional file 3). The pathways most strongly represented by mapped unigenes were “metabolic pathways” (3459 unigenes), “biosynthesis of secondary metabolites” (1920 unigenes), and “plant-pathogen interaction” (1110 unigenes), indicating that these were the active pathways in R. trigyna. Finally, 38197 unigenes were annotated in the transcriptome of R. trigyna, of which 35459 were generated by BLAST analysis and 2738 were annotated using the ESTscan program.

10The analysis showed that 32697 unigenes were up-regulated and 31997 were down-regulated upon treatment of R. trigyna seedlings with salt solution. In total, 5032 DEGs were identified in both libraries, including 2370 up- and 2662 down-regulated unigenes (Figure 5). Among those DEGs responding to salt-stress treatment, 3719 were annotated and 1313 had no BLAST hits.

11Enrichment analysis was conducted to clarify the biological functions of the identified DEGs. The results indicated that 1947 DEGs were enriched in 27 GO terms, using a corrected P-value ≤ 0.05 as the threshold (Additional file 4). Among these GO categories, “oxidoreductase activity” (171 DEGs), “catalytic activity” (742 DEGs) and “response to stress” (155 DEGs) were significantly enriched among DEGs compared with the whole transcriptome background. In total, 2086 DEGs were enriched in 33 metabolic pathways (q-value ≤ 0.05) (Additional file 4). The most enriched metabolic pathways were “metabolic pathways” (483 DEGs), “biosynthesis of secondary metabolites” (351 DEGs), and “plant-pathogen interaction” (198 DEGs). The pathways of phenylpropanoid biosynthesis, flavonoid biosynthesis, stilbenoid, diarylheptanoid, and gingerol biosynthesis were very active, while those associated with ascorbate, aldarate, and nitrogen metabolism were strongly suppressed in salt-treated plants (Figure 6, Additional file 5).

12Among the unigenes, we identified 135, 72, and 25 genes related to regulation of K+ uptake, H+ pumping, and Na+ efflux, respectively (Figure 7, Table 2, Additional file 6). Clearly, genes associated with K+-transport formed the largest proportions of screened genes, suggesting that these genes played important role in reestablishing Na+/K+ homeostasis in salt-stressed R. trigyna. Among 135 K+ transport genes, 42 were DEGs, including 24 up- and 18 down-regulated genes in response to salt-stress. KUP (unigene14170, unigene26051) and CNGC (unigene3965, unigene7183) transcripts were present at high levels in both libraries, but were down-regulated by nearly 2-fold under salt-stress. Seven HKT1 genes, which are responsible for Na+ influx, were salt-responsive, of which three (unigene4890, unigene27586, and unigene20631) were suppressed by more than 2-fold. We identified 16, 44, and 12 genes encoding plasma membrane H+-ATPases (PM-H+-ATPases), vacuolar H+-ATPases (V-H+-ATPases) and H+-pyrophosphatases (V-H+-PPases), respectively. Out of 44 V-ATPase homologs, four showed down-regulated transcription under salt treatment while the other 40 did not. Transcripts of six PM-H+-ATPase homologs were highly abundant in both libraries. Two PM-H+-ATPases (unigene15650 and unigene248) were up-regulated by salt, showing increases in their Reads Per Kilobase of exon model per Million mapped reads (RPKM) values from 143.1 to 376.9 and from 79.8 to 240.3, respectively. Five PPases showed high transcript levels in both libraries, but their transcriptions did not change under salt-stress. Among 35 Na+ efflux genes, one PM-Na+/H+ antiporter gene (unigene798 or SOS1B) showed moderate transcript levels and two showed low transcript abundance. Transcription of unigene798 was slightly up-regulated by salt-stress (0.3-fold increase in abundance). Among 12 genes encoding V-Na+/H+ antiporter family proteins, four (unigene20634, unigene5272, unigene8445, and unigene752) were highly abundant. The transcript levels of unigene20634 and unigene5272 were up-regulated by nearly 2-fold, whereas those of unigene8445 and unigene752 were unchanged under salt-stress. There were abundant Nha (multisubunit Na+/H+ antiporter, unigene16859) transcripts in both libraries (RPKM values of 223.0 in C21 and 311.2 in T43).

13In total, 298 unigenes were predicted to encode enzymes related to reactive oxygen species (ROS) scavenging (Table 3, Additional file 7). There were 130, 84, 63, and 7 unigenes categorized into the GSH-ascorbate cycle, the peroxiredoxin/thioredoxin (PrxR/Trx) pathway, the glutathione peroxidase (GPX) pathway, and the catalase (CAT) pathway, respectively. The largest group was genes encoding thioredoxins (Trxs), glutathione S-transferases (GSTs), glutaredoxins (GLRs), and PODs (Figure 7). We detected 14 genes encoding SODs; three were up-regulated and 11 were down-regulated by salt treatment. In the GSH-ascorbate ROS removal pathway, seven ascorbate peroxidase (APX) genes (unigene44406, unigene47495, unigene17786, unigene13095, unigene57465, unigene25879, unigene10954), two GLR genes (unigene33573, unigene7889), one monodehydroascorbate reductase gene (unigene33304), one dehydroascorbate reductase gene (unigene42162), and two glutathione reductase genes (unigene8218, unigene3408) showed relatively abundant transcript levels (RPKM values of > 200). We identified five PrxRs, which play a role in the PrxR/Trx pathway in the antioxidant defense system. Two PrxR genes (unigene63275 and unigene61053) showed high transcript levels in both control and salt-treated plants (RPKM values of > 200). Out of 79 Trxs, five transcripts (unigene12493, unigene17142, unigene7940, unigene34328, and unigene2658) were expressed at high levels in both libraries. In the GPX pathway, there were 7 unigenes encoding GPXs but 53 unigenes encoding GSTs. Ten GST showed increased transcription under salt-stress, with RPKM values ranging from 102.6 to 477.2 in the T43 library but from 62.6 to 313.4 in the C21 library. In addition, 11 genes belonging to the Tau family of GSTs showed increased transcription under salt stress.

14PCR amplification showed that all qPCR primers produced only single fragments of the expected lengths (155–350 bp, a 100% success rate). The qPCR results for 30 selected unigenes showed general agreement with their transcript abundance changes determined by RNA-seq, suggesting the reliability of the transcriptomic profiling data. However, moderate discrepancies between the expression levels and RPKM values were observed in ten unigenes, i.e., unigene7220, unigene20645, unigene5040, unigene65175, unigene55260, unigene65219, unigene42056, unigene19588, unigene41893 and unigene39 (Figure 8). On the other hand, amplifications of longer sequences (> 500 bp) had a lower success rate. Among the 30 unigenes, 26 unigenes yielded fragments > 500 bp and the other four failed to produce amplification products, yielding a success rate of 86.7%. All positive clones from the validation studies were sequenced by Sanger sequencing, and the results completely confirmed the Illumina results.

15For many non-model species, there is no background genomic information available for researchers to conduct comprehensive investigations into the genetic mechanisms underlying their unique features. Therefore, the newly developed NGS technique has been used widely to explore genomic solutions to important physiological questions. Gene annotation is an important element of NGS in which biological information is attached to the predicted genes or unigenes. A high proportion of unigenes with high-confidence BLASTX similarity to protein sequences from annotated gene catalogs of other plant species is considered to indicate the integrity of transcript sequences assembled from Illumina short-read data. This assumption was also verified by the results generated in the present study. In R. trigyna, 54.27% of 65340 unigenes were annotated by BLAST analysis and functional bioinformatics analyses (e.g., GO, Swiss-Prot, and KEGG). Overall, the top five species with BLAST hits to annotated unigenes were Arabidopsis thaliana, Oryza sativa, Arabidopsis lyrata, Populus trichocarpa, and Vitis vinifera, species for which the annotations of their genomes are comprehensive and largely accepted. This suggested that the sequences of the R. trigyna unigenes generated in the present study were assembled and annotated correctly.

16The lengths of the assembled sequences are crucial in determining the level of significance of a BLAST match. Out of all the assembled unigenes, 66.79% had lengths between 300 and 500 bp, among which 37.1% had significant BLAST hits in the public databases. There were 1224 unigenes showing strong homology with the sequences hit in the database (E-value < 1.0E-50). Among 21697 unigenes with sequence lengths > 500 bp, 87.0% had significant BLAST scores and 59.2% showed E-values less than 1.0E-50. Out of all the unannotated unigenes (28991), 70.65% were longer than 500 bp, demonstrating that the lack of annotation for these unigenes was not because of a shorter sequence length but because of a genuine lack of hits to sequences in the database. Therefore, we can speculate that these unannotated unigenes represent a specific genetic resource of R. trigyna, which warrants further investigation.

17In this study, we also used other strategies to enhance the effectiveness of the short reads assembly, apart from the use of bioinformatics methods. These included preparation of high-quality RNA samples, removal of dirty raw reads, BLAST assembly of unigenes using multiple databases, and the large sample population for sequencing. First, RNA was isolated from sterile R. trigyna seedlings to minimize the risk of contamination by foreign RNase. Second, three seedlings were used to extract RNA samples, not only to reduce sample bias, but also to ensure comprehensive coverage of the R. trigyna transcriptome. Last, two sequencing libraries (C21 and T43) were merged to generate longer sequences and to increase the sequencing depth. These strategies were an effective way to improve the quality of assembly and annotation of assembled unigenes.

18In summary, an extensive and diverse expressed gene catalog, representing a large proportion of the transcribed genes in R. trigyna, was successfully sampled in the present study. The gene catalog provided a comprehensive understanding of the gene transcription profiles of R. trigyna, and laid a solid foundation for further study of salt-tolerance mechanisms and identification of novel genes in this species.

19Gene transcription and/or expression is often compared among different developmental stages, among different plant organs, or among plants under different growth conditions. In this study, many genes showing transcriptional changes under salt stress were identified by comparing NaCl-treated seedlings with untreated controls. There were 64694 unigenes showing differences in transcript abundance, and 5032 were defined as DEGs using the thresholds of false discovery rate (FDR) ≤ 0.001 and |log2Ratio| ≥ 1. GO clustering analysis suggested the potential biological functions of these DEGs. For example, many DEGs were enriched in GO terms such as “oxidoreductase activity”, “catalytic activity” and “response to stress”. This information will be useful to elucidate salt-tolerance mechanisms and to find new salt-stress-related genes specific to R. trigyna. KEGG enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways involving DEGs. Surprisingly, genes in the phenylpropanoid biosynthesis and flavonoid biosynthesis pathways were transcribed at much high levels under salt stress than in the control in R. trigyna, which was not expected from the results of previous studies. This result suggested that genes in these two pathways may play vital roles in resisting oxidation and maintaining membrane integrity.

20An outstanding advantage of NGS is the detection of transcripts present in low copy numbers. Among the 5032 identified DEGs, almost half were detected as low-abundance transcripts that were up-regulated after exposure to salt-stress. For example, only 14 reads from C21 could be mapped to unigene24841, encoding flavanone 3-hydroxylase (F3H), compared with 470 reads from T43. Nearly 30-fold up-regulation suggested that this gene may act as a key element in the salt-stress response. Overall, transcriptome comparison analysis provided sufficient information to study the salt-responsive mechanisms and genes related to salt-tolerance related in R. trigyna.

21Ion transport is a crucial element in response to salt stress in plants. This is particularly true in halophytes. It has been shown that halophytes are able to grow under extreme salinity conditions because of their anatomical and morphological adaptations and/or their avoidance mechanisms. R. trigyna is a representative recretohalophyte with unique morphological characteristics, such as salt glands and succulent leaves, allowing it to adapt to the salinized conditions of the Ordos desert. Our previous study showed that the salt glands of this species functionally excrete salt ions under normal and salt-stressed conditions. However, the amount of Na+ and Cl- excreted significantly increased under salt treatment. In addition, R. trigyna shows a strong ability to uptake K+ from barren soil. Therefore, it is possible that the ion transport mechanisms and their gene expression patterns in this plant differ from those of other plant species to some extent.

22Once Na+ is taken up into the cell, ATPases (PM-H+-ATPases, V-H+-ATPases, and V-H+-PPases) are induced to create the driving force for Na+ transport. This results in not only extrusion of Na+ into the external environment by Na+/H+ antiporters in the plasma membrane, but also in compartmentalization into vacuoles by tonoplast Na+/H+ antiporters, which are essential for reestablishing cellular ion homeostasis in salt-stressed plants. In this study, five V-H+-PPases were detected in R. trigyna, and were proposed to generate a proton electrochemical gradient. Taking the morphological characteristics of this plant into account, the largest share of the driving force generated by such enzymes is likely to be consumed by the succulent leaves and the salt glands, which include two mesophyll-like collecting cells with obvious vacuoles, and six secretory cells. All mesophyll cells in leaves need V-H+-ATPases and V-H+-PPases to generate an electrochemical gradient to compartmentalize excess Na+ ions into the tonoplast, as a mechanism to cope with physiological drought conditions and ion toxicity. When Na+ is excreted via secretory cells, both V-H+-ATPases and V-H+-PPases are required to produce an H+ electrochemical gradient to pump Na+ into collecting cells. PM-H+-ATPases generate sufficient driving force to exchange Na+ ions out of the cells. In addition, higher H+-pumping activity may be indispensable to maintain the high concentration of K+ in salt-stressed R. trigyna. The activation of PM-H+-ATPases in salinized Populus euphratica led to repolarization or hyperpolarization of the plasma membrane and thus decreased NaCl-induced K+ loss through depolarization-activated outward rectifying K+ channels (DA-KOCs). In R. trigyna, 10 out of 17 RtKOCs showed high homology with those of Populus, demonstrating similar functions of these genes. The analysis of transcripts related to K+ uptake showed that the outstanding potassium uptake capacity of R. trigyna was probably conferred by enhanced KUP/HAK (potassium transporter) and AKT K+ uptake systems, among which the number of up-regulated unigenes (39) far exceeded that of down-regulated unigenes (16) under salt stress. In summary, proton electrochemical gradient-dependent K+ maintenance, ions secretion, compartmentalization, and enhanced K+ uptake systems may represent important components in reestablishing ion homeostasis in R. trigyna.

23On the other hand, we detected 12 NHXs that were probably responsible for Na+ sequestration. The significantly up-regulated unigene20634 and 5272 showed high homology to the Citrus reticulata NHX1 and a Salsola komarovii Na+/H+ antiporter, respectively. The slightly down-regulated unigene8445 and 752 had a highly significant BLAST hit to the Na+/H+ antiporter of M. crystallinum and the NHX1 of Tetragonia tetragonioides, respectively. These four unigenes also showed strong homology to the NHX2 of A. thaliana, suggesting that these AtNHX2-like proteins play an important role in avoiding or mitigating the deleterious effects of high Na+ levels in the cytosol or in regulating intravacuolar K+ and pH, which has been demonstrated in Arabidopsis. Our prediction is also supported to some degree by the studies of succulent leaves of the halophytes Suaeda salsa and Halostachys caspica.

24In R. trigyna, we identified only one moderately expressed and up-regulated PM Na+/H+ antiporter (SOS1B, unigene798), named RtSOS1. The transcript profile of this gene was consistent with those of ThSOS1 and PeSOS1, which remained relatively constant or were slightly up-regulated under salt-tress. Escherichia coli complementation experiments with PeNhaD1 rescued salt-exposed bacteria by lowering salt accumulation, suggesting that the P. euphratica gene PeNhaD1 also functions as an Na+⁄H+ antiporter. Interestingly, a highly abundant transcript (unigene16859, RtNha1) encoding an Nha protein was identified in our dataset. This may play a role in Na+ exclusion in R. trigyna. At present, it is still unknown how the moderately up-regulated RtSOS1 achieves effective efflux of excess salt ions under extreme salt-stress, and whether any other gene products are also involved in this process. This requires further investigations in the near future.

25Salt stress causes a rapid increase in ROS, including superoxide radicals (·O2-), hydrogen peroxide (H2O2), and hydroxyl radicals (·OH), which perturb cellular redox homeostasis and result in oxidative damage to many cellular components and structures. According to our previous studies, ·O2- increased rapidly in R. trigyna seedlings treated with 400 mM NaCl, as did the concentrations of other antioxidants, for example, GSH. At the same time, increased activities of antioxidant enzymes (SOD, POD) were also detected. This suggested that R. trigyna probably has a similar ROS scavenging system to that of other plant species, and that this system can be enhanced to increase its antioxidative ability. As the first line of defense against oxidative damage, SODs are usually induced by salinity to rapidly dismutate ·O2- into oxygen and H2O2, which is subsequently removed through various pathways. In this study, 14 SODs were identified. Two of them (unigene23472, unigene8340) were transcribed at high levels, and showed up-regulated transcription under salt stress. Therefore, they were probably associated with enhanced SOD activity and responsible for converting the increased ·O2- into H2O2. On the other hand, increased GSH could be consumed in the ascorbate-GSH cycle and the GPX cycle. In the ascorbate-GSH cycle, the oxidized ascorbate (i.e., monodehydroascorbate) can be converted into dehydroascorbate (DHA); DHA is then reduced to ascorbate at the expense of GSH, yielding oxidized GSH (i.e., glutathione disulfide). In the GPX pathway, GPX can reduce H2O2 to the corresponding hydroxyl compounds using GSH. However, to fulfill its function as an antioxidant, GSH must be catalyzed by GSTs, which have GPX activity and can use GSH to reduce organic hydroperoxides of fatty acids and nucleic acids to the corresponding monohydroxy alcohols. Interestingly, in our dataset, there were 35 salt-induced GSTs genes that likely participated in such catalytic processes, especially the 11 Tau family GSTs. This finding was similar to the results of other studies on salt-stressed plants. Therefore, the interaction between GSH and GST may be an important factor in the ROS scavenging system of R. trigyna. In addition, we identified 84 genes encoding components of the PrxR/Trx defense system, which employs a thiol-based catalytic mechanism to reduce H2O2 and is regenerated using Trxs as electron donors.

26This study describes a platform in which publicly available information was matched to the R. trigyna transcriptome using NGS technology. The substantially assembled sequences represented a considerable portion of the transcriptome of this plant. Transcriptome comparison identified the transcripts and metabolic pathways that play significant roles in the response to salt stress. These results will prompt studies on the molecular mechanisms of salt resistance, facilitate the analysis of expression profiles of genes related to salt tolerance, and accelerate the discovery of specific stress-resistance genes in R. trigyna. The information provided here can also further extend the knowledge on the salt tolerance of halophytes that survive in high sodic stress in areas such as the semi-desert saline areas in Ordos, Inner Mongolia, China.

27Seeds of R. trigyna were collected in September 2009 from the Eastern Alxa-Western Ordos area in Inner Mongolia, China. Intact plump seeds were selected, immersed in 10% sodium hypochlorite for 15 min, and then rinsed three times with sterilized double-distilled water. The seeds were germinated in a 150-ml conical flask containing 40 ml MS medium in the dark for 72 h. Germinated seeds were then grown on the same medium at 25°C under 70% relative humidity and a 16-h light/8-h dark photoperiod for 15 d. When seedlings were approximately 10 cm high, they were transferred into a test tube containing 50 ml ½-strength Hoagland’s medium, and further cultured for another 4 weeks with a change of medium every 2 days. Three healthy seedlings of a similar size were selected for further analysis. Prior to NaCl treatment, leaf samples were collected for RNA extraction from each seedling as controls. The roots were then immersed in Hoagland’s medium containing NaCl. The NaCl concentration was increased stepwise from 100 mM to 400 mM at a rate of 100 mM per step every 8 h. Leaves were collected before each increase in NaCl concentration. Leaves were also collected after the seedlings were treated in 100 mM NaCl for 0.5 h, 1 h, 2 h, 4 h, and 8 h. The leaf samples were immediately snap-frozen in liquid nitrogen and stored at −80°C until analysis.

28Total RNA was extracted with Plant Plus RNA reagent (DP437, Tiangen, Beijing) according to the manufacturer’s instructions. The extracted RNA was treated with RNase-free DNase I (TaKaRa Bio Inc., Otsu, Shiga, Japan) for 45 min at 37°C to remove residual DNA. RNA quality was evaluated using the Agilent 2100 Bioanalyzer with a minimum RNA integrity number value of 8. In this study, two sequencing libraries were prepared, the C21 library from control samples and the T43 library from salt-stressed samples.

29Beads with Oligo (dT) (Illumina) were used to isolate poly (A) + RNA from 20 μg of each RNA pool. Fragmentation buffer (Ambion, Austin, TX, USA) was added to break mRNA into short fragments. First-strand cDNA was synthesized using N6 random hexamers (Illumina Inc., San Diego, CA, USA) and SuperScript® III reverse transcriptase (Invitrogen, Grand Island, NY, USA). The second-strand cDNA was synthesized with DNA polymerase I (Invitrogen, Grand Island, NY, USA). Short fragments were purified with a QIAquick PCR extraction kit (Qiagen, Hilden, Germany) and resolved with EB buffer for end-repair and poly (A) addition, and then linked to sequencing adapters. Agarose gel electrophoresis and PCR amplification were used to select suitable fragments. Two cDNA libraries were constructed and sequenced on the Illumina HiSeq™ 2000 platform. The sequences were base-called and quality-checked by the Illumina data processing pipeline.

30The raw reads were filtered to obtain high-quality clean reads prior to assembly. This was performed by removing adaptor sequences, duplicated sequences, reads containing more than 5% “N” (i.e., ambiguous bases in reads), and reads in which more than 50% of the bases showed a Q-value ≤ 5. Transcriptome de novo assembly was carried out with the short reads assembly program, SOAPdenovo, with the default settings, except that the K-mer value was 25-mer. The different K-mer sizes were assessed; the 29-mer yielded the best assembly for the desired application and so it was used to construct the de Bruijn graph. Contigs without ambiguous base reads were obtained by conjoining K-mers in an unambiguous path. The reads were then mapped to contigs to construct scaffolds with paired-end information. Paired-end reads were used again for gap filling of scaffolds to obtain unigenes. To evaluate the depth of sequence coverage, all usable reads were realigned to the unigenes with SOAPaligner. To minimize sequence redundancy, the scaffolds were clustered using the GICT. The clustering output was passed to the sequence clustering software CAP3 assembler for multiple alignment and generation of consensus sequences. The sequences that did not reach the threshold set and did not fall into any assembly were assigned as singletons.

31To assign predicted gene descriptions for the assembled unigenes, they were aligned against the plant protein dataset of nr, Swiss-Prot/Uniprot protein database, and COG databases, respectively, using BLASTX with a significance threshold of E-value ≤ 10-5. To identify the best BLASTX hits from the alignments, putative gene names, ‘CDS’, and predicted proteins of the corresponding assembled sequences were produced. The orientation of sequences was also derived from BLAST annotations. For sequences aside from those obtained from BLAST searches, the EST Scan program was used to predict the ‘CDS’ and its orientation. Functional categorization by GO terms was performed based on the best BLASTX hits from the nr database using BLAST2GO software according to molecular function, biological process, and cellular component ontologies with an E-value threshold of 10-5. The data were statistically analyzed using WEGO software. The pathway assignments were carried out by sequence searches against the KEGG database, also using the BLASTX algorithm with an E-value threshold of 10-5.

32To identify DEGs between the two samples, the transcript abundance of all assembled unigenes were analyzed using the RPKM method. With reference to “the significance of digital gene expression profiles”, a rigorous algorithm has been developed to identify DEGs between two samples. The probability of gene A being transcribed equally between two samples was calculated using the following formula:

33where the P-value corresponds to the differential gene expression test; N1 and N2 represent the total clean tag number of C21 and T43 samples, respectively; and x and y denote the tag numbers of the gene of interest present in C21 and T43, respectively.

34The FDR is a method to determine the threshold P-value in multiple tests and analysis through manipulating the FDR value. FDR ≤ 0.001 and the absolute value of log2Ratio ≥ 1 were used as the threshold to identify DEG. Functional enrichment analysis including GO and KEGG were performed using the following ultra-geometric test to find which DEGs were significantly enriched in GO terms (P-value ≤ 0.05) and metabolic pathways (q-value ≤ 0.05) compared with the whole transcriptome background.

35To experimentally validate the transcriptional abundance results from sequencing and computational analysis, 30 unigenes (10 up-regulated, 10 down-regulated, and 10 unigenes with no significant changes in transcript abundance) were selected for RT-PCR and qPCR analysis. Reverse transcription reactions were performed using SuperScript® III Reverse Transcriptase (Invitrogen, Grand Island, NY, USA) with approximately 5 μg total RNA following the manufacturer’s instructions. Primers for RT-PCR and qPCR were designed using Primer Premier 5 and Beacon Designer 7.0 software, respectively (shown in Additional file 8). β-Actin was used as the internal control gene. qPCR was performed on a Qiagen Rotor-gene Q realtime PCR platform (Qiagen, Hilden, Germany) using SYBR-Green real-time PCR mix (TaKaRa Bio Inc., Otsu, Shiga, Japan) to detect transcript abundance. The amplification was achieved by the following PCR protocol: first denaturation 95°C for 30 s, then 40 cycles of denaturation at 95°C for 5 s, annealing and extension at 55°C for 30 s. The relative expression levels of the selected unigenes normalized to β-Actin was calculated using 2-∆∆Ct method. All reactions were performed with three replicates, and data were analyzed using Rotor-gene Q series software.

36NGS: Next-Generation Sequencing; DEG: Differentially Expressed Genes; RPKM: Reads Per Kilobase of exon model per Million mapped reads; FDR: False Discovery Rate; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; COG: Clusters of Orthologous Groups; AKT: Potassium Channel; HKT: High-Affinity K+ Transporter; CNGC: Cyclic Nucleotide-Gated Channel; KUP: K+ Ion Transmembrane Transporter; HAK: High Affinity K+ Transporter; CHX: Cation H+ Exchanger; KEA: K+ Efflux Antiporter; KOC: Outward Rectifying K+ Channel; SOS1: Salt Overly Sensitive 1; NHX: Na+/H+ Exchanger; NhaP: Nhap-Type Na+/H+ Antiporter; V-ATPase: V-H+-ATPase, vacuolar H+-ATPase; V-PPase: V-H+-PPase, H+-Pyrophosphatase; P-ATPase: PM-H+-ATPase, Plasma Membrane H+-ATPase; ROS: Reactive Oxygen Species; GSH: Glutathione; ·O2-: Superoxide Radicals, H2O2, Hydrogen Peroxide; ·OH: Hydroxyl Radicals; GLR: Glutaredoxin; APX: Ascorbate Peroxidase; MDAR: Monodehydroascorbate Reductase; DHAR: Dehydroascorbate Reductase; GR: Glutathione Reductase; GST: Glutathione S-Transferase; GPX: Glutathione peroxidase; POD: Peroxidases; GLP: Germin-Like Protein; CAT: Catalase; PEX11: Peroxisomal Biogenesis Factor11; Trx: Thioredoxin; PrxR: peroxiredoxin; SOD: Superoxide Dismutase; A. lyrata: Arabidopsis lyrata subsp. lyrata; A. thaliana: Arabidopsis thaliana; A. hypogaea: Arachis hypogaea; B. gymnorhiza: Bruguiera gymnorhiza; C. sinensis: Camellia sinensis; C. chinense: Capsicum chinense; C. quinoa: Chenopodium quinoa; C. limon: Citrus limon; C. unshiu: Citrus unshiu; C. reticulata: Citrus reticulata; C. sativus: Cucumis sativus; D. carota: Daucus carota; E. guineensis: Elaeis guineensis; E. camaldulensis: Eucalyptus camaldulensis; G. max: Glycine max; G. hirsutum: Gossypium hirsutum; H. caspica: Halostachys caspica; H. brasiliensis: Hevea brasiliensis; H. vulgare: Hordeum vulgare; I. batatas: Ipomoea batatas; I. lacteal: Iris lacteal; J. curcas: Jatropha curcas; K. foliatum: Kalidium foliatum; L. bicolor: Limonium bicolor; L. chinensis: Litchi chinensis; M. pumila: Malus pumila; M. pusilla: Malva pusilla; M. truncatula: Medicago truncatula; M. crystallinum: Mesembryanthemum crystallinum; N. nucifera: Nelumbo nucifera; N. alata: Nicotiana alata; N. benthamiana: Nicotiana benthamiana; N. tabacum: Nicotiana tabacum; O. sativa: Oryza sativa; P. acutifolius: Phaseolus acutifolius; P. sativum: Pisum sativum; P. euphratica: Populus euphratica; P. trichocarpa: Populus trichocarpa; P. juliflora: Prosopis juliflora; P. persica: Prunus persica; P. tenuiflora: Puccinellia tenuiflora; P. communis: Pyrus communis; R. australe: Rheum australe; R. communis: Ricinus communis; S. komarovii: Salsola komarovii; S. lycopersicum: Solanum lycopersicum; S. nigrum: Solanum nigrum; S. tuberosum: Solanum tuberosum; S. oleracea: Spinacia oleracea; S. salsa: Suaeda salsa; T. androssowii: Tamarix androssowii; T. hispida: Tamarix hispida; T. tetragonioides: Tetragonia tetragonioides; T. octandrum: Trichostigma octandrum; V. radiate: Vigna radiate; V. unguiculata: Vigna unguiculata; V. vinifera: Vitis vinifera; Z. aethiopica: Zantedeschia aethiopica; Z. mays: Zea mays; Z. jujuba: Ziziphus jujuba; Z. marina: Zostera marina.

37The authors declare that they have no competing interests.

38ZHD cultivated plants; performed RNA extraction, RT-PCR, and qPCR; analyzed, interpreted, and submitted data; and drafted the manuscript. LLZ and ZQ participated in the design of the study and data analysis. JW and ZG participated in plant management and data submission to the database. SBW participated in analysis and interpretation of data, and critically revised the manuscript. YCW conceived of the study, participated in its design and data interpretation, and revised the manuscript critically. All authors read and approved the final manuscript.