Transcriptomic profiles of Dunaliella salina in response to hypersaline stress

Background Dunaliella salina is a good model organism for studying salt stress. In order to have a global understanding of the expression profiles of Dunaliella salina in response to hypersaline stress, we performed quantitative transcriptomic analysis of Dunaliella salina under hypersaline stress (2.5 M NaCl) of different time duration by the second and third generation sequencing method. Results Functional enrichment of the up-regulated genes was used to analyze the expression profiles. The enrichment of photosynthesis was observed, accompanied by enrichments of carbon fixation, pigment biosynthetic process and heme biosynthetic process, which also imply the enhancement of photosynthesis. Genes responsible for starch hydrolysis and glycerol synthesis were significantly up-regulated. The enrichment of biosynthesis of unsaturated fatty acids implies the plasma membrane undergoes changes in desaturation pattern. The enrichment of endocytosis implies the degradation of plasma membrane and might help the synthesis of new glycerophospholipid with unsaturated fatty acids. Co-enrichments of protein synthesis and degradation imply a higher protein turnover rate. The enrichments of spliceosome and protein processing in endoplasmic reticulum imply the enhancement of regulations at post-transcriptional and post-translational level. No up-regulation of any Na+ or Cl− channels or transporters was detected, which implies that the extra exclusion of the ions by membrane transporters is possibly not needed. Voltage gated Na+ and Cl− channels, mechanosensitive ion channel are possible signal receptors of salt stress, and Ca2+ and MAP kinase pathways might play a role in signal transduction. Conclusion At global transcriptomic level, the response of Dunaliella salina to hypersaline stress is a systematic work, possibly involving enhancements of photosynthesis, carbon fixation, and heme biosynthetic process, acceleration of protein turnover, spliceosome, protein processing in endoplasmic reticulum, and endocytosis, as well as degradation of starch, synthesis of glycerol, membrane lipid desaturation. Altogether, the changes of these biological processes occurred at trancriptomic level will help understand how a new intracellular balance achieved in Dunaliella salina to adapt to hypersaline environment, which are worth being confirmed at the physiological levels.

glycerol is converted to dihydroxyacetone by glycerol dehydrogenase, and then dihydroxyacetone is converted to DHAP by dihydroxyacetone kinase [3]. As the key enzyme in the pathway, GPDH was extensively studied [4,5]. Third, the Na + /H + antiporter activity was detected in plasma membrane and was thought to function as exclusion of Na + in vivo [6,7].
In twenty-first century, proteomic methods were used to understand the molecular mechanism of salt tolerance at omics level. Proteins such as transferrin, carbonic anhydrases, Na + /H + antiporter, fatty acid elongase, GPDH, small GTP-binding protein and tubulin were found upregulated significantly under salt stress. These proteins can be classified in carbon assimilation, energy production, transporters, signal transduction, protein synthesis and cell defense [8,9]. However, due to the limitation of the two-dimensional electrophoresis, the information obtained from this technique is limited [8,9], with detected number of differently expressed proteins below 100, of which only about 60% can be annotated.
Compared with proteomic approaches, transcriptomic methods are more reproducible, sensitive with higher genome coverage. A transcrtiptome of 17,845 transcripts was reported when Dunaliella tertiolecta was investigated to identify the pathways and genes involved in lipid synthesis under nitrogen stress, which covers about 97% of the core eukaryotic genes (CEGs) [10][11][12]. Hong et al. reported the transcriptome of Dunaliella salina at different phases of their growth cycle (30d, 80d, 120d), but no transcriptome under salt stress was reported [13]. Alkayal reported the expressed sequence tag (EST) profiling of Dunaliella salina after 5 h of hypersaline shock, in which a transcriptome of 1401 unique transcripts was reported and the annotated transcripts can be classified into protein synthesis, energy, primary metabolism and protein fate [14]. However, no transcriptome before salt stress was generated, so there was no comparison to present the underlying changes during this shock period. In order to have a better understanding of how Dunaliella salina responds to hypersaline shock at transcriptomic level, the second and the third generation sequencing were used to generate the transcriptome of Dunaliella salina at different duration time under stress. Because intracellular glycerol synthesis accomplished in about 2 h after hypersaline shock [15,16], we reported the transcriptomic profiles on time duration of 0.5-h, 1h and 2-h under hypersaline stress and the profiles were compared with those before stress.

