Identification of Salt Stress Responding Genes Using Transcriptome Analysis in Green Alga Chlamydomonas reinhardtii

Salinity is one of the most important abiotic stresses threatening plant growth and agricultural productivity worldwide. In green alga Chlamydomonas reinhardtii, physiological evidence indicates that saline stress increases intracellular peroxide levels and inhibits photosynthetic-electron flow. However, understanding the genetic underpinnings of salt-responding traits in plantae remains a daunting challenge. In this study, the transcriptome analysis of short-term acclimation to salt stress (200 mM NaCl for 24 h) was performed in C. reinhardtii. A total of 10,635 unigenes were identified as being differently expressed by RNA-seq, including 5920 up- and 4715 down-regulated unigenes. A series of molecular cues were screened for salt stress response, including maintaining the lipid homeostasis by regulating phosphatidic acid, acetate being used as an alternative source of energy for solving impairment of photosynthesis, and enhancement of glycolysis metabolism to decrease the carbohydrate accumulation in cells. Our results may help understand the molecular and genetic underpinnings of salt stress responses in green alga C. reinhardtii.


Introduction
Salinity is one of the most important abiotic stresses threatening agricultural productivity worldwide. Although plants have gradually evolved a series of adaptive molecular, physiology and biochemistry processes to respond to salinity stress, it could threaten 30% of cultivable soils by 2050 [1,2]. Understanding the molecular machineries of salt stress response in model plants of basal taxa, such as green algae, may contribute to finding the evolutionary cues of abiotic stress response in plants and developing salt-resistant crops with additional salt-responding traits [2][3][4][5][6][7][8][9].
Salt stress causes diverse impacts on plant growth by disturbing the osmotic/ionic balance and eliciting Na + toxicity [9,10]. Under aquatic saline stress, a series of physical and biochemical processes are recruited by algae to respond to the damage caused by osmotic and ionic stresses, such as photosynthesis inhibition, macromolecular compound synthesis and homeostasis adjustment [6,[10][11][12][13][14]. It has been reported that salt stress leads to decreased photosynthetic efficiency [15,16] which influences chlorophyll content in plant leaves [17,18]. In green algae, salt stress remarkably influences the structure and functions of the photosynthetic apparatus in Scenedesmus obliquus [19] and reduces the maximum quantum yield of photosystem II (PSII) in Dunaliella salina [20]. In alga Botryococcus braunii, metabolism of lutein was significantly enhanced under stress conditions [12].
Chlamydomonas reinhardtii is a free-living freshwater alga with unicellular vegetative cell. Previous studies exposed the C. reinhardtii strain 21 gr and CC-503 to salt stress and demonstrated the physiological and metabolic processes impacted by ionic toxicity and osmotic stress caused by salt damage [6][7][8]16,21]. Vega [22] demonstrated that 200 mM NaCl in the culture medium was highly toxic for C. reinhardtii productivity. The addition of NaCl immediately blocked the photosynthetic activity of the alga which partially recovered, after 1 h of treatment, remaining high during the following 24 h. However, after 24 h treatment with NaCl 200 mM, the intracellular catalase activity of the alga reached a 20-fold higher level than in the control cells. The physiological data indicate that saline stress induces in C. reinhardtii an increase of intracellular peroxide, which parallels a significant inhibition of the photosynthetic-electron flow. However, the related machineries of up-stream regulating and the triggering of appropriate cellular and physiological responses to cope with stress circumstances are still largely unknown.
Keeping this in mind, the work presented here was carried out to explore the saline stress-responding mechanisms of C. reinhardtii by transcriptome sequencing of strains GY-D55 wild type. The aim of this study was to identify dys-regulated genes in C. reinhardtii cells under salt stress by RNA-seq, screen physiological and biochemical cues by gene ontology (GO) terms and MapMan functional enrichment analyses, and investigate the physiological adaptions and cellular regulatory networks for salt stress responding.

Transcriptome Profiling of C. reinhardtii
After sequencing with the Illumina HiSeq X platform, a total of 56,438,218, 72,853,712, 47,551,786, 56,962,722, 52,926,804 and 55,998,748 high-quality pair-end reads were obtained from three control and three salt stress treated samples of C. reinhardtii (Table 1), respectively. De novo transcriptome assembly generated 91,242 unigenes, with an average length of 2691 nt and N50 of 4554. On average, 90.66% of the reads from six samples were mapped to the reference genome ( Table 1). The assembled transcriptome information of C. reinhardtii is shown in Supplementary Figure S1.

