Metabolic Analyses of Nitrogen Fixation in the Soybean Microsymbiont Sinorhizobium fredii Using Constraint-Based Modeling

Formulation of metabolic network reconstruction of the symbiotic form of S. fredii CCBAU45436.A metabolic reconstruction (iCC541) was developed to represent the symbiotic form of the microsymbiont S. fredii CCBAU45436. Metabolic reconstructions previously generated (19, 24) were used as a starting point, with 303 reactions selected in total from both models. Reactions used as the template are indicated in Table S1 in the supplemental material. The draft model was manually curated using information from literature, the S. fredii strain pathway database in Kyoto Encyclopedia of Genes and Genomes (KEGG) (28), automated reconstruction databases (29), and other available databases to ensure that it captures the biochemical and physiological knowledge available on nitrogen fixation. A schematic representation of the nutrient and cofactor exchanges between the host plant and the bacteroid is shown in Fig. 1.

• Open in new tab
• Download powerpoint

FIG 1

Schematic representation of main metabolic pathways and transport systems in soybean microsymbiont S. fredii CCBAU45436. Carbon sources in the plant cytosol are transported across the symbiosome and bacteroid membranes in exchange for fixed nitrogen. Malate and succinate enter the TCA cycle to be metabolized. Transport of cofactors required for nitrogen fixation is indicated. Dct, C₄-dicarboxylate transport system; ala-L, alanine; asp-L, aspartate; OAA, oxaloacetate; AKG, 2-oxoglutarate; G3P, glyceraldehyde-3-phosphate; PHB, poly-3-hydroxybutyrate.

A symbiosis reaction was defined to integrate the essential and characteristic elements of the soybean-rhizobia symbiosis association. Specifically, the legume-rhizobium symbiosis determines the accumulation of metabolites inside the microsymbionts related to carbon and energy storage. S. fredii is able to nodulate different legumes and form determinate or indeterminate nodules depending on the host species (30). In particular, poly-3-hydroxybutyrate (PHB), glycogen, and lipids were included in the symbiosis reaction, as S. fredii has been shown to accumulate these compounds in determinate soybean nodules (31–34). The synthesis of storage compounds is important for managing the carbon source and reducing power supplied by the host plant. However, the role of carbon storage compounds in the symbiotic process differs depending on the rhizobium strain (12, 20). Furthermore, bacteroid development and nitrogen fixation are dependent on the amino acids supplied by the host plant (35, 36). Therefore, the uptake and export of various amino acids such as glutamate, alanine, aspartate, and arginine between the bacteroid and host plant were also considered in the symbiosis reaction. Bacteroids isolated from soybean nodules possess an active uptake system for glutamate supplied by the plant cytosol (37, 38). Experimental evidence shows that alanine and aspartate can be supplied to the plant by the bacteroid when malate is provided as a carbon source (39). Meanwhile, arginine also plays a role in the symbiotic process, since the decarboxylation products of arginine are the precursors of polyamines (40, 41). Legume root nodules can contain high levels of polyamines, and it plays a role in the regulation of nutrient exchange (42, 43). Transcriptome data have shown that genes associated with polyamine biosynthesis are upregulated in the symbiosis stage (44). We have also included valine and isoleucine in modeling the symbiotic capacity of S. fredii, because their biosynthetic pathways have been indicated as essential for effective symbiosis in Sinorhizobium and Rhizobium (45), and two transposon Tn5-induced isoleucine and valine auxotrophs of S. fredii HH303 were unable to form effective nodules on soybean as well (46). Moreover, nitrogenous base production has also been demonstrated to be essential for the formation of effective soybean nodules, and S. fredii purine and pyrimidine mutants form ineffective nodules in their hosts (42, 43, 46–49). Therefore, purine and pyrimidine metabolism also constitute a component of the reconstruction. However, symbiosis-defective mutants were not genetically characterized, and further experiments are required to explain this phenotype.