Data quality and sequences annotation
To obtain high quality sequence data, total RNAs of high quality were extracted (not shown). After the second generation sequencing, each library gave high quality clean reads with Q20% ranging from 97.21 to 98 with error rate about 0.01% (Additional file 1: Table S1). The GC content is about 56.5%, which is close to Dunaliella salina (CCAP19/18) reported [17] . The number of the named "full length transcripts" generated from the third generation sequencing was 43,864, ranging from 242 to 8978 bp in length with a mean length of 1009 bp and median length of 918 bp. About 80% transcripts of them are in the length range of 500 to 2000 bp (Additional file 1: Table S2). By ORF analysis, among the 43,864 transcripts, 35,175 transcripts are classified into coding sequences and 8689 transcripts are classified into long non-coding sequences. Among the 35,175 coding sequences, 29,071 sequences are annotated and 6104 sequences cannot be annotated so far. In order to estimate the coverage of the transcriptome, transcripts hit the same gene (the same sequence ID) in Nr, Nt or SwissProt database are defined as the splice variants generated by alternative splicing from a single gene. By this method isoenzymes and artificially spliced sequences are also excluded. Finally 9256 individual genes from the 29,071 transcripts are generated. Genome sequencing of Dunaliella salina (CCAP19/18) and Chlamydomonas reinhardtii predicted 16,697 and 17,741 loci containing protein-coding transcripts respectively [17,18]. Compared with the predicted numbers of gene loci of the two green algae, 9256 is a rather high number, since many genes aren't expressed and their mRNAs can't be detected. Furthermore, approximately 87.1% of the core eukaryotic genes (CEGs) were identified from the 9256 individual genes by sequence similarity search which suggests a rather high coverage of the Dunaliella salina transcriptome.

General pattern of gene expression
Based on gene expression value, clustering analysis was performed (Additional file 2: Figure S1), we can see the similarities of the expression patterns of the samples with good repeatability in the same group (the same stress time).
While compared with the 0-h of stress (no salt stress was applied), the number of differentially expressed genes increased with the increasing of stress duration time (Fig. 1). The number of up-regulated genes increases from 569 on 0.5-h of stress to 915 on 1-h of stress, and then to 3071 on 2-h of stress. On the other hand, the number of down-regulated genes increases from 513 on 0.5-h stress to 810 on 1-h stress, and then to 2580 on 2-h stress.
In order to have an overall understanding of the upregulated genes under salt stress, functional enrichments were performed by GO (gene ontology) ( Table 1). On 0.5h of stress, carboxylic acid biosynthetic process, cellular lipid metabolic process, carbohydrate metabolic process, response to temperature stimulus, photosynthesis (light harvesting), photosynthesis (light reaction), cofactor metabolic process, pigment biosynthetic process, and tetrapyrrole biosynthetic process are significantly enriched. On 1h of stress, protein folding and DNA replication are included in the list of significantly enriched biological processes, cellular lipid metabolic process and response to temperature stimulus are enriched but not statistically significant, while photosynthesis is excluded due to rapid decreasing of gene number (Table 2). On 2-h of stress, new terms such as macromolecule modification, cellular catabolic process, cell redox homeostasis, reproductive process, and ferrous iron transport are significantly enriched, while transcription (DNA-templated) is enriched, but not statistically significant. The terms enriched on 1-h of stress, such as carboxylic acid metabolic process, cellular lipid metabolic process, carbohydrate metabolic process, response to temperature stimulus, cofactor metabolic process, protein folding, and DNA replication, are also enriched and show a quick increasing of the gene numbers compared with that of 1-h of stress. These biological processes are not statistically significant due to the rapid increasing of the number of the total upregulated genes, but they are still worth focusing on. In general, the significantly enriched biological processes can be classified into photosynthesis, carbohydrate metabolism, lipid metabolism, and amino acids and protein metabolism. We focused on analyzing these biological processes in the following sections.
On the other hand, the functional enrichment of the down-regulated genes was also performed by GO (Additional file 1: Table S3). On 0.5-h of stress, no terms were significantly enriched, but carbohydrate binding and protein binding were worth focusing on since the numbers of down-regulated genes involved are large. On 1-h of stress, DNA metabolic process, protein binding, cytoskeleton, glycoprotein biosynthetic process, glycosaminoglycan biosynthetic process, and dynein complex were significantly enriched. On 2-h of stress, more GO terms were significantly enriched beside the terms enriched on 1-h of stress, these terms include transferase activity, protein modification process, regulation of RNA biosynthetic process, response to nitrate, inorganic anion transport, lipid transport, DNA integration, autophagy, and GTPase activator activity. From the point of gene numbers, we can see that the down-regulated genes are mainly involved in protein binding, transferase activity, protein modification process, DNA metabolic process, regulation of RNA biosynthetic process, and cytoskeleton. These terms are also important for understanding the hypersaline stress of Dunaliella salina, however, this paper only focuses on the analysis of the terms enriched by the up-regulated genes.

