A Network Biology Approach to Denitrification in Pseudomonas aeruginosa

Pseudomonas aeruginosa is a metabolically flexible member of the Gammaproteobacteria. Under anaerobic conditions and the presence of nitrate, P. aeruginosa can perform (complete) denitrification, a respiratory process of dissimilatory nitrate reduction to nitrogen gas via nitrite (NO 2), nitric oxide (NO) and nitrous oxide (N 2 O). This study focuses on understanding the influence of environmental conditions on bacterial denitrification performance, using a mathematical model of a metabolic network in P. aeruginosa. To our knowledge, this is the first mathematical model of denitrification for this bacterium. Analysis of the long-term behavior of the network under changing concentration levels of oxygen (O 2), nitrate (NO 3), and phosphate (PO 4) suggests that PO 4 concentration strongly affects denitrification performance. The model provides three predictions on denitrification activity of P. aeruginosa under various environmental conditions, and these predictions are either experimentally validated or supported by pertinent biological literature. One motivation for this study is to capture the effect of PO 4 on a denitrification metabolic network of P. aeruginosa in order to shed light on mechanisms for greenhouse gas N 2 O accumulation during seasonal oxygen depletion in aquatic environments such as Lake Erie (Laurentian Great Lakes, USA). Simulating the microbial production of greenhouse gases in anaerobic aquatic systems such as Lake Erie allows a deeper understanding of the contributing environmental effects that will inform studies on, and remediation strategies for, other hypoxic sites worldwide.


Introduction
Denitrification is a facultative anaerobic process in which nitrate is utilized as an alternative terminal electron receptor and dissimilatory nitrate is reduced to nitrogen gas via nitrogen oxides [1][2][3].
Since denitrification is one of the few pathways for producing atmospheric N 2 , it is a major component of the nitrogen cycle [4]. Denitrification occurs in several habitats such as soils, lakes, rivers and oceans [5]. Nitrogen fluxes from marine systems to the atmosphere are between 25 × 10 9 and 179 × 10 9 kilograms per year via microbial denitrification [6]. Pseudomonas aeruginosa, a facultative ubiquitous, and metabolically flexible member of the Gammaproteobacteria, can perform (complete) denitrification under anaerobic conditions and the presence of nitrate. Complete denitrification consists of four sequential steps to reduce nitrate (NO 3 ) to dinitrogen (N 2 ) via nitrite (NO 2 ), nitric oxide (NO), and nitrous oxide (N 2 O), and each step of the pathway is catalyzed by (denitrification) enzymes such as nitrate reductase (nar), nitrite reductase (nir), nitric oxide reductase (nor), and nitrous oxide reductase (nos). The identification and transcriptional control of denitrification genes encoding nar, nir, nor and nos has been largely established. Transcription is dependent on a hierarchy of the FNR-like Crp family transcription factors Anr and Dnr, the two-component system NarXL, and the CbbQ family protein NirQ [7,8], summarized in [4], allowing for experimental validation of N 2 O yield as environmental parameters change.
We have built a combined gene regulatory and metabolic network for the denitrification pathway in Pseudomonas aeruginosa PAO1, a well-studied denitrifier strain (Fig. 1). With this study, we hope to shed light on the environmental factors contributing to greenhouse gas N 2 O accumulation, of particular interest in Lake Erie (Laurentian Great Lakes, USA). Environments such as Lake Erie experience seasonal periods of hypoxic conditions favorable for denitrification, and the endemic microbial community regulates expression of alternative respiratory pathways to adapt to low oxygen (O 2 ) tension.
We are interested in using the model to investigate the effect of PO 4 on the denitrification performance of P. aeruginosa under anaerobic conditions with high NO 3 . Although there are several studies on regulation of denitrification by kinetic mathematical modeling approaches (e.g. [9][10][11][12]), these attempts are not enough to cover the phenomenon at different levels [2]. One of the challenges in building kinetic mathematical models of networks, such as systems of differential equations, is that many of the needed parameters are either not known or unmeasurable. Furthermore, for large networks, kinetic models are difficult to analyze mathematically. Therefore, we take a qualitative approach to model denitrification distinct from the quantitative denitrification models attempted previously. We use a discrete model framework that provides coarse-grained information about the temporal biochemical output of the network in response to environmental conditions. This framework captures attractors (and their biological correspondence, phenotypes) yet it does not render any measurements of time or concentration. In particular, we prefer a time-discrete and multi-state deterministic framework, Polynomial Dynamical System (PDS) [13], to model our denitrification network in Pseudomonas aeruginosa.

