Lipid accumulation of Chlorella sp. TLD6B from the Taklimakan Desert under salt stress

Algal strain and culture conditions

The experimental algae were collected from the Taklimakan Desert (37°36′N, 80°23′E, 1,254 m altitude, 3.5 °C average surface temperature, 0.4 M average NaCl salt concentration, 71.2 W m⁻² net radiation at 16:00) in Xinjiang Province, China. Field experiments were approved by the College of Life Sciences of Shihezi University (project number: 2011-GB-200900-G). The algae were isolated in the laboratory from a mixed culture of Taklimakan Desert soil samples and were identified as Chlorella sp. TLD6B (Wang et al., 2014). They were maintained in an autotrophic culture system at the College of Life Sciences, Shihezi University, China. Bold basal medium (BBM; Nichols, 1973) was autoclaved, and algal cells in the logarithmic growth phase were inoculated into Erlenmeyer flasks containing 500 mL BBM medium at a ratio of 20% (V/V). The initial cell density was controlled at an OD₆₈₀ of approximately 0.1, and the preparation of BBM medium is described in Table S1. Light was provided by white fluorescent lamps at an intensity of 72 µmol photons m⁻² s⁻¹ with a 12h/12 h light-dark cycle, and the temperature was 23 °C. Flasks were shaken three times a day at regular intervals and cultured for 24 d (Mu et al., 2016).

Salt stress treatments

Algal cells in the logarithmic growth phase were inoculated into Erlenmeyer flasks containing 300 mL BBM medium, and 0.0, 0.1, 0.2, 0.4, 0.6, or 0.8 M NaCl was added to supplement the medium to a final volume of 500 mL. The 0.0 NaCl treatment was the control, and there were three replicates of each treatment. Algal culture was performed as described above. Samples were obtained every 6 d after inoculation to measure physiological and biochemical indexes such as biomass, carbohydrate concentration, and total lipid concentration.

Measurement of physiological and biochemical indexes

Culture solution (100 mL) was added to 10-mL centrifuge tubes at 4 °C, and the tubes were centrifuged at 8,000 rpm for 15 min. The supernatant was discarded, and two mL of deionized water was added. The cell sediments were washed three times by vortexing (5–10 s each time), then divided into two equal parts. One part was used for dry weight measurement, and the other was immediately frozen in liquid nitrogen and stored at −80 °C for further analysis.

Optical density at 680 nm (OD₆₈₀) was determined using a UV-5800H UV–VIS spectrophotometer (Metash, Shanghai, China), and a growth curve was drawn (Gong et al., 2013). The centrifuged Chlorella sp. TLD6B sediment was dried at 80 °C for 24 h to measure the dry weight (Wan et al., 2015). The specific growth rate was calculated using the method of Wang et al. (2016). (1)${\displaystyle \mu =\frac{lg{\mathrm{N}}_1-lg{\mathrm{N}}_0}{{\mathrm{t}}_1-{\mathrm{t}}_0}\times 3.322}$

where µis the specific growth rate, mg L⁻¹ d⁻¹; N₀ is the biomass from the previous sampling date, mg L^-1; N₁ is the current biomass, mg L^-1; and t₀ and t₁ are the culture times, t₁ −t₀ = 6 d.

The phenol-sulfuric acid method was used to measure carbohydrate content (Dubois et al., 1956). The absorbance was measured at 490 nm on an iMark Microplate Reader (BIO-RAD, USA). Different concentrations of glucose were used to create a standard curve (y = 8.55x −0.07, where x is the standard concentration (µg mL ⁻¹) and y is the absorbance value), from which the glucose concentration X_C of the sample could be calculated. X_C (µg mL⁻¹) was defined as shown in equation Eq. (2), and carbohydrate content (CC, %) was defined as shown in Eq. (3). (2)${\displaystyle {\mathrm{X}}_{\mathrm{C}}=\frac{\left(△\mathrm{A}+0.07\right)}{8.55}}$ (3)${\displaystyle \mathrm{CC}=\frac{{\mathrm{X}}_{\mathrm{C}}\times \mathrm{V}\times 100\text{\%}}{{\mathrm{V}}_{\mathrm{C}}\times \mathrm{W}}}$ where △A is the absorbance at 490 nm, V is the total volume of the sample (mL), V_C is the volume of the subsample used for measurement (µL), and W is the dry weight of the sample (g L⁻¹).

