Development of a Chlamydomonas reinhardtii metabolic network dynamic model to describe distinct phenotypes occurring at different CO2 levels

The increase in atmospheric CO2 due to anthropogenic activities is generating climate change, which has resulted in a subsequent rise in global temperatures with severe environmental impacts. Biological mitigation has been considered as an alternative for environmental remediation and reduction of greenhouse gases in the atmosphere. In fact, the use of easily adapted photosynthetic organisms able to fix CO2 with low-cost operation is revealing its high potential for industry. Among those organism, the algae Chlamydomonas reinhardtii have gain special attention as a model organism for studying CO2 fixation, biomass accumulation and bioenergy production upon exposure to several environmental conditions. In the present study, we studied the Chlamydomonas response to different CO2 levels by comparing metabolomics and transcriptomics data with the predicted results from our new-improved genomic-scale metabolic model. For this, we used in silico methods at steady dynamic state varying the levels of CO2. Our main goal was to improve our capacity for predicting metabolic routes involved in biomass accumulation. The improved genomic-scale metabolic model presented in this study was shown to be phenotypically accurate, predictive, and a significant improvement over previously reported models. Our model consists of 3726 reactions and 2436 metabolites, and lacks any thermodynamically infeasible cycles. It was shown to be highly sensitive to environmental changes under both steady-state and dynamic conditions. As additional constraints, our dynamic model involved kinetic parameters associated with substrate consumption at different growth conditions (i.e., low CO2-heterotrophic and high CO2-mixotrophic). Our results suggest that cells growing at high CO2 (i.e., photoautotrophic and mixotrophic conditions) have an increased capability for biomass production. In addition, we have observed that ATP production also seems to be an important limiting factor for growth under the conditions tested. Our experimental data (metabolomics and transcriptomics) and the results predicted by our model clearly suggest a differential behavior between low CO2-heterotrophic and high CO2-mixotrophic growth conditions. The data presented in the current study contributes to better dissect the biological response of C. reinhardtii, as a dynamic entity, to environmental and genetic changes. These findings are of great interest given the biotechnological potential of this microalga for CO2 fixation, biomass accumulation, and bioenergy production.