Additionally, symbiotic cofactors are key for the activity of nitrogenase and were therefore added to the symbiosis reaction. Nitrogenase is the enzyme responsible for the reduction of atmospheric nitrogen supplied by the host plant to ammonia. The nitrogenase complex consists of two enzymes: dinitrogenase reductase and dinitrogenase (5, 10). Dinitrogenase reductase is a dimeric iron (Fe) protein and dinitrogenase is a tetrameric iron-molybdenum (FeMo) protein. Nitrogenase activities depend on the availability of energy (ATP), low-potential electrons (iron-sulfur clusters), and microaerobic conditions (5, 50). One of these cofactors is homocitrate, which is a component of the iron-molybdenum cofactor of nitrogenase. From the genome sequence of CCBAU45436, no homocitrate synthase has been annotated; thus, the host plant has to provide homocitrate to the bacteroid to synthesize nitrogenase (51). MucR1 has been identified as a transcription regulator associated with the expression of ion transporters linked to phosphate, molybdenum, zinc, iron, and sulfur. Symbiotic defects were observed in S. fredii mucR1 mutants (52). Additional evidence highlights the role of phosphate in the nitrogen fixation process (33). After nitrogen, phosphate is the most limiting macronutrient, and a high demand for it has been reported in the nodules of nitrogen-fixing legumes (53). Heme and nicotinate are also synthesized by bacteroids, which are supplied to the houseplant for the synthesis of leghemoglobin to maintain the microaerobic condition for nitrogen fixation (36, 46, 47). Other cofactors include pyridoxine (19, 24), magnesium (54), and glutathione (55, 56). The symbiosis equation is a combination of the components described above and it is provided in Table S1. The element balance was checked, and the chemical formulae of the symbiotic product was determined to ensure a molecular weight (MW) of 1 g mmol⁻¹ (57).

Evaluation of iCC541 reconstruction.iCC541 describes 538 reactions encoded by 541 genes as well as 508 metabolites. Table 1 summarizes the main properties of iCC541. iCC541 was evaluated using three methods. First, the reconstruction was validated using experimental data reported in the literature to ensure the model is able to capture unique characteristics of the nitrogen fixation capacity of S. fredii. The conditions tested and captured by the model can be found in Table 2. Published experimental data represent the current knowledge status of an organism, and the validation procedure tests the biological accuracy of the reconstruction. Strain-specific metabolic capabilities included carbon sources and phenotypes of mutant strains.

Additionally, the iCC541 formulation was assessed using a standardized genome-scale metabolic model test suite called memote (26). This metabolic model test approach provides a series of comparisons to benchmark metabolic models and facilitates quality control during model reconstruction. Scores are assigned to properties such as consistency (stoichiometry, mass and charge balance, and metabolite connectivity) and annotation (database cross-references and ontology terms). Scores represent the degree of completeness with respect to the tested property, with a maximum achievable score of 100%. An overall score of 70% was assigned to iCC541 by memote (Fig. 2a). Consistency, metabolites, reactions, genes, and systems biology ontology (sbo) annotations achieved scores of 85%, 68%, 69%, 33%, and 62%, respectively (Fig. 2b). The subsections under consistency, namely, stoichiometric consistency, mass balance, charge balance, and metabolite connectivity, were assigned scores of 100%, 99.2%, 96.9%, and 100%, respectively (Fig. 2c). A maximum score in each of these subsections is desirable to ensure reliable model predictions. These scores are comparable to those for well-curated models of Escherichia coli (E. coli core model and iAF1260b) and Saccharomyces cerevisiae (iMM904) from the Biochemical, Genetic, and Genomic (BiGG) knowledge base (58). iCC541 has significantly improved quality in terms of model consistency and annotations compared to that of the published Sinorhizobium meliloti reconstructions, iGD1575 and iHZ565 (25, 59), which were also tested by memote in this study using the published Systems Biology Markup Language (SBML) files (Fig. 2). The full report returned by memote for iCC541 and the two S. meliloti models can be found in Fig. S1.

• Open in new tab
• Download powerpoint

FIG 2

