Metabolic traits specific for lipid-overproducing strain of Mucor circinelloides WJ11 identified by genome-scale modeling approach

Introduction

Due to the world population growth and the industrial revolution, the value addition of agricultural materials and residues has been attributed to the sustainable production of biobased products. The prices of petroleum-derived products have gradually increased, and concern about the depletion of crude oil reservoirs has stimulated the use of alternative sources to replace the fossil-based products (Zinoviev et al., 2010; Williams et al., 2009; Dellomonaco, Fava & Gonzalez, 2010). In particular, the dedicated biobased oleochemicals and functional lipids with industrial applications (Salimon, Salih & Yousif, 2012; Vanhercke et al., 2013; Hatti-Kaul et al., 2007) are of a great interest. Oleaginous microorganisms are currently of interest as cell chassis for pathway manipulation because they exhibit phenotypes with advantages in the cultivation process and lipid production. Indeed, the particular microbes have been used in several industrial sectors, including feed, food, and biofuels (Thevenieau & Nicaud, 2013; Mata, Martins & Caetano, 2010). Of these, Mucor circinelloides is a well-known oleaginous fungus that accumulates lipids at high levels, particularly under certain culture conditions (Xia et al., 2011). The main storage forms in M. circinelloides are lipid bodies (LBs) or lipid particles (LPs), in which triacylglycerol (TAGLY) is predominant. Moreover, M. circinelloides is capable of synthesizing the nutritionally important polyunsaturated fatty acid (PUFA), γ-linolenic acid (GLA, C18:3 n-6; cis 6, 9, 12-octadecatrienoic acid), which has beneficial effects on human and animal health (Zurier, 1998).

The genome-scale metabolic model (GEM) is one of the computational tools used to predict metabolic behaviors and it is usually employed in systems biology in conjunction with other modern technologies, such as gene editing and synthetic biology (Thiele & Palsson, 2010). Of the oleaginous strains, the GEM of the yeast Yarrowia lipolytica was reconstructed, aiming to enhance the production of biodiesels and other valuable products (Pan & Hua, 2012). In oleaginous fungi, the reconstruction of the GEM of Mortierella alpina was implemented to investigate the metabolic characteristics for enhancing the production of lipids rich in arachidonic acid (Ye et al., 2015). Furthermore, Vongsangnak et al. (2016) also demonstrated the empowering of GEM for dissecting the growth behavior of M. circinelloides strain CBS 277.49 on various nutrient sources through the comparative analysis of the three GEMs of M. circinelloides (iWV1213), M. alpina (iCY1106) (Ye et al., 2015) and Y. lipolytica (iYL619_PCP) (Pan & Hua, 2012).

Among the oleaginous strains, M. circinelloides strain WJ11 has been identified as a promising strain for the overproduction of lipid-derived products based on its ability to accumulate lipids up to 36% of its dry cell weight (DCW) under nitrogen-deficiency conditions, which is higher than that has been reported in the reference strain CBS 277.49 (15% lipid of DCW) (Tang et al., 2017). As a consequence, the cellular metabolic mechanisms governing the oleaginicity of this particular strain have been addressed through several approaches, such as multilevel omics analysis (Tang et al., 2016; Tang et al., 2015a; Tang et al., 2015b; Tang et al., 2017). Recently, the constructed genome-scale metabolic network of M. circinelloides WJ11 was employed to explain its metabolic routes, focusing on lipid metabolism and carotenoid biosynthesis (Vongsangnak et al., 2018). To gain a more precise scaffold for the predictive analysis of the metabolic control involved in the oleaginicity of the WJ11 strain, a functional GEM of M. circinelloides WJ11 was developed in this work. Briefly, the metabolic network of M. circinelloides WJ11 (Vongsangnak et al., 2018) was initially used as a scaffold for improving gene annotation through metabolic reconstruction. For metabolic modeling and analysis, flux balance analysis (FBA), phenotype phase plane analysis (PhPP), fast flux variability analysis (fastFVA), and uniform sampling with coordinate hit-and-run with rounding (CHRR) were afterwards executed for observing cellular phenotypes of M. circinelloides WJ11. Finally, a comparison of the models between the lipid-overproducing strain WJ11 and the reference strain CBS277.49 of M. circinelloides was performed. This work demonstrated the efficiency and reliability of the GEM in describing growth behavior and specific metabolic traits of lipid overproduction in the M. circinelloides strain WJ11, which is applicable in pathway manipulation for cell optimization relevant to desired products with industrial applications.