INTRODUCTION
Atmospheric CO 2 concentration has significantly increased over the last 100 years as a result of numerous anthropogenic activities.This has strongly contributed to climate change and global warming leading to a severe environmental crisis (Sayre, 2010).In the recent years, the use of photosynthetic organisms able to convert atmospheric CO 2 into metabolites and other compounds of interest arose as an important alternative to mitigate such problem.In addition, CO 2 -fixing organism have been also shown to be a suitable platform for the production of inexpensive bulk chemicals with potential commercial interest using the atmospheric CO 2 as a feedstock (Boyle & Morgan, 2009).
Among CO 2 -fixing organisms, aquatic photosynthetic microorganisms have gained increased attention since they perform almost 50% of global photosynthesis.Even when the diffusion of CO 2 in aqueous solution has been reported to be a thousand times slower as compared to its diffusion on air (Moroney & Somanchi, 1999), several aquatic photosynthetic microorganisms have evolved CO 2 concentrating mechanisms (CCMs) to adapt to their environment.These mechanisms increase photosynthetic productivity by scaling up the levels of inorganic carbon over the environmental levels of CO 2 (Moroney & Ynalvez, 2007).One of these organisms is the algae Chlamydomonas reinhardtii, which is considered a promising organism with a high potential for CO 2 fixation and bioenergy production.This well-known microorganism has been used as a model to study the cellular and biochemical response to a large number of environmental conditions including nutrient deprivation (Harris, 2001).For example, it has been shown that lipid accumulation is favored once C. reinhardtii is grown under conditions of nitrogen starvation (Scranton et al., 2015).Moreover, it has been shown that modifications in the culture conditions of C. reinhardtii significantly influence its final biomass (Rosenberg et al., 2008).Interestingly, this algae has been successfully photoautotrophically, heterotrophically, and mixotrophically cultivated (Rupprecht, 2009).
As explained above, C. reinhardtii have evolved CCMs in order to increase intracellular CO 2 concentrations by using active C i -transport, thereby facilitating the rate of photosynthetic CO 2 fixation even at low concentrations of external C i .Given that this microalga must overcome the 1000-fold slower diffusion of CO 2 in water than in air, the active transport and accumulation of C i , either as CO 2 or as HCO 3 − , plays a critical role (Moroney & Ynalvez, 2007).Carbonic anhydrases (CAs) has been shown to play a major role in the CCM because they catalyze the interconversion of CO 2 and HCO 3 − , thereby maintaining and accelerating transport across the plasma membrane and to chloroplasts (Fang et al., 2012).
Significant advances have been made in order to elucidate the Chlamydomonas' CCM, including not only components of C i uptake systems, but additional regulatory factors (Yamano & Fukuzawa, 2009).For example, the expression of these C i transporters was shown to be regulated by several transcription factors that are responsive to changes in the environmental CO 2 concentration (Winck, Páez Melo & González Barrios, 2013).Under CO 2 limiting conditions, cells capture more C i by altering the expression of thousands of genes, which, in turn, may be involved in the acceleration and enhancement of C i acquisition.In other words, low CO 2 concentrations has been reported to induce the C. reinhardtii CCM, which have facilitated the identification of several additional CCMrelated genes by determining their expression level under limiting CO 2 (lower than 0.05%) as compared to high CO 2 (1% to 5%) (Fang et al., 2012).Moreover, additional factors such as the cooperation between light and carbon signaling have been shown to be also necessary for the modulation of CCM-related genes (Winck, Páez Melo & González Barrios, 2013).In this regard, our group has previously identified a vast range of genes and proteins that integrate carbon-related mechanisms (Winck, Páez Melo & González Barrios, 2013).
C. reinhardtii shows three CO 2 acclimation states: high concentration (>5% [CO 2 ]), low concentration (0.04%-0.03% [CO 2 ]), and very low concentration (<0.01%[CO 2 ]) (Fang et al., 2012).Therefore, an understanding of the metabolic pathways involved in the acquisition and accumulation of CO 2 over these stages is fundamental for the identification of metabolites that contribute to acclimation, cellular growth, and biomass production in response to changes in environmental CO 2 concentrations (Boyle & Morgan, 2009).In this regard, in silico analysis to model these metabolic pathways and its behaviors under variable conditions became a relevant strategy.Experimental validation of these models improves their predictive capabilities because not all metabolic pathways are necessarily activated in response to a specific cellular response or adaptation (Schmidt et al., 2010).Therefore, experimental data are essential for describing the substrates and products present in the cell.In this regard, the information obtained through experimental approaches from OMICS (i.e., transcriptomics, proteomics, and metabolomics) can be integrated into a reconstructed metabolic network model by the inclusion of new constraints as described by Kauffman, Prakash & Edwards (2003).Application of these restrictions reduces the solution space of possible phenotypes, thereby allowing a more accurate reconstruction and characterization of the metabolic network (Blazier & Papin, 2012).All these strategies can improve our capacity to predict the cellular behavior and identify pathways involved in biomass accumulation based on the observation of effects in the outputs.This is crucial to develop new strategies for improving the ability of C. reinhardtii to acquire and accumulate carbon (Winck, Páez Melo & González Barrios, 2013).
This work arose from our previously developed network for C. reinhardtii (Winck et al., 2016) and compares metabolomics and transcriptomics data with results predicted by our reconstructed genome-scale model (GSM) at steady-state and dynamic conditions at different CO 2 levels.Our main goal was to improve the capacity to predict metabolic routes for biomass accumulation in C. reinhardtii.For this purpose, we first curated the reconstructed model for subsequent analysis mainly of pathways sensitive to CO 2 and genes related to carbon uptake and accumulation applying a steady-state and dynamic model.Then, we generated a phenotypically accurate and predictive computational model for C. reinhardtii, which represents a major advance over previous models (Imam et al., 2015;Winck et al., 2016).
Our predicted results improved our understanding of the biological response of C. reinhardtii, as a dynamic entity, to environmental changes.So far, almost all physical networks reported in the literature have been examined under a single static condition.Nevertheless, biological systems are highly dynamic and static approaches have failed to identify and integrate parameters that are condition specific (Ideker & Krogan, 2012).Considering the effects of CO 2 concentration from a metabolic dynamic perspective is of great significance given the well-known biotechnological potential of C. reinhardtii for CO 2 fixation, biomass accumulation, and bioenergy production.