Results
The denitrification network consists of molecules, proteins and genes all of which can play an important role in the denitrification process in Pseudomonas aeruginosa. Fig. 1 illustrates a static representation of the variables and their regulations. The blue circular nodes are molecules (O 2 , PO 4 , NO 3 , NO 2 , NO, N 2 O, N 2 ), the yellow rectangular nodes are proteins (PhoRB, PhoPQ, PmrA, Anr, Dnr, NarXL, NirQ) and the pink hexagonal nodes are genes (nar, nir, nor, nos) in the network. The large gray rectangle represents the bacterial cell. The regulatory edges between the nodes are either upregulation/activation (green solid arrows) or downregulation/ inhibition (red dashed arrows). The pathway begins with the phosphate-sensing two component regulatory system PhoRB [14]. PhoRB, the main PO 4 sensor activating the pho regulon, has been recently shown to be a regulator of PhoPQ transcription in the gammaproteobacterium Escherichia coli [15]. In light of the fact that Pseudomonas aeruginosa possesses a similar regulatory system to PhoRB in E. coli [16], it is appropriate to label the PO 4 -sensing regulatory protein as PhoRB in the denitrification network. In this case, the red dashed arrow from PO 4 to PhoRB means that the availability of phosphate, PO 4 , reduces PhoRB function, and the green arrow from PhoRB to PhoPQ means that PhoPQ is activated by PhoRB. Thus, the availability of PO 4 downregulates PhoPQ via PhoRB. The green solid arrow from Anr (anaerobic regulation of arginine deiminase and nitrate reduction) to NarXL and the green solid arrow from NO 3 to the arrow between Anr and NarXL indicate that Anr activates NarXL in the presence of NO 3 . In the same setting, PhoPQ inhibits the expression of PmrA [17]. Low oxygen (O 2 ) tension, which is the major initial signal to turn on the denitrification pathway, can be sensed by Anr [1]. Under anaerobic conditions, Anr primarily promotes Dnr (dissimilatory nitrate respiration regulator) transcription [4]. The effect of Anr on Dnr can be amplified by NarXL [8]. The mechanism of inhibitory effect of PmrA on Dnr [17] is not known, so we assumed that the effect of Anr on Dnr can be reduced by PmrA. The regulatory protein NirQ, which can be activated by NarXL or Dnr, regulates nir and nor coordinately to keep the level of NO low because of toxicity of NO [4]. A NO 3 -responding regulatory protein, NarXL, directly activates nar, and indirectly activates nir and nor via NirQ [4]. The main regulator of the system, Dnr, controls the expression of all denitrification genes (nar, nir, nor, nos) in the presence of NO [18]. Of particular note is the influence of the two-component system PhoPQ on PmrA expression and, subsequently, Dnr expression [17], suggesting that phosphorus (P) availability influences denitrification gene expression (see Fig. 1). This is particularly relevant, since linkages between anaerobic Fe(III) reduction and P release adsorbed to FeOOH in sediments have been recognized for many years [19,20], and recently documented in Lake Erie by stable isotope methods [21].
The actual mechanisms of the relationships in the denitrification network ( Fig. 1) may be quite complex and involve several intermediates. Thus, the network does not represent a biochemical reaction network, for instance, but rather captures the regulatory logic driving the network in a similar way that a circuit diagram explains the function of a circuit board. In the network (Fig. 1), O 2 , PO 4 and NO 3 are external parameters and the remaining nodes are variables. In the discrete setting that is used to model the denitrification network, each node (e.g. an external parameter O 2 or a variable nos) can take up to three states (low, medium, high), and time is implicit and progresses in discrete steps. Our interest lies in perturbation of the external parameters and their effect on the long-term behavior of the variables in the system. S1 Table indicates the discretization values (low/medium/high) for external parameters and nitrogen oxides. Such values incorporate appropriate ranges of long-term nutrient and seasonal oxygen concentrations for Lake Erie [22,23].
The denitrification network is an open system; it exchanges molecules with the outside environment and responds to external stimuli [24]. The molecule NO 3 enters the bacterium and N 2 exits the system once the system is triggered by low O 2 . The model predicts the long-term behavior of the denitrification pathway under various environmental conditions and these predictions are either supported by the literature or validated by experimental results. Fig. 2 indicates the (predicted) attractors of the system under some possible configuration of the external parameters. There are two conditions that we did not focus on. The low NO 3 and low PO 4 condition and the low NO 3 and high PO 4 condition, while possible, are less likely in freshwaters based on a worldwide survey of lakes revealing that N:P stoichiometric ratios average above the ideal Redfield ratio of 16 [25]. Besides, these conditions would be less relevant to current conditions in Lake Erie, for example, as current measurements of nitrate concentrations (averaging 14μM) typically exceed the Km (Michaelis-Menten constant) for nitrate-dependent denitrification in Pseudomonas spp. (for more information, see [26,27]). However, a high P, high NO 3 condition can arise in lakes affected by agricultural nutrient inputs and deposition of P in sediments.
• Prediction 1: If the concentration levels of O 2 and PO 4 are low, and NO 3 is high, then it is a perfect condition for complete denitrification to N 2 . The model suggests that all variables in the network except PmrA are expected to be high and the bacterium reduces NO 3 to N 2 via nitrogen oxides. This prediction is supported by the following studies [1,4,8]. In this condition, Anr senses low O 2 and activates NarXL in the presence of NO 3 [4]. Since the effect of Anr on Dnr is amplified by NarXL but is not reduced by PmrA under low PO 4 conditions, Dnr is highly expressed [8]. The main regulator of the system, Dnr, promotes activation of all denitrification genes (nar, nir, nor, nos), so NO 3 is reduced to N 2 via NO 2 , NO and N 2 O [1].
•   [24]. Phenotypes, biological interpretations of the long-term behavior (steady states), of the system under various environmental conditions can be found in Table 2. Based  on the steady state analysis above, the Pseudomonas network model predicts that elevated PO 4 , hypothesized to increase under hypoxia, acts to modulate the transcriptional network to limit nos gene expression. Thus, the physiological output under this condition will be an increased yield of N 2 O relative to N 2 . Given the prediction that increased PO 4 will influence the N 2 O yield, our experimental results thus far indicate that PO 4 availability modestly, but significantly increases N 2 O yield in this model species (ANOVA p = 0.012; Table 1). While other studies have suggested linkages between N 2 O accumulation and factors such as nosZ vs. nirS/K abundance [29,30], nirS (heme dependent nitrite reductase) genetic diversity [31], or soil pH [32], the data presented here are the first to suggest a role for PO 4 in regulating the denitrification pathway. Given the elevated PO 4 release from FeOOH complexes following sedimentary anaerobic Fe(III) reduction [19,20], hypoxia may yield a high P, high NO 3 condition that enhances N 2 O production.