Photosynthesis
On the 0.5-h of stress, photosynthesis-light reaction and photosynthesis-light harvesting are significantly enriched by GO analysis on up-regulated genes, which implies the enhancement of photosynthesis. In time course, most genes are highly expressed on 0.5-h, decreased a little on 1-h, and then decreased to low levels even lower than that of 0-h. The expression pattern is like a pulse style and most peaks of gene expression are induced on or before 0.5-h of stress (Fig. 2). Many genes of Chlorophyll a-b binding proteins show pulse expression patterns, such as Chlorophyll a-b binding protein of LHCII type I, Chlorophyll a-b binding protein type 1 member F3, Chlorophyll a-b binding protein P4, and Chlorophyll a-b binding protein CP29 et al.. Some of the genes show high expression on 2-h of stress, including ATPdependent zinc metalloprotease FTSH 2, Photosystem II repair protein PSB27-H1, D-amino-acid transaminase, and Photosystem II protein D1. A few genes show a decreasing of expression, including Protein TIC 20-II, Oxygen-evolving enhancer protein, and DNA-binding  (Fig. 2). Chlorophyll biosynthetic process is also enriched, which indicates the synthesis of photosynthetic pigments to enhance photosynthesis ( Table 2). This is consistent with previous study [9]. With the stress going on, the gene numbers of photosynthesis-light reaction and photosynthesis-light harvesting decreased (Table 2), while the gene number of carbon fixation constantly increased, from 12 on 0.5-h to 25 on 1-h, and to 39 on 2-h of stress (Table 2, Additional file 2: Figure S2), key genes such as carbonic anhydrase and rubisco activase are significantly upregulated (Additional file 1: Table S4). Compared with the decreased gene number of photosynthesis-light reaction and photosynthesis-light harvesting, the constantly increased gene number of carbon fixation indicates that these biological processes may be controlled by different signaling pathways.
With the stress going on, the gene number of chlorophyll biosynthetic process decreased, while the gene number of tetrapyrrole biosynthetic process remained stable and the gene number of heme biosynthetic process kept increasing ( Table 2). The increasing of gene number of heme biosynthetic process and the decreasing of gene number of chlorophyll biosynthetic process together resulted in the stableness of gene number of tetrapyrrole biosynthetic process since the latter is the father term of the former two. This is consistent with the result of heat-map analysis, of which some genes show pulse expression pattern, these genes are clustered to chlorophyll biosynthetic process, while some genes show high expression values on 2-h of stress, these genes are clustered to heme biosynthetic process (Additional file 2: Figure S3). The significant enrichment of tetrapyrrole biosynthetic process and heme biosynthetic process on 0.5-h and 1-h of stress are very interesting. In plants and algae, tetrapyrroles are plastid signals demonstrated to regulate nuclear gene expression [19][20][21][22]. Heme signaling also appears to play a role in starch biosynthesis and drought tolerance in plants [23,24]. We see the constant increasing of gene number of heme biosynthetic process with the increasing of stress time, while large amount of signal molecules are usually not needed, so the constant increasing gene number of heme synthesis could be for the synthesis of heme-containing enzymes, such as catalase and ascorbate peroxidase, which play important roles in detoxification of reactive oxygen species (ROS) [25]. Consistently, the expression of ascorbate peroxidase is up-regulated and also confirmed by qPCR (Additional file 1: Table S4, Additional file 3).

Glycolysis and glycerol synthesis
Glycolysis is significantly enriched by KEGG Pathway analysis on up-regulated genes. The up-regulations of PGM, GPI, the rate-limiting enzyme PFK1 (6-phosphofructokinase 1, 2.7.1.11), and fructose-bisphosphate aldolase were seen on 0.5-h of stress (Additional file 1: Table S4), which implies alphpa-D-Glucose-1p from hydrolysis of starch goes to glycolysis (Fig. 4). Interestingly,  Our data show that the Dunaliella salina specific didomain glycerol-3-phosphate dehydrogenase (DsGPDH) can convert DHAP (an intermediate of glycolysis) to glycerol directly [27]. So the glycerol synthesis pathway of Dunaliella salina can be drawn based on the genes from   (Fig. 4). The new glycerol synthesis pathway includes only 6 enzyme-catalyzing steps to convert starch to glycerol, the six enzymes are all up-regulated on 0.5-h of stress. The up-regulation of TPI can be explained as that more enzyme is needed to converting GADP to DHAP for glycerol synthesis. Consistently, DsGPDH is up-regulated 0.5, 0.9 and 1.0 fold on 0.5-h, 1-h, and 2-h of stress respectively (Additional file 1: Table S4).

