Contextualized genome-scale model unveils high-order metabolic effects of the specific growth rate and oxygenation level in recombinant Pichia pastoris.

Pichia pastoris is recognized as a biotechnological workhorse for recombinant protein expression. The metabolic performance of this microorganism depends on genetic makeup and culture conditions, amongst which the specific growth rate and oxygenation level are critical. Despite their importance, only their individual effects have been assessed so far, and thus their combined effects and metabolic consequences still remain to be elucidated. In this work, we present a comprehensive framework for revealing high-order (i.e., individual and combined) 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 relationship 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.

Introduction

The methylotrophic yeast Pichia pastoris (syn. Komagataella phaffii) is a popular cell factory for heterologous protein production due to its convenient physiology and the availability of genetic tools for its manipulation (Gasser and Mattanovich, 2018; Peña et al., 2018; Theron et al., 2018). Amongst its features, P. pastoris displays a Crabtree-negative growth phenotype, leading to reduced synthesis of by-products under aerobic conditions (Baumann et al., 2010; Calik et al., 2015). In addition, the availability of strong promoters for heterologous expression, added to its natural capacity to perform post-translational modifications, has made this yeast the platform of choice for producing a range of recombinant proteins up to grams per liter titers (Ciplys et al., 2015; Hasslacher et al., 1997; Heyland et al., 2010; Wang et al., 2001). Despite these achievements, efficient heterologous protein production in P. pastoris remains challenging, as poorly tuned protein overexpression can affect relevant cellular processes, such as protein folding and secretion (Delic et al., 2014; Gasser et al., 2007; Love et al., 2012). Moreover, codon usage level (Hu et al., 2013; Xiang et al., 2016), promoter selection (Prielhofer et al., 2013), as well as culture medium composition (Heyland et al., 2011) and operational conditions (Cos et al., 2006; Maurer et al., 2006) may also play major roles on process performance. In particular, the operational conditions have gained increasing attention as they are known to introduce substantial variability in the process, significantly affecting the recombinant protein secretion (Looser et al., 2015).

High recombinant protein expression in P. pastoris relies on the use of strong promoters, like pAOX1 (promoter from alcohol oxidase I encoding gene) and pGAP (promoter from glyceraldehyde-3-phosphate dehydrogenase encoding gene). While pAOX1 offers strong inducible expression with methanol – thereby enabling uncoupling fast growth from production –, pGAP provides comparable constitutive expression (Peña et al., 2018). P. pastoris cultures incur in high oxygen consumption and heat production during methanol oxidation, and hence, its use poses major challenges for large-scale protein production (Mattanovich et al., 2014). Once a suitable expression system has been chosen, the next step is to optimize culture conditions to achieve the target productivity. Factors such as temperature, pH, osmolality, specific growth rate (μ) and dissolved oxygen (DO) are critical for the efficient operation of the culture, and their influence on protein production and culture performance has been individually assessed (Baumann et al., 2008; Charoenrat et al., 2005; Dragosits et al., 2009, 2010; Garcia-Ortega et al., 2017; Heyland et al., 2010; Maurer et al., 2006). Although several studies have reviewed the relationships between protein production and growth (refer to Looser et al. (2015) for a comprehensive review), and how DO impacts the yeast’s physiology (Adelantado et al., 2017; Baumann et al., 2010; Garcia-Ortega et al., 2017), current studies fail to evaluate both the individual and combined (high-order) effects of these operational parameters on the metabolic performance of P. pastoris. Notably, such evaluation is not trivial as both, an appropriate experimental design and modelling framework, are required to appropriately interpreting the results.