Functional Annotations of Unigenes
Similarity searches were performed to annotate unigenes against different databases using BLASTX. For C. reinhardtii, 65,679 (71.98%) unigenes were annotated in at least one database ( Figure 1C and Supplementary Figure S1). A total of 52,884 (57.96%) and 58,062 (63.63%) unigenes showed similarity to sequences in NR and PFAM databases with an E-value threshold of 1 × 10 −5 . About 58,651 (64.28%) unigenes were annotated in the GO database by Blast2GO v2.5 with an E-value cutoff of 1 × 10 −6 ( Figure 1C and Supplementary Figure S1). Unigenes of the C. reinhardtii were assigned to C. reinhardtii and A. thaliana gene IDs for GO annotation mapping and TFs/PKs perdition. By sequence alignment, a total of 48,158 unigenes were aligned to C. reinhardtii PLAZA genome genes. A total of 54,509 unigenes were assigned to TAIR10 locus IDs by BLASTP with an E-value cutoff of 1 × 10 −5 and classified into GO categories for GO analysis (Supplementary Table S1).  Table S1). Red and green dots represent DEGs, blue dots indicate genes that were not differentially expressed. In total, 10,635 unigenes were identified as DEGs (padj < 0.05) between S_200 and C_0, including 5920 upregulated genes and 4715 downregulated genes.

Differently Expressed Genes (DEGs) Calculation
To evaluate the relative level of gene expression in C. reinhardtii under control or salt stress treatment, the FPKM values were calculated based on the uniquely mapped reads. The FPKM distributions of unigenes in six samples are shown in Supplementary Figure 2. The FPKM value for genes detected in six samples ranged from 0 to 40,486.05, with mean value of 7.08. By comparative analysis, a part of the genes was observed to be differently expressed in 200 Mm NaCl treated samples: 5920 unigenes were calculated as up-regulated in salt treated samples and 4715 filtered as down-regulated genes with the cutoff of padj < 0.05 and |log2(foldchange)| > 1 (Supplementary Table S2). Red and green dots represent DEGs, blue dots indicate genes that were not differentially expressed. In total, 10,635 unigenes were identified as DEGs (padj < 0.05) between S_200 and C_0, including 5920 upregulated genes and 4715 downregulated genes.

Differently Expressed Genes (DEGs) Calculation
To evaluate the relative level of gene expression in C. reinhardtii under control or salt stress treatment, the FPKM values were calculated based on the uniquely mapped reads. The FPKM distributions of unigenes in six samples are shown in Supplementary Figure S2. The FPKM value for genes detected in six samples ranged from 0 to 40,486.05, with mean value of 7.08. By comparative analysis, a part of the genes was observed to be differently expressed in 200 Mm NaCl treated samples: 5920 unigenes were calculated as up-regulated in salt treated samples and 4715 filtered as down-regulated genes with the cutoff of padj < 0.05 and |log2(foldchange)| > 1 (Supplementary Table S2).

MapMan Enrichment of DEGs
A more specific comparison of metabolic and regulatory pathways was conducted using MapMan. A total of 5920 up-and 4715 down-regulated genes were assigned to 1334 and 1050 homologs in Arabidopsis thaliana, respectively. Consequently, these uniquely expressed genes were mapped to 797 pathways by MapMan, of which, 22 pathways were filtered enriched by the dysregulated genes with the cutoff p-value < 0.05 ( Figure 3A; Supplementary Table S4). The expression of genes implicated in "TCA/org. transformation.TCA", "Tetrapyrrole synthesis", "Starch" and "Sucrose" were over-expressed in C. reinhardtii, while those genes involved in "PS.lightreaction", "PS.lightreaction.photosystem I" and "PS.lightreaction.photosystem I.LHC-I" were down-regulated in C. reinhardtii during salt stress responding ( Figure 3A).