Lipid metabolism
On 0.5-h of stress, fatty acid biosynthesis is significantly enriched by KEGG Pathway analysis on up-regulated genes. The number of genes involved in this pathway remains stable on 1-h of stress, and decreases a little on 2h of stress. While the gene number of fatty acid degradation and biosynthesis of unsaturated fatty acids increase constantly ( On 2-h of stress, key genes involved in biosynthesis of unsaturated fatty acids are also significantly up-regulated, including delta-7 desaturase, omega-3 fatty acid desaturase and delta-12 desaturase (omega-6 fatty acid desaturase), acyl-[acyl-carrier-protein] desaturase and acyl-CoA oxidase, implying the synthesis of unsaturated fatty acids (Additional file 2: Figure S4). The genes classified into fatty acid degradation by KEGG include long-chain acyl-CoA synthetase, acyl-CoA oxidase, alcohol dehydrogenase class-3, acetyl-CoA c-acetyltransferase, and glyoxysomal fatty acid beta-oxidation multifunctional protein MFP-a. Among them, long-chain acyl-CoA synthetase and acyl-CoA oxidase are also classified into fatty acid biosynthesis, but the up-regulation of the latter three genes implies the degradation of fatty acid. On the other hand, the number of up-regulated genes involved in Glycerophospholipid and glycerolipid metabolism increase constantly with the increasing of stress time (Table 4). Genes catalyzing the formation of phospatidylcholine (lecithin), phosphatidylethanolamine, phosphatidyl-L-serine, phosphatidylglycerol, and phosphatidyl-D-myo-inositol are significantly up-regulated on 2-h of stress (Additional file 1: Table S4), indicating the synthesis of glycerophospholipid on 2-h of stress. Study reported that fatty acids of Dunaliella salina plasma  [28]. So it is possible that the new synthesized unsaturated fatty acids might be incorporated into the new glycerophospholipids.
Interestingly, endocytosis, which can absorb and start the degradation of cell membrane [29], is significantly enriched on 2-h of stress (Table 4). When confronting salt stress, Dunaliella cells begin to shrink to cause excessive cell membrane. It is possible that the excessive membrane is degraded into building blocks for synthesis of new glycerophospholipid with unsaturated fatty acids, or further degraded for fuel energy metabolism.

Amino acids and protein metabolism
Several amino acid metabolism pathways are enriched by KEGG Pathway analysis on up-regulated genes (Table 5), including some high abundant amino acids such as glutamine, glutamic acid, serine, alanine, proline, aspartic acid and asparagine. These high abundant amino acids were reported to be synthesized during abiotic stress to act as compatible osmolytes, precursors for secondary metabolites, or storage forms of organic nitrogen [30]. Several amino acid biosynthesis pathways are also enriched especially on 2-h of stress, including phenylalanine, tyrosine and tryptophan biosynthesis, lysine biosynthesis, valine, leucine and isoleucine biosynthesis, and histidine metabolism (Table 5). These are low abundant amino acids and the synthesis of these amino acids is energetically costly and requires multiple reaction steps [30]. With the stress time increasing, the number of genes increased in all the pathways except arginine biosynthesis which remains stable, the increasing is very slowly or stable from 0.5-h of stress to 1-h of stress and sharply from 1-h of stress to 2-h of stress ( Table 5). The sharp increase of gene numbers implies that the synthesis of these low abundant amino acids at a significantly higher rate on 2-h of stress. Ribosome biogenesis and transcription are also significantly  enriched on 2-h of stress (Table 6). Genes involved in the two biological processes are highly expressed on 2-h of stress (Additional file 2: Figure S5 and S6). These together implies the synthesis of new protein on 2-h of stress at a significantly higher rate, compared with that of 0.5-h and 1-h of stress. Spliceosome is significantly enriched on 1-h and 2-h of stress (Table 6; Additional file 2: Figure S7). Post-transcriptional regulation mediated by spliceosome is an important way of gene regulation and was reported to involve in responses to various abiotic stress including salt and temperature stress [31,32]. Alternative splicing is a fast way to generate new mRNAs (splice variants) using the existed pre-mRNAs, and then these splice variants can be translated into different proteins. Protein processing in endoplasmic reticulum is also significantly enriched on 1h and 2-h of stress (Table 6; Additional file 2: Figure S8), which implies the regulation at post-translational level.
On the other hand, ubiquitin mediated proteolysis is significantly enriched on 2-h of stress, which implies the degradation of proteins (Table 6; Additional file 2: Figure S9). The enrichments of protein synthesis and degradation at the same time implies a higher protein turnover rate on 2-h of stress. Some proteins are destined to degradation, such as misfolded proteins caused by stress [33], while new proteins are needed to adapt to the hypersaline environmental, such as chaperones (heat shock proteins) and reactive oxygen species scavenging enzymes (superoxide dismutase, catalase, and peroxiredoxins) [34]. Consistently, the expressions of these genes are up-regulated (Additional file 1: Table S4; Additional file 3).

Ion homeostasis
Surprisedly, the expression of all the possible Na + and Cl − transporters or channels are not up-regulated up to 2-h of stress (Table 7), which implies that the enhancement of exclusion of Na + or Cl − are not needed, probably because the extracellular ions can't go into Dunaliella cell easily. When confronting salt stress, Dunaliella cells undergo rapid shrinking by efflux of water through aquaporin [1]. The cell shrinkage and efflux of water result in the increased solutes concentrations in cytoplasm, and at the same time the glycerol synthesis starts, which contributes to the osmotic balance between cell membrane. From this point of view, if Dunaliella cells can restrict ion permission, the enhanced exclusion of the ions by membrane transporters seems no need, so we can't see the up-regulation of these ion transporters.
Some other metal ion transporters are also not upregulated on 0.5-h and 1-h of stress (Additional file 1: Table S5), while some metal ion transporters such as magnesium, zinc, and boron transporters are upregulated on 2-h of stress. It is possible that the upregulation of these transporters is to uptake metal ion for protein synthesis since Dunaliella cells began to synthesize large amounts of protein on 2-h of stress (Table 6).

Salt stress sensing and signal transduction
Kinetic of glycerol synthesis in Dunaliella shows that the synthesis can be detected within minutes and is independent of protein synthesis [35], which indicates that the expression of the signal receptors and signal transducers might not be changed upon signal stimulation. Therefore, using the GO enrichments of the up-regulated genes to show the involved signaling pathways may not reflect the real truth. So far, very little has been known about the sensing and signal transduction cascades of salt stress in Dunaliella. For this reason, we tried to search the possible cell surface receptors based on the sequence annotation data. Two possible receptors are found: a putative Gprotien coupled seven-transmembrane receptor (GPCR) and a receptor-like serine/threonine protein kinase. Small GTP-binding proteins were identified in D. salina plasma membrane by Mass Spectrometry and were thought to involve in sensing and signaling of salt stress in D. salina [8]. However, the ligands of the two receptors are less likely to be Na + or Cl − . Due to the properties of salt stress, we focused on ion channels. Genes encoding sodium channel, chloride channel and mechanosensitive ion channel are found in the sequence data. One gene is annotated as a sodium channel, which shows sequence similarity with Sodium channel protein 60E from Drosophila melanogaster. Sodium channel protein 60E is a voltage-gated sodium channel, and the opening and closing of which is triggered by changing the ion concentration, and hence the charge gradient, between the two sides of the cell membrane [36]. In rabbit chemoreceptor cells, the activation of voltagegated sodium channels can cause cell membrane depolarization, which will result in Ca 2+ influx by activation of voltage-dependent Ca 2+ channels, so the signal is switched to intracellular Ca 2+ signal [37,38]. No sequences encoding voltage-dependent Ca 2+ channel are found in the sequence data. It cannot be concluded that there are no voltage-dependent Ca 2+ channels in Table 7 The expressions of the possible sodium and chloride ion transporters or channels identified in the transcriptiome transporting ATPase and calcium uniporter are found in the sequence data. Furthermore, sequences encoding Calcium dependent kinases, Ca 2+ /calmodulin-dependent protein phosphatase, Calcium sensing receptor, calmodulin, Calcium/calmodulin-dependent 3′, 5′-cyclic nucleotide phosphodiesterase are also found in the sequence data. These indicate that Ca 2+ plays an important role in signal transduction. Two genes are annotated as voltage gated chloride channel. Since the salt stress was applied by only increasing the concentration of NaCl in medium, we wonder if the chloride channel can sense the salt stress by the influx of Cl − . WNK kinase (serine/threonine-protein kinase) is reported to function as a direct chloride sensor, and several genes encoding WNK are also found in the transcriptome. WNK are activated by reduction in intracellular Cl − concentrations and inactivated by binding of a Cl − . Chloride binding inhibited the autophosphorylation of WNK1, thereby inhibiting kinase activity [39,40]. WNK is also reported to inhibit sodium [41,42] and potassium channel [43,44]. Except for these, how the influx of Cl − triggers intracellular pathways are almost unknown.
Mechanosensitive ion channels will respond to the mechanical deformation, which includes changes in the tension, thickness, or curvature, of the membrane. Mechanosensitive channels respond to membrane tension by altering their conformation between an open state and a closed state [45,46]. When confronted with salt stress, Dunaliella cells shrink rapidly, thus the membrane deformation could cause the opening of mechanosensitive ion channels. One gene showing significant sequence similarity with the mechanosensitive ion channel protein 5 of Arabidopsis thaliana is found; however we don't know whether it is a calcium channel or other ion channel. The applying of salt stress can trigger the voltage-gated sodium channel, the voltage-gated chloride channel and also the mechanosensitive ion channel, it is hard to tell which one is the very first signal receptor. Dunaliella cells undergo rapid shrinkage by efflux of water through aquaporin [1]. So the mechanosensitive ion channel seems unlikely to be the first receptor, but might be the enhancer.
A previous paper showing that MAPK signaling cascade may be involved in the salt stress of Dunaliella by western blot analysis [47]. Many genes on MAPK signaling cascade are also found in the sequence data, including mitogen-activated protein kinase (MAPK), Mitogenactivated protein kinase kinase (MAP 2 K), Mitogenactivated protein kinase kinase kinase (MAP 3 K), Mitogen-activated protein kinase kinase kinase kinase (MAP 4 K), MAP kinase phosphatase 1, MAP kinase phosphatase 5, MAP kinase phosphatase 6, and dual specificity MAP kinase phosphatase, which indicates that MAPK signaling pathway may be involved in salt stress signal transduction of Dunaliella salina.

Discussion
Most of the qPCR results shew up-regulations on at least one of the three time points of hypersaline stress, which is consistent with the NGS results except catalase and acetyl-CoA carboxylase. No up-regulations of catalase and acetyl-CoA carboxylase were detected by qPCR (Additional file 3). Acetyl-CoA carboxylase, which converts cytosol acetyl-CoA into malonyl CoA, the first committed step in the synthesis of fatty acids, is a multisubunit enzyme in the chloroplasts of most plants and algae [48]. The up-regulation of acetyl-CoA carboxylase beta subunit was detected by qPCR. The up-regulation of its beta subunit could possibly imply the enhancement of the holoenzyme activity.
Dunaliella is a good model organism for studying salt tolerance. Omics methods are particularly appropriate to analyze the stress response of Dunaliella at a global level, because Dunaliella is a unicellular alga and there will be no problem caused by different cell types. Previous omics studies can't give a global understanding of the response process of Dunaliella to salt stress due to limited discovery of genes involved or unappropriated time duration of stress [13,14]. Here we reported the transcriptomic profiles on 0.5-h, 1-h and 2-h of hypersaline stress and compared the profiles with that before stress. Our results imply that the response of Dunaliella to hypersaline stress is almost finished on 2-h stress. The reasons are that nuclear origin of replication recognition complex and reproductive process are significantly enriched on 2-h stress, which implies that the algae cell starts to reproduce on 2-h of stress and cells already adapted the new environment. This is in agreement with the previous studies reporting that intracellular glycerol synthesis accomplished and cells recovered their original volume in about 2 h after hypersaline stress [15,16].
When the hypersaline stress is applied, many of the genes classified into photosynthesis show a pulse expression pattern (quick increasing and then quick decreasing), implies that hypersaline stress induces enhancement of photosynthesis, but no constant enhancement of photosynthesis is observed. Carbon fixation from photosynthesis has its upper limit under a given condition. Study reported that only a small part of carbon source of glycerol synthesis comes from CO 2 fixation of photosynthesis, and others are from starch breakdown [49]. However, constant increasing of carbon fixation is observed, this may be induced by NaHCO 3 in the high salt medium when stress is applied. Constant increase of pigment biosynthetic process including heme biosynthetic process is also observed. Pigments such as carotenoids play a role in response to abiotic stress [50], and heme could be used to synthesize heme-containing enzymes, such as catalase and ascorbate peroxidase, which play important roles in detoxification of reactive oxygen species (ROS) [25]. Starch is the main carbon source of glycerol synthesis when hypersaline stress is applied [50]. Consistently, starch and sucrose metabolism is significantly enriched. The new drawn glycerol synthesis pathway includes only 6 enzyme-catalyzing steps to convert starch to glycerol, which is theoretically faster than previous reported pathways [9,51]. Synthesis of unsaturated fatty acid is consistent with previous report [28]. Interestingly, endocytosis is significantly enriched, together with the enrichment of Glycerophospholipid metabolism. We speculate that endocytosis could cause the degradation of excessive plasma membrane and help the synthesis of new glycerophospholipid with unsaturated fatty acids.
From amino acids and protein metabolism, we found an indication of sharp increasing of protein synthesis rate on 2-h stress based on the significant enrichment of ribosome biogenesis and transcription. Enhancement of regulation on protein synthesis is also implied by the significant enrichment of spliceosome and protein processing in endoplasmic reticulum. While, the significant enrichment of ubiquitin mediated proteolysis implies the degradation of proteins. Taken together, accelerated protein turnover rate on 2-h of stress is implied.
It is very interesting that no sodium and chloride channels or transporters are up-regulated during the period of 2-h of stress. No up-regulation of Na + /H + antiporter was detected by qPCR (Additional file 3). This is in consistent with the reports of Na + /H + antiporter, which regards Na + /H + antiporter as the main player to eliminate Na + during salt stress [7,52]. We can't see the up-regulation of Na + /H + antiporter, maybe the activity of Na + /H + antiporter is regulated at protein level instead of mRNA level. However, if Dunaliella cell can restrict Na + permission, the exclusion of the ions by membrane transporters seems no need.
The time setting of stress gives us good results. We can see the constant increasing of the number of upregulated genes with the increasing of stress time. Many of the results of functional enrichment are consistent with the previous studies, such as the enhancement of photosynthesis and carbon fixation [9], the degradation of starch for glycerol synthesis [49], and the desaturation of membrane lipid [28] when hypersaline stress is applied.
We also have some new findings, such as the enhancement of heme biosynthetic process and endocytosis, the accelerated protein turnover rate, the enhancement of spliceosome and protein processing in endoplasmic reticulum in response to hypersaline stress. These new findings shed light on the mechanism of salt tolerance of photosynthetic plants and may help improving salt tolerance of crops by genetic manipulation.
Finally, we like to acknowledge the limitation of this study, that is, all the analyses are based on the transcriptomic data, the changes at mRNA level. Although protein levels are usually positively correlated with the mRNA levels, there are post-transcriptional and translational regulations that may affect the correlation. From this point of view, it is necessary to further figure out how a new intracellular balance at the physiological levels achieved in Dunaliella salina to adapt to hypersaline environment.

Conclusion
At global transcriptomic level, the response of Dunaliella salina to hypersaline stress is systematic, possibly involving enhancement of photosynthesis, constant increasing of carbon fixation and heme biosynthetic process, degradation of starch, synthesis of glycerol, membrane lipid desaturation, accelerated protein turnover, enhancement of spliceosome, protein processing in endoplasmic reticulum, and endocytosis. The changes of these biological processes will help to understand the achieved new balance adapting to the hypersaline environment.  [53]. Total 1200 mL algae in log phase (about 10 6 cells/mL) were divided into 12 bottles, 100 mL per bottle. Three bottles were used for control; the other 9 bottles were used for hyper-salt stress treatment. For hyper-salt stress treatment, each of the 9 bottles was added directly with 100 mL of high salt medium (containing 4 M NaCl) to form a final NaCl concentration of 2.5 M and was further cultured in light at 20°C. Algae cells under stress were harvested on time of 0.5-h, 1-h, and 2-h for total RNA extraction. The control algae cells were also harvested for total RNA extraction. In all, there are 4 time points (4 groups), 0 h, 0.5 h, 1 h and 2 h, each time point with three replicates. In total 12 RNA samples were prepared from 12 algae samples.

RNA extraction
Total RNA was extracted with TRIzol (Invitrogen, USA) by following the user manual. RNA degradation and contamination was monitored on 1% agarose gels. RNA purity and integrity was checked by NanoPhotometer (IMPLEN, Germany) and RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, USA). RNA concentration was measured by Qubit® 2.0 Flurometer together with Qubit® RNA Assay Kit (Life Technologies, USA).