In this study, we present an integrated framework for mapping, at the genome–scale, the individual and combined metabolic effects of key operational parameters – here μ and DO – during recombinant protein production by P. pastoris, under glucose-limited conditions in continuous cultures. As a case study, we analyzed the metabolic behavior of a recombinant strain producing the sweet-tasting, low-calorie protein thaumatin. This protein has 207 amino acid residues and 8 disulfide bonds (Illingworth et al., 1989), which are critical for its sweet taste (Masuda et al., 2016) and are considered the main reason behind the low titers achieved so far (Moralejo et al., 2001) (~ 100 mg L−1 in high-density cell cultures (Masuda et al., 2010)). Folding of recombinant proteins with many disulfide bounds is both difficult and costly, as it requires a high supply of NAD(P)H cofactors which can affect redox homeostasis and lead to negative physiological responses like the Unfolded Protein Response (UPR) and Endoplasmic-Reticulum-Associated Degradation (ERAD) (Gasser et al., 2007; Puxbaum et al., 2015). Thus, understanding the effects of μ and DO that have a major metabolic impact is critical for optimizing heterologous protein production in P. pastoris. To capture the individual and combined effects of μ and DO on the metabolic performance of this strain, we implemented a fractional experimental design to determine the statistically significant effects from a reduced set of experiments, and later integrated these with a tailored Genome-Scale Metabolic Model (GSMM) to interpret and further analyze the metabolic and energetic consequences. GSMMs describe the set of biochemical reactions that represent the entire metabolic potential of an organism, and as such, are excellent tools to understand the physiological state of the cell under different environmental conditions (Price et al., 2004). Our results revealed opposite operational regions for efficient growth and protein production, as well as a tighter energetic constraint on P. pastoris growth under glucose-limited, low DO conditions.

Materials and Methods

Plasmid construction and strain transformation

The thaumatin gene – including its natural pre-region secretion signal – was synthesized by Genscript (Piscataway, NJ, USA) and was codon-optimized for expression in P. pastoris. The commercial recombinant vectors pGAPZB (Invitrogen, Carlsbad, CA, USA) were employed for the construction of the transformation plasmids using Gibson assembly (Gibson et al., 2009), and later used to integrate the thaumatin gene into the yeast’s genome under the constitutive GAP promoter. The sequences of the employed primers are shown in Table 1. Amplification of DNA fragments for Gibson assembly was carried out by PCR using Phusion High-Fidelity DNA Polymerase (ThermoFisher, Waltham, MA, USA). All PCR products were treated with DpnI restriction enzyme to remove original vector residues and later were purified by gel extraction, using the Qiaquick Gel Extraction kit (Qiagen, Hilden, Germany). Finally, chemically competent E. coli TOP 10 cells were transformed using the sought construct. These cells were grown at 37 °C in low salt-LB medium, containing 25 μg mL−1 zeocin for selection of clones transformed with pGAPZB-TAU vector.

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

P. pastoris wild-type strain GS115 (Invitrogen, Carlsbad, CA, USA) was used as a host strain throughout this study, which was transformed using an in-house-built vector to revert its histidine auxotrophy (refer to Supplementary Text S1). AvrII was employed to linearize the transformation vector, which was introduced by electroporation into the competent cells, as described by Gasser et al. (2006). Both plasmids and transformations were verified by DNA sequencing (Macrogen Inc., Seoul, Korea).

Cell cultivation

Continuous cultures were started from pre-inocula grown overnight at 30 °C and 150 rpm in 200-mL shake flasks, containing YPG medium with 100 μg mL−1 zeocin. Prior to the inoculation of the bioreactors, each inoculum was centrifuged at 5000 rpm for 5 min and resuspended in fresh culture medium without trace elements. Chemostat cultures were performed in 2-L benchtop Biostat B bioreactors (Sartorius AG, Göttingen, Germany) using 1 L of working volume.

P. pastoris was grown in glucose-limited chemostats under the culture conditions described by Adelantado et al. (2017) with only minor differences (see Supplementary Text S2). For DO control, a mixture of air, pure nitrogen and pure oxygen was used. Each gas was independently controlled using three mass flow meters (Bronkhorst High-Tech, Ruurlo, Netherlands) until the desired DO levels were reached (see below). The pH, agitation, and temperature were set respectively to 5.0, 700 rpm, and 25 °C in all conditions. Samples were taken after four and five residence times, when the culture reached steady state.