Genome-scale metabolic network reconstruction
GSM reconstruction involved the following steps: (i) Identification of metabolic pathways and membrane transport reactions based on the genomic annotation, enzyme homology, and experimental observations.This step used databases such as KEGG, ModelSEED, TCD, and TransportDB, among others; (ii) Development of equations for biomass components based on physiological data or comparison between similar species; (iii) Identification of incomplete pathways and transport reactions of missing metabolites through semi-automated methods, such as reverse engineering, or automated methods, such as GapFind and GapFilling (Kumar, Dasika & Maranas, 2007;Senger & Papoutsakis, 2008;Serrano-Bermúdez, 2016).In this work, we first improved our previously developed network for C. reinhardtii (Winck et al., 2016), which required the extraction of all the available metabolic information (e.g., reactions, compounds, pathways, genes, enzymes) from the following public databases: Plant Metabolic Network (ChlamyCyc 5.0, 2015), KEGG (Kanehisa Laboratories, 2016), BioCyc (BioCyc Database Collection, 2016), and ModelSEED (ModelSEED, 2015).In addition, the iCre1355 model (Imam et al., 2015) was used to enhance our reconstruction.Manual curation was required through the comparison between models in order to verify the presence of novel substrates, products, and reactions by comparing several models.Manual gap-filling was carried out to fill additional gaps using literature and experimental data from different databases.Thus, manual curation was used to resolve issues related to reaction stoichiometries of reactions not found in the databases, or genes coming from extensive bioinformatic analysis not yet submitted.Reactions directionality were determined as reported by Mavrovouniotis and co-workers (Mavrovouniotis, 1990).Each metabolite in our GSM model may have one or more designations denoting its cellular compartment, which allowed the inclusion of exchange reactions designed to connect novel metabolites in their cellular compartments.

Chlamydomonas reinhardtii kinetic modeling
Kinetic parameters associated with substrate consumption at different growth conditions (i.e., 0% CO 2 -heterotrophic and high CO 2 -mixotrophic) were used as additional constraints in the dynamic model.Then, based on the growth conditions, we constructed two Chlamydomonas kinetic models: (i) heterotrophic growth and (ii) mixotrophic growth.

Heterotrophic dynamic model for cell growth and substrate consumption
In this work, we used the Haldane model, which includes terms to describe the effects of a nutrient on the growth rate at high concentrations (Zhang, Chen & Johns, 1999;Lee, Jalalizadeh & Zhang, 2015).
where X is algal density or concentration and m is the maintenance energy coefficient.

Mixotrophic dynamic model for cell growth and substrate consumption
The behavior of microalgal cultures was described with a set of nonlinear algebraicdifferential equations deduced from mass balance considerations for both liquid and gaseous phases, assuming well-mixed conditions, and following the Haldane model (Tebbani et al., 2013).
out is the CO 2 molar fraction in the output gas, and K 1 and K 2 are equilibrium constants.

Constraint-based modeling
Flux Balance Analysis (FBA) FBA assumes that metabolic networks achieve a steady state that is constrained by the stoichiometry and that they keep on self-adjusting and optimizing in a changing external environment.This method then uses a mathematical model to simulate the network self-adjustment (Kauffman, Prakash & Edwards, 2003;Chen, Wang & Wei, 2010).Here, we reconstructed our model following the same parameters described in our previous work (Winck et al., 2016), but for all new tested conditions (i.e., photoautotrophic, heterotrophic, and mixotrophic conditions).
Then, CO 2 fluxes were changed in order to determine the impact of the different tested CO 2 concentrations supply at the steady state of cell metabolism.The biomass function was exactly the same as described in our previous study (Winck et al., 2016).The LBs and UBs of the reactions were fixed according to standard values, the criteria of the free Gibb's energy at 27 • C and pH 7, and the database information.A maximization of the biomass reaction flux was performed, and the optimization problem was resolved to employ the GAMS R (General Algebraic Modeling System: https://www.gams.com/)distribution 23.9.5 using FBA by solving the linear programming problem with constraints set as follows: where S is the stoichiometric matrix, v is the flux vector, LB is the lower bound, UB is the upper bound, c is a vector of zeros that sets the objective function, m is the number of metabolites, n is the number of reactions, and pos( 1) is the objective function (biomass) (Gianchandani, Chavali & Papin, 2010).Finally, those reactions that changed at different CO 2 levels were identified by FBA.In that way, we were able to determine candidate metabolic pathways of relevance for a CO 2 response.

