A diazotrophy-ammoniotrophy dual growth model for the sulfate reducing bacterium Desulfovibrio vulgaris var. Hildenborough

Sulfate reducing bacteria (SRB) comprise one of the few prokaryotic groups in which biological nitrogen fixation (BNF) is common. Recent studies have highlighted SRB roles in N cycling, particularly in oligotrophic coastal and benthic environments where they could contribute significantly to N input. Most studies of SRB have focused on sulfur cycling and SRB growth models have primarily aimed at understanding the effects of electron sources, with N usually provided as fixed-N (nitrate, ammonium). Mechanistic links between SRB nitrogen-fixing metabolism and growth are not well understood, particularly in environments where fixed-N fluctuates. Here, we investigate diazotrophic growth of the model sulfate reducer Desulfovibrio vulgaris var. Hildenborough under anaerobic heterotrophic conditions and contrasting N availabilities using a simple cellular model with dual ammoniotrophic and diazotrophic modes. The model was calibrated using batch culture experiments with varying initial ammonium concentrations (0–3000 µM) and acetylene reduction assays of BNF activity. The model confirmed the preferential usage of ammonium over BNF for growth and successfully reproduces experimental data, with notably clear bi-phasic growth curves showing an initial ammoniotrophic phase followed by onset of BNF. Our model enables quantification of the energetic cost of each N acquisition strategy and indicates the existence of a BNF-specific limiting phenomenon, not directly linked to micronutrient (Mo, Fe, Ni) concentration, by-products (hydrogen, hydrogen sulfide), or fundamental model metabolic parameters (death rate, electron acceptor stoichiometry). By providing quantitative predictions of environment and metabolism, this study contributes to a better understanding of anaerobic heterotrophic diazotrophs in environments with fluctuating N conditions.


Introduction
Nitrogen (N) is the fourth most abundant element in biological systems, and its availability often constrains primary production in terrestrial [17,37,57] and oceanic ecosystems [18,55,60]. A phylogenetically diverse group of prokaryotes have evolved the specialized ability to acquire new N from the large atmospheric N 2 reservoir using a process termed biological nitrogen fixation (BNF) catalyzed by the enzyme nitrogenase. The ecological advantage provided by BNF has allowed N fixers to successfully colonize a diversity of environments characterized by low N availability where they help support primary production. Over the last two centuries, the development of human societies has substantially impacted N cycling, with reactive forms of N entering natural environments at higher rates, widely alleviating N limitation but often favoring invasive species [48], and contributing to long term greenhouse gas production [38,7]. A better understanding of how increases in N availability will influence the activity and adaptability of N-fixing organisms, particularly in carbon (C)-rich ecosystems like fresh water and marine sediments, is important to understand their resilience and contribution to ecosystem services in the long term under anthropic pressure.
Sulfate reducing bacteria (SRB), most of which are deltaproteobacteria, have well-accepted roles in benthic biogeochemistry [21,3,4,41,44] based on their use of sulfate as the terminal electron acceptor for oxidative transformations of organic carbon and hydrogen (H 2 ). Widespread findings of diverse SRB-like nif genes, which encode for the most ubiquitous form of nitrogenase, within sediments with relatively high fixed-N availability have been surprising, prompting questions on the role and regulating controls on benthic BNF [36,40,41].
One fundamental reason why fixed-N is generally thought to be preferred over BNF lies in the higher metabolic cost of BNF [15,22]. For instance, NH 4 + transport by the most abundant amtB transporter could cost only 1 ATP, the cost of one ion transport, or perhaps even no energy expenditure if transport occurs in the deprotonated form, NH 3 [31,32,61]. In contrast, the chemical bond between the two N atoms in N 2 is among the strongest, and 16 moles of ATP are required for the biochemical conversion of one mole of N 2 into 2 moles of NH 4 + , with the obligate production of 1 mol of H 2 as a by-product [50,51]. Based on the direct costs of N acquisition, ammonium is thus possibly ∼ 8 times less costly a source of N than BNF. Additionally, indirect metabolic costs such as biosynthesis and enzyme maintenance in N and related pathways can substantially impact the overall cost of N acquisition for an organism [23]. Further, diazotrophs also commonly possess transporters to access "cheaper" forms of nitrogen such as ammonium (NH 4 + ) [22,35].
Yet, several studies in benthic environments have reported the presence of nif genes, coding for the canonical Mo-nitrogenase, in areas with relatively high bulk levels of ammonium, which have raised questions on the possible contributions of BNF to N input in the presence of NH 4 + and on the spatio-temporal complexity of benthic N cycling [36]. In support of the traditional conceptualization of BNF as a less-preferred N source, recent findings showed clear downregulation of benthic BNF activity in the presence of ammonium [11,13]. However, measured NH 4 + thresholds for BNF inhibition in naturally heterogeneous benthic sediments were relatively high (< 30 µM NH 4 + , [13], 0-20 µM NH 4 + , [11]) as compared to cellular thresholds in liquid culture (< 2 µM, [11]). This difference suggests that sediment bio-physical complexity (e.g., nutrient diffusion, remineralization, flow heterogeneity) could significantly obscure our ability to observe "true" sensitivities in various sediment types. Similarly, a mechanistic understanding of how in-situ heterogeneity influences biological activities has been limited by our capacities for measuring key phenomena at microbiologically-relevant time and spatial scales. By disentangling biology from physicochemical factors, quantitative growth models could help provide a refined mechanism of physical, chemical and metabolic interactions in real ecosystems [14,30].
In order to better quantify and predict metabolic differences in nitrogen acquisition by anaerobic benthic diazotrophs, we developed a simple cellular growth model for the anaerobic heterotroph Desulfovibrio vulgaris var. Hildenborough (DvH) [25], a well-known sulfate reducer, which resolves both NH 4 + uptake and BNF as N sources. This is the first dual N source model for anaerobic heterotrophic N fixers. We parameterized the model using measurement data from the growth of DvH in batch cultures with increasing initial [NH 4 + ] ranging from the background (∼10 µM) to 3000 µM and literature values of ammonium transporter enzyme kinetics [10,35]. Finally, we discuss the implications of our findings for understanding N flows in benthic environments and the potential applications and limitations of the proposed model. By providing quantitative mechanistic predictions of environment and metabolism, this study contributes to a better understanding of anaerobic diazotrophs in spatiotemporally changing sediment and soil environments.

