(C) PLOS One
This story was originally published by PLOS One and is unaltered.
. . . . . . . . . .



Genome-scale metabolic modeling reveals metabolic trade-offs associated with lipid production in Rhodotorula toruloides [1]

['Alīna Reķēna', 'Department Of Chemistry', 'Biotechnology', 'Tallinn University Of Technology', 'Tallinn', 'Marina J. Pinheiro', 'Department Of Materials', 'Bioprocess Engineering', 'University Of Campinas', 'Campinas']

Date: 2023-05

Rhodotorula toruloides is a non-conventional, oleaginous yeast able to naturally accumulate high amounts of microbial lipids. Constraint-based modeling of R. toruloides has been mainly focused on the comparison of experimentally measured and model predicted growth rates, while the intracellular flux patterns have been analyzed on a rather general level. Hence, the intrinsic metabolic properties of R. toruloides that make lipid synthesis possible are not thoroughly understood. At the same time, the lack of diverse physiological data sets has often been the bottleneck to predict accurate fluxes. In this study, we collected detailed physiology data sets of R. toruloides while growing on glucose, xylose, and acetate as the sole carbon source in chemically defined medium. Regardless of the carbon source, the growth was divided into two phases from which proteomic and lipidomic data were collected. Complemental physiological parameters were collected in these two phases and altogether implemented into metabolic models. Simulated intracellular flux patterns demonstrated the role of phosphoketolase in the generation of acetyl-CoA, one of the main precursors during lipid biosynthesis, while the role of ATP citrate lyase was not confirmed. Metabolic modeling on xylose as a carbon substrate was greatly improved by the detection of chirality of D-arabinitol, which together with D-ribulose were involved in an alternative xylose assimilation pathway. Further, flux patterns pointed to metabolic trade-offs associated with NADPH allocation between nitrogen assimilation and lipid biosynthetic pathways, which was linked to large-scale differences in protein and lipid content. This work includes the first extensive multi-condition analysis of R. toruloides using enzyme-constrained models and quantitative proteomics. Further, more precise k cat values should extend the application of the newly developed enzyme-constrained models that are publicly available for future studies.

Transition towards a biobased, circular economy to reduce the industrial dependence on fossil-based resources requires new technologies. One of the options is to convert available biomass feedstocks into valuable chemicals using microbes as biocatalysts. Rhodotorula toruloides is a nonpathogenic, nonconventional yeast that has recently emerged as one of the most promising yeasts for sustainable production of chemicals and fuels due to its natural ability to synthesize large amounts of lipids. However, its unique metabolic properties are not yet fully understood. We have computationally predicted metabolic fluxes in R. toruloides while growing in economically viable growth conditions inducing lipid accumulation and analyzed them together with absolute proteome quantification. Our holistic approach has highlighted metabolic pathways important for lipid biosynthesis and revealed metabolic trade-offs associated with NADPH allocation during lipogenesis. In addition, our work highlighted the necessity for accurate computational approaches in characterizing enzymatic kinetic properties that would improve the metabolic studies of R. toruloides.

In the present study, we created condition-specific enzyme-constrained genome-scale metabolic models of R. toruloides, ecRhtoGEMs, and used them to predict intracellular fluxes. Flux bounds to constrain the model were obtained from bioreactor (1 L) experiments with yeast cultivation in chemically defined medium, with three carbon sources studied individually—glucose, xylose and acetate. These very detailed physiological data sets enabled us to precisely characterize metabolism at exponential growth and lipid accumulation phase. In all conditions, we performed mass spectrometry (MS) based absolute proteome quantification. Also, biomass macromolecular composition in regard to lipids and proteins was determined, including lipid profiling by gas chromatography (GC) analysis. Using this data, we generated 6 different versions of the R. toruloides model with enzyme constraints and biomass composition specificity, where we were able to demonstrate trade-offs in NADPH requirements for the cells growing exponentially versus in nitrogen limitation. To our knowledge, this is the first proteomics analysis of acetate-grown R. toruloides cells and the first detailed GEM analysis combined with proteome analysis of acetate and glucose conditions for this strain.