Flux Variability Analysis (FVA)
Linear program (LP) problems arising in FBA can have multiple solutions with the exact same optimal value for the objective function and still satisfy all the constraints.The number of these solutions determines the dimensionality of the alternative optimal space.These variables can be identified using flux variability analysis (FVA), where the flux of each reaction in the network is maximized and minimized, one at a time, while the optimal value is obtained from FBA (Mahadevan & Schilling, 2003;Maranas & Zomorrodi, 2016).
Here, the optimization problem was as follows: maximize(and minimize) z = v j subject to FVA was then used to identify blocked reactions in a metabolic network under a given environmental and growth condition.Blocked reactions were reactions with both a minimum and maximum flux value equal to zero under the examined condition.The relevance of identifying blocked reactions was: (i) that they could be used as a basis for improving the quality of metabolic network reconstruction by helping to identify gaps, or (ii) that they could significantly reduce the search space of FBA by pre-setting their value to zero (Maranas & Zomorrodi, 2016).

Thermodynamically infeasible cycles (TICs)
The loop law states that around a metabolic loop all the thermodynamic driving forces must add up to zero.In other words, a net flux cannot exist around a closed cycle in a network at steady state (Schellenberger, Lewis & Palsson, 2011).By using the FVA method a loopless analysis was performed, which consisted of a sequential maximization and minimization of each reaction in the model.Those reactions where the range (max v i -min v j ) differed by >10 −6 were defined as loop reactions (Schellenberger, Lewis & Palsson, 2011).The loopless condition used to identify all cycles can be written as the condition where G is a vector of energies for each reaction.All loops can be expressed as a linear combination of the null basis of S matrix in the form v = Nxα, where N = null(S), and αare weights.As described above, the loopless condition generates an LP problem.The driving force of each reaction was then indicated by the vector of continuous variables (G i ) .This quantity can be thought of as being analogous to the G of each reaction.As a loop is entirely defined by the direction of the flux distribution (Schellenberger, Lewis & Palsson, 2011), if NxG i = 0, no loop must be present.This methodology consisted of the detection of reactions with maximum and/or minimum possible values on fluxes (1,000 and −1,000 (mmol/g*h)) using the FVA.Subsequently, these reactions and the metabolites involved are expressed in matrix form, and the null space of the matrix is calculated.Each thermodynamically infeasible cycle (TIC) is thus detected with the reactions involved.The TICs are solved through manual revision of reversibility or irreversibility of the reactions involved, removal of repeated reactions, or as a last option, the blocking of some reactions.Finally, this procedure is repeated until the null space is zero; thus, the GSM reconstructed model is cured (Maranas & Zomorrodi, 2016).Here, the optimization problem (FVA) was resolved using GAMS R (General Algebraic Modeling System) distribution 23.9.5.
Computation for obtaining the null space of S was performed in MATLAB R version 2015b (MathWorks, Natick, MA, USA) followed by manual correction of reversibility in the involved reactions according to the databases sources.

Chlamydomonas reinhardtii metabolic network dynamic simulation
A variant of FBA, called dynamic flux balance analysis (DFBA), allows a description of the transience of metabolism.DFBA can be implemented in two ways: by a dynamic optimization approach (DOA), which requires solving a nonlinear program (NLP), and by a static optimization approach (SOA), which requires solving an LP (Gianchandani, Chavali & Papin, 2010).On one hand, DOA involves optimization over the entire time period of interest to obtain time profiles of fluxes and metabolite levels (Mahadevan, Edwards & Doyle, 2002).On the other hand, SOA distributes the time into different time intervals and solve the immediate optimization problem at the beginning of each of these intervals.This allows to predict the fluxes of each time point, which is followed by integration over the interval to compute species concentrations over time (Mahadevan, Edwards & Doyle, 2002;Gianchandani, Chavali & Papin, 2010).Here, we performed the dynamic model using SOA that involved kinetic parameters as additional constraints, and the optimization problem was resolved using GAMS R (General Algebraic Modeling System) distribution 23.9.5 by setting the following: subject to where z i is the extracellular metabolite concentration (based on the components of the culture media), z i,0 and X 0 are initial conditions for extracellular metabolite and algal concentrations, respectively, X is the algal concentration, µ is the specific growth rate, c j is the reaction weight, ĉ (z,v) is a vector of nonlinear constraints (i.e., kinetics parameters associated with substrate consumption), t 0 and t f are the initial and the final times, respectively, and G is the number of intervals in which the time is discretized.