Diazotrophic -ammoniotrophic dual growth model (DAD-GM)
The DAD-GM is an idealized model that allows a population of bacterial cells to grow on both ammonium and atmospheric dinitrogen as N sources (Fig. 1). In the model, the growth of new biomass (Cell), expressed as ODmL, (i.e., the amount of biomass in 1 milliliter of culture at a density of 1 unit of OD 600 , optical density at 600 nm), is determined by the acquisition of nitrogen (N) by ammonium transporters and nitrogenase enzyme activities following Eqs. 1  were calibrated experimentally (see Table 2 and below). (1) The main features of batch culture experiments with DvH under sub-replete initial ammonium concentrations ([NH 4 + ] initial ) (< 2300 µM, 30 mM of Pyruvate, 10 mM Sulfate, [11]) are summarized in Fig. 1A. Cultures exhibited a clear bi-phasic growth, with the initial growth phase supported by ammonium uptake (first ∼ 75 h.) followed by a second, diazotrophic growth phase. The two growth phases were separated by a lag phase (> 24hrs., Fig. 1A), indicating the diazotrophic regime could be modelled independently from the ammoniotrophic regime as a second, separate exponential growth curve. Hence, in our model formulation, the total active cell number (Cell tot , expressed in ODmL) can be separated into two cell groups with distinct functional characteristics. One group consists of idealized specialized N-fixing cells (Cell BNF , expressed in ODmL), which contain active NH 4 + transporters and a fixed and optimal number of nitrogenase enzyme required to support maximal BNF activity. ] > 3000 µM) experimental data. At each time increment (dt), Cell BNF and Cell NH4 pools acquire nitrogen according to the pathway described above and time-dependent enzyme-specific activities for nitrogen fixation and ammonium transport (see Table 1). The amount of each specialized biomass type (i.e., Cell BNF and Cell NH4 ) is increased using cell-specific cellular N quotas (Q N,BNF and Q N,NH4 , both expressed in µmol N .ODmL -1 , for diazotrophic and ammoniotrophic cells, respectively), and the composition (NH 4 + , SO 4 2-, and pyruvate) of the surrounding medium is then balanced to account for the increase in biomass and by removing stochiometric amounts of sulfate and pyruvate according to measured and theoretical stoichiometries ( Fig. 1 and Table 1).
To study the dynamics of N fixation activity following NH 4 + addition, we used acetylene reduction assays (ARA), a common method to measure nitrogenase activity [24]. To allow the model to integrate data from ARA experiments, we accounted for the use of acetylene as a nitrogenase substrate in diazotrophic cells. We assumed total nitrogenase activity was partitioned between acetylene reduction and nitrogen reduction according to nitrogenase acetylene saturation level (Sat = [Ac]/([Ac]+K M,Ac ), expressed in %, where [Ac] is acetylene headspace concentration in %v/v, and K M,Ac the Michaelis-Menten constant of nitrogenase for acetylene, in %v/v). Hence, during ARA experiments, diazotrophic cells acquire nitrogen to support the increase of the Cell BNF pool at only (1-Sat) % of the measured nitrogenase activity v BNF (see Table 1). According to measured delays in ammonium inhibition [11], we incorporated in the model a 30 min delay between NH 4 + addition and the occurrence of the transition between diazotrophy and ammoniotrophy.
In a first set of model simulations, the DAD-GM was applied to reproduce batch culture growth data using the initial media conditions of our laboratory experiments (see below) with increasing Table 1 Mass balance equations used to model the growth of DvH on ammonium (NH 4 + ) and by biological nitrogen fixation (BNF) following a time (t) increment of dt (i.e., from t to t + dt). Index i refers to either BNF or NH 4 + . Variables a The nitrogenase acetylene saturation (Sat, in %), the decay rate (r death , in hrs -1 ), and pathway specific nitrogen cellular quota requirement (Q N,i, in mol N .ODmL -1 ) are defined in the text (see also Fig. 1 A and 2 B, and Table 2). Enzyme specific activities (v BNF and v NH4 ) are defined in Eqs. 1 and 2. V is the media volume, expressed in mL. Cellular biomass was expressed as OD x mL, or the number of cells found in 1 mL of medium at OD 600 = 1.
[NH 4 + ] initial (0-3000 µM, 30 steps). The DAD-GM was iterated for a total simulated time of 300 h with a time step (dt) of 0.5 h (600 steps). We performed a second type of simulation to investigate growth, BNF activity, and ammonium uptake during the addition of increasing [NH 4 + ] (0-3000 µM) to actively growing diazotrophic culture under acetylene reduction conditions. To facilitate comparison with laboratory experiments, where volume could change with media sampling, results for cell number are expressed as cell density (OD 600

Batch culture growth experiments and biomass measurements
We conducted growth experiments containing different initial concentrations of ammonium (NH 4 + ), which ranged from background concentration (< 10 µM) to 3000 µM. Ammonium was added as ammonium chloride to the medium from 10-to 100-fold stock solutions prepared in aseptic culture medium to reach the required final initial medium concentration. Growth experiments were initiated at an initial OD 600 ∼ 0.05 by needle addition of 100-300 mL of a starter culture grown under the same conditions for at least 10 generations. Cells were grown at 30 o C in the dark, with tubes placed obliquely in an orbital shaker at a constant agitation speed of 130 rpm. Culture cell density was followed by turbidimetry at 600 nm using a Thermo Fisher Scientific (Waltham, MA, USA) GENESYS spectrophotometer equipped with a tube holder. Each experiment was performed in at least triplicate. The typical results of our growth experiments are summarized and explained in Fig. 1A.

Acetylene reduction assay experiments
Data from experiments using the Acetylene Reduction Assay (ARA) described in a prior publication [11] were used to calibrate nitrogenase maximal specific activity (v BNF,MAX ) and to validate ethylene production, growth, and ammonium uptake under the ARA condition. Briefly, the experiments consist of diazotrophic cultures of DvH in pyruvate (30 mM) or lactate (40 mM) with replete sulfate as the electron acceptor (20 mM). Once cultures reached an OD 660 ∼ 0.08 ± 0.01, the ARA experiments were started by adding acetylene in the headspace. Since acetylene directly inhibits nitrogen fixation, we used two different acetylene concentrations depending on the experimental needs. During experiments designed for calibrating nitrogenase maximal specific activity, we employed a close to saturating acetylene condition for molybdenum nitrogenase according to available kinetic data (5.9% v/v, [10]) and pyruvate, a fast-growing substrate, to ensure high rates of ethylene production and thus reduced measurement uncertainty. Conversely, we used a lower acetylene concentration (1.7% v/v) in a second type of experiment with lactate, a slower growing substrate, to reduce the acetylene inhibition of nitrogen fixation, ensure maximal cell growth, and allow sufficient time to conduct sampling at frequent timepoints. In all experiments, once we confirmed ethylene production was active (after ∼ 3 h), ammonium was introduced as a single addition into actively fixing cultures (final concentration 10-3000 µM) and headspace aliquots from each culture were sampled for the next 100-120 h and analyzed for their ethylene content (see [11] for additional details).

Calibration of diazotrophic growth phase in the model with Acetylene Reduction Assays
Nitrogenase maximal specific activity (v BNF,MAX ) was estimated in pyruvate experiments conducted with ∼5.9 %v/v acetylene headspace concentration ([Ac]), assuming an enzyme saturation (Sat) of 60% (with a K M , Ac ∼ 0.04 atm, [12]) and a theoretical R ratio (i.e., the ratio of acetylene to dinitrogen reduction rate) of 3.2 for the Mo nitrogenase [2]. Values were averaged across replicates for each timepoint, leaving out conditions in which ammonium was present in the media and the last time point when growth reached stationary phase, for a total of 16 individual cultures and 51 measurements over four timepoints (40 ± 18 nmol N .hr -1 . ODml -1 , range 25-65 nmol N .hr -1 . ODml -1 , SI Dataset S1, VBNF tab).

Measurement of sulfate utilization during growth
The quantity of sulfate required during diazotrophic and ammoniotrophic growth in this study was estimated by measuring dissolved hydrogen sulfide concentrations during and at the end of growth experiments using the methylene blue method [8]. At multiple times along the experiment, at the end of visible growth, (i.e., maximal biomass yield), and during the observed decay phase (∼48 hrs. after maximal yield), aliquots of medium were sampled anaerobically using needle and syringe, filter sterilized, and preserved in a 2% zinc acetate solution at 4 o C for < 7 days until measurement. The extinction coefficient was obtained from external calibrations. The calibration (slope=0.028 ± 0.0003 mmol.L -1 .cm -1 , range=0.1-6 µ M, dilution rate 250X, n = 5) was obtained using sodium sulfide nonahydrate (Sigma Aldrich, Allentown, PA, USA) prepared in a serum bottle sealed with a 22 mm blue butyl septum. A selected crystal of appropriate mass was properly cleaned with MilliQ® water (Sigma Millipore, West Springfield, MA, USA), thoroughly dried using absorbent paper, weighed, and dissolved in a known volume of MilliQ® water. Total H 2 S (H 2 S tot ) produced (a proxy for total SO 4 2used during growth) was calculated by correcting HSaq for partition between chemical species (S 2-, HS -, and H 2 S, pKa 1 =7) and neglecting dissolved H 2 S (g) using Henry's constant in water Hcc (water/air) = 0.0036 [47], at 30 o C and pH 7.1 (final correction factor applied: HSaq / H 2 S tot = 0.44). For each replicate, the highest H 2 S tot value measured was used to calculate the mean value of the condition.

Model parametrization using experimental data
We estimated the key parameters of the model from experimental data using regression analyses between biomass production and sulfate usage ( Fig. 2A), biomass decay in relation to biomass density (Fig. 2B), and sulfate and nitrogen utilization during growth (Fig. 2C&D) Fig. 2A and Table 2). We estimated r death (in hr -1 ), the rate at which the biomass was removed from the media due to senescence or death (see Fig. 1A Fig. 2D).