Scores and consistency assigned to iCC541, iHZ565, and iGD1575 by memote. (a) Total scores. (b) Scores by main categories: consistency and annotations (systems biology ontology [sbo], reactions, metabolites, and genes). (c) Consistency across subsections.

Finally, the simulation performance of iCC541 was compared against the S. meliloti models. Inconsistencies reported by memote were solved in order to facilitate cross-comparisons between models. For that purpose, mass and charge imbalances and energy-generating cycles (ATP) (60) were fixed for the S. meliloti models. The symbiosis reactions were scaled to make sure the total mass of the material through the symbiosis reactions was the same in all three models in order to compare symbiosis yields (symbiosis product in grams per millimole carbon) (Table 3). The symbiosis reaction flux was maximized under microaerobic conditions for each model: 0.0107, 0.0145, and 0.0130 g mmol⁻¹ were obtained for iCC541, iHZ565, and iGD1575, respectively. Similar yields relative to nitrogen uptake were observed for iCC541 and iGD1575: 0.0747 and 0.0906 g mmol⁻¹ nitrogen, respectively, whereas iHZ565 predicts a higher yield with 0.1774 g mmol⁻¹. C/N ratios obtained under standard nitrogen fixation conditions (defined according to literature) were comparable in iCC541 and iGD1575, with iHZ565 showing a lower ratio. The number of share reactions, metabolites, and genes between iCC541 and the S. meliloti models are indicated in Table 4.

Symbiotic pathway utilization analyses.In symbiotic nitrogen fixation nodules, the host plant provides C₄-dicarboxylic acids (malate or succinate) to bacteroids as the major energy source, while the microsymbiont returns organic fixation products to the host. The primary carbon source, malate, is transported into bacteroids to be metabolized through the tricarboxylic acid (TCA) cycle to provide reducing equivalents, energy, and precursors for biosynthetic pathways (12, 61, 62). The TCA cycle together with oxidative phosphorylation generates the energy and conditions required for the nitrogen fixation process. Oxygen levels must be maintained at the minimum to avoid inhibition of the nitrogenase activity, since it acts as a regulator of nitrogen fixation genes (nif and fix genes) (12, 20, 63). To examine the extent to which the model can capture these metabolic characteristics, flux balance analysis (FBA) was performed by maximizing the symbiosis reaction with malate as the primary carbon source and succinate and inositol as secondary carbon sources, with a small amount of oxygen available (maximum uptake rates being constrained by experimental data) (see Table S1). Figure 3 describes the flux distribution predicted by the model. In detail, the simulated exchange profile of the model shows the active uptake of glutamate and cofactors (including homocitrate, molybdenum, thiamine, cobalamin, zinc, iron, myo-inositol, phosphate, magnesium, and sulfate) as expected from the formulation of the symbiosis reaction (62, 64, 65). Oxygen is utilized by cytochrome oxidase to generate ATP. Ammonia, alanine, and aspartate are transported to the plant cell cytosol (39, 87). In addition, the predicted internal flux distribution describes the metabolic route of malate. Malate is actively transported into the cell and converted to pyruvate and CO₂ via malate dehydrogenase. Then, pyruvate is decarboxylated by pyruvate dehydrogenase to form acetyl coenzyme A (acetyl-CoA) and enters the TCA cycle. The full oxidative TCA cycle is utilized by the bacteroid during the process, according to the model predictions (nonzero fluxes confirmed by flux variability analysis). Enzymes participating in the TCA cycle appear to have an active role in the nitrogen fixation process in S. fredii. TCA cycle mutants of S. meliloti, Rhizobium tropici, and Rhizobium leguminosarum have the common phenotype of being unable to fix nitrogen (68–70). The model, predicting the use of malate to produce energy and products for symbiosis through the TCA cycle, is therefore able to capture this important aspect of rhizobium-legume symbiosis. However, the role of the TCA cycle varies among rhizobia, with high metabolic plasticity observed in Bradyrhizobium japonicum (71, 72). PHB biosynthesis and gluconeogenesis pathways are active according to the simulation results and the definition of the symbiosis reaction. The use of these pathways is validated by the fact that both of these carbon storage compounds have been found to be present in determinate nodules of soybean (31–34). The model also predicts a low but nonzero flux through some reactions of the pentose phosphate pathway. This pathway generates ribose 5-phosphate, which is a precursor of purines and pyrimidines. Reduced symbiotic properties have been observed in pyrimidine and purine mutants of S. fredii (42, 43, 46–49). According to our reconstruction formulation and in silico results, the ammonia exported to the plant cell cytosol originates from the ammonia produced by the nitrogenase reaction (Fig. 3).