Total biomass measurement
Cells were grown in both heterotrophic (0% CO 2 and acetic acid as a carbon source) and mixotrophic (10% CO 2 and acetic acid as carbon sources) growth conditions.Then, the total biomass was determined based on dry weights following the method described in our previous work (Winck et al., 2016).

Gene expression analysis by quantitative RT-PCR
Genes were selected based on our previously reported sensitivity analysis and included genes related to CCM, the Calvin Cycle, glycolysis/gluconeogenesis and glyoxylate/dicarboxylate metabolism (Winck et al., 2016).We performed a new sensitivity analysis using FBA and DFBA to contrast with the previous results and to describe other pathways involved using expression data (i.e., transcriptomic data).Cell pellets from 45 mL cell culture aliquots were collected on the last cultivation day in both growth conditions (heterotrophic and mixotrophic).Then, the total RNA was extracted by using the RNeasy Plant Mini Kit (QIAGEN 2012) and subsequently treated with TURBO DNAase (Ambion and Life Technologies, 2012) to remove remaining DNA.Synthesis of cDNA was performed using the iScript cDNA Synthesis Kit (BioRad, Hercules, CA, USA) according to the manufacturer's instructions.The actin gene [ACT] was used as a reference gene.All the cDNA samples were amplified in 96-well plates in Agilent Technologies Stratagene Mx3005P system.Real-time PCR was performed used the same conditions as described in Winck et al. (2016).For each gene, a relative expression ratio was determined as previously described by Schmittgen and Livak (Schmittgen & Livak, 2008).PCR efficiencies were determined as previously described (Winck et al., 2016).

Metabolomics analysis
Cell pellets from a culture in stationary phase were collected from the 0% CO 2heterotrophic and high CO 2 -mixotrophic cultures and metabolites were extracted and treated as previously described in Winck et al. (2016).Metabolites derivatization, GC-TOF-MS analysis, chromatogram analysis, peak detection was also performed used the methods previously described (Winck et al., 2016).

Significantly improved and robust genome-scale metabolic network model of C. reinhardtii
Based on our previously reconstructed network and biomass reactions (Winck et al., 2016), we initially performed a curation process by reverse engineering.This allows us to identify some inconsistencies, such as disruptions in the balance of proton (H + ) charges, disruptions in the balance of water consumption and production in some reactions, the oxygen requirement for growth, and the presence of the same reaction twice but one of them was the reverse reaction, among other factors.Initially, our model consisted of 3554 reactions and 2,342 metabolites, in contrast to iCre1355, which has 2,394 reactions and 1,845 metabolites (Imam et al., 2015).A cross-referencing process with iCre1355 allowed us to identify 172 new reactions, including reactions related to exchange, transport, amino acid synthesis, oxidative phosphorylation, and fatty acid biosynthesis, among others.Manual examination of KEGG and ModelSEED databases was used to link some products by applying gap filling.The reactions were compartmentalized into ten compartments: extracellular, cytosol, chloroplast, mitochondria, nucleus, eyespot, flagellum, glyoxysome, thylakoid lumen, and Golgi apparatus.Finally, we generated a new comprehensive genomescale network reconstruction for C. reinhardtii, which consisted of 3,726 reactions and 2,436 metabolites (Data S1, S2 and S3).
The main characteristics of our GSM model are summarized in Table 1.This model exhibits increases of 5% in the number of reactions over our previous model, and 56% in the number of reactions over the iCre1355 model, as well as increases of 4% and 32% in the number of metabolites over the previous model and the iCre1355 model, respectively.The new model we are presenting here also has 1,416 unique metabolites, which allows us to encompass a broader range of metabolic functions.As shown in Fig. 1, all reactions presented in the new reconstructed model were distributed across several compartments, including transport and exchange reactions.Compartmentalization was necessary to guarantee the location and viability of reactions as reported in databases and literature (Petersen et al., 2011;Winck et al., 2016).The interconnections in the current metabolic network are shown in Fig. 2. As in the previously developed model and the iCre1355 model, the reactions in the current model are distributed over nine intracellular compartments and the extracellular space, with the majority of reactions localized to the cytosol, chloroplast, and mitochondria.Interestingly, the chloroplast was shown to be highly relevant in light-driven metabolism due to its relative centrality (Figs. 1 and 2).In fact, the chloroplast represents 25% of the total reactions in the network.These compartments involve crucial pathways for photoautotrophic growth, including photosynthesis, chlorophyll synthesis, and the mechanism for phototaxis (Chang et al., 2011).