Hiseq libraries construction for the second generation sequencing
The 12 RNA samples were used to construct 12 Hiseq libraries, one library per RNA sample. Three microgram total RNA per RNA sample was used to prepare the sequencing library with NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, USA) by following the manufacturer's instruction. Briefly, mRNA was purified by using poly-T magnetic beads and was fragmented to perform first-strand cDNA by using M-MuLV reverse transcriptase and random hexamer-primers. Second-strand cDNA was synthesized by using DNA polymerase I and RNase H. The cDNA was end-repaired and A-tailed, and NEB-Next Adaptor with hairpin loop structure was linked. After size selection with AMPure XP system (Beckman Coulter, Beverly, USA) and PCR amplification with Phusion High-Fidelity DNA polymerase, samples were sequenced on the Illumina Hi-Seq 2000 system, generating paired-end (PE) reads with a length of 150 bp.

PacBio library construction for the third generation sequencing
The only reported Dunaliella salina (CCAP19/18) genome shows low sequence similarity with our algae strain (CCAP 19/30) and can't be used as the reference genome for mRNA sequence assembly and annotation, so the third generation sequencing was applied to generate the full length mRNA sequences. One microgram mRNA (generated by equally mixing mRNAs from the 9 stress samples and the 3 control samples) were used for first-strand cDNA synthesis with SMARTer PCR cDNA Synthesis Kit (Clontech Cat. 634,926). Double-strand cDNA was amplified with PrimeSTAR GXL DNA Polymerase (Clontech R050B). The > 4 kb cDNA fraction was generated with the BluePippin Size Selection System and mixed with the cDNA before fractionation at 1:1 ratio. The mixed cDNA was used to construct a combined SMRTbell Libray for sequencing. The library was subsequently sequenced on the PacBio Sequel System platform using SMRT Cell 1 M v2.

