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

Daniela Alejandra Mora Salguero; Miguel Fernández-Niño; Luis Miguel Serrano-Bermúdez; David O. Páez Melo; Flavia V. Winck; Camila Caldana; Andrés Fernando González Barrios

doi:10.7717/peerj.5528

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

Daniela Alejandra Mora Salguero¹, Miguel Fernández-Niño¹, Luis Miguel Serrano-Bermúdez², David O. Páez Melo¹, Flavia V. Winck³, Camila Caldana^4,5, Andrés Fernando González Barrios ¹

1Grupo de Diseño de Productos y Procesos (GDPP), Department of Chemical Engineering, Universidad de Los Andes, Bogotá, Colombia

2Bioprocesses and Bioprospecting Group, Universidad Nacional de Colombia, Bogotá D.C., Colombia

3Laboratory of Regulatory Systems Biology, Department of Biochemistry, Institute of Chemistry, Universidade de São Paulo, São Paulo, Brazil

4Brazilian Bioethanol Science and Technology Laboratory, Brazilian Center for Research in Energy and Materials, Campinas, Brazil

5Max Planck Partner Group, Brazilian Bioethanol Science and Technology Laboratory, Brazilian Center for Research in Energy and Materials, Campinas, Brazil

DOI: 10.7717/peerj.5528

Published: 2018-09-03
Accepted: 2018-08-06
Received: 2018-02-10

Academic Editor: Michael Henson

Subject Areas: Computational Biology, Climate Change Biology
Keywords: Chlamydomonas reinhardtii, CO₂ fixation, Metabolomics, Biomass, Transcriptomics, Biotechnology, Metabolic network, Dynamic model

Copyright: © 2018 Mora Salguero et al.
Licence: This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.

Cite this article: Mora Salguero DA, Fernández-Niño M, Serrano-Bermúdez LM, Páez Melo DO, Winck FV, Caldana C, González Barrios AF. 2018. Development of a Chlamydomonas reinhardtii metabolic network dynamic model to describe distinct phenotypes occurring at different CO₂ levels. PeerJ 6:e5528 https://doi.org/10.7717/peerj.5528

Abstract

The increase in atmospheric CO₂ 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 CO₂ 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 CO₂ fixation, biomass accumulation and bioenergy production upon exposure to several environmental conditions. In the present study, we studied the Chlamydomonas response to different CO₂ 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 CO₂. 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 CO₂-heterotrophic and high CO₂-mixotrophic). Our results suggest that cells growing at high CO₂ (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 CO₂-heterotrophic and high CO₂-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 CO₂ fixation, biomass accumulation, and bioenergy production.

Introduction

Atmospheric CO₂ 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₂ into metabolites and other compounds of interest arose as an important alternative to mitigate such problem. In addition, CO₂-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₂ as a feedstock (Boyle & Morgan, 2009).

Among CO₂-fixing organisms, aquatic photosynthetic microorganisms have gained increased attention since they perform almost 50% of global photosynthesis. Even when the diffusion of CO₂ 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₂ 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₂ (Moroney & Ynalvez, 2007). One of these organisms is the algae Chlamydomonas reinhardtii, which is considered a promising organism with a high potential for CO₂ 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₂ concentrations by using active C_i-transport, thereby facilitating the rate of photosynthetic CO₂ fixation even at low concentrations of external C_i. Given that this microalga must overcome the 1000-fold slower diffusion of CO₂ in water than in air, the active transport and accumulation of C_i, either as CO₂ or as HCO₃⁻, 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₂ and HCO₃⁻, 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₂ concentration (Winck, Páez Melo & González Barrios, 2013). Under CO₂ 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₂ concentrations has been reported to induce the C. reinhardtii CCM, which have facilitated the identification of several additional CCM-related genes by determining their expression level under limiting CO₂ (lower than 0.05%) as compared to high CO₂ (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₂ acclimation states: high concentration (>5% [CO₂]), low concentration (0.04%–0.03% [CO₂]), and very low concentration (<0.01% [CO₂]) (Fang et al., 2012). Therefore, an understanding of the metabolic pathways involved in the acquisition and accumulation of CO₂ 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₂ 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₂ 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₂ 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₂ concentration from a metabolic dynamic perspective is of great significance given the well-known biotechnological potential of C. reinhardtii for CO₂ fixation, biomass accumulation, and bioenergy production.

Materials and Methods

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₂-heterotrophic and high CO₂-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). $μ = μ_{m a x} \frac{[A c]}{K_{s} + [A c] + \frac{{[A c]}^{2}}{K_{i}}}$ $\frac{d X}{d t} = μ_{m a x} \frac{[A c]}{K_{s} + [A c] + \frac{{[A c]}^{2}}{K_{i}}} X$ $- \frac{d [A c]}{d t} = \frac{1}{Y_{x}} \frac{d X}{d t} + m X$ where X is algal density or concentration [g∕L], [Ac] is substrate (acetate) concentration [g∕L], t is time [h], μ_m is maximum specific growth rate [1∕h], K_s is the acetate saturation constant [g∕L], K_i acetate inhibition constant [g∕L], Y_x is yield constant on acetate $(\frac{g A l a g e p r o d u c e d}{g o f n u t r i e n t}) [g ∕ g]$ , 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 algebraic-differential 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).