Cell growth and substrate consumption behavior
Cell growth and substrate consumption profiles were obtained in batch culture under heterotrophic and mixotrophic conditions in order to utilize them as restrictions for the dynamic model (Figs. 3 and 4).pH variations ranging from 6.3 and 7.5 were chosen, based on previous kinetic modeling for the heterotrophic growth of C. reinhardtii in batch  cultures using acetate as carbon source and ammonium as nitrogen source (Zhang, Chen & Johns, 1999).A lag phase of around 10 hours was observed and acetate was completely exhausted after 70 hours.The resulting cell concentration at heterotrophic conditions was around 0.55 g/L for the model and 0.58 g/L in our experimental data, which was in the same range and similar to the values reported by others (0.44 g/L) (Zhang, Chen & Johns, 1999).Likewise, the response of dissolved inorganic carbon concentration agreed with previous measurements of CO 2 fluxes in Chlamydomonas under similar conditions (Badger, Palmqvist & Yu, 1994).The increasing cell concentration under mixotrophic conditions was due to the presence of both carbon sources (acetate and CO 2 ) and the continuous input of CO 2 into the media.

The curated genome-scale metabolic network model of C. reinhardtii is highly sensitive under both steady-state and dynamic conditions
The guarantee of a robust and sensitive model to predict metabolic fluxes and to use in subsequent dynamic simulations required the correction of all TICs, along with a balancing of H + charge and water consumption in the entire network.A total of 187 loops were identified and many cycles were corrected and eliminated to finally generate a metabolic network curated without any thermodynamically infeasible cycle.The final reconstructed model was used along with the FBA to predict the flux distribution in C. reinhardtii under photoautotrophic (using light and CO 2 as energy and carbon sources), heterotrophic (using acetate as carbon source), and mixotrophic (using both acetate and CO 2 as carbon sources, and light) growth conditions.This new sensitivity analysis implemented FBA and DFBA to describe new pathways using transcriptomic data, which provided a reference for optimizing biomass production upon different growth conditions.The predictions were performed with maximization of biomass.The reactions we have found showed significant variation depending on the CO 2 concentration and confirmed effects on different metabolic pathways, such as glycerolipid synthesis, amino acid synthesis, glycolysis, carbohydrate catabolism, and carbon fixation, among others (Data S2).We assessed the ability of this model to predict growth rates by performing simulations under photoautotrophic, heterotrophic, and mixotrophic conditions.Fluxes of carbon dioxide, acetate, and light were adjusted according to the experimental evidence (Chang et al., 2011;Chaiboonchoe et al., 2016) to predict growth rates in the model (Table 2).Our results showed a significant coherence between conditions when compare the experimental and the in silico growth rates, with the lowest growth rate under heterotrophic conditions and the highest under the mixotrophic condition, which is in agreement with several previous reports (Moon et al., 2013;Imam et al., 2015).The growth rates predicted by FBA in terms of the CO 2 and acetate input fluxes are shown in Fig. S1.The model predicts optimal values of acetate input fluxes close to 100 mmol/gDW h, without restriction of ATP production in the mitochondria or the chloroplasts.We have found that ATP production is an important restriction for growth that is dependent on environmental conditions.Response surfaces of growth rates predicted by FBA with respect to fluxes of acetate and CO 2 and the reactions of ATP production in the mitochondria and chloroplasts are shown in Fig. 5. Simulation of the different growth conditions required blockade  of ATP production in mitochondria under photoautotrophic conditions, blockade of ATP production in the chloroplast under heterotrophic conditions, and activity of both reactions under mixotrophic conditions.FBA suggested that at all tested conditions (i.e., heterotrophic, photoautotrophic, and mixotrophic conditions), the metabolic pathways related to photorespiration persist active, which fulfilled the model constraints.
The magnitude of the fluxes for glycine transport from chloroplasts to mitochondria, serine transport from mitochondria to chloroplasts, and production of hydroxy pyruvate, glycerate, and 3-phosphoglycerate in the chloroplast showed positive values.These routes contribute to the maximization of biomass production because they are related to photorespiration processes (Kebeish et al., 2007).FBA revealed an increase in the flux for metabolic pathways happening in the chloroplast and mitochondrial transport systems in response to CO 2 input variations as previously reported (Winck et al., 2016).Moreover, our results suggested that cells at high CO 2 (photoautotrophic and mixotrophic conditions) have the ability to increase their biomass production.

