Metabolic control of nitrogen fixation in rhizobium-legume symbioses

Catabolism of dicarboxylates at low oxygen drives N2 fixation and promotes ammonia secretion in rhizobium-legume symbioses.


INTRODUCTION
Biological nitrogen fixation provides 50 to 70 Tg of bioavailable nitrogen in agricultural systems per year (1) and sustains global food security. The most efficient contribution to biologically fixed nitrogen is from symbioses between legumes and rhizobia (2), which are soil bacteria that induce formation of nodules on plant roots. Inside nodules, rhizobia differentiate into bacteroids that reduce atmospheric N 2 into ammonia for secretion to the plant host in exchange for dicarboxylates, primarily succinate and malate (3,4). Succinate and malate are typically metabolized via malic enzyme and pyruvate dehydrogenase, yielding acetyl-coenzyme A that can be oxidized in the tricarboxylic acid (TCA) cycle (5). Whether bacteroids need a complete TCA cycle remains unclear, as it is essential for Rhizobium and Sinorhizobium species, but mutants of Bradyrhizobium japonicum lacking 2-oxoglutarate dehydrogenase activity achieve wild-type levels of nitrogen fixation on a per-bacteroid basis (6,7). Flux through the TCA cycle produces reduced nicotinamide adenine dinucleotide (NADH) and flavin adenine dinucleotide (FADH 2 ), which provide electrons for nitrogen fixation and adenosine triphosphate (ATP) generation (8). Despite rhizobia requiring oxygen for ATP synthesis, oxygen levels in nodules are only 10 to 40 nM (9), which is a requirement for highly oxygen-sensitive nitrogenase enzyme activity (10). The challenge of balancing carbon allocation under these conditions is evidenced by the synthesis of lipids and carbon polymers, such as polyhydroxybutyrate (PHB), which indicate imbalances in nutrient supply (8,11). PHB and lipids have been suggested to play a role as carbon and redox sinks for bacteroids, although their accumulation is variable between rhizobial strains and their role remains to be elucidated (8,12,13).
The defining distinction between nitrogen fixation by rhizobial bacteroids compared to free-living bacteria is the secretion of fixed ammonia to the plant. However, there is no known metabolic mechanism forcing secretion of fixed nitrogen to the plant instead of assimilation by the bacteroid. In addition to the main secretion product ammonia, a substantial portion of fixed nitrogen is apparently secreted in the form of alanine and aspartate (14,15). At the same time, ammonia assimilation by the glutamine synthetaseglutamine oxoglutarate aminotransferase (GS-GOGAT) system, which is active in free-living rhizobia, is down-regulated during the symbiosis (16)(17)(18).
Because of the complexity of bacteroid metabolism, several studies have used computational approaches to investigate the symbiosis. Notably, metabolic models for various rhizobial species have been reconstructed (19)(20)(21)(22), and some reconstructions have been integrated with models for the host plant (23,24). Metabolic models describe all enzymatic and transport reactions encoded in an organism's genome and enable the simulation of metabolic flux distributions under defined environmental conditions (25,26). The most common approach for analyzing metabolic models is flux balance analysis (FBA), where an objective function reflecting the evolutionary goal of the organism under investigation is optimized to calculate a flux distribution (27). Maximum cellular growth is the most commonly used objective function and has been found to reproduce experimental results for strains adapted to growth under laboratory conditions (28), but this is clearly not applicable in the case of growth-arrested bacteroids. In silico studies of bacteroid metabolism, which have so far mostly applied standard FBA methods, defined an objective reaction comprising nitrogen export to the plant and synthesis of storage compounds (19)(20)(21)(22).
The result of FBA calculations is a single flux distribution, which is often not a unique solution for the optimization problem (27). In contrast, methods such as elementary flux mode enumeration describe all minimal sets of steady-state fluxes through the metabolic network (29). While this approach is more computationally expensive than FBA and currently not feasible for genome-scale metabolic networks (30), it provides a more comprehensive and unbiased description of metabolism that does not rely on an artificially defined objective. Addressing the problem of combinatorial explosion during elementary flux mode enumeration, elementary conversion modes (ECMs) have been proposed as an alternative approach to compute all possible stoichiometries between input and output metabolites (31). While ECMs do not provide information on the underlying metabolic pathways for a specific conversion, they can still capture all metabolic capabilities of an organism. The feasibility of ECM enumeration for genome-scale metabolic networks when focusing on subsets of metabolites has recently been demonstrated (32).
Despite decades of research efforts, a comprehensive view of the links between central carbon and nitrogen metabolism in bacteroids is missing. Most experimental studies have focused on individual metabolic pathways, and computational models have used artificial objective functions that may not be relevant in the natural system.
In this study, we combine experimental and metabolic modeling approaches to explain fundamental features of bacteroid metabolism. Using proteome, transcriptome, and gene essentiality data, we reconstruct a model of metabolic pathways active during nitrogen fixation in Rhizobium leguminosarum bv. viciae. We implement modeling strategies that circumvent the limitations of traditional FBA to explain the importance of the experimentally observed storage polymer synthesis and amino acid secretion in bacteroids. We further investigate the role of the TCA cycle during nitrogen fixation and validate model predictions by 13 C metabolic flux analysis of Azorhizobium caulinodans. Our model provides insights into the fundamental constraints on rhizobial metabolism during symbiotic nitrogen fixation. An improved understanding of metabolic processes in bacteroids is of central importance for ongoing efforts in optimizing existing symbioses and engineering novel plantmicrobe interactions for sustainable agriculture (2,33).