$\frac{d X}{d t} = μ_{m a x} \frac{[C_{D I C}]}{K_{s C_{D I C}} + [C_{D I C}] + \frac{{[C_{D I C}]}^{2}}{K_{i C_{D I C}}}} * \frac{[A c]}{K_{s A C} + [A c] + \frac{{[A c]}^{2}}{K_{i A C}}} X$ $\frac{d [C_{D I C}]}{d t} = - \frac{1}{M_{x}} \frac{d X}{d t} + N_{C O_{2}}$ $\frac{d [A c]}{d t} = - \frac{1}{Y_{x}} \frac{d X}{d t} + m X$ $N_{C O_{2}} = k_{L} a_{C O_{2}} (\frac{P}{H_{C O_{2}}} y_{o u t}^{C O_{2}} - [C O_{2}])$ $[C O_{2}] = \frac{[C_{D I C}]}{(1 + K_{1} 1 0^{p H} + K_{1} K_{2} 1 0^{2 p H})}$ where [C_DIC] is dissolved inorganic carbon concentration [mol∕L], K_{sC_DIC} is the C_DIC saturation constant [mol∕L], K_{iC_DIC} is the C_DIC inhibition constant [mol∕L], M_x is the C-mole mass, N_CO₂ is the volumetric mass transfer rate for CO₂ $[\frac{m o l}{m^{3} * h}]$ , k_La_CO₂ is the volumetric mass transfer coefficient [1∕h], P is the pressure $[P a]$ , H_CO₂ is the CO₂Henry constant [Pa∗m³∕mol], $y_{o u t}^{C O_{2}}$ is the CO₂ molar fraction in the output gas, and K₁ and K₂ 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₂ fluxes were changed in order to determine the impact of the different tested CO₂ 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^® (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: $m a x_{z} c^{T} v (1)$ $subject to S v = 0$ $L B \leq v \leq U B$ $c^{T} \in R^{n} | c^{T} = [000 \dots 1 \dots 000] p o s (1) = B i o m a s s r e a c t i o n$ $v \in R^{n}$ $S \in R^{m x n}$ LB ∈ Rⁿ, UB ∈ Rⁿ

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 $p o s (1)$ is the objective function (biomass) (Gianchandani, Chavali & Papin, 2010). Finally, those reactions that changed at different CO₂ levels were identified by FBA. In that way, we were able to determine candidate metabolic pathways of relevance for a CO₂ 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$ $\sum_{j ϵ J} S_{i j} v_{j} = 0, \forall i ϵ I$ $L B_{j} \leq v_{j} \leq U B_{j}, \forall j ϵ J$ $v_{b i o m a s s} = f v_{b i o m a s s}^{m a x}$ $v_{j} ϵ R, \forall j ϵ J .$

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⁻⁶ 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^{^®} (General Algebraic Modeling System) distribution 23.9.5. Computation for obtaining the null space of S was performed in MATLAB^® 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^{^®} (General Algebraic Modeling System) distribution 23.9.5 by setting the following: $m a x \sum_{j = 1}^{N} c_{j} \cdot v_{j} (t)$ $subject to$ $z_{i} (t + Δ t) = z_{i} (t) - \sum_{j = 1}^{N} S_{i j} \cdot v_{j} (t) \cdot X (t) \cdot Δ t \forall i ϵ 1, \dots, M_{e x t r a c e l l u l a r}$ $X (t + Δ t) = X (t) + μ \cdot X (t) \cdot Δ t$ $\sum_{j = 1}^{N} S_{i j} \cdot v_{j} = 0 \forall i ϵ 1, \dots, M_{e x t r a c e l l u l a r}$ $v_{j}^{m i n} < v_{j} (t) < v_{j}^{m a x} \forall j ϵ 1, \dots, N$ $\hat{c} (z_{i} (t), v_{j} (t)) \leq 0$ $z_{i} (t) \geq 0 z_{i} (t_{0}) = z_{i, 0} \forall i ϵ 1, \dots, M_{e x t r a c e l l u l a r}$ $X (t) \geq 0 X (t_{0}) = X_{0}$ $Δ t = \frac{t_{f} - t_{0}}{G} \forall G ϵ 0 \dots G$ $\forall t_{g} ϵ [t_{0}, t_{f}]$ where z_i is the extracellular metabolite concentration (based on the components of the culture media), z_i,0 and X₀ 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, $\hat{c} (z, v)$ is a vector of nonlinear constraints (i.e., kinetics parameters associated with substrate consumption), t₀ and t_f are the initial and the final times, respectively, and G is the number of intervals in which the time is discretized.

Cell strain and culture conditions

The Chlamydomonas reinhardtii strain CC503 cw92mt+ (Chlamydomonas Resource Center University of Minnesota, USA) was cultivated in tris-acetate-phosphate (TAP) medium (Anderson, 2005): phosphate solution (1 ml) (K₂HPO₄, 28.8 g 100 ml⁻¹; KH₂PO₄, 14.4 g 100 ml⁻¹), salts (25 ml) (NH₄Cl, 15g L⁻¹; MgSO₄ 7H₂O, 4g L⁻¹; CaCl₂ 2H₂O, 2g L⁻¹), trace elements (1 ml) (H₃BO₃,1.14 g 100 ml⁻¹; ZnSO₄ 7H₂O, 2.2 g 100 ml⁻¹; MnCl₂ 4H₂O, 0.5 g 100 ml⁻¹; FeSO₄ 7H₂O, 0.5 g 100 ml⁻¹; CoCl₂ 6H₂O, 0.16 g 100 ml⁻¹; CuSO₄ 5H₂O, 0.16 g 100 ml⁻¹; (NH₄)₆Mo₇O₂₄ 4H₂O, 0.11 g 100 ml⁻¹; EDTA 5 g 100 ml⁻¹), acetic acid (1 ml), and tris (2.42 g/L) at 27 °C in heterotrophic and mixotrophic growth conditions, under different CO₂ concentrations [0% and 10% CO₂, respectively], in a R’ALF Plus bioreactor (Bioengineering, Inc., Cambridge, MA, USA) under continuous illumination. Cell growth, bioreactor experiments and counting of cells were performed as previously described in Winck et al. (2016).

Total biomass measurement

Cells were grown in both heterotrophic (0% CO₂ and acetic acid as a carbon source) and mixotrophic (10% CO₂ 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₂-heterotrophic and high CO₂-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).

Results

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 genome-scale 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).

