Environment-coupled models of leaf metabolism

The plant leaf is the main photosynthetic organ where light energy is converted into chemical energy to support the maintenance and growth of the leaf and the rest of the plant. The leaf is also the site of gas-exchange (O₂, CO₂ and water vapour) through the stomata and different forms of photosynthesis and leaf architecture — C3, C4 and crassulacean acid metabolism (CAM) — have evolved to cope with different environments. However, given its large surface area the plant leaf is also vulnerable to pathogen attacks. Due to these various factors, the leaf is a main determinant of plant growth and health. Therefore, gaining a better understanding of leaf–environment interactions is key to developing strategies to better equip our crop plants for productivity requirements in current and future climates. Computational models of leaf metabolism can aid this quest by explaining metabolic, architectural, and evolutionary aspects to guide the engineering of more productive and resistant crop species. In this review, I focus on flux-balance models of leaf metabolism and highlight their applications to study plant–environment interactions.

Flux-balance models rely on the stoichiometry of the metabolic network under consideration and an optimality assumption based on which metabolic steady-state fluxes can be predicted [1,2]. In the past decade, plant stoichiometric network models have been reconstructed for a large range of species and to varying extent, covering cell suspension models, multi-cell type, and whole-plant models [3–6]. Several pipelines for the automated reconstruction of plant metabolic network models are available [7–10]. The optimality assumption or ‘objective function’ depends on the metabolic system under consideration and reflects its function in the plant. For instance, for a growing leaf the objective function could be the maximization of biomass production and for a mature leaf the export of photosynthetic assimilates to the phloem [11]. Other optimality criteria have also been applied, such as the minimization of internal fluxes or photon usage efficiency [12,13]. The latter objectives are often chosen when metabolic outputs, such as growth rates, have been experimentally determined and can be applied as constraints. Additionally, external constraints, such as limited nutrient or light availability, can be modelled by constraining the respective uptake fluxes. Depending on these constraints the model then predicts optimal metabolic steady-state flux modes which can then be analyzed with respect to changes between conditions and complement experimental observations. See Figure 1 for a schematic representation.

As an extension to this, various data-integrative approaches have been proposed to model condition-specific metabolic fluxes or to extract context-specific metabolic networks. Most of these approaches employ transcript, protein or metabolite abundances to constraint the flux boundaries of the respective reactions [14–17].

While any flux-balance model of leaf metabolism is inherently coupled to the environment via exchange reactions, such as light influx, CO₂ and O₂ exchange, and nutrient uptake; the focus here will be on studies that investigate the leaf's metabolic response to changes in those parameters or which integrate data collected from leaf material exposed to different environments. The highlighted studies are divided into two groups — environmental constraint-driven and data-integrative — based on how the environment-specificity is achieved. Moreover, the studies are organized according to the modelled interaction, i.e. abiotic or biotic. A schematic overview of the various modelling approaches is shown in Figure 2 and summarized in Table 1.

In one of the earliest studies in 2013, Poolman et al. [18] analyzed a genome-scale model representing a developing rice leaf cell to identify changes in reaction fluxes in response to changes in light influx. One of their findings put forward a beneficial role for the energy-consuming photorespiratory pathway under supraoptimal light conditions. The photorespiratory pathway recycles toxic intermediates which are formed due to the dual activity of ribulose-1,5-bisphosphate carboxylase-oxygenase (rubisco), the main enzyme to convert atmospheric carbon dioxide to energy-rich molecules. With increasing temperature rubisco's ratio of carboxylation to oxygenation (v_C/v_O) shifts towards O₂-fixation. This side reaction forms byproducts which are then recycled in the energy-wasting process of photorespiration. The model predicted that the photorespiratory pathway was spontaneously activated at supraoptimal light intensities and thereby dissipated excess light energy to prevent cell damage or overreduction.

A year later the same authors extended this study by exploring flux rearrangements during the transition from a growing to a mature leaf [11]. This was modelled as a transition from synthesizing leaf biomass to export of nutrients to the phloem as a function of light intensity. The authors found that the predicted flux patterns were relatively insensitive to the specific metabolic function of the leaf and attributed flux differences to changes in energy and redox requirements for the different metabolic outputs.