• Open in new tab
• Download powerpoint

FIG 3

Flux distribution during symbiotic nitrogen fixation in S. fredii. Reactions are colored according to a scale of purple (flux of 0.04 mmol g⁻¹ [dry weight] h⁻¹) to red (flux of 1.8 mmol g⁻¹ [dry weight] h⁻¹ or greater). Reactions carrying zero flux or less than 0.04 mmol g⁻¹ (dry weight) h⁻¹ are colored gray. Colored arrowheads indicate the directionality of the flux. Flux metabolic networks were constructed using the Web tool Escher.

Identification of genes involved in symbiosis.Next, the model was used to identify symbiotic genes under standard nitrogen fixation conditions. As defined by the gene-protein-reaction (GPR) association, a gene was classified as a symbiotic gene if one or more reactions associated with that gene need to carry a nonzero flux to satisfy the symbiotic objective function. Results of simulations were categorized according to the effect of the gene deletion compared to that of the wild-type strain (see Table S2). According to the simulations, 141 of 541 genes (26%) were predicted to affect the symbiotic process, wherein 121 genes were classified as essential and 20 genes as having a partial effect, where their deletions resulted in a reduction in the capacity of S. fredii to fix atmospheric nitrogen. Most of these genes (84%) were located in the chromosome, while the others were present in the symbiosis plasmid, accessory plasmids, and chromid (Table S2). Coexpression patterns in the replicon of S. fredii strains have been observed under both free-living and symbiotic conditions (17).

