Toward mechanistic modeling and rational engineering of plant respiration

Abstract Plant respiration not only provides energy to support all cellular processes, including biomass production, but also plays a major role in the global carbon cycle. Therefore, modulation of plant respiration can be used to both increase the plant yield and mitigate the effects of global climate change. Mechanistic modeling of plant respiration at sufficient biochemical detail can provide key insights for rational engineering of this process. Yet, despite its importance, plant respiration has attracted considerably less modeling effort in comparison to photosynthesis. In this update review, we highlight the advances made in modeling of plant respiration, emphasizing the gradual but important change from phenomenological to models based on first principles. We also provide a detailed account of the existing resources that can contribute to resolving the challenges in modeling plant respiration. These resources point at tangible improvements in the representation of cellular processes that contribute to CO2 evolution and consideration of kinetic properties of underlying enzymes to facilitate mechanistic modeling. The update review emphasizes the need to couple biochemical models of respiration with models of acclimation and adaptation of respiration for their effective usage in guiding breeding efforts and improving terrestrial biosphere models tailored to future climate scenarios.


Modelling CO2 assimilation and respiration in genome-scale metabolic models
To obtain a physiologically meaningful flux distribution for a photoautotrophic growth, the ratio of fluxes for the oxygenation to carboxylation reactions of RuBisCO ( , ) must be constrained; this is needed to avoid zero flux through photorespiration (except for organisms with carbon concentrating mechanisms, where photorespiration is negligible). Since we attempted to get an overview of photosynthetic and respiratory capabilities based on genome-scale metabolic models of 15 species, we derived a single ratio for and that represents the average over many measurements. To this end, we first obtained values for the specificity of RuBisCO for CO2 over O2, ⁄ , from three different studies (Parry et al., 1989;Zhu et al., 1992;Hermida-Carrera et al., 2016). These values were converted from (unitless) mol/mol to (unitless) bar/bar by using a factor determined from the results of Walker and colleagues (Walker et al., 2013), where ⁄ was given in mol/mol and bar/bar. Finally, the value for the oxygenation to carboxylation ratio, , was calculated using intercellular partial pressures for CO2 ( ) and O2 ( ) given in Farquhar et al. (1980) (Farquhar et al., 1980): In the optimization programs, the ratio between and (correspond to and ) was fixed within two standard deviations ( ) of the distribution of values: − 2 ≤ ≤ + 2 .
The average value for was 0.43 and was 0.06, so the ratio of oxygenation to carboxylation flux was allowed to range between 0.31 and 0.55.
First, we performed flux balance analysis (FBA) with additional constrains on the ratio of and and contraints from loopless FBA (Schellenberger et al., 2011). These additional constraints were applied in three optimization problems, which essentially differed in their objectives. For all simulations, the reversible reactions in the models were split into two irreversible reactions.

∈ {0,1}
∈ ℝ ∈ is the stoichiometric matrix, is the part of that only contains internal reactions i.e., excludes import/export reactions, and denotes the vector of fluxes through all reactions in the model. The binary variable indicated for each internal reaction , whether its direction is forward or backward. The vector contains the Gibbs free energy proxies of the reactions (Schellenberger et al., 2011).
In models that do not contain a single biomass reaction, rather than maximizing , we maximize the sum of fluxes of reactions that export biomass components.

Maximization of CO2 release
To obtain the maximum possible CO2 release, we maximized the flux through an artificial reaction that draws CO2 from the system, while guaranteeing for at least 99% of the optimal objective value determined in problem I.
constraints from problem I where ∑ * is the optimal value of problem I. The variable 2 denotes a reaction that exports CO2 from the system.

Maximization of CO2 assimilation
To estimate the maximum possible CO2 assimilation, all previously existing CO2 uptake or release reactions were blocked and the difference between an artificial CO2 uptake and release reaction was maximized. Again, the sum of biomass reaction(s) was fixed to 99% of the optimal value (problem I). Here, 2 represents a reaction that produces CO2 and thus imports it into the system. All three optimizations were followed by minimization of the first nom of the flux distribution (i.e. the sum of the absolute values of fluxes). The simulations were carried out using the COBRA toolbox function optimizeCbModel with options corresponding to the additional constraints mentioned above. All optimization problems were solved using the Gurobi solver (Gurobi Optimization, 2021) with a feasibility tolerance of 10 −9 .
Calculation of net CO2 assimilation and carbon use efficiency from day models from predicted

flux distributions
The values for net CO2 assimilation ( ) and carbon use efficiency were obtained from the solution to problem I. Net assimilation of CO2 was determined by subtracting the sums of fluxes of all CO2consuming reactions from all CO2-producing reactions. Carbon use efficiency (CUE) was calculated by = 1 − 0.5 + .
In this formula, corresponds to day respiration, which was determined by We note that this definition differs from the classical, given by = 1 − net carbon gain gross carbon assimilation , that involves the net CO2 assimilation, day respiration, and night respiration over a period of time. From the specification, it is expected that over a day, the two formulations result in proportional values.

