Carotenoids Overproduction in Dunaliella Sp.: Transcriptional Changes and New Insights through Lycopene β Cyclase Regulation

Dunaliella is a green microalga known for its ability to produce high levels of carotenoids under well-defined growing conditions. Molecular responses to the simultaneous effect of increasing salinity, light intensity and decrease of nitrogen availability were investigated in terms of their effect on different metabolic pathways (isoprenoids synthesis, glycolysis, carbohydrate use, etc.) by following the transcriptional regulation of enolase (ENO), 1-deoxy-D-xylulose 5-phosphate synthase (DXS), lycopene β-cyclase (LCYB), carotene globule protein (CGP), chloroplast-localized heat shock protein (HSP70), and chloroplast ribulose phosphate-3-epimerase (RPE) genes. The intracellular production of carotenoid was increased five times in stressed Dunaliella cells compared to those grown in an unstressed condition. At transcriptional levels, ENO implicated in glycolysis, and revealing about polysaccharides degradation, showed a two-stage response during the first 72 h. Genes directly involved in β-carotene accumulation, namely, CGP and LCYB, revealed the most important increase by about 54 and 10 folds, respectively. In silico sequence analysis, along with 3D modeling studies, were performed to identify possible posttranslational modifications of CGP and LCYB proteins. Our results described, for the first time, their probable regulation by sumoylation covalent attachment as well as the presence of expressed SUMO (small ubiquitin-related modifier) protein in Dunaliella sp.


Introduction
The genus Dunaliella is composed of 24 unicellular species, uninucleate [1], enclosed by a glycocalyx [2] and without a rigid polysaccharide wall, allowing a rapid response to hypo or hyperosmotic conditions [3]. The halotolerant Dunaliella salina, the holotype of this genus, is characterized by an exceptional capacity to survive under some environmental stress by producing added value compounds as glycerol, polyunsaturated fatty acids, and a high percentage of β-carotene reaching 10% to 14% of dry weight [2].
The majority of carotenoids are C 40 isoprenoids with color ranging from yellow to red and are divided into carotenes and xanthophylls [4]. The β-carotene molecule is an orange compound, comprising a long conjugated chain of eight isoprene units with a center of symmetry and a β-cyclohexene ring at each end [5]. It can be found in four possible stereoisomers: 9-cis β-carotene,

Microalga and Growth Conditions
Dunaliella sp. was isolated from the Sebkha of Sidi El Hani (Sousse, Tunisia), having 97% of similarity with D. salina, D. quartolecta, and D. polymorpha [34]. It was cultivated in the F/2 medium [35] based on artificial seawater (ASW). Cells were grown either under unstressed condition (DSC) with a salinity (NaCl concentration) equal to 0.46 M (27 g·L −1 ), a light intensity of 80 µmol photons m −2 s −1 , and a Nitrogen source (0.88 mM NaNO 3 , 9.89 mM KNO 3 ) or stressed condition (DSS) with a salinity equal to 2 M (116.88 g·L −1 ), a light intensity of 540 µmol photons m −2 s −1 , and a Nitrogen source (1 mM NaNO 3 , 0 mM KNO 3 ). For this purpose, an exponential growing starting culture was used to inoculate both control and stressed culture at day 0. Cultures were maintained at 25 • C under continuous illumination by a fluorescent lamp.

Pigments Extraction
Pigments were extracted from 2 mL of cell cultures. To this pellet, 2 mL of 96% ethanol was added, and the mixture was sonicated at 40 kHz over 20 min and then incubated overnight at 4 • C in the dark. Cells debris was removed by centrifugation (10,000× g, 10 min) to recover pigments in the supernatant. Total chlorophyll and carotenoid amounts were measured by spectrophotometry using equations given by Wellburn and Lichtenhaler [36]. As stated in a previous study, we assumed that the total carotenoid content could be a good indicator of the β-carotene content as it is the main carotenoid in the pigment composition of D. salina [37].

Flow Cytometry Analysis
To study the intrinsic and extrinsic characteristics of the microalgae cell, flow cytometry analysis was used, as described by Dammak et al. [38]. All cytometric measurements were done using a BD LSR (Becton Dickinson, Franklin Lakes, NJ, USA) equipped with two lasers: argon (λex488 nm) and helium-neon (λex633 nm). Data were processed with CellQuest Pro software. The diagrams presented in this document were processed online (www.cytobank.org). Argon probe was used as it is associated with the emission channels FL1 (FITC), FL2 (PE), and FL3 (PerCP), which emit at the wavelengths 530/30, 575/26, and 682/13 nm, respectively. The signal intensities of fluorescence were calculated from the mean fluorescence intensity of 10,000 events.