DISCUSSION
The main goal of this study was to understand the biological response of C. reinhardtii to environmental changes in CO 2 concentrations from a metabolic dynamics perspective.For this reason, we performed a sensitivity analysis to provide a reference for optimizing biomass production in response to variations in heterotrophic, photoautotrophic, and mixotrophic growth conditions by changing the carbon sources available in the extracellular medium.
We also wanted to test the effects of changing the initial concentration of inorganic compounds, such as photons and water, in the extracellular medium to promote the photosynthesis process and biomass production.Some reports on the energy requirements for photoautotrophic growth of C. reinhardtii indicated that, at high light supply, the biomass yields decrease due to light saturation effects (Kliphuis et al., 2012).Light saturation effects also have been shown to be depended on DIC availability (Chang et al., 2011).The optimal value between low and high light supply still has to be elucidated, which represents an important future work.The dynamic predictions using SOA agreed with the experimental data of biomass production where, at day three of cultivation, we obtained a biomass concentration of 1.3 g/L.The current reconstructed model used the same achieved values of the SOA, for less cultivation time, with the improved biomass production showing consumption of both carbon sources.Our model predicted mixotrophic growth accurately: acetate is used initially, but as acetate is consumed and its concentration is reduced, the CO 2 consumption starts (Fig. 6).Our model was able to predict the behavior of variations in the synthesis of some amino acids (lysine, glycine, leucine, tyrosine, isoleucine, valine, and phenylalanine); some of these are biomass precursors and showed the most prominent alterations in response to high CO 2 concentrations (Fig. 7).In the same way, our dynamic model predicts the fluxes of important transport reactions that involved relevant cofactors and substrates of the TCA cycle and carbon fixation mechanism.These were activated pathways mainly associated with the chloroplast and mitochondria and included hydroxypyruvate, 3-phosphoglycerate, ribulose 1,5-biphosphate, aspartate, and phosphoenolpyruvate production.All these processes may be related to photorespiration and showed positive values in our current reconstructed model.All these routes were shown to enhance biomass production under mixotrophic conditions.It is important to mention that CO 2 will only increase the biomass if the cells are not limited for other nutrients like N or P.

