Abstract While nitrogen fertilizers are widely used in agricultural production, their application incurs significant environmental and energetic costs. In contrast, some crops are less dependent on these fertilizers because they engage in symbioses with rhizobia, nitrogen-fixing bacteria provide ammonium to the plant in exchange for carbon. However, the carbon cost associated with nitrogen fixation can negatively impact crop yields. Improving the efficiency of this metabolic process could alleviate this impact on crop productivity. Mathematical models can help us quantitatively explore metabolic behavior and identify potential targets for metabolic engineering. In this work, we developed a kinetic model of determinate root nodule metabolism, where this symbiotic exchange of carbon from the plant and nitrogen from the bacteria occurs. We used this model to evaluate how the predicted metabolic behavior differs between inefficient and efficient nodules, and to identify potential engineering targets for improving nitrogen fixation efficiency and rate. We show that the enzymes phosphoenolpyruvate carboxylase and pyruvate kinase have significant influence on the predicted rate and efficiency of nitrogen fixation, especially when their expression is varied in combination with oxidative Pentose Phosphate Pathway enzymes like glucose-6-phosphate dehydrogenase and 6-phosphogluconolactonase. The model predicts that pairing a 3-fold decrease in glucose-6-phosphate dehydrogenase activity along with either a 3-fold increase in phosphoenolpyruvate carboxylase activity or decrease in pyruvate kinase activity could increase nitrogen fixation rate by 5.51% while improving nitrogen fixation efficiency by 7.74%.

Engineering a novel N2-fixing rhizobia symbiosis in cereal crops is a strategy being pursued to improve agricultural sustainability. However, if such a symbiosis were introduced, it would have to be economically viable in the context of plants' existing nutrient acquisition strategies, including existing symbioses with arbuscular mycorrhizal fungi (AMF) that most plants already engage in. It is important to understand how the metabolic costs and benefits from these symbioses with overlapping functions might impact plant growth when evaluating the potential benefits of this engineering strategy. To address this, we developed metabolic models describing how the relative growth rate of Zea mays is impacted by the AMF Rhizophagus irregularis and a hypothetical N2-fixing symbiosis with Bradyrhizobium diazoefficiens in isolation and in tandem. The metabolic models of the plant-AMF symbiosis and plant-AMF-rhizobia symbiosis are the first of their kind. To validate the AMF component of our model, we conducted a field evaluation comparing AMF-compatible and mutant AMF-incompatible maize hybrids. The empirically measured AMF-mediated growth benefit agreed well with model predictions. Our model of the rhizobium symbiosis predicted that the lower N content of cereal crops makes the growth penalty associated with acquiring nitrogen from rhizobia smaller than in legumes. Finally, the model of the plant-AMF-rhizobia symbiosis predicted positive synergies between rhizobia and AMF under nutrient-limited conditions but negative synergies under phosphorus-replete conditions. This indicates that these bioengineering strategies could improve cereal crop yields and may achieve greater gains in tandem, but soil nutrient levels and plant nitrogen requirements should be considered.

Photorespiration is the second largest carbon flux in most leaves and is integrated into metabolism broadly including one-carbon (C1) metabolism. Photorespiratory intermediates such as serine and others may serve as sources of C1 units, but it is unclear to what degree this happens in vivo, whether altered photorespiration changes flux to C1 metabolism, and if so through which intermediates. To clarify these questions, we quantified carbon flux from photorespiration to C1 metabolism using 13CO2 labelling and isotopically non-stationary metabolic flux analysis in Arabidopsis thaliana under different O2 concentrations which modulate photorespiration. The results revealed that 5.8

Engineering N2-fixing rhizobia symbiosis in cereals and enhancing association with arbuscular mycorrhizal fungi (AMF) are two strategies being pursued to improve agricultural sustainability. However, the metabolic fluxes and consequent costs and benefits of these strategies in a cereal crop have not been examined. We developed a metabolic model describing how the relative growth rate of Zea mays is impacted by hypothetical association with the rhizobium Bradyrhizobium diazoefficiens and the AMF Rhizophagus irregularis. We found that the lower N content of cereal crops makes the relative growth rate cost associated with acquiring nitrogen from N2-fixing rhizobia smaller than in legumes. Our model also predicts positive synergies between rhizobia and AMF under nutrient-limited conditions but negative synergies under nutrient, particularly phosphorus, replete conditions. These findings indicate that these bioengineering strategies could improve cereal crop yields and may achieve greater gains in tandem, but soil nutrient status of target sites should be considered. ### Competing Interest Statement The authors have declared no competing interest.

