Transient Accumulation of NO2- and N2O during Denitrification Explained by Assuming Cell Diversification by Stochastic Transcription of Denitrification Genes.

Denitrifying bacteria accumulate [Formula: see text], NO, and N2O, the amounts depending on transcriptional regulation of core denitrification genes in response to O2-limiting conditions. The genes include nar, nir, nor and nosZ, encoding [Formula: see text]-, [Formula: see text]-, NO- and N2O reductase, respectively. We previously constructed a dynamic model to simulate growth and respiration in batch cultures of Paracoccus denitrificans. The observed denitrification kinetics were adequately simulated by assuming a stochastic initiation of nir-transcription in each cell with an extremely low probability (0.5% h-1), leading to product- and substrate-induced transcription of nir and nor, respectively, via NO. Thus, the model predicted cell diversification: after O2 depletion, only a small fraction was able to grow by reducing [Formula: see text]. Here we have extended the model to simulate batch cultivation with [Formula: see text], i.e., [Formula: see text], NO, N2O, and N2 kinetics, measured in a novel experiment including frequent measurements of [Formula: see text]. Pa. denitrificans reduced practically all [Formula: see text] to [Formula: see text] before initiating gas production. The [Formula: see text] production is adequately simulated by assuming stochastic nar-transcription, as that for nirS, but with a higher probability (0.035 h-1) and initiating at a higher O2 concentration. Our model assumes that all cells express nosZ, thus predicting that a majority of cells have only N2O-reductase (A), while a minority (B) has [Formula: see text]-, NO- and N2O-reductase. Population B has a higher cell-specific respiration rate than A because the latter can only use N2O produced by B. Thus, the ratio [Formula: see text] is low immediately after O2 depletion, but increases throughout the anoxic phase because B grows faster than A. As a result, the model predicts initially low but gradually increasing N2O concentration throughout the anoxic phase, as observed. The modelled cell diversification neatly explains the observed denitrification kinetics and transient intermediate accumulations. The result has major implications for understanding the relationship between genotype and phenotype in denitrification research.

Introduction

The dissimilative reduction of nitrate () to nitrite (), nitric oxide (NO), nitrous oxide (N2O), and finally to N2 (denitrification) is an indispensable process in the nitrogen cycle, returning N to the atmosphere as N2. However, denitrification significantly leaks the gaseous intermediates NO and N2O, both with serious consequences for the environment. N2O catalyses depletion of the stratospheric ozone [1] and causes global warming, contributing ~10% to the anthropogenic climate forcing [2]. Data suggests that since the 1950s, the atmospheric N2O has been increasing, and before being photolysed in the stratosphere, the gas persists for an average ~120 years in the troposphere [3]. ~70% of global N2O emissions are tentatively attributed to microbial nitrification and denitrification in soils [4], where denitrification, generally, is considered a more dominant source [5].

To mitigate N2O emissions, we need to understand the physiology of denitrifiers

To devise robust strategies for mitigating global N2O emissions, a good understanding of its primary source is imperative, i.e., genetics, physiology, and regulatory biology of denitrifiers. Any knowledge of the environmental controllers of N2O is incomplete without understanding the causal relationships of such controllers at the physiological level [6].

The biogeochemical models developed for understanding the ecosystem controls of denitrification and N2O emissions treat the denitrifying community of soils and sediments as a single homogenous unit with certain characteristic responses to O2 and concentrations [6,7]. Natural denitrifying communities, however, are mixtures of organisms with widely different denitrification regulatory phenotypes [8]. The regulatory response of such mixtures is not necessarily equal to the ‘sum of its components’ because there will be interactions, not the least, via the intermediates NO and . Hence, it is probably a mission impossible to predict the regulatory responses of complex communities based on their phenotypic composition. Nevertheless, investigations of the regulation in model organisms like Pa. denitrificans provide us with essential concepts, enhancing our ability to understand the regulatory responses of mixed communities and to generate meaningful hypotheses. Thus, future biogeochemical models of N2O and NO emissions are expected to have more explicit simulations of the regulatory networks involved, and a first attempt has recently been published [9].

Simulating the cell diversification in response to impending anoxia to analyse its implications for , N2, and N2O kinetics

Dynamic modelling has been used to a limited extent to analyse various denitrification phenotypes; for example, to analyse and reduction and gas-kinetic data for individual strains [10] and mixtures of selected phenotypes [11]; to model the consequence of competition for electrons between denitrification reductases [12,13]; to investigate the control of O2 on denitrification enzymes and inhibition of cytochrome c oxidase by NO in Agrobacterium tumefaciens [14]; and to examine the effect of copper availability on N2O reduction in Paracoccus denitrificans [15]. In our previous model [16], we simulated O2 and N2 kinetics from batch incubations of Pa. denitrificans [8,17] to test if a postulated cell diversification, driven by stochastic initiation of nirS, could explain the N2 production kinetics in -supplemented media. The available data also contained -supplemented treatments but and were not monitored, and the experiment provided no information about the N2O kinetics, except that the concentrations were extremely low (below the detection limit of the thermal conductivity detector used). Recently, a neat dataset was generated from batch incubations supplemented with , with frequent measurements of and a more sensitive detection of N2O by an electron capture detector [18]. That encouraged us to extend our previous model and simulate the cell diversification during transition from oxic to anoxic conditions, targeting the regulation of Nar and cNor/NosZ (N2O emissions) in Pa. denitrificans.

Regulatory network of denitrification in Paracoccus denitrificans

Pa. denitrificans is a facultative anaerobe capable of reducing all the way to N2:

In response to impending anoxic conditions, the organism sustains respiratory metabolism by producing the membrane-bound cytoplasmic nitrate reductase (Nar), cytochrome cd nitrite reductase (NirS), cytochrome c dependent nitric oxide reductase (cNor), and nitrous oxide reductase (NosZ). Transcription of the genes encoding these reductases (narG, nirS, norBC, and nosZ, respectively) are regulated by the FNR-type proteins FnrP, NarR, and NNR. FnrP contains a 4Fe-4S cluster for sensing O2, and NNR harbours a NO-sensing haem; NarR, however, is poorly characterised and is most likely a -sensor [19-21]. All these sensors remain inactive during aerobic growth conditions [19].