Calculation of the carbon molar fraction in the biomass reaction
First, the mass fraction of carbon atoms ( , ) for a metabolite (substrate) in the biomass reaction was calculated by: where denotes the molecular weight of the metabolite, , represents the number of carbon atoms in the sum formula of , and is the molecular weight of carbon.
Second, the weight of each metabolite, , was calculated by: where , denotes the stoichiometric coefficient of metabolite in the biomass reaction with unit mmol per gram dry weight (gDW).
Next, the weight of carbon per metabolite weight ( , ) was calculated using the mass fractions of carbon atoms calculated above: Finally, the molar fraction of carbon in the biomass reaction was determined by dividing the sum of carbon weights, divided by the molecular weight of carbon ( [ ]): = 1000 ⋅ , .

Day-and night-specific modelling of respiration
The AraCore model (Arnold and Nikoloski, 2014) was updated by integrating changes in reaction lower and upper bounds (Arnold et al., 2015). Moreover, the ratio between carboxylation flux and oxygenation flux was constrained for the day-specific model as described above in the beginning of the document.
To enable flux through the biomass reaction in the night-specific model, we made the following that differ from the changes in the cited publication. Flux distributions for both models were obtained by solving problem I, followed by minimization of the sum of all fluxes. In the night-specific model, the predicted growth rate was limited to the optimal value for the relative growth rate of the day-specific model. The values for growth respiration ( ) were obtained by summing up all predicted flux values of CO2-producing reactions, except for transport and import reactions ("Tr_*" and "Im_*").

Reaction ID Reaction equation Lower Bound Upper bound Comment
The obtained values for were then scaled to the product of the predicted relative growth rate ( ) and the molar fraction of carbon in the biomass reaction ( ), as described above:

Modelling of single-and double knock-outs and their effect on respiration
To this ends, the AraCore model was used with constraints on the lower and upper bound of the ratio between fluxes through oxygenation and carboxylation reactions. Single and double knock-outs were simulated by blocking the individual reactions of pairs of reactions. Flux distributions were predicted by solving problem I. The predicted relative growth rates and values were then used to calculate the value for ′ , as described above.

Minimization of growth respiration
The AraCore model was used with constraints on the ratio between oxygenation and carboxylation fluxes as described above. All reversible reactions in the model were split into two irreversible reactions.
The following linear fractional program was then solved, which minimizes the ration between fluxes through CO2-producing reactions (  Supplemental Figure S1. Flux of CO2-producing reactions in the AraCore day-and night-specific models. The AraCore model was updated to a day-and night-specific model (Arnold and Nikoloski, 2014;Arnold et al., 2015; see Supplemental Methods for more detailed information). The relative growth rate for both simulations was kept equal to the optimal value predicted using the day-specific model. Flux distributions were predicted using parsimonious flux balance analysis. The flux through CO2-producing reactions is shown which are contained in both the day-and night-specific model, (Fig.  1 for graphical representation). The abbreviations c, m, and h indicate the compartments, where the reactions are located (c: cytosol, m: mitochondrion, h: chloroplast). gDW denotes gram dry weight.
Supplemental Figure S2. Gene expression of respiratory enzymes across different plant tissues. Transcript abundances across 30 tissues were obtained for the enzymes described in Fig. 1 (Mergner et al., 2020 Supplemental Figure S3. Carbon molar fraction in selected models of photosynthetic eukaryotes.
The molar fraction of carbon was calculated as described in the Supplementary Methods. The values were scaled by the sum of metabolite weights that are substrates in the biomass reaction to make the value comparable between models. This value should theoretically be 1 g/gDW but this is not the case for all models. The calculation was performed on all models shown in Fig. 5, but here only we only those models are represented for which the calculation of the carbon molar fraction could be performed. gDW denotes gram dry weight.

Supplemental Figure S4. Changes in nominal and scaled growth respiration in single and double reaction knock-outs.
We predicted flux distributions of all possible single and double knock-outs of reactions in the AraCore model (Arnold and Nikoloski, 2014). The sum of fluxes through CO2-producing reactions was used to approximate nominal growth respiration, . The values for were classified as lower (orange) or higher (blue) compared to nominal of the wild-type model (lower left triangle). The upper right triangle shows the ratio of scaled growth respiration, , between the simulated knock-out and the wild-type model. The values for were obtained by dividing the nominal by the product of predicted relative growth rate and the molar fraction of carbon in the biomass reaction. White color indicates predicted lethal knock-outs. The heatmap only shows the top ten percent of reactions with respect to their average change in ⁄ . The average ⁄ values per reaction across all combinations (n=549) is shown on the right with the same colour coding (i.e. lower (orange) or higher (blue) compared to nominal of the wild-type model). Error bars denote standard deviation. The abbreviations c, m, p, and h indicate the compartments, where the reactions are located (c: cytosol, m: mitochondrion, p: peroxisome, h: chloroplast). The underlying modelling is explained in the Supplemental Methods.  Table S2. Ranking of CO2-producing reactions for the day-and night-specific AraCore adaptations. The column "Flux in pFBA" shows fluxes obtained from solving problem I (see Supplementary Methods) followed by minimization of the sum of fluxes. The last column contains the same fluxes, scaled by the product of the predicted relative growth rate and the carbon molar fraction in the biomass reaction. The abbreviations c, m, p, and h indicate the compartments, where the reactions are located (c: cytosol, m: mitochondrion, p: peroxisome, h: chloroplast). gDW denotes gram dry weight.

Reaction ID Reaction name EC Numbers
Reaction subsystem / pathway