Experimentally Validated Reconstruction and Analysis of a Genome-Scale Metabolic Model of an Anaerobic Neocallimastigomycota Fungus

The genome of N. lanati is rich in carbohydrate-active enzymes and metabolically similar to other anaerobic gut fungi.Given the large repeat-rich genomes inherent to anaerobic gut fungi (12), PacBio sequencing was used to obtain a high-quality genome of the isolate N. lanati; this isolate was sourced from a fecal pellet of a sheep obtained from an enclosure at the Santa Barbara Zoo (the Index Fungorum identification number is IF557810). While the genome of this fungus is large (Table 1), it is the second least fragmented of all 5 published gut fungal genomes (see Table S1 in the supplemental material). The N. lanati genome encodes a rich array of carbohydrate-active enzymes (CAZymes) in numbers similar to those reported from other gut fungal genomes (Table S1) (3, 4, 7). In total, 1,788 CAZymes were identified in the genome, of which, 1,253 were expressed in the transcriptome. Like other anaerobic gut fungi, N. lanati deploys both complexed (cellulosomes) and uncomplexed CAZymes through its rhizoidal network (Table 1; Fig. 1), both of which contribute to the decomposition of lignocellulose. Figure S1 shows the breakdown of CAZyme domains identified in the genome of N. lanati. This, in combination with its relatively high growth rate on defined M2 media, suggests that N. lanati is a good model anaerobic gut fungus.

• Open in new tab
• Download powerpoint

FIG 1

The morphology of Neocallimastix lanati aids in the decomposition of unpretreated lignocellulose by disrupting the lignocellulosic plant biomass to increase the surface area available for enzymatic attack. A micrograph of a mature N. lanati sporangium growing on corn stover in M2 medium after 3 days of growth at 39°C. The filamentous rhizoidal network is used to increase the surface area for its lignocellulolytic enzymes that decompose the lignocellulosic corn stover into its fermentable sugar constituents.

Despite advances in sequencing and annotation, a large number of putative gut fungal genes remain completely unannotated (∼48% of the 27,677 predicted genes of N. lanati) (Table S1), which is consistent with previous genomic annotations in this clade. These unannotated genes contributed to gaps found in the draft reconstructed metabolism of N. lanati. A comparative genomic analysis within the primary metabolism across all high-quality publicly available gut fungal genomes (Anaeromyces robustus, Neocallimastix californiae, Pecoramyces ruminatium, and Piromyces finnis) revealed that the gut fungi are metabolically similar. Of the 1,023 unique EC numbers identified across these 5 genomes, fewer than ∼3% are unique to each isolate (Fig. 2). This suggests that gut fungi share a similar primary metabolism. Thus, some metabolic gaps in the draft reconstruction could be filled by searching for genes in N. lanati that are homologous to those encoded, and annotated, in the genomes of the other gut fungal isolates. In this way, key enzymes in the biosynthesis pathways of arginine, asparagine, biotin, riboflavin, lipids, and fatty acids were identified and included in the metabolic reconstruction of N. lanati. In total, 35 gaps in the primary metabolic pathways were identified and annotated in this manner, as noted in the confidence score and homologous gene annotation fields in the model.

• Open in new tab
• Download powerpoint

FIG 2