Cheung et al. [19] used the effect of light intensity on leaf metabolism together with a cost-weighted flux-balance formulation to study alternative metabolic pathways in Arabidopsis. The authors simulated C-limited conditions by fixing cellular output and maintenance costs and varied photon influx. They then simulated flux distributions for sets of randomly chosen flux weighting factors using flux minimization as the objective function. Their results emphasized the potential contribution of alternative fluxes in rebalancing and consuming ATP and NADPH, such as the chlorophyll and xanthophyll pigment cycles and various futile cycles, in different environmental and physiological conditions.

Later, Chatterjee et al. [20] used a different rice leaf model to study its responses to different light intensities and varying v_C/v_O ratios of rubisco, thereby mimicking increased photorespiration under drought stress, normal and suppressed photorespiration. Their comprehensive environment-scan combined with cost-weighted flux-balance analysis revealed different metabolic modes for maintaining redox and ATP balance including flux rearrangements across compartments and shifts in transport reactions. These findings highlighted alternative routes and metabolic flexibility for stress adaptation in leaf metabolism.

Plant species, such as maize, sugarcane, and sorghum perform C4 photosynthesis. This specialized form of photosynthesis involves a carbon-concentration mechanism in which initial CO₂ fixation and internal re-fixation by rubisco are separated in two adjacent cell types — the bundle sheath and mesophyll cells. This separation increases both nitrogen (N) and water use efficiency compared with C3 plants. To investigate the metabolic impact of N availability in maize Simons et al. [21] presented a C4 model capturing the interactions between bundle sheath and mesophyll cells. The authors used condition-specific biomass compositions and explored the effect of two knockdown mutants of glutamine synthetase (GS) — an enzyme essential for N assimilation. The authors used transcriptomic and proteomic data to introduce regulatory constraints that represented the metabolism of wild type and GS mutants under N-complete and -deficient conditions. This way, between nine and 100 reaction fluxes per condition were constrained. Using these condition-specific models they determined flux-sum ranges as a flux measure for the reactions associated with either the production or consumption of a given metabolite. They achieved up to 90% accuracy when comparing their model predictions with measured metabolite levels. Furthermore, the authors identified a set of genes for further testing which the model's predictions had related to changes in biomass formation in the considered conditions.

In the same year, Arnold and Nikoloski [22] presented an Arabidopsis core model equipped with three biomass compositions to simulate carbon-limited, N-limited, and optimal growth conditions. This core model and two previous Arabidopsis models [23,24] were later used to in silico analyze the effect of N supply (nitrate/ammonium) on individual amino acid synthesis costs under autotrophic/heterotrophic and day/night growth conditions [25]. The authors quantified the synthesis costs of amino acids in terms of ATP demand and found these costs to be highly dependent on the environment — most amino acid costs for night conditions were higher than those for heterotrophic day conditions and the costs for autotrophic conditions were higher than the costs for heterotrophic conditions. The study further confirmed $NH4+$ uptake to be cheaper than $NO3−$ uptake to meet the plant's N demand due to the extra cost for nitrate reduction.

Shaw and Cheung used their previously developed diel Arabidopsis [26] model to generate a sophisticated dynamic multi-tissue, day–night modelling framework in which they studied the effect of both N-rich and -limiting conditions on Arabidopsis leaf and root growth [27]. Their framework demonstrated the power of integrated whole-plant analysis to uncover optimal growth strategies. More specifically, their model showed that it is energetically most efficient to store nitrate taken up into the root during the night in the vacuole and to then transport and fix it in the leaf during the day when light energy is available.

Studies that investigated the effect of different CO₂ levels on leaf metabolism and particularly photorespiration include models of rice [20,28], tomato [29] and generic models of CAM photosynthesis [30]. Typically, constraints on rubisco's v_C/v_O ratio were implemented to model the leaf's metabolic response to changing CO₂ levels or temperature (as v_C/v_O is temperature-dependent, see also section ‘Response to different light intensities’).