Data processing and transcriptome assembly
Firstly, circular consensus sequences (CCS) were generated from Pacbio raw reads. Then, full length non chimera (FLNC) reads were generated from CCS. FLNC reads should have sequence regarding to both 5′ and 3′ Clontech kit primers as well as a poly-A tail signal preceding the 3′ primer. FLNC reads were clustered by ICE (Iterative isoform-clustering) and a consensus sequence was generated from each cluster. The consensus sequences were corrected by non-full length sequences to generate the polished consensus sequences.
Hiseq raw reads were processed to generate Hiseq clean reads by removing adapter sequences, excluding reads contains > 50% nucleotides with Qphred ≤20, and reads with ambiguous bases N. Hiseq clean reads were used to correct the polished consensus sequences obtained from the third generation sequencing. The correction process was done by Proovread [54]. The final corrected unique full length sequences were used for annotation and gene expression analysis.

Coding sequence analysis and gene annotation
Transcripts were firstly submitted to Nr and Swissprot protein database to analyze the potential open reading frames (ORF), then the transcripts failed to hit any proteins in the two databases were submitted to the software ESTSCAN (3.0.3) [55] to analyze the potential ORF, finally the transcripts failed the ESTSCAN are defined as non-coding sequences. The corrected unique full length sequences were annotated against 7 databases, Nr (NCBI's non-redundant protein sequences), Nt (NCBI's nucleotide sequences), Pfam (protein family), KOG/COG (clusters of orthologous groups of proteins/ eukaryotic ortholog groups), Swiss-Prot, KEGG (Kyoto encyclopedia of genes and genomes) and GO (gene ontology) database. The softwares used are diamond v0.8.22 for Nr (e-value = 1e-5), KOG/COG (1e-3) and Swiss-Prot (e-value = 1e-5); NCBI blast 2.2.28 for Nt (evalue = 1e-5); HMMER 3.0 for Pfam (e-value = 0.01); KEGG automatic annotation server (KAAS) for KEGG (e-value = 1e-10), Blast2GO v2.5 for GO (e-value = 1e-6). When available, the functional annotation from Swiss-Prot was preferred because of its high accuracy.