Anaerobic gut fungi have very similar genetic metabolic potentials, suggesting that metabolic gaps can be filled by looking for homologous genes found in the other sequenced isolates. Each Venn diagram was generated by inspecting the intersection of the annotated EC numbers contained in the genome of each fungus for each metabolic module. Overlapping regions imply that those isolates share the EC assignments contained in each of the metabolic modules. The EC numbers contained in each module are based the KEGG database (60) (see supplemental data 4 in the iNlan20 GitHub repository available at https://github.com/stelmo/iNlan20 for the list of modules encompassing each Venn diagram), while the EC assignments for each fungus are based on the JGI and bidirectional annotation data as described in Materials and Methods.

The curated metabolic model of N. lanati captures the carbon, amino acid, vitamin, fatty acid, nucleotide, and lipid metabolism.Based on the draft metabolic reconstruction of N. lanati, a manually curated genome-scale metabolic model of N. lanati was built (iNlan20) according to an established protocol for generating high-quality reconstructions (26). The model contains 1,023 reactions, 816 metabolites, and 1,018 genes distributed across 3 compartments (hydrogenosome, cytoplasm, and extracellular space). Where possible, experimental data were used to curate the model. Materials and Methods details specifics on the curation process, as well as experiments used to construct the biomass objective function of the model. Briefly, Table 2 shows the experimentally measured macromolecular components of N. lanati that were used to construct the biomass objective function for the genome-scale metabolic model. Further simplifying assumptions were made to construct the specific biomass objective function used in iNlan20. The insoluble carbohydrate component of the biomass was assumed to be solely chitin (27), and the amino acid composition of the protein component of the biomass was assumed to follow the amino acid distribution of the predicted genes (i.e., the predicted proteome). Similarly, the nucleotide composition was assumed to follow the composition of the genome (for the DNA nucleotides) and the transcriptome (for the RNA nucleotides). The lipid component was assumed to be composed of myristic, palmitic, and stearic acids, which were found to be the major fatty acid components of the lipid fraction of N. lanati (see Fig. S2). The growth-associated and non-growth-associated maintenance (GAM and NGAM, respectively) functions were estimated using experimental data (see Fig. S3).

The model focuses on the primary metabolism but includes CAZymes as generalized cellulase and hemicellulase reactions. Of the 791 metabolic genes included in the model, 216 do not have gene assignments—reflecting the understudied nature of gut fungal metabolism. Despite this, the model is stoichiometrically consistent as well as mass and charged balanced (see Memote report or data set S1 in the iNlan20 GitHub repository available at https://github.com/stelmo/iNlan20). Table S2 explains the confidence rating assigned to the genes associated with each reaction in the model. In the energy-generating pathways, particular attention was paid to modeling the hydrogenosome (a mitochondrion-like organelle that functions completely anaerobically), which is discussed in greater detail in the following sections. More generally, the Embden-Meyerhof-Parnas variant of glycolysis is present in N. lanati as well as pathways for mixed-acid fermentation (succinate, acetate, lactate, formate, and ethanol), which are typically found in anaerobic gut bacteria (28). Interestingly, it was found that N. lanati possesses both the NAD⁺ and NADP⁺ variants of glyceraldehyde-3-phosphate dehydrogenase in glycolysis, with the latter used to conserve energy as NADPH instead of ATP. The pentose phosphate pathway of N. lanati is incomplete, with glucose-6-phosphate dehydrogenase and 6-phosphogluconate missing. These reactions regenerate NADPH and possibly explain the presence of the NADP⁺ variant of glyceraldehyde-3-phosphate dehydrogenase as a compensating mechanism (29, 30). Despite these missing genes, the model is still capable of producing nucleotide precursors. The xylose isomerase pathway is also present in N. lanati, as has been found in other sequenced anaerobic gut fungi (21).

The major components (amino acids, nucleotides, vitamins, fatty acids, and lipids) of the anabolic metabolism of N. lanati were found to be present, in agreement with its ability to grow in M2 medium without these components added. Specifically, the complete biosynthesis pathways for all the proteogenic amino acids and the modeled fatty acids were found. Most of the canonical vitamin and cofactor (vitamin B₅, vitamin B₆, riboflavin, and thiamine) biosynthesis pathways were found to be complete, with the exception of folate, where no synthesis mechanism of 4-aminobenzoate was found. The heme and biotin biosynthesis pathways appeared to be incomplete; however, the latter vitamin was not found to be essential (experimentally) in defined medium, suggesting that the pathway is either poorly annotated or that the fungus does not require it for growth. Finally, the model recapitulates the experimentally observed growth rate in defined medium using only the measured flux of glucose (1.5 mmol/g [dry weight]/h) as an input constraint (flux balance analysis predicted μ = 0.044 h⁻¹ versus an experimentally measured μ = 0.045 ± 0.003 h⁻¹).

iNlan20 includes an expanded model of the hydrogenosomal metabolism.Anaerobic gut fungi possess a variant of the hydrogenosome, with the core set of enzymes that catalyze the conversion of malate and pyruvate to acetate, H₂ and formate already identified, as shown in Fig. 3 (2, 14, 16, 31). However, the metabolic pathway leading to H₂ production is not resolved, with literature suggesting either pyruvate ferredoxin oxidoreductase (PFO) or pyruvate formate lyase (PFL) as possible routes (Fig. 3). Both enzymes were identified in the genome and transcriptome and are thus included in the model of the hydrogenosome.

• Open in new tab
• Download powerpoint

FIG 3

An expanded model of the hydrogenosome is included in the model based on genomic annotation, literature, and predicted localization data (14–16). Core hydrogenosome enzymes are colored in blue, while speculative enzymes are shown in black. PFL, pyruvate formate lyase; PFO, pyruvate ferredoxin oxidoreductase; Ac, acetate; SucCoA, succinyl coenzyme A; CoA, coenzyme A; AcCoA, acetyl coenzyme A; Frdx, ferredoxin.

By combining literature sources, gene annotation, transcriptomic expression, and subcellular localization data, we have included additional pathways in the model of the hydrogenosome for N. lanati (Table 3 and Fig. 3). Of note is the inclusion of a putative ATP synthase (32, 33), which was previously speculated to be present in other anaerobic gut fungal isolates (7, 16). Additionally, we also found evidence that complex 2 of the mitochondrial electron transport chain is present: homologs to all four subunits were found to be expressed and localized to the hydrogenosome (complex 2, subunits A, B, C, and D). We could not find any homologs of the membrane-bound subunits of complex 1 or the ATP synthase in the N. lanati genome, as was also reported previously for other anaerobic gut fungi (34). A possible explanation for this is that many of the membrane-bound subunits of the electron transport chain are encoded by mitochondrial DNA, which the fungal hydrogenosome appears to have lost. Despite this, it is perplexing that no homologs of the membrane-bound subunits of complex 1 were found, since these are used to shuttle electrons between the two complexes in the inner membrane of the mitochondrion, and it remains unclear how complex 2 could function without them. However, homologs of the soluble subunits of complex 1, nuoF and nuoE, are highly expressed relative to the other core enzymes of the hydrogenosome (Table 3). The presence of the soluble subunits, coupled with the absence of the membrane-associated subunits of complex 1, has also been observed in the hydrogenosomes of, e.g., Trichomonas vaginalis (35, 36) and Nyctotherus ovalis (37). This raises two possibilities. First, that N. lanati possesses a proton-pumping mechanism, but the missing membrane-bound genes are unannotated. While the membrane-bound subunits of complex 1 and the ATP synthase appear to be missing, we do find preliminary evidence of a pH gradient inside the hydrogenosome (see Fig. S4). This provides some support for a mechanism, possibly involving complex 1, 2, and the ATP synthase, that makes use of this gradient to generate energy, as found in other hydrogenosomes (32), if these genes are present but unannotated. Second, it is possible that its hydrogenase associates with the identified nuoF-like subunit of complex 1 to form a bifurcating hydrogenase, as has been speculated to occur in the hydrogenosome of T. vaginalis (32, 35). Indeed, we find high-homology sequences in the N. lanati genome to all three of the bifurcating hydrogenase subunits that have been enzymatically characterized in Thermotoga maritima (38) (see data set S2 in the iNlan20 GitHub repository available at https://github.com/stelmo/iNlan20). In this case, the function of the membrane-bound subunits of complex 2 and the soluble subunits of the ATP synthase are unknown. In either case, further experimental work is needed to clarify these questions.

Taken together, our expanded model of the hydrogenosome includes the core enzymes previously reported in other gut fungal species as well as a speculative bifurcating hydrogenase, ATP synthase, and proton-pumping module composed of the complex 1 and complex 2 enzymes identified in the N. lanati genome (similar to what has been found in other H₂-producing mitochondria [32]). Given the speculative nature of the proton-pumping mechanism, the ATP synthase, and the bifurcating hydrogenase, these reactions are constrained to carry zero flux in the base case model. These constraints were modified to investigate the consequences of this extended hydrogenosomal metabolism, as discussed later. Additionally, it was previously suggested that a hydrogen dehydrogenase [NAD(P)⁺ + H₂ ↔ H⁺ + NAD(P)H] operates in the reverse direction in the hydrogenosome (7, 14, 15). Consequently, this hydrogen dehydrogenase simultaneously produces H₂ and prevents the accumulation of NAD(P)H produced by the malic enzyme in the hydrogenosome. However, in this direction, the reaction is energetically very unfavorable, with a ΔG of ≈34 ± 5.9 kJ/mol assuming physiologically realistic conditions. Therefore, the flux bounds of this reaction in the hydrogenosome were set to reflect the assumption that the hydrogen dehydrogenase only carries flux in the forward, energetically favorable, direction.

iNlan20 accurately predicts substrate utilization and in vivo fluxes.The curated model was validated using a combination of growth curve, extracellular metabolite, and metabolic flux analysis (MFA) data. Substrate utilization tests were performed on 36 different carbon sources, focusing on metabolites that are present in the digestive tract environment of the anaerobic gut fungi (Table 4). The qualitative prediction accuracy of the model for the substrate utilization and vitamin essentiality validation tests is 89% (Matthews correlation coefficient, 0.79). Interestingly, despite the apparent presence of a full xylose isomerase pathway, N. lanati did not grow using xylose as its sole carbon source, as has been found in other gut fungi (21). In this case, the model’s predictions were incorrect. It was previously suggested that transport limitations may cause this issue (21). However, a xylose transporter was identified in the genome, suggesting that cellular regulation might explain this discrepancy better. Vitamin essentiality tests were also conducted (Table 4). It was found that both heme and 4-aminobenzoate were essential for growth, in agreement with the model’s predictions. In other gut fungi, heme has also been found to be essential (39), suggesting that its de novo biosynthesis pathway may be absent across the clade. It was found that only cysteine could be used as a sulfur source. However, it is not clear if this is a nutritional requirement, since every other reducing agent tested (Na₂S, 2-mercaptoethanol, and dithiothreitol) appeared to be toxic to the fungus. Since cysteine was used to ensure anaerobicity of the medium, we could not test nitrogen source utilization.

Metabolic flux analysis (MFA) was also used to experimentally verify the predicted intracellular fluxes of the GEM. A 1,2-¹³C-labeled glucose tracer was used in conjunction with a carbon atom transition model built from the N. lanati metabolic reconstruction. For the MFA model, metabolic degeneracy caused by the ability of the hydrogenosome to metabolize both malate and pyruvate resulted in large bounds on the fluxes involving these metabolites. To circumvent this, the MFA model was constrained to only import pyruvate into the hydrogenosome, based on previous observations (14). Extracellular metabolic product measurements (ethanol, formate, H₂, acetate, succinate, and lactate) were also used to constrain the MFA model. This resulted in accurate internal metabolic flux measurements based on a statistically significant fit between measured and simulated proteinogenic amino acid labeling patterns (Fig. 4). These ¹³C measured fluxes were then compared to the fluxes predicted using the GEM with independently measured metabolite flux constraints (Table 5). Parsimonious flux-based analysis (pFBA) was then used to find unique flux predictions. Using these constraints, the coefficient of determination between the pFBA and MFA simulation was found to be 0.98 (linear regression fit P < 0.01) (see Fig. S5). This suggests that the constrained metabolic model accurately predicts the steady-state measured intracellular fluxes of N. lanati.

• Open in new tab
• Download powerpoint

FIG 4

The genome-scale metabolic model accurately predicts the in vivo carbon metabolism of N. lanati. Experimentally determined MFA fluxes and predicted pFBA fluxes (top and bottom, respectively) for glycolysis, the TCA cycle, and the hydrogenosome of N. lanati. Error estimates denote one standard deviation from the reported mean for the MFA measurements. Three serially passaged [1,2-¹³C]glucose tracer experiments, grown in M2 medium at 39°C and harvested during exponential phase, were used to measure the in vivo fluxes (see Materials and Methods for more details and the model).

The core hydrogenosome metabolism uses PFO to produce hydrogen, but PFL carries the most flux.There remains uncertainty regarding the presence of pyruvate ferredoxin oxidoreductase (PFO) and its relative importance in the hydrogenosomal metabolism of anaerobic fungi. Earlier enzymatic characterization of hydrogenosomal proteins in Neocallimastigomycota suggested that PFO is the primary route for H₂ production through an associated ferredoxin hydrogenase, as found in the hydrogenosomes of other organisms (16, 31, 32, 40). However, more recent studies suggest that PFO is either absent or of only marginal importance in the gut fungal hydrogenosomal metabolism (14, 15). These later studies suggest that pyruvate formate lyase (PFL), which was likely acquired through horizontal gene transfer from bacteria (41), is significantly more active than PFO. It has been suggested that hydrogen evolution occurs through a hydrogen dehydrogenase working in an energetically infeasible reverse direction (7, 14). Both PFO and PFL were identified in all published gut fungal genomes as well as in N. lanati (Table 3). The model was used to reconcile the role and relative importance of these two enzymes to hydrogenosome function under steady-state growth conditions.

Due to the reaction stoichiometry of PFL, the molar ratio of formate to acetate and ethanol produced is expected to approach unity (1:1) if PFL is metabolically dominant. This is because ethanol and acetate are only produced from PFL (via acetyl coenzyme A) in the cytosol (15). Since PFO only produces acetate, and not formate, the ratio of formate to acetate and ethanol will not be unity if PFO carries significant metabolic flux. Figure S6A shows that the experimentally measured molar ratios of formate to acetate and ethanol are not significantly different (using the unequal variance test, P < 0.05) from unity for N. lanati. This is in agreement with earlier metabolite measurements for Piromyces sp. strain E2, suggesting that PFL is dominant (14). Figure S6B shows that the unconstrained model predicts a wide range of possible ratios, reflecting the metabolic degeneracy of the carbon metabolism of N. lanati. Since there is no energetic cost associated with using PFO versus PFL (both produce one ATP molecule per pyruvate) (Fig. 3), the model predicts that both could be used to maximize ATP production in the hydrogenosome. However, external metabolite flux measurements show only modest H₂ production (Table 5), suggesting that cellular regulation may play a role in diverting flux to PFL instead of PFO. This can also be seen in the relative expression difference between PFL and PFO (an order of magnitude difference between them) in Table 3. When the model is constrained by the measured metabolite fluxes shown in Table 5, the range of possible ratios is reduced to those observed experimentally (Fig. S6B). Since PFO is the only (known) energetically feasible way to produce H₂, this result is not surprising. Using this constraint, the model suggests that PFL carries the most flux in the hydrogenosome, but that PFO is used to produce H₂.

Electron bifurcation and proton pumping may form part of the hydrogenosomal metabolism.Electron bifurcation is an energy conservation mechanism that can be used to drive thermodynamically unfavorable reactions by coupling endergonic and exergonic reactions through an enzyme complex (42). In the simplest case, this phenomenon is used by anaerobes to increase the yield of ATP through their carbon metabolism by using H₂ as an electron sink for the recycling of NADH to NAD⁺ (42). In the case of the anaerobic gut fungi, the hydrogenosome can be used to generate 2 extra moles of ATP for every mole of glucose that enters glycolysis. However, not all the glycolytic flux can be diverted to the hydrogenosome, because NAD⁺ needs to be regenerated from the NADH that is produced by glycolysis to maintain cellular redox balance. As mentioned before, NAD⁺ is unlikely to be produced by the hydrogen dehydrogenase since the redox potential of NADH/NAD⁺ is too electropositive to reduce H⁺ directly (32). On the other hand, the ferredoxin-based hydrogenase included in the model only recycles the oxidized ferredoxin, produced by PFO, to reduced ferredoxin and does not impact the NADH/NAD⁺ pools in the cell. Sequence homology suggests that the N. lanati hydrogenosome could potentially house a bifurcating hydrogenase, which would couple the reduction of H⁺ to the oxidation of NADH through the ferredoxin produced by PFO. This hydrogenase enzyme complex would allow more flux to be channeled into the hydrogenosome for energy production, since the hydrogenase would now generate NAD⁺ as well as H₂ (Fig. 3).

These findings suggest that there is a significant energetic advantage associated with possessing a bifurcating hydrogenase. The model captures this benefit by predicting a 16% increase in growth rate associated with the use of the bifurcating hydrogenase, as opposed to the ferredoxin hydrogenase (μ = 0.051 h⁻¹ versus μ = 0.044 h⁻¹, respectively). The model also dictates that the production of NAD⁺ shifts from the cytosol to the hydrogenosome when the bifurcating hydrogenase is used, as shown in Fig. 5. However, this requires metabolic flux to be diverted from PFL to PFO in the hydrogenosome. Consequently, significantly more H₂ is predicted to be produced, which is not observed experimentally (Table S3). This discrepancy could be due to metabolic regulation that is unaccounted for in the GEM. Further experimental work needs to be conducted to investigate the potential presence of the bifurcating hydrogenase in the anaerobic gut fungal hydrogenosome.

• Open in new tab
• Download powerpoint

FIG 5

The effects of including additional reactions in the hydrogenosome on NAD⁺ and ATP production show that the putative bifurcating hydrogenase has a large positive effect on NAD⁺ and ATP generation. Conversely, the putative ATP synthase has a negligible effect on both. The increase in growth rate caused by the putative bifurcating hydrogenase is due to NAD⁺ being regenerated in the hydrogenosome; this allows more flux to be channeled into the organelle, which in turn produces more ATP. Flux sampling was used to determine the fluxes associated with NAD⁺ and ATP production in each metabolic configuration. The base case model only includes the ferredoxin hydrogenase. The ferredoxin hydrogenase was replaced by a bifurcating hydrogenase to analyze its effect on the model. Finally, a complex 1, 2, and ATP synthase module was added to the base case model to investigate the consequences of this expanded metabolism. The model was constrained to produce biomass at 90% of the maximum yield; subsequently, 2,000 samples were drawn from each case. The average production of each metabolite in the hydrogenosome and cytosol are shown.

Given that the hydrogenosome in anaerobic fungi is relatively understudied yet related to the mitochondrion (32, 43), we assumed that their metabolite trafficking machinery is similar. Specifically, we assumed that the hydrogenosome has an inner membrane that is impermeable to H⁺, like those of mitochondria (44). This implies that H⁺ can only enter and leave the hydrogenosome through the action of transporters. The H⁺ balance in the hydrogenosome has a direct effect on the ability of the putative ATP synthase to produce ATP. However, the impact of the ATP synthase on ATP production was found to be small, with it only supporting small fluxes (∼3% of the glucose flux into the model) (Fig. 5). This suggests that the putative hydrogenosomal proton gradient (Fig. S4) may not be important for the generation of ATP, as is also suggested by the low expression of the ATP synthase complex subunits (Table 3), or that the proton gradient mechanism is not yet fully understood. Without further experimental evidence, it remains an open question whether the identified ATP synthase components are localized to the hydrogenosome or to some other subcellular organelle (34). Furthermore, the mechanism by which the putative complex 2 operates without the membrane-bound subunits of complex 1 remains to be determined.

Metabolic degeneracy is related to the redox balance.As is typical of unconstrained GEMs, the modeled gut fungal metabolism displays significant degeneracy, as shown by the high degree of flux variability (Table S4). The degeneracy is primarily due to the ability of N. lanati to regulate how NAD⁺ is regenerated through its mixed acid fermentation pathways, i.e., through a combination of lactate dehydrogenase, acetaldehyde dehydrogenase, and alcohol dehydrogenase. Interestingly, the relative mean error between the predicted flux distributions and experimental measurements of the fermentation products is much more sensitive to constraints placed on acetate production than any other single measured external metabolite flux, as shown in Fig. 6. Likewise, constraints placed on lactic acid flux also narrow the deviation of the predicted flux distributions. This effect is due to the different yields of NAD⁺ that can be achieved per mole of pyruvate depending on which mixed-acid fermentation pathway, or combination thereof, is constrained. Without the bifurcating hydrogenase, H₂ production does not significantly impact the overall redox balance of the cell. This possibly explains why its constraint has the smallest effect on the flux variability predicted by the model and may allow the cell to fine tune its metabolism to suit the environmental needs, e.g., sugar availability, by up- or downregulating the flux channeled to the hydrogenosome (2).

• Open in new tab
• Download powerpoint

FIG 6

The absolute relative error between the model predictions and the experimentally measured values suggest that constraining the flux of acetate production has the biggest impact on the model’s accuracy. The flux of acetate (Ac), ethanol (EtOH), formate (For), H₂, and lactate (Lac) was constrained, individually, to their observed ranges (variables on the x axis). The resultant predicted fluxes of these metabolites (generated by sampling 2,000 possible solutions where the biomass objective function was within 90% of its optimal value and subject to the respective additional constraints as shown in the figure) were then compared to the experimental observations as shown in the legend.