Many genetic studies have been conducted to examine the legume-rhizobium symbiosis, providing important insights into the symbiotic capacity of rhizobia. A detailed comparative genomic study of rhizobia that nodulate soybean was used to validate the predicted symbiotic genes from our model (16). Genes predicted to affect the symbiotic process were compared to a list of known symbiosis genes. The analysis matched 50 of the 141 predicted symbiotic genes to previously identified genes in the literature. Based on our model, mutations in these genes had a negative impact on the symbiotic phenotypes of rhizobial strains. For example, symbiotic genes related to oxidative phosphorylation were identified. According to the simulation in our model, oxidative phosphorylation must be active in vivo to effectively fix nitrogen, as previously identified by pathway analysis (63). A complete list of these genes and the references providing the corresponding validating evidence are available in Table S3. Among the 91 remaining predicted symbiotic genes, 58 are orthologs to essential genes identified by iHZ565 (24). These genes have already been listed in the Database of Essential Genes (DEG) (73) and therefore have been shown to perform core functions that are required under both free-living and symbiotic conditions. The remaining 33 genes have not been shown to be essential to the symbiotic process. Their locations and predicted functions are provided in Table S2. Most of these genes are linked to the purine, pyrimidine, TCA, and amino acid pathways. One example is the gene AB395_00003474 (sucB). sucB and sucA encode subunits of the α-ketoglutarate dehydrogenase which catalyzes the conversion of α-ketoglutarate to succinyl-CoA in the TCA cycle. R. leguminosarum mutants form ineffective nodules (69). However, B. japonicum sucA mutants only showed a delay in the development of bacteroids, but their nitrogen fixation rates were similar to those of the wild type (72). On the other hand, AB395_00003483, which encodes a succinate dehydrogenase flavoprotein subunit, was predicted here to be a symbiotic gene. A gene encoding succinate dehydrogenase A (sdhA) has been identified as essential in Bradyrhizobium sp. strain ORS278 (74). Another symbiotic gene predicted by our model, AB395_00004794, is involved in the isoleucine and valine biosynthesis pathway. Tn5 experiments have successfully identified this pathway as being related to the symbiotic process in S. fredii HH303, but the genes have not been characterized (46). The predicted gene is an ilv gene that encodes a dihydroxy-acid dehydratase. A BLASTP search of the amino acid sequence of the dihydroxy-acid dehydratases from S. meliloti 1021 (NP_384214, NP_436655, and NP_386934) against the S. fredii genome revealed identities of 89.37%, 64%, and 31.67%, respectively, with AB395_00004794 at the amino acid level. Mutations of ilv genes have been reported to result in phenotypes ranging from no nodule formation (Nod⁻) (35) to being unable to fix nitrogen (Nod delayed, Fix⁻) (75). In addition, genes related to the transport of cofactors such as sulfate were also predicted as novel symbiotic genes by our model. Additional experiments are required to confirm the functions and relevance of these genes in the TCA cycle, branched-amino-acid biosynthesis, and transport systems to the symbiotic processes in S. fredii. As stated above, myo-inositol uptake was active according to the model predictions. Simulations showed a slight reduction in the fixation rate, and shadow price analyses indicated it was a rate-limiting substrate. This is qualitatively but not quantitatively consistent with a previous finding that a mutation in the myo-inositol dehydrogenase gene (idhA) caused a much more significant reduction in the nitrogen fixation rate (66). Then, false-negative cases were evaluated. These cases represent discrepancies with experimentally observed phenotypes and indicate potential errors or knowledge gaps in the model. The symbiotic genes predictions exhibited 45 (8%) false negatives. The list of these genes can be found in Table S4. Analysis of these genes revealed active roles in processes such as nodulation, bacteroid development, microoxic respiration regulation, osmoregulation, and the inositol pathway. Regulation processes and steps previous to the symbiotic stage (e.g., root nodule formation) were not included in the current scope of the model, which explains this result. Additional false-negative cases occurred with knockout of genes involved in inositol catabolism, including iolA (AB395_0000388), iolC (AB395_00002655), and iolD (AB395_00002654) (76, 77). Shadow price analyses identified inositol as a rate-limiting step as stated above. pstC (AB395_0000141) and pit (AB395_00003896) involved in phosphate transport were classified as nonessential (78). FBA may use any transport reaction in the network to fix nitrogen under the standard nitrogen fixation condition. Additional constraints are required to simulate phosphate transport via high- and low-affinity systems.