Discussion
In an aquatic system, oxygen dissolves in water to be available to living aerobic organisms. Hypoxia is the phenomenon of dissolved oxygen below 4mgO 2 per liter. Common reasons for hypoxia include aerobic respiration of decaying algal biomass from bloom events. Such blooms are fueled by increased availability of N and P due to anthropogenic inputs such as agricultural runoff and industrial pollutants [33]. The linkage between high nutrient (N, P) loads and N losses (N 2 and N 2 O) through dissimilatory anaerobic processes was described recently [34]. Hypoxic (low-oxygen) areas, so-called dead zones, often occur in several large bodies of water affected by human activity, including Lake Erie, which is of particular interest. Establishing a better understanding of the nutrient cycling of Lake Erie has quite wide ranging socioeconomic impacts on its recreational area and economy, primarily fisheries. Through denitrification, dead zones lead to microbial production of the greenhouse gas nitrous oxide (N 2 O), which plays a crucial role in ozone layer depletion and climate change. Simulating the microbial production of greenhouse gases in anaerobic aquatic systems such as Lake Erie allows a deeper understanding of the contributing environmental effects that will inform studies on, and remediation strategies for, other hypoxic sites worldwide. During hypoxia, the denitrification rate in Lake Erie is about 150μmolN 2 m −2 h −1 [35]. In addition to oxygen, the intersections of the nitrogen cycle with other geochemical cycles may be important factors influencing denitrification and nitrogen (N) sinks in aquatic systems. In particular, the increased availability of phosphorus (P) has been shown to dictate the rate of nitrogen removal in aquatic systems [34]. Indeed, the transcriptional regulatory network developed for P. aeruginosa indicates that bioavailable phosphate (PO 4 ) is an environmental factor that should be considered. The bacterium Pseudomonas aeruginosa is an example of an abundant microbe in aquatic systems [36], and analysis of Lake Erie metagenomic data sets reveals abundant pseudomonads capable of denitrification (Unpublished data, DOE-JGI). This study describes a computational model of a denitrification network of this bacterium to capture the effect of PO 4 on its denitrification performance in order to shed light on greenhouse gas N 2 O accumulation during oxygen depletion. To our knowledge, this is the first mathematical model of denitrification for this bacterium. Transcription is dependent on a hierarchy of the FNR-like Crp family transcription factors Anr and Dnr, the two-component system NarXL, and the CbbQ family protein NirQ [7,8,37], allowing for experimental measurement of N 2 O as external (environmental) parameters change. The model was constructed based on the pertinent biological literature. Model predictions either agree with current published results or are validated by experimentation. The new biology that our model discovers is that PO 4 availability strongly affects the denitrification activity of P. aeruginosa under anaerobic conditions and the presence of nitrate; high PO 4 can cause less N 2 O reduction to N 2 during denitrification. The data presented here are the first to suggest a role for PO 4 in regulating the denitrification pathway in Pseudomonas aeruginosa.
Current efforts will be expanded to determine how PO 4 affects greenhouse gas N 2 O accumulation during denitrification in P. aeruginosa. According to the model, the activation of Dnr by Anr or the activation of nos in the presence of NO by Dnr can be prevented by high PO 4 . These hypotheses will be tested utilizing quantitative reverse transcriptase PCR (qRT-PCR) to determine Dnr, norB (nitric oxide reductase large subunit gene) and nosZ (encoding nitrous oxide reductase) transcript levels in denitrifying cultures grown in increasing P. Synergistic interactions between individual members of population of Pseudomonas aeruginosa may need to be taken into account and incorporated to the model. For instance, Toyofuku and his colleagues stated that denitrification performance of P. aeruginosa does not only depend upon activation of denitrification genes (nar, nir, nor, nos) but also cell-cell communications under denitrifying conditions [38].
The model described here works well for cultured Pseudomonas, and the next step is to test natural complex microbial communities from different denitrification sites. The effects of PO 4 on N 2 O production will be tested in mesocosms of hypoxic Lake Erie water samples to see if the model described here predicts the community as a whole. By testing the model on environmental samples in mesocosms from Lake Erie and elsewhere, the study can likely be applied broadly to other marine dead zones such as those that routinely occur in the Gulf of Mexico.