Transcription of denitrification genes in Pa. denitrificans

FnrP and NarR facilitate a product-induced transcription of the nar genes, and NNR facilitates a product-induced transcription of the nirS genes (Fig 1, see P1 and P2): Low oxygen concentration ([O2]) activates the self-regulating FnrP, which induces nar transcription in coaction with NarR. The self-regulating NarR was previously assumed to be activated by either or [21], but a recent proteomics study indicates that is the activator [19]. Thus once a cell starts producing traces of , nar expression becomes autocatalytic. Transcription of nirS is induced by NNR, which is apparently inactivated by O2 [22,23], but under anoxic/micro-oxic conditions, NNR is activated by NO. Thus, once traces of NO are produced, the expression of nirS also becomes autocatalytic [19,20]. In contrast, nor transcription is substrate (NO) induced via NNR while nosZ is equally induced by NNR or FnrP [24]. High concentrations of NO may constrain nar transcription by inactivating FnrP [25] and, like O2, render NosZ dysfunctional by inactivating the CuZ subunit of the reductase [26], but these observations are ignored in our model because Pa. denitrificans restricts [NO] to very low levels.

Fig 1

Regulatory network of denitrification in Pa. denitrificans.

The network is driven by four core enzyme-complexes: Nar (transmembrane nitrate reductase encoded by the narG gene), NirS (cytochrome cd nitrite reductase encoded by nirS), cNor (NO reductase encoded by norBC), and NosZ (N2O reductase encoded by nosZ). When anoxia is imminent, the low [O2] is sensed by FnrP, which in some interplay with NarR induces nar transcription. NarR is activated by ; thus once a cell starts producing traces of , nar expression becomes autocatalytic (see P1). Transcription of nirS is induced by NNR, which is activated under anoxic/micro-oxic conditions by NO; thus once traces of NO are produced, the expression of nirS also becomes autocatalytic (see P2) [20]. The activated P2 will also induce nor and nosZ transcription via NNR. The transcription of nosZ, however, can also be induced equally and independently by FnrP [24]. Micromolar concentrations of NO may inactivate both FnrP [25] and NosZ [26]. These observations, however, are ignored for our modelling because Pa. denitrificans restricts NO to nanomolar levels.

Entrapment of cells in anoxia: The underlying hypothesis and modelling

Denitrification proteome, once produced in response to an anoxic spell, is likely to linger within the cells under subsequent oxic conditions, ready to be used if anoxia recurs. But the proteome will be diluted by aerobic growth because the transcription of denitrification genes is inactivated under oxic conditions [20]. Hence, a population growing through many generations under fully oxic conditions is expected to undertake de novo synthesis of denitrification enzymes when confronted with anoxia. Batch cultivations of such aerobically raised Pa. denitrificans provided indirect evidence for a novel claim that, in response to anoxia, only a small fraction of the incubated population is able to produce denitrification proteome [8,17,27,28]. Our dynamic modelling of Bergaust et al.’s [17] -supplemented incubations corroborated this, suggesting that a probabilistic function (specific probability = 0.005 h-1) resulting in the recruitment of 3.8–16.1% of all cells to denitrification is adequate to explain the measured N2 kinetics [16].

Our model was based on the hypothesis that the entrapment of a large fraction in anoxia is due to a low probability of initiating nirS transcription, which in response to O2 depletion is possibly mediated through a minute pool of intact NNR, crosstalk with other factors (such as FnrP), unspecific reduction of to NO by Nar, and/or through non-biologically formed traces of NO found in a -supplemented medium. Regardless of the exact mechanism(s), once nirS transcription is initiated, the positive feedback via NO/NNR (Fig 1, see P2) would allow the product of a single transcript of nirS to induce a subsequent burst of nirS transcription. The activated positive feedback will also help induce nor and nosZ transcription via NNR, rapidly transforming a cell into a full-fledged denitrifier. We further hypothesised that recruitment to denitrification will only be possible as long as a minimum of O2 is available because, since Pa. denitrificans is non-fermentative, the synthesis of first molecules of NirS will depend on energy from aerobic respiration.

The above hypothesis was modelled by segregating the culture into two pools (subpopulations): one for the cells without (ND−) and the other with denitrification enzymes (ND+). Initially, all cells were ND−, growing by consuming O2. As [O2] fell below a certain threshold, ND− recruited to ND+ with a constant probability (h-1), assumed to be that of the nirS transcriptional activation, and the recruitment halted as O2 was completely exhausted, assuming lack of energy (ATP) for enzyme synthesis.

Underlying assumptions and aims of the present modelling

The present model is an extension of that developed in Hassan et al. [16]. Here we have divided the respiring culture into four pools (Fig 2A):

Fig 2

A stock and flow diagram illustrating the model’s structure.

A. Cell diversification and growth; B. O2 kinetics; C. Denitrification kinetics. The squares represent state variables, the circles the rate of change of the state variables, the edges (thicker arrows) depict flows into or out of the state variables, the shaded ovals auxiliary variables, and the arrows portray mutual dependencies between the variables. All feedback relationships among the three model sectors could not be shown; however, for illustration the feedback relationships of one sub-population (Z−) are shown (dashed arrows). Within each square (state variable), t0 refers to the initial value.

Z−: cells without Nar, NirS, and cNor

ZNa: cells with Nar

ZNaNi: cells with Nar, NirS, and cNor

ZNi: cells with NirS and cNor

All these subpopulations are assumed to scavenge O2 (if present) and produce NosZ in response to impending anoxia. The latter because the nosZ genes are readily induced by the O2-sensor FnrP [24].

The Z− pool (Fig 2A) contains the inoculum that grows by aerobic respiration. As [O2] falls below a critical threshold [empirically determined, 18], the cells within Z− are assumed to start synthesising Nar with a certain probability and populate the ZNa pool. The aim here is to investigate whether, like for nirS, the initiation of nar transcription (by some combined activity of FnrP and NarR) can also be explained as a probabilistic phenomenon, quickly differentiating a cell into a full-fledge scavenger through product () induced transcription via NarR (Fig 1, see P1). If so, we were interested to estimate what fraction of the cells is required to adequately simulate the measured data ( production), aiming at scrutinising the general assumption that all cells in batch cultures produce Nar in response to impending anoxia.