A better understanding of how different metabolic pathways contribute to lipid accumulation under different substrates would help to design better metabolic engineering strategies. GEMs can be a powerful and helpful tool in metabolic studies, if their predictive power is good. Enzyme-constrained GEMs integrate additional constraints on enzyme capacity and their total abundances (as thoroughly reviewed by Chen and Nielsen 2021 [ 16 ]). Phenomenological constraint is imposed on metabolic flux (v; mmol/gDCW/h), formulated as enzyme kinetics ( Eq 1 ) (1) where E is protein abundance (mmol/gDCW) and k cat is the enzyme’s turnover number (1/s), provided with an upper limit on individual or total protein abundances The integration of enzymatic constraints in S. cerevisiae significantly improved phenotype prediction [ 17 ]. The strength of proteome constraints has also been demonstrated by predicting overflow metabolism in E. coli [ 18 ] and metabolic shift in arginine catabolism in L. lactis [ 19 ]. A similar coarse-grained approach that allowed the prediction of maximal growth without constraining the model with any exchange fluxes in S. cerevisiae was demonstrated by applying a global thermodynamics constraint [ 20 ].

A holistic view on metabolism can be provided by genome-scale metabolic models (GEMs). GEMs are metabolic networks reconstructed from a genome sequence of a specific organism. They contain all known biochemical reactions of the organism. GEMs allow the calculation of metabolic fluxes that represent activity of different metabolic pathways under specified conditions, e.g., an uptake of a particular carbon source. GEMs of R. toruloides were built based on the genome sequences of strains NP11 [ 15 ] and IFO 0880 [ 11 ]. Flux balance analysis predicted that up to 87% of NADPH was regenerated from xylose through the oxidative part of pentose phosphate pathway (oxPPP) [ 1 ]. Phosphoketolase was suggested as the main supplier of acetyl-CoA during lipogenesis in xylose-grown cells [ 1 ]. On the other hand, TCA cycle related enzymes were suggested for NADPH production on acetate-grown cells [ 7 ], demonstrating that metabolic operations can vary significantly with the carbon source uptake. Models have also been used to study metabolism during cell growth on glucose [ 11 ] and glycerol [ 7 ].

It has been reported that on xylose R. toruloides is growing 2 to 3 times slower compared to glucose [ 8 ], but the underlying mechanisms are yet to be identified. In our previous proteomics study [ 1 ], we discovered from proteomics quantification that xylulokinase, encoded in the genome as the third step in the currently known xylose assimilation pathway, is not present in the proteomic data set, suggesting potential limitation in xylose metabolism. Later, a similar finding was reported by Jagtap et al. 2021 [ 3 ] and Kim et al. 2021 [ 9 ] using a different R. toruloides strain, IFO 0880. The latter proposed an alternative xylose assimilation pathway for this species.

Metabolic pathways producing intracellular metabolite acetyl-CoA and a cofactor NADPH in R. toruloides have been the main focus of metabolic studies due to their central role in lipid biosynthesis. Fatty acids, which mainly accumulate in the form of triacylglycerols (TAGs), are produced via the sequence of four enzymatic reactions that require 1 ATP and 2 NADPH molecules per 1 acetyl-CoA added to the fatty acid chain [ 13 ]. To study lipid metabolism in R. toruloides, previous studies have taken the systems biology approach. Genome sequencing and high-throughput multi-omics analysis facilitated the reconstruction of the metabolic networks. Based on a genome sequence of R. toruloides strain NP11, the first metabolic network of R. toruloides included its central carbon metabolism and lipid biosynthetic pathways [ 10 ]. R. toruloides possesses several enzymatic pathways that differ from the model yeast Saccharomyces cerevisiae and which specifically facilitate the generation of lipid precursors. The key differences included the synthesis of acetyl-CoA from citrate by ATP citrate lyase (ACL), synthesis of acetyl-CoA from xylulose 5-phosphate by phosphoketolase (XPK), and the conversion of S-malate into pyruvate by malic enzyme (ME) that provides for NADPH [ 10 , 14 ]. Proteomics analysis has suggested NADPH regeneration primarily through the pentose phosphate pathway on xylose and glucose, but the role of malic enzyme is not clearly understood [ 8 , 10 ]. The role of XPK in the generation of acetyl-CoA has not been acknowledged previously, whereas ACL has been demonstrated to be upregulated during lipid accumulation [ 10 ], especially in presence of xylose [ 8 ].