Quantification of gene expression level
Briefly, Hiseq clean reads from each replicate were mapped to the corrected unique full-length sequences using HTSeq v0.6.1 [56]. And then FPKM (expected number of Fragments per Kilobase of transcript sequence per Millions base pairs sequenced) of each gene was calculated based on the length of the gene and reads count mapped to this gene. FPKM considers the effect of sequencing depth and gene length at the same time, and is currently the most commonly used method for estimating gene expression levels [57]. Clustering analysis of gene expression value was performed by SOM method [58].

Differential expression assessment
Differential expression analysis of two groups was performed by using the DESeq R package with the read counts from HTSeq [59]. DESeq provides statistical routines for determining differential expression in digital gene expression data, using a model based on the negative binomial distribution. The resulting P values were adjusted using the Benjamini and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted P-value < 0.05 found by DESeq were assigned as differentially expressed.

Gene enrichment analysis
Gene Ontology (GO) terms were assigned to each genes based on the BLAST2GO [60]. Gene Ontology (GO) enrichment analysis of differentially expressed genes were implemented by the GOseq R package [61], in which gene length bias was corrected. GO terms with corrected P value < 0.05 were considered significantly enriched by differential expressed genes.
KEGG pathways were assigned to each genes based on the KAAS. KOBAS software [62] was used to test the statistical enrichment of differential expression genes in KEGG pathways. The enrichment analysis was tested using a hypergeometric test at a significance cutoff of 0.05 false discovery rate (FDR). Pathways were assigned as significantly enriched if they had a FDR below 5%.

