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

Abstract

Anaerobic gut fungi in the phylum Neocallimastigomycota typically inhabit the digestive tracts of large mammalian herbivores, where they play an integral role in the decomposition of raw lignocellulose into its constitutive sugar monomers. However, quantitative tools to study their physiology are lacking, partially due to their complex and unresolved metabolism that includes the largely uncharacterized fungal hydrogenosome. Modern omics approaches combined with metabolic modeling can be used to establish an understanding of gut fungal metabolism and develop targeted engineering strategies to harness their degradation capabilities for lignocellulosic bioprocessing. Here, we introduce a high-quality genome of the anaerobic fungus Neocallimastix lanati from which we constructed the first genome-scale metabolic model of an anaerobic fungus. Relative to its size (200 Mbp, sequenced at 62× depth), it is the least fragmented publicly available gut fungal genome to date. Of the 1,788 lignocellulolytic enzymes annotated in the genome, 585 are associated with the fungal cellulosome, underscoring the powerful lignocellulolytic potential of N. lanati The genome-scale metabolic model captures the primary metabolism of N. lanati and accurately predicts experimentally validated substrate utilization requirements. Additionally, metabolic flux predictions are verified by ¹³C metabolic flux analysis, demonstrating that the model faithfully describes the underlying fungal metabolism. Furthermore, the model clarifies key aspects of the hydrogenosomal metabolism and can be used as a platform to quantitatively study these biotechnologically important yet poorly understood early-branching fungi.IMPORTANCE Recent genomic analyses have revealed that anaerobic gut fungi possess both the largest number and highest diversity of lignocellulolytic enzymes of all sequenced fungi, explaining their ability to decompose lignocellulosic substrates, e.g., agricultural waste, into fermentable sugars. Despite their potential, the development of engineering methods for these organisms has been slow due to their complex life cycle, understudied metabolism, and challenging anaerobic culture requirements. Currently, there is no framework that can be used to combine multi-omic data sets to understand their physiology. Here, we introduce a high-quality PacBio-sequenced genome of the anaerobic gut fungus Neocallimastix lanati Beyond identifying a trove of lignocellulolytic enzymes, we use this genome to construct the first genome-scale metabolic model of an anaerobic gut fungus. The model is experimentally validated and sheds light on unresolved metabolic features common to gut fungi. Model-guided analysis will pave the way for deepening our understanding of anaerobic gut fungi and provides a systematic framework to guide strain engineering efforts of these organisms for biotechnological use.

Keywords: 13C metabolic flux analysis; Neocallimastigomycota; Neocallimastix lanati; anaerobes; anaerobic fungi; flux balance analysis; genome-scale metabolic model; nonmodel fungus.

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.

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.

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.

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

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.

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.