The lipid content was measured by the colorimetric sulfo-phospho-vanillin method as described by Mishra et al. (2014) with some modifications (Li et al., 2012). The absorbance was measured at 528 nm, and the corresponding linoleic acid concentration (µg mL⁻¹) was calculated from the standard curve (y = 0.0005x + 0.008). X_L (µg mL⁻¹) was defined as shown in Eq. (4), and the lipid content (LC, %) was defined as shown in Eq. (5). (4)${\displaystyle {\mathrm{X}}_{\mathrm{L}}=\frac{\left(△\mathrm{A}-0.008\right)}{0.0005}}$ (5)${\displaystyle \mathrm{LC}=\frac{{\mathrm{X}}_{\mathrm{L}}\times \mathrm{V}\times 100\text{\%}}{{\mathrm{V}}_{\mathrm{L}}\times \mathrm{W}}}$ where △A is the absorbance at 528 nm, V is the total volume of the extract (mL), V_L is the volume of the subsample used for measurement (μL), and W is the dry weight of the sample (g L⁻¹).

RNA extraction and sequencing

Based on the results of physiological measurements under salt stress, we selected Chlorella sp. TLD6B cells grown under low salt stress (0.1 M NaCl addition) and high salt stress (0.8 M NaCl addition) for 18 d for use in transcriptome sequencing and analysis. Samples collected from the control group, the 0.1 M NaCl treatment group, and the 0.8 M NaCl treatment group were labeled CK18, Nacl1, and Nacl2, respectively. Each treatment was replicated two times. Samples were centrifuged at 8000 rpm for 8 min at 4 °C. After removal of the supernatants, Chlorella sp. TLD6B cells were washed with distilled water three times and stored at −80 °C for RNA extraction.

Total RNA was isolated from all samples using the TRIzol reagent (Invitrogen, USA) according to the manufacturer’s instructions. Total RNA was sent to Beijing Novogene Technology Co., Ltd. for RNA quality assessment, library construction, and sequencing (150-bp paired-end reads) on the Illumina NovaSeq 6000 platform. Raw sequence reads have been deposited into the NCBI GEO under accession number GSE162916.

Transcriptome analysis

The total raw reads were processed using FASTQC software with default parameters. High-quality clean reads were obtained by removing reads that contained adapters, reads that contained more than 10% unknown nucleotides (N), and low-quality reads that contained more than 50% low-quality bases (Q_phred ≤ 20). At the same time, the GC content, Q20, Q30, and sequence repetition level of the clean reads were calculated. Transcriptome de novo assembly was performed using Trinity v2.4.0 software (Grabherr et al., 2011) with the parameter Kmer = 25. Corset v1.05 (Davidson & Oshlack, 2014) was used with default parameters to cluster the transcript sequences in order to remove redundant sequences and obtain unigene sequence sets for subsequent analysis. The longest transcript of each subcomponent was used as the unigene for functional annotation. All assembled unigenes were finally annotated by comparing their sequences against the NR (NCBI non-redundant protein sequences, e-value = 10⁻⁵), NT (NCBI nucleotide sequences, e-value = 10⁻⁵), SwissProt (e-value = 10⁻⁵), GO (Gene Ontology, e-value = 10⁻⁶), KEGG (Kyoto Encyclopedia of Genes and Genomes, e-value = 10⁻¹⁰), PFAM (protein family, e-value =0.01), and KOG/COG (euKaryotic Ortholog Groups/ Clusters of Orthologous Groups of Proteins, e-value = 10⁻³) databases. For NR, SwissProt, and KOG/COG annotations, Diamond v0.8.22 software was used. For NT annotation, Blastn (NCBI blast 2.2.28+) was used. PFAM protein family alignments were performed using the HMMER 3.0 package. GO annotations were obtained based on the annotation results from NR and PFAM using Blast2GO v2.5 software (Götz et al., 2008). The KEGG Automatic Annotation Server of KAAS was used for KEGG pathway analysis. The clean reads were aligned to the Trinity transcripts using Bowtie 2 with default parameters (Langmead & Salzberg, 2012). Following alignment, raw read counts for each transcript were derived using RSEM v1.2.15 (Li & Dewey, 2011) with default parameters and then normalized to FPKM (Fragments Per Kilobase of transcript per Million mapped reads) (Trapnell et al., 2010). The read counts were used as input for the DESeq R package (1.10.1) (Anders & Huber, 2010) to identify differentially expressed genes (DEGs) using the thresholds adjusted P value (padj) <0.05 and —log₂ FoldChange— ≥ 1. DEGs were mapped to terms in the GO and KEGG databases for functional and pathway analysis. The GOseq R package (Young et al., 2010) was used for GO enrichment analysis, and GO terms with corrected P value <0.05 were defined as significantly enriched in the DEG set. Significantly enriched KEGG pathways were identified in the DEGs using KOBAS v2.0.12 with a corrected P value <0.05 (Mao et al., 2005). The GO and KEGG enrichment maps were generated using the Beijing Novogene Technology Co., Ltd. server. The two figures were merged using Adobe Photoshop CS5 (Adobe Inc.; San Jose, CA, USA).