Next, when [O2] is further depleted to another critical threshold [18], the Z− and ZNa cells are assumed to initiate nirS transcription with a low per hour probability and, thereby, populate the ZNi and ZNaNi pools, respectively. As explained above for our previous model, NirS + cNor production is assumed to be a) coordinated because the transcription of both nirS and nor is induced by NO via the NO-sensor NNR (Fig 1), and b) stochastic because the initial transcription of nirS (paving the way for the autocatalytic expression of NirS and substrate-induced nor transcription) happens in the absence of NO or at too low [NO] to be sensed by NNR.

Synthesis of denitrification enzymes requires energy, which all the subpopulations can obtain by respiration only. Hence, the initiation of the autocatalytic expression of nar and nirS (i.e., recruitment to ZNa and ZNaNi/ZNi, respectively, Fig 2A) depends on the availability of the relevant terminal e--acceptor(s) above a critical concentration to sustain a minimum of respiration. For Z−, the only relevant e--acceptors are O2 and the traces of N2O produced by ZNi and ZNaNi. The same applies For ZNa, but in addition, this subpopulation can also obtain energy by reducing , if present. In our previous model [16], we assumed that recruitment to denitrification was sustained by energy from O2-respiration only; not because we simulated -supplemented treatments, and not by N2O because we naively assumed that the pool of this e--acceptor was insignificant (N2O concentrations were below the detection limit of the system used for those experiments). However, the present model assumes that the recruitment from Z− to ZNa and Z− to ZNi is sustained by both O2- and N2O-reduction, and the recruitment from ZNa to ZNaNi is sustained by O2-, N2O- and -reduction, when above a critical minimum (). The default value for was set to an arbitrary low value (= 0.44% of the maximum e--flow rate to O2), and we have investigated the consequences of increasing, decreasing, and setting = 0.

The expressions of nar and nirS + nor (recruitments to ZNa and ZNaNi/ZNi, respectively, Fig 2A) are modelled as instantaneous discrete-events in each cell, thus ignoring the time-lag from the initiation of gene transcription till the cell is fully equipped with the reductase(s) in question. That is because the lag observed between the emergence of denitrification gene transcripts and the subsequent gas products suggests that the synthesis of denitrification enzymes takes less than half an hour [17,18], which is negligible for our purposes here.

The main purpose of the present modelling is to investigate if a full-fledged model, including all four functional denitrification reductases, could adequately simulate the observed kinetics and stoichiometry of denitrification products [18]. These cultures reduced all available to prior to the onset of gas production and accumulated traces of N2O throughout the anoxic phase, as illustrated in S1 Fig In particular, we were interested to investigate the kinetics, controlled by nar- and nirS transcription, and to test if the peculiar N2O kinetics (low, but increasing concentrations throughout the anoxic phase) could be explained by our modelled cell diversification.

Materials and Methods

An overview of the modelled experiment

Batch incubation

Qu [18] incubated Pa. denitrificans (DSM-413) at 20°C using 50 mL Sistrom’s [29] medium in 120 mL gas-tight vials. Either succinate or butyrate (5 mM) was used as the main carbon source, enough to secure consumption of all available e--acceptors. After distribution of the medium, each vial was loaded with a magnetic stirring bar, sterilised through autoclaving, supplemented with 2 mM KNO3, and was tightly sealed. To remove O2 and N2 from the headspace, the headspace air was evacuated and replaced by helium (He) through several cycles of evacuation and He-filling (He-washing). Some vials were supplemented with oxygen to reach 7 vol.% O2 in headspace (treatment designated 7% O2). The remaining vials received no O2 (designated 0% O2, although there were traces of O2 present despite the He washing). For each treatment (i.e., C source and initial O2), there were three replicates, and each vial was inoculated with 2.2×108 aerobically grown cells.

and gas measurement

Gases (CO2, O2, NO, N2O, and N2) were monitored by frequent sampling of the headspace, using an improved version of the robotised incubation system [30]. In short, the system draws gas samples from the headspace (peristaltic pumping) via the septum pierced by a needle, filling three loops used for injecting samples into the two GC columns and the chemiluminescence NO analyser. The sample drawn is replaced by He (reversing the peristaltic pump), thus securing ~1 atm pressure. The primary improvements of the new system are a more sensitive detection of N2O (by an electron capture detector), lower sampling volumes (~1 mL), and lower leaks of O2 and N2 through the sampling system (4 nmol O2 and 12 nmol N2 per sampling, which is ~20% of that for the old system).

To extract samples for measuring without tampering the original vials, identical (parallel) vials were prepared for each treatment. Using sterile syringes, samples of 0.1 mL were regularly drawn from the liquid-phase of the parallel vials and immediately analysed for .

Results for one of the treatments are shown in S1 Fig, illustrating the complete reduction of to prior to the onset of significant N-gas production. In previous experiments [17], N2O concentrations were below the detection limit of the system, but thanks to the new system, the N2O kinetics were monitored with a reasonable precision.

The model

The model is constructed in Vensim DSS 6.2 Double Precision (Ventana Systems, inc. http://vensim.com/) using techniques from the field of system dynamics [31].

Cell diversification and growth

The respiring population is divided into four subpopulations, according to their reductases (Fig 2A): 1) Z−: cells without Nar, NirS, and cNor; 2) ZNa: cells with Nar; 3) ZNaNi: cells with Nar, NirS, and cNor; and 4) ZNi: cells with NirS and cNor. All the subpopulations are assumed to equally respire O2, if present, and express nosZ in response to oxygen depletion [24]. Z− contains the inoculum (= 2.2×108 cells) that grows by aerobic respiration. As O2 is depleted, the Z− cells populate the other pools by producing Nar and/or NirS + cNor.

The recruitment from Z− to ZNa (RNa, Fig 2A) takes place first:

(cells h-1)

where r (O ,N O) is a conditional specific probability (h-1) for any Z− cell to initiate nar transcription (quickly transforming it into a scavenger through autocatalytic gene expression, see Fig 1, P1):

(h-1)