KEGG Enrichment of DEGs
To gain a deeper insight into the regulation of photosynthesis underlying salt stress response, down-regulated unigenes involved in "photosynthesis" KEGG pathways (ko00195) were mapped and shown in Figure 2B. Orthologs of 44 genes annotated in this pathway were filtered as down-regulated in the NaCl treated samples in the green alga, such as, photosystem II oxygen-evolving enhancer protein

The Differentially Expressed TFs and PKs
Among the expressed unigenes, 2050 and 1624 sequences were assigned to 45 TF families and 78 PK families, respectively (Supplementary Table S5). Of the TF families, MYB family had the largest number of upregulated genes (16 unigenes), including MYB109 ortholog unigenes Cluster-2749.35807 (L 2 fc = 2.5798) and Cluster-2749.70085 (L 2 fc = 1.3722). In contrast, SET family had the largest number of downregulated genes (16 unigenes). Of the PKs families, TKL-Cr-3 family was uncovered to contain the largest number of upregulated genes. By comparison, CAMK_CDPK and Group-Cr-2 family contained the largest number of downregulated unigenes (Supplementary Table S5).

Real-Time Quantitative PCR Validation
To verify the RNA-seq results, an alternative strategy was selected for the upregulated unigenes. In total, five over-expressed unigenes were randomly selected for validation by qRT-PCR using the same RNA samples that were used for RNA-seq. Primers were designed to span exon-exon junctions (see Supplementary Table S6 and Figure S3). In most cases, the gene expression trends were similar between these two methods; the result is shown in Figure 3. The ortholog of cytosolic small heat shock protein encoding genes HSP17.6A, Cluster-2749.57700, which was detected by RNA-Seq as up-regulating genes in the salt treated samples (L 2 fc = 9.76), was also detected to be significantly over-expressed by qRT-PCR method (Figure 3).

Real-Time Quantitative PCR Validation
To verify the RNA-seq results, an alternative strategy was selected for the upregulated unigenes. In total, five over-expressed unigenes were randomly selected for validation by qRT-PCR using the same RNA samples that were used for RNA-seq. Primers were designed to span exon-exon junctions (see Supplementary Table S6 and Figure S3). In most cases, the gene expression trends were similar between these two methods; the result is shown in Figure 3. The ortholog of cytosolic small heat shock protein encoding genes HSP17.6A, Cluster-2749.57700, which was detected by RNA-Seq as up-regulating genes in the salt treated samples (L2fc = 9.76), was also detected to be significantly over-expressed by qRT-PCR method ( Figure 3).