In the future, plants may encounter increased light and elevated CO 2 levels. How consequent alterations in photosynthetic rates will impact fluxes in photosynthetic carbon metabolism remains uncertain. Respiration in light ( R L ) is pivotal in plant carbon balance and a key parameter in photosynthesis models. Understanding the dynamics of photosynthetic metabolism and R L under varying environmental conditions is essential for optimizing plant growth and agricultural productivity. However, measuring R L under high light and high CO 2 (HLHC) conditions poses challenges using traditional gas exchange methods. In this study, we employed isotopically nonstationary metabolic flux analysis (INST-MFA) to estimate R L and investigate photosynthetic carbon flux, unveiling nuanced adjustments in Camelina sativa under HLHC. Despite numerous flux alterations in HLHC, R L remained stable. HLHC affects several factors influencing R L , such as starch and sucrose partitioning, v o / v c ratio, triose phosphate partitioning, and hexose kinase activity. Analysis of A/C i curve operational points reveals that HLHC's major changes primarily stem from CO 2 suppressing photorespiration. Integration of these fluxes into a simplified model predicts changes in CBC labeling under HLHC. This study extends our prior discovery that incomplete CBC labeling is due to unlabeled carbon reimported during R L , offering insights into manipulating labeling through adjustments in photosynthetic rates.

The thermoacidophilic red alga Cyanidioschyzon merolae survives its challenging environment likely in part by operating a carbon-concentrating mechanism (CCM). Here, we demonstrated that C. merolae's cellular affinity for CO2 is stronger than the affinity of its rubisco for CO2. This finding provided additional evidence that C. merolae operates a CCM while lacking the structures and functions characteristic of CCMs in other organisms. To test how such a CCM could function, we created a mathematical compartmental model of a simple CCM, distinct from those we have seen previously described in detail. The results of our modeling supported the feasibility of this proposed minimal and non-canonical CCM in C. merolae. To facilitate the robust modeling of this process, we measured and incorporated physiological and enzymatic parameters into the model. Additionally, we trained a surrogate machine-learning model to emulate the mechanistic model and characterized the effects of model parameters on key outputs. This parameter exploration enabled us to identify model features that influenced whether the model met the experimentally derived criteria for functional carbon concentration and efficient energy usage. Such parameters included cytosolic pH, bicarbonate pumping cost and kinetics, cell radius, carboxylation velocity, number of thylakoid membranes, and CO2 membrane permeability. Our exploration thus suggested that a non-canonical CCM could exist in C. merolae and illuminated the essential features generally necessary for CCMs to function. The thermoacidophilic red alga Cyanidioschyzon merolae possesses a unique carbon-concentrating mechanism (CCM), and a model illuminates the essential features generally required for CCM function.

Carbon-concentrating mechanisms (CCMs) have evolved numerous times in photosynthetic organisms. They elevate the concentration of CO2 around the carbon-fixing enzyme rubisco, thereby increasing CO2 assimilatory flux and reducing photorespiration. Biophysical CCMs, like the pyrenoid-based CCM (PCCM) of Chlamydomonas reinhardtii or carboxysome systems of cyanobacteria, are common in aquatic photosynthetic microbes, but in land plants appear only among the hornworts. To predict the likely efficiency of biophysical CCMs in C3 plants, we used spatially resolved reaction-diffusion models to predict rubisco saturation and light use efficiency. We found that the energy efficiency of adding individual CCM components to a C3 land plant is highly dependent on the permeability of lipid membranes to CO2, with values in the range reported in the literature that are higher than those used in previous modeling studies resulting in low light use efficiency. Adding a complete PCCM into the leaf cells of a C3 land plant was predicted to boost net CO2 fixation, but at higher energetic costs than those incurred by photorespiratory losses without a CCM. Two notable exceptions were when substomatal CO2 levels are as low as those found in land plants that already use biochemical CCMs and when gas exchange is limited, such as with hornworts, making the use of a biophysical CCM necessary to achieve net positive CO2 fixation under atmospheric CO2 levels. This provides an explanation for the uniqueness of hornworts' CCM among land plants and the evolution of pyrenoids multiple times.

13C-Metabolic Flux Analysis (13C-MFA) and Flux Balance Analysis (FBA) are widely used to investigate the operation of biochemical networks in both biological and biotechnological research. Both of these methods use metabolic reaction network models of metabolism operating at steady state, so that reaction rates (fluxes) and the levels of metabolic intermediates are constrained to be invariant. They provide estimated (MFA) or predicted (FBA) values of the fluxes through the network in vivo, which cannot be measured directly. A number of approaches have been taken to test the reliability of estimates and predictions from constraint-based methods and to decide on and/or discriminate between alternative model architectures. Despite advances in other areas of the statistical evaluation of metabolic models, validation and model selection methods have been underappreciated and underexplored. We review the history and state-of-the-art in constraint-based metabolic model validation and model selection. Applications and limitations of the χ2-test of goodness-of-fit, the most widely used quantitative validation and selection approach in 13C-MFA, are discussed, and complementary and alternative forms of validation and selection are proposed. A combined model validation and selection framework for 13C-MFA incorporating metabolite pool size information that leverages new developments in the field is presented and advocated for. Finally, we discuss how the adoption of robust validation and selection procedures can enhance confidence in constraint-based modeling as a whole and ultimately facilitate more widespread use of FBA in biotechnology in particular.