where rNa (h-1) is a constant specific probability for a cell to initiate nar transcription once O2 concentration in the aqueous-phase ([O2]aq, mol L-1) falls below a critical concentration ([O2]), empirically determined as the [O2]aq (= 4.75×10−5 mol L-1) at the outset of accumulation in the medium [18]. The second condition for a cell to produce first molecules of Nar is a minimum of e--flow to an e--acceptor (, mol e- cell-1 h-1), assumed to generate minimum ATP required for protein synthesis. and (mol e- cell-1 h-1) are the cell-specific velocities of e--flow to O2 and N2O, respectively. The latter is weighed down by 0.5 because mole ATP per mole e- transferred to /NOx is lower for denitrification than for aerobic respiration [17,20]. For a Z− cell, and are not considered here, since such a cell is assumed to have no NirS and cNor.

The fraction of the cells that successfully produces Nar (FNa) is calculated based on the integral of the recruitment (Eq 1):

(dimensionless)

where t is the time-window available for the recruitment. In theory, t is the time-period when (Eq 2). Since the e--flow to N2O started after all had been reduced to (S1 Fig), the recruitment based on would be inconsequential for the simulated (and measured) kinetics. Therefore, to calculate the functional FNa actually responsible for producing , we ignored the N2O-sustained recruitment, thus considering t to be the time when .

Next, the cells within ZNa and Z− are recruited to ZNaNi and ZNi (RNaNi and RNi, respectively, Fig 2A), as they are assumed to stochastically initiate nirS transcription, paving the way for NO/NNR mediated autocatalytic expression of nirS + nor (Fig 1). In principle, the rates of both these recruitments are modelled as that of the recruitment from Z− to ZNa (Eqs 1 and 2): a) Both trigger as O2 falls below another critical concentration ([O2]), low enough to activate NNR to induce nirS transcription; [O2] (= 1.16×10−5 mol L-1) is empirically determined as the O2 concentration at the outset of NO accumulation [18]. b) Both continue as long as a minimum of e--flow to the relevant terminal e--acceptor is possible, sustaining the respiratory metabolism to generate ATP for protein synthesis:

(cells h-1)

(h-1)

where rNi is a constant specific probability (h-1) for the initiation of nirS transcription. and are multiplied with 0.5 for the same reasons as described for Eq 2.

The recruitment from Z− to ZNi (RNi, Fig 2A) is modelled as a product of Z− and a conditional specific probability, r (O ,N O), which is different from Eq 5 only in that is omitted, since Z− do not possess Nar:

(cells h-1)

(h-1)

The fraction that successfully produced NirS + cNor (FNi) is calculated based on the integral of RNaNi and RNi:

(dimensionless)

where t is the duration of the recruitment from ZNa to ZNaNi, i.e., when (Eqs 4 and 5), FNa is the fraction recruited to the pool of Nar positive cells (ZNa, Eq 3), and t is the duration of the recruitment from Z− to ZNi, i.e., when (Eqs 6 and 7).

Each of the populations will grow depending on the rates of e--flow to the various e--acceptors they are able to use:

(cells h-1)

where (cells mol-1 e- to X = O2 or /NOx) is the growth yield determined under the actual experimental conditions, and (mol e- cell-1 h-1) is the cell-specific velocity of e--flow to X (O2 or /NOx), which depends on the concentration of the e--acceptor (see Eqs 17, 20 and 28). For and , a restricted velocity () is used so that when electrons flow to O2, , and simultaneously, the total ve− per cell does not exceed the maximum electrons that the TCA cycle () can deliver per hour (see Eqs 21 and 22).

O2 kinetics

O2 is initially present in the headspace (, mol, initialised according to the experiment, see Table 1) but is transported to the liquid-phase () due to its consumption therein (Fig 2B). The transport rate () is modelled according to Molstad et al. [30]:

Table 1

Simulated experiment [18].

Batch	C-source	O_2__g(t_0) (vol.%)*	NO3−(t_0) (mM)	Replicates
1	Butyrate	~0	2	3
2	Butyrate	7	2	3
3	Succinate	~0	2	3
4	Succinate	7	2	3

*Target values for initial O2 concentrations in the headspace (vol.%). ~0 means that the intended concentration should be zero, but there were detectable traces of O2, despite several cycles of evacuation and He-flushing of the headspace.

(mol h-1)

where kt (L h-1) is the empirically determined coefficient for the transport of gas between the headspace and the liquid, (mol L-1 atm-1) is the solubility of O2 in water at 20°C, (= [O2]g × R × T, atm) is the partial pressure of O2 in the headspace, and [O2]aq (mol L-1) is the O2 concentration in the liquid .

In addition, the model simulates the changes in due to sampling. The robotised incubation system used monitors gas concentrations by sampling the headspace, where each sampling alters the concentrations in a predictable manner: a fraction of is removed and replaced by He (dilution), but the sampling also results in a marginal leakage of O2 through the tubing and membranes in the injection system. The net change in (ΔO2(S)) as a result of each sampling is calculated as:

(mol h-1)

where (mol vial-1) is the O2-leakage into the headspace, D (dilution) is the fraction of each headspace gas removed and replaced by equal amount of He, and ts (h) is the time taken to complete each sampling. ΔO2(S) is negative if is high and marginally positive at very low oxygen concentrations.

O2 in the liquid-phase (, mol, Fig 2B) is initialised by assuming equilibrium with at the time of inoculation (). The dynamics of are modelled as a function of transport between the headspace and the liquid (, Eq 13) and its reduction rate (, mol h-1):

(mol h-1)

where Z−, ZNa, ZNaNi, and ZNi (cells) are all the sub-populations present (described above); thus, we assume that all cells have the same potential to consume O2. (mol cell-1 h-1) is the cell-specific velocity of O2 consumption, obtained by the velocity of e--flow to O2 , where is modelled as a Michaelis-Menten function of oxygen concentration:

(mol e- cell-1 h-1)

where (mol e- cell-1 h-1) is the maximum velocity of e--flow to O2 per cell (determined under the actual experimental conditions), [O2]aq (mol L-1) is the O2 concentration in the liquid-phase, and (mol L-1) is the half-saturation constant for O2 reduction.

Denitrification kinetics

The and pools (mol, Fig 2C) are initialised according to the experiment (Table 1; = 0). The kinetics of these nitrogen oxyanions () are modelled as:

(mol h-1)