Nitrogen fixation capacity of S. fredii bacteroids in wild versus in cultivated soybean.Finally, transcriptome profiles were used to determine the nitrogen fixation capacity of CCBAU45436 with two different soybean cultivars, the wild soybean accession W05 and the cultivated soybean accession C08. Condition-specific models were built from iCC541 using the E-Flux approach (79). Gene expressions were used to constrain reaction flux bounds to assess the impact of the host plant on nitrogen fixation. Tables S5 and S6 list the flux ranges for each condition and raw expression data, respectively. Flux distributions are shown in Fig. 4. Exchange profiles of the models showed carbon and essential nutrient requirements. Additionally, certain secretion patterns overlap one another. As expected, malate is transported inside the cell and pyruvate is generated by the reaction catalyzed by malic enzyme. Pyruvate is further decarboxylated by pyruvate dehydrogenase to form acetyl-CoA. In both soybean cultivars, soybean bacteroid diverts acetyl-CoA into the TCA cycle and into the production of the lipid-like polymer PHB, both of which are pathways known for electron allocation in nitrogen-fixing bacteroids (80). Although similitudes are observed in the flux distributions, there are also differences in the use of some reactions and/or their magnitudes (Fig. 4c). To fully assess the impact on nitrogen fixation capacity, yields were calculated based on the model predictions (Table 5); 0.0032 and 0.0097 g symbiosis product per mmol of carbon uptake were obtained in the symbiosis with W05 and C08, respectively. Additionally, there was a difference in the flux of fixed ammonia with 0.0038 and 0.0275 mmol g⁻¹(dry weight) h⁻¹ transported from the bacteroid to W05 and C08, respectively. Simulations predicted a higher nitrogen fixation capacity of CCBAU45436 with C08 than with W05. C/N ratios of 0.1435 and 0.0426 were calculated for the C08 and W05 symbiotic conditions, respectively. The ammonia produced by nitrogen fixation in bacteroids is transported back to the host plant and assimilated into glutamate and glutamine by glutamine synthetase and glutamate synthase (20). In determinate nodules, such as those in soybean roots, glutamine enters the purine biosynthetic pathway to be further converted to ureides, allantoin, and allantoic acid (65, 81). Ureides, products of nitrogen fixation, are then transported from the root to the shoot. Thus, a higher ureides and nitrogen content should be expected in the root nodules of C08 compared to those of W05 when soybean is inoculated with CCBAU45436. Ureides and nitrogen content in root nodules are commonly used as an indicator of nitrogen fixation capacity and were previously investigated by our group using the same symbiotic conditions, with W05 and C08 being inoculated with S. fredii CCBAU45436. The results corroborate this model prediction (21). Experimental observations showed cultivated soybean outperformed wild soybeans on all the tested nitrogen fixation-related traits, in particular, the total nitrogen and total ureides accumulation (21). An imbalance in carbon and nitrogen metabolism in wild soybean was suggested as the cause of the lower capacity for nitrogen fixation. This experimental evidence validates the model predictions and indicates a consequence of the process of domestication and human selection in soybean. Furthermore, a simple evolutionary sensitivity analysis was performed to identify possible reactions (Fig. 4c) responsible for the higher performance of S. fredii with the cultivated soybean, by starting from the W05 condition (i.e., the model with flux bounds inferred from the E-Flux method used in the previous analysis), iteratively imposing bounds for the C08 condition randomly and fixing the bounds that increase nitrogen fixation until no further increase could be achieved (Fig. 5a). Across 10,000 simulations, different activity levels in the inositol pathway, cytochrome oxidase reaction, and iron transport were consistently identified to be the determinants for the predicted phenotypes (Fig. 5b). The inositol pathway is a lower-yielding pathway having a marginal effect on improving nitrogen fixation under the W05 condition (increased by 20%), while the cytochrome oxidase reaction appears to be a major bottleneck. However, a gene encoding an ATP-binding cassette (ABC)-type metal ion transport system (AB395_00003450) seemed to limit the nitrogen fixation rate under the C08 condition because of its lower expression level under the C08 condition than the W05 condition. Potential explanations for this include the existence of regulatory or resource allocation constraints that are not accounted for by the model. The proportionality assumption between maximum enzymatic activities and transcript levels applied to the model can serve as a proxy for estimating relative activities across conditions in the absence of additional data, but it is probably an oversimplification and may underestimate the activity of the metal ion transport reaction under the C08 condition or vice versa. Further in vitro and in vivo testing of these potential targets will help resolve this. Overall, a fine tuning of the energy, reducing factors, and regulations of gene expression are required to optimize the nitrogen fixation process.

• Open in new tab
• Download powerpoint

• Open in new tab
• Download powerpoint

FIG 4

Nitrogen fixation capacity of S. fredii CCBAU45436 in wild soybean accession W05 (a) and cultivated soybean accession C08 (b). (c) The difference between the flux distributions of C08 and W05. Flux metabolic networks were constructed using the Web tool Escher. Reactions are colored according to a scale of purple (median flux or greater) to red (maximum flux). Colored arrowheads indicate the directionality of the flux. Reactions carrying zero flux or less than the median are colored gray.

• Open in new tab
• Download powerpoint

FIG 5

Evolutionary sensitivity analysis. (a) Symbiosis fluxes were determined by simulations starting from the W05 condition where bounds of the C08 condition were randomly imposed. (b) Reactions were identified where S. fredii inoculated in cultivated soybean outperformed that in wild soybean in nitrogen fixation.