The modeling of rates of biochemical reactions-fluxes-in metabolic networks is widely used for both basic biological research and biotechnological applications. A number of different modeling methods have been developed to estimate and predict fluxes, including kinetic and constraint-based (Metabolic Flux Analysis and flux balance analysis) approaches. Although different resources exist for teaching these methods individually, to-date no resources have been developed to teach these approaches in an integrative way that equips learners with an understanding of each modeling paradigm, how they relate to one another, and the information that can be gleaned from each. We have developed a series of modeling simulations in Python to teach kinetic modeling, metabolic control analysis, 13C-metabolic flux analysis, and flux balance analysis. These simulations are presented in a series of interactive notebooks with guided lesson plans and associated lecture notes. Learners assimilate key principles using models of simple metabolic networks by running simulations, generating and using data, and making and validating predictions about the effects of modifying model parameters. We used these simulations as the hands-on computer laboratory component of a four-day metabolic modeling workshop and participant survey results showed improvements in learners' self-assessed competence and confidence in understanding and applying metabolic modeling techniques after having attended the workshop. The resources provided can be incorporated in their entirety or individually into courses and workshops on bioengineering and metabolic modeling at the undergraduate, graduate, or postgraduate level.

Motivation The accurate prediction of complex phenotypes such as metabolic fluxes in living systems is a grand challenge for systems biology and central to efficiently identifying biotechnological interventions that can address pressing industrial needs. The application of gene expression data to improve the accuracy of metabolic flux predictions using mechanistic modeling methods such as Flux Balance Analysis (FBA) has not been previously demonstrated in multi-tissue systems, despite their biotechnological importance. We hypothesized that a method for generating metabolic flux predictions informed by relative expression levels between tissues would improve prediction accuracy. Results Relative gene expression levels derived from multiple transcriptomic and proteomic datasets were integrated into Flux Balance Analysis predictions of a multi-tissue, diel model of Arabidopsis thaliana’s central metabolism. This integration dramatically improved the agreement of flux predictions with experimentally based flux maps from 13C Metabolic Flux Analysis (MFA) compared with a standard parsimonious FBA approach. Disagreement between FBA predictions and MFA flux maps, as measured by weighted averaged percent error values, dropped from between 169-180% and 94-103% in high light and low light conditions, respectively, to between 10-12% and 9-11%, depending on the gene expression dataset used. The incorporation of gene expression data into the modeling process also substantially altered the predicted carbon and energy economy of the plant. Availability Code is available from https://github.com/Gibberella/ArabidopsisGeneExpressionWeights Contact yairhill@msu.edu

Since they emerged approximately 125 million years ago, flowering plants have evolved to dominate the terrestrial landscape and survive in the most inhospitable environments on earth. At their core, these adaptations have been shaped by changes in numerous, interconnected pathways and genes that collectively give rise to emergent biological phenomena. Linking gene expression to morphological outcomes remains a grand challenge in biology, and new approaches are needed to begin to address this gap. Here, we implemented topological data analysis (TDA) to summarize the high dimensionality and noisiness of gene expression data using lens functions that delineate plant tissue and stress responses. Using this framework, we created a topological representation of the shape of gene expression across plant evolution, development, and environment for the phylogenetically diverse flowering plants. The TDA-based Mapper graphs form a well-defined gradient of tissues from leaves to seeds, or from healthy to stressed samples, depending on the lens function. This suggests that there are distinct and conserved expression patterns across angiosperms that delineate different tissue types or responses to biotic and abiotic stresses. Genes that correlate with the tissue lens function are enriched in central processes such as photosynthetic, growth and development, housekeeping, or stress responses. Together, our results highlight the power of TDA for analyzing complex biological data and reveal a core expression backbone that defines plant form and function.