where (mol h-1) is the reduction rate, ZNa + ZNaNi (cells) is the total number of cells with Nar, ZNaNi + ZNi (cells) is the total NirS active population, and (mol cell-1 h-1) is the cell-specific velocity of consumption, obtained by the velocity of e--flow to . The latter is modelled as a Michaelis-Menten function of concentration:

(mol e- cell-1 h-1)

where (mol e- cell-1 h-1) is the maximum velocity of e--flow to per cell (determined under the actual experimental conditions), (mol L-1) is the concentration in the aqueous-phase, and (mol L-1) is the half-saturation constant for reduction.

The velocity of and consumption had to be restricted () to ensure that when electrons flow to O2, and simultaneously, the total ve− per cell does not exceed an estimated maximum delivery of electrons from the TCA cycle (). In competition for electrons, O2 is prioritised [20], followed by and , respectively [18]:

(mol e- cell-1 h-1)

where is the realised e--flow to , limited either by available or the availability of electrons (due to competition with O2); is the realised e--flow to . Such competition for electrons was not implemented for and because at the onset of NO-, N2O- and N2 production, the total velocity of e--flow to all available e--acceptors (as predicted by the enzyme kinetics alone) never exceeded .

Gas consumption and production takes place in the aqueous phase, but the gases are transported between aqua and the headspace depending on their concentrations in the two phases. Each gas in aqua, Xaq (molN, Fig 2C), is modelled as a function of production, consumption (not applicable to N2), and the net transport, where N2Oaq and are initialised with zero, and NOaq is initialised with a negligible 1×10−25 mol to avoid division by zero (in Eq 28).

(molN h-1)

where (molN h-1) is the relevant /NOx reduction rate, and represents the gas transport rate between aqua and the headspace (Eq 29; N.B. < 0 for the net transport from aqua to the headspace).

The reduction of NO to N2O (RrNO) and N2O to N2 () is modelled likewise as a function of the number of relevant cells and the velocity of e--flow to NO and N2O (mol e- cell-1 h-1), respectively:

(molN h-1)

where v NO and are obtained by the velocity of e--flow to NO and N2O, respectively . is modelled as a Michaelis-Menten function of [N2O]aq, similarly as that of O2, , and (Eqs 17 and 20), but is modelled assuming a cooperative binding of two NO molecules with cNor to form N2O [32]:

(mol cell-1 h-1)

where (mol e- cell-1 h-1) is the empirically determined maximum velocity of e--flow to NO per cell, [NO]aq (mol L-1) is the NO concentration in the liquid-phase, and K1NO & K2NO (mol L-1) are the equilibrium dissociation constants for the cNor/NO- and cNor/(NO)2 complex, respectively.

The transport of NO, N2O, and N2between the liquid and the headspace (Eqs 23–25) is modelled as:

(molN h-1)

where kt is the empirically determined coefficient for the transport of each gas between the headspace and the liquid, kH(N) (molN L-1 atm-1) is the solubility of NO, N2O, or N2 in water at 20°C, PN (= [N]g×R×T, atm) is the partial pressure of each gas in the headspace, and [N]aq (mol L-1) represents the concentration of each gas in the liquid-phase.

The amount of NO and N2O in the headspace (, molN, Fig 2C) is a function of transport (Eq 29) and the disturbance by gas sampling. The latter is simulated as discrete events at time-points given as input to the model (equivalent to the sampling times in the experiment):

(molN h-1)

where ΔNOx(S) is the net change in the amount of (molN), D (dilution) is the fraction of each gas removed and replaced by equal amount of He, and ts (h) is the time taken to complete each sampling. For N2, the model ignores the sampling loss because the N2 production data to be compared with the model output are already corrected for the sampling disturbance [30]. Thus, the model estimates somewhat higher N2 concentrations than that experienced by the organisms, which is acceptable, since the concentration of N2 is unlikely to have consequences for the metabolism.

Parameterisation

Most of the parameter values used in the model are well established in the literature (see Table 2); however, uncertain parameters include , , , and .

Table 2

Model parameters.

		Description	Value	Units	Reference
Butyrate treatments
vemaxTCA−		Max. cell-specific rate of e--delivery from the TCA cycle	1×1014	mol e- cell-1 h-1	[18]
vemaxO_2−		The maximum cell-specific velocity of e--flow to O2	4.22×10−15	mol e- cell-1 h-1	Optimisation
vemaxNO3−−		The maximum cell-specific velocity of e--flow to NO3−	1×10−14	mol e- cell-1 h-1	[18]
vemaxNO2−−		The maximum cell-specific velocity of e--flow to NO2−	2.65×10−15	mol e- cell-1 h-1	[18]
vemin−		The min. velocity of e--flow to O2/NOx−/NOx required for protein synthesis (ATP)	1.87×10−17	mol e- cell-1 h-1	Assumption
Ye_O_2−		The growth yield per mole of electrons transferred to O2	2.74×1013	cells (mol e-)-1	[18]
Ye_NO_x−		The growth yield per mole e- to NO3−, NO2−, NO, or N2O	1.12×1013	cells (mol e-)-1	[18]
Succinate treatments
vemaxTCA−		Max. cell-specific rate of e--delivery from the TCA cycle	9.34×10−15	mol e- cell-1 h-1	[18]
vemaxO_2−		The maximum cell-specific velocity of e--flow to O2	4.42×10−15	mol e- cell-1 h-1	[18]
vemaxNO3−−		The maximum cell-specific velocity of e--flow to NO3−	9.34×10−15	mol e- cell-1 h-1	[18]
vemaxNO2−−		The maximum cell-specific velocity of e--flow to NO2−	2.01×10−15	mol e- cell-1 h-1	[18]
vemin−		The minimum velocity of e--flow to O2/NOx−/NOx required for protein synthesis (ATP)	1.95×10−17	mol e- cell-1 h-1	Assumption
Ye_O_2−		The growth yield per mole of electrons transferred to O2	4.97×1013	cells (mol e-)-1	[18]
Ye_NO_x−		The growth yield per mole e- to NO3−, NO2−, NO, or N2O	1.52×1013	cells (mol e-)-1	[18]
Parameters common to both succinate and butyrate treatments
[O2]na	The [O2] in aqua below which Nar production triggers		5.95×10−5	mol L-1	[18]
[O2]ni	The [O2] in aqua below which NirS production triggers		9.75×10−6	mol L-1	[18]
rNa	The specific-probability for Nar production		0.035	h-1	Optimisation
rNi	The specific-probability for NirS production		0.004	h-1	Optimisation
vemaxNO−	The maximum cell-specific velocity of e--flow to NO		3.56×10−15	mol e- cell-1 h-1	[33]
vemaxN_2O−	The maximum cell-specific velocity of e--flow to N2O		5.5×10−15	mol e- cell-1 h-1	[24]
K_mO_2	The half-saturation constant for O2 reduction		2.25×10−7	mol L-1	Optimisation
K_mNO3−	The half-saturation constant for NO3− reduction		5×10−6	mol L-1	[34,35]
K_mNO2−	The half-saturation constant for NO2− reduction		4.13×10−6	mol L-1	[36,37]
K1NO	The equilibrium dissociation constant for cNor/NO complex		8×10−14	mol L-1	[33]
K2NO	The equilibrium dissociation constant for cNor/(NO)2 complex		34×10−9	mol L-1	[33]
K_mN_2O	The half-saturation constant for N2O reduction		5.93×10−7	mol N2O-N L-1	Optimisation
D	Dilution (due to sampling): fraction of gas replaced by He		0.013–0.016	–	[18]
k_H(O_2)	Solubility of O2 in water at 20°C		0.0014	mol L-1 atm-1	[38]
kH(NO)	Solubility of NO at 20°C		0.0021	mol L-1 atm-1	[30]
k_H(N_2O)	Solubility of N2O at 20°C		0.056	mol N2O-N L-1 atm-1	[38]
k_H(N_2)	Solubility of N2 at 20°C		0.00035	mol N2-N L-1 atm-1	[38]
kt	The coeff. for gas transport between headspace and liquid		3.6	L vial-1 h-1	Measured
O2leak	O2 leakage into the vial during each sampling		2.92×10−9	mol	Measured
R	Universal gas constant		0.083	L atm K-1 mol-1	–
T	Temperature		293.15	K	[18]
ts	The time taken to complete each sampling		0.017	h	[30]
Volg	Headspace volume		0.07	L	[18]
Volaq	Aqueous-phase volume		0.05	L	[18]