Metabolomic and transcriptomic data
We performed a canonical discriminant evaluation of the experimental data to compare with our dynamic modeling results.The canonical discriminant analysis allows evaluation of the presence of significant differences between groups of a set of measured variables in order to classify them.Figure S2 presents the canonical discriminant analysis grouping by growth condition alone.Figure S3 shows the canonical discriminant functions analysis obtained after grouping by both time and growth condition.A substantial discrimination was revealed by the opposite ranges of scores (Fig. S2), suggesting that growth condition is a function that discriminates in a proper way.In Fig. S3, the centroids show the CO 2 groups in culture (0%-heterotrophic, 10%-mixotrophic) considering the metabolic abundance.Metabolites and abundance were shown to be significantly affected by changes in CO 2 concentration level.Our dynamic modeling agreed with these results since it showed a significant increase of flux in metabolic routes including the TCA cycle, glycolysis/gluconeogenesis, and amino acid biosynthesis at high CO 2 concentrations, as confirmed by the experimental abundance of metabolites involved in these pathways.
Our data suggested that high CO 2 concentration may activate mechanisms that are able to control carbon fixation.
The dynamic model predicted the profile only for the extracellular metabolites, so this required obtaining the fluxes of the reactions that mainly contributed to the experimentally observed metabolic abundance in order to make appropriate comparisons.The metabolomic data primarily showed the presence of amino acids and few sugars; therefore, to compare these data with the results predicted by our model, we primarily analyzed reactions involved in amino acid biosynthesis.Figure 7 shows the predicted reaction fluxes associated with amino acid biosynthesis, where a differential behavior was clearly observed between growth conditions (high values for reaction fluxes at high CO 2 concentrations and low values for reaction fluxes at low CO 2 concentrations).This is in agreement with the experimental results of metabolic abundance that showed opposing behavior between the two CO 2 concentrations.This suggests that our proposed metabolic network model is highly sensitive and predictive of environmental changes in both steady-state and dynamic conditions.
The gene expression analysis performed by RT-PCR allowed the identification and relative quantification of genes involved in the carbon concentrating mechanism, TCA cycle, glycolysis, and gluconeogenesis.The expression level of those genes encoding carbonic anhydrases is shown in Fig. 8A.Here, we analyzed these genes since the previous evidence was indicating that they were of relevance in carbon uptake and biomass accumulation (Winck et al., 2016).Our results showed that gene transcripts for CAH1, CAH4, CAH5, and LCIA were more abundant at high CO 2 concentrations, in agreement with previous reports of carbon uptake mechanisms in C. reinhardtii (Brueggeman et al., 2012).Fig. 8B presents the expression levels of genes involved in glycolysis/gluconeogenesis and some of TCA cycle.Genes encoding for glycerate kinase (GLYK), ATP synthase alpha subunit (ATP1A), fructose 1,6-biphosphatase (FBP1) showed more abundant gene transcripts at high CO 2 concentration.On the other hand, genes encoding for pyruvate decarboxylase (PDC3), NAD-dependent malate dehydrogenase (MDH3) did not show a differential behavior.These results are in agreement with our steady-state and dynamic simulations, which showed abundant gene transcripts in crucial genes involved mainly in carbon uptake mechanisms and the TCA cycle.Furthermore, they were related to the variation in differential growth conditions, which was confirmed in the fluxes of these reactions (including transport and exchange reactions) in our current reconstruction model (Fig. 9 and Data S3).
Our predicted results led us to consider new target genes that own the crucial plasticity needed to support growth adaptability.Figure 9 shows the reaction fluxes predicted for other target genes associated with the carbon concentrating mechanism, photosynthesis, glycolysis/gluconeogenesis, and the TCA cycle.It is important to note that we assumed that the reaction flux is directly related to the transcript abundance, given the enzyme activity, so differential behavior is clearly observed between growth conditions.Our predicted results suggest that transcripts for GAP3 (glyceraldehyde-3-phosphate dehydrogenase), RBCS1 (ribulose-bisphosphate carboxylase), PGK (phosphoglycerate kinase), PGI1 (glucose-6-phosphate isomerase), and PGH1 (enolase) must be more abundant at high CO 2 concentrations.In the same way, our predicted results showed high and unique expression of genes related to photosynthesis under photoautotrophic and mixotrophic growth conditions.We have found that ATP production is an important restriction for growth, which depends on environmental conditions, in agreement with previous flux balance analysis in metabolic networks (Serrano-Bermúdez et al., 2017).This result suggests that our proposed metabolic network model could give complementary and high-value perspectives on these transcriptional changes.

CONCLUSIONS
In this work, we have significantly improved the previous models, obtaining a new robust genome-scale metabolic network model for the microalga Chlamydomonas reinhardtii.
Our current model is appropriately curated and lacks any thermodynamically infeasible cycles, and it is highly sensitive to environmental changes in both steady-state and dynamic conditions, with kinetic parameters associated with substrate consumption at different growth conditions (0% CO 2 -heterotrophic and high CO 2 -mixotrophic).The experimental results agreed with steady-state and dynamic simulations, showing abundant gene transcripts in crucial genes involved in carbon uptake mechanisms and the TCA cycle.Our proposed metabolic network model can, therefore, provide complementary and high-value perspectives on transcriptional changes, metabolic abundance, and biomass production under different environmental conditions.Finally, up to now, no reports exist regarding dynamic flux balance analysis in C. reinhardtii.Consequently, our results and predictions can improve the understanding of the biological response of C. reinhardtii, as a dynamic entity, to environmental changes.This is of great interest given the biotechnological potential of this organism for CO 2 fixation, biomass accumulation, and bioenergy production.
Model predictions with FBA allows the description of different metabolic states regarding the variations in the growth conditions.In addition, the dynamic approaches have revealed a robust model with a high predictive capacity, which is a useful tool to further analyze a priori effects of several perturbations on the culture media, limiting factors, metabolic pathways, potential mutations, biomass production, lipids accumulation, and so on.Thus, our model represents a useful contribution to the field since it will allow to propose, analyze and evaluate different conditions affecting biological systems from a dynamic point of view.This is of relevance as biological systems are highly dynamic entities affected by a wide variety of changes in the environment, as a part of an adaptation process.

Figure 1
Figure 1 Subcellular localization of all reactions in the current reconstructed model (including exchange and transport reactions).Full-size DOI: 10.7717/peerj.5528/fig-1