Data-based reconstruction of a bacteroid metabolic model
We reconstructed a metabolic model for pea bacteroids of R. leguminosarum bv. viciae 3841 using bacteroid-specific experimental data. First, we quantified the proteome of unlabeled bacteroids relative to 15 N-labeled free-living bacteria (Fig. 1, figs. S1 and S2, Supplementary Text, and data S1). In addition, genes up-regulated in transcriptional datasets of bacteroids ( fig. S3) (34) and genes identified as specifically essential for symbiosis formation by insertion sequencing (INSeq) (35) were included in the model. The final core metabolic network named iCS323 contained 323 genes, 237 metabolites, and 299 reactions, 207 of which are metabolic (excluding transport, demand, and sink reactions) (table S1 and data S2 to S4). Of the 207 metabolic reactions, 177 (86%) are supported by experimental evidence from at least one of the bacteroid-specific datasets. iCS323 was evaluated using MEMOTE (36), confirming the absence of stoichiometrically balanced cycles, orphan metabolites, and deadend reactions (data S5).
The main pathways in the final model were central carbon metabolism (TCA cycle, gluconeogenesis, and pentose phosphate pathway), amino acid metabolism, and carbon polymer synthesis. As an initial validation, malate, succinate, and oxygen uptake were constrained according to the flux boundaries defined in a modeling study of Sinorhizobium meliloti (21), and standard FBA maximizing nitrogenase activity was performed. The obtained flux distribution captured key features of bacteroid metabolism, including use of the TCA cycle, pyruvate synthesis via malic enzyme, and minor activity of gluconeogenesis ( fig. S4) (5,8,37).
When comparing gene essentiality predictions to mutant phenotypes identified by INSeq (35), we found that 280 (87%) of the genes in the model had an in silico phenotype that agreed with the experimental evidence (data S6). Genes were defined as essential in silico if their deletion prevented flux through the nitrogenase reaction in the model. Of the 31 false negatives (in silico nonessential genes that were essential according to INSeq), 9 were involved in sugar metabolism. In agreement with the results of a recent modeling study of S. meliloti (24), this may suggest that rhizobia differentiating into bacteroids have access to sugars as a carbon source, whereas nitrogen-fixing bacteroids do not. In addition, four genes involved in serine metabolism were incorrectly predicted to be  nonessential. A possible explanation is the role of serine as a precursor for cysteine biosynthesis. Cysteine plays a role in the synthesis of iron-sulfur clusters for the nitrogenase enzyme (38), which is not explicitly accounted for in the model. Pyruvate kinase (PykA) was further predicted to be nonessential in disagreement with the INSeq data. A pykA mutant had higher nitrogenase activity than the wild type in plants harvested after 28 days (5), indicating that the INSeq phenotype may be a result of developmental delay and the model correctly identifies the gene as nonessential. Overall, our model shows good predictive quality for gene essentiality. It is important to note that the perfect agreement of in silico predictions with the experimental data is not expected. This is due to our model being limited to central metabolic pathways in bacteroids, neglecting, for example, the synthesis of nucleotides for DNA replication. In addition, genes that are determined to be essential for bacteroids by INSeq may actually cause a growth defect at earlier stages of symbiosis formation and are not necessarily required for nitrogen fixation itself.