Carbon-neutral production of valuable bioproducts is critical to sustainable development but remains limited by the slow engineering of photosynthetic organisms. Improving existing synthetic biology tools to engineer model organisms to fix carbon dioxide is one route to overcoming the limitations of photosynthetic organisms. In this work, we describe a pipeline that enabled the deletion of a conditionally essential gene from the Shewanella oneidensis MR-1 genome. S. oneidensis is a simple bacterial host that could be used for electricity-driven conversion of carbon dioxide in the future with further genetic engineering. We used Flux Balance Analysis (FBA) to model carbon and energy flows in central metabolism and assess the effects of single and double gene deletions. We modeled the growth of deletion strains under several alternative conditions to identify substrates that restore viability to an otherwise lethal gene knockout. These predictions were tested in vivo using a Mobile-CRISPRi gene knockdown system. The information learned from FBA and knockdown experiments informed our strategy for gene deletion, allowing us to successfully delete an "expected essential" gene, gpmA. FBA predicted, knockdown experiments supported, and deletion confirmed that the "essential" gene gpmA is not needed for survival, dependent on the medium used. Removal of gpmA is a first step toward driving electrode-powered CO2 fixation via RuBisCO. This work demonstrates the potential for broadening the scope of genetic engineering in S. oneidensis as a synthetic biology chassis. By combining computational analysis with a CRISPRi knockdown system in this way, one can systematically assess the impact of conditionally essential genes and use this knowledge to generate mutations previously thought unachievable.

Significance Photosynthesis metabolites are quickly labeled when 13CO2 is fed to leaves, but the time course of labeling reveals additional contributing processes involved in the metabolic dynamics of photosynthesis. The existence of three such processes is demonstrated, and a metabolic flux model is developed to explore and characterize them. The model is consistent with a slow return of carbon from cytosolic and vacuolar sugars into the Calvin–Benson cycle through the oxidative pentose phosphate pathway. Our results provide insight into how carbon assimilation is integrated into the metabolic network of photosynthetic cells with implications for global carbon fluxes.

Since they emerged ~125 million years ago, flowering plants have evolved to dominate the terrestrial landscape and survive in the most inhospitable environments on earth. At their core, these adaptations have been shaped by changes in numerous, interconnected pathways and genes that collectively give rise to emergent biological phenomena. Linking gene expression to morphological outcomes remains a grand challenge in biology, and new approaches are needed to begin to address this gap. Here, we implemented topological data analysis (TDA) to summarize the high dimensionality and noisiness of gene expression data using lens functions that delineate plant tissue and stress responses. Using this framework, we created a topological representation of the shape of gene expression across plant evolution, development, and environment for the phylogenetically diverse flowering plants. The TDA-based Mapper graphs form a well-defined gradient of tissues from leaves to seeds, or from healthy to stressed samples, depending on the lens function. This suggests there are distinct and conserved expression patterns across angiosperms that delineate different tissue types or responses to biotic and abiotic stresses. Genes that correlate with the tissue lens function are enriched in central processes such as photosynthetic, growth and development, housekeeping, or stress responses. Together, our results highlight the power of TDA for analyzing complex biological data and reveal a core expression backbone that defines plant form and function. Significance statement A grand challenge in biology is to link gene expression to phenotypes across evolution, development, and the environment, but efforts have been hindered by biological complexity and dataset heterogeneity. Here, we implemented topological data analysis across thousands of gene expression datasets in phylogenetically diverse flowering plants. We created a topological representation of gene expression across plants and observed well-defined gradients of tissues from leaves to seeds, or from healthy to environmentally stressed. Using this framework, we identified a core and deeply conserved expression backbone that defines plant form and function, with key patterns that delineate plant tissues, abiotic, and biotic stresses. Our results highlight the power of topological approaches for analyzing complex biological datasets.

Epichloë species (Clavicipitaceae, Ascomycota) are endophytic symbionts of many cool-season grasses. Many interactions between Epichloë and their host grasses contribute to plant growth promotion, protection from many pathogens and insect pests, and tolerance to drought stress. Resistance to insect herbivores by endophytes associated with Hordeum species has been previously shown to vary depending on the endophyte-grass-insect combination. We explored the genetic and chemotypic diversity of endophytes present in wild Hordeum species. We analyzed seeds of Hordeum bogdanii, H. brevisubulatum, and H. comosum obtained from the US Department of Agriculture's (USDA) National Plant Germplasm System (NPGS), of which some have been reported as endophyte-infected. Using polymerase chain reaction (PCR) with primers specific to Epichloë species, we were able to identify endophytes in seeds from 17 of the 56 Plant Introduction (PI) lines, of which only 9 lines yielded viable seed. Phylogenetic analyses of housekeeping, alkaloid biosynthesis, and mating type genes suggest that the endophytes of the infected PI lines separate into five taxa: Epichloë bromicola, Epichloë tembladerae, and three unnamed interspecific hybrid species. One PI line contained an endophyte that is considered a new taxonomic group, Epichloë sp. HboTG-3 (H. bogdanii Taxonomic Group 3). Phylogenetic analyses of the interspecific hybrid endophytes from H. bogdanii and H. brevisubulatum indicate that these taxa all have an E. bromicola allele but the second allele varies. We verified in planta alkaloid production from the five genotypes yielding viable seed. Morphological characteristics of the isolates from the viable Hordeum species were analyzed for their features in culture and in planta. In the latter, we observed epiphyllous growth and in some cases sporulation on leaves of infected plants.