(Eq 17). Pa. denitrificans has three haem-copper terminal oxidoreductases [39] with ranging from nM to µM [40,41], so we decided to estimate the parameter value by optimising for the low [O2] treatments data. Vensim was used for the optimisation, where = 2.25×10−7 neatly simulated the O2 depletion for both the succinate- and butyrate-supplemented treatments.

In vitro studies of NosZ from Pa. denitrificans estimate the values for = 5 μM at 22°C and pH 7.1 [42] and 6.7 μM at 25°C and pH 7.1 [43]. When our model was simulated with in this range, given our empirically estimated [24], the simulated N2O reached concentrations much higher than that measured (see Results/Discussion). A more adequate parameter value (= 0.6 μM) was found by optimising in Vensim. The value is within the range determined for soil bacterial communities [44].

(Eq 17) could be estimated using the empirically determined cell yield per mole of electrons to O2 (, cells per mol e-) and the maximum specific growth rate (μ, h-1): . We are confident about the yields for the two C-substrates used, but the empirically determined μ for the butyrate treatments is suspiciously low (= 0.067 h-1), providing = 2.45×10−15 mol e- cell-1 h-1. Simulations with this value grossly underestimated the rate of O2 depletion compared to measured, which forced us to estimate the parameter value by optimisation: = 4.42×10−15 and 4.22×10−15 mol e- cell-1 h-1 for the succinate- and butyrate treatments, respectively. These values give μ = 0.22 and 0.12 h-1, respectively: for the succinate treatments, the value is very close to that empirically determined (= 0.2 h-1); for the butyrate treatments, the value seems more realistic than 0.067 h-1.

(Eqs 2, 5 and 7) is the per cell velocity of e--flow to O2 () assumed to generate minimum ATP required for synthesising the initial molecules of denitrification enzymes. Since we lack any empirical or other estimations for this parameter, it is arbitrarily assumed to be the when [O2]aq reaches 1 nM. At this concentration, is determined by the Michaelis-Menten equation , using and given above. The values obtained for the succinate- and butyrate-supplemented treatments = 1.96×10−17 and 1.87×10−17 mol e- cell-1 h-1, respectively, which for both the cases is 0.44% of . To investigate the impact of on the model behaviour (rNa and rNi, Eqs 1, 2, 4, 5, 6 and 7), sensitivity analyses were performed by simulating the model with corresponding to [O2]aq = 5×10−9, 5×10−10, and 0 mol L- 1 (see Results/Discussion).

Results/Discussion

Low probabilistic initiation of nar transcription, resulting in the fraction of the population with Nar < 100%

To test the assumption of a single homogeneous population with all cells producing Nar in response to O2 depletion, we simulated the model with the specific probability for a Z− cell to initiate nar transcription (rNa) = 4 h-1. This resulted in 98% of the cells possessing Nar within an hour (see Eqs 1–3). Evidence suggests that less than half an hour is required to synthesise denitrification enzymes [17,18], but an hour’s time is assumed here to allow margin for error. The results show that, for all the treatments, the simulated production (mol vial-1) grossly overestimates that measured (Fig 3).

Fig 3

Comparison of measured and simulated accumulation assuming definitive versus stochastic initiation of nar transcription.

To test the assumption of a single homogeneous population with almost all cells expressing nar in response to O2 depletion, we forced our model to achieve 98% Nar-positive cells (ZNa) within an hour by setting the specific-probability of initiating nar transcription (rNa) = 4 h-1. This resulted in grossly overestimated rates of accumulation for all treatments (grey curves). In contrast, we simulated the model with rNa = 0.035 h-1 obtained through optimisation, resulting in a reasonable agreement with measurements for all treatments (except for an apparent time frameshift for the Butyrate, 7% O2 treatment).

To find a reasonable parameter value, we optimised rNa for the 0% O2 treatments, so that the simulated production matches that measured. The results (Table 3) suggest that a low probabilistic initiation of nar transcription (average rNa = 0.035 h-1) is adequate to simulate the measured kinetics (Fig 3). In the Butyrate, 7% O2 treatment (Fig 3B), the simulated starts earlier, but the rate of accumulation is similar to that measured.

Table 3