Table 1:

Comparison of components in the models.

Components	Previous model	iCre1355	Current model
Reactions	3,554	2,394	3,726
Transport	479	426	498
Exchange	72	55	91
Compartments	10	10	10
Metabolites	2,342	1,845	2,436
Unique metabolites	1,382	1,133	1,416
Cytosol	955	730	1,008
Chloroplast	626	478	638
Mitochondria	388	308	423

DOI: 10.7717/peerj.5528/table-1

Figure 1: Subcellular localization of all reactions in the current reconstructed model (including exchange and transport reactions).

Download full-size image

DOI: 10.7717/peerj.5528/fig-1

Map Chlamydomonas reinhardtii. Metabolic Network model. — Figure 2: Map *Chlamydomonas reinhardtii.* Metabolic Network model.
(Nodes: Metabolites, Edges: Reactions). Green, cytosol; light blue, chloroplast; yellow, Mitochondria; maroon, glyoxysome; orange, flagellum; red, thylakoid lumen; pink, Extracellular; dark blue, Nucleus; violet, golgi apparatus; gray, eyespot. This graphic map was performed in the visualization software Cytoscape (http://www.cytoscape.org/).

Download full-size image

DOI: 10.7717/peerj.5528/fig-2

Figure 3: Profiles predicted by kinetic modeling at heterotrophic growth conditions (A) acetate consumption and (B) cell concentration.

Download full-size image

DOI: 10.7717/peerj.5528/fig-3

Time courses of (A) dissolved inorganic carbon concentration, (B) cell concentration and (C) acetate consumption in mixotrophic growth conditions of C. reinhartii in batch culture. — Figure 4: Time courses of (A) dissolved inorganic carbon concentration, (B) cell concentration and (C) acetate consumption in mixotrophic growth conditions of *C. reinhartii* in batch culture.

Download full-size image

DOI: 10.7717/peerj.5528/fig-4

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₂ 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₂) and the continuous input of CO₂ 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₂ as energy and carbon sources), heterotrophic (using acetate as carbon source), and mixotrophic (using both acetate and CO₂ 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₂ 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₂ 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₂ 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.

Figure 5: Phenotypic phase planes.
Response surfaces of growth rates predicted by FBA in function to (A) Heterotrophic growth: Flux acetate and flux of reaction R1612 (ATP production in mitochondria). (B) Photoautotrophic growth: Flux CO2 and flux reaction R3576 (ATP production in chloroplast).

Download full-size image

DOI: 10.7717/peerj.5528/fig-5

Table 2:

Growth rates predicted, and input fluxes of CO2 and acetate upon photoautotrophic, heterotrophic, and mixotrophic conditions.

Components	Photoautotrophic	Heterotrophic	Mixotrophic
CO₂ Flux (mmol/gDW h)	−12	85,73	−12
Acetate Flux (mmol/gDW h)	0	−10	−10
µmax (1/h)	0,313	0,522	0,835

DOI: 10.7717/peerj.5528/table-2

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₂ input variations as previously reported (Winck et al., 2016). Moreover, our results suggested that cells at high CO₂ (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₂ 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.

Figure 6: Profiles predicted by SOA of (A) biomass production and (B) acetate consumption over time.

Download full-size image

DOI: 10.7717/peerj.5528/fig-6

Figure 7: Reactions fluxes predicted associated to amino acids biosynthesis.
Valine, leucine and isoleucine biosynthesis: (ACAS 2ahbut) acetolactate synthase, (IMDH) 3- isopropylmalate dehydrogenase, (IMDHT) 3-isopropylmalate dehydratase, (IPMS) 2- isopropylmalate synthase, (KARI) ketol-acid reductoisomerase, (MOPC) 4-Methyl-2- oxopentanoate conversion. Phenylalanine, tyrosine and tryptophan biosynthesis: (ANTPPT) pyrophosphate phosphoribosyl-transferase, (ATNS nh4) Anthranilate synthase, (CHRM) Chorismate mutase, (CHRS) Chorismate synthase, (DHQS) 3-dehydroquinate synthase, (DPHS) 3-deoxy-7-phosphoheptulonate synthase, (DQDH) 3-deydroquinate dehydratase, (G3PL) D-glyceraldehyde-3-phosphate lyase, (IGPS) indole-3-glycerol-phosphate synthase. Lysine biosynthesis: (DAPAT) diaminopimelate aminotransferase, (DAPDC) diaminopimelate decarboxylase, (DAPEP) diaminopimelate epimerase, (DHDPS) dihydrodipicolinate synthase, (DHPR) dihydrodipicolinate reductase.

Download full-size image

DOI: 10.7717/peerj.5528/fig-7

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₂ 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₂ 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₂ 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₂ groups in culture (0%-heterotrophic, 10%-mixotrophic) considering the metabolic abundance. Metabolites and abundance were shown to be significantly affected by changes in CO₂ 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₂ concentrations, as confirmed by the experimental abundance of metabolites involved in these pathways. Our data suggested that high CO₂ 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₂ concentrations and low values for reaction fluxes at low CO₂ concentrations). This is in agreement with the experimental results of metabolic abundance that showed opposing behavior between the two CO₂ 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₂ 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₂ 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).