Characterization of bacteroid metabolism using ECMs
Previously published modeling studies of bacteroid metabolism used FBA and maximized a lumped objective reaction comprising ammonia and amino acid export as well as storage polymer synthesis (19)(20)(21)(22). While this approach constrains flux distributions to reflect experimentally observed phenotypes, it also introduces an artificial stoichiometric coupling between the objective metabolites, which precludes the investigation of changes in carbon and nitrogen allocation depending on nutrient availability. To characterize bacteroid metabolism with a minimum number of preset assumptions, we used ecmtool (32) to enumerate ECMs of the metabolic network ( Fig. 2 and data S7). In the context of ECMs, the overall transformation of nutrients into secreted metabolites and cell components is called a conversion. The set of possible conversions is determined by an organism's metabolism, which is described by a metabolic model. Enumerating an exhaustive set of minimal conversions (ECMs) thus provides a complete overview of all metabolic capabilities in terms of input-output stoichiometry. This approach gives a comprehensive description of metabolism without assuming optimality with respect to a specific objective.
Conversion modes were found for all expected products (ammonia, alanine, aspartate, and storage compounds) with a minimal set of inputs, mainly malate/succinate, oxygen, and N 2 (data S8). The presence of conversion modes for storage polymer synthesis without nitrogenase activity suggests that the two processes are not inherently linked, which agrees with reports of storage polymer accumulation before the onset of nitrogen fixation (12,39). Potential co-catabolism of an amino acid in addition to malate/succinate was investigated by analyzing conversion modes with each of the proteinogenic amino acids or -aminobutyric acid (GABA) as an additional input. Conversion modes were found for arginine, cysteine, GABA, glutamine, glutamate, glycine, serine, and threonine. All conversions involving arginine, cysteine, or serine had a negative net nitrogen output and were therefore considered unlikely to be biologically relevant. For the remaining conversion modes with a positive net nitrogen output, no benefit was found for any amino acid in terms of oxygen demand or carbon cost ( fig. S5). We therefore limited our analysis to conversion modes using GABA, which is known to be metabolized in pea bacteroids (40).
Oxygen demand per carbon uptake was decreased for all conversion modes that produced storage polymers compared to conversions that did not ( Fig. 3A and fig. S6). Carbon storage polymers thus function as carbon and redox sinks under oxygen-limiting conditions, enabling electrons to be partitioned to polymer synthesis rather than oxygen as a terminal electron acceptor. The oxygen uptake relative to fixed N 2 was increased for conversion modes generating glycogen or lipids (Fig. 3B), consistent with the ATP requirement for producing these storage compounds, which would add to the ATP demand of the nitrogenase reaction. In addition, conversion modes generating carbon polymers increased the carbon cost per nitrogen supplied to the plant (Fig. 3C). This aligns with observations that plant nodule cells accumulate starch when occupied by bacteroid glycogen synthase mutants (12); i.e., the plant has excess carbon when bacteroid polymer synthesis is restricted.
The carbon cost determined for conversion modes without polymer production is similar to the theoretical cost of nitrogen fixation (2.5 g of carbon per gram of nitrogen), whereas that for PHB-and some palmitate-and glycerolipid-producing conversion modes is close to the experimentally observed value (8 g of carbon per gram of nitrogen) (2). Given that the plant host regulates the supply of nutrients such as carbon and oxygen in response to nitrogen output (41,42), this would potentially limit excess diversion of carbon into polymers.
Conversions without carbon polymer production could be subgrouped into those generating alanine and/or aspartate and those only producing ammonia, where amino acid secretion reduced oxygen demand per carbon uptake (Fig. 3A) and oxygen demand per nitrogen output ( fig. S6). Alanine dehydrogenase catalyzes the NADH-dependent synthesis of alanine from pyruvate and ammonia, making this pathway an oxygen-independent carbon sink for NAD + regeneration as well as a sink for carbon and protons. Synthesis of alanine in particular enables bacteroids to maintain nitrogen export to the plant in a low-oxygen environment and explains the mixed secretion of ammonia and amino acids observed experimentally (14). Amino acid secretion could be favored over polymer synthesis under certain conditions because it allows for removal of carbon from the bacteroid rather than intracellular accumulation. The same trends for storage polymer synthesis and amino acid secretion were observed for ECMs of a model for Sinorhizobium fredii bacteroids ( fig. S7) (19), indicating that these principles are likely to govern symbiotic metabolism across rhizobial species.