Experimental design

Experimental design is a statistical tool that is used to systematically examine the effects of two or more factors that shape a determined response function. Among the different possible designs, Doehlert design displays higher efficiency and flexibility, as it allows fitting complex response models with less data than complete designs, and without requiring extreme points difficult to acquire (Doehlert, 1970; Hanrahan and Lu, 2006).

In this study, the Doehlert matrix (Doehlert, 1970) was constructed to evaluate the individual and combined effects of μ and DO on the macroscopic growth and production parameters of the thaumatin-producing P. pastoris strain. The experimental domain was defined over 0.05–0.15 h−1 for μ with five dilution rate levels; 0.05, 0.075, 0.1, 0.125 and 0.15 h−1 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 M−1 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 polynomial (Eq. (1)) was fitted to determine the response of the dependent variable as a function of the different factors, i.e. μ and DO.

Fig. 1

Overview of the experimental design and GSMM contextualization workflow employed in this study. A) Doehlert fractional experimental design for evaluating the metabolic impact of the specific growth rate (μ) and dissolved oxygen (DO) level on the metabolic performance of P. pastoris. This design considers 1 central point and 6 extreme points over 3 DO levels (4%, 30% and 56% oxygen saturation), and 5 μ levels (0.05, 0.075, 0.1, 0.125 and 0.15 h−1). B) Manual GSMM contextualization workflow for flux simulations. The latest P. pastoris GSMM (iMT1026 v3.0) was incrementally constrained for accurately describing metabolic fluxes under glucose-limited conditions using reported experimental data and flux feasibility criteria.

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 describes an interaction term (combined effect). Statistical calculations and significant tests (ANOVA) were performed using STATGRAPHICS Plus 5.1 (Statpoint Technologies, Inc., VA, USA).

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 Tomàs-Gamisans et al. (2016). Extracellular thaumatin was quantified in triplicate using Thaumatin ELISA Kit (CellTrend, Luckenwalde, Germany), according to manufacturer instructions. Finally, Dry Cell Weight (DCW) was quantified in triplicate, following the method of Jordà et al. (2012).

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; Tomàs-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. (2016).

In order to achieve accurate and meaningful predictions underpinned by realistic flux distributions, we set up a sequential refinement workflow starting from the latest and most curated P. pastoris GSMM – the iMT1026 v3.0 (Tomàs-Gamisans et al., 2018) –, where intracellular reactions were incrementally constrained (both activities and directionalities), based on reported fluxomics, proteomics and transcriptomics data under appropriate experimental conditions (Clasquin et al., 2011; Gasser et al., 2007; Rußmayer et al., 2015; Zhang et al., 2017), and the imposition of the loopless flux condition (elimination of thermodynamically infeasible closed cycles) (Saa and Nielsen, 2016a, 2016b) (Fig. 1B). In this way, the refined metabolic model generated realistic intracellular fluxes, i.e. consistent with reported metabolic pathway utilization behaviors and flux patterns and, most importantly, accurately predicted observed specific growth and exchange rates. Additional details about the contextualization and refinement workflow can be found in Supplementary Text S3 and Tables S1–S3. The contextualized constrained metabolic model is available in Supplementary File S1 and has been deposited in the BioModels database (Chelliah et al., 2015) under MODEL1910040001 identifier. A MATLAB script for generating the contextualized models is also available in Supplementary File S1.

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.