Figure 8: Experimental gene expression profiling.
(A) Carbon concentrating mechanism: Carbonic anhydrases (CAH1, CAH3, CAH4, CAH5, CAH7, CAH8, CAH9, CAH10, CAH11, CAH12), Carbon concentrating mechanism 1 (CCM1), Aldolase reductases (LCIA, LCI28) (B) Glycolysis/Gluconeogenesis and TCA cycle: Alcohol dehydrogenase/acetaldehyde dehydrogenase (ADH1), ATP synthases (ATP1A, ATP2) citrate synthase (CIS1), fructose-1,6-bisphosphatase (FBP1), glycerate kinase (GLYK), NAD-dependent malate dehydrogenase (MDH3), pyruvate decarboxylase (PDC3), phosphoglycerate kinase (PGK1), triose phosphate isomerase (TPIC).

Download full-size image

DOI: 10.7717/peerj.5528/fig-8

Figure 9: Reactions fluxes predicted for target genes associated with (A) carbon concentrating mechanism (B) photosynthesis, (C) glycolysis/gluconeogenesis and TCA cycle.
Alanine aminotransferase (AAT1), aspartate aminotransferase (AST1), fumarate hydratase (FUM1), glyceraldehyde-3-phosphate (GAP3), isocitrate dehydrogenase (IDH3), malate dehydrogenases (MME3, MDH3), phosphoglycerate kinases (PGK, PGKh, PGKm), phosphoenolpyruvate carboxykinase (PCK1), phosphoribulokinase (PRK1), pyruvate kinase (PYK1), ribulose-bisphosphate carboxylase (RBCS1), succinyl-CoA ligase ADPforming (SCL1), triosephosphate isomerase (TPI1), cytochrome b6 (PETM), cyclic Electron Flow (FDX1), ferredoxin reductase (FNR1), photosystem I (LHCA1), photosystem II (CPLD45), glyceraldehyde 3-phosphate dehydrogenase (GAP1), phosphoglucomutases (GPM1, GPM2), glucose-6-phosphate isomerase (PGI1), phosphoglycerate mutase (PGM2), succinate dehydrogenase ubiquinone (SDH2).