Validation of calibration parameters
We confirmed the validity of the model parameters and output using independently measured variables as well as theoretical and literature values. Our independent estimates for SO 4  We directly estimated Q N,BNF from the ARA experiments, by dividing the total amount of N fixed during the experiments in each individual culture (dmol N , assuming 40%, or 1-60% saturation, for enzyme activity toward N fixation, SI Dataset S1, QNBNF tab) by the total biomass created corrected for ammonium contribution (dODmL). The resulting value was at the higher end of the calculated estimate using v BNF,max (1.62 + 0.14 vs. a range of 0.6-1.5 umol N . ODmL -1 , see Table S1, SI Dataset S1, QNBNF tab).
Finally, we estimated the amount of growth inhibition to DvH due to use of the 5.9%v/v acetylene in ARA experiments by calculating dln (OD)/dt for every timepoint and each individual culture, removing the first timepoint and conditions under which ammonium was present. The average growth rate under 5.9%v/v acetylene led to a value of 0.018 ± 0.007 hr -1 (SI Dataset S1, Mu_ARA tab), or about ∼ 41% of the diazotrophic growth rate under non-ARA conditions. This result supports our assumption of a 60:40 partition of nitrogenase total activity between acetylene reduction and nitrogen fixation and the validity of the K M,Ac value used here (i.e., K M,Ac = 4%v/v, Table 2).