Computational Methods
Our network consists of two different sub-networks (metabolic and gene regulatory) and consequently different time scales. From a discrete modeling perspective, this issue can be tackled or ignored only if the long-term behavior of the system is of interest. One could address this issue either (1) using a stochastic framework such as Stochastic Discrete Dynamical System (SDDS) [39] if how fast/slow the reactions are in the network are known/inferred out of a time-course experimental data or (2) introducing time delays by an asynchronous update schedule. Due to inadequate information on the reaction rates, we do not focus on a stochastic framework. Even with a fully asynchronous update schedule, the attractors are preserved for each configuration of external parameters; however, this asynchronous update schedule requires more time steps to reach a steady state than a synchronous update schedule does. Since an asynchronous update schedule provides us more on transient behavior of the system and we are interested in long-term behavior of the system, we prefer to use a deterministic framework with a synchronous update schedule, Polynomial Dynamical System (PDS), which allows us to model regulatory networks over a finite field [13].
Definition 1 Let x 1 , x 2 , . . ., x n be variables which can take values in finite fields X 1 , X 2 , . . ., X n respectively. Let X = X 1 × Á Á Á × X n be the Cartesian product. For each i = 1, 2, . . ., n, we define f i : X ! X i which is an update function that describes the regulation of x i through interaction with other variables in the system. A Polynomial Dynamical System is a collection of n update functions f ¼ ðf 1 ; f 2 ; . . . ; f n Þ : X !X In the model, all external parameters (O 2 , PO 4 , NO 3 ) and some variables (PhoPQ, PmrA, Anr, NarXL) are Boolean (low or high), and other variables are ternary (low, medium or high). There are 14 variables, each of which is labeled for the mathematical formulation. Table 3 indicates the variables, their discretization, update rules and the literature evidence that support these update rules. Inflow substances (i.e. external parameters: O 2 , PO 4 , NO 3 ) in this model give inputs to variables and are involved in the update rules. They do not have update rules because not only they do not have regulators but also we are interested in analyzing the longterm behavior of the model under different configurations of them. The model has only one outflow substance, N 2 , whose regulation depends upon the greenhouse gas N 2 O and its reductase, nos.
Based on the literature, we formulate the regulation of the variables with MIN, MAX and NOT, which correspond to AND, OR and NOT in a Boolean setting. The following are examples for how the update rules are decided: • An update rule of NarXL can be defined as "MIN (Anr, NO 3 )" because NarXL is activated by Anr only in the presence of NO 3 , i.e. both Anr and NO 3 need to be high for NarXL regulation.
• An update rule of nir can be labeled as "MAX (NirQ, MIN(Dnr, NO))" due to the fact that nir is activated by NirQ or Dnr in the presence of NO.
• An update rule of PhoRB can be "NOT (PO 4 )" since PO 4 downregulates PhoRB, i.e. one is low when another is high.
From the update rules in Table 3, for each network variable, we constructed a corresponding transition table, which describes how a specific variable responds to different configurations of their regulators. Although the regulations for most variables can be formulated by MIN, MAX and/or NOT, the regulations of a few variables are very close to some formulation but not quite. For the sake of consistency with biology, we decided to slightly modify the transition table of Dnr, NirQ, nar and NO 2 , whose update rules are marked with an asterix ( Ã ) in Table 3. The transition tables of these variables and more explanation on why the changes were necessary can be found in S2 Table, S3 Table, S4 Table and S5 Table respectively. Table 3. Summary of the model variables, their discretization, update rules and supportive argument. The update rules with an asterix (*) means this update rule is very close to the biological correspondence but not quite. The transition tables of the variables having update rules with an asterix (*) can be found in the Supplementary material.