Verification of RNA-seq gene expression by quantitative real-time PCR

Samples of total RNA (3 µg) were taken from the same total RNA used for RNA sequencing and used for the synthesis of first-strand cDNA with the M-MLV reverse transcription kit (Takara, Japan) according to the manufacturer’s instructions. The cDNA obtained by reverse transcription was diluted 10-fold for use. Primer 5 software was used to design gene-specific real-time quantitative PCR (qRT–PCR) primers. qRT–PCR was performed using the Roche Light Cycler 480 system (Roche, Rotkreuz, Switzerland) and the SYBR Premix ExTaq kit (Takara, Japan) with a 10-µL PCR reaction system. The thermal cycle procedure was as follows: pre-incubation at 95 ° C for 5 min, followed by 40 cycles of 95 ° C for 15 s, 60 ° C for 15 s, and 72 ° C for 20 s. The 2^−ΔΔCt method was used to calculate relative gene expression using the GAPDH gene as a reference (Livak & Schmittgen, 2001). Pearson correlations between RNA-seq and qRT–PCR expression fold changes were calculated using SPSS 20.0. Treatment differences were assessed by analysis of variance using SPSS 20.0 (IBM Inc., Chicago, USA).

Data analysis

All physiological data are presented as the mean ± standard deviation of three replicates. SPSS 20.0 was used to perform one-way analysis of variance (ANOVA) and correlation analysis, and Tukey’s test was used to identify differences among individual treatments (P < 0.05). Correlations between biomass and physiological indexes after different durations of salt stress were obtained using the bivariate module in the correlate analysis tool of SPSS 20.0. Graphs were constructed using SigmaPlot 12.5 software.

Effects of different salt stress levels on the physiological characteristics of Chlorella sp. TLD6B

Biomass increased under low salt stress and decreased under high salt stress, and Chlorella sp. TLD6B entered the logarithmic growth phase on the sixth day of culture. The OD₆₈₀ measurement increased rapidly before day 18, but slowed by day 24, entering a stable phase. The OD₆₈₀ measurement gradually decreased as the concentration of NaCl added to the medium increased from 0.1 M to 0.8 M (Fig. 1A and Table S2). The algae were able to tolerate high levels of salt (0.6 M and 0.8 M NaCl addition), but their biomass was very low (Fig. 1B and Table S2). Biomass was higher in the 0.1 M and 0.2 M NaCl addition treatments than in the control, but the OD₆₈₀ measurement decreased with increasing NaCl concentrations in the remaining treatments and also decreased with time. In general, compared with the control, the 0.1 M and 0.2 M NaCl addition treatments significantly inhibited the OD₆₈₀ but increased biomass. Biomass increased markedly after 12 days of NaCl stress in the 0.2 M NaCl addition treatment, whereas addition of 0.4–0.8 M NaCl significantly inhibited both OD₆₈₀ and biomass.

Carbohydrate content decreased as NaCl concentration increased: with the exception of the 0.2 M NaCl addition treatment, the higher the NaCl concentration in the culture solution, the lower the algal carbohydrate content (Fig. 1C and Table S2). However, lipid content and biomass were improved by the addition of up to 0.2 M NaCl. Lipid content as a function of NaCl concentration was directly proportional to algal biomass, and Chlorella sp. TLD6B exhibited the highest lipid contents of 29.95% and 31.27% after the addition of 0.1 M and 0.2 M NaCl, respectively (Figs. 1B and 1D and Table S2). When NaCl concentrations greater than 0.2 M were added, lipid content decreased significantly as salinity increased (Fig. 1D and Table S2).