Evaluation of model goodness-of-fit
To assess the performance of the DAD-GM model, we evaluated the goodness of fit using a linear regression of the model output with experimental data for each [NH 4 + ] init condition (Fig. 3) and for each measurement time (Fig. 4&5). For each figure and each variable, we report the goodness-of-fit as the adjusted R-squared of the regression (R 2 ).

DAD-GM simulations of DvH growth under varying levels of ammonium repleteness
To evaluate DAD Growth Model performance, we first compared modeled and measured growth rates, biomass yields (reported as maximal biomass density in 10 mL cultures), and maximum H 2 S tot produced in batch cultures of DvH at increasing [NH 4 + ] init (Fig. 3). Modeled growth rates for the ammoniotrophic and diazotrophic phase were calculated using an approach similar to that applied to laboratory data (i.e., linear regression of the logarithm of the total active cell, Fig. 3A), and we extracted the maximal biomass yield for ammoniotrophic growth (i.e., in relation to the first growth plateau, see Fig. 1A) and complete growth (maximum of Cell tot ) (Fig. 3B).
Maximum produced [H 2 S] tot in the model was also compared to measured values at the end of visible growth (i.e., maximal biomass yield) and in multiple cases 48 h later (Fig. 3C). The model and experimental data for growth rate and final and NH 4 plateau biomass yield showed close agreement with each other (R 2 > 0.85), and modeled [H 2 S] tot values reproduce well the average measured [H 2 S] tot , albeit demonstrating higher variability (R 2 = 0.41).
Assuming the growth experiment was started with fully competent N 2 -trophic cells at an initial biomass density (set to 0.002 OD 600 (i.e., for 10 mL culture, Cell BNF,0 = 0.02 ODmL and Cell NH4,0 =0 ODmL), simulations yielded a good fit between the measured and modeled diazotrophic and ammoniotrophic growth rate for [NH 4 + ] init ranging from 0 to 3000 µM (Fig. 3A, R 2 > 0.85). Changing values of K M,NH4 between 0 and 20 µM, within the range of published K M values for diazotrophs (0-20 µM, [35]), produced similar results. These K M values are of similar magnitude to the measured threshold value for nitrogenase inhibition found by [10] for DvH in pure culture and in slurry incubation ([NH 4 + ] threshold ∼10 µM) and used as K i,NH4 used in this study (Table 2). Moreover, the Monod constant for ammonium K s,NH4 , estimated using non-linear regression of our simulated growth rate data, was found to be K s,NH4 = 32 µM ± 4 µM, slightly higher than the preset value of K M,NH4 of 20 µM within the model (see Fig. 3A).
In addition to the good performance of the model for quantifying cellular parameters and medium final composition (R 2 > 0.85 for growth rate, > 0.95 for biomass density, and ∼ 0.4 for [H 2 S] tot ), the DAD Growth Model closely reproduces (0.56 < R 2 < 0.96) biomass increase with time for increasing [NH 4 + ] init across three independent replicates and using the exact same starting condition (i.e., the same initial biomass quantity of ammoniotrophic and diazotrophic cells). Good visual agreement between modeled and measured data is particularly evident during ammoniotrophic growth phases (see Fig. 4 H-N), but there are some disparities for the diazotrophic growth phases, particularly at low initial [NH 4 + ], which we mostly attributed to fluctuation in the length of the first (see Fig. 4M) or second lag-phase between NH 4 + -trophic and N 2 -trophic growth (see

DAD-GM simulations of DvH diazotrophic growth response to ammonium addition
To further evaluate model capability, we tested whether the DAD-GM could reproduce growth dynamics after ammonium addition to diazotrophically growing DvH, as recorded using ARA experiments in which ammonium was added to growing cultures at various concentrations (0-3000 µM, Fig. 5).
The modeled ethylene production (Fig. 5A, R 2 =0.86), growth curves (Fig. 5B, R 2 =0.92), and dynamics of N fixation inhibition   ([11]). We also applied the same DAD-GM to evaluate a similar experiment using lactate instead of pyruvate as the electron donor (40 mM), and a change in the electron donor: sulfate ratio to fit the stoichiometry of lactate (2:1 instead of 4:1, [43]), and the growth rate of the diazotrophic and ammoniotrophic phase (µ BNF,lac =0.022 hr -1 and µ NH4,lac =0.044 hr -1 ). In both the modeled and measured data, ethylene production, culture growth and ammonium consumption showed similar patterns (Fig. 5 D, E & F, R 2 > 0.9). Of note, there was a 5-10 hrs delay in modeled and measured ethylene production (Fig. 5D). Conversely, the model shows a slower and delayed NH 4 + uptake than our measurements for the 800 µM NH 4 + addition experiment (Fig. 5F). In both cases, the delay was not observed in biomass data (compare Fig. 5 D, E, and F).
In addition, contrary to the model assumption, there was no visible decay phase in the lactate condition (Fig. 5E).

DAD-GM identification of N fixation specific growth limitation
To be able to model the decrease in maximum H 2 S tot produced in experiments with [NH 4 + ] init between 500 and 1500 µM (Fig. 3C), it was necessary to arbitrary impose a limiting step preventing BNF activity in the model. Indeed, under nitrogen deplete conditions, measurements indicate cultures appear to stop growing once total biomass (i.e., ammonium and BNF-supported) reached a density of OD 600 ∼ 0.4 (compare Fig. 3B and C and see Fig. 4). We could not identified the metabolic process underpinning this limitation. It was thus introduced in the model as a condition on the total amount of dead biomass: "no BNF activity once total dead biomass density is over x OD 600 ", and was tuned independently in simulations of experiments with and without acetylene to achieve the best fit between model and data (0.03 and 0.065 OD dead with and without acetylene, respectively).

Metabolic perspectives on cost of nitrogen acquisition: NH 4 + versus N fixation
The measured experimental data, supported by the overall good goodness-of-fit of the DAD-GM (R 2 > 0.85 for growth rates, > 0.95 for biomass density yields, ∼ 0.4 for [H 2 S] tot, > 0.56 for growth curves, > 0.85 for NH 4 + addition experiments), enable comparison of the metabolic cost of diazotrophy vs. ammoniotrophy in DvH ( Table 3). Estimates of ATP cost per unit of biomass can be calculated using direct measurements of sulfate usage per biomass produced (18.9 and 9.7 µmol SO4 . ODmL -1 , Table 3) for BNF and NH 4 -supported growth, respectively. Assuming 1 mol of SO 4 2used per mol ATP produced [5], diazotrophy-derived biomass requires twice as much ATP as NH 4 + -derived biomass (18.9 vs. 9.7 µmol ATP . ODmL -1 ). Assuming 5.5 × 10 8 cell.ODmL -1 (as found at 600 nm for Desulfovibrio desulfuricans, [59]), ∼ 20 billion and 10 billion ATP molecules are required for one cell made under diazotrophic and ammoniotrophic conditions, respectively. In comparison, Escherichia coli, a non-diazotrophic organism of similar shape but smaller than DvH, requires ∼12-20 billion ATP to synthesize one cell (range for various growth conditions, [53]), indicating DvH could have a relatively lower biomass cost than this well-known model organism.
When the overall metabolic cost was reported per mol of N using the N biomass requirement (Q N,BNF = 0.9 and Q N,NH4 = 2.8 mol N .ODmL -1 , Table 2), we calculated an average cost for BNF ∼ 21 ± 3 mol ATP .mol N-BNF -1 (range between 13 and 32 mol ATP .mol N-BNF -1 ), almost three times larger than the direct biochemical cost of 8 mol ATP .mol N-BNF -1 for nitrogen fixation (i.e., 16 mol ATP per mol of N 2(g) reduced, [50]). Direct and indirect costs associated with BNF growth, of a similar order of magnitude, have been found in cyanobacteria and aerobic heterotrophs, due to large energy expenditure associated with oxygen protection mechanisms [23,27]. The indirect cost for BNF under anaerobic conditions is similarly quite substantial and cannot be due to oxygen, but likely illustrates fundamental indirect costs associated with the synthesis of complex enzymatic machinery required to undertake BNF, such as the higher N content of nitrogenase enzyme versus typical ammonium transporters (∼ 5:1 N:protein BNF/NH4 ), other enzymes required for nitrogenase co-factors maturation and nitrogenase regulation, or the requirement for essential metal uptake (Mo, Fe) [1], which could become more costly under anaerobic and sulfidic conditions due to the lower solubility of some metal chemical forms. The value for Q N,BNF used to estimate these costs, calculated from growth rate and nitrogen specific activity ranging from 0.6 to  1.5 µmol N .ODmL -1 , is likely an underestimate of the true net cellular N quota for N 2 -trophic derived biomass. First, another estimate using accumulated N-fix and total biomass increase during the entire ARA experiment yields a slightly higher value at 1.6 µmol N .ODmL -1 , with a corresponding metabolic cost of 12 mol ATP .mol N-BNF . Second, both estimates for Q N,BNF were obtained under strong N limitation due to the 60% inhibition of N fixation by acetylene, as evidenced by the lower observed growth rates during the experiments (0.018 v. 0.044 hr -1 ). So while our parameters allow good in silico reproduction of the ARA experiments, it is likely that cells decrease their N cellular quota to help cope with nitrogenase interactions with non-N 2 substrate. In addition, our model does not consider fluctuations in biomass per cell or cell size between the different conditions, which could result from a lower N quota (i.e., change in cell. ODmL -1 ] in media. Panel (C) shows BNF contribution to biomass growth, as calculated in [10]. In panel (E), "No ARA" condition refers to growth without acetylene in the headspace. Uncertainty bands (filled and dotted) represent model output for the range of initial biomass amounts used in the laboratory experiments (OD 600 = 0.07-0.09).

Table 3
Comparison of metabolic costs and growth outcomes between ammoniotrophic (NH4) and diazotrophic (BNF) growth of DvH, an anaerobic heterotrophic diazotroph. values). Finally, electron flux to nitrogenase was assumed to be separate between N reduction and AR in a constant sum, but as two binding sites for acetylene exist [33], it is likely that the sum of assumed BNF and ARA rates is not equal to total BNF activity in the absence of acetylene. In addition, the competition between acetylene and nitrogen is dependent on the ratio of the two sub-units of nitrogenase (nitrogenase reductase and dinitrogenase), which could vary in vivo [12]. Hence, it is likely that the actual Q N,BNF is at the higher end of our estimate (< 1 µmol N .ODmL -1 ) and that cost of BNF is < 20 mol ATP .mol N . Notwithstanding these study's limitations, our results still indicate a significant indirect cost of 4-12 additional ATP.mol N -1 for BNF under anaerobic conditions.
Our estimate for Q N,NH4 (2.83 ± 0.04 µmol N .ODmL -1 ) is very close to the value found in the literature (2.95 µmol N .ODmL -1 , [42]). The net cost of ammonium uptake is 3.3 ± 0.3 mol ATP .mol N-NH4 -1 , three times higher than the direct cost of ion transport through the membrane [31]. The mechanism of NH 4 + transport through the cell membrane has been highly debated over the last 20 years [32,34,58,61]. One possibility would be that NH 4 + is transported as a charged ammonium ion requiring the cell to pump out one ion at the cost of 1 ATP to retain the electrochemical potential of the cell. The other possibility is that NH 4 + transport occurs in an electro-neutral form as ammonia (NH 3 ), with possibly no ATP cost [32]. Assuming a cost of 1 ATP for NH 4 + uptake, the remaining indirect metabolic cost for biosynthesis and maintenance of ammoniotrophic cells would be between 2 and 3 ATP per mol of biomass N, remaining lower than our lowest estimate for the indirect cost of BNF (4 ATP.mol N -1 ).
We attribute the two-fold difference between the ratio of ATP cost and observed growth rate for the two acquisition strategies (mol ATP .mol N,BNF /mol ATP .mol N,NH4 = 7 vs. µ NH4 / µ BNF = 2.7, Table 3) to the difference in cellular N content between the two strategies (Q N =0.9 vs 2.8 mol N .ODmL -1 for BNF and NH 4 , respectively, ratio NH 4 / BNF = 3, Table 3). This stochiometric flexibility, in which BNF has a lower Q N (most likely between 1 and 2 mol N . ODmL -1 ) would also explain the relatively higher maximal biomass yield ratio of diazotrophy vs. ammoniotrophy when compared to the ratio of growth rate and overall cost of BNF (Table 3). A calibration factor of 5.5 x10 8 cell.ODmL -1 (as estimated at 600 nm for Desulfovibrio desulfuricans, [59]) would yield values of 2 and 5 fmol N .cell -1 for N 2 -trophy and NH 4 -trophy, respectively. These values are in the range of Q N found in the literature for pure culture non-diazotrophic organisms (2-7 fmol N .cell -1 , [19,56]), but are 3-4 times lower than values found for another model N fixer, Azotobacter vinelandii, grown under aerobic diazotrophic conditions (∼17 fmol N .cell -1 , [1]). However, we note that Azotobacter vinelandii growth rate is > 5 times faster than DvH under similar temperature conditions (0.04 vs. 0.2 hr -1 at 30 and 25 o C, respectively, [1]). Hence, the high cellular quota of A. vinelandii can be attributed to the higher energy yield from oxygenic respiration (4-20 time more energy than sulfate respiration). This indicates that oxygenic respiration could allow A. vinelandii, an organism of similar size and shape to DvH, to afford higher N biomass content, even considering the large amount of energy spent in O 2 protection strategies, such as wasteful respiration or production of a low permeability alginate layer [27,46].

Model performance and limitations
Our model is able to resolve the main qualitative and quantitative features of heterotrophic anaerobic growth (e.g., growth rates, biomass yield, ammonium consumption) under a wide range of [NH 4 + ] init conditions using only a single set of initial parameters. The DAD-GM adds to the body of diazotrophic models that cover a variety of metabolisms, from aerobic heterotrophs to phototrophs [29,30,28,45]. One of the key simplifying features of the DAD-GM is that its parameterization relies on activity measurements (e.g., nitrogen uptake or nitrogen fixation, sulfide production), which are easier to obtain than growth rate in environmental samples and would allow the model framework to be adapted to function in more chemically and physically complex conditions (e.g., sediment, soil cores).
The primary disparity between our modeled and measured data occurs in the diazotrophic regime of the model, particularly when reproducing growth curve dynamics. Modeled diazotrophic growth is generally faster than observations (see Fig. 4 A, D, and G), which show delayed or flattened growth. In addition, growth rates, biomass density yields, and decay rates in laboratory experiments were generally more variable during the diazotrophic phase than for ammoniotrophy (see standard deviations in Table 2, Fig. 3B, and Fig. 4 A,C,D, and F). There are several explanations that could contribute to these phenomena. Firstly, while all cells would be able to use ammonium at the beginning of the experiment, the initial amount of active nitrogenase enzymes in cells could fluctuate. Our measurements suggest there could be generally less of an initial active nitrogenase pool than what we arbitrary assumed for the model (Cell BNF,0 = 0.002). Thus, a better temporal fit between experimental data and model output could be obtained by adjusting Cell BNF,0 parameters for each experiment, or even for each replicate. Other explanations could stem from nitrogenases extreme sensitivity to oxygen [20] or possible competition between sulfate and oxygen as final electron acceptor [25], both of which could slow down some cultures and result in lower maximal biomass yields due to the associated metabolic cost of protective measures [27] (see SI Supplementary Material Fig. S1).
There is evidence of decoupling between N acquisition and biomass build up, as seen during the 5-10 hrs delay between measured and modelled [NH 4 + ] uptake following addition that does not have much effect on the biomass growth (compare Fig. 5 E and F). This could be explained by the existence of a legacy intracellular N pool, which would support current growth. This phenomenon would warrant further characterization before being incorporated into the model.

Nature of the putative BNF-limiting phenomenon and its environmental significance
The plateau of maximum biomass density at OD 600 = 0.38-0.4 obtained for [NH 4 + ] init between 0 and 1000 µM was an interesting and unexpected feature of our experimental results (Fig. 3B). Specifically, at [NH 4 + ] init = 1000 µM, and using the 2-fold measured difference in maximal biomass yield between ammoniotrophy and diazotrophy as found under purely diazotrophic and ammoniotrophic conditions (0.38 and 0.83 OD 600 , respectively, Table 2), a biomass density of OD 600 ∼ 0.6 would be expected. However, the consistent findings of a lower OD 600 yield indicate complexities associated with the dual growth mode induced by sub-replete NH 4 + conditions. Notably, this biomass plateau occured at the same time as a decrease in maximal H 2 S tot production for [NH 4 + ] init between 300 and 2000 µM (Fig. 3C), reaching its lowest value for [NH 4 + ] init = 1000 µM, indicating that not all pyruvate and sulfate were used under these conditions. We changed model parameters to better understand the processes involved in the yield plateau at sub-replete NH 4 + . While we were able to reproduce the maximal biomass yield plateau by changing the magnitude of biomass decay, it was not possible to reproduce the measured decreases in maximal H 2 S tot concentration without the introduction of an arbitrary imposed limiting step specific to BNF that we only allowed to operate after a certain amount of biomass had been produced, or a certain amount of the existing biomass had died (compare SI Supplementary Material Fig. S2, Fig. S3, and Fig. 3C). While a change in the Pyr: SO 4 2ratio could have explained our results, we found a significant decrease in acetate, yielded from the incomplete oxydation of pyruvate, between [NH 4 + ] init = 0 and 1500 µM and a significant increase between 1500 and 3000 µM (SI Supplementary Material Fig S4). This indicates that stoichiometric flexibility, if present, could not account for all the observed features of our data (i.e., biomass yield plateau and maximal [H 2 S] tot decrease), and points to the existence of a phenomenon that would specifically prevent diazotrophic growth in batch culture. ] init of 0 and 1000 µM. We tested similarly a two-fold increase in Ni, required for hydrogenase, and the removal of the initial background of 2% v/v H 2 , a direct inhibitor of BNF, using 100% N 2 (SI Supplementary Material Fig. S5). In another experiment, we tested the addition of an initial H 2 S concentration of 4 mM (as Na 2 S), a strong metal complexing agent under anaerobic conditions [54], and the substitution of Fe +2 for a chelated and presumably less prone to precipitate form of Fe at the same concentration (Fe-EDTA, 2.5 µM) (SI Supplementary Material Fig. S6). We found no improvements in maximal measured yields or in growth rates at [NH 4 ] init = 0 µM, indicating neither direct metal availabilities nor the high concentration of H 2 S found at the end of growth or the inhibition of nitrogenase by H 2 were obviously involved in this BNF limitation. In addition, in case of nutrient limitation, metabolic flexibility would generally allow growth to continue at a slower rate in a linear fashion (i.e., each subsequent doubling of cell halves the intracellular quantity of the limiting nutrient and reduces instantaneous growth rate), which we did not observed (Fig. 4). Toxicity from H 2 S could also be ruled out, as diazotrophy was able to support biomass production in N depleted cultures up to the stochiometric 7.5 mM H 2 S required to fully consume the initial 30 mM of pyruvate, and H 2 S produced at [NH 4 + ] init of 1000 µM was lower than this maximal (4 mM, Fig. 3C). Other potential explanations include continuous, slow diazotrophic growth or a long lag-phase overlaid by a faster decay phase. However, sequential measurement of H 2 S for [NH 4 + ] init of 1000 µM up to 100 hrs after the culture reached its maximal biomass yield showed very little increase in H 2 S production, and the maximal measured H 2 S value was used in Fig. 3C. As we were able to model the dip in H 2 S by imposing a constraint on total decayed biomass, we suggest that this phenomenon is either linked to the direct production of an inhibiting substance by DvH, or the indirect interaction of some substances (e.g., metals, metabolites) with unknown forms of organic matter [26,9]. Further experiments are required to decipher the exact nature of BNF-specific limitation under the tested culture conditions.