Specific-probability of nar and nirS transcriptional initiation (rNa and rNi, respectively) estimated for each treatment by optimisation (best match between the simulated and measured data).

Batch	C-source	Treatment*: O2 (vol.%), NO3− (mM)	Optimal rNa (h-1)	Optimal rNi (h-1)
1	Butyrate	~0, 2	0.041	0.005
2	Butyrate	7, 2	–	0.004
3	Succinate	~0, 2	0.030	0.005
4	Succinate	7, 2	–	0.003
			Avg. = 0.035	Avg. = 0.004

*Treatment refers to the C-source, initial oxygen concentration in the headspace (measured as headspace-vol.%), and initial concentration in the medium (mM).

Once O2 falls below a certain threshold, the production of Nar is assumed to trigger with rNa = 0.035 h-1 and last until a minimum of respiration is sustained by the e--flow to O2 and N2O ( and ), assumed to fulfil the ATP needs for Nar production (Eqs 1 and 2). But the production of Nar sustained by was inconsequential for simulating the measured production, since was already exhausted when N2O started accumulating (i.e., when > 0). For this reason, the fraction that produced Nar (FNa, Eq 3 and Table 4) is calculated as functional (= 0.23–0.43) and theoretical (= 0.56–0.81), where the first is the fraction actually responsible for production (sustained by ), but the latter also incorporates the fraction that produced Nar after the exhaustion of (sustained by as well as ). The rationale behind calculating the theoretical FNa is the empirical data indicating that Nar transcription is not turned off in response to depletion [18]. Although our model cannot test the theoretical FNa, but the functional FNa suggests that, contrary to the common assumption, the measured kinetics can be neatly explained by only 23–43.3% of the population producing Nar in response to O2 depletion.

Table 4

The fraction of the population with Nar (FNa) and NirS (FNi) estimated based on the optimal specific-probability of nar and nirS transcriptional initiation (rNa and rNi), respectively.

Batch	C-source	O2 (vol.%), NO3− (mM)	Functional FNa * (unitless)	Theoretical FNa * (unitless)	FNi (unitless)
1	Butyrate	~0, 2	0.433	0.813	0.221
2	Butyrate	7, 2	0.343	0.656	0.088
3	Succinate	~0, 2	0.357	0.803	0.206
4	Succinate	7, 2	0.230	0.564	0.077

*Functional FNa is the fraction of cells expressing Nar while is still present, while Theoretical FNa is the fraction expressing Nar when including the theoretical recruitment after depletion (supported by energy from N2O reduction).

Very low probabilistic initiation of nirS transcription

When we optimised the specific probability of nirS transcriptional activation (rNi, see Eqs 4, 5, 6 and 7) to fit the measured data, the average rNi = 0.004 h-1 (Table 3) adequately simulated the measured depletion and N2 accumulation (Fig 4). The recruitment to denitrification lasted for 19.5–47.3 h, i.e., the time when [O2] was below a critical concentration and the velocity of e--flow to O2 and the relevant /NOx remained above a critical minimum (Eqs 4, 5, 6 and 7). The resulting fraction recruited to denitrification (FNi, see Eq 8 and Table 4) was 0.08–0.18, the bulk of which depended on the e--flow to and N2O (instead of aerobic respiration).

Fig 4

Comparison of measured and simulated data assuming stochastic initiation of nirS transcription.

Each panel compares the measured depletion (sub-panel) and N2 accumulation (main panel; n = 3–4) with simulations. The simulations are carried out with an optimised specific-probability of nirS transcriptional initiation (average rNi = 0.004 h-1, Eqs 4, 5, 6 and 7), allowing 7.7–22.1% of the population to produce NirS + cNor (Eq 8) during the available time-window (= 19.5–47.3 h).

To test whether the measured data could be explained without the recruitment sustained by and N2O respiration, we also simulated the model with the recruitment as a function of O2 alone and re-optimised rNi, which on average increased to 0.012 h-1 (providing FNi = 0.083–0.35). This was expected since O2 is exhausted rather quickly, shrinking the time-window available for the recruitment. Comparatively, these simulations were less satisfactory: using the average rNi = 0.012 h-1 generally resulted in larger deviations than for the default simulations (S2 Fig), and the optimal rNi for individual treatments varied grossly (50% higher values for the ~0% O2 treatments than for the 7% O2 treatments). This contrasts the default simulations, where the optimal rNi values for individual treatments were quite similar.

Sensitivity of rNa and rNi to

Recruitment to denitrification (both nar and nirS transcription) is assumed to continue only as long as the combined e--flow to O2, and N2O is greater than (Eqs 1, 2, 4, 5, 6 and 7). To test the model’s sensitivity to this parameter, we estimated rNa and rNi by optimisation for different values of , relative to the default value = 1.95×10−17 mol e- cell-1 h-1. For all cases, the model was able to adequately simulate the measured N2 kinetics by moderate adjustments of rNa and rNi. Table 5 shows the average optimal values of rNa and rNi, obtained by fitting the simulated N2 kinetics to the data, for different values of . S3 Fig shows adequate simulations of the measured N2 kinetics assuming = 0, with optimised rNa = 0.033 h-1 and rNi = 0.0033 h-1. Thus, although assuming > 0 appears logical, it is not necessary to explain the measured data.

Table 5

Estimated rNa and rNi, depending on a.

vemin− (mol e- cell h-1)	Optimal rNa (h-1)	Optimal rNi (h-1)
5 × Default *	0.041	0.0062
Default	0.035	0.0041
0.5 × Default	0.034	0.0035
0	0.033	0.0033

*Refers to the default value = 1.95×10−17 mol e- cell-1 h-1.

N2O kinetics

To simulate the N2O kinetics, we used = 5.5×10−15 mol e- cell-1 h-1, empirically determined under similar experimental conditions as simulated here [24], and adopted the literature values for [= 5 and 7 μM 42,43, respectively]. But with = 5 μM, the model predicted N2O accumulation ~10–20 times higher than measured for the ~0% and ~2–3 times higher for the 7% O2 treatments (Fig 5). This forced us to simulate the model with the parameter value estimated by optimisation, providing the average = 0.6 μM.

Fig 5

Comparison of the measured N2O with that simulated.