Materials & Methods

Results

Growth simulation and validation of iNI1159

Using FBA through the MATLAB environment (see ‘Materials & Methods’), the specific growth rates of iNI1159 on a glucose-containing medium were predicted. When the glucose uptake rates were fixed at 2.518 (Zhao et al., 2015) and 1.807 (Tang et al., 2016) mmol gDW⁻¹ h⁻¹ and the other environmental conditions, e.g., oxygen and nitrogen uptake rates were unconstrained, the results showed that the specific growth rates were consistent with the experimental data, in which the deviation values were 0.05% and 9.19%, respectively as shown in Table 1 (see File S5). These results suggest that the model is efficient for further analysis using the PhPP, fastFVA and CHRR approaches. The effects of nitrogen and glucose uptake rates on the cell growth and lipid production rates of M. circinelloides WJ11 were then investigated by PhPP analysis. As shown in Fig. 1A, the glucose uptake rate, rather than the nitrogen uptake rate, affected the biomass production rate. Although the nitrogen uptake rate was important for the biomass production rate, a low biomass growth rate was obtained under conditions with a high nitrogen uptake rate. In other words, this result suggested that the growth rate was linearly dependent on the glucose uptake rate when the nitrogen source was in excess. However, a progressive increase in the glucose uptake could not enhance the biomass growth rate if the nitrogen was limited. These results further indicated that the ratio of carbon and nitrogen contents might contribute to the growth rate. In addition, the results of the PhPP analysis indicated the effect of nitrogen and glucose uptake rates on the lipid production of WJ11 as shown in Fig. 1B. When the TAGLY production rate was set as the objective function and the biomass growth rate was fixed at 0.1671 h⁻¹, the TAGLY production rate was directly proportional to the glucose uptake rate, but inversely proportional to the nitrogen uptake rate. This means that high glucose and low nitrogen conditions resulted in enhanced lipid accumulation as previously described (Tang et al., 2016). This suggests that the lipid production rate was induced by an increase in the glucose uptake rate and a decrease in the nitrogen uptake rate, which is the cultivation condition for triggering the lipid accumulation in oleaginous strains (Botham & Ratledge, 1979; Ratledge & Wynn, 2002; Ratledge, 2004).

Major metabolic reactions associated with lipid accumulation under nitrogen-limited conditions