Carbohydrate content gradually increased early in cultivation (days 1–12) but decreased rapidly as total lipid content increased after day 12. At day 24, carbohydrate contents were 0.74-, 0.88-, 0.62-, 0.53-, and 0.38-fold that of the control in the 0.1, 0.2, 0.4, 0.6, and 0.8 M NaCl addition treatments, respectively (Fig. 1C and Table S2). By contrast, the lipid contents of the same treatments were 1.37-, 1.43-, 1.20-, 1.15-, and 1.13-fold that of the control (Fig. 1D and Table S2). At the same time, the stability of the specific growth rate of cells under low salt stress (0.1 and 0.2 M NaCl addition) was higher than that under medium and high salt stress (0.4–0.8 M NaCl addition), resulting in low lipid productivity under high salt stress (Figs. 1E and 1F and Table S2). This is consistent with the change in lipid content.

The biomass of Chlorella sp. TLD6B was significantly correlated with its physiological parameters on day 18 of salt stress. Cell density (OD₆₈₀) and carbohydrate content were significantly correlated with biomass on day 18, but the correlation between biomass and total lipid content was highest (r =0.79) (Table 1). Therefore, Chlorella cells under low salt stress (0.1 M NaCl addition) and high salt stress (0.8 M NaCl addition) on day 18 were selected for transcriptome sequencing to further study the mechanisms of lipid accumulation under salt stress.

RNA-sequencing, de novo transcriptome assembly, and functional annotation

RNA-seq was performed on the mRNA extracted from the two NaCl treatments and the control, generating 963,078,184 raw reads. A total of 947,225,244 clean reads were obtained after filtering on base quality score and read length. The GC percentage of the clean reads was nearly 66.0%, and the Q20 was greater than 96% (Table S3). Trinity was used to create a de novo transcriptome assembly from the high-quality clean reads, producing 219,577 transcripts with an average length of 1,394 bp. 155,503 non-redundant unigenes were assembled; their length ranged from 200 to 23,825 bp with an average length of 1,842 bp (Tables S4 and S5).

All assembled unigenes were functionally annotated by searching against multiple public databases. Of the 155,503 unigenes, 111,238, 60,439, 84,173, 105,026, 109,098, and 51,484 were successfully annotated in the NR, NT, SwissProt, PFAM, GO, and KOG databases, respectively. Unigenes were also functionally annotated by searching against the KEGG database. A total of 47,032 (30.24%) unigenes had an annotated function in the KEGG database, and 131,623 (84.64%) unigenes were successfully annotated in at least one database (Table 2).

DEGs in response to salt stress

There were a total of 16,402 DEGs across the low salinity (0.1 M NaCl addition, Nacl1) and high salinity (0.8 M NaCl addition, Nacl2) treatments, based on thresholds of P < 0.05 and —log₂FC— ≥ 1 (Fig. 2A and Table S6). Of these, 1,326 were up-regulated and 2,304 were down-regulated under low salinity, and 9,776 were up-regulated and 5,798 were down-regulated under high salinity. Eight hundred twenty-eight unigenes were differentially expressed only under low salinity, 12,772 were differentially expressed only under high salinity, and 2,802 were differentially expressed under both conditions (Fig. 2B and Table S6). Of the 16,402 total DEGs, 82.9% were differentially expressed under only the low or the high salinity treatment.

qRT–PCR verification of RNA-seq data

Nine DEGs in Nacl1 vs. CK18 (the 0.1 M NaCl treatment compared to the control on day 18) and Nacl2 vs. CK18 (the 0.8 M NaCl treatment compared to the control on day 18) were selected for qRT–PCR verification. These unigenes were mainly involved in fatty acid synthesis, TAG synthesis, the glycolysis/gluconeogenesis pathway, the pentose phosphate pathway, and energy metabolism of the mitochondrial respiratory chain. Seven were up-regulated and two were down-regulated. The qRT–PCR primers listed in Table S7 were used to verify the RNA-seq data (Table S7), and the results are shown in Fig. 3A and Table S7. The expression patterns of nine unigenes were similar in the qRT–PCR and RNA-seq data, although there was some deviation in fold change values between the two methods. The qRT–PCR results were highly correlated with the RNA-seq results (r =0.890, r² =0.791), indicating that the RNA-seq data were generally accurate and reliable (Fig. 3B).

GO and KEGG pathway enrichment analysis of the DEGs

