1
SCIENTIFIC REPORTS | (2019) 9:12638 | https://doi.org/10.1038/s41598-019-49079-z
www.nature.com/scientificreports
Modeling the Interplay between
Photosynthesis, CO
2
Fixation, and
the Quinone Pool in a Purple Non-
Sulfur Bacterium
Adil Alsiyabi, Cheryl M. Immethun & Rajib Saha
Rhodopseudomonas palustris CGA009 is a purple non-sulfur bacterium that can x carbon dioxide
(CO
2
) and nitrogen or break down organic compounds for its carbon and nitrogen requirements. Light,
inorganic, and organic compounds can all be used for its source of energy. Excess electrons produced
during its metabolic processes can be exploited to produce hydrogen gas or biodegradable polyesters. A
genome-scale metabolic model of the bacterium was reconstructed to study the interactions between
photosynthesis, CO
2
xation, and the redox state of the quinone pool. A comparison of model-predicted
ux values with available Metabolic Flux Analysis (MFA) uxes yielded predicted errors of 5–19% across
four dierent growth substrates. The model predicted the presence of an unidentied sink responsible
for the oxidation of excess quinols generated by the TCA cycle. Furthermore, light-dependent energy
production was found to be highly dependent on the quinol oxidation rate. Finally, the extent of CO
2
xation was predicted to be dependent on the amount of ATP generated through the electron transport
chain, with excess ATP going toward the energy-demanding Calvin-Benson-Bassham (CBB) pathway.
Based on this analysis, it is hypothesized that the quinone redox state acts as a feed-forward controller
of the CBB pathway, signaling the amount of ATP available.
Purple non-sulfur bacteria (PNSB) are considered to be among the most metabolically versatile groups of bacte-
ria
1,2
. Within this class, Rhodopseudomonas palustris CGA009 (hereaer R. palustris) demonstrates this elasticity
through its ability to survive in a myriad of diverse environmental conditions
3
. It can grow either aerobically or
anaerobically, utilize organic (heterotrophic) or inorganic (autotrophic) carbon sources, and exploit light to obtain
energy when growing anaerobically
3
. Several interesting features have been observed in this bacterium, such as
its consumption of fatty acids, dicarboxylic acids, and aromatic compounds including lignin breakdown products
(LBPs)
4–6
. It is also one of two known bacteria that can express three unique nitrogenases, each with a dierent
transition-metal cofactor
7
. Furthermore, this metabolically versatile strain’s genome encodes the aerobic and anaer-
obic pathways for three of the four known strategies that microbes use to break down aromatic compounds, such
as LBPs
8
. Harnessing R. palustris’ unique metabolic versatilities for the conversion of plant biomass to value-added
products, such as polyhydroxybutyrate (PHB)
9
, n-butanol
10
, and hydrogen
11,12
, has garnered increasing interest.
However, lack of a systems-level understanding of how the bacterium’s complex web of metabolic modules operates
in response to environmental changes is hindering the development of the PNSB as a biochemical chassis.
Several studies conducted on R. palustris showed that in addition to the Calvin-Benson-Bassham (CBB) cycle’s
role of carbon assimilation during autotrophic growth, the pathway plays a major role in maintaining redox bal-
ance under heterotrophic conditions
10,12–14
. It was shown that heterotrophic growth of the PNSB on substrates that
are more reduced than biomass, such as LBPs, is dependent on the availability of an electron sink
13
. CO
2
-xation
using the enzyme ribulose-1,5-biphosphate carboxylase/oxygenase (RuBisCO), nitrogen-xation through the
enzyme nitrogenase
12
, and supplementation with an electron acceptor (e.g., trimethylamine-N-oxide (TMAO))
15
all prevent the inhibitory accumulation of excess reducing agents. erefore, the use of CO
2
as a redox balancing
strategy for the conversion of plant biomass to value-added products is an attractive approach that could increase
protability while improving sustainability. However, the complex interplay between the electrons supplied by
the catabolism of dierent carbon sources, CO
2
xation, and the cyclic electron ow during photosynthesis is not
fully understood, thus diminishing the ability to engineer this promising bacterium.
Department of Chemical and Biomolecular Engineering, University of Nebraska-Lincoln, Lincoln, NE, USA.
Correspondence and requests for materials should be addressed to R.S. (email: rsaha2@unl.edu)
Received: 29 May 2019
Accepted: 19 August 2019
Published: xx xx xxxx
OPEN
2
SCIENTIFIC REPORTS | (2019) 9:12638 | https://doi.org/10.1038/s41598-019-49079-z
www.nature.com/scientificreports
www.nature.com/scientificreports/
A Genome-Scale Metabolic Model (GSMM) provides a mathematical representation of an organism’s met-
abolic functionalities
16,17
by translating the repertoire of biochemical transformations into a stoichiometric
matrix
18
. Due to the underdetermined nature of metabolic networks, optimization tools are used to predict reac-
tion rates for a pre-specied objective function, such as the maximization of biomass
19
. One of the most common
optimization tools used to model metabolism is Flux Balance Analysis (FBA). FBA performs a pseudo-steady
state mass balance for each metabolite in the network to predict the maximum growth rate and corresponding
reaction uxes during the cell’s exponential growth phase
20–24
. Due to the high dimensionality of the network,
other tools such as Flux Variability Analysis (FVA) are used to determine the sensitivity of growth rate as a func-
tion of each reaction ux
25
. Finally, a modied FBA formulation can be used to predict the set of essential genes
under a specied growth condition
26
. us far, a limited number of small-scale metabolic reconstructions have
been developed for PNSB, examining either the central carbon metabolism
27
or the electron transport chain
28
.
However, these models are limited in scope, as they consider less than 4% of the organism’s metabolic function-
ality and are therefore incapable of capturing system-wide interactions between dierent metabolic modules.
Very recently, a GSMM of the bacterium was reconstructed and used to test an array of cellular objectives during
phototrophic growth. Anaerobic growth on acetate, benzoate, and 4-hydroxybenzoate was simulated using eight
dierent biologically relevant objective functions
29
. e model predicted that the organism primarily optimized
for growth, ATP production, and metabolic eciency
29
. However, the model could be improved further by inte-
grating recently annotated metabolic pathways for lignin monomer degradation
30
, as well as making use of exper-
imental data on gene essentiality
31
and metabolic ux analysis for growth under dierent carbon sources
13,14
to
validate and rene the network.
In this work, a GSMM of R. palustris (iRpa940) was constructed to model the bacterium’s metabolic func-
tionality under dierent environmental conditions. e model was used to simulate growth on dierent carbon
sources and showed excellent agreement with experimentally measured uxes
13,14
. Gene essentiality analysis was
also performed for aerobic and anaerobic growth on acetate. e predicted essential genes were compared with
available trans-mutagenesis data
31
, and an accuracy of 84% was achieved. Aer the model indicated the presence
of an unidentied quinol sink, in silico simulations were combined with published in vivo ux measurements
13,14
to study the eect (and the extent) of the quinone redox state on cellular growth, electron transport rate, and CO
2
xation. It was observed that an increase in the quinol oxidation rate resulted in an increase in the electron trans-
port rate, and therefore ATP generation. ese results suggest that redox state acts as a feed-forward controller of
the highly energy-demanding CBB cycle by regulating the rate of light-generated ATP. Overall, an understanding
of the metabolic control points of this interconnected system constitutes the rst step towards engineering strains
capable of more eciently harnessing photosynthetic energy and rerouting this energy towards bio-production
and lignin valorization.
Methods
Model reconstruction. A dra model was rst generated in KBase
32
based on R. palustris’ genome (down-
loaded from the NCBI database on 04/12/2018). KBase uses annotated features in the genome to construct a
list of reactions associated with genes in the organism. Previously published work of the bacterium’s metabolic
network
27
was used to manually curate pathways from the central carbon metabolism and to ensure correct
cofactor usage and gene association. is resulted in an expanded network of high-condence reactions, all asso-
ciated with genes in R. palustris. Experimentally measured concentrations of biomass components are available
for R. palustris when grown on acetate
13
, and were used to develop the biomass equation (see Supplementary
File 1). To minimize the addition of low-condence reactions during gap-lling, the process was broken down
into two steps. First, a subset of high-condence reactions from a recently published genome-scale model of R.
palustris
29
was added to the dra model. Here, high-condence reactions are dened to be the reactions that are
associated with at least onepublished source of annotation. At the end of this step, the majority of the gaps in the
network that precluded the production of biomass existed in partially incomplete linear pathways. erefore, the
ModelSEED database
33
was used to ll the gaps in the network, and a biomass producing model was generated in
KBase
32
. In addition, annotated metabolic pathways for the breakdown of multiple aromatic compounds includ-
ing lignin breakdown products were found in literature
30
and in organism-specic biochemical databases
34,35
,
and were subsequently added to the model (see Fig.S1 in Supplementary File2). Finally, annotated R. palustris
genes were mined from three databases (KEGG
34
, BioCyc
35
, and UniProt
36
) to validate the Gene-Protein-Reaction
(GPR) associations established in the model and to include GPR relationships for reactions added during the
gap-lling process (see Supplementary File 3).
Model simulations. Parsimonious Flux Balance Analysis (pFBA)
37
was used to simulate growth under dif-
ferent environmental conditions. pFBA is analogous to FBA but adds a second objective that minimizes the
sum of all reaction uxes. e two objectives were reformulated into one function through objective tilting
38
as
displayed below.
Maximize vv
subjectto
Sv iI
0 0001
0
(1)
biomass
jJv
j
jJ
ij j
biomass
∑
∑
−.
⋅=∀∈
∈−
∈
≤≤ ∀∈LB vUBjJ
(2)
jj j
3
SCIENTIFIC REPORTS | (2019) 9:12638 | https://doi.org/10.1038/s41598-019-49079-z
www.nature.com/scientificreports
www.nature.com/scientificreports/
where I and J are the sets of metabolites and reactions in the model, respectively. S
ij
is the stoichiometric coe-
cient of metabolite i in reaction j and v
j
is the ux value of reaction j. Parameters LB
j
and UB
j
denote the minimum
and maximum allowable uxes for reaction j, respectively. v
biomass
is the ux of the biomass reaction which mimics
the cellular growth rate.
Model validation. Metabolic Flux Analysis
39,40
(MFA) measurements foranaerobic growth on acetate
13
,
fumarate
14
, succinate
14
, and butyrate
14
were compared with model predicted uxes. Model accuracy for each
growth condition was calculated by taking the sum of percent errors between pFBA-predicted and MFA values
(see Supplementary File 4 for anexample). In addition, R. palustris’ essential genes, determined experimentally
for aerobic growth on acetate
31
, were used to validate the essential genes predicted by the model. Gene essentiality
was predicted in the model by sequentially knocking out each reaction and determining the resulting eect on
the biomass reaction rate
26
. If a reaction knockout resulted in a predicted growth rate that was less than 10% of
the wild type growth rate, the reaction was considered essential
41,42
. Reaction GPRs were then used to map the list
of essential reactions to essential genes. Finally, the list of experimentally determined essential metabolic genes
31
were compared with model predicted essential genes to determine the specicity and sensitivity of the predictions
(see Supplementary File 5).
Results and Discussion
Model Reconstruction and validation. A summary of the iRpa940 model’s major statistics is shown in
Fig.1A. Overall, the 940 genes associated with 1393 model reactions account for 62% of the genes involved in
energy metabolism, biosynthesis, carbon & nitrogen metabolism, and cellular processes in R. palustris’ genome
3
.
Figure1B shows the relative molar abundance of each macromolecular class in R. palustris
13
. is data was used to
calculate the stoichiometric coecients of components in the model’s biomass equation (see Methods). us, an
initial high-condence model containing 540 genes and 915 reactions with no orphan reactions was constructed.
e gap-lling procedure was carried out next in KBase
32
using reactions from the ModelSEED database
33
. Out
of the 478 reactions added during gap-lling, 368 were annotated using information from organism-specic data-
bases (see Methods). A breakdown of the number of GPR relationships established from each annotation source
is shown in Fig.1C. is resulted in the addition of 328 annotated and 110 unannotated (orphan) reactions. e
inclusion of these reactions was necessary to ensure biomass production.
pFBA was used to simulate growth on a number of dierent carbon sources, including carboxylic acids
(acetate, fumarate, succinate and butyrate) and lignin monomers. pFBA is analogous to FBA but adds an outer
objective that minimizes the sum of all reaction uxes (see Methods). is is justied by the assumption that
cells synthesize the minimum amount of cellular machinery required to maintain the maximal growth rate
37
.
Simulating growth using pFBA has two main advantages over FBA. First, pFBA avoids unrealistic ux predic-
tions for reactions participating in thermodynamically infeasible cycles (TICs)
43
. TICs are usually removed from
GSMMs to avoid false predictions; however, when analyzing highly connected networks like that of R. palustris,
Figure 1. Summary of the iRpa940 model statistics and validation. (A) Overall model statistics. (B) Model
biomass component compositions. (C) Sources of gene annotation. (D) Gene essentiality analysis results. G:
Growth (non-essential gene), NG: No Growth (essential gene).
4
SCIENTIFIC REPORTS | (2019) 9:12638 | https://doi.org/10.1038/s41598-019-49079-z
www.nature.com/scientificreports
www.nature.com/scientificreports/
removing these cycles can lead to the model missing certain functionalities and metabolic modes utilized by the
organism. pFBA avoids these false predictions by the additional constraint that reaction uxes should be mini-
mized. Second, the pFBA formulation results in a signicantly reduced set of optimal solutions compared to FBA.
Flux Balance Analysis usually results in a large number of alternate optimal solutions (especially in highly con-
nected networks), most of which are not biologically relevant, and can therefore lead to false conclusions
44
. pFBA’s
additional objective greatly restricts the solution space and leads to more biologically insightful conclusions
37
.
In silico gene essentiality analysis identied 368 essential reactions, out of which 249 were associated with gene
annotations in the model. ese essential reactions were then compared with in vivo gene essentiality data for
aerobic growth on acetate
31
to check the model accuracy (Fig.1D). e calculated sensitivity and false negative
rate (FNR) were consistent with recently published GSMMs
45,46
. Moreover, given that this is a non-model organ-
ism with no well-characterized close relatives, high-condence annotation was not available for the less-studied
pathways. erefore, an automated pipeline like GrowMatch
47
could not be implemented with justiable accuracy
to further improve essentiality predictions.
The eect of the quinone pool on light uptake, carbon dioxide xation, and growth. During
initial phototrophic growth simulations, growth on any of the four carbon sources (acetate, fumarate, succi-
nate, and butyrate) was observed to be hindered due to the accumulation of excess quinols formed in the TCA
cycle. Flux analysis of the electron transport chain (ETC) revealed that the rate of quinol oxidization through the
cytochrome bc1 complex was equivalent to the rate of quinone reduction in the Reaction center (RC). is result
is consistent with previous studies in PNSB
28
, and is necessary for steady-state ow of electrons through the cycli-
cal chain. Furthermore, previous studies on the activity of the ETC concluded that the thermodynamically unfa-
vorable process of reverse electron transfer through NADH dehydrogenase had very low activity compared to the
rate through the RC
28,48
. erefore, this reaction could not account for the oxidation of the excess quinols pro-
duced in the TCA cycle. Since no other high-condence reaction was found to consume quinols in R. palustris,
a quinol “sink reaction” was added to the iRpa940 model. Sink reactions are oen incorporated into metabolic
models when a metabolite is known to be produced during metabolism but for which no means of consumption
have been identied
49
, or to describe the accumulation of a storage compound
49
(e.g. glycogen). Furthermore,
recent experimental work with R. palustris TIE-1 reported the presence of an unidentied quinol-oxidizing reac-
tion that had not been accounted for previously
48
, giving further support to this prediction.
Figure 2. Eect of the Quinol sink rate on: (A) Light uptake rate, (B) Growth rate, (C) Carbon source uptake
rate, and (D) Carbon xation rate for growth on four carbon sources. ace: acetate, but: butyrate, suc: succinate,
fum: fumarate.
Carbon
source
QH
2
oxidation rate
(mmol/gDW/hr)
Q reduction rate
a
(mmol/gDW/hr)
Electron transport
rate (dmol/gDW/hr)
CO
2
xation rate
(mmol/gDW/hr)
% CO
2
xed
b
Net CO
2
excretion
rate (mmol/gDW/hr)
Acetate 52.5 39.1 5.3 29.7 73.2 10.9
Butyrate 54.9 37.4 5.4 57.8 — −18.6
c
Succinate 49.8 49.2 3.0 35.6 50.6 34.8
Fumarate 0 0 2.3 17.3 25.1 51.5
Table 1. Predicted reaction rates for growth on four dierent carbon sources.
a
e rate of quinone reduction
in the TCA cycle.
b
e rate of CO
2
xation divided by the rate of total CO
2
produced.
c
CO
2
was supplied in the
media during growth on butyrate.
5
SCIENTIFIC REPORTS | (2019) 9:12638 | https://doi.org/10.1038/s41598-019-49079-z
www.nature.com/scientificreports
www.nature.com/scientificreports/
pFBA simulations were conducted under dierent quinol sink rates to qualitatively predict how changes in
the quinone redox state aected the rest of the metabolic network. e quinol sink reaction was treated as a
parameter in the model and pFBA simulations were conducted at varying quinol oxidation (sink) rates to deter-
mine how light uptake (i.e. Electron Transport Rate or ETR), growth, and CO
2
xation are aected by changes
in the quinone redox state (Fig.2). Carbon uptake was restricted to a maximum value of 100 mmol/gDW/hr
for acetate and 50 mmol/gDW/hr for fumarate, succinate, and butyrate to ensure the same number of carbons
were being taken up. MFA values were scaled to the same carbon uptake rates
13,14
. For growth on butyrate, the
supplementation of CO
2
is required for growth, as the substrate is more reduced than biomass and requires an
electron sink
14
. e media was supplied with CO
2
at a maximum uptake rate of 32.1 mmol/gDW/hr to match
MFA observations. Since steady-state GSMMs cannot capture metabolite concentrations, the redox state cannot
be quantied directly. Instead, the qualitative behavior of the redox state was predicted by varying the rate of the
quinol sink. As the quinol oxidation rate increases, the quinone pool becomes more oxidized. Using experimental
MFA data
13,14
, the quinol oxidation rate was predicted for each of the four substrates (Table1). ese values were
calculated by minimizing the sum of errors between the in silico-generated pFBA uxes and the MFA ux values.
e table also shows the quinone reduction rate through the TCA cycle for each carbon source. e percentage
of CO
2
xed was dened as the rate of CO
2
xation divided by the total rate of CO
2
produced metabolically.
Figure3 shows the resulting ux predictions obtained at the predicted quinol oxidation rates for growth on
acetate (Fig.3A), and the calculated percent errors of these predictions for each carbon substrate (Fig.3B). A
Figure 3. Comparison of model-predicted vs MFA-generated ux values for reactions involved in central
carbon metabolism. (A) Metabolic ux map showing reaction rates for growth on acetate (B) Percentage error
between model predictions and MFA ux values for growth on four carbon sources.