Role of oxygen limitation in shaping bacteroid metabolism
We next characterized the network response to different carbon and oxygen availability when optimum nitrogenase activity is maintained. Because uptake fluxes of bacteroids are difficult to determine experimentally, we performed phenotype phase plane analysis (43). Avoiding artificial biases of carbon and nitrogen allocation, nitrogenase activity was evaluated without maximizing the production of storage compounds. Four feasible regions with distinct metabolic behavior were identified, with phase I characterized by carbon limitation and increasing oxygen limitation from phase II to most limited in phase IV ( Fig. 4A and fig. S8).
Flux variability analysis showed that with increasing oxygen limitation, flux through TCA cycle enzymes, as represented by citrate synthase, decreased (Fig. 4B). Carbon was increasingly channeled into pyruvate, which caused accumulation of PHB and lipids as well as alanine production, with all nitrogen being secreted in the form of alanine instead of ammonia in phase IV. The tendency to increase alanine secretion at low oxygen levels has recently been shown experimentally in B. japonicum (44). Alanine secretion is thus important to sustain bacteroid metabolism under oxygen limitation, which probably occurs in the natural system (45). This is consistent with a 20% reduction in dry weight of peas inoculated with alanine dehydrogenase mutants of R. leguminosarum that no longer secrete alanine (46).
Independently of maximum nitrogenase activity, general network properties were investigated using ensemble-evolutionary FBA (47). The results supported increased PHB and lipid synthesis under oxygen-limiting conditions and decreased activity of the TCA cycle (Fig. 4, C and D). The trend for increased alanine synthesis was not observed, indicating that it is linked to nitrogen fixation rather than being an inherent network property. The identified shifts in metabolic behavior support the role of alanine, PHB, and lipids as sinks for carbon when oxygen is limiting, and usage of the TCA cycle becomes disadvantageous owing to the accumulation of reduced electron carriers (8), which agrees with the conversion mode analysis.

Metabolic constraints on ammonia assimilation by bacteroids
The predicted down-regulation of the TCA cycle would reduce the availability of 2-oxoglutarate required by the GS-GOGAT pathway, which is the sole pathway for ammonia assimilation coupled to growth in rhizobia. It should be noted that rhizobial GS-GOGAT mutants cannot grow on ammonia as a nitrogen source, and alanine dehydrogenase, for example, cannot substitute for GS-GOGAT (18). Increased ammonia assimilation into glutamate by bacteroids would require increased TCA cycle activity and hence increase oxygen demand as indicated by a significant positive correlation between oxygen uptake and GS activity in the conversion mode analysis (Fig. 5A). To assess the impact of ammonia assimilation on the metabolic fluxes in bacteroids, we forced flux through a demand reaction for glutamate in the model. In this scenario, flux through the decarboxylating arm of the TCA cycle had to be maintained under oxygen limitation to supply 2-oxoglutarate for ammonia assimilation, which caused an increased oxygen demand (Fig. 5B).
Our modeling results indicate the importance of the carbon source provided to bacteroids to enforce use of the TCA cycle as the main catabolic pathway. It is currently unclear why plants provide bacteroids with C 4 -dicarboxylates, especially considering that photosynthate is transported to nodules as sucrose (48). We therefore compared the effects of using sucrose instead of malate as a carbon source in silico. For a given carbon uptake rate, the model predicted higher nitrogenase activity for sucrose compared to malate (Fig. 5C) in agreement with an integrated plant-bacteroid model for S. meliloti (24). Furthermore, less oxygen was needed for sucrose ( Fig. 5B) (or glucose; fig. S9) catabolism, which is in accordance with experimentally determined oxygen uptake rates of free-living R. leguminosarum (Fig. 6A). Catabolism of arabinose, a sugar metabolized via 2-oxoglutarate, induced an oxygen demand similar to growth on succinate. This suggests that catabolism of TCA cycle intermediates generally creates a high oxygen demand and thus causes a more severe growth impairment compared to glucose under low oxygen conditions (Fig. 6C). However, the NADH/NAD + ratio was similar for growth on glucose and arabinose but significantly higher for succinate catabolism (Fig. 6B). This could be explained not only by 2-oxoglutarate supply from arabinose catabolism partly obviating use of the decarboxylating TCA cycle arm but also by differences in the regulation of sugar versus dicarboxylate uptake and metabolism.
Overall, these findings indicate that supply of dicarboxylates such as succinate to bacteroids creates both a high oxygen demand and a highly reduced redox ratio, which is important for supply of electrons to the nitrogenase enzyme (8). The absolute value of the redox ratio did not change significantly for different oxygen concentrations and was also not correlated with the growth rate (Fig. 6C). A slightly higher NADH/NAD + ratio was measured for cultures grown on succinate at 21% oxygen compared to 1% oxygen. However, this is most likely due to minor effects of sample preparation, which is very time-sensitive in the case of fast-growing cultures with high respiration rates. Thus, the redox state of R. leguminosarum is mainly determined by the carbon source, implying that the nature and quantity of the supplied carbon source are of central importance for defining metabolic fluxes in bacteroids, with succinate providing a higher redox value than glucose. Plants provide dicarboxylates as a carbon source to bacteroids although they are less efficient at supporting nitrogen fixation and increase oxygen demand relative to sucrose in the oxygen-limited nodule.