To further explore their functions, the DEGs from the two salt treatments were annotated with GO terms and KEGG pathways for functional and pathway enrichment analyses (Table S6). Among the top 30 most highly enriched GO terms in the 828 specific DEGs from the low salinity treatment, one molecular function (MF) GO term was significantly enriched (corrected P value <0.05): transferase activity, transferring one-carbon groups (GO:0016741). Among the top 30 most highly enriched GO terms in the 2,802 common DEGs, 10 biological process (BP) terms and 20 MF terms were significantly enriched (corrected P value <0.05). In the BP category, the three most significantly enriched GO terms were G protein-coupled receptor signaling pathway (GO:0007186), defense response to other organism (GO:0098542), and response to external biotic stimulus (GO:0043207). The three most significantly enriched MF GO terms were peptidase activity (GO:0008233), peptidase activity, acting on L-amino acid peptides (GO:0070011), and endopeptidase activity (GO:0004175). The most strongly and significantly enriched BP GO term in the 12,772 specific DEGs from the high salinity treatment was transcription initiation from RNA polymerase II promoter (GO:0006367). The three most significantly enriched cellular component (CC) GO terms were RNA polymerase II transcription factor complex (GO:0090575), nuclear transcription factor complex (GO:0044798), and transcription factor TFIIA complex (GO:0005672). The most significantly enriched MF GO term was secondary active transmembrane transporter activity (GO:0015291) (Fig. 4 and Table S8). There were no significantly enriched KEGG pathways in the low salinity-specific DEGs. Among the top 20 KEGG pathways enriched in the common DEGs, seven were significantly enriched (corrected P value <0.05): plant hormone signal transduction (ko04075), arginine and proline metabolism (ko00330), glutathione metabolism (ko00480), betalain biosynthesis (ko00965), sulfur metabolism (ko00920), cysteine and methionine metabolism (ko00270), and nitrogen metabolism (ko00910). The 12,772 high salinity-specific DEGs were significantly enriched in mismatch repair (ko03430), DNA replication (ko03030), nucleotide excision repair (ko03420), proteasome (ko03050), photosynthesis (ko00195), porphyrin and chlorophyll metabolism (ko00860), and galactose metabolism (ko00052) (Fig. 5 and Table S9). Consistent with the KEGG analysis, GO enrichment analysis showed that common DEGs under low and high salinity were involved primarily in stress response, defense response, and membrane structure.

Transcriptional expression of unigenes involved in fatty acid and TAG biosynthesis pathways

Fatty acid biosynthesis

Acetyl-CoA, the precursor for fatty acid synthesis in the cytosol, derives from pyruvate (PYR) produced in glycolysis by pyruvate kinase (PK). Genes encoding two rate-limiting enzymes in glycolysis, hexokinase (HXK) and PK, were up-regulated under both the high and low salinity treatments, and 6-phosphofructokinase (PFK) was up-regulated under the high salinity treatment (Fig. 6 and Table S10). These changes would have promoted the accumulation of acetyl-CoA and energy.

The expression of key genes of the de novo fatty acid synthesis pathway also changed in response to salinity treatment. For example, the acetyl CoA carboxylase (ACCase) genes BC and BCCP were up-regulated on day 18 under both salinity treatments. The α-CT gene (Cluster-31803.29741) was up-regulated 2.7-fold under high salinity (Fig. 7A and Table S11), and the malonyl-CoA:ACP transacylase (MCTA) gene (Cluster-31803.71595) was also up-regulated only under high salinity (Fig. 7B and Table S11). Two 3-oxoacyl-[acyl carrier protein] synthase II (KAS II) genes (Cluster-31803.17167 and Cluster-31803.96146) were up-regulated under both salinity treatments, but these differences were not significant (Fig. 7C and Table S11). One ω-3 fatty acid desaturase (FAD) gene (Cluster-31803.83020) was down-regulated 4.5-fold under high salinity, and the other FAD gene (Cluster-31803.96909) was up-regulated 1.5-fold under low salinity (Fig. 7D and Table S11). A stearoyl ACP desaturase (SAD) gene (Cluster-31803.15683) was down-regulated under low salinity, whereas three SAD genes were up-regulated under high salinity (Cluster-31803.73844, Cluster-31803.73842, and Cluster-31803.73850) (Fig. 7E and Table S11).

Expression of key unigenes related to lipid synthesis pathways

To further explore mechanisms of lipid accumulation in Chlorella sp. TLD6B under salt stress, we compared the expression of genes associated with the metabolic pathways of starch-lipid biosynthesis. Under high salinity, a single gene encoding the starch-degrading enzyme α-amylase (AMY3, Cluster-31803.62982) was up-regulated 4.6-fold (Fig. 6 and Table S10).

Furthermore, genes encoding ACL, PDC/PDH, ME, G6PD, and PEPC were primarily up-regulated. All of these genes showed greater expression under high salinity than under low salinity (Fig. 6 and Table S10). The upregulation of these genes whose products function in central carbon metabolism may directly or indirectly promote the accumulation of fatty acids.