R. toruloides is a red basidiomycota known for its ability to accumulate high amounts of intracellular lipids [ 1 ] and consume different carbon substrates [ 2 , 3 ]. It has been studied for its ability to consume complex biomass substrates, including from the lignocellulosic origin [ 4 – 6 ] that would make it interesting for a biorefinery concept. However, studies aimed at fundamental investigation of R. toruloides metabolism have been mainly conducted using a single carbon source as substrate, such as xylose [ 1 , 7 – 9 ], glucose [ 8 – 11 ], glycerol [ 7 ], acetate [ 7 ], L-arabinose and p-coumarate [ 9 ], in a chemically defined mineral medium and occasionally rich cultivation medium (YP) [ 3 ]. Secondary nutrient limitation induces lipid accumulation [ 12 ]. In nitrogen limitation, 65% of lipids of dry cell weight were reached in a batch cultivation regime [ 1 ].

Results

Relation between ribosomal content, growth rate and translation We used absolute quantification of the proteome and the ribosomal content to calculate the rate of protein synthesis per ribosome, also known as ribosome efficiency or protein translation rate (for instructions see S3 Table). The ribosome of R. toruloides strain NP11 was characterized by 178 structurally distinct proteins reported in Uniprot.org, from which 147 were identified in CCT 7815 strain and quantified (S2 Table). The calculated translation rates varied from 0.8 to 6.6 aa/s (Fig 2C), which was very similar as observed in S. cerevisiae (between 2.8 and 10 aa/s) [28]. Among the 6 conditions analyzed, we observed a linear correlation between the translation rate and specific growth rate μ (R2 = 0.99, p-value < 0.001). The mass-wise ribosome content of proteome (g/g_protein) (Fig 2C) had no such distinct correlation with the μ (R2 = 0.68, p-value = 0.043). Interestingly, the lowest ribosome content in proteome was detected during growth on acetate as compared to other substrates.