Each main panel (A–D) compares the measured N2O (single vial results) with the default simulation using the parameter values given in Table 2, i.e., = 0.6 μM (estimated through optimisation) and = 5.5×10−15 mol e- cell-1 h-1 [24]. In contrast, each inserted panel shows the simulated N2O assuming 1) N2O consumption only by the cells producing N2O (ZNaNi + ZNi), and 2) the literature value for = 5 μM [42]. The results show that the default simulation best explains the measured N2O kinetics, assuming its production by a small fraction (ZNaNi + ZNi) and consumption by the entire population (Z− + ZNa+ ZNaNi + ZNi).

The measured N2O shows a conspicuous increase throughout the entire active denitrification period, and this phenomenon is neatly captured by the model. The reason for this model prediction is that the number of N2O producing cells (ZNaNi + ZNi, Fig 2A) is low to begin with compared to the number of N2O consuming cells (Z− + ZNa+ ZNaNi + ZNi), but the fraction of N2O producers will increase during the anoxic phase for two reasons: one is the recruitment to ZNaNi & ZNi, another is the fact that the model predicts approximately three times faster cell-specific growth rate for ZNaNi & ZNi than for Z− & ZNa ( is identical for all groups, while and are both zero for Z− & ZNa but for ZNaNi & ZNi, it holds that ≈ > . To illustrate this phenomenon, we ran the model, assuming that the Z− & ZNa cells had no N2O reductase, resulting in a) constant N2O concentration throughout the entire anoxic phase and b) much higher N2O concentrations than measured (Fig 5). The overestimation is a trivial result, easily avoidable by increasing or decreasing moderately. However, the prediction of a constant N2O concentration is clearly in conflict with the experimental data, and no parameterisation could force the model to reproduce this phenomenon other than the differential expression of denitrification genes.

Hence, although there is room for further refinements, our default assumption regarding differential expression of NirS and NosZ explains the observed N2O kinetics: 1) abrupt initial accumulation to very low levels due to recruitment of relatively small numbers to the N2O producing pools (ZNaNi & ZNi), and 2) increasing N2O concentration due to recruitment and faster cell-specific growth of ZNaNi & ZNi than that of the cells only consuming N2O (Z− + ZNa).

This modelling exercise sheds some light on the possible role of regulatory biology of denitrification in controlling N2O emissions from soils. If all cells in soils had the same regulatory phenotype as Pa. denitrificans, their emission of N2O would probably be miniscule, and soils could easily become strong net sinks for N2O because the majority of cells would be ‘truncated denitrifiers’ with only N2O reductase expressed. It remains to be tested, however, if the regulatory phenotype of Pa. denitrificans is a rare or a common phenomenon among full-fledged denitrifiers. We foresee that further exploration of denitrification phenotypes will unravel a plethora of response patterns.

Conclusion

Using dynamic modelling, we have demonstrated that the denitrification kinetics in Pa. denitrificans can be adequately explained by assuming low probabilistic transcriptional activation of the nar and nirS genes and a subsequent autocatalytic expression of the enzymes. Such autocatalytic gene expressions are common in prokaryotes, rendering a population heterogeneous because of the stochastic initiation of gene transcription, with a low probability [45]. For N2O kinetics, our hypothesis was that a) the gas is produced by a fraction of the incubated population that is able to initiate nirS transcription with a certain probability, leading to a coordinated expression of nirS + nor via NO [16], and b) N2O is consumed by the entire population because, in response to anoxia, nosZ is readily induced by FnrP [24]. Our model corroborated this hypothesis by reasonably simulating the N2O kinetics with the specific-probability of nirS transcriptional activation = 0.004 h-1, resulting in 7.7–22.1% of the population producing NirS + cNor (hence N2O), but all cells producing NosZ (hence equally consuming N2O).

The folder contains the dynamic model used in this study ‘Hassan_et_al_2015_Pa._denitrificans.mdl’.

The model requires Vensim (Double Precision), which is available at http://vensim.com/free-download/. The zip folder also contains files with the empirical data; these files are automatically loaded into the model when it is run.

(ZIP)

Click here for additional data file.

Pa. denitrificans gas and kinetics.

Typical gas kinetics (O2, NO, N2O, N2) and accumulation in Pa. denitrificans during the transition from aerobic respiration to denitrification; batch cultures, n = 3; 20°C; Sistrom’s medium; 2 mM KNO3 and 7 vol% initial O2 in the headspace. All the available (100 μmol vial-1) was recovered as before the onset of N-gas production. In previous experiments [17], N2O concentrations were below the detection limit of the system, but thanks to a new system with electron capture detector, the N2O kinetics were monitored with reasonable precision. Adapted from [18].

(TIF)

Click here for additional data file.

Comparison of measured and simulated data assuming stochastic initiation of nirS transcription with aerobic respiration being the only energy source for producing NirS + cNor.

In each panel, the measured depletion (sub-panel) and N2 accumulation (main panel; n = 3–4) are compared with simulations. The simulations here are to be compared with the default simulations (Fig 4), which were run assuming that the coordinated NirS + cNor production (via nirS transcriptional activation) is sustained by the energy generated by O2 as well as and/or N2O reduction. The default simulations provided an average specific-probability of nirS transcriptional activation (rNi) = 0.004 h-1 (Eqs 4, 5, 6 and 7) by optimisation, allowing 7.7–22.1% of the population to produce NirS + cNor (Eq 8) in 19.5–47.3 h. To match the measured data here, the average rNi had to be raised to 0.012 h-1, since the time available for the enzyme synthesis shrank (= 3.5–16 h) due to a rapid exhaustion of O2. Comparatively, the assumption that the ATP from and/or N2O reduction should help cells produce denitrification enzymes seems more plausible and provide better agreement with the measured data.

(TIF)

Click here for additional data file.

Measured vs. simulated N2 kinetics assuming = 0.

The default simulations are carried out assuming that for a cell to produce first molecules of Nar and NirS, a minimum of e--flow to an available e--acceptor (, mol e- cell-1 h-1) is necessary to generate a minimum of ATP required for protein synthesis (Eqs 1, 2, 4, 5, 6 and 7). Although assuming > 0 seems logical, the measured N2 kinetics are adequately simulated here with = 0. This shows that the assumption is not necessary to explain the measured data.

(TIF)

Click here for additional data file.