Temporal and spatial aspects of nitrogen growth modes
The long transition time between the end of ammoniotrophic growth for ammonium adapted cells and the start of an exponential diazotrophic phase contrasts with the extremely fast resumption of BNF activity in BNF-adapted cultures upon depletion of any added NH 4 + , both in pure cultures and in benthic sediment slurries [11] These results indicate that the initial transition from ammoniotrophy to diazotrophy occurs under extreme N limitation, in which the build up of the nitrogenase enzyme pool could be itself strongly limited by the acquisition of new N. This interpretation is further supported by the rather sharp termination of the ammoniotrophic growth phase, without much of a slow-down before the complete exhaustion of NH 4 + (see Fig. 4), the fast ammonium uptake following addition (see Fig. 5F), and the generally high affinity of the NH 4 + transporter in the literature [35] and required in our model to fit our data (K M,NH4 < 20 µM). Accordingly, this however suggests that DvH is unable to sense ammonium depletion to onset nitrogenase synthesis prior to reaching extremely low [NH 4 + ] (< 1 µM). Under extreme N deplete condition, the second diazotrophic exponential phase can only be explained as either (i) the exponential self-reproduction from newly-fixed N of the initial pool of nitrogenase (i.e., from the inoculation) carried over by each cell following growth dilution during the ammoniotrophic phase, or (ii) by the replication of the first cells able to synthesize new nitrogenase enzyme (i.e., cell specialization), which would then have a strong competitive advantage to take over the population of cells. Thus, cells that could retain a pool of nitrogenase enzymes in a fluctuating environment in which [NH 4 + ] remains low (0-20 µM) can strongly benefit from their initial resource investment. Such a phenomenon likely contributes to results from phenotypic heterogeneity experiments carried on with N fixers, which shows that single-cell BNF rates prior to N depletion correlates with their growth rates post limitation [49]. Conversely, constant exposure to high N loadings followed by N depletion would prevent DvH from quickly shifting to diazotrophic growth in natural ecosystems.
The above considerations support segregating the total cells into ammoniotrophic and diazotrophic-competent cells in the model to efficiently reproduce macroscopic population level results. Indeed, for a unit of biomass created from both NH 4 + uptake and BNF, it is not possible at the timescale at which experiments take place to macroscopically distinguish the biomass being created sequentially (i.e., first NH 4 + uptake, then BNF), simultaneously in two different cells, or simultaneously in the same cells. This simple model would be a good candidate to include in more complex large-scale biogeochemical models, and is easily adaptable for spatially explicit models to help further our understanding of aspects related to marine biogeochemistry (e.g., particulate biogeochemistry and niche expansion, [6,39]) or to microbial ecology (e.g., metabolic heterogeneity within a population at the single cell level, [49,62]).