Discussion
Salinity is one of the major environmental factors threatening crop productivity and plant growth worldwide [2,9,42]. Due to the complexity of abiotic stress-responding processes, although several hundreds of salt-responding genes have been reported in plants, understanding the genetic underpinnings of salt-responding traits in plantae remains a daunting challenge. The model alga C. reinhardtii, which contains one large cup-shaped chloroplast, has the ability to adapt rapidly to changing environmental conditions, such as high salinity, via the generation of novel traits [8,14,43,44]. Given previous results from analysis of salt stress in C. reinhardtii and other plants, we analyzed the Illumina RNA-seq data from this alga grown in BG-11 medium with the addition of 200 mM NaCl and analyzed in triplicate after 24 h of incubation [16,22]. In this study, a total of 5920 and 4715 unigenes were identified as up-and down-regulated genes in C. reinhardtii under salt stress by RNA-seq. Our study found some molecular cues for reducing the negative effects due to ionic/osmotic toxicity and photosynthesis impairment under saline conditions in C. reinhardtii.
Previous studies discovered that the cell density of C. reinhardtii cells obviously reduced when stressed by NaCl [

Discussion
Salinity is one of the major environmental factors threatening crop productivity and plant growth worldwide [2,9,42]. Due to the complexity of abiotic stress-responding processes, although several hundreds of salt-responding genes have been reported in plants, understanding the genetic underpinnings of salt-responding traits in plantae remains a daunting challenge. The model alga C. reinhardtii, which contains one large cup-shaped chloroplast, has the ability to adapt rapidly to changing environmental conditions, such as high salinity, via the generation of novel traits [8,14,43,44]. Given previous results from analysis of salt stress in C. reinhardtii and other plants, we analyzed the Illumina RNA-seq data from this alga grown in BG-11 medium with the addition of 200 mM NaCl and analyzed in triplicate after 24 h of incubation [16,22]. In this study, a total of 5920 and 4715 unigenes were identified as up-and down-regulated genes in C. reinhardtii under salt stress by RNA-seq. Our study found some molecular cues for reducing the negative effects due to ionic/osmotic toxicity and photosynthesis impairment under saline conditions in C. reinhardtii.
Previous studies discovered that the cell density of C. reinhardtii cells obviously reduced when stressed by NaCl [8,14,16,21,22]. Neelam et al. demonstrated that at the morphological level, 150 or 200 mM NaCl salt stress led to palmelloid morphology, flagellar resorption, reduction in cell size, and slower growth rate in C. reinhardtii [21]. It should be noted that dead and dying cells have dys-regulated mRNA and contribute to transcript levels under saline stress. In our study, programmed cell death (PCD) in the C. reinhardtii cell was found with PCD-regulating proteins being significantly up-regulated, e.g., condensin complex subunit (Cluster-2749.11751: L 2 fc = 9.368; Cluster-2749.12889: L 2 fc = 7.766), sucrose-phosphatase 1 (Cluster-2749.35394: L 2 fc = 4.304) and stress tolerance related fibrillin family member (Cluster-2749.70284: L 2 fc = 1.920).
Saline stress leads to the overproduction of reactive oxygen species (ROS) in plants which are highly reactive and toxic and cause damage to lipids, carbohydrates, proteins and DNA which ultimately results in oxidative stress [8,9,14,45]. The accumulation of ROS also influences the expression of a number of genes and therefore controls many processes, such as growth, cell cycle, PCD, secondary stress responses and systemic signaling [8,9,14,45]. The excess Na + and oxidative stress in the intracellular or extracellular environment activates the acytoplasmic Ca 2+ signal pathway for regulating an osmotic adjustment or homeostasis regulating of salt stress responses [24,29,39,[46][47][48][49][50][51]. In our study, calcium-related pathway in the C. reinhardtii cell was found with several calcium ion binding proteins being significantly upregulated, e.g., peroxygenase 3 (Cluster-2749.55812: L 2 fc = 10.431; Cluster-2749.59997: L 2 fc = 7.680) and calreticulin (Cluster-2749.35394: L 2 fc = 3.082).
The short-term (within 48 h) acclimation to salt stress in C. reinhardtii involves activation of phospholipid signaling, leading to the accumulation of phosphatidic acid (PA), which is a lipid second messenger in plant and animal systems [52][53][54]. In the case of C. reinhardtii, incubation in 150 mM NaCl leads to a three-to four-fold rise of PA levels within minutes [52,55]. Lysophosphatidic acid (LPA) has also been shown to accumulate in this alga under salt stress, with the dose-dependent response reaching a maximum at 300 mM NaCl [55,56]. In this study, soluble lysophosphatidic acid acyltransferase (Cluster-2749.8269: L 2 fc = 9.126; Cluster-2749.9895: L 2 fc = 8.720) was found to be significantly up-regulated in salt stress treated samples, which indicated the potential role of this gene in maintaining the lipid homeostasis by regulating PA under saline stress [55]. Further, analysis of glycerophospholipid metabolism pathways showed that the alga cells had significant up-regulation of FAD (flavin adenine dinucleotide)-dependent oxidoreductase family protein (Cluster-2749.52046; L 2 fc = 2.695) that involves storing lipid catabolism and glycerol assimilation, and in glycerol-3-phosphate shuttle, which transports reduced power from cytosol to mitochondrion [8]. This suggests that the intracellular glycerol pool in C. reinhardtii cells likely increased as a response to salt stress, similar to what has been shown for the green alga Dunaliella tertiolecta [57,58].
Requirement of energy to maintain ion homeostasis is the major metabolic impact of salt stress. The reduction of oxidative stress and osmotic stress, and the up-regulation of heatshock proteins were speculated to aid protein renaturation and recover homeostasis [59][60][61]. In this study, the stress response is apparent in the C. reinhardtii cells with significant up-regulation of genes involved in oxidative/osmotic stress reduction process including glyceraldehyde-3-phosphate dehydrogenase C subunit 1 (Cluster-2749.27769: L 2 fc = 1.930) and fumarase 1 (Cluster-2749.35832: L 2 fc = 8.306). In bacterium Escherichia coli, trehalose is synthesized as a compatible solute and enables cells to exclude toxic cations and to acclimate to high concentrations of salt in the growth medium [62]. For maize, trehalose has helped to reduce the negative effects of saline stress as an osmoprotectant [63]. In our study, enzymes involved in trehalose synthesis significantly up-regulated, e.g., trehalose-6-phosphatase synthase S8 (Cluster-2749.17684: L 2 fc = 6.453) and trehalose-6-phosphate synthase (Cluster-2749.61951: L 2 fc = 1.123). These results indicated the potential underpinnings for these to maintain homeostasis in C. reinhardtii under saline conditions. In plants, saline stress generally causes ion injury and osmotic stress, which interferes with numerous biochemical and physiological processes, including energy metabolism pathways such as photosynthesis [26,36,64,65] and photorespiration [8]. Previous pigment analyses have demonstrated that photosystem I-light harvesting complexes (LHCs) are damaged by ROS at high salt conditions, and PSII proteins involved in oxygen evolution are impaired [21,45]. In our study, impairment of photosynthesis in the C. reinhardtii cell population was found, with several photosystem I-light harvesting complex (LHC) proteins being significantly down-regulated ( Figure 2B), e.g., photosystem I light harvesting complex gene LHCA2 (Cluster-2749.32743: L 2 fc = −4.74; Cluster-2749.52511: L 2 fc = −3.28), LHCA3 (Cluster-2749.43129: L 2 fc = −6.583) and LHCA5 (Cluster-2749.40312: L2fc = −11.375; Cluster-2749.34085: L 2 fc = −6.553). Further, we found most of the chloroplast encoded transcripts (e.g., PsaA, B, C, J, M) in photosystem I (PSI) were relatively unchanged in level while the nuclear genes (e.g., PsaD, E, G, F, H) down-regulated under saline conditions ( Figure 2B). Existing studies have demonstrated the usage of acetate in the medium as alternative source of energy to compensate for the lowered efficiency in photosynthesis [66]. Consistent with this view, we found that acetyl-CoA synthetase (Cluster-2749.60516: L 2 fc = 5.144; Cluster-2749.25511: L 2 fc = 2.495), which combines acetate and CoA to form acetyl-CoA, was significantly up-regulated in the alga cells under saline conditions. In this study, a significant down-regulation was found in a key enzyme of the glyoxylate cycle-isocitrate lyase (ICL, [Cluster-2749.51492; L 2 fc = −3.119]) [8,45,66,67]-which catalyzes the cleavage of isocitrate to succinate and glyoxylate. Together with malate synthase, ICL bypasses the two decarboxylation steps of the tricarboxylic acid cycle (TCA cycle) [8]. The spatio-temporal expression patterns of genes suggest that in alga cells acetyl-CoA is introduced into energy generation pathways for salt stress responses.
Glycolysis is considered to play an important role in plant development and adaptation to multiple abiotic stresses, such as cold, salt, and drought. It is the key respiratory pathway for generating ATP and carbohydrates metabolites [50,[68][69][70][71][72][73]. In our work, salt stress significantly increased the expression of genes participating in the metabolism of main carbohydrates, such as starch, sucrose, soluble sugar and glucose (Figure 2A). For example, 31 genes of "glycolytic process" (GO:0006096) over-expressed during salt stress responding, including plastidic pyruvate kinase PKP-ALPHA (Cluster-2749.14688: L 2 fc = 8.53) and PKP-BETA1 (Cluster-2749.26182: L 2 fc = 3.68). This is consistent with Zhong et al. [68], who reported salt stress significantly increased the main carbohydrate contents of cucumber leaves [53]. Carbohydrates are involved not only in osmotic adjustment, but also can be used as protective agents for homeostasis regulating during salt stress tolerance [24,30,39,48,69,70,[74][75][76][77]. Given that salt injury caused the destruction of photosynthesis, which might inhibit transport of carbohydrate and accumulate excess starch or sucrose, we speculate C. reinhardtii enhanced glycolysis metabolism to decrease carbohydrate accumulation in cells, which would promote the respiratory metabolism and mitochondrial electron transport, thus reducing the effects of ionic toxicity and osmotic stress caused by salt damage.

Chlamydomonas Material Preparation, Salt Stress Treatment and RNA Extraction
The C. reinhardtii strain GY-D55 wild type from LeadingTec (Shanghai, China) were grown in 150 mL of BG11 media, and placed on a shaking table with 120 rpm and maintained at light (16 h)/dark (8 h) at 23 • C, with an illumination of 100 µmol m −2 ·s −1 . The density of cell cultures was determined by using the blood cell counting plate, with each value being the means of 6 repeats. Under this condition, C. reinhardtii cells were grown in BG11 for 14 d.
The methods published by Zhao [16] and Vega [22] were referenced for NaCl treatment in this study. A total of 50 mL medium with 800 mM NaCl was added to the 150 mL culture medium on a shaking table for finishing 200 mM NaCl treatment, the added NaCl was rapidly diluted, and then the pH value was adjusted to 7.0. A parallel set of cells that were unexposed to NaCl stress conditions and cultured in medium served as the experimental control. A total of 50 mL medium without NaCl was added into the control group. Each treatment had 3 repeats. For 24 h, 200 mM NaCl treatment significantly affected the cellular physiology of the alga, such as its photosynthetic and intracellular catalase activity; in this study, the culture time for C. reinhardtii under salt stress was 24 h.
After 24 h, 100 mL cell culture medium was extracted from the NaCl treated and control culture bottles, respectively. The collected cells were centrifuged at 3000× g for 5 min, and the collected cells were resuspended in 25 mL RNAlater (Ambion, Shanghai, China) solution for RNA extraction. The cells of each repeat were mixed and total RNAs were extracted separately using the TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) following the manufacturer's procedures. RNA quality was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA) and the NanoDrop 2000 spectrophotometer (Thermo Scientific, Wilmington, NC, USA).

Illumina Library Construction and Sequencing
A total amount of 1.5 µg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext ® Ultra TM RNA Library Prep Kit for Illumina ® (NEB, San Diego, CA, USA) by following manufacturer's procedures, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. The random hexamer primer and M-MuLV Reverse Transcriptase (RNase H − ) were used to synthesize the first strand cDNA and the DNA Polymerase I and RNase H were used for second strand cDNA synthesis. Fragments of 150~200 bp cDNA were purified with the AMPure XP system (Beckman Coulter, Beverly, MA, USA). Then, 3 µL USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 • C for 15 min followed by 5 min at 95 • C. Then, PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. Ten cycles were used for PCR enrichment. Finally, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia, San Diego, CA, USA) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq X platform (Illumina, San Diego, CA, USA), according to the manufacturer's procedures. All genetic data have been submitted to the NCBI Sequence Read Archive (SRA) database (https://www.ncbi.nlm.nih.gov/sra), SRA accession: PRJNA490089.

De Novo Transcriptome Assembling and Unigene Annotation
RNA sequencing and de novo transcriptome assembling were conducted to create reference sequence libraries for C. reinhardtii. The RNA sample of each repeat was sequenced separately. cDNA library construction and Illumina pair-end 150 pb sequencing (PE150) were performed at Novogene Co., Ltd. (Shanghai, China), according to instructions provided by Illumina Inc. Reads containing adapter, ploy-N and low-quality reads were removed from raw data for obtaining clean reads. The filtered high-quality reads were used for transcriptome assembling by the Trinity software with default parameters [78]. Clean datasets of 6 samples were pooled for de novo assembling and comprehensive sequence library construction. The Basic Local Alignment Search Tool (BLAST) searches of de novo assembled sequences against public databases (NR, NT, Swiss-Prot, Pfam, KOG/COG, Swiss-Prot, KEGG Ortholog database and Gene Ontology) with an E-value threshold of 10 −10 were used for unigenes' annotation.

Calculation and Comparison of Unigene Expression
The independent transcripts libraries of 3 repeats under NaCl treatment conditions and 3 under control conditions were generated for C. reinhardtii by a PE150 sequencing analysis. The clean reads were aligned to the de novo assembled transcriptome and estimated by the RSEM [79] method. Gene expression levels were calculated by the fragment per kilobase of exon model per million mapped reads (FPKM) method. DESeq2 [80] was used to compare the expression levels between NaCl treated and control samples with an cutoff of adjusted p-value (padj) < 0.05 and |log2(foldchange)| > 1.

GO, KEGG and MapMan Annotation and Enrichment
The GO enrichment analysis for DEGs of C. reinhardtii was performed by the topGO package of R. KEGG [82] is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS [83] software to test the statistical enrichment of differential expression genes in KEGG pathways. MapMan (version 3.5.1 R2) [84] was also used to annotate the DEGs onto metabolic pathways. The DEGs of C. reinhardtii unigene IDs were transferred to the Arabidopsis Information Resource (TAIR) locus IDs during the MapMan analysis.

Real-Time Quantitative PCR (qRT-PCR) Verification
Real-time quantitative PCR (qRT-PCR) was performed to verify the expression patterns revealed by the RNA-seq analysis. The purified RNA of samples under salt stress and control conditions were treated with DNaseI and converted to cDNA using the PrimeScript RT Reagent Kit with gDNA Eraser (Takara, Dalian, China) according to the manufacturer's procedures. Five up-regulated unigenes in C. reinhardtii were selected for the qRT-PCR assay, including Cluster-2749.49004 (ortholog of HSP81-2), Cluster-2749.57700 (ortholog of HSP17.6A), Cluster-2749.55812 (ortholog of RD20), Cluster-2749.36436 (ortholog of RGP), and Cluster-2749.35807 (ortholog of MYB109). Gene-specific qRT-PCR primers (18-20 bp) (Table S6) were designed using Premier 5.0 software. qPCR was performed using SYBR Green qPCR Master Mix (DBI, Ludwigshafen, Germany) in ABI7500 Real-Time PCR System (ABI, Waltham, MA, USA). Three replicates were performed, and the amplicons were used for melting curve analysis to evaluate the amplification specificity. Relative gene expression was quantified using the 2 −(∆∆Ct) method [85]. Ortholog of the A. thaliana housekeeping GTP binding Elongation factor Tu family member AT5G60390 in C. reinhardtii (Cluster-2749.43263) was used to normalize the amount of template cDNA added in each reaction.

Conclusions
We performed a transcriptome analysis of short-term acclimation to salt stress (200 mM NaCl for 24 h) in C. reinhardtii. In total, 10,635 unigenes were identified as differentially expressed in C. reinhardtii under salt stress by RNA-seq, including 5920 that were up-and 4715 that were down-regulated. A series of molecular cues were screened by GO terms, MapMan and KEGG functional enrichment analyses, which were identified as potential mechanisms for salt stress responses. These mainly include maintaining the lipid homeostasis by regulating phosphatidic acid, acetate being used as an alternative source of energy for solving impairment of photosynthesis and enhancement of glycolysis metabolism to decrease the carbohydrate accumulation in cells. Our results may help understand the molecular and genetic underpinnings of salt responding traits in green alga C. reinhardtii.
Supplementary Materials: The following are available online at http://www.mdpi.com/1422-0067/19/11/3359/s1, Figure S1: The assembled transcriptome information of C. reinhardtii, Figure S2: The FPKM density distribution of C. reinhardtii, Table S1: Unigenes of C. reinhardtii annotated in A. thaliana genome by BLASTX analysis, Table S2: Information of the 5920 up-and 4715 down-regulated unigenes in C. reinhardtii, Table S3: The Gene Ontology (GO) enrichment results of the dys-regulated genes in C. reinhardtii, Table S4: The MapMan pathways enrichment results of the dys-regulated genes in C. reinhardtii, Table S5: Information of the differently expressed transcription factors (TFs) and protein kinases (PKs), Table S6: Information of the qRT-PCR primers.
Author Contributions: L.Z. and N.W. conceived and designed the study. L.Z., N.W., Z.Q., M.L., X.Z. performed the data collection and analysis. L.Z. and N.W. wrote the paper. L.Z., S.F. and X.Z. reviewed and edited the manuscript. All authors read and approved the manuscript.
Funding: This work was supported by National Natural Science Foundation of China (31800185) and A Project of Shandong Province Higher Educational Science and Technology Program (J18KA147).

Conflicts of Interest:
The authors declare no conflict of interest.