Metabolic Engineering Communications 9 (2019) e00103 Contents lists available at ScienceDirect Metabolic Engineering Communications journal homepage: www.elsevier.com/locate/mec Contextualized genome-scale model unveils high-order metabolic effects of the specific growth rate and oxygenation level in recombinant Pichia pastoris Paulina Torres a, Pedro A. Saa a, Joan Albiol b, Pau Ferrer b, Eduardo Agosin a, * a Department of Chemical and Bioprocess Engineering, School of Engineering, Pontificia Universidad Catolica de Chile, Vicu~ na Mackenna, 4860, Santiago, Chile b Department of Chemical, Biological, and Environmental Engineering, Universitat Autonoma de Barcelona, 08193, Bellaterra (Cerdanyola del Valles), Catalonia, Spain A R T I C L E I N F O A B S T R A C T Keywords: Pichia pastoris is recognized as a biotechnological workhorse for recombinant protein expression. The metabolic Pichia pastoris performance of this microorganism depends on genetic makeup and culture conditions, amongst which the Experimental design specific growth rate and oxygenation level are critical. Despite their importance, only their individual effects have Metabolic modelling been assessed so far, and thus their combined effects and metabolic consequences still remain to be elucidated. In Dissolved oxygen Recombinant protein this work, we present a comprehensive framework for revealing high-order (i.e., individual and combined) Thaumatin metabolic effects of the above parameters in glucose-limited continuous cultures of P. pastoris, using thaumatin production as a case study. Specifically, we employed a rational experimental design to calculate statistically significant metabolic effects from multiple chemostat data, which were later contextualized using a refined and highly predictive genome-scale metabolic model of this yeast under the simulated conditions. Our results revealed a negative effect of the oxygenation on the specific product formation rate (thaumatin), and a positive effect on the biomass yield. Notably, we identified a novel positive combined effect of both the specific growth rate and oxygenation level on the specific product formation rate. Finally, model predictions indicated an opposite rela- tionship between the oxygenation level and the growth-associated maintenance energy (GAME) requirement, suggesting a linear GAME decrease of 0.56 mmol ATP/gDCW per each 1% increase in oxygenation level, which translated into a 44% higher metabolic cost under low oxygenation compared to high oxygenation. Overall, this work provides a systematic framework for mapping high-order metabolic effects of different culture parameters on the performance of a microbial cell factory. Particularly in this case, it provided valuable insights about optimal operational conditions for protein production in P. pastoris. 1. Introduction liter titers (Ciplys et al., 2015; Hasslacher et al., 1997; Heyland et al., 2010; Wang et al., 2001). Despite these achievements, efficient heter- The methylotrophic yeast Pichia pastoris (syn. Komagataella phaffii) is ologous protein production in P. pastoris remains challenging, as poorly a popular cell factory for heterologous protein production due to its tuned protein overexpression can affect relevant cellular processes, such convenient physiology and the availability of genetic tools for its as protein folding and secretion (Delic et al., 2014; Gasser et al., 2007; manipulation (Gasser and Mattanovich, 2018; Pe~ na et al., 2018; Theron Love et al., 2012). Moreover, codon usage level (Hu et al., 2013; Xiang et al., 2018). Amongst its features, P. pastoris displays a et al., 2016), promoter selection (Prielhofer et al., 2013), as well as Crabtree-negative growth phenotype, leading to reduced synthesis of culture medium composition (Heyland et al., 2011) and operational by-products under aerobic conditions (Baumann et al., 2010; Calik conditions (Cos et al., 2006; Maurer et al., 2006) may also play major et al., 2015). In addition, the availability of strong promoters for het- roles on process performance. In particular, the operational conditions erologous expression, added to its natural capacity to perform have gained increasing attention as they are known to introduce sub- post-translational modifications, has made this yeast the platform of stantial variability in the process, significantly affecting the recombi- choice for producing a range of recombinant proteins up to grams per nant protein secretion (Looser et al., 2015). Abbreviations: GSMM, genome-scale metabolic model; DO, dissolved oxygen; GAME, growth-associated maintenance energy; NGAME, non-growth-associated maintenance energy; FBA, flux balance analysis; ll-FBA, loopless flux balance analysis; FVA, flux variability analysis. * Corresponding author. E-mail address: agosin@ing.puc.cl (E. Agosin). https://doi.org/10.1016/j.mec.2019.e00103 Received 15 July 2019; Received in revised form 4 October 2019; Accepted 8 October 2019 2214-0301/© 2019 The Authors. Published by Elsevier B.V. on behalf of International Metabolic Engineering Society. This is an open access article under the CC BY- NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Table 1 Primers used in this study. Primer Sequence Thau-pGAPZB-F CTATTTCGAAACGAGGAATTCACGATGGCCGCTACCACCTGC Thau-pGAPZB-R GAACAGTCATGTCTAAGGCTAAAACTCACTCGTCTTCCAACTCCAAGGCTGTA pGAPZB-Thau-F TACAGCCTTGGAGTTGGAAGACGAGTGAGTTTTAGCCTTAGACATGACTGTTC pGAPZB-Thau-R GCAGGTGGTAGCGGCCATCGTGAATTCCTCGTTTCGAAATAG High recombinant protein expression in P. pastoris relies on the use 2. Materials and Methods of strong promoters, like pAOX1 (promoter from alcohol oxidase I encoding gene) and pGAP (promoter from glyceraldehyde-3-phosphate 2.1. Plasmid construction and strain transformation dehydrogenase encoding gene). While pAOX1 offers strong inducible expression with methanol – thereby enabling uncoupling fast growth The thaumatin gene – including its natural pre-region secretion signal from production –, pGAP provides comparable constitutive expression – was synthesized by Genscript (Piscataway, NJ, USA) and was codon- (Pe~ na et al., 2018). P. pastoris cultures incur in high oxygen consump- optimized for expression in P. pastoris. The commercial recombinant tion and heat production during methanol oxidation, and hence, its use vectors pGAPZB (Invitrogen, Carlsbad, CA, USA) were employed for the poses major challenges for large-scale protein production (Mattanovich construction of the transformation plasmids using Gibson assembly et al., 2014). Once a suitable expression system has been chosen, the (Gibson et al., 2009), and later used to integrate the thaumatin gene into next step is to optimize culture conditions to achieve the target pro- the yeast’s genome under the constitutive GAP promoter. The sequences ductivity. Factors such as temperature, pH, osmolality, specific growth of the employed primers are shown in Table 1. Amplification of DNA rate (μ) and dissolved oxygen (DO) are critical for the efficient operation fragments for Gibson assembly was carried out by PCR using Phusion of the culture, and their influence on protein production and culture High-Fidelity DNA Polymerase (ThermoFisher, Waltham, MA, USA). All performance has been individually assessed (Baumann et al., 2008; PCR products were treated with DpnI restriction enzyme to remove Charoenrat et al., 2005; Dragosits et al., 2009, 2010; Garcia-Ortega original vector residues and later were purified by gel extraction, using et al., 2017; Heyland et al., 2010; Maurer et al., 2006). Although several the Qiaquick Gel Extraction kit (Qiagen, Hilden, Germany). Finally, studies have reviewed the relationships between protein production and chemically competent E. coli TOP 10 cells were transformed using the growth (refer to Looser et al. (2015) for a comprehensive review), and sought construct. These cells were grown at 37  C in low salt-LB medium, how DO impacts the yeast’s physiology (Adelantado et al., 2017; Bau- containing 25 μg mL1 zeocin for selection of clones transformed with mann et al., 2010; Garcia-Ortega et al., 2017), current studies fail to pGAPZB-TAU vector. evaluate both the individual and combined (high-order) effects of these P. pastoris wild-type strain GS115 (Invitrogen, Carlsbad, CA, USA) operational parameters on the metabolic performance of P. pastoris. was used as a host strain throughout this study, which was transformed Notably, such evaluation is not trivial as both, an appropriate experi- using an in-house-built vector to revert its histidine auxotrophy (refer to mental design and modelling framework, are required to appropriately Supplementary Text S1). AvrII was employed to linearize the trans- interpreting the results. formation vector, which was introduced by electroporation into the In this study, we present an integrated framework for mapping, at the competent cells, as described by Gasser et al. (2006). Both plasmids and genome–scale, the individual and combined metabolic effects of key transformations were verified by DNA sequencing (Macrogen Inc., Seoul, operational parameters – here μ and DO – during recombinant protein Korea). production by P. pastoris, under glucose-limited conditions in continuous cultures. As a case study, we analyzed the metabolic behavior of a re- 2.2. Cell cultivation combinant strain producing the sweet-tasting, low-calorie protein thau- matin. This protein has 207 amino acid residues and 8 disulfide bonds Continuous cultures were started from pre-inocula grown overnight (Illingworth et al., 1989), which are critical for its sweet taste (Masuda at 30  C and 150 rpm in 200-mL shake flasks, containing YPG medium et al., 2016) and are considered the main reason behind the low titers with 100 μg mL1 zeocin. Prior to the inoculation of the bioreactors, each achieved so far (Moralejo et al., 2001) (~ 100 mg L1 in high-density cell inoculum was centrifuged at 5000 rpm for 5 min and resuspended in cultures (Masuda et al., 2010)). Folding of recombinant proteins with fresh culture medium without trace elements. Chemostat cultures were many disulfide bounds is both difficult and costly, as it requires a high performed in 2-L benchtop Biostat B bioreactors (Sartorius AG, supply of NAD(P)H cofactors which can affect redox homeostasis and G€ottingen, Germany) using 1 L of working volume. lead to negative physiological responses like the Unfolded Protein P. pastoris was grown in glucose-limited chemostats under the culture Response (UPR) and Endoplasmic-Reticulum-Associated Degradation conditions described by Adelantado et al. (2017) with only minor dif- (ERAD) (Gasser et al., 2007; Puxbaum et al., 2015). Thus, understanding ferences (see Supplementary Text S2). For DO control, a mixture of air, the effects of μ and DO that have a major metabolic impact is critical for pure nitrogen and pure oxygen was used. Each gas was independently optimizing heterologous protein production in P. pastoris. To capture the controlled using three mass flow meters (Bronkhorst High-Tech, Ruurlo, individual and combined effects of μ and DO on the metabolic perfor- Netherlands) until the desired DO levels were reached (see below). The mance of this strain, we implemented a fractional experimental design pH, agitation, and temperature were set respectively to 5.0, 700 rpm, and to determine the statistically significant effects from a reduced set of 25  C in all conditions. Samples were taken after four and five residence experiments, and later integrated these with a tailored Genome-Scale times, when the culture reached steady state. Metabolic Model (GSMM) to interpret and further analyze the meta- bolic and energetic consequences. GSMMs describe the set of 2.3. Experimental design biochemical reactions that represent the entire metabolic potential of an organism, and as such, are excellent tools to understand the physio- Experimental design is a statistical tool that is used to systematically logical state of the cell under different environmental conditions (Price examine the effects of two or more factors that shape a determined et al., 2004). Our results revealed opposite operational regions for response function. Among the different possible designs, Doehlert design efficient growth and protein production, as well as a tighter energetic displays higher efficiency and flexibility, as it allows fitting complex constraint on P. pastoris growth under glucose-limited, low DO response models with less data than complete designs, and without conditions. requiring extreme points difficult to acquire (Doehlert, 1970; Hanrahan 2 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 h1 for μ with five dilution rate levels; 0.05, 0.075, 0.1, 0.125 and 0.15 h1 and over 0–60% DO saturation with three oxygenation levels; low (4%, L condition), medium (30%, M condition) and excess (56%, E condition) (Fig. 1A). The latter DO levels correspond to approximately 10.4 (L), 77.7 (M) and 145.1 μM (E) DO concentrations assuming a Henry’s constant of 957.7 atm M1 at 25  C for a fermentation broth (Saa et al., 2012). This design samples seven experimental conditions in total: six external conditions performed in duplicate and in random order, and a central point performed in triplicate. Finally, a second-degree poly- nomial (Eq. (1)) was fitted to determine the response of the dependent variable as a function of the different factors, i.e. μ and DO. Z ¼ b0 þ b1 X1 þ b2 X2 þ b11 X 21 þ b12 X1 X2 þ b22 X 22 (1) In Eq. (1), Z denotes the predicted response, X1 and X2 represent the independent variables, b0 is the mean value of the response, b1 and b2 are linear (main) effects, b11 and b22 represent squared effects, and b12 de- scribes an interaction term (combined effect). Statistical calculations and significant tests (ANOVA) were performed using STATGRAPHICS Plus 5.1 (Statpoint Technologies, Inc., VA, USA). 2.4. Extracellular metabolite and biomass quantification Extracellular glucose, arabitol and ethanol were analyzed by HPLC. Samples of 2 mL each were taken and centrifuged at 10,000 rpm and 4  C for 5 min. The supernatant was collected and filtered using 0.45 μm-pore nitrocellulose filters (Merck Millipore, Carrigtwohill, Ireland). HPLC analysis was performed in duplicate in a HP 1050 liquid chromatography system (Dionex Corporation, Sunnyvale, USA) as described by Tomas-- Gamisans et al. (2016). Extracellular thaumatin was quantified in trip- licate using Thaumatin ELISA Kit (CellTrend, Luckenwalde, Germany), according to manufacturer instructions. Finally, Dry Cell Weight (DCW) was quantified in triplicate, following the method of Jorda et al. (2012). 2.5. Contextualization and refinement of the genome-scale metabolic model for accurate prediction of glucose-limited oxygen-sufficient growth Several GSMMs have been made available for modelling P. pastoris metabolism and optimizing heterologous protein production under different conditions (Caspeta et al., 2012; Irani et al., 2016; Sohn et al., 2010; Tomas-Gamisans et al., 2018; Ye et al., 2017). These models are typically validated by comparing experimental and predicted specific growth and exchange rates obtained using Flux Balance Analysis (FBA) (Orth et al., 2010). However, this is often insufficient to validate GSMMs, as identical FBA predictions may be underpinned by completely different intracellular flux distributions. As such, consideration of available intracellular flux data (e.g., measured by 13C fluxomics), is an important (often neglected) refinement step for accurately capturing the metabolic behavior for a simulated phenotype. An excellent example of the importance of the above has been recently presented by Pereira et al. Fig. 1. Overview of the experimental design and GSMM contextualization (2016). workflow employed in this study. A) Doehlert fractional experimental design for In order to achieve accurate and meaningful predictions underpinned evaluating the metabolic impact of the specific growth rate (μ) and dissolved by realistic flux distributions, we set up a sequential refinement workflow oxygen (DO) level on the metabolic performance of P. pastoris. This design starting from the latest and most curated P. pastoris GSMM – the iMT1026 considers 1 central point and 6 extreme points over 3 DO levels (4%, 30% and v3.0 (Tomas-Gamisans et al., 2018) –, where intracellular reactions were 56% oxygen saturation), and 5 μ levels (0.05, 0.075, 0.1, 0.125 and 0.15 h1). incrementally constrained (both activities and directionalities), based on B) Manual GSMM contextualization workflow for flux simulations. The latest P. pastoris GSMM (iMT1026 v3.0) was incrementally constrained for accurately reported fluxomics, proteomics and transcriptomics data under appro- describing metabolic fluxes under glucose-limited conditions using reported priate experimental conditions (Clasquin et al., 2011; Gasser et al., 2007; experimental data and flux feasibility criteria. Rußmayer et al., 2015; Zhang et al., 2017), and the imposition of the loopless flux condition (elimination of thermodynamically infeasible and Lu, 2006). closed cycles) (Saa and Nielsen, 2016a, 2016b) (Fig. 1B). In this way, the In this study, the Doehlert matrix (Doehlert, 1970) was constructed to refined metabolic model generated realistic intracellular fluxes, i.e. evaluate the individual and combined effects of μ and DO on the macro- consistent with reported metabolic pathway utilization behaviors and scopic growth and production parameters of the thaumatin-producing flux patterns and, most importantly, accurately predicted observed spe- P. pastoris strain. The experimental domain was defined over 0.05–0.15 cific growth and exchange rates. Additional details about the contextu- alization and refinement workflow can be found in Supplementary Text 3 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Table 2 Metabolic performance parameters of recombinant P. pastoris in glucose-limited chemostats under different μ and DO conditions.a. Conditionb μ (h1) qS (mmol g1 1 DCW h ) qCO2 (mmol g1 1 DCW h ) qO2 (mmol g1 1 DCW h ) qP (mg g1 1 DCW h ) YS,X (gDCW g1) YS,P (mg g1) 1 (L) 0.072  0.003 0.79  0.038 1.84  0.101 1.84  0.097 0.05  0.005 0.51  0.002 0.27  0.042 2 (L) 0.119  0.006 1.25  0.031 2.77  0.017 2.79  0.317 0.04  0.003 0.53  0.011 0.16  0.016 3 (M) 0.048  0.002 0.48  0.029 1.03  0.051 1.05  0.068 0.03  0.001 0.55  0.008 0.40  0.019 4 (M) 0.101  0.001 0.99  0.009 1.98  0.016 1.97  0.076 0.03  0.001 0.57  0.002 0.19  0.001 5 (M) 0.150  0.002 1.50  0.021 3.02  0.052 2.96  0.099 0.035  0.001 0.57  0.002 0.13  0.002 6 (E) 0.072  0.001 0.68  0.007 1.21  0.004 1.30  0.005 0.03  0.001 0.59  0.004 0.26  0.007 7 (E) 0.120  0.001 1.14  0.014 2.18  0.147 2.31  0.129 0.03  0.001 0.58  0.004 0.18  0.009 a Results represent the average and standard error of two or three independent experiments (biological replicates) for the extreme and the central points, respectively. In each case, the carbon balance closed satisfactorily and reconciled no less than 94.7% of the consumed carbon source (see Supplementary Table S5). b The experimental conditions are abbreviated L for low oxygenation (4% dissolved oxygen saturation), M for medium oxygenation (30% dissolved oxygen satu- ration), and E for excess oxygenation (56% dissolved oxygen saturation). S3 and Tables S1–S3. The contextualized constrained metabolic model is available in Supplementary File S1 and has been deposited in the Bio- Models database (Chelliah et al., 2015) under MODEL1910040001 identifier. A MATLAB script for generating the contextualized models is also available in Supplementary File S1. 2.6. Estimation of energetic requirements for growth GSMMs have cellular ATP requirement divided into Growth- Associated Maintenance Energy (GAME) and Non-Growth-Associated Maintenance Energy (NGAME). The former (GAME) represents the (milli)moles of ATP required per gram of biomass (YX,ATP) – represented by the corresponding stoichiometric coefficient in the biomass equation; whereas the latter (NGAME) describes the ATP consumption (mATP) required for cellular maintenance (Pirt, 1982). These parameters together quantify the total ATP requirement per gram of dry cell weight per hour (qATP) for growing under a defined condition and follow the relation shown below. qATP ¼ YX;ATP  μ þ mATP (2) So far, the energetic impact of the oxygenation on P. pastoris meta- bolism has not been fully determined, as reported experimental designs are inadequate to properly assess its effects, e.g. chemostat data has only been reported at extreme DO conditions and at a single μ of 0.1 h1 (Adelantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., 2017). In contrast, the experimental design proposed here enables ex- amination of the energetic parameters that best fit the experimental data (growth and production rates at different DO and μ conditions), shedding light on the energetic impact of the DO level. To determine the Fig. 2. Response surfaces of P. pastoris macroscopic culture parameters under latter, the qATP of the refined GSMM was maximized in all the experi- different DO and μ conditions. The fitted response surfaces were for the mental conditions using the loopless FBA (ll-FBA) algorithm (Saa and (encoded) biomass yield (A), specific thaumatin production rate (B) and thau- Nielsen, 2016a) – i.e., a flux estimator that avoids thermodynamically matin yield (C) are, respectively: YS;X ¼ 0:565 þ 0:036  DO  0:02  μ  DO  infeasible loops. In all the simulations, the macromolecular biomass 0:017  DO2 (R2 ¼ 0.948), qP ¼ 0:033  0:004  DO þ 0:008  μ  DO (R2 ¼ 0.68), composition employed was based on available data from glucose-limited and YS;P ¼ 0:186  0:127  μ þ 0:08  μ2 (R2 ¼ 0.959), where μ ¼ 0:05 þ 0:05  cultures (Carnicer et al., 2009), which is available in iMT1026 v3.0. For ðμ þ1Þ ; μ 2 ½ 1; 1 and DO ¼ 4 þ 30  ðDO þ 0:8666Þ ; DO 2 ½  0:8666; the qATP estimation, the mATP reaction was chosen as the objective, the 0:8666. Further details about the regression and statistical effects can be found GAME was set to zero, and the growth and exchange rates were fixed to in Supplementary Table S6. the observed values. Once qATP was computed, the GAME was fitted to the different DO conditions using Eq. (2) after setting NGAME to the 3. Results most recent measured value (Rebnegger et al., 2016). Another tested model considered fitting NGAME for different DO conditions and setting 3.1. Effects of the specific growth rate and oxygenation level on the GAME to the GSMM value of 72 mmol ATP g1 DCW (Tom as-Gamisans metabolic performance of P. pastoris under glucose-limited conditions et al., 2016). At this point, three models were hypothesized and fitted: one where YX,ATP and NGAME were the same for all DO conditions; Growth and metabolic production parameters of the recombinant another, where only YX,ATP was fitted to each DO condition; and a final P. pastoris strain were measured in seven chemostat cultures (conditions), model, where only NGAME was fitted to each DO condition. Finally, the under different μ and oxygenation conditions. In each case, biomass and Akaike’s Information Criterion for low number of data points (AICc) was thaumatin yields (YS,X and YS,P, respectively), specific thaumatin and applied to compare the goodness-of-fit of each model (Hurvich and Tsai, carbon dioxide production rates (qP and qCO2, respectively), and specific 1989). The model with the smallest AICc was considered the best glucose and oxygen consumption rates (qS and qO2, respectively) were choice. estimated (Table 2). In all conditions, ethanol and arabitol were not detected, pointing to a purely respiratory metabolism, which was later 4 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Fig. 3. Evaluation of contextualized P. pastoris GSMM for metabolic flux prediction under glucose-limited conditions. A) Step-by-step contextualization process of the P. pastoris iMT1026 v3.0 GSMM for describing growth under glucose-limited conditions. B) Evaluation of maximum specific growth rate predictive performance of the initial versus contextualized GSMM. C) Comparison of the predicted intracellular flux through cytoplasmic reactions of the initial (slope ¼ 0.911, R2 ¼ 0.783) and contextualized GSMM (slope ¼ 1.047, R2 ¼ 0.972) against experimental data under appropriate growth conditions. D) Comparison of the predicted intracellular flux through mitochondrial reactions of the initial (slope ¼ 0.112, R2 ¼ 0.076) and contextualized GSMM (slope ¼ 0.721, R2 ¼ 0.867) against experimental data under the same conditions used in C). For more information about the employed experimental conditions refer to Materials and Methods. confirmed with measured RQ values (Supplementary Table S5). sequential refinement workflow. In a first step, the reaction di- Oxygenation exerted the greatest impact on both, YS,X and qP, rectionalities and activities were modified based on literature data. At whereas μ had the main effect on YS,P. In the case of YS,X, there were this point, it is important to note that iMT1026 v3.0 supports growth on notable differences between different oxygenation levels at a given μ glycerol and methanol through the introduction of accessory pathways in (e.g., L and E conditions showed a 14% relative difference for μ ¼ 0.075 different compartments (e.g., peroxisomal pentose phosphate pathway, h1), but not vice versa (< 4.7% relative difference for a fixed oxygen- PPP) (Tomas-Gamisans et al., 2018). The activity of such pathways is, ation level and different μ). Overall, YS,X varied from 0.507 to 0.589 gDCW however, absent under glucose-limited conditions, as a previous report g1, attaining the lowest yield at μ ¼ 0.075 h1 in the L condition. Sta- suggests (Rußmayer et al., 2015). Fail to block these pathways led to the tistical analysis showed that DO had the main positive and significant simulation of unrealistically high glycolytic fluxes (Supplementary effect (p-value<0.001, Supplementary Table S6), although the interac- Fig. S1). Also, the cytosolic NADP-dependent isocitrate dehydrogenase tive and quadratic terms displayed also significant (both negative) in- was not blocked, as previously suggested (Pereira et al., 2016; Saitua fluences (p-value<0.01 and p-value<0.05, respectively; refer to et al., 2017), as there is no reported evidence for the lack of activity in Supplementary Table S6 for details). The YS,X regression explained a P. pastoris. Another modification to the model considered fixing the substantial part of the variability in the data (R2 ¼ 0.948) and, more sedoheptulose 1,7-bisphosphate D-glyceraldehyde-3-phosphate-lyase importantly, it suggested high oxygenation and low μ as the most reaction in the direction of riboneogenesis (Clasquin et al., 2011), which convenient conditions for efficient growth (Fig. 2A). improved considerably the flux predictions in the non-oxidative branch In the case of the thaumatin production parameters qP and YS,P, they of the PPP, as compared to above 13C labeling flux data. Cytosolic showed different behaviors. While YS,P was almost unaffected by DO and acetyl-CoA synthesis was forced to follow the canonical pathway through displayed a negative correlation with μ (e.g., varying from 0.129 to pyruvate by blocking the methylmalonate-semialdehyde dehydrogenase. 0.399 mg thaumatin g1 in the high and low μ conditions, respectively), The directionality of the mitochondrial malate dehydrogenase was fixed qP was largely influenced by DO (negatively); and to a lesser extent by μ in the direction of malate oxidation, consistent with an inactive glyox- (positively) (Table 2 and Fig. 2B). In this case, the regression analysis ylate cycle under various carbon sources (Sola et al., 2004, 2007) and explained a lower fraction of the variability of the data (R2 ¼ 0.68) and more recent 13C flux data in aerobic, glucose-limited cultures (Sola et al., indicated a significant negative effect of DO (p-value<0.05; Supple- 2004; Zhang et al., 2017). Finally, the activity of the malic enzyme was mentary Table S6) with a positive interaction term (p-value<0.05, Sup- blocked as it carries a neglectable flux under the hypoxic and normoxic plementary Table S6). Here, the magnitudes of both, the individual and conditions in glucose-limited cell cultures (Baumann et al., 2010; Zhang combined effects, were quite similar (see Supplementary Table S6). In the et al., 2017). Altogether, correction of these reactions substantially case of the YS,P, significant individual (negative) and quadratic (positive) improved mitochondrial flux predictions, as compared to available 13C effects of μ were suggested by the regression analysis (in both cases p- flux measurements (see below). value<0.001, and R2 ¼ 0.959, see Supplementary Table S6). Altogether, Integration of transcriptomic and proteomic data under appropriate regions of both, low μ and oxygenation, favor thaumatin production and experimental conditions (Clasquin et al., 2011; Gasser et al., 2007; yield (Fig. 2B–C). Rußmayer et al., 2015; Zhang et al., 2017), and directionalities based on reaction thermodynamics, reduced the number of reversible reactions by 25 and blocked 23 initially active reactions (Fig. 3A). Furthermore, 3.2. GSMM contextualization and refinement subsequent thermodynamic feasibility analysis, based on imposition of the loopless flux condition, yielded irreversible 215 reactions previously The prediction power of iMT1026 v3.0 was improved through a 5 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 tighter as compared to the initial iMT1026 v3.0 (Supplementary Fig. S3). The final list with constraints employed can be found in Supplementary Table S3. 3.3. GAME evaluation under oxygen-sufficient, glucose-limited conditions The contextualized iMT1026 v3.0 was employed to evaluate the en- ergetic effects of the DO level on the P. pastoris GAME. To this end, the GAME and NGAME parameters were refitted to the experimental data of this study and allowed to vary depending on the DO level (refer to Esti- mation of energetic requirements for growth for details). Our results showed that the model with a variable GAME as a function of the DO level ex- plains a higher amount of the data variability (94.9% compared to 87.1% and 88.6%), and has a stronger statistical support, compared to the alternative (Fig. 4A). Notably, this model displayed a decreasing GAME with DO level, starting from 95 mmol ATP g1 DCW for the L condition, and descending to 75 and 66 mmol ATP g1 DCW for the M and E conditions, respectively (44% higher cost between L and E) (Fig. 4B). These results indicate a decreasing growth cost with increasing DO level. 3.4. Simulation of intracellular fluxes under extreme DO levels As the DO level had the greatest effect on the metabolic parameters of P. pastoris, the intracellular flux distributions under extreme DO condi- tions were simulated at a fixed μ level (μ ¼ 0.075 h1). The absolute glycolytic flux increased by approximately 15% in the L condition, although no substantial carbon redistribution was predicted in cytosolic fluxes, as observed by the almost identical glycolysis, PPP and biosyn- thetic relative fluxes (Fig. 5A). Both the TCA cycle and ATP synthase relative fluxes increased proportionally with qCO2 and qO2 (Fig. 5A and Table 2). Although there is also an important difference between the pyruvate export between both conditions, there is no compensation by other by-products, i.e. ethanol or arabitol (Fig. 5A). In both conditions, the glyoxylate cycle is inactive consistent with a canonical operation of the TCA cycle. Surprisingly, the absolute TCA flux under the L condition is approx. 42% greater than in the E condition, in agreement with a greater ATP production (Fig. 5A). Interestingly, analysis of the feasible flux ranges revealed slightly tighter metabolic constraints under the E condition, despite the increased energetic growth cost of the L condition (i.e., higher YX,ATP). These results are, however, consistent with the Fig. 4. Evaluation of the DO influence on the GAME and NGAME requirements of P. pastoris under glucose-limited conditions. A) Comparison of the specific increased glucose uptake rate in the L condition, which provides more growth rate predictions of the model with a fixed GAME for all DO conditions flexibility for carbon allocation (Fig. 5B). (slope ¼ 0.752, R2 ¼ 0.871), the model with a variable GAME for each DO con- dition (slope ¼ 0.844, R2 ¼ 0.949), and the model with a variable NGAME for 4. Discussion each DO condition (slope ¼ 0.890, R2 ¼ 0.886) against the experimental growth data of this study. The model with a variable GAME for each DO condition has 4.1. Oxygenation exerts an individual negative effect and a synergistic stronger statistical support as shown by the lower AICc value (-158.9 versus positive effect with the specific growth rate on thaumatin production -146.8 and -154.7). B) Predicted GAME values for different DO conditions (95 mmol ATP g1 1 DCW for the L condition, 75 mmol ATP gDCW for the M condition, The specific thaumatin production rate was influenced negatively by and 66 mmol ATP g1 DCW for the E conditions). the DO level, and positively by μ through the interaction with DO (Fig. 2). The latter effect is consistent with previous reports of a positive corre- deemed reversible (Fig. 3A). Altogether, these measures disabled infea- lation between μ and protein secretion rates at fixed medium DO level, sible cycles responsible for the flux prediction inaccuracies. Next, although the estimated effect was not so strong as previously reported adjustment of the NGAME to 0.55 mmol ATP g1 DCW h 1 (assuming a yield (Buchetics et al., 2011; Maurer et al., 2006). Recent studies have showed of 32 mol of ATP per mol of glucose (Tomas-Gamisans et al., 2016), and that specific protein production in P. pastoris can also depend on the using the recent experimental maintenance value of 3.1 mg glucose g1 DCW nature of the produced protein, expression system and oxidative condi- h1 measured in glucose-limited retentostats (Rebnegger et al., 2016), tions (Gasser et al., 2006; Liu et al., 2013; Maccani et al., 2014; Uchiyama increased the prediction fidelity of iMT1026 v3.0 by ca. 86% for the and Shioya, 1999). Indeed, despite having constitutive protein expres- observed μ – from approx. 7.0% average relative error to 4.2% (Fig. 3B). sion, non-growth-associated protein production patterns may arise due to Moreover, the agreement between predicted and measured intracellular global stress responses, high protein turnover and limitations in endo- flux data improved considerably – explained flux variability increased plasmic reticulum processing, among others (Liu et al., 2013). The latter from 78.3% to 97.2% in the cytosol (Fig. 3C) and from 7.6% to 86.7% in points to protein expression, processing and/or secretion as likely – and the mitochondria (Fig. 3D). We note that iMT1026 v3.0 uses a biomass increasingly more recognized – bottlenecks (Looser et al., 2015). macromolecular composition with ashes, although correction by the The calculated individual DO effect on protein production agrees well estimated ashes content yielded almost identical results (Supplementary with previous studies using the pGAP-based expression system, reporting Fig. S2). Finally, the flux ranges of the contextualized model were slightly higher protein secretion under hypoxic conditions at μ ¼ 0.1 h1 6 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Fig. 5. Flux distribution analysis of P. pastoris metabolism under extreme DO conditions. A) Predicted metabolic flux distributions under extreme DO conditions (4% and 56% oxygen saturation) at μ ¼ 0.075 h1. Parsimonious FBA (Lewis et al., 2010) was employed to obtain a representative flux distribution under each condition. B) Feasible flux ranges for main central carbon reactions of P. pastoris under extreme oxygenation conditions. Flux ranges were calculated using loopless FVA (Saa and Nielsen, 2016a) as described in the Materials and Methods section. Abbreviations: HEX, hexokinase; PGI, glucose-6-phosphate isomerase; TPI, triose-phosphate isomerase; GAPD, glyceraldehyde-3-phosphate dehydrogenase; ENO, enolase; PYK, pyruvate kinase; G6PDH2, glucose-6-phosphate dehydrogenase; PGL, 6-phospho- gluconolactonase; GND, phosphogluconate dehydrogenase; RPI, ribose-5-phosphate isomerase; TKT1(2), transketolase 1(2); PYRt2m, pyruvate transport (mito- chondrial); PDH1a, pyruvate dehydrogenase; CSm, citrate synthase (mitochondrial); ACONTm, aconitate hydratase (mitochondrial); ICDHxm, isocitrate dehydrogenase NAD-dependent (mitochondrial); ICDHym, isocitrate dehydrogenase NADP-dependent (mitochondrial); AKGDbm, oxoglutarate dehydrogenase dihydrolipoamide S-succinyltransferase (mitochondrial); SUCOASm, oxoglutarate dehydrogenase dihydrolipoamide S-succinyltransferase (mitochondrial); MDHm, malate dehydrogenase (mitochondrial); PC, pyruvate carboxylase; PYRDC, pyruvate decarboxylase; PGMT, phosphoglucomutase. (Adelantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., optimal operational regions for achieving efficient growth and increased 2017). Likely causes of the latter considered stronger transcriptional productivity, respectively. While maximum qP is predicted to be at either induction of glycolysis, modification of membrane composition, and high or low DO and μ, optimal YS,X is strictly confined to the high DO and upregulation of UPR genes (Adelantado et al., 2017; Baumann et al., low μ region (Fig. 6). The former case is somewhat different from what 2010). The first cause is particularly supported by our data and flux has been reported (Hensing et al., 1995), specifically in the case of high μ simulations as the absolute glycolytic fluxes were higher under the L (> 0.2 h1) and DO (> 40%). A closer analysis of the statistical effects of condition, at a fixed μ of 0.075 h1 (Table 2 and Fig. 5A). In fact, μ and DO revealed, however, that μ only has a significant positive in- increased transcriptional levels of glycolytic genes under low DO levels fluence at very high DO conditions, and not vice versa (Fig. 2B). Meta- have been suggested as a plausible cause for higher protein production bolic studies are typically performed at low-medium oxygenation (< under the (glycolytic) GAP promoter (Baumann et al., 2010). Impor- 30%), which may explain why this synergistic effect has gone unnoticed tantly, our analysis revealed that the DO negative effect can be poten- for other proteins. Previous studies assessing protein secretion under the tially counteracted at high μ, yielding a more complex picture of the pGAP/TDH3 system reported an increase with μ at constant (low-- optimal protein production frontier (see below). medium) DO levels (Khasa et al., 2007; Maurer et al., 2006), which overall agrees with our results (Fig. 6). In the case of the DO effect, previous data showed that TDH3 expression – and hence pGAP-mediated 4.2. Optimal growth and thaumatin production lie in opposite operational protein expression – is upregulated at very low DO levels (i.e., hypoxic) regions (Baumann et al., 2010), which is consistent with the trend observed in our data. Importantly, this may help to explain the increase in qP Identification of the response surface of the macroscopic growth (i.e., (Table 2) even under higher energetic growth requirements (Fig. 4B). YS,X) and production (i.e., qP) parameters enables determination of the 7 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 glucose-limited conditions for various specific growth rates and constant DO level (Heyland et al., 2011), the TCA flux cycle is dramatically reduced under high DO concentrations (Fig. 5A). The relative ATP excess from the TCA under low DO can be then used to offset a seemingly higher ATP growth cost without great flux differences in other pathways like glycolysis and PPP (Fig. 5A). Altogether, our results suggest that the DO level is a critical operational parameter that has significant consequences on major central metabolic processes, as well as on the specific protein production performance of P. pastoris. 5. Conclusions Determination of the optimal culture conditions for microbial growth and production lies at the core of the optimization of any biotechno- logical process. Metabolic analysis of the impact of key culture parame- ters is challenging, as they typically require many experiments whose results are difficult to properly interpret. To tackle this limitation, we presented here a comprehensive framework for revealing high-order metabolic effects integrating rigorous design-of-experiments methodol- ogies with genome-scale metabolic modelling. Application of this Fig. 6. Phase plot of optimal P. pastoris growth and thaumatin production as a framework to the recombinant production of thaumatin in P. pastoris in function the specific growth rate and dissolved oxygen. The optimal conditions oxygen-sufficient, glucose-limited chemostats enabled mapping the for efficient growth (dashed line) and increased thaumatin production (contin- growth and production landscape of this yeast under different culture uous line) are shown in thick black lines, whereas suboptimal conditions are conditions. Furthermore, integration of the culture data with a high- represented with thin lines using the same symbols as before. As shown in the quality contextualized genome-scale metabolic model, revealed an un- figure, optimal growth lie opposite to optimal production conditions. suspected relationship between the DO level and the growth-associated ATP requirement, suggesting higher metabolic growth costs under rela- The latter highlights that recombinant protein production depends on a tively low oxygenation. Particularly for the case study presented here, complex interplay between the expression system, produced protein, and this information will be useful for designing more rational fermentation culture conditions (Looser et al., 2015). While the metabolic conse- strategies enabling uncoupling fast cellular growth from high protein quences of different DO levels have been reasonably well characterized in production in P. pastoris. this study, the underpinning mechanisms by which protein production are, likely, more complex and remain to be elucidated in future studies. Acknowledgements Finally, in the case of YS,X, our results indicated a positive correlation with DO that is consistent with the available data at constant μ (Ade- This project was funded by CORFO (Chile) grant No. 11CEII-9568 and lantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., 2017). MINECO/FEDER (Spain, EU) grant No. CTQ2016-74959-R. PT was sup- Overall, our model indicates that optimal growth and protein production ported by the Doctoral National Program (2014) grant No. 21140759 lie in opposite operational regions, highlighting the natural trade-off from CONICYT (Chile). PF and JA are supported by the Reference between these metabolic tasks. Network in Biotechnology (XRB) and grant No. 2014-SGR-452 from the Predictive GSMM reveals higher energetic growth requirements of Generalitat de Catalunya. Finally, we thank Alexandra Lobos and Vasni P. pastoris at lower DO levels. Zavaleta for technical assistance. The introduction of proper constraints to GSMMs based on experi- mental data and feasibility criteria is critical to achieve reliable flux Appendix A. Supplementary data predictions under the simulated conditions. Here, we have employed a rigorous step-by-step curation process where the predictive capabilities Supplementary data to this article can be found online at https://doi. of the GSMM were evaluated and improved (Fig. 3), which later enabled org/10.1016/j.mec.2019.e00103. us to evaluate the energetic impact of DO on P. pastoris growth. Adjust- ment of the GAME for different DO improved significantly the growth predictions of the GSMM (Fig. 4A), revealing an important energetic References effect of this variable on P. pastoris metabolism. Notably, the fitted GAME Adelantado, N., Tarazona, P., Grillitsch, K., García-Ortega, X., Monforte, S., Valero, F., values showed excellent agreement with available data (e.g., the pre- Feussner, I., Daum, G., Ferrer, P., 2017. The effect of hypoxia on the lipidome of dicted GAME under the M condition was 75 mmol ATP g1 DCW, similar to recombinant Pichia pastoris. Microb. Cell Factories 16, 86. the reported values of 70.5 and 72 mmol ATP g1 DCW under comparable DO Baumann, K., Carnicer, M., Dragosits, M., Graf, A.B., Stadlmann, J., Jouhten, P., Maaheimo, H., Gasser, B., Albiol, J., Mattanovich, D., Ferrer, P., 2010. A multi-level conditions (Caspeta et al., 2012; Tomas-Gamisans et al., 2016), which study of recombinant Pichia pastoris in different oxygen conditions. BMC Syst. Biol. prompted us to analyze the DO impact on the metabolism of P. pastoris. 4, 141. Our results showed a linear GAME decrease of 0.56 mmol ATP g1 DCW per Baumann, K., Maurer, M., Dragosits, M., Cos, O., Ferrer, P., Mattanovich, D., 2008. Hypoxic fed-batch cultivation of Pichia pastoris increases specific and volumetric each 1% increase in DO (Fig. 4B), which describes an important increase productivity of recombinant proteins. Biotechnol. Bioeng. 100, 177–183. in the growth efficiency at high DO levels (up to approx. 16% YS,X Buchetics, M., Dragosits, M., Maurer, M., Rebnegger, C., Porro, D., Sauer, M., Gasser, B., decrease between L and E, Table 2). Our results also indicate that the Mattanovich, D., 2011. Reverse engineering of protein secretion by uncoupling of cell extra metabolic capacity is not transferred to thaumatin production cycle phases from growth. Biotechnol. Bioeng. 108, 2403–2412. Calik, P., Ata, O., Gunes, H., Massahi, A., Boy, E., Keskin, A., Ozturk, S., Zerze, G.H., (Table 2) or associated to a substantial carbon redistribution (Fig. 5A), Ozdamar, T.H., 2015. Recombinant protein production in Pichia pastoris under which suggests DO plays a more complex regulatory role, likely related to glyceraldehyde-3-phosphate dehydrogenase promoter: from carbon source changes in cellular composition and/or transcriptional regulation (Bau- metabolism to bioreactor operation parameters. Biochem. Eng. J. 95, 20–36. Carnicer, M., Baumann, K., Toplitz, I., Sanchez-Ferrando, F., Mattanovich, D., Ferrer, P., mann et al., 2010; Gasser et al., 2007). Most notably, our results indicate Albiol, J., 2009. Macromolecular and elemental composition analysis and that, as opposed to previous studies under oxygen-sufficient, extracellular metabolite balances of Pichia pastoris growing at different oxygen levels. Microb. Cell Factories 8, 65. 8 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Caspeta, L., Shoaie, S., Agren, R., Nookaew, I., Nielsen, J., 2012. Genome-scale metabolic Love, K.R., Politano, T.J., Panagiotou, V., Jiang, B., Stadheim, T.A., Love, J.C., 2012. reconstructions of Pichia stipitis and Pichia pastoris and in silico evaluation of their Systematic single-cell analysis of Pichia pastoris reveals secretory capacity limits potentials. BMC Syst. Biol. 6, 24. productivity. PLoS One 7, e37915. Charoenrat, T., Ketudat-Cairns, M., Stendahl-Andersen, H., Jahic, M., Enfors, S.O., 2005. Maccani, A., Landes, N., Stadlmayr, G., Maresch, D., Leitner, C., Maurer, M., Gasser, B., Oxygen-limited fed-batch process: an alternative control for Pichia pastoris Ernst, W., Kunert, R., Mattanovich, D., 2014. Pichia pastoris secretes recombinant recombinant protein processes. Bioproc. Biosyst. Eng. 27, 399–406. proteins less efficiently than Chinese hamster ovary cells but allows higher space- Chelliah, V., Juty, N., Ajmera, I., Ali, R., Dumousseau, M., Glont, M., Hucka, M., time yields for less complex proteins. Biotechnol. J. 9, 526–537. Jalowicki, G., Keating, S., Knight-Schrijver, V., Lloret-Villas, A., Natarajan, K.N., Masuda, T., Ide, N., Ohta, K., Kitabatake, N., 2010. High-yield secretion of the Pettit, J.B., Rodriguez, N., Schubert, M., Wimalaratne, S.M., Zhao, Y.Y., recombinant sweet-tasting protein thaumatin I. Food Sci. Technol. Res. 16, 585–592. Hermjakob, H., Le Novere, N., Laibe, C., 2015. BioModels: ten-year anniversary. Masuda, T., Ohta, K., Ojiro, N., Murata, K., Mikami, B., Tani, F., Temussi, P.A., Nucleic Acids Res. 43, D542–D548. Kitabatake, N., 2016. A hypersweet protein: removal of the specific negative charge Ciplys, E., Zitkus, E., Gold, L.I., Daubriac, J., Pavlides, S.C., Hojrup, P., Houen, G., at Asp21 enhances thaumatin sweetness. Sci Rep-Uk 6, 20255. Wang, W.A., Michalak, M., Slibinskas, R., 2015. High-level secretion of native Mattanovich, D., Jungo, C., Wenger, J., Dabros, M., Maurer, M., 2014. Yeast suspension recombinant human calreticulin in yeast. Microb. Cell Factories 14, 165. culture. In: Meyer, H., Schmidhalter, D.R. (Eds.), Industrial Scale Suspension Culture Clasquin, M.F., Melamud, E., Singer, A., Gooding, J.R., Xu, X., Dong, A., Cui, H., of Living Cells. Campagna, S.R., Savchenko, A., Yakunin, A.F., Rabinowitz, J.D., Caudy, A.A., 2011. Maurer, M., Kühleitner, M., Gasser, B., Mattanovich, D., 2006. Versatile modeling and Riboneogenesis in yeast. Cell 145, 969–980. optimization of fed batch processes for the production of secreted heterologous Cos, O., Ramon, R., Montesinos, J.L., Valero, F., 2006. Operational strategies, monitoring proteins with Pichia pastoris. Microb. Cell Factories 5, 37-37. and control of heterologous protein production in the methylotrophic yeast Pichia Moralejo, F.J., Watson, A.J., Jeenes, D.J., Archer, D.B., Martin, J.F., 2001. A defined level pastoris under different promoters: a review. Microb. Cell Factories 5. of protein disulfide isomerase expression is required for optimal secretion of Delic, M., Gongrich, R., Mattanovich, D., Gasser, B., 2014. Engineering of protein folding thaumatin by Aspegillus awamori. Mol. Genet. Genom. 266, 246–253. and secretion-strategies to overcome bottlenecks for efficient production of Orth, J.D., Thiele, I., Palsson, B.O., 2010. What is flux balance analysis? Nat. Biotechnol. recombinant proteins. Antioxidants Redox Signal. 21, 414–437. 28, 245–248. Doehlert, D.H., 1970. Uniform shell designs. J. R Stat. Soc. Ser C (Appl Stat.) 19, 231–239. Pe~na, D.A., Gasser, B., Zanghellini, J., Steiger, M.G., Mattanovich, D., 2018. Metabolic Dragosits, M., Stadlmann, J., Albiol, J., Baumann, K., Maurer, M., Gasser, B., Sauer, M., engineering of Pichia pastoris. Metab. Eng. 50, 2–15. Altmann, F., Ferrer, P., Mattanovich, D., 2009. The effect of temperature on the Pereira, R., Nielsen, J., Rocha, I., 2016. Improving the flux distributions simulated with proteome of recombinant Pichia pastoris. J. Proteome Res. 8, 1380–1392. genome-scale metabolic models of Saccharomyces cerevisiae. Metab. Eng. Commu. 3, Dragosits, M., Stadlmann, J., Graf, A., Gasser, B., Maurer, M., Sauer, M., Kreil, D.P., 153–163. Altmann, F., Mattanovich, D., 2010. The response to unfolded protein is involved in Pirt, S.J., 1982. Maintenance energy: a general model for energy-limited and energy- osmotolerance of Pichia pastoris. BMC Genomics 11, 207. sufficient growth. Arch. Microbiol. 133, 300–302. Garcia-Ortega, X., Valero, F., Montesinos-Seguí, J.L., 2017. Physiological state as Price, N.D., Reed, J.L., Palsson, B.O., 2004. Genome-scale models of microbial cells: transferable operating criterion to improve recombinant protein production in Pichia evaluating the consequences of constraints. Nat. Rev. Microbiol. 2, 886–897. pastoris through oxygen limitation. J. Chem. Technol. Biotechnol. 92, 2573–2582. Prielhofer, R., Maurer, M., Klein, J., Wenger, J., Kiziak, C., Gasser, B., Mattanovich, D., Gasser, B., Mattanovich, D., 2018. A yeast for all seasons - is Pichia pastoris a suitable 2013. Induction without methanol: novel regulated promoters enable high-level chassis organism for future bioproduction? FEMS Microbiol. Lett. 365. expression in Pichia pastoris. Microb. Cell Factories 12, 5-5. Gasser, B., Maurer, M., Gach, J., Kunert, R., Mattanovich, D., 2006. Engineering of Pichia Puxbaum, V., Mattanovich, D., Gasser, B., 2015. Quo vadis? The challenges of pastoris for improved production of antibody fragments. Biotechnol. Bioeng. 94, recombinant protein folding and secretion in Pichia pastoris. Appl. Microbiol. 353–361. Biotechnol. 99, 2925–2938. Gasser, B., Maurer, M., Rautio, J., Sauer, M., Bhattacharyya, A., Saloheimo, M., Rebnegger, C., Vos, T., Graf, A.B., Valli, M., Pronk, J.T., Daran-Lapujade, P., Penttil€a, M., Mattanovich, D., 2007. Monitoring of transcriptional regulation in Mattanovich, D., 2016. Pichia pastoris exhibits high viability and a low maintenance Pichia pastoris under protein production conditions. BMC Genomics 8, 179-179. energy requirement at near-zero specific growth rates. Appl. Environ. Microbiol. 82, Gibson, D.G., Young, L., Chuang, R.-Y., Venter, J.C., Hutchison Iii, C.A., Smith, H.O., 4570. 2009. Enzymatic assembly of DNA molecules up to several hundred kilobases. Nat. Rußmayer, H., Buchetics, M., Gruber, C., Valli, M., Grillitsch, K., Modarres, G., Methods 6, 343. Guerrasio, R., Klavins, K., Neubauer, S., Drexler, H., Steiger, M., Troyer, C., Al Hanrahan, G., Lu, K., 2006. Application of factorial and response surface methodology in Chalabi, A., Krebiehl, G., Sonntag, D., Zellnig, G., Daum, G., Graf, A.B., Altmann, F., modern experimental design and optimization. Crit. Rev. Anal. Chem. 36, 141–151. Koellensperger, G., Hann, S., Sauer, M., Mattanovich, D., Gasser, B., 2015. Systems- Hasslacher, M., Schall, M., Hayn, M., Bona, R., Rumbold, K., Luckl, J., Griengl, H., level organization of yeast methylotrophic lifestyle. BMC Biol. 13, 80. Kohlwein, S.D., Schwab, H., 1997. High-level intracellular expression of Saa, P.A., Moenne, M.I., Perez-Correa, J.R., Agosin, E., 2012. Modeling oxygen dissolution hydroxynitrile lyase from the tropical rubber tree Hevea brasiliensis in microbial and biological uptake during pulse oxygen additions in oenological fermentations. hosts. Protein Expr. Purif. 11, 61–71. Bioproc. Biosyst. Eng. 35, 1167–1178. Hensing, M.C., Rouwenhorst, R.J., Heijnen, J.J., van Dijken, J.P., Pronk, J.T., 1995. Saa, P.A., Nielsen, L.K., 2016a. Fast-SNP: a fast matrix pre-processing algorithm for Physiological and technological aspects of large-scale heterologous-protein efficient loopless flux optimization of metabolic models. Bioinformatics 32, production with yeasts. Antonie Leeuwenhoek 67, 261–279. 3807–3814. Heyland, J., Fu, J., Blank, L.M., Schmid, A., 2010. Quantitative physiology of Pichia Saa, P.A., Nielsen, L.K., 2016b. ll-ACHRB: a scalable algorithm for sampling the feasible pastoris during glucose-limited high-cell density fed-batch cultivation for solution space of metabolic networks. Bioinformatics 32, 2330–2337. recombinant protein production. Biotechnol. Bioeng. 107, 357–368. Saitua, F., Torres, P., Perez-Correa, J.R., Agosin, E., 2017. Dynamic genome-scale Heyland, J., Fu, J., Blank, L.M., Schmid, A., 2011. Carbon metabolism limits recombinant metabolic modeling of the yeast Pichia pastoris. BMC Systems Biology 11, 27. protein production in Pichia pastoris. Biotechnol. Bioeng. 108, 1942–1953. Sohn, S.B., Graf, A.B., Kim, T.Y., Gasser, B., Maurer, M., Ferrer, P., Mattanovich, D., Hu, H., Gao, J., He, J., Yu, B., Zheng, P., Huang, Z., Mao, X., Yu, J., Han, G., Chen, D., Lee, S.Y., 2010. Genome-scale metabolic model of methylotrophic yeast Pichia 2013. Codon optimization significantly improves the expression level of a keratinase pastoris and its use for in silico analysis of heterologous protein production. gene in Pichia pastoris. PLoS One 8, e58393. Biotechnol. J. 5, 705–715. Hurvich, C.M., Tsai, C.L., 1989. Regression and time-series model selection in small Sola, A., Jouhten, P., Maaheimo, H., Sanchez-Ferrando, F., Szyperski, T., Ferrer, P., 2007. samples. Biometrika 76, 297–307. Metabolic flux profiling of Pichia pastoris grown on glycerol/methanol mixtures in Illingworth, C., Larson, G., Hellekant, G., 1989. Secretion of the sweet-tasting plant chemostat cultures at low and high dilution rates. Microbiology (Read., Engl.) 153, protein thaumatin by Streptomyces lividans. J. Ind. Microbiol. 4, 37–42. 281–290. Irani, Z.A., Kerkhoven, E.J., Shojaosadati, S.A., Nielsen, J., 2016. Genome-scale metabolic Sola, A., Maaheimo, H., Ylonen, K., Ferrer, P., Szyperski, T., 2004. Amino acid model of Pichia pastoris with native and humanized glycosylation of recombinant biosynthesis and metabolic flux profiling of Pichia pastoris. Eur. J. Biochem. 271, proteins. Biotechnol. Bioeng. 113, 961–969. 2462–2470. Jorda, J., Jouhten, P., C amara, E., Maaheimo, H., Albiol, J., Ferrer, P., 2012. Metabolic Theron, C.W., Berrios, J., Delvigne, F., Fickers, P., 2018. Integrating metabolic modeling flux profiling of recombinant protein secreting Pichia pastoris growing on glucose: and population heterogeneity analysis into optimizing recombinant protein methanol mixtures. Microb. Cell Factories 11, 57. production by Komagataella (Pichia) pastoris. Appl. Microbiol. Biotechnol. 102, Khasa, Y.P., Khushoo, A., Srivastava, L., Mukherjee, K.J., 2007. Kinetic studies of 63–80. constitutive human granulocyte-macrophage colony stimulating factor (hGM-CSF) Tomas-Gamisans, M., Ferrer, P., Albiol, J., 2016. Integration and validation of the expression in continuous culture of Pichia pastoris. Biotechnology Letters 29, genome-scale metabolic models of Pichia pastoris: a comprehensive update of protein 1903–1908. glycosylation pathways, lipid and energy metabolism. PLoS One 11, e0148031. Lewis, N.E., Hixson, K.K., Conrad, T.M., Lerman, J.A., Charusanti, P., Polpitiya, A.D., Tomas-Gamisans, M., Ferrer, P., Albiol, J., 2018. Fine-tuning the P. pastoris iMT1026 Adkins, J.N., Schramm, G., Purvine, S.O., Lopez-Ferrer, D., Weitz, K.K., Eils, R., genome-scale metabolic model for improved prediction of growth on methanol or K€onig, R., Smith, R.D., Palsson, B.Ø., 2010. Omic data from evolved E. coli are glycerol as sole carbon sources. Microb. Biotechno. 11, 224–237. consistent with computed optimal growth from genome-scale models. Mol. Syst. Biol. Uchiyama, K., Shioya, S., 1999. Modeling and optimization of alpha-amylase production 6, 390-390. in a recombinant yeast fed-batch culture taking account of the cell cycle population Liu, Z., Hou, J., Martinez, J.L., Petranovic, D., Nielsen, J., 2013. Correlation of cell growth distribution. J. Biotechnol. 71, 133–141. and heterologous protein production by Saccharomyces cerevisiae. Appl. Microbiol. Wang, Y., Liang, Z.H., Zhang, Y.S., Yao, S.Y., Xu, Y.G., Tang, Y.H., Zhu, S.Q., Cui, D.F., Biotechnol. 97, 8955–8962. Feng, Y.M., 2001. Human insulin from a precursor overexpressed in the Looser, V., Bruhlmann, B., Bumbak, F., Stenger, C., Costa, M., Camattari, A., Fotiadis, D., methylotrophic yeast Pichia pastoris and a simple procedure for purifying the Kovar, K., 2015. Cultivation strategies to enhance productivity of Pichia pastoris: a expression product. Biotechnol. Bioeng. 73, 74–79. review. Biotechnol. Adv. 33, 1177–1193. 9 P. Torres et al. Metabolic Engineering Communications 9 (2019) e00103 Xiang, L., Wang, Q., Zhou, Y., Yin, L., Zhang, G., Ma, Y., 2016. High-level expression of a Zhang, M., Yu, X.W., Xu, Y., Jouhten, P., Swapna, G.V.T., Glaser, R.W., Hunt, J.F., ZEN-detoxifying gene by codon optimization and biobrick in Pichia pastoris. Montelione, G.T., Maaheimo, H., Szyperski, T., 2017. 13 C metabolic flux profiling of Microbiol. Res. 193, 48–56. Pichia pastoris grown in aerobic batch cultures on glucose revealed high relative Ye, R., Huang, M., Lu, H., Qian, J., Lin, W., Chu, J., Zhuang, Y., Zhang, S., 2017. anabolic use of TCA cycle and limited incorporation of provided precursors of Comprehensive reconstruction and evaluation of Pichia pastoris genome-scale branched-chain amino acids. FEBS J. 284, 3100–3113. metabolic model that accounts for 1243 ORFs. Bioresourc. Bioprocess. 4, 22. 10