Conclusion
The Diazotrophic-Ammoniotrophic Dual Growth Model was able to accurately reproduce the growth dynamics of the anaerobic heterotroph sulfate reducer Desulfovibrio vulgaris Hildenbourough and its response to ammonium addition with two different substrates (pyruvate and lactate), using the same single set of parameters, and the additional adjustment of one parameter related to a diazotrophic-specific growth limitation and of the substrate sulfate usage stoichiometry. Our model was able to accurately reproduce the data of laboratory experiments, from growth rate, maximal biomass yield, sulfate consumption, and growth dynamics, as well as headspace ethylene increase and ammonium consumption during ARA experiments. In the future, the DAD-GM, which is built on simple metabolic relationships for ammoniotrophy and diazotrophy at the cellular level, represents an ideal model to study the behavior of N fixers under fluctuating environments in which ammonium may not be fully replete (e.g., due to diffusion of ammonium in sediment, or low frequency of fresh media supply in chemostat), and could easily be adapted to a spatially explicit or agent based models to investigate the microscale phenomena in physically and chemically complex environments. Additionally, our results demonstrate the existence of a putative limiting phenomenon that selectively hinders the activity of nitrogenase but not ammonium uptake under laboratory conditions. Should a similar phenomenon occur in nature, the potential of benthic ecosystems to obtain new nitrogen under N limitation could be reduced.

Declaration of Competing Interest
The authors declare no competitive interests.