Of particular sophistication is a study by Mallmann et al. which combined a kinetic model of photosynthesis [31] and a stoichiometric model [32] to explore the evolution of C4 photosynthesis in the genus Flaveria [33]. Here, environment-specificity of the model was achieved by using the output of the kinetic model of C3–C4 intermediate photosynthesis to constrain key C4 parameters, such as net CO₂ uptake, rubisco's v_C/v_O ratio in mesophyll and bundle sheath cells, CO₂ leakage from the bundle sheath, PEP-carboxylase activity in the mesophyll, the activity of NADP-malic enzyme (ME) in the bundle sheath, plasmodesmatal flux of glycine and serine, and decarboxylation by the glycine decarboxylase complex in the genome-scale stoichiometric model of C4 photosynthesis. This way, the authors found that C2 photosynthesis, an intermediate state where CO₂ is relocated from the mesophyll to the bundle sheath cells via the photorespiratory intermediates glycolate and glycerate, caused a N missbalance and that rebalancing N metabolism was a driving force for the evolution of C4 photosynthesis.

In a related study, Blätke and Bräutigam [34] studied the evolutionary trajectory of C4 metabolism by examining selective pressures that potentially have led to the occurrence of this C-fixation mechanism in a two-celled stoichiometric model representing the C4-typical mesophyll-bundle sheath cell Kranz anatomy. In contrast with earlier studies of C4 metabolism [21,32,33] in their model C4 photosynthesis was not enforced by fixing a priori defined flux patterns between the two cell types but the C4 syndrome was allowed to emerge under a set of different input constraints. The authors employed an artifice in which they approximated the CO₂ concentration-dependent changes in rubisco's v_C/v_O ratio by modelling two rubisco populations in the bundle sheath cells — the native rubisco, which performs both carboxylation and oxygenation at a typical C3 plant ratio of 3 : 1 [35], and a carbon-concentration mechanism-dependent rubisco population which only catalyzed the carboxylation of ribulose 1,5-bisphosphate and used exclusively CO₂ released in the bundle sheath cells. Using this setup, they tested several scenarios and found that high photorespiration and N limitation drove the emergence of C4 flux patterns in the model. The model also predicted that light availability and distribution across the leaf section, i.e. between mesophyll and bundle sheath cells could play a role in the evolutionary choices of possible decarboxylation enzymes NAD-ME, NADP-ME, and PEP-carboxykinase. Finally, the model predicted that C2 photosynthesis was optimal under particular conditions and supported the hypothesis of C2 photosynthesis as a stable intermediate state between C3 and C4 photosynthesis.

In our recent work, we focussed on studying the tradeoff between water-saving and leaf productivity and the analysis of alternative CAM-like flux modes in a C3 metabolic network [36]. In contrast with C4 photosynthesis, in CAM photosynthesis initial and re-fixation of CO₂ are not spatially but temporally separated between day and night. This way CO₂ can be taken up through the stomata at colder and more humid night-time hours, thereby minimizing water loss through transpiration. To capture this behaviour, we coupled a biophysical model of gas–water exchange to a time-resolved diel leaf metabolic model. This way, we were able to model the effect of three main abiotic parameters — temperature, relative humidity, and light intensity — on leaf metabolism. We used this environment-coupled model to study the emergence of water-saving flux modes on the Pareto frontier for water-saving and leaf productivity. Our analysis revealed that vacuolar storage capacity is a main determinant of the extent of CAM. Additionally, we identified the mitochondrial enzyme isocitrate dehydrogenase (ICDH) and an isocitrate-citrate-proline-2-oxo-glutarate cycle as a potential contributor to initial carbon fixation at night. Model analysis across a wide range of environmental parameters showed that the water-saving potential of CAM-like flux modes strongly depends on the environment and that under certain conditions CAM with night-time carbon fixation by ICDH could reach 11% total water saving for the conditions tested.