So far, the energetic impact of the oxygenation on P. pastoris metabolism 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 h−1 (Adelantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., 2017). In contrast, the experimental design proposed here enables examination 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 latter, the qATP of the refined GSMM was maximized in all the experimental conditions using the loopless FBA (ll-FBA) algorithm (Saa and Nielsen, 2016a) – i.e., a flux estimator that avoids thermodynamically infeasible loops. In all the simulations, the macromolecular biomass composition employed was based on available data from glucose-limited cultures (Carnicer et al., 2009), which is available in iMT1026 v3.0. For the qATP estimation, the mATP reaction was chosen as the objective, the GAME was set to zero, and the growth and exchange rates were fixed to the observed values. Once qATP was computed, the GAME was fitted to the different DO conditions using Eq. (2) after setting NGAME to the most recent measured value (Rebnegger et al., 2016). Another tested model considered fitting NGAME for different DO conditions and setting GAME to the GSMM value of 72 mmol ATP gDCW−1 (Tomàs-Gamisans 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; another, where only YX,ATP was fitted to each DO condition; and a final model, where only NGAME was fitted to each DO condition. Finally, the Akaike’s Information Criterion for low number of data points (AICc) was applied to compare the goodness-of-fit of each model (Hurvich and Tsai, 1989). The model with the smallest AICc was considered the best choice.

Results

Effects of the specific growth rate and oxygenation level on the metabolic performance of P. pastoris under glucose-limited conditions

Growth and metabolic production parameters of the recombinant P. pastoris strain were measured in seven chemostat cultures (conditions), under different μ and oxygenation conditions. In each case, biomass and thaumatin yields (YS,X and YS,P, respectively), specific thaumatin and carbon dioxide production rates (qP and qCO2, respectively), and specific glucose and oxygen consumption rates (qS and qO2, respectively) were estimated (Table 2). In all conditions, ethanol and arabitol were not detected, pointing to a purely respiratory metabolism, which was later confirmed with measured RQ values (Supplementary Table S5).

Table 2

Metabolic performance parameters of recombinant P. pastoris in glucose-limited chemostats under different μ and DO conditions.a.

Conditionb	μ (h−1)	qS (mmol gDCW−1 h−1)	qCO2 (mmol gDCW−1 h−1)	qO2 (mmol gDCW−1 h−1)	qP (mg gDCW−1 h−1)	YS,X (gDCW g−1)	YS,P (mg g−1)
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

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

The experimental conditions are abbreviated L for low oxygenation (4% dissolved oxygen saturation), M for medium oxygenation (30% dissolved oxygen saturation), and E for excess oxygenation (56% dissolved oxygen saturation).

Oxygenation exerted the greatest impact on both, YS,X and qP, whereas μ had the main effect on YS,P. In the case of YS,X, there were notable differences between different oxygenation levels at a given μ (e.g., L and E conditions showed a 14% relative difference for μ = 0.075 h−1), but not vice versa (< 4.7% relative difference for a fixed oxygenation level and different μ). Overall, YS,X varied from 0.507 to 0.589 gDCW g−1, attaining the lowest yield at μ = 0.075 h−1 in the L condition. Statistical analysis showed that DO had the main positive and significant effect (p-value<0.001, Supplementary Table S6), although the interactive and quadratic terms displayed also significant (both negative) influences (p-value<0.01 and p-value<0.05, respectively; refer to Supplementary Table S6 for details). The YS,X regression explained a substantial part of the variability in the data (R2 = 0.948) and, more importantly, it suggested high oxygenation and low μ as the most convenient conditions for efficient growth (Fig. 2A).

Fig. 2

Response surfaces of P. pastoris macroscopic culture parameters under different DO and μ conditions. The fitted response surfaces were for the (encoded) biomass yield (A), specific thaumatin production rate (B) and thaumatin yield (C) are, respectively: (R2 = 0.948), (R2 = 0.68), and (R2 = 0.959), where and . Further details about the regression and statistical effects can be found in Supplementary Table S6.

In the case of the thaumatin production parameters qP and YS,P, they showed different behaviors. While YS,P was almost unaffected by DO and displayed a negative correlation with μ (e.g., varying from 0.129 to 0.399 mg thaumatin g−1 in the high and low μ conditions, respectively), qP was largely influenced by DO (negatively); and to a lesser extent by μ (positively) (Table 2 and Fig. 2B). In this case, the regression analysis explained a lower fraction of the variability of the data (R2 = 0.68) and indicated a significant negative effect of DO (p-value<0.05; Supplementary Table S6) with a positive interaction term (p-value<0.05, Supplementary Table S6). Here, the magnitudes of both, the individual and combined effects, were quite similar (see Supplementary Table S6). In the case of the YS,P, significant individual (negative) and quadratic (positive) effects of μ were suggested by the regression analysis (in both cases p-value<0.001, and R2 = 0.959, see Supplementary Table S6). Altogether, regions of both, low μ and oxygenation, favor thaumatin production and yield (Fig. 2B–C).

GSMM contextualization and refinement

The prediction power of iMT1026 v3.0 was improved through a sequential refinement workflow. In a first step, the reaction directionalities and activities were modified based on literature data. At this point, it is important to note that iMT1026 v3.0 supports growth on glycerol and methanol through the introduction of accessory pathways in different compartments (e.g., peroxisomal pentose phosphate pathway, PPP) (Tomàs-Gamisans et al., 2018). The activity of such pathways is, however, absent under glucose-limited conditions, as a previous report suggests (Rußmayer et al., 2015). Fail to block these pathways led to the simulation of unrealistically high glycolytic fluxes (Supplementary Fig. S1). Also, the cytosolic NADP-dependent isocitrate dehydrogenase was not blocked, as previously suggested (Pereira et al., 2016; Saitua et al., 2017), as there is no reported evidence for the lack of activity in P. pastoris. Another modification to the model considered fixing the sedoheptulose 1,7-bisphosphate D-glyceraldehyde-3-phosphate-lyase reaction in the direction of riboneogenesis (Clasquin et al., 2011), which improved considerably the flux predictions in the non-oxidative branch of the PPP, as compared to above 13C labeling flux data. Cytosolic acetyl-CoA synthesis was forced to follow the canonical pathway through pyruvate by blocking the methylmalonate-semialdehyde dehydrogenase. The directionality of the mitochondrial malate dehydrogenase was fixed in the direction of malate oxidation, consistent with an inactive glyoxylate cycle under various carbon sources (Sola et al., 2004, 2007) and more recent 13C flux data in aerobic, glucose-limited cultures (Sola et al., 2004; Zhang et al., 2017). Finally, the activity of the malic enzyme was blocked as it carries a neglectable flux under the hypoxic and normoxic conditions in glucose-limited cell cultures (Baumann et al., 2010; Zhang et al., 2017). Altogether, correction of these reactions substantially improved mitochondrial flux predictions, as compared to available 13C flux measurements (see below).

Integration of transcriptomic and proteomic data under appropriate experimental conditions (Clasquin et al., 2011; Gasser et al., 2007; 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, subsequent thermodynamic feasibility analysis, based on imposition of the loopless flux condition, yielded irreversible 215 reactions previously deemed reversible (Fig. 3A). Altogether, these measures disabled infeasible cycles responsible for the flux prediction inaccuracies. Next, adjustment of the NGAME to 0.55 mmol ATP gDCW−1 h−1 (assuming a yield of 32 mol of ATP per mol of glucose (Tomàs-Gamisans et al., 2016), and using the recent experimental maintenance value of 3.1 mg glucose gDCW−1 h−1 measured in glucose-limited retentostats (Rebnegger et al., 2016), increased the prediction fidelity of iMT1026 v3.0 by ca. 86% for the observed μ – from approx. 7.0% average relative error to 4.2% (Fig. 3B). Moreover, the agreement between predicted and measured intracellular flux data improved considerably – explained flux variability increased from 78.3% to 97.2% in the cytosol (Fig. 3C) and from 7.6% to 86.7% in the mitochondria (Fig. 3D). We note that iMT1026 v3.0 uses a biomass macromolecular composition with ashes, although correction by the estimated ashes content yielded almost identical results (Supplementary Fig. S2). Finally, the flux ranges of the contextualized model were slightly 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.

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.

GAME evaluation under oxygen-sufficient, glucose-limited conditions

The contextualized iMT1026 v3.0 was employed to evaluate the energetic 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 Estimation of energetic requirements for growth for details). Our results showed that the model with a variable GAME as a function of the DO level explains 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 gDCW−1 for the L condition, and descending to 75 and 66 mmol ATP gDCW−1 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.

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 growth rate predictions of the model with a fixed GAME for all DO conditions (slope = 0.752, R2 = 0.871), the model with a variable GAME for each DO condition (slope = 0.844, R2 = 0.949), and the model with a variable NGAME for 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 stronger statistical support as shown by the lower AICc value (-158.9 versus -146.8 and -154.7). B) Predicted GAME values for different DO conditions (95 mmol ATP gDCW−1 for the L condition, 75 mmol ATP gDCW−1 for the M condition, and 66 mmol ATP gDCW−1 for the E conditions).

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 conditions were simulated at a fixed μ level (μ = 0.075 h−1). 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 biosynthetic 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 increased glucose uptake rate in the L condition, which provides more flexibility for carbon allocation (Fig. 5B).

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 h−1. 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-phosphogluconolactonase; GND, phosphogluconate dehydrogenase; RPI, ribose-5-phosphate isomerase; TKT1(2), transketolase 1(2); PYRt2m, pyruvate transport (mitochondrial); 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.

Discussion

Oxygenation exerts an individual negative effect and a synergistic positive effect with the specific growth rate on thaumatin production

The specific thaumatin production rate was influenced negatively by the DO level, and positively by μ through the interaction with DO (Fig. 2). The latter effect is consistent with previous reports of a positive correlation between μ and protein secretion rates at fixed medium DO level, although the estimated effect was not so strong as previously reported (Buchetics et al., 2011; Maurer et al., 2006). Recent studies have showed that specific protein production in P. pastoris can also depend on the nature of the produced protein, expression system and oxidative conditions (Gasser et al., 2006; Liu et al., 2013; Maccani et al., 2014; Uchiyama and Shioya, 1999). Indeed, despite having constitutive protein expression, non-growth-associated protein production patterns may arise due to global stress responses, high protein turnover and limitations in endoplasmic reticulum processing, among others (Liu et al., 2013). The latter points to protein expression, processing and/or secretion as likely – and increasingly more recognized – bottlenecks (Looser et al., 2015).

The calculated individual DO effect on protein production agrees well with previous studies using the pGAP-based expression system, reporting higher protein secretion under hypoxic conditions at μ = 0.1 h−1 (Adelantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., 2017). Likely causes of the latter considered stronger transcriptional induction of glycolysis, modification of membrane composition, and upregulation of UPR genes (Adelantado et al., 2017; Baumann et al., 2010). The first cause is particularly supported by our data and flux simulations as the absolute glycolytic fluxes were higher under the L condition, at a fixed μ of 0.075 h−1 (Table 2 and Fig. 5A). In fact, increased transcriptional levels of glycolytic genes under low DO levels have been suggested as a plausible cause for higher protein production under the (glycolytic) GAP promoter (Baumann et al., 2010). Importantly, our analysis revealed that the DO negative effect can be potentially counteracted at high μ, yielding a more complex picture of the optimal protein production frontier (see below).

Optimal growth and thaumatin production lie in opposite operational regions

Identification of the response surface of the macroscopic growth (i.e., YS,X) and production (i.e., qP) parameters enables determination of the optimal operational regions for achieving efficient growth and increased productivity, respectively. While maximum qP is predicted to be at either high or low DO and μ, optimal YS,X is strictly confined to the high DO and low μ region (Fig. 6). The former case is somewhat different from what has been reported (Hensing et al., 1995), specifically in the case of high μ (> 0.2 h−1) and DO (> 40%). A closer analysis of the statistical effects of μ and DO revealed, however, that μ only has a significant positive influence at very high DO conditions, and not vice versa (Fig. 2B). Metabolic studies are typically performed at low-medium oxygenation (< 30%), which may explain why this synergistic effect has gone unnoticed for other proteins. Previous studies assessing protein secretion under the pGAP/TDH3 system reported an increase with μ at constant (low-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 protein expression – is upregulated at very low DO levels (i.e., hypoxic) (Baumann et al., 2010), which is consistent with the trend observed in our data. Importantly, this may help to explain the increase in qP (Table 2) even under higher energetic growth requirements (Fig. 4B). The latter highlights that recombinant protein production depends on a complex interplay between the expression system, produced protein, and culture conditions (Looser et al., 2015). While the metabolic consequences of different DO levels have been reasonably well characterized in this study, the underpinning mechanisms by which protein production are, likely, more complex and remain to be elucidated in future studies.

Fig. 6

Phase plot of optimal P. pastoris growth and thaumatin production as a function the specific growth rate and dissolved oxygen. The optimal conditions for efficient growth (dashed line) and increased thaumatin production (continuous line) are shown in thick black lines, whereas suboptimal conditions are represented with thin lines using the same symbols as before. As shown in the figure, optimal growth lie opposite to optimal production conditions.

Finally, in the case of YS,X, our results indicated a positive correlation with DO that is consistent with the available data at constant μ (Adelantado et al., 2017; Baumann et al., 2008; Garcia-Ortega et al., 2017). Overall, our model indicates that optimal growth and protein production lie in opposite operational regions, highlighting the natural trade-off between these metabolic tasks.

Predictive GSMM reveals higher energetic growth requirements of P. pastoris at lower DO levels.

The introduction of proper constraints to GSMMs based on experimental data and feasibility criteria is critical to achieve reliable flux predictions under the simulated conditions. Here, we have employed a rigorous step-by-step curation process where the predictive capabilities of the GSMM were evaluated and improved (Fig. 3), which later enabled us to evaluate the energetic impact of DO on P. pastoris growth. Adjustment of the GAME for different DO improved significantly the growth predictions of the GSMM (Fig. 4A), revealing an important energetic effect of this variable on P. pastoris metabolism. Notably, the fitted GAME values showed excellent agreement with available data (e.g., the predicted GAME under the M condition was 75 mmol ATP gDCW−1, similar to the reported values of 70.5 and 72 mmol ATP gDCW−1 under comparable DO conditions (Caspeta et al., 2012; Tomàs-Gamisans et al., 2016), which prompted us to analyze the DO impact on the metabolism of P. pastoris. Our results showed a linear GAME decrease of 0.56 mmol ATP gDCW−1 per each 1% increase in DO (Fig. 4B), which describes an important increase in the growth efficiency at high DO levels (up to approx. 16% YS,X decrease between L and E, Table 2). Our results also indicate that the extra metabolic capacity is not transferred to thaumatin production (Table 2) or associated to a substantial carbon redistribution (Fig. 5A), which suggests DO plays a more complex regulatory role, likely related to changes in cellular composition and/or transcriptional regulation (Baumann et al., 2010; Gasser et al., 2007). Most notably, our results indicate that, as opposed to previous studies under oxygen-sufficient, 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.

Conclusions

Determination of the optimal culture conditions for microbial growth and production lies at the core of the optimization of any biotechnological process. Metabolic analysis of the impact of key culture parameters 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 methodologies with genome-scale metabolic modelling. Application of this framework to the recombinant production of thaumatin in P. pastoris in oxygen-sufficient, glucose-limited chemostats enabled mapping the growth and production landscape of this yeast under different culture conditions. Furthermore, integration of the culture data with a high-quality contextualized genome-scale metabolic model, revealed an unsuspected relationship between the DO level and the growth-associated ATP requirement, suggesting higher metabolic growth costs under relatively low oxygenation. Particularly for the case study presented here, this information will be useful for designing more rational fermentation strategies enabling uncoupling fast cellular growth from high protein production in P. pastoris.