RNA Extraction
Stressed and non-stressed cells at different stages of incubation (0, 4, 24, 48, and 72 h) were harvested by centrifugation at 1800× g for 10 min. RNA was extracted using TRIzol ® reagent (Invitrogen, Carlsbad, CA, USA) and further purified using column and wash buffer of the RNeasy Mini Kit (Qiagen, Hilden, Germany). RNA concentrations and quality were determined by the measurement of absorbance (260 nm/280 nm) using a NanoDrop 2000 Spectrophotometer (Thermo Scientific, Waltham, MA, USA) and agarose gel electrophoresis. Approximately 1 µg of the resulting total RNA was used for the synthesis of cDNA using the PrimeScript™ RT reagent Kit with gDNA Eraser (TaKaRa, Shiga, Japan) according to the manufacturer's protocol.

Primers Design and Validation, Reference Gene Selection
Specific primer pairs for alpha-tubulin (α-TUB), beta-tubulin (β-TUB), Enolase (ENO), 1-deoxy-D-xylulose 5-phosphate synthase (DXS), lycopene β-cyclase (LCYB), chloroplast-localized heat shock protein (HSP70), and chloroplast ribulose phosphate-3-epimerase (RPE) were designed based on the sequences of the corresponding genes in the D. salina genome deposited in GenBank. Whereas, primer pairs for carotene globule protein (CGP) amplification were designed based on the sequence of the corresponding gene in Dunaliella bardawil genome, as the sequence among Dunaliella salina is not available. All primer sequences were designed using Primer-BLAST from NCBI and Primer 3 Input 4.0. Primer pairs used for Actin (Act) and 18S rRNA are those cited by [39,40], respectively. Primer set and amplicon characteristics are represented in Table 1. The designed primers obeyed to the following criteria: a length between 18-24 nucleotides, an amplicon length from 100 to 200 bp, chosen temperature (Tm) around 60 • C, and avoiding secondary structures and self-and cross-annealing. The specificity of the designed primer was verified by the blast search function (http://www.ncbi.nlm.nih.gov). Furthermore, the optimal concentrations were set after testing three different primer concentrations ranging from 100 to 200 nM to avoid non-specific products and/or primer-dimer artifacts.
Gene coding for Act, 18S RNA, α-TUB, and β-TUB was used for the screening of the most stable gene in all experiment conditions as they are known to be classic plant reference gene candidates [41]. Quantification cycle (Cq) values were imported into the comprehensive web-based tool RefFinder (http://150.216.56.64/referencegene.php?type=reference), which integrates commonly used reference gene evaluation programs together, including geNorm, NormFinder, Bestkeeper, and the comparative delta Ct methods [41,42]. Results elaborated by RefFinder were confirmed by GeNorm, which is available on excel file [43].

qPCR Conditions
qPCR was performed using the Step One QuantiFast SYBR Green PCR Kit (Qiagen, Germany) according to the manufacturer's manual using the LightCycler 480 II (Roche, Basel, Switzerland). Cycling parameters were an initial denaturation at 95 • C for 5 min, 40 cycles at 95 • C for 10 s, and 60 • C for 30 s. A step of fusion by heating from 60 • C to 95 • C was added to verify the amplification specificity of each primer pair. Data acquisition and analysis were performed by LightCycler ® 480 Software, Version 1.5 (Roche). Relative gene expression was calculated using reference genes α-TUB. The Cq values were corrected with the estimated PCR efficiency, which was evaluated using a four-point standard curve consisting of a 2-fold dilution series. All reactions were performed in triplicate. The relative quantification of the stressed sample at different times versus non-stressed sample at T0 was calibrated, and the fold changes were then calculated using the following Equation (1) [44]: Calculated fold changes were analyzed statistically using XLSTAT 2017.5.47467 version (Addinsoft). One-way analysis of variance (ANOVA) followed by Tukey (honest significant difference) HSD test and Dunnett test were applied.

Retrieval of Protein Sequences and In Silico Analysis
The amino acid sequences of LCYB and CGP from Dunaliella sp. were retrieved from the protein database of the National Center for Biotechnology Information (NCBI, http://www.ncbi. nlm.nih.gov/protein/) under the accession numbers of ACA34344.1 and AFU62958.1, respectively. The nucleotide sequence of SUMO (small ubiquitin-related modifier) was searched by blast against the Expressed Sequence Tag (EST) database of Dunaliella sp. Then, the open reading frame (ORF) sequences were generated from the corresponding cDNA (GenBank: HO848823.1) using ORFfinder from NCBI, and the obtained proteins were blasted against the non-redundant protein database. The protein sequence corresponding to ORF2 presented 52.87% and 49.47% with SUMO protein from two Chlorophyta Ostreococcus tauri and Monoraphidium neglectum, respectively. The motifs search was performed using the motif search tools from KEGG (Kyoto Encyclopedia of Genes and Genomes) (https://www.genome.jp/tools/motif/). Prediction of posttranslational modifications was performed with SIB (Swiss Institute of Bioinformatics) bioinformatics tools from the ExPASy bioinformatics resource portal (www.expasy. org). Sumoylation was further assessed using SUMOplot (http://www.abgent.com/tool/sumoplot). Multiple sequence alignment of LCYB sequences was performed with Clustalw using default parameters [45] and formatted with Espript 3.0 [46]. TMHMM server (http://www.cbs.dtu.dk/services/ TMHMM) [47] was used to predict transmembrane helices.

Three-Dimensional Protein Modeling
Homology modeling was performed with the SWISS-MODEL server (http://www.expasy.org/ swissmod/). For the LCYB model, the crystal structure of geranylgeranyl reductase (PDB ID 4OPL) was used as a template since no structure of LCYB was found in the database [16]. The modeled structure was further validated by performing the Ramachandran map analysis using the PROCHECK server [48]. However, the generated model was truncated in its N-terminal and C-terminal ends due to the non-total coverage sequence with the template. For SUMO protein, the solution structure of SUMO from Trypanosoma brucei was used as a template (PDB ID 2k8h.1), and the generated model was equally validated.
The models were viewed using The Discovery Studio software [Dassault Systemes BIOVIA, Discovery Studio Modeling Environment, Release 4.5, 2015] and The PyMol Molecular Graphics System v1.2r3 pre (Schrödinger, LLC, New York, NY, USA). Possible active sites of LCYB were determined using the CASTp web server [49].

Flavin Adenine Dinucleotide (FAD)-Docking
The FAD docking was performed using the Galaxy web server, where parameters are set as default. GalaxyDock performs conformational space annealing (CSA) global optimization to find the optimal binding pose of a ligand both in the rigid-receptor mode and the flexible-receptor mode [50].

Influence of Stressed Condition on Cell Growth
Dunaliella sp. was grown in F/2 medium under unstressed (DSC) and stressed (DSS) conditions. In the latter case, a rise in salinity and light intensity was applied against a decrease in nitrogen source. The follow-up of the growth among 29 days showed that DSS maintained the same growth rate as DSC during the first 6 days of culture ( Figure 1a). Subsequently, DSC continued to grow during the follow-up while DSS appeared to stop dividing.

Pigments Accumulation
The production of pigments by DSC and DSS during one month of culture was followed by biochemical quantification. For DSC, the increase in pigment concentration over the time was due to the increase in cell numbers and not to an intracellular over-accumulation as we noticed that the concentrations presented in µg/mL increased over time (Figure 1b) whereas those presented in pg/cell were almost constant ( Figure 1c). Whereas for DSS, the metabolism of cells was totally oriented to carotenoid accumulation. Chlorophyll production measured in the culture (µg/mL) and on the intracellular side (µg/cell) were maintained at a basal level (Figure 1b,c). Carotenoids production in culture and intracellularly increased simultaneously (Figure 1b,c).
The accumulation of carotenoids by DSS was greater than that in DSC. For example, the number of carotenoids produced was 3.38 ± 0.33 and 17.96 ± 0.67 pg/cell in DSC and DSS, respectively, on day 17, which represents a five-fold increase. The ratio of carotenoids/chlorophylls in DSS culture vs. DSC culture was highly increased, reaching 18 times at the end of the exponential phase.
Dunaliella sp. was grown in F/2 medium under unstressed (DSC) and stressed (DSS) conditions. In the latter case, a rise in salinity and light intensity was applied against a decrease in nitrogen source. The follow-up of the growth among 29 days showed that DSS maintained the same growth rate as DSC during the first 6 days of culture ( Figure 1a). Subsequently, DSC continued to grow during the follow-up while DSS appeared to stop dividing.

Flow Cytometry Analysis
Cell populations were selected based on cell size FSC (Forward Scatter) and granulation SSC (Side Scatter) (Figure 2a). For DSC, there were two subpopulations with different granularity characteristics: a major sub-population (65%) with high SSC and a sub-population (35%), which was less grainy. While for DSS, a unique cell population was present resembling the major subpopulation of DSC. emission signal FL3 (682/13 nm) was more adequate to detect the fluorescence of chlorophyll a as its fluorescence emission spectrum is characterized by a major peak at 683/20 nm, attributable to the activity of PSII [51]. Stressed cells had a reduced fluorescence emission in FL3 corresponding to a regression of intracellular content of chlorophyll measured by biochemical methods. In contrast, an important increase of fluorescence emission was detected through FL1 and FL2 channels, which was likely due to the overproduction of carotenoids in stressed cells compared to unstressed ones [52,53].  The mean fluorescence intensities of the stressed and unstressed cells were also determined ( Figure 2b). Fluorescence was detected on the channel FL1, FL2, and FL3. The red fluorescence emission signal FL3 (682/13 nm) was more adequate to detect the fluorescence of chlorophyll a as its fluorescence emission spectrum is characterized by a major peak at 683/20 nm, attributable to the activity of PSII [51]. Stressed cells had a reduced fluorescence emission in FL3 corresponding to a regression of intracellular content of chlorophyll measured by biochemical methods. In contrast, an important increase of fluorescence emission was detected through FL1 and FL2 channels, which was likely due to the overproduction of carotenoids in stressed cells compared to unstressed ones [52,53].

Evaluation of Reference Gene Stability
The comprehensive ranking generated by RefFinder for the total experimental conditions revealed that alpha-tubulin (α-Tub) was the most stable gene followed by the actin gene (Act), while 18S RNA was the least stable gene (Supplementary Figure S1a). The analysis via GeNorm showed an expression stability value (M) less than 1.5 for all tested genes, which indicated that they were all accepted as reference gene candidates (Supplementary Figure S1b). Nonetheless, GeNorm classified α-Tub and Act as the best stable reference genes since they had lower M values compared to those of 18S RNA.

Gene Expression Profile under Stressed Conditions
A follow-up of the expression of some genes involved in the carotene accumulation process during the first 72 h of stress time was processed by quantitative PCR. The first gene followed was ENO involved in the glycolysis cycle and giving an idea about carbon origin (polysaccharide degradation, photosynthesis, etc.) and destiny (isoprenoids synthesis, carbohydrate synthesis, etc.) in the cell. Its transcript level decreased 0.48 times during the first 24 h of stress. Then, it tended to re-stabilize after 72 h by having a relative expression level close to 1 (Figure 3a) and without a significant difference in its expression with unstressed cells (Figure 3g  . Bands intensity correlates with values found by qPCR quantification. Linear regression between the relative expression of LCYB and CGP, showing a linear relationship between the expressions of these two genes (h). An asterisk indicates a significant difference between gene expressions at different stressed times, according to the one-way ANOVA test followed by Tukey test. * p < 0.05; ** p < 0.005; *** p < 0.001. Data are expressed as means ± SD (n = 3). Model, confidence interval (Mean 95%), confidence interval (obs. 95%).   . Bands intensity correlates with values found by qPCR quantification. Linear regression between the relative expression of LCYB and CGP, showing a linear relationship between the expressions of these two genes (h). An asterisk indicates a significant difference between gene expressions at different stressed times, according to the one-way ANOVA test followed by Tukey test. * p < 0.05; ** p < 0.005; *** p < 0.001. Data are expressed as means ± SD (n = 3).
confidence interval (Mean 95%),  . Bands intensity correlates with values found by qPCR quantification. Linear regression between the relative expression of LCYB and CGP, showing a linear relationship between the expressions of these two genes (h). An asterisk indicates a significant difference between gene expressions at different stressed times, according to the one-way ANOVA test followed by Tukey test. * p < 0.05; ** p < 0.005; *** p < 0.001. Data are expressed as means ± SD (n = 3).
The two genes LCYB and CGP, which are directly implicated in carotene accumulation, showed the most important rise. Data showed that their expression was increased since 4 h of stress time (Figure 3c,d) at an exponential pace. For the LCYB gene, its expression was found to be about 3-fold higher in stressed culture than the unstressed one during the first 72 h (Figure 3g). However, for CGP, the expression level in the stressed cell was very high; it was improved by almost 11-fold at 24 h (Figure 3d). There was a drastic increase in their expression between 24 h to 72 h of stress conditions: up to 10-fold for LCYB and to 54-fold for CGP (Figure 3c,d). We noticed also a positive correlation in their mRNA expression (R 2 = 0.997).
The less important increase was registered in HSP70B mRNA expression level, where the increase reached 3 and 2.27 times at 72 h in stressed and non-stressed conditions, respectively (Figure 3e,g). Chloroplast ribulose phosphate-3-epimerase (RPE) transcript level showed no significant variation during the first 72 h of stress (Figure 3f,g).

In Silico Analysis
The following in silico analysis was focused on the two genes showing the most important changes in their transcriptional levels. The motif search results showed that CGP proteins presented a SOUL heme-binding protein motif (PF04832) ranged from residue 283 to residue 350, while LCYB possessed FAD-binding motifs (PF00890), ranging from residue 123 to residue 154.
Prediction of posttranslational modifications (PTMs) for CGP showed four putative protein kinase C phosphorylation sites distributed in the N terminal part, two probable N-myristoylation sites, and two sumoylation sites (Table 2). On the other hand, LCYB was found to have eight probable protein kinase C phosphorylation sites mostly distributed in its N-terminal part, one N-glycosylation site (residue 51), six N-myristoylation sites, and one high probability sumoylation site, as well as two sumoylation interaction motifs ( Table 2).

Enzyme Modeling
An investigation of CGB, LCYB, and SUMO (small ubiquitin-related modifier) protein 3D structures was undertaken to understand or find possible molecular determinants related to their up-regulation. However, only LCYB and SUMO protein were successfully modeled, since no structural PDB files having an acceptable identity with CGP were found in the database. Thus the following part was essentially focused on LCYB molecular study and its probable interaction with SUMO protein in Dunaliella.

Alignment Analysis
The amino acid sequence of LCYB from D. salina (accession number: ACA34344.1) was compared with those of known β-lycopene cyclase homologous proteins (Figure 4). D. salina LCYB showed high similarity with many β-lycopene cyclases from other microalgae, e.g., 68% identity with that from Haematococcus lacustris (AAO64977.1), 65% with that from Chlamydomonas eustigma (GAX83195.1). In contrast, less identity was observed with plant LCYB (50% with that from Arabidopsis thaliana (AAA81880.1)) and very low identity with that from cyanobacteria (33% with LYBC from Acaryochloris marina (WP_012166204.1). The sequence alignment also revealed that D. salina LCYB exhibited N-terminal and C-terminal insertions when compared to plant and cyanobacteria LCYB in addition to an internal insertion of 15 residues corresponding to the α6 helix in the generated 3D model. Several conserved regions in the protein sequences already identified in other lycopene cyclases from bacteria and plants were also found in D. salina LCYB [18]. Amongst these conserved sequences, PTFLYAM and FLEETSL, highlighted in red in Figure 4, were involved in the formation of β12 and β13 strands. To further understand the role of these conserved regions in the enzyme activity, we proposed further investigation of the 3D model.

Overall 3D Model of D. salina LCYB
Because no structure has been reported in the LCYB family, the crystal structure of geranylgeranyl reductase from a Sulfolobus acidocaldarius strain, which displays a 27% sequence identity with D. salina LCYB, was used as a template, as described previously by [16].
The overall structure of the enzyme could be divided into two distinct domains: a catalytic domain involving a FAD binding site, as well as a substrate-binding site and a transmembrane domain ( Figure 5). The catalytic domain belongs to the Rossmann-fold family structure, consisting of a core height-stranded β-sheet flanked by four α-helices. In our case, the catalytic domain also comprised of residues Gly128-Ser129-Gly130-Pro131-Ser132-Gly133, known as a G-x-G-x-x-G motif, for nucleotide-binding [54].  Because no structure has been reported in the LCYB family, the crystal structure of geranylgeranyl reductase from a Sulfolobus acidocaldarius strain, which displays a 27% sequence identity with D. salina LCYB, was used as a template, as described previously by [16].  Calculations of the size and shape of binding sites performed using the CAST program ( Figure  6) demonstrated the presence of a small cavity for FAD-binding (radius probe was set to 2), with a total volume of 139.5 Å 2 , and a larger cavity for the substrate binding (radius probe was set to 1.4), with a total volume of 480.14 Å 2 and an interior length of about 35 Å, able to accommodate a large and hydrophobic substrate like lycopene.  The second domain consisted of only three large helices covering the total length of the enzyme, which are suspected to constituting the transmembrane part of the enzyme, as predicted by the TMHMM server.
Calculations of the size and shape of binding sites performed using the CAST program ( Figure 6) demonstrated the presence of a small cavity for FAD-binding (radius probe was set to 2), with a total volume of 139.5 Å 2 , and a larger cavity for the substrate binding (radius probe was set to 1.4), with a total volume of 480.14 Å 2 and an interior length of about 35 Å, able to accommodate a large and hydrophobic substrate like lycopene. Calculations of the size and shape of binding sites performed using the CAST program ( Figure  6) demonstrated the presence of a small cavity for FAD-binding (radius probe was set to 2), with a total volume of 139.5 Å 2 , and a larger cavity for the substrate binding (radius probe was set to 1.4), with a total volume of 480.14 Å 2 and an interior length of about 35 Å, able to accommodate a large and hydrophobic substrate like lycopene.

Discussion
Dunaliella is a green microalga, well-known by its high production of carotenoids under specific conditions. In this context, we characterized carotenoids production among an isolated strain of Dunaliella sp. cultured under two conditions: unstressed and stressed. The growth rate was influenced by the stress condition, causing a decrease in the DSS growth rate after a week. This was probably due to the depletion of the nitrogen source initially used in the medium. Similar results were also found by Sánchez-Estudillo et al. [56], where D. salina grown under nitrogen limitation maintained the same growth rate as those under nitrogen sufficiency during 8 days and then decreased to the half in about 20 days. The decrease in cell number was accompanied rather with an increase in intracellular carotenoids accumulation and a decrease of chlorophyll amount. Lv et al. in 2016 [21] studied the effect of nutrition depletion on carotenoid production that increased from approximately 2 pg/cell in complete medium to 9 pg/cell in nitrogen-depleted medium on day 15. This observation is in agreement with a previous study showing that the limitation of nitrogen availability accompanied by an increase of salinity and light intensity induced carotenoid accumulation [34]. The increase of salinity engenders also an increase in glycerol and polysaccharide synthesis accompanied by starch degradation [57][58][59].
For a better characterization of microalgae cells, they were subjected to flow cytometry analysis. Among DSC, we noticed the presence of two populations, which could be due to different evolution cells stage. The analysis of cells at the end of the exponential phase did not permit to display the increase in size and volume of DSS cells versus DSC. Davidi et al. [32] showed that during the first day of environmental stress, cells accumulated lipid globules, causing an increase of cell granularity. After two days, this granularity declined as the stress level would be reduced due to the increase in

Discussion
Dunaliella is a green microalga, well-known by its high production of carotenoids under specific conditions. In this context, we characterized carotenoids production among an isolated strain of Dunaliella sp. cultured under two conditions: unstressed and stressed. The growth rate was influenced by the stress condition, causing a decrease in the DSS growth rate after a week. This was probably due to the depletion of the nitrogen source initially used in the medium. Similar results were also found by Sánchez-Estudillo et al. [56], where D. salina grown under nitrogen limitation maintained the same growth rate as those under nitrogen sufficiency during 8 days and then decreased to the half in about 20 days. The decrease in cell number was accompanied rather with an increase in intracellular carotenoids accumulation and a decrease of chlorophyll amount. Lv et al. in 2016 [21] studied the effect of nutrition depletion on carotenoid production that increased from approximately 2 pg/cell in complete medium to 9 pg/cell in nitrogen-depleted medium on day 15. This observation is in agreement with a previous study showing that the limitation of nitrogen availability accompanied by an increase of salinity and light intensity induced carotenoid accumulation [34]. The increase of salinity engenders also an increase in glycerol and polysaccharide synthesis accompanied by starch degradation [57][58][59].
For a better characterization of microalgae cells, they were subjected to flow cytometry analysis. Among DSC, we noticed the presence of two populations, which could be due to different evolution cells stage. The analysis of cells at the end of the exponential phase did not permit to display the increase in size and volume of DSS cells versus DSC. Davidi et al. [32] showed that during the first day of environmental stress, cells accumulated lipid globules, causing an increase of cell granularity. After two days, this granularity declined as the stress level would be reduced due to the increase in β-carotene accumulation, and the cell's granularity would become almost like the unstressed one [60]. Based on autofluorescence proprieties of chlorophylls and carotenoids pigment, we confirmed the deviation of Dunaliella metabolism to carotenoids accumulation. In view of this change observed at the phenotypic level, we decided to analyze its response during the first three days of stress, at the transcriptional level of ENO, DXS, LCYB, CGP, HSP70, and RPE genes. For this purpose, the first step was the determination of the best reference gene, which appeared to be the α-Tub. Two other studies have also reported the stability of the α-Tub gene under severe conditions, one of which was conducted with Dunaliella sp. grown under high light stress [61]. Similarly, a second study realized with Syntrichia caninervis under abiotic stress; the α-Tub gene was also found to be stable amongst 15 reference genes, and 18S was the less unstable one [41].
Enolase is a ubiquitous enzyme that catalyzes the penultimate reversible step of glycolysis, which is the unique dehydration of D-2-phosphoglycerate (2-PGA) to phosphoenolpyruvate (PEP) [62]. Its regulation is a key step in the orientation of the organism's metabolism to glycerol or tricarboxylic acid formation [63]. In our study, we found that the transcript level of ENO declined during the first 24 h of stress, which was in accordance with results founds by Xia et al., reporting a decrease in mRNA expression and ENO activity during 12 h of hyperosmolarity stress [59]. The expression of ENO tended to increase and restabilize at 72 h of culture with a relative expression level close to 1 ( Figure 3a) and similar to that of DSC (Figure 3h). The same profile was found by Ruan et al., where under hypersaline stress (3 M), the expression level of ENO protein and its activity decreased after 3 h of stress among D. salina and re-increased at 6 h [62]. The authors explained that the down-regulation of ENO inhibited the transformation of 2-PGA to pyruvate, which decreased the amount of available carbon for the tricarboxylic acid cycle and the formation of PEP. Dunaliella is indeed capable of surviving under hyperosmotic stress, thanks to the rapid accumulation of osmoprotectants like glycerol. In this condition, starch reserved in the chloroplast will be degraded into glucose to be converted to glycerol [2]. As a result, 2-PGA and then DHAP (dihydroxyacetone phosphate) will be accumulated in the middle of the glycolysis pathway. The DHAP is then transformed via GPDH (glycerol-3-phosphate dehydrogenase) and GPP (glycerol-3-phosphate phosphatase) to glycerol. The restabilization of ENO gene expression after 72 h allowed a sufficient supply of carbon to PEP synthesis and, subsequently, of the pyruvate, essential for the synthesis of isoprenoids. This hypothesis was in accordance with that of Ramos et al. [2], supposing that response to salt stress goes through two stages where, in the first step, it accumulates glycerol and glycine betaine, and in a later time, it produces neutral lipid and carotenoid.
The pyruvate and G3P (glyceraldehyde-3-phosphate) will be converted to DXP (deoxyxylulose 5-phosphate) by DXS, catalyzing the first reaction in MEP (methylerythritol 4-phosphate) pathway [1]. Therefore, it is implicated in the synthesis of isoprenoids, such as carotenoids, chlorophylls, volatiles [64], and squalene [65]. The regulation of the MEP pathway required more clarification, resulting in the disagreement for the regulation of DXS. Lv et al. [21] found that the transcript level of DXS increased under nitrogen deprivation to provide more substrate to carotene synthesis. While Sánchez-Estudillo et al. [56] reported that the mRNA level of DXS decreased under nitrogen limitation on account of growth decrease, low photosynthetic activity, and chloroplast development. Our results showed a DXS increase at 24 h, but it trended to stabilize after 72 h with lower expression level than non-stressed cells. The discrepancy in result between publications suggested the possibility of the existence of more than one isoform of DXS having different expression profiles and functions like among the microalgae, Botryococcus braunii, and plants (tomato, Arabidopsis, rice, maize, citrus fruit) [64,66]. Besides, the expression of DXS is subject to translational, posttranslational, and posttranscriptional modulations [67,68]. In fact, it was demonstrated that DXS expression was regulated by feedback exercised by isopentenyl diphosphate (IDP) and Dimethylallyl diphosphate (DMADP) [67,69]. It is important to recall the assumption that D. salina lost the mevalonate (MVA) pathway and uses only the MEP pathway to provide isoprenoids precursors [2]. Also, Sánchez-Estudillo and his collaborators supposed that carotenogenesis, dependent on the regulation of enzymes, is located in the final step of the MEP pathway [56]. This supposition was confirmed by the significant rise at 72 h of the expression level of LCYB, which is involved in the cyclization of linear lycopene to β-carotene. Thus, the accumulation of this isoprene is directly related to the activity of this enzyme and the molecular regulation and control of its gene [20]. Our result was in accordance with those of Ramos et al. [18]. In D. bardawil, LCYB was regulated in salt stress, thanks to two salt regulated elements (SRE), having GT1GMSCAM4 sequence in the promoter and a GT-rich region in its first intron [70]. The promotor sequence of LCYB also contains a light-inducible element (GATABOX, CGCGBOXAT, SORLIP1AT) and W-boxes, a binding domain for the WRKY transcriptional factor [71]. WRKY transcription factor family is of high interest since it is involved in diverse biotic/abiotic stress response and developmental/physiological processes [70]. Moreover, our data showed that the increase in LCYB was concomitant with a rise in a lipid-associated protein CGP. It has a crucial role in maintaining the stability of βC-plastoglobules as it is localized at the periphery of the chloroplast [31][32][33], preventing the coalescence of the β-carotene globule, thanks to its hydrophilic layer, which covers the hydrophobic pigment core [33]. The gene coding for CGP is present in four different copies in D. salina 'Teodoresco' [32]. High light and/or nitrogen deficiency induces the accumulation of CGP in D. bardawil in parallel with β-carotene [33]. In accordance with our results, an increase in CGP mRNA expression was noticed during the first hours of nitrogen deprivation of D. bardawil [32]. This increase was associated with an accumulation of triacylglycerol (TAG) and β-carotene [72]. D. bardawil mobilizes TAG from the cytoplasm into the chloroplasts to form more plastoglobules [32]. In the non-carotenogenesis condition, the absence of these lipids globules was noticed.
Heat shock protein (HSP70B) is a 70 kDa protein. It is a member of ATP-dependent molecular chaperones families, highly conserved and encoded by a single gene in green algae [73]. The increase in HSP70 expression lasting for 72 h of stress time can be explained by a cell requirement for chaperone molecule. A similar result was obtained by Yokthongwattana et al., 2001 [74], where the level of HSP70 gene expression increased immediately under irradiance stress by shifting cells from 75 to 2500 µmol photon. m −2 ·s −1 among D. salina. Authors revealed that cells lost D1 protein by photodamage; thus, it would be replaced by the accumulation of HSP70B, which could afford conformational protection to the dissembled PSII core complex. In fact, protein misfolding is a process under oxidative and environmental stress, leading to an increase in molecular chaperone and protease expression [75].
The RPE catalyzes the reversible epimerization reaction of D-ribulose-5-phosphate into D-xylulose-5-phosphate in the Calvin cycle. During these first three days of stress, the level of its mRNA was significantly less expressed in stressed condition than in unstressed one after 72 h of stress. Similarly, but at the proteomic level, the amount of this enzyme was also decreased among Dunaliella under stress [34] and among Nannochloropsis oceanica under nitrogen deprivation [76]. This means that the photosynthetic efficiency is affected by nitrogen amount [76], but, at the same time, the down-regulation of this enzyme allows the accumulation of xylulose-5P, which probably would serve as a substrate for G3P production in the pentose phosphate cycle [34].
In order to identify if these transcriptional changes go along with posttranslational ones, bioinformatics tools were used to analyze protein sequences of the two genes (CGP and LCYB) having the most increased transcriptional levels. Prediction of PTMs showed numerous possible phosphorylation, palmitoylation, and myristoylation sites, as well as sumoylation sites. Thes results corroborated the fine regulation of the two enzymes after their fold and the major regulatory role they have in modulating the ratio of carotenoid biosynthesis. Unlike plants, data concerning posttranslational regulation modalities of CGP and LCYB were not well studied for microalgae. In addition, the presence of SUMO protein in Dunaliella was not reported previously.
Posttranslational modification by SUMO is a major regulating mechanism for eukaryotic protein functions. SUMO becomes commonly covalently attached to specific targets (lysine residue) through a succession of molecular "handoffs" relating multienzyme cascades consisting of an E1 activating enzyme, an E2 conjugating enzyme, and an E3 ligase. Attachment of SUMO is known to control target functions, such as protein-protein interactions, protein-DNA interactions, and subcellular localization.
Sumo regulates equally many important processes, such as signaling, transcription, DNA repair, and other stress responses. In plants, it was demonstrated that sumoylation plays an important role in reacting to stress and pathogens [77]. For microalga, three SUMO homologs and SUMO-related proteins were identified in C reinhardii and were established to be essential for the survival of the strain under stress conditions [78]. Concerning the regulation of LCYB by sumoylation, recent research reported the detection of a recognition sequence for sumoylation in LCYB from Bixa orellana that could modulate enzyme activity or localization [79]. To the best of our knowledge, this is the first report for a possible regulation by the sumoylation of LCYB in Dunaliella sp., and to further investigate the attachment modalities between the two proteins, we tried to carry out a 3D modeling analysis. For this purpose, a crystal structure of geranylgeranyl reductase from a Sulfolobus acidocaldarius strain was used since it displays an acceptable sequence identity/coverage, and it was used in the previous reports of LCYB models. The overall structure was divided into two distinct domains: a catalytic domain involving a FAD-binding site, as well as a substrate-binding site and a transmembrane domain. In fact, several studies have shown that lycopene cyclases use FAD as a cofactor [80] and that many of the enzymes involved in lycopene and carotene transformations are membrane-bound [81]. We found that the binding site was not far from the membrane helices and that the FAD-binding domain existed in the diagonally opposite side of the enzyme (Figure 6a). The inspection of the putative SUMO interacting residue (K271), as well as sumoylation interaction motifs, showed that they were both located on the same side and were very close to each other and the FAD-binding site (Figure 7b). This observation could be interesting in understanding the regulation mechanism exerted by SUMO on LCYB, and the overall work is expected to help the future investigation of transcriptional/translational regulation studies in Dunaliella sp.

Conclusions
In this study, Dunaliella sp. was submitted to stress conditions by varying simultaneously the salinity, the light intensity, and the nitrogen availability. Mentioned conditions result in the deviation of the metabolism for the carotenoid accumulation. The latter could substitute chlorophyll with the advantage of preserving cell viability until roughly the end of the exponential phase. Flow cytometry analyses revealed that stressed and unstressed cells exhibited differences in their intrinsic and extrinsic characteristics. At the transcriptional level, the expression of genes, closely related to the synthesis and the accumulation of β-carotene, increased significantly. On the other hand, the inspection of the three-dimensional model of lycopene-β-cyclase, a key enzyme in the carotenoid synthesis pathway, highlighted possible further regulation by posttranslational modification, especially sumoylation.
Supplementary Materials: The following are available online at http://www.mdpi.com/2076-3417/9/24/5389/s1, Figure S1: Geomean of ranking values of the candidate reference genes based on the RefFinder algorithm (a); Average expression stability M values of remaining control genes based on GeNorm algorithm (b). The mean Cq of three independent replicate experiments was used to generate Data. Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.