The approaches highlighted so far all have in common that the environment-specificity was achieved by constraining input parameters, such as light intensity, CO₂ or N uptake, the light gradient in the leaf cross section, the gas–water exchange or few reaction fluxes based on experimental data. Most of the fluxes were left unconstrained and could vary between a generic lower and upper boundary. In contrast, several approaches that rely on the integration of Omics data and typically constrain several hundreds to thousands of reaction fluxes have been applied to plant metabolic models and resulted in environment-specific models and flux mode predictions. The majority of the methods employ transcriptomics data to constrain flux boundaries and can be grouped based on whether they employ binarized (e.g. PROM [37], iMAT [37,38], MADE [39], GIMME [40], AdaM [41]) or continuous gene-expression levels (e.g. E-Flux [42]) to constrain fluxes through the respective reactions.

In an earlier study, we used an Arabidopsis metabolic network model accounting for both primary and secondary metabolism [43] and applied a transcriptomics data integrative approach to study changes in leaf metabolism in response to changes in temperature and light intensities [44,45]. The data integration was achieved by using a modified version of the E-Flux approach and enabled us to define three optimization-based indices to characterize different aspects of metabolic pathway behaviour in the context of the entire network. Our study highlighted pathways that showed differential behaviour with respect to a null model and their involvement in the different acclimation responses.

Lakshmanan et al. [46] investigated the effect of four different light treatments and darkness on rice leaf metabolism by combining a rice genome-scale model with transcriptome profiles. The authors performed flux sampling based on the E-Flux method and analyzed the up- and down-regulation of individual metabolic pathways. Additionally, they compared their modelling results to previously measured metabolomics data [47]. Using these combined approaches, they found that photosynthesis and secondary metabolism were up-regulated in blue light and reserve carbohydrates degradation was pronounced in the dark. The analysis further identified phytohormones, such as abscisate, ethylene, gibberellin, and jasmonate as keymarkers of light-mediated regulation and helped elucidate the transcriptional control of red and blue light signals.

Another data-integrative study explored Arabidopsis's response to low and elevated CO₂ by generating condition-specific models based on the integration of transcriptomics data using the iMat approach. Liu et al. showed that the simulated CO₂-fixation at different CO₂ concentrations was consistent with measured CO₂-assimilation curves. Additionally, they predicted post-transcriptionally regulated reactions across different CO₂ concentrations and that low CO₂ stress requires stronger metabolic adjustment than elevated CO₂ conditions.

While the above-outlined approaches incorporated time-resolved and/or compartment-specific data, Bogart and Myers [48] used spatially resolved data to study the transition from source to sink tissue along the leaf gradient. Their analysis predicted metabolic fluxes by integrating spatially resolved transcriptomics and enzyme activity data with a genome-scale model of maize leaf metabolism. Additionally, they modelled the non-linear relationships between CO₂ and O₂ levels and reaction rates in the C4 system. This resulted in a non-linearly constrained model describing mesophyll and bundle sheath cells in 15 segments of the developing maize leaf. The analysis successfully recaptured results from radiolabeling experiments and the observed base-to-tip transition between C-importing and -exporting tissue along the leaf axis. Thus, this approach laid the foundation for studying the response of various heterogeneous metabolic systems to environmental and biochemical perturbations by means of flux-balance analysis.

Nägele and Weckwerth [49] proposed a metabolomics data integrative approach which they applied to study metabolite compartmentation in leaves of different Arabidopsis accessions exposed to low temperature. The authors used a compartmentalized Arabidopsis network reconstruction [43] and compartment-specific metabolomics data to generate a reduced metabolic interaction matrix, i.e. a matrix that contains only measurable metabolites and their respective compartmentation. The extracted subcellular metabolic network comprised over 500 metabolic intermediates and interactions and was integrated with metabolite covariances of different Arabidopsis thaliana accessions. Analysis of these data-integrated networks revealed differences in the regulation of carbohydrate compartmentation in the three investigated accessions and highlighted differential regulation of the interconversion and transport of vacuolar glucose and fructose.

In a followup study to the above-mentioned analysis of Arabidopsis's response to changes in temperature and light intensities [44] we tested the informative value of the transcriptomics data-integrative approach with respect to inferences that could be made on the level of metabolites [50]. We observed that substrates in pathways which our transcript-data integrative study had classified as important for the adjustment processes showed less temporal fluctuations. Moreover, these pathways had on average fewer substrates than the average pathway investigated. This led us to the conclusion that the observed substrate robustness is an inducible genetic mechanism, both depending on the metabolic network structure and the specific environmental condition.