Integrating fluxomic and proteomic analysis using an enzyme-constrained genome-scale model Genome-scale models allow an in silico simulation of intracellular flux patterns in accordance with exchange fluxes obtained from cultivation experiments. To improve the predictive power and consider the capacity constraints imposed by enzymatic catalytic capacities and their protein levels, we developed an enzyme-constrained GEM using the GECKO Toolbox [17]. In lieu of a strain-specific model, we used the NP11-based GEM [15] to represent the CCT 7815 strain used in this study. The genome of its parental strain CCT 0783 possesses two versions of the same gene, one presenting >90% identity and the other version presenting >70% identity to the genome of haploid strain NP11 [29]. We integrated individual protein concentrations with their corresponding catalytic activities (k cat ) in the model to constrain individual metabolic fluxes. We created separate models for exp and Nlim phases on xylose (X), glucose (G) and acetate (A), respectively. Hence, 6 different versions of the proteome constrained model with modified biomass composition, fatty acid profiles and flux bounds from the experimental data were constructed. Proteome constraints included the concentrations of 773 different enzymes across all conditions (S4A Fig and S4 Table), which were applied to 1515 metabolic reactions (30% of all reactions) (S5 Table). The coverage of these constraints was greatly improved by manually assigning EC numbers to 461 R. toruloides enzymes (S4 Table), which enabled GECKO Toolbox to assign their k cat values. At first, BRENDA was queried for exact matching reaction, substrate and organism. But as kinetic parameter data for non-model organisms such as R. toruloides were not readily available, GECKO Toolbox step-wise relaxes the stringency when matching EC number, organism and substrate, to assign reasonable estimates of k cat values [30]. Mass-wise, the proteome constraints of the measured fraction of enzymes covered between 14% (Gexp) to 25% (ANlim) of the quantified proteome (S4B Fig). Aside from enzyme concentrations, proteome constraints contained 535 unique k cat values automatically queried from the BRENDA database (S3 Dataset). Models, data sets and scripts are hosted on a dedicated Github repository ecRhtoGEM (www.github.com/alinarekena/ecRhtoGEM). Next we used Flux Balance Analysis [31] to simulate intracellular flux patterns and random sampling of the solution space [32] with 2000 sampling iterations to evaluate flux variability (S4 Dataset). To constrain the set of feasible solutions during sampling, we fixed the upper bound and lower bounds on the observed exchange fluxes, ATP hydrolysis (non-growth related maintenance) and protein pool exchange (see Methods). The average flux variability, estimated as a percentage of SD divided by the flux values, was 19% (median value of all conditions) (S4 Dataset). From the simulated flux values, we calculated apparent enzyme catalytic activities (k app ) as ratio of model-predicted fluxes and measured protein concentration. k app represents the apparent in vivo enzyme turnover which drives the biological processes in the environment, in contrast to k cat representing maximum enzyme capacity. As we used model-predicted fluxes that were constrained by k cat values in the ec-model, the k app values that we obtained cannot be higher than the k cat value, which means that we cannot capture any potential in vivo enzyme activity enhancement effect. Regardless, in case of high k app values, high reaction rates are catalyzed by low protein concentration, and vice versa. This study is the first report on the in silico k app values in R. toruloides. Calculated k app values in all growth conditions are available in S3 Dataset. Vast majority of k app values were in the range from 0.1 to 100 (s-1) (S5 Fig), which is in the range of “average enzyme” k cat of 10 s-1 reported by Bar-Even et al. [33]. Some of the lowest k app values in acetate condition were associated with fatty acid degradation and beta oxidation metabolic pathways. We found that during the Nlim phase, when lipid accumulation occurs, the number of enzymes with relatively low k app values (0.1 to 1 s-1) was increased (S5 Fig). This reflects the fact that absolute fluxes decreased more than protein concentrations during the Nlim phase in comparison to exp growth phase, suggesting that for many reactions downregulation of the enzyme did not affect its reaction rate directly.