To identify major metabolic reactions possibly responsible for increased lipid production behavior under nitrogen-limited conditions, fastFVA followed by relative flux change analysis was applied (see ‘Materials & Methods’). Table 3 shows that the top 20 metabolic reactions with high values of relative flux change in the metabolism of carbohydrates, lipids and certain amino acids. The values of relative flux change for all metabolic reactions of iNI1159 when nitrogen-limited conditions are given in File S6. For carbohydrate metabolism, most of the metabolic reactions were in the glycolysis pathway, such as pyruvate oxidation catalyzed by the pyruvate dehydrogenase complex (EC: 1.2.4.1: 2.3.1.12: 1.8.1.4), and the conversion reaction catalyzed by phosphoglycerate kinase to yield D-glyceraldehyde 3-phosphate (EC: 2.7.2.3). Moreover, the flux of a metabolic reaction in the pentose phosphate pathway, which was catalyzed by ribulose-phosphate 3-epimerase (EC: 5.1.3.1) for conversion of intermediate substance from D-ribulose 5-phosphate to D-xylulose 5-phosphate, was also altered. For lipid metabolism, it was found that some reactions responsible for fatty acyl supply were markedly influent, which included the reaction catalyzed by acetyl-CoA carboxylase (EC: 6.4.1.2) for converting acetyl-CoA to malonyl-CoA, the oxidation reaction of acetoacetyl-CoA to yield acetyl-CoA catalyzed by acetoacetyl-CoA thiolase (EC: 2.3.1.9), and the reaction in biosynthesis of lanosterol to produce acetyl-CoA catalyzed by hydroxymethylglutaryl CoA synthase (EC: 2.3.3.10). Moreover, the reaction in nitrogen metabolism to convert carbon dioxide into carbonic acid catalyzed by carbonic anhydrase (EC: 4.2.1.1) was observed. Interestingly, it was detected that some amino acid reactions with relative flux change were the reaction catalyzed by threonine dehydratase to pyruvate and ammonia (EC: 4.3.1.19), the reaction catalyzed by cystathionine beta-synthase (EC: 4.2.1.22) to produce L-cystathionine, the reaction catalyzed by cystathionine beta-lyase (EC: 4.4.1.8) to generate pyruvate and ammonia, the reaction catalyzed by alanine aminotransferase to yield alanine (EC: 2.6.1.2), and the reaction responsible for converting L-serine to pyruvate and ammonia that is catalyzed by L-serine dehydratase (EC: 4.3.1.17).

Table 3:

List of top 20 metabolic reactions with high flux change in relation to the lipid production identified in iNI1159 model.

EC number	Metabolic reaction	Subsystem	Relative flux change
4.1.1.1	PYR[c] ->CO2[c] + ACAL[c]	Pyruvate metabolism	21.08
1.2.4.1: 2.3.1.12: 1.8.1.4	NAD[m] + PYR[m] + COA[m] ->NADH[m] + ACCOA[m] + CO2[m]	Glycolysis	21.08
6.4.1.2	ATP[c] + HCO3[c] + ACCOA[c] ->ADP[c] + PI[c] + MALCOA[c]	Fatty acid biosynthesis	20.16
2.7.1.40	ADP[c] + PEP[c] ->ATP[c] + PYR[c]	Glycolysis	15.56
1.2.1.12	T3P1[c] + PI[c] + NAD[c] <=>13PDG[c] + NADH[c]	Glycolysis	15.56
2.7.2.3	ADP[c] + 13PDG[c] <=>ATP[c] + 3PG[c]	Glycolysis	15.56
5.4.2.1	3PG[c] <=>2PG[c]	Glycolysis	15.56
4.2.1.11	2PG[c] <=>PEP[c] + H2O[c]	Glycolysis	15.56
4.2.1.1	CO2[c] + H2O[c] ->HCO3[c] + H_PO[c]	Nitrogen metabolism	14.64
2.3.3.10	COA[c] + H3MCOA[c] <=>ACCOA[c] + AACCOA[c]	Biosynthesis of lanosterol	14.34
2.3.1.9	COA[c] + AACCOA[c] <=>2 ACCOA[c]	Fatty acid oxidation	14.34
4.3.1.19	SER[c] <=>PYR[c] + NH3[c]	Glycine, serine and threonine metabolism	11.04
4.2.1.22	SER[m] + HCYS[m] ->LLCT[m]	Glycine, serine and threonine metabolism	11.04
4.4.1.8	H2O[m] + LLCT[m] ->PYR[m] + NH3[m] + HCYS[m]	Methionine metabolism	11.04
4.1.2.13	FDP[c] <=>T3P2[c] + T3P1[c]	Glycolysis	11.04
2.6.1.2	PYR[c] + GLU[c] <=>AKG[c] + ALA[c]	Alanine/aspartate and asparagine metabolism	11.04
5.1.3.1	RL5P[c] <=>XUL5P[c]	Pentose phosphate pathway	11.04
2.7.1.11	ATP[c] + F6P[c] ->ADP[c] + FDP[c]	Glycolysis	11.04
4.3.1.17	SER[c] ->PYR[c] + NH3[c]	Cysteine metabolism	11.04
5.3.1.1	T3P2[c] <=>T3P1[c]	Glycolysis	10.04