Index Variables Discretization Update Rules
Literature Evidence Besides, if the variable takes three states (low, medium, high), the current state of the variable is included its own transition table. This does not mean autoregulation/self-regulation; but it is to prevent the variable from jumps between the low (0) state and the high (2) state at the next time step. In other words, including the current state of a ternary variable in its transition table provides a smooth transition among its own states. On the other hand, such jumps cannot occur in a Boolean variable.
After constructing a transition table for each variable x i , an update function can be obtained by interpolating its transition table using the polynomial form: After having all update functions (see S1 Text), we computed the basin of attraction of the whole system under the environmental conditions of interest (see Fig. 2). For model construction and steady state analysis, we used customized Ruby and Perl scripts, which are a part of the source code of Analysis of Dynamic Algebraic Models (ADAM, available at http://adam. plantsimlab.org/), a free of charge web-tool to analyze the dynamics of discrete biological systems [41].

Experimental Methods
Pseudomonas aeruginosa PAO1 cultures were grown in stoppered 20mL serum vials containing glucose minimal medium [42] supplemented with 110mM glucose and 16mM nitrate (NO 3 ). Phosphate (PO 4 ) concentration varied from 1.0mM to 7.5mM, and triplicate culture vials were sampled for headspace gases at 24h and 72h post-inoculation. Gases were dispensed into evacuated exetainers and assayed for nitrous oxide by gas chromatography. Gas production was normalized to cell counts obtained by flow cytometry of culture fluids.
Supporting Information S1