Growth on xylose Next, we explored R. toruloides metabolism during growth on xylose. Xylose is metabolized by xylose reductase (XR, NADPH-dependent), which reduces xylose to xylitol, xylitol dehydrogenase (XDH) and xylulokinase (XK), and further assimilated into central carbon metabolism via transketolase (TKT1-2) or XPK pathway. The expression of XK was not detected on proteome level in any of the conditions studied, suggesting an alternative pathway to the known fungal xylose assimilation pathway (Fig 3). The experimental detection of D-arabinitol isoform suggested the conversion of D-xylulose to D-arabinitol. This mechanism was supported by the presence of two genes in the R. toruloides genome encoding D-arabinitol dehydrogenase, RHTO_07844 and RHTO_07702. Only protein RHTO_07844 was detected in our proteomics analysis (1913 μg/g_protein) (S10 Fig), suggesting its role as D-arabinitol 4-dehydrogenase (DAD-4), converting D-xylulose to D-arabinitol. Arabinitol dehydrogenase could also convert arabinitol to ribulose (D-arabinitol 2-dehydrogenase) [34]. L-xylulose reductase (LXR) of fungal A. monospora has been reported to reversibly convert D-ribulose to D-arabinitol [35]. In support of this mechanism, protein levels of L-xylulose reductase (EC 1.1.1.10, RHTO_00373) were 10-fold upregulated during growth on xylose versus other substrates. Therefore, RHTO_00373 was selected as D-arabinitol 2-dehydrogenase (DAD-2) (converting D-arabinitol to D-ribulose). Arabinitol dehydrogenase is known to use NADH as cofactor [34]. LXR is mostly known for NADP(+)/NADPH specificity [35]. D-ribulose can enter the non-oxidative part of PPP via phosphorylation by D-ribulokinase (RK). An equivalent pathway was recently reported by Kim et al. 2021 [9]. One gene in R. toruloides IFO 0880 GEM (version 4.0) was annotated as D-ribulokinase (ID 14368) and we used it to identify potential RK in NP11 strain, which is more similar to the strain CCT 7815 used in this study [29]. Gene RHTO_00950 was identified as an ortholog of protein ID 14368 by a BLAST search, which found a match with 98.5% identity. Interestingly, in strain IFO 0880 orthologs of both genes RHTO_07844 and RHTO_07702 were identified as DAD-2 and DAD-4, respectively, and both were using NAD/NADH as the cofactor [9]. While it is not known, which cofactor of DAD-2/LXR enzyme is operational in R. toruloides strain CCT 7815, we analyzed the fluxes of both possible scenarios. Both simulation results were highly similar, with a difference in where NADPH was regenerated. The alternative pathway through DAD was preferred even when XK was not constrained with proteome. In a scenario when DAD-2/LXR was NADP-dependent, during both exp and Nlim growth phases between 46–49% of carbon derived from xylose was directed via glucose 6-phosphate isomerase (GPI) in a reverse direction to the glycolytic flux. In a combination with that, 42% of carbon was directed via the oxPPP and returned to the Ru5P branching point, indicating that a loop associated with NADPH recycling is taking place. Alternatively, up to 88% of xylose-derived carbon was directed via oxPPP (S7 Dataset). In the first scenario, ZWF and GND provided more NADPH than LXR/DAD-2 during the exp phase (S12 Fig). During the Nlim phase, when the yield of NADPH slightly increased (S13 Fig), the flux of ZWF and GND remained unchanged, while the flux of LXR/DAD-2 increased (Fig 4A). XR consumed at least 2-fold more NADPH than any other NADP(+)-dependent enzyme during both growth phases. However during Nlim, more NADPH consumed by FAS1-2 was spent on lipid biosynthesis (S12 Fig). From the proteomics analysis, the concentrations of enzymes involved in the xylose pathway were 1.1 to 1.6-fold downregulated during Nlim phase versus exp phase (Fig 4B), consistent with the decrease in xylose uptake rate (S1 Table). Lower concentration of RK (644±8 μg/g_protein) compared to other enzymes involved in xylose assimilation was measured (S10 Fig), suggesting enzyme limitation in the XK bypass pathway. At relatively low protein levels, high flux through RK resulted in relatively higher k app values (S11 Fig). Aside from enzymes directly involved in xylose assimilation, the intracellular flux patterns on xylose were the closest to growth on glucose, in comparison to growth on acetate. The flux of XPK was upregulated 1.7-fold between 13% to 22% during the Nlim phase (Fig 4A), which was also the main source of acetyl-CoA during lipogenesis. At the Nlim phase, the yields of ATP and NADH significantly increased (S8 and S9 Figs). The additional mitochondrial NADH during the Nlim phase was provided via internal cycling of MDH1 (S14 Fig).