Metabolic flux analysis of A. caulinodans supports model predictions
To validate core findings of our model predictions experimentally, we performed 13   We chose this rhizobial strain for its ability to perform nitrogen fixation in both free-living conditions and in symbiosis with Sesbania rostrata. When comparing non-diazotrophic growth of A. caulinodans at different oxygen levels, PHB synthesis increased threefold under microaerobic (3% oxygen) conditions compared to aerobic growth and further increased in bacteroids ( Fig. 7 and data S9). This agrees with the predicted importance of storage polymer synthesis for balancing carbon allocation under oxygen-limited conditions. TCA cycle fluxes increased during free-living diazotrophic growth compared to non-diazotrophic growth, which is consistent with an increased demand for reductant for nitrogen fixation as well as for cell growth. However, an overall down-regulation of the TCA cycle was observed when comparing bacteroid metabolism to non-diazotrophic growth under microaerobic conditions. The highest relative downregulation was observed for the decarboxylating arm of the TCA cycle between citrate and succinate. As predicted by our model, bacteroids would thus have limited capability for 2-oxoglutarate synthesis and would consequently be restricted for ammonia assimilation into glutamate. In agreement with these findings, proteome and transcriptome data showed a strong down-regulation of isocitrate dehydrogenase and the GS-GOGAT system in bacteroids (figs. S2 and S3).

DISCUSSION
In this study, we combined experimental and computational methods to provide explanations for fundamental properties of bacteroid metabolism, which have so far mostly been investigated separately.
Our modeling results indicate that oxygen limitation is the driving factor behind the observed synthesis of carbon polymers (8,12) and induces the secretion of alanine, in addition to ammonia. The level of alanine secretion is dependent on the ratio of carbon and oxygen supply, which agrees with experimental studies (44). Rhizobial nitrogen fixation is fueled by dicarboxylates (3), requires a lowoxygen environment to protect nitrogenase (10), and involves downregulation of ammonia assimilation into glutamate in bacteroids (18). Both model predictions and metabolic flux analysis of A. caulinodans bacteroids indicate that these factors are interconnected: Catabolism of dicarboxylates induces a highly reduced redox poise and creates a high oxygen demand. In the low-oxygen environment of the nodule, this promotes down-regulation of the decarboxylating TCA cycle arm, which would decrease ammonia assimilation into glutamate by bacteroids. This agrees with glutamate levels being 20-fold lower in bacteroids compared to free-living rhizobia (8). In addition to promoting ammonia secretion to the plant, limited ammonia assimilation into glutamate would contribute to the growth arrest of bacteroids. As glutamate is the transamination donor for most other amino acids, this could explain the dependence of bacteroids on supply of amino acids by the plant (49,50). Intriguingly, the increased oxygen stress on bacteroids and their limited ability to assimilate ammonia into glutamate will force ammonia secretion and is therefore an important mechanism for enforcing a mutualistic relationship with the host plant. Furthermore, the ample supply of reducing equivalents would then be available to satisfy the high electron demand of the nitrogenase reaction, while carbon polymer synthesis and secretion of alanine, in particular, balance the allocation of carbon and regeneration of electron carriers. High oxygen consumption rates induced by dicarboxylate catabolism also contribute to lowering local oxygen levels in bacteroids, which would additionally protect the oxygen-sensitive nitrogenase enzyme. The proposed metabolic principles provide a general framework for understanding the constraints on bacteroid metabolism that favor ammonia (and amino acid) secretion and contribute to the growth arrest of bacteroids. However, rhizobium-legume symbioses are highly evolved mutualisms. As a result, many factors, such as nodule cysteine-rich peptides, play a role in inducing terminal differentiation of some bacteroids (51) and ammonia secretion by bacteroids is further forced by transcriptional regulation of genes involved in nitrogen metabolism (e.g., Ntr) (52).
Previous modeling studies of intermicrobial interactions found anoxic conditions to promote mutualistic interactions and increase the diversity of secreted metabolites (53,54), and oxygen availability was determined to be a better indicator of secreted metabolites than species identity (54). It is therefore possible that metabolic principles similar to those outlined in this study have a broader significance in governing metabolite exchanges in interspecies interactions.