DOI: 10.7717/peerj.7015/table-3

Comparative analysis between iNI1159 and iWV1213 for guiding key metabolic routes involved in lipid accumulation

Because of the colinear evolution between the genomes of the lipid-overproducing and reference strains (Tang et al., 2015a), iNI1159 and iWV1213 undoubtedly shared almost common metabolic reactions in their metabolic networks as shown in Fig. 2. Interestingly, it was observed that the numbers of reactions in carbohydrate, amino acid and cofactor metabolisms between iNI1159 and iWV1213 were different. To gain a better understanding of the metabolic behaviors between the WJ11 and reference CBS277.49 models, the flux distributions in both models were determined and compared based on fastFVA and CHRR. Upon setting the glucose uptake rate to 2.518 mmol gDW ⁻¹h⁻¹ and biomass growth as the objective function, FBA was performed and showed that the specific growth rates of iNI1159 and iWV1213 were 0.1671 h ⁻¹(see Table 1) and 0.2593 h⁻¹ (see File S5), respectively. All solutions of the exchange reactions obtained from FBA were then fixed to perform either fastFVA or CHRR for determining flux distributions. The possible range of solutions obtained from fastFVA is shown in File S7. The fastFVA solutions were then divided into three categories, including “always-active reactions”, “sometimes-active reactions”, and “never-active reactions”. The reactions that had the minimum and maximum flux of nonzero values with the same sign were referred to as “always-active reactions”. The reactions with the minimum and maximum flux values spanning zero were called “sometimes-active reactions”. The “never-active reactions” were the reactions with both the minimum and maximum fluxes equal to zero. As a result, Fig. 3 shows that iNI1159 contained 252 always-active reactions (19% of total reactions), 152 sometimes-active reactions (11% of total reactions), and 951 never-active reactions (70% of total reactions). For iWV1213, 240 always-active reactions (18% of total reactions), 175 sometimes- active reactions (13% of total reactions), and 912 never-active reactions (69% of total reactions) were observed. Detailed information is provided in File S8 and File S9. As expected, most common always-active reactions shared between iNI1159 and iWV1213 were the reactions catalyzed by hexokinase (EC: 2.7.1.1) through pyruvate kinase (EC: 2.7.1.40), and some reactions in pentose phosphate pathway catalyzed by ribulose-phosphate 3-epimerase (EC: 5.1.3.1) through phosphogluconate dehydrogenase (EC: 1.1.1.44). The most common sometimes-active reactions were involved in the tricarboxylic acid cycle. However, the metabolic fluxes between these two strains were different, particularly in lipid metabolism, including a set of reactions catalyzed by the enzymes involved in glycerophospholipid and glycerolipid metabolism, such as lysophosphatidic acid acyltransferase (EC: 2.3.1.51), phosphatidate phosphatase (EC: 3.1.3.4), and the enzymes involved fatty acids biosynthesis (acetyl-CoA carboxylase (EC: 6.4.1.2), fatty-acyl-CoA synthase (EC: 2.3.1.86), and fatty acid elongase (EC: 2.3.1.199), as shown in Fig. 4. As a result, the iNI1159 had the unique always-active reactions in lipid metabolism, whereas iWV1213 had the unique sometimes-active reactions in lipid metabolism (File S10). Moreover, it was found that the flux distributions over the central metabolic pathways between iNI1159 and iWV1213 are different indicated by dotted arrows (Jaccard index = 0) as shown in Fig. 5. Besides, the average flux values obtained from CHRR are shown, i.e., top figures representing fluxes of iNI1159 and bottom figures representing fluxes of iWV1213 (Fig. 5). The fluxes between the two strains were found to be significantly different according to the Student’s t-test (File S11). It was noted that iNI1159 possessed higher metabolic fluxes in lipid metabolism, such as lanosterol, zymosterol, glycerolipid and fatty acids biosynthesis than iWV1213. In contrast, the iWV1213 showed higher flux distribution in the metabolism of carbohydrates (i.e., pentose phosphate pathway) and amino acids (arginine, cysteine, glycine and threonine metabolism), which were relevant to cell growth or biomass production. The detailed information of Jaccard index, average fluxes and their standard deviations as well as statistical differences (p-value) for the entire networks of the two models are provided in File S11.

Discussion

GEM is an effective prediction tool by using constraint-based analysis integrated with linear programming algorithm. In this work, the developed GEM of M. circinelloides WJ11 expanded the metabolic information in the lipid biosynthetic pathway by incorporating the TAGLY production and exchange reaction (Table 2), which were a representative of lipids accumulated in LPs. The iNI1159 model could predict the simulated growth rate with low percent deviation (Table 1 and Fig. 1A). Interestingly, the lipid-accumulation stage of oleaginous strains, which is generally triggered by the carbon-excess and nitrogen-limited condition (Lamers et al., 2016), could be defined by the model (Fig. 1B). These results suggested that the developed model could be used for investigating the metabolic flux at particular conditions. As a consequence, the metabolic reactions affecting the lipid-accumulation process were identified by simulating the culture condition, in which the carbon source was fixed, and the nitrogen source was a variable or in turn the C/N ratios were varied. The top 20 metabolic reactions with high flux change in relation to the lipid production flux were found in iNI1159 (Table 3) when simulating the culture condition under the glucose uptake rate (2.518 mmol gDW⁻¹ h⁻¹) and nitrogen uptake rate (2–10 mmol gDW⁻¹ h⁻¹) (see File S6). Of them, some metabolic reactions have been previously identified as the most connected metabolites in the oleaginous strains, Y. lipolytica and M. alpina, such as precursors (i.e., acetyl-CoA) in lipid biosynthesis and current metabolites (i.e., ATP and NADPH) (Pan & Hua, 2012; Vongsangnak et al., 2013). Particularly, the increased fluxes of some reactions in the central carbon metabolism of iNI1159, such as glycolysis and pentose phosphate pathway might promote the supply of acetyl-CoA and NADPH, which are basic precursors for the lipid biosynthesis similarly to the previous findings in GEMs of the oleaginous yeast (Kerkhoven et al., 2016). These results are also coincided with the previous report of proteomic study in this oleaginous fungus, in which most enzymes in the glycolysis of M. circinelloides WJ11 were up-regulated under nitrogen exhaustion (Tang et al., 2016). As a result of the uniform random sampling analysis in i NI1159 (Fig. 5), the increased fluxes in acetyl-CoA generation of WJ11 strain by the catalytic functions of pyruvate dehydrogenase, citrate synthase, ATP:citrate lyase (ACL) and acetyl-CoA synthetase were further linked to significant increases of the fluxes toward the biosynthesis of lipids, particularly TAGLY, squalene (SQL) and steryl ester (ERGOSE or SE). These findings could explain the high lipid accumulation in M. circinelloides strain WJ11 by the fact that TAGLY and SE are main components of storage lipids in form of lipid body (LB) or lipid particle (LP). Recently, the efficient strategy for rational improvement of SQL production in Y. lipolytica by engineering the acetyl-CoA metabolism based on genome-scale metabolic model has also been reported (Huang et al., 2018).

Not surprisingly, all possible metabolic reactions in lipid metabolism, e.g., fatty acids biosynthesis and fatty acid oxidation, were among the top 20 metabolic reactions with high flux change in relation to the lipid production rate. Not only the carbohydrate and lipid metabolisms, but some reactions in amino acid metabolism appeared to contribute to lipid overproduction under nitrogen depletion. Most of them were involved in pyruvate and ammonia production (Table 3). Of them, some amino acid pathways might be targeted for genetic engineering to improve the pyruvate production, which is an important precursor for fatty acids biosynthesis. It has been reported that the redirecting carbon flux from amino acids to lipids was found in Y. lipolytica during nitrogen limitation by integrative analysis of multilevel omics data through the iYali4 GEM (Kerkhoven et al., 2016).

Regarding the fastFVA, some metabolic reactions were differently displayed even though iNI1159 was developed using iWV1213 as a draft network. The common always-active reactions found in the glycolysis and pentose phosphate pathway of the two models indicated a core metabolic route for generating acetyl-CoA and NADPH, which are key precursors for fatty acids biosynthesis in common eukaryotic cells. It is noteworthy that the sometimes-active reactions in tricarboxylic acid cycle also shared among these models, which are similar to the previous findings (Kavšček et al., 2015; Vongsangnak et al., 2016). Interestingly, fastFVA indicated that fluxes in some lipid metabolic reactions categorizing either sometimes-active or always-active were different between these two models, such as glycerophospholipid and glycerolipid metabolism, as well as fatty acids biosynthesis. It could be implied that iNI1159 provided high lipid accumulation, whereas iWV1213 provided low lipid accumulation due to more sometimes-active in lipid metabolism (Fig. 4). These findings recommend that sometimes-active reactions should be further focused in more details due to these metabolic reactions may associate to both growth and lipid accumulation phases.

When compared the flux distributions between the lipid-overproducing and reference strains, both fastFVA and CHRR provided similar results. The Jaccard index based on fastFVA as well as mean fluxes obtained from CHRR indicated that i NI1159 had higher fluxes in the glycolysis and lipid metabolism while iWV1213 showed higher fluxes in the pentose phosphate, carbohydrate and amino acid metabolisms (Fig. 5). The achieved results can possibly be explained by the difference in stoichiometric coefficients of biomass equations between these two models (Table 4). Comparatively, iWV1213 had higher values of stoichiometric coefficients for protein leading to greater fluxes in amino acid metabolism, whereas iNI1159 had higher values of stoichiometric coefficients for lipid leading to greater fluxes in lipid metabolism. Notably, changes in the metabolic fluxes across amino acid metabolism of these two models were clearly observed. For examples, the metabolic fluxes of aspartate, glycine and cysteine metabolisms of iWV1213 were higher than those of iNI1159. Possibly, these might be alternative metabolic routes for leveraging between biomass and lipid production.

Taken together, the developed model of iNI1159 in this study was supported by the related experimental results, particularly in terms of physiological responses during lipid accumulation stage, where low growth rate and high lipid production rate was observed (Zhao et al., 2015; Tang et al., 2016). As a result, the iNI1159 model showed a tendency to produce cellular lipid at high level, whereas iWV1213 displayed the high biomass production (Fig. 5).

Conclusions

The developed iNI1159 empowering the prediction of biomass and lipid production of M. circinelloides WJ11 could be exploited to explain its metabolic phenotypes through changes in metabolic flux distribution. The metabolic reactions in amino acid metabolisms influencing the lipid-accumulation process in WJ11 were identified in the nitrogen-limited conditions. The high relative flux change in carbohydrate, lipid, and amino acids metabolisms identified in the lipid-overproducing WJ11 strain by fastFVA might contribute to overflow of metabolic fluxes to lipid accumulation. Comparative metabolic flux distributions using Jaccard index and uniform random sampling clearly demonstrated distinct flux flows in the central metabolic pathways between the lipid-overproducing strain WJ11 (iNI1159) and reference strain CBS277.49 (iWV1213), enabling analysis of metabolic control involved in the oleaginicity. Thus, this iNI1159 model offers a promising scaffold for metabolic engineering at precise targets to achieve an optimized strain for industrial purpose.

Supplemental Information