Growth on acetate Lastly, we explored R. toruloides metabolism during growth on acetate. Acetate can cross the plasma membrane to enter the cells via simple or facilitated diffusion, but at pH below neutral (< pH6) the diffusion of the undissociated form of the acid induces the stress response or causes negative effect on metabolic activity [36]. In R. toruloides, two permeases have been found upregulated during growth on acetate-based rich medium in comparison to glucose-based rich medium [3], suggesting that facilitated diffusion is taking place. Once inside the cells, acetate is assimilated via ACS that directly provides acetyl-CoA (Fig 3), one of the main precursors for lipid biosynthesis. From acetyl-CoA branching point, the flux is channeled into the central metabolic pathways via isocitrate lyase (ICL1-2) and malate synthase (MLS), which are predicted to be located in cytosol, but no experimental evidence is available. Metabolic model predicted that at acetyl-CoA branching point, 18% of carbon from acetate during exp growth phase was directed to lipid biosynthesis via ACC, while the majority of carbon (51%) entered glyoxylate shunt. In addition, a significant amount of carbon from acetyl-CoA (29%) was channeled via CRC carrier, which was predicted to have a minor activity on glucose condition. The CRC route was preferred over the PDH pathway towards mitochondrial acetyl-CoA (MLS-ME-PDH), which channeled only 18% of carbon from acetate at exp phase. Metabolic model predicted 5% of carbon from acetate excreted as succinate from the glyoxylate shunt, in addition to 2% of carbon excreted as citrate, which was confirmed by HPLC. During the Nlim phase, the main fluxes demonstrated different regulation (Fig 4A). The increase in flux via CRC (1.4-fold) reflects that more carbon entered the TCA cycle during the Nlim phase. Interestingly, the flux of ACC was downregulated 3.1-fold at Nlim compared to the exp growth phase. Using the rate of lipid production, which decreased during Nlim phase, it can be explained that the lipid production in absolute amounts was higher during the exp phase to sustain the growth together with moderate lipid production (S1 Table). On acetate, fluxes of the TCA cycle were the lowest, while measured protein levels were the highest among all conditions analyzed (Fig 4A and 4B). At Nlim phase, 28% of carbon from acetate was predicted to be excreted as OAA, while the biomass yield decreased (S1 Table). Flux levels of the TCA cycle indicated that an internal cycling of carbon similar as in other conditions was taking place (S6 Fig). It involved different transporter proteins,—the citrate-oxoglutarate (CTP) and succinate-fumarate (SFC) transport -, which allow channeling of the flux from the TCA cycle to glyoxylate shunt (Fig 3). During the Nlim phase, ATP turnover, produced entirely via ETC (S15 Fig), and NADH turnover, produced almost entirely via the TCA cycle (S16 Fig), both decreased (S8 and S9 Figs), unlike observed in glucose or xylose conditions, where ca 80% of the ATP originated from ETC, while the rest came mainly from glycolysis. Aside from enzymes directly involved in the TCA cycle, cytosolic ME was the sole supplier of NADPH during the Nlim phase only in acetate condition (S17 Fig). This is supported also by the measured protein levels of ME, which were significantly higher under acetate conditions, although the absolute levels of ME were relatively low under all studied conditions (189±3 μg/g_protein) (Figs 4B and S6). Only during the growth on acetate the NADPH yield decreased during the Nlim phase (S13 Fig). During the Nlim phase, more NADPH was consumed by FAS and spent on lipid biosynthesis (S17 Fig). We also observed few significant changes in fluxes of enzymes involved in gluconeogenesis, which is an important pathway during growth on acetate to provide xylose phosphate-based precursors for ribonucleotide synthesis. The normalized flux towards gluconeogenesis, channeled via MDH2, carried 11% of carbon exp growth phase. It may reflect the fact that PEPCK, the first enzyme in the gluconeogenesis pathway, consumed ATP, but we found that PEPCK was consuming only 2.4% of ATP during exp phase (S15 Fig). The concentration of PEPCK (627 μg/g_protein) (S10 Fig) and its k cat value (38 s-1) (S2 Dataset) were low, suggesting that PEPCK could have been a rate-limiting step of gluconeogenesis during the exp phase. From proteomics analysis, the concentration of enzymes involved in fatty acid beta oxidation (RHTO_04957, RHTO_00300, RHTO_02848, RHTO_07118, RHTO_00476) at higher concentrations (from 163±9 to 531±7 μg/g_protein) as compared to cells grown on other substrates at exp phase (S2 Dataset), suggesting this pathway might be more active in R. toruloides during growth on acetate.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011009

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/

You are viewing proxied material from magical.fish. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.