Real-time PCR (qPCR)
In order to confirm the expression results from RNAsequencing, we perform qPCR analysis of the key genes from each functional group. The qPCR kit was from TaKaRa (DRR081A, TaKaRa Biotechnology, and Dalian, China). The primers used in qPCR were in the supplementary materials (Additional file 3). The reference gene used in qPCR was 18S rRNA. The relative expression quantification was calculated by the 2 -ΔΔCt method.
Additional file 1: Table S1. The data qualities of the libraries of the second-generation sequencing. Table S2. Transcripts properties from the third-generation sequencing. Table S3. GO terms enriched from the down-regulated genes. Table S4. Up-regulation of some key genes on different stress time compare with that before stress. Table S5. The expressions of the possible metal ion (excluding sodium and chloride ion) transporters or channels identified in the transcriptome.
Additional file 2: Figure S1. Clustering analysis of the differentially expressed genes in each sample. Figure S2. Heat-map of carbon fixation. Figure S3. heat-map of tetrapyrrole biosynthetic process. Figure S4. Heat-map of fatty acid biosynthesis. Figure S5. Heat-map of ribosome biogenesis. Figure S6. Heat-map of transcription. Figure S7. Heat-map of spliceosome. Figure S8. Heat-map of protein processing in endoplasmic reticulum. Figure S9. Heat-map of ubiquitin mediated proteolysis.