In an attempt to model flux re-routing upon external perturbations we developed a metabolomics data-integrative approach and applied it to predict reactions and pathways with altered fluxes in wild type Arabidopsis leaves and four photorespiratory mutants undergoing high-to-low CO₂ acclimation [51,52]. The approach relied on relative metabolite abundances for at least two metabolic states and employed a mass-action like representation of differential fluxes. We found that the observed flux alterations in the knock-out mutants were mainly attributed to ATP synthase and the photosystem I, and fluxes through reactions from the Calvin–Benson–Bassham cycle, proline biosynthesis, N and redox metabolism and glycolysis. The study highlighted the power of integrated differential flux analysis to complement labor-intensive flux measurements which are usually restricted to small-scale networks.

In recent years, several studies have focussed on modelling plant-microbe interactions by analyzing, both beneficial symbiont relationships, such as interactions between the roots and N fixing soil bacteria [53,54] and plant–pathogen interactions. Two studies have focused on modelling the plant leaf–pathogen interactions between the solanaceous species potato (Solanum tuberosum) and tomato (Solanum lycopersicum) and the late blight causing oomycete Phytophthora infestans.

Botero et al. aimed at gaining a better understanding of the molecular basis of the previously identified decrease in photosynthetic activity of infected potato plants [55,56]. To this end, the authors reconstructed a genome-scale model of potato and integrated gene-expression data from three time-points after pathogen infection [57] using the maximization of biomass as an objective function. They predicted decreased photophosphorylation, the activity of the Calvin–Benson–Bassham cycle and starch synthesis as well as transiently increased photorespiration during the host–pathogen interaction and compared their results to experimental observations from the literature.

In contrast with this host-centric study, Rodenburg et al. studied the nutrient flux from the host to the pathogen during the infection [59]. Thus, the authors coupled a genome-scale model of tomato leaf metabolism [29] and P. infestans [58] by allowing all common metabolites to be exchanged between the plant and the pathogen. Using this combined plant-pathogen model they analyzed four scenarios in which they maximized biomass production of the pathogen and determined (i) the minimal set of nutrients P. infestans needed to import from the tomato leaf, (ii) and (iii) the minimum and maximum number of reactions used by P. infestans, (iv) a combination of (i) and (ii) where they determined the minimal nutrient uptake combined with minimal usage of pathogen reactions. The authors reasoned that the four investigated scenarios represented extreme cases and that the pool of imported nutrients would likely be a combination of the metabolite pools predicted in these scenarios. Additionally, the authors integrated dual-transcriptome time series data of a full late blight infection cycle on tomato leaves. Analysis of these contextualized models indicated that, as the infection progresses, P. infestans performed less de novo synthesis of metabolites and more scavenging from the tomato plant.

• Environment-coupled flux-balance models of leaf metabolism are a versatile means for multi-omics data integration and interpretation [17,60], they have closed gaps in our understanding of plant metabolism and led to the suggestion of engineering strategies to enhance plant's performance [36].

• Yet, one should be aware that flux-balance models are best suited to model plant systems in optimal growth conditions for which an objective function as well as a cellular maintenance costs can be readily defined. For non-optimal growth conditions determining an adequate objective function is less obvious and cellular maintenance costs might vary substantially and could account for a significant portion of the total energy budget [13].

• In the future, we will see continued efforts to integrate the here discussed environment-coupled flux-balance models into multi-scale modelling frameworks which can cover computational models of gene regulation, protein synthesis, metabolic pathways, and plant architecture up to the level of ecosystem models [61–63]. With the advancement and automatization of plant computational modelling I anticipate these models to play a crucial role in future agriculture research by guiding the development of crop species with improved yield and resistance in current and future climates.

The author declares that there are no competing interests associated with this manuscript.

N.T. did not receive specific funding for this work.

CAM

crassulacean acid metabolism

GS

glutamine synthetase

ICDH

isocitrate dehydrogenase

ME

malic enzyme