General modeling procedures
Standard FBA calculations were performed in MATLAB R2019b (Mathworks) using scripts from the COBRA Toolbox v3.0 (55) and the Gurobi 8.0.1 solver (www.gurobi.com). When using the optimizeCbModel function, the Taxicab norm was minimized to avoid loops in the calculated flux distributions. The functions phenotypePhasePlane and fluxVariability were used for phenotype phase plane analysis and flux variability analysis, respectively. For flux variability analysis, loopless solutions were calculated that allowed for maximum nitrogenase activity. In silico gene deletion analysis was performed with the function singleGeneDeletion using the FBA method. A gene was considered essential when its deletion prevented flux through the nitrogenase reaction (data S6). All MATLAB scripts are available on Zenodo (https://doi.org/10.5281/ zenodo.4911594). Flux maps were created using Escher (56).

ECM analysis
ECMs were calculated using ecmtool (https://github.com/tjclement/ ecmtool) (32). Boundary conditions for calculating conversion modes were based on experimental studies and are detailed in data S8. To limit our analysis to biologically meaningful scenarios, only conversion modes with a positive cost lower than 40 g of carbon per gram of nitrogen were considered. We further restricted our analysis to conversion modes using one amino acid at most as an input. The  full set of conversion modes calculated with succinate, malate, and amino acids as inputs is available in data S7. For determining the correlation of oxygen demand and GS activity, glutamate was set as an additional output and a virtual metabolite was added to the GS reaction to track the flux through this reaction. ECM analysis for iCC541, a metabolic model of S. fredii bacteroids (19), was performed according to similar principles. To decrease computation time, cofactors were hidden in this analysis (data S8).

Statistical analysis
Nonparametric Spearman correlation and analysis of variance (ANOVA) followed by Tukey's multiple comparisons test were performed in GraphPad Prism 8.4.3. P < 0.05 was considered statistically significant.

Flux balance analysis
Constraints for FBA-based computations were defined similar to those for conversion mode analysis. For all proteinogenic amino acids except for asparagine, alanine, and aspartate, demand fluxes were constrained to values between 0.01 and 0.05 to mimic low levels of protein synthesis in bacteroids. The upper bound for alanine and aspartate demand reactions was left unconstrained because those amino acids are secreted by bacteroids. All constraints are detailed in data S8.
Ensemble-evolutionary FBA was performed as previously described (47,57). Briefly, 50,000 objective functions containing a random number of model reactions that were assigned random weights were generated (see fig. S10 for determination of ensemble size). Flux distributions for all objective functions were then calculated using FBA for a fixed malate uptake rate of 4 flux units and varying oxygen uptake rate. Constraints for ensemble-evolutionary FBA were defined as for standard FBA, but with a minimum flux of 0.01 for the nitrogenase reaction and no upper bound on all amino acid demand reactions.

Model reconstruction
A database of gene-protein-reaction associations was derived from genome annotation of R. leguminosarum bv. viciae 3841 (Rlv3841). First, an orthology-based reconstruction was obtained using AuReMe (58). Because both a genome-scale model [iGD1575b (21,59)] and a highly curated core model [iGD726 (59)] for the rhizobial strain S. meliloti 1021 are available, these were used as templates. Reconstruction in AuReMe was performed for both S. meliloti models, resulting in draft networks containing 1304 and 613 reactions for iGD1575b and iGD726, respectively. Where the same reaction occurred in both reconstructions, the gene-protein-reaction association from the iGD726-based reconstruction was selected as it can be expected to be more accurate. To account for gene functions annotated in the Rlv3841 genome but not present in S. meliloti, a second draft metabolic model was obtained from KBase (60). Briefly, the Rlv3841 genome was annotated with Prokka (v1.12) and the "Build Metabolic Model" function was run without gap-filling. This generated a network with 607 reactions, 98 of which were not contained in the reaction list obtained from the AuReMe reconstructions. A model of bacteroid metabolism in Rlv3841 was then built by manually selecting reactions from the draft reconstructions that are catalyzed by enzymes detected in the bacteroid proteome, as well as those associated with up-regulated (34) and essential genes for bacteroids (35). Gene-protein-reaction associations were refined by comparison to gene essentiality data for free-living Rlv3841 (61). Gaps were filled on the basis of literature evidence and guided by the Kyoto Encyclopedia of Genes and Genomes database (62). Because the focus of this study was on central metabolic pathways, biosynthesis of some cofactors was excluded and GTP (guanosine triphosphate) was replaced with ATP where applicable. Some linear pathways, namely, most reactions associated with lipid biosynthesis and myo-inositol catabolism, were summarized in lumped reactions. Further details on pathways included in the model and the definition of exchange fluxes are provided in the Supplementary Text.

Bacterial strains and culture conditions
Oxygen consumption and NADH/NAD + ratios were measured for Rlv3841. Cultures were grown in universal minimal salts (UMS) (61) or acid minimal salts (AMS) (63) medium at 28°C. UMS was supplemented with different carbon and nitrogen sources at the following final concentrations: succinate, 20 mM; arabinose, 20 mM; glucose, 10 mM; and NH 4 Cl, 10 mM. Cultivations were performed in 250-ml Erlenmeyer flasks with an initial filling volume of 50 ml. Low-oxygen cultivations were performed in a glove box (Belle Technology) with the atmosphere adjusted to the desired oxygen concentration by flushing with nitrogen gas.
A. caulinodans ORS571 was grown in UMS medium supplemented with 10 mM succinate and 0.3 mM nicotinic acid at 37°C. For diazotrophic growth, cultures were continuously sparged with a gas mixture containing 97% N 2 and 3% O 2 , which had been determined to be the optimal oxygen level for diazotrophic growth.

Measurement of NADH/NAD + ratios
NADH/NAD + ratios were determined using the NAD/NADH-Glo Assay (Promega) according to the manufacturer's instructions. Cells were harvested during exponential phase [OD 600 (optical density at 600 nm) = 0.4 to 0.6] for all measurements.

Measurement of oxygen consumption rates
Oxygen consumption rates were measured for liquid cultures in early exponential phase. Cultures were diluted to OD 600 ~ 0.15, and a 25-ml glass vial containing an OxyDot was filled to the top with liquid culture and sealed. Measurements of oxygen concentration were taken every 15 s using the OxySense 325i system while the culture was stirred with a magnetic stirrer. The data were analyzed with the OxySense Gen III software and oxygen consumption was calculated from measurements obtained between 18 and 15% oxygen concentration.

Proteomics
Bacteroids were obtained from nodules of pea (Pisum sativum cv. Avola) plants inoculated with Rlv3841 and grown in an illuminated environment-controlled growth room as previously described (8). Plants were harvested at 28 days after inoculation, with bacteroids extracted from excised root nodules by double Percoll gradient purification (14). Bacteroids were isolated from nodules obtained from a total of 21 plants and processed as three independent replicates.
For the free-living cultures, Rlv3841 was grown in AMS medium with 20 mM succinate and 10 mM 15 NH 4 Cl. Cultures were grown to late log phase (OD 600 ~ 0.8) at 28°C on a gyratory shaker at 250 rpm and subcultured into fresh AMS six times, as preliminary trials showed that this yielded >99% 15 N incorporation into cell proteins.
From the sixth subculture, three separate AMS cultures were inoculated and harvested at mid-log phase (OD 600 ~ 0.4). Cells were harvested by centrifugation, washed in AMS, and stored at −80°C. Bacteroid and cell pellets were resuspended separately in 10 mM HEPES (pH 7.2) buffer and an aliquot of each was taken and lysed by two rounds on a FastPrep FP120 Ribolyser (Bio 101/Savant) at a setting of 6.5 for 30 s, with samples kept on ice for 5 min between each round. The protein content of the extracted aliquots was then determined by Bradford assay, where bovine serum albumin (Pierce) was used to generate a standard curve. These values were then used to mix equivalent proportions of unextracted bacteroid and cells from the original suspensions to yield at least 200 g/ml of total combined protein. Mixed bacteroid and cell samples were then ribolyzed as described above, another Bradford determination was performed to confirm the protein concentration, and then the equivalent of 50 g of protein was extracted by methanol-chloroformwater precipitation (64).
The protein pellets from the three replicates were dissolved in SDS-gel sample loading buffer, heated at 80°C for 10 min, and loaded onto a Novex gel (10% bis-tris SDS gel, Life Technologies, Carlsbad, CA). After separation and staining with InstantBlue (Expedeon, Harston, UK), the gel lanes were cut into 12 to 15 slices that were washed, reduced and alkylated, and treated with trypsin according to standard procedures. After digestion, peptides were extracted with 5% formic acid/50% acetonitrile, dried, and redissolved in 0.1% trifluoroacetic acid.
Liquid chromatography-tandem mass spectrometry analysis of all gel fractions was performed on an LTQ-Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA) coupled with a nanoACQUITY UPLC system (Waters, Manchester, UK). Sample aliquots were loaded onto a trap column (Symmetry R C18, 5 m, 180 m by 20 mm; Waters), and the peptides were then separated on an analytical column (BEH C18, 1.7 m, 75 m by 250 mm; Waters) and infused into the mass spectrometer via a 10-m SilicaTip nanospray emitter (New Objective, Woburn, MA) attached to a nanospray interface (Proxeon, Odense, Denmark).
For separation, the following gradient of solvents A (0.1% formic acid in water) and B (0.1% formic acid in acetonitrile) was used at a flow rate of 250 nl min −1 : solvent B at start: 0%; 0 to 3 min: linear ramp to 5% B; 3 to 56 min: ramp to 40% B; 56 to 62 min: ramp to 85% B; 85% B kept for 3 min followed by 100% A for 20 min for reequilibration.
The mass spectrometer was operated in positive ion mode at a capillary temperature of 200°C. The source voltage and focusing voltages were tuned for the transmission of MRFA peptide [mass/ charge ratio (m/z) 524] (Sigma-Aldrich, St. Louis, MO). Datadependent analysis was carried out in Orbitrap-IT parallel mode using CID fragmentation on the five most abundant ions in each cycle. The full scan mass spectrometry was performed in the Orbitrap at a resolution of 60,000 over the range m/z 350 to 1800. For the CID-MS2, the mono-isotopic 2+ and 3+ charged precursors were selected with an isolation width of 2 Da. MS2 was triggered by a minimal signal of 10 3 with an AGC target of 3 × 10 4 ions and 150-ms maximum scan time using the chromatography function for peak apex detection. Collision energy was 35, and dynamic exclusion was set to 1 count and 60 s with a mass window of ±20 parts per million (ppm). Mass spectrometry scans were saved in profile mode while MS2 scans were saved in centroid mode.
The analysis of all gel fractions from three replicates resulted in 52 raw files (with replicate 2 run twice). Data were processed using Mascot Distiller 2.7 and Mascot Server 2.7 through Mascot Daemon 2.7 (Matrix Science, London, UK). Peak lists generated by Mascot Distiller were used for a database search using Mascot Server on the GeneDB_Rleguminosarum_Proteins fasta database from www.sanger. ac.uk/resources/downloads/bacteria/rhizobium-leguminosarum.html (May 2020, 7144 entries). For protein annotation, data from https:// rhizosphere.org/lab-page/molecular-tools/genomes/rlv3841-genome/ (May 2020, 7288 entries) were used. A small database containing common contaminants (MaxQuant contaminants 2017, 250 entries) was included in the search.
Each raw file was processed and searched separately using the enzyme trypsin with two missed cleavages, a precursor mass tolerance of 10 ppm, and a fragment mass tolerance of 0.6 Da. Carbamidomethyl (C) was set as a fixed modification, and oxidation (M), deamidation (N, Q), and acetylation (protNterm) were set as variable modifications. Intensities for light and heavy 15 N-labeled peptides were extracted using a Mascot Server 15 N metabolic quantitation method with the following parameters: 99.2% labeling (as determined by Mascot Distiller), Simpson's integration, isotope match rho = 0.6, XIC threshold = 0.1, isolation threshold = 0.5, peptide expect threshold = 0.05, outlier removal auto, and normalization none. The experimental design with three replicates and corresponding fractions was generated in the quantitation table exported via Mascot Daemon. The resulting expression table was used for ratio calculation and statistical analysis in RStudio 1.2.5033 with R version 3.6.3 (65). Potential contaminants were removed, and the table was filtered for proteins quantified in all three replicates. Light and heavy protein intensities were log 10 -transformed, and ratios were calculated as light/heavy (bacteroids/Rhizobium cells) for each replicate. Those ratios were used for statistical testing using the limma eBayes function in R. The final ratio was calculated as median from the replicates. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (66) partner repository with the dataset identifier PXD019467.

Metabolic flux analysis
Metabolic flux analysis of A. caulinodans ORS571 was performed as previously described (8). Briefly, A. caulinodans was grown in UMS medium containing 20% [ 13 C 4 ]-succinate and 80% unlabeled sodium succinate. Bacteroids were isolated from the root nodules of S. rostrata as previously described (67). Isolated bacteroids were incubated in UMS medium supplemented with 10 mM 20% [ 13 C 4 ]-succinate for 24 hours at ≤1% O 2 . Nitrogenase activity of free-living cultures and isolated bacteroids during labeling experiments was determined by an acetylene reduction assay (8). Free-living cultures were harvested in late exponential phase and bacteroids after 24 hours for metabolite extraction and gas chromatography-mass spectrometry analysis. Mass spectrometry data processing and isotopomer analysis of protein-derived amino acids and hydroxybutyrate obtained by hydrolysis of PHB was done using methods described previously (8,68).
Metabolic modeling was performed with INCA (Isotopomer Network Compartmental Analysis) using an iterative procedure (69). A complete description of the model, including the network carbon atom transitions, and net flux data are provided in data S9. The model along with the mass spectrometry measurements was simulated to obtain an optimized flux pattern in the network, followed by statistical validation of the flux maps. The fluxes were estimated relative to the succinate uptake flux fixed at 100 units. The assessment of goodness of fit for the flux maps with statistically valid sum of squared residuals was performed by the comparison of simulated mass spectrometry measurements with that of experimentally measured values. To assess the precision of flux estimates, parameter continuation was performed in INCA to calculate the lower and upper bounds of the 95% confidence interval for the flux estimates. Comparison of fluxes estimated for the free-living bacteria and bacteroids were based on the confidence intervals (upper and lower limits) for a specific flux. A flux determined under two different conditions was deemed to be significantly different if the confidence intervals for the flux under the two conditions did not overlap.

SUPPLEMENTARY MATERIALS
Supplementary material for this article is available at http://advances.sciencemag.org/cgi/ content/full/7/31/eabh2433/DC1 View/request a protocol for this paper from Bio-protocol.