Download full-size image

DOI: 10.7717/peerj.5528/fig-9

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₂ 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₂-heterotrophic and high CO₂-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₂ 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.

Supplemental Information

Response surface of growth rates predicted by FBA in function of CO2 and acetate fluxes [mmol/gDW h]

DOI: 10.7717/peerj.5528/supp-1

Download

Canonical discriminant analysis of metabolomic data (A) heterotrophic growth condition (0% CO2) (B) mixotrophic growth condition (10% CO2)

DOI: 10.7717/peerj.5528/supp-2

Download

Canonical discriminant functions analysis of metabolomic data (A) grouping by time (B) grouping by growth condition and time

DOI: 10.7717/peerj.5528/supp-3

Download

New comprehensive genome-scale network reconstruction for C. reinhardtii, which consisted of 3726 reactions and 2436 metabolites

DOI: 10.7717/peerj.5528/supp-4

Download

New comprehensive genome-scale network reconstruction for C. reinhardtii

DOI: 10.7717/peerj.5528/supp-5

Download

Effects on different metabolic pathways, such as glycerolipid synthesis, amino acid synthesis, glycolysis, carbohydrate catabolism, and carbon fixation, among others

DOI: 10.7717/peerj.5528/supp-6

Download

[1] Anderson RA. 2005. Algal culturing techniques (1st edition). Amsterdam: Elsevier. 589

[2] Badger MR, Palmqvist K, Yu J-W. 1994. Measurement of CO₂ and HCO₃-fluxes in cyanobacteria and microalgae during steady-state photosynthesis. Physiologia Plantarum 90(3):529-536

[3] BioCyc Database Collection. 2016. Summary of Chlamydomonas reinhardtii. Version 20. 1 data

[4] Blazier AS, Papin JA. 2012. Integration of expression data in genome-scale metabolic network reconstructions. Frontiers in Physiology 3 Article 299

[5] Boyle N, Morgan J. 2009. Flux balance analysis of primary metabolism in Chlamydomonas reinhardtii. BMC Systems Biology 3(1):4

[6] Brueggeman AJ, Gangadharaiah DS, Cserhati MF, Casero D, Weeks DP, Ladunga I. 2012. Activation of the carbon concentrating mechanism by CO₂ deprivation coincides with massive transcriptional restructuring in Chlamydomonas reinhardtii. The Plant Cell 24:1860-1875

[7] Chaiboonchoe A, Ghamsari L, Dohai B, Ng P, Khraiwesh B, Jaiswal A, Jijakli K, Koussa J, Nelson DR, Cai H, Yang X, Chang RL, Papin J, Yu H, Balaji S, Salehi-Ashtiani K. 2016. Systems level analysis of the Chlamydomonas reinhardtii metabolic network reveals variability in evolutionary co-conservation. Molecular BioSystems 12:2394-2407

[8] Chang RL, Ghamsari L, Manichaikul A, Hom EF, Balaji S, Fu W, Shen Y, Hao T, Palsson BØ, Salehi-Ashtiani K, Papin JA. 2011. Metabolic network reconstruction of Chlamydomonas offers insight into light-driven algal metabolism. Molecular Systems Biology 7(518) Article 518

[9] Chen Q, Wang Z, Wei DQ. 2010. Progress in the applications of flux analysis of metabolic networks. Chinese Science Bulletin 55(22):2315-2322

[10] ChlamyCyc 5.0. 2015. Plant metabolic network. software

[11] Fang W, Si Y, Douglass S, Casero D, Merchant SS, Pellegrini M, Ladunga I, Liu P, Spalding MH. 2012. Transcriptome- wide changes in Chlamydomonas reinhardtii gene expression regulated by carbon dioxide and the CO₂-concentrating mechanism regulator CIA5/CCM1. The Plant Cell 24(5):1876-1893

[12] Gianchandani EP, Chavali AK, Papin JA. 2010. The application of flux balance analysis in systems biology. Wiley Interdisciplinary Reviews: Systems Biology and Medicine 2(3):372-382

[13] Harris EH. 2001. Chlamydomonas as a model organism. Annual Review of Plant Physiology and Plant Molecular Biology 52(1):363-406

[14] Ideker T, Krogan NJ. 2012. Differential network biology. Molecular Systems Biology 8:565

[15] Imam S, Schäuble S, Valenzuela J, Lõpez García De Lomana A, Carter W, Price NE, Baliga NS. 2015. A refined genome-scale reconstruction of Chlamydomonas metabolism provides a platform for systems-level analyses. The Plant Journal 84(6):1239-1256

[16] Kanehisa Laboratories. 2016. KEGG: Kyoto Encyclopedia of Genes and Genomes.

[17] Kauffman KJ, Prakash P, Edwards JS. 2003. Advances in flux balance analysis. Current Opinion in Biotechnology 14(5):491-496

[18] Kebeish R, Niessen M, Thiruveedhi K, Bari R, Hirsch H, Rosenkranz R, Stäbler N, Schönfeld B, Kreuzaler F, Peterhänsel C. 2007. Chloroplastic photorespiratory bypass increases photosynthesis and biomass production in Arabidopsis thaliana. Nature Biotechnology 25:593-599

[19] Kliphuis AMJ, Klok AJ, Martens DE, Lamers PP, Janssen M, Wijffels RH. 2012. Metabolic modeling of Chlamydomonas reinhardtii: energy requirements for photoautotrophic growth and maintenance. Journal of Applied Phycology 24(2):253-266

[20] Kumar VS, Dasika M, Maranas CD. 2007. Optimization based automated curation of metabolic reconstructions. BMC Bioinformatics 8:212

[21] Lee E, Jalalizadeh M, Zhang Q. 2015. Growth kinetic models for microalgae cultivation: a review. Algal Res 12:497-512

[22] Mahadevan R, Edwards JS, Doyle FJ. 2002. Dynamic flux balance analysis of diauxic growth in Escherichia coli. Biophysical Journal 83(3):1331-1340

[23] Mahadevan R, Schilling CH. 2003. The effects of alternate optimal solutions in constraint-based genome-scale metabolic models. Metabolic Engineering 5(4):264-276

[24] Maranas C, Zomorrodi A. 2016. Optimization methods in metabolic networks. Hoboken: Wiley.

[25] Mavrovouniotis ML. 1990. Group contributions for estimating standard Gibbs energies of formation of biochemical compounds in aqueous solution. Biotechnology and Bioengineering 36(10):1070-1082

[26] ModelSEED. 2015. Model seed genome-scale metabolic models.

[27] Moon M, Kim CW, Park WK, Yoo G, Choi YE, Yang JW. 2013. MixotrophicGrowth with acetate or volatile fatty acids maximizes growth and lipid production in Chlamydomonas reinhardtii. Algal Research 2(4):352-357

[28] Moroney JV, Somanchi A. 1999. How do algae concentrate CO₂ to increase the efficiency of photosynthetic carbon fixation? Plant Physiology 119(1):9-16

[29] Moroney JV, Ynalvez RA. 2007. Proposed carbon dioxide concentrating mechanism in Chlamydomonas reinhardtii. Eukaryotic Cell 6(8):1251-1259

[30] Petersen TN, Brunak S, Von Heijne G, Nielsen H. 2011. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nature Methods 8(10):785-786

[31] Rosenberg JN, Oyler GA, Wilkinson L, Betenbaugh MJ. 2008. A green light for engineered algae: redirecting metabolism to fuel a biotechnology revolution. Current Opinion in Biotechnology 19(5):430-436

[32] Rupprecht J. 2009. From systems biology to fuel—Chlamydomonas reinhardtii as a model for a systems biology approach to improve biohydrogen production. Journal of Biotechnology 142(1):10-20

[33] Sayre R. 2010. Microalgae: the potential for carbon capture. BioScience 60(9):722-727

[34] Schellenberger J, Lewis NE, Palsson BO. 2011. Elimination of thermodynamically infeasible loops in steady-state metabolic models. Biophysical Journal 100:544-553

[35] Schmidt BJ, Lin-Schmidt X, Chamberlin A, Salehi-Ashtiani K, Papin JA. 2010. Metabolic systems analysis to advance algal biotechnology. Biotechnol. J 5(7):660-670

[36] Schmittgen TD, Livak KJ. 2008. Analyzing real-time PCR data by the comparative CT method. Nature Protocols 3(6):1101-1108

[37] Scranton MA, Ostrand JT, Fields FJ, Mayfield SP. 2015. Chlamydomonas as a model for biofuels and bio-products production. The Plant Journal 82(3):523-531

[38] Senger RS, Papoutsakis ET. 2008. Genome-scale model for Clostridium acetobutylicum: Part I. Metabolic network resolution and analysis. Biotechnology and Bioengineering 101(5):1036-1052

[39] Serrano-Bermúdez LM. 2016. Análisis de balance de flujo dinámico de la producción de 1, 3—propanodiol a partir de Clostridium sp. Bogota: Universidad Nacional de Colombia.

[40] Serrano-Bermúdez LM, González Barrios AF, Maranas CD, Montoya D. 2017. Clostridium butyricum maximizes growth while minimizing enzyme usage and atp production: metabolic flux distribution of a strain cultured in glycerol. BMC Systems Biology 11(1):58

[41] Tebbani S, Titica M, Caraman S, Boillereaux L. 2013. Estimation of Chlamydomonas reinhardtii growth in a torus photobioreactor. In: 12th IFAC symposium on computer applications in biotechnology, 16–18, December 2013. The International Federation of Automatic Control. Mumbai.