A low supply of biologically available nitrogen (N) often constrains primary productivity across a wide range of ecosystems. This constraint—termed N limitation—is demonstrated by the fact that additions of N often increase rates of productivity, both alone and interactively with other elements. Moreover, N additions enhance productivity in more terrestrial ecosystems and do so to a greater extent than do additions of any other nutrient (Elser et al. 2007; LeBauer and Treseder 2008; Fay et al. 2015; Hogberg et al. 2017). The widespread importance of N limitation begs the question “What controls the availability and cycling of N?”

This is not a question that answers itself, because absent supplementation by fertilization, anthropogenic atmospheric deposition, or similar mostly human-driven processes, N supply does not exert independent control over terrestrial ecosystems. Rather, the pools and fluxes of N are consequences of multiple other ecosystem processes operating on multiple time scales and hierarchical levels of control. These controls range from proximate ones that may themselves be influenced by other controlling processes, to ultimate, independent controls that may be remote from immediate influence on N cycling. Following the framework of Jenny (1980), terrestrial ecosystems can be conceptualized as being controlled (ultimately) by a few independent “state factors”, the most important of which Jenny identified as climate, relief or topography, parent material, organisms (defined as the species pool that could be present in a site), time, and human activity. In this framework, state factors ultimately control the development and properties of soils and the composition and functioning of biological communities (the species actually present on a site, and their dynamics). The absence of N from this list of controls (except as a consequence of human activity) is striking, given the widespread occurrence of N limitation.

The existence and importance of N limitation has long been a mystery in biogeochemistry, in both marine and especially terrestrial ecosystems (Howarth 1988; Vitousek and Howarth 1991; Howarth and Marino 2006). We focus on terrestrial ecosystems here, because organisms with the capacity to fix N biologically are ubiquitous there, and symbiotic N fixers in particular have the capacity to add enough N to quickly overcome any deficit in N (as well as possessing what appears to be a substantial competitive advantage in N-limited ecosystems). Publications that demonstrate that N is limiting, and even those that show it is not limiting in particular situations, far outnumber those that ask how N can be limiting. Two exceptions are publications by Vitousek and Howarth (1991) and Menge (2011), and other publications by these authors. These publications suggest that the co-occurrence of two conditions is necessary for the existence of N limitation that is more than marginal and/or ephemeral. First, there must be constraints to symbiotic N fixation that keep it from responding fully to N deficiency. Second, there must be losses of N from the ecosystem that cannot be prevented by organisms within an ecosystem, even where the growth of those organisms normally is limited by N. The first condition is widely appreciated (Vitousek and Howarth 1991; Menge et al. 2008); the second condition is equally important because even low rates of N inputs (such as atmospheric deposition of fixed N in unpolluted regions) would otherwise accumulate, eventually to the point that even the slow turnover of a large pool of N would suffice to meet the needs of organisms and thus to reverse N limitation. Menge et al. (2008) pointed out this second condition in their analysis of constraints to symbiotic N fixation, reporting that their model showed that large N losses that could not be controlled by biota would favor the growth of symbiotic N fixers. They emphasized the continuous process of leaching losses of unavailable forms of dissolved organic N (DON) (Hedin et al. 1995). Their focus was on controls of symbiotic N fixation, while ours is on N limitation, so we evaluate the co-occurrence of uncontrollable N losses (which we define as either losses of N in forms that cannot be acquired by N-limited organisms or losses of forms of N that are biologically available at times and/or places where those forms are present in excess of organism’s requirements) and constraints to N fixation (as did Menge 2011), in our case with a focus on losses of N and their controls. We show here that temporal variability in rainfall, as is expected (and observed) to increase with ongoing greenhouse-gas-driven climate change, can also cause uncontrollable losses of N and so could drive the development of N limitation (where there are sufficient constraints to symbiotic N fixation). In addition, losses of N caused by disturbance (including harvest, fire, and windthrow) are similarly uncontrollable and so can drive the development of N limitation (Vitousek 2003).

We use a simple “toy” model of N dynamics, one that includes symbiotic biological N fixation (BNF), implemented along a climate gradient to evaluate the nature and causes of N limitation to primary productivity. A “toy” model is defined as a simplified set of objects and equations that can be used to understand a mechanism. Although the model we use is relatively complex, we call it a toy model because its purpose is to explore a general principle rather than to capture the quantitative details of a particular system. In the sense of Rastetter (2017), it is a model designed to further understanding rather than a model designed to yield numbers. We find that developing a model like this is useful for testing our understanding of the implications of biogeochemical processes, because such a model forces us to write equations that express our understanding of ecosystem processes, and running the model shows whether our understanding (when pulled together in the toy model) actually yields anything like the patterns we observe in the field. This model was developed for a climate gradient in Hawaiʻi; the model and a sensitivity analysis of it and the model’s output compared to the pattern of 15N natural abundance and of a peak in biologically available soil phosphorus (P) that is caused by biological uplift in the mesic portion of that gradient are summarized in an excerpt from Vitousek et al. (2021) associated with this publication as on-line supplementary material. The most recent use of the model (Vitousek et al. 2021) sought to understand the consequences of biogeochemical processes on this particular gradient. Here we use this “toy” model to evaluate the general question (though grounded in a particular place) of what permits the development of N limitation to primary productivity, especially as an aid to understanding the pathways and implications of losses of N that cannot be prevented by N-limited organisms.

Methods

Toy Model. We wrote the program that runs the toy model in Matlab version 2018a (Mathworks Inc.). We used three versions of the basic model here, all of which are provided in supplementary material. The three are newprodall, which was close to the complete model reported in Vitousek et al. (2021) (with the modifications described below); newprodnofixcr, which is a simplified and readily modified version of the basic model without deep soils and with an easy way to turn off BNF, this model is run across the full climate gradient, as is newprodall; and tcnewprodnofixcr, an equally simplified version of the model designed to evaluate ecosystem dynamics over time at a single level of mean rainfall. For model runs across the gradient, we used the programs newprodall and newprodnofixcr, and looped through the levels of mean precipitation reported for that gradient (from 2 to 27 cm/month, from Giambelluca et al. 2013]), and then executed an inner loop in time (from 2 to 60,000 months) at each level of mean precipitation. For some model runs, we used tcnewprodnofixcr at a selected level of mean precipitation at which N limits primary productivity and the accumulation of plant biomass (in the model), and ran the model over time at that selected mean precipitation. The complete program newprodall is built upon earlier versions of the toy model that evaluated controls of symbiotic BNF (Vitousek and Field 1999, 2001), that were then modified to evaluate the weathering of soil minerals in Vitousek et al. (2016), and modified again to include a deep soil layer (conceptualized as 50–200 cm deep) and to apply to a climate gradient in Vitousek et al. (2021). The structure of the program used here was like that in Vitousek et al. (2021) with a few modifications, the most important of which was that we changed the cost of N acquisition by a non N fixer to be a consistent function of the quantity of biologically available N in the soil (from all sources, including N carried over from previous times, N mineralization, N deposition, and sometimes N fertilization) rather than being a function of the fraction of the biologically available N that was present in the soil at the beginning of that time step (Fig. 1). This change is logical and likely more realistic, and simulations across the climate gradient demonstrated that it yielded patterns similar to those in the older versions. Other changes from the program used in Vitousek et al. (2021) are different background levels of temporal variability in precipitation, designed to reproduce observations of temporal variation in precipitation in a 30+ year observational record (P. von Holt, personal communication) collected by Ponoholo Ranch (adjacent to the climate gradient). Also, we returned to the formulation used in Vitousek and Field (1999) for the effects of shading by a non-N-fixing plant on the productivity of an N-fixing plant, and (like the 1999 paper) we included the possibility of differentially greater grazing on the N- and P-rich N fixing plant. Finally, we corrected an error in the formulation of P mineralization which had had the effect of creating unforced temporal variability in P dynamics, we increased the coefficient for ammonia volatilization from a maximum of 10% of the pool per month to 20% of the pool, we decreased the coefficients for maximum DON and DOP leaching from 0.5% to 0.05% of their organic pools in soil (respectively) to 0.2% and 0.01%, and we decreased the coefficient for leaching of phosphate (PO4) in soil from a maximum of 10% of the pool per month to 1%; these changes are logical and consistent with observations in Hawaiian ecosystems (Neff et al. 2000) and with the greater mobility of N than of P in soils generally.

Fig. 1
figure 1

Conceptual model for the costs of N acquisition (characterized as the proportion of plant NPP spent on acquiring N) by a non N fixer and so for the transition to production by an N fixer. which we assume occurs when the cost of N acquisition by the non N fixer exceeds the (constant) cost of N acquisition by the N fixer. a The formulation used originally (from Vitousek and Field 1999) and in subsequent iterations of the model (including Vitousek et al. 2021). Here the x-axis is the N in available forms in the soil; it becomes costlier for the non N fixer to acquire N when half of that original N has been acquired. When the cost of N acquisition by the non N fixer exceeds the cost of N fixation by the N fixer, the N fixer alone grows. With this formulation, the supply of N in the ecosystem will equilibrate over time to point B. b The formulation used in this paper. Here the x axis is N supply in the soil (from all sources, including carryover of N from earlier time steps, N mineralization, atmospheric deposition, and fertilization), and the cost of N acquisition by the non N fixer increases steeply at a non fixation threshold (nonfixthresh) (set here to 2% of the N supply required for maximum production by the non N fixer [MaxprodN]/50). With this formulation, the N fixer will grow and add N to the system as long as the supply of N is less than MaxprodN + nonfixthresh, so the supply of N in the ecosystem will increase to that point (and decrease above that point)

Full size image

The model program is described in detail in Vitousek et al. (2021) and in the excerpt from that paper included here as on-line supplementary material; the structure of the model is summarized as Fig. 2. Briefly, we simulate productivity and biomass of a non N fixing plant species and a symbiotic N fixing perennial plant species. We did not simulate different plant tissues such as leaves, roots, or woody tissue, and so we do not specify the type of ecosystem simulated (other than that the vegetation is perennial). Biomass turnover in the model is constant, so (at equilibrium) biomass corresponds to time-averaged net primary production. The non N fixer and the N fixer differ in their stoichiometry, with the N fixer richer in N and P relative to C, in comparison with the non N fixer. The model embodies the assumption that the N fixer fixes N obligately (Menge et al. 2009, 2015) and derives all of its N from fixation. Where there is N available for uptake, the non N fixer obtains it more cheaply (energetically) than the fixer can fix N, so the non N fixer has absolute priority to light, water, and P—as long as N is available to it. The N fixer will only grow (and fix N) if the non N fixer runs out of N and there is still light, water, and P available to support the growth of the N fixer. Even with this strong constraint to growth and N fixation by the N fixer, we find that N limitation is marginal and ephemeral—unless there are additional constraints imposed by features such as the P-rich stoichiometry of the N fixer, differential effects of shading that reduce growth of the N fixer, and/or preferential grazing on the N- and P-rich N fixer (Vitousek and Field 1999).

Fig. 2
figure 2

The structure of the complete model in the program newprodall, from Vitousek et al. (2021)

Full size image

For the analysis of N limitation, other important features of the toy model are its treatment of temporal variability in precipitation and of disturbance. In the absence of such temporal variability or disturbance, the model quickly comes to equilibrium, with a strong threshold in simulated soil and ecosystem properties where precipitation equals evapotranspiration; with temporal variability in precipitation included, this threshold is less sharp, and it migrates over time to progressively drier positions along a precipitation gradient over time, in a reasonable pattern (Vitousek et al. 2016, 2021). Precipitation variability is particularly important in that increasing precipitation variability (together with the related occurrence of extreme precipitation events—floods and droughts) is one of the most robust predictions (and observations) for consequences of greenhouse-gas-driven anthropogenic climate change (Kharin et al. 2007; Pendergrass et al. 2017). In the toy model, we represent precipitation and its variability with a random component to month-to-month variation. For the gradient, we simulate precipitation with a mean level for that position on the gradient; variability is calculated as:

Simulated Precipitation (for that time step, assumed and adjusted to be a month) = [Mean Monthly Precipitation (for that position on the gradient)] + [Mean Monthly Precipitation * (a constant > 1) * (rand (a random number evenly distributed from 0 to 1) – (a constant > 0.5 and < 1)]. If the result of this calculation is a negative number, we set precipitation to 0 for that month.

The constants were selected (and varied) to yield a progressive increase in the variability of monthly precipitation in progressively drier sites, as is observed in practice at nearby Ponoholo Ranch. We chose constants that would yield a range in precipitation variability with the appropriate mean precipitation for that position on the gradient. Modeled and observed coefficients of variation (standard deviation/mean) for precipitation across the climate gradient are shown in Fig. 3.

Fig. 3
figure 3

Observed (solid symbols) versus simulated (hollow symbols) Coefficient of Variation (CV) (CV = standard deviation/mean) for precipitation on a climate gradient in Leeward Kohala, Hawaiʻi. Both observed and modeled mean precipitation was taken from values for those sites from the rainfall atlas of Giambelluca et al. (2013); the precipitation that we measured was lower (Marshall et al. 2017), but (for comparison with other work) we used the atlas values for both observed and modelled precipitation here. Observed CVs from a 32 year record of measurements on Ponoholo Ranch (Pono von Holt, personal communication), parallel to the climate gradient

Full size image

To model disturbance, we considered one frequency of disturbance (every 500 months) and three types of disturbance that loosely correspond to harvest, fire, and windthrow. To model harvest, we removed all of the biomass (and the nutrients it contained) after 45,000 months and every 500 months thereafter. For fire, we used the same timing and assumed that all of the C and N contained in biomass was lost to the ecosystem (volatilized) with each fire, while all of the P was transferred to the mineral P pool in the soil. This assumption exaggerates the difference between N and P in responses to fire (Ewel et al. 1981), but it does so in the right direction. Finally, to simulate windthrow we again used the same timing, and assumed that all of the biomass and the N and P it contained were added to soil organic pools with each disturbance Simulated windthrow is particularly interesting in that the disturbance itself does not involve any loss of nutrients, though there are losses of N and P following disturbance.

For analyses across a precipitation gradient, we used an outer loop in mean monthly precipitation and an inner loop in time (hereafter we refer to the time steps as months). We started each rainfall level with no time step-to-time step variation in precipitation to allow rapid equilibration, then (where such variation was included) we introduced random variation (as described above) after 12,000 months. When we introduced added N (or P) to evaluate nutrient limitation, we did so after 45,000 months. In any case, we averaged (and report) results for the last 10,000 months of simulation. For analyses at a particular rainfall, we followed the same procedure, generally reporting the time course of responses for the last 10,000 months. These time steps and run lengths were selected as a compromise between the simulated development of features of interest and the time it took to run the program.

Observations. We compared toy model outputs with observations from a precipitation gradient in Kohala, Hawaiʻi (Vitousek et al. 2021). Measures of most soil properties on that gradient are summarized in Vitousek and Chadwick (2013); N dynamics, including δ15N, are summarized in von Sperber et al. (2017) and Burnett et al. (in press). We assumed that most losses of nitrate are via leaching rather than denitrification, consistent with the nitrate 15 N results in von Sperber et al. (2017) (and discussed in on-line supplementary material)—but for the development of N limitation under temporal variation in precipitation the accumulation of nitrate is most important, not the pathway of losses of that nitrate. For precipitation and its variability, we drew upon 32 years of observations across a parallel gradient on neighboring Ponoholo Ranch (von Holt, personal communication).

Results and discussion

Simulated N and P limitation to the NPP of a non N fixer across the climate gradient, evaluated for the full model using the program newprodall and including BNF and a deep soil layer, is shown in Fig. 4. The productivity of the simulated non N fixer is limited by both N and P between ~ 8 and ~ 19 cm/month and above ~ 25 cm/month of precipitation, and limited by N (alone) between ~ 19 and ~ 25 cm/month of simulated precipitation. (P limitation to simulated NPP occurs over less of the gradient than is reported in Vitousek et al. 2021) because of the changes described above in modeling the P cycle). The maximum extent of N limitation (maximum difference between NPP of the non N fixer with versus without added N) occurs at 20 cm/month of simulated precipitation. Since it is long-term response to simulated fertilization that is reported here (fertilization with 10 units of N or P per month − enough to ensure that neither N nor P limited NPP − began after 45,000 months, and results are averaged from 50,000 to 60,000 months), it is possible that where both N and P limit the growth of the non N fixer that added P stimulates the growth of the N fixer which then adds N to the system, enhancing growth of the non N fixer.

Fig. 4
figure 4

Effects of simulated additions of N and P on the simulated NPP (in arbitrary units that are close to kg/ha) of a non N fixer across the climate gradient, with all controllable and uncontrollable losses of nitrogen active and with Biological N Fixation (BNF) included, from the program newprodall

Full size image

We then used the related program newprodnofixcr to evaluate nutrient limitation both with and without N fixation, with no losses of dissolved organic nutrients and no temporal variation in simulated precipitation or disturbance (Fig. 5). In this situation (without BNF or any uncontrollable losses of nutrient elements, but with N losses as nitrate or ammonia), neither N nor P limited production of the non N fixer above 9 cm/month of simulated precipitation because weathering of P-containing minerals and atmospheric deposition of N accumulated sufficient nutrients to offset limitation. Below 10 cm/month of simulated precipitation, N supply limited NPP of the non N fixer because insufficient N accumulated from atmospheric N deposition to offset limitation, in a long-term transient phenomenon. We demonstrated this process of accumulation by running the model (here tcnewprodnofixcr) without BNF or uncontrollable N losses for 120,000 months rather than the standard 60,000 months, with simulated N fertilization after 105,000 months, and results averaged from 110,000 to 120,000 months (Fig. S1). Under these conditions, enough N accumulated to offset N limitation at 8 (not shown) and 9 cm/month of simulated precipitation (Fig. S1), as well as at 12 cm/month (as had been true after 60,000 months—Fig. 5).

Fig. 5
figure 5

Simulated N limitation to simulated NPP across the climate gradient with no uncontrollable losses of N, from the program newprodnofixcr. a No BNF. NPP of a non N fixer without (solid black line) and with (thick gray line) simulated N fertilization. b BNF active. NPP of a non N fixer (solid black line, hidden behind the summed line) and an N fixer (dashed black line), of an N-fertilized non N fixer without BNF (dotted line), and the sum of NPP by a non N fixer plus 1.5 times NPP of an N fixer. The factor of 1.5 represents the greater concentration of N in the N fixer, and so it puts the NPP of both plants on a comparable basis. This sum is very slightly greater than is NPP by the N-fertilized non N fixer alone, demonstrating that where BNF is active it can prevent N limitation under these conditions. As Fig. S1 demonstrates, the N limitation observed in at least some of the low-rainfall sites here is a transient phenomenon

Full size image

We then evaluated each of the potential pathways of biologically uncontrollable N losses that we identified (DON leaching, temporal variability in precipitation, ecosystem-level disturbance) by adding these processes into the programs newprodnofixcr and/or tcnewprodnofixcr individually, with and without BNF, and determining where and how each process could drive N losses that organisms could not control. We evaluated losses of DON first (Fig. 6a, b). Without temporal variation in precipitation, there is no leaching loss of water or nutrients at or below the point where potential evapotranspiration (PET) equals precipitation (PPT) (here 15 cm/month of simulated PPT); leaching losses of water (and so DON losses) increase progressively where PPT exceeds PET in wetter sites. Accordingly, in the absence of BNF, NPP of the non-fixer becomes progressively (and profoundly) more N-limited in wetter portions of the gradient, as is demonstrated by the response to simulated N fertilization there (Fig. 6a). With BNF active (Fig. 6b), NPP of the N fixer increases in wetter portions of the gradient, and the N that the fixer adds to the system maintains higher NPP than occurs in the absence of BNF (Fig. 6a). BNF also adds N that prevents N limitation in drier sites (except for the driest sites, where BNF is suppressed within the program as described in Vitousek et al. (2021) and the excerpt from that paper on-line to match the low abundance of symbiotic N fixers observed in these dry sites (Fig. 6a,b). The ongoing uncontrollable losses of DON drive ongoing N limitation to the NPP of the non N fixer (and so ongoing BNF by the N fixer) in wetter sites; we wonder if one of the reasons for the widespread abundance of symbiotic N fixing trees in wet tropical forests (Hedin et al. 2009) is the high rainfall in many of those systems, with its attendant DON losses. With BNF active (Fig. 6b), the sum of NPP by a non N fixer plus 1.5 times the NPP of an N fixer [1.5 to account for the N-rich stoichiometry of the N fixer—McKey (1994)—and so to put NPP of non N fixers and N fixers on a comparable basis] is very slightly greater than the NPP of an N-fertilized non N fixer, showing that the existence of N limitation here requires barriers to BNF as well as losses of DON.

Fig. 6
figure 6

Consequences for simulated N limitation of including leaching losses of dissolved organic N (DON), as well as losses of nitrate and ammonia. a No BNF. NPP of a non N fixer without (solid black line) and with (thick gray line) simulated N fertilization. b BNF active. NPP of a non N fixer (solid black line) and an N fixer (dashed black line), of an N-fertilized non N fixer without BNF (dotted line), and the sum of NPP by a non N fixer plus 1.5 times NPP of an N fixer. The factor of 1.5 represents the greater concentration of N in the N fixer, and so it puts the NPP of both plants on a comparable basis. This sum is very slightly greater than is NPP by the N-fertilized non N fixer alone, demonstrating that where BNF is active it can prevent N limitation under these conditions

Full size image

Next we evaluated the consequences of temporal variability in precipitation. Vitousek et al. (2021) suggested that variability in precipitation could drive an imbalance between the supply of N to N-limited organisms and demand for that N. Losses of nitrate (normally a form of N for which losses can be controlled by N-limited organisms) could occur by either nitrate leaching or denitrification when the supply of N was in excess of demand, and those losses could mean that there would not be enough biologically available N for maximum growth at other times, such that N could be limiting, if there are also constraints to biological N fixation. This asynchrony between N supply and demand, leading to N losses, had been reported in semi-arid and arid ecosystems by Dijkstra et al. 2012, Nielson and Ball (2015), and others. We modified the program newprodnofixcr to include the level of temporal variation in precipitation modeled (and observed) across the climate gradient (Fig. 3), and evaluated N limitation to NPP caused by precipitation variability with and without BNF, and with only putatively controllable losses of ammonium or nitrate active (Fig. 7a, b). Precipitation variability at the levels observed on the climate gradient caused N limitation to the non N fixer at all levels of precipitation below ~ 22 cm/month of simulated precipitation (not just in semi-arid or arid ecosystems) (Fig. 7a), as demonstrated by the response to simulated N fertilization there; when we made BNF active, it increased the NPP of the non N fixer, but not to the level that simulated N fertilization did (Fig. 7b). Indeed, the NPP of the non N fixer plus 1.5 times that of the N fixer was less than that of the non N fixer with simulated N fertilization (Fig. 7a, b) in the precipitation range from ~ 12 to ~ 22 cm/month.

Fig. 7
figure 7

Consequences of temporal variability in precipitation across the climate gradient with losses of ammonia and nitrate but not losses of DON included. a No BNF. NPP of a non N fixer without (solid black line) and with (thick gray line) simulated N fertilization. b BNF active. NPP of a non N fixer (solid black line) and an N fixer (dashed black line), of an N-fertilized non N fixer without BNF (dotted line), and the sum of NPP by a non N fixer plus 1.5 times NPP of an N fixer. The factor of 1.5 represents the greater concentration of N in the N fixer, and so it puts the NPP of both plants on a comparable basis. This sum is less than is NPP by the N-fertilized non N fixer in the range from ~ 12 to ~ 22 cm/month of simulated precipitation, suggesting that precipitation variability constrains the responsiveness of BNF to N deficiency there

Full size image

We evaluated the potential importance of precipitation variability in driving N limitation further by examining different intensities of precipitation variability at 20 cm/month of precipitation, the point at which N limitation to NPP was maximized (Fig. 4). To carry out this analysis, we ran the program tcnewprodnofixcr at 5 different intensities of temporal variability in precipitation, with the middle intensity corresponding to the observed temporal variation at 20 cm/month of simulated PPT on the climate gradient, two lower levels of temporal variability, including no temporal variability, and two higher levels of temporal variability in precipitation. For each level of variability, we compared the biomass of a non N fixer without added nutrients to its biomass with added N to visualize the extent of N limitation. (We used biomass rather than NPP here because biomass is less variable than NPP and so the effects of temporal variation are more observable, but as described above biomass is a many month integration of NPP.) Results of this analysis are shown in Fig. 8a–e. Nitrogen limitation was absent where there was no temporal variation in precipitation (Fig. 8a), but increasing levels of temporal variation in simulated precipitation in the absence of BNF increased N limitation and decreased the biomass of a non N fixer progressively (Fig. 8b–e).

Fig. 8
figure 8

Consequences of 5 intensities of temporal variability in precipitation at a simulated precipitation of 20 cm/month, the level of simulated precipitation at which N limitation was greatest (Fig. 4). These simulations include losses of ammonia and nitrate but not of DON, with no BNF. Also, the simulated effects of fertilizing with more than enough N to saturate plant demand (thereby reversing any N limitation). Biomass of a non N fixing plant is shown here (solid black line = without simulated N fertilization, thicker gray line = with simulated N fertilization). a No temporal variability in precipitation; there is no effect of simulated N fertilization under these conditions. b Low level of temporal variability in precipitation. Here the Coefficient of Variation (CV) (standard deviation/mean) is 0.29. c Intermediate level of temporal variation in precipitation, equivalent to that observed at 20 cm/month precipitation in Kohala, Hawaiʻi (Fig. 3); CV is 0.76. d High level of temporal variation in precipitation, CV is 1.18. e Very high level of temporal variation in precipitation, CV is 1.64

Full size image

The results in Fig. 8 provide strong evidence that temporal variation in precipitation like that anticipated (and observed) to occur with ongoing anthropogenic climate change can drive biotically uncontrollable losses of N, even in relatively high-rainfall systems, and that it is already doing so at current levels of temporal variation (especially in low-rainfall sites where such variability is greater). We tested the mechanism underlying this result by plotting the demand for N by non N fixers (the amount of N that would support maximum productivity by the non N fixer, given the climate conditions in that time step), the supply of N (the amount of N available for uptake, from all sources), leaching of water, and leaching of nitrate from the system with the observed level of temporal variability in precipitation at 20 cm/month of simulated precipitation, using the program tcnewprodnofixcr (Fig. S2) beginning 12,500 months into the simulation, shortly after temporal variation in precipitation began. Most losses of nitrate occurred when there was leaching loss of water at a time when the supply of N was greater than demand for N. We then evaluated the consequences of these losses of N by turning off both losses of nitrate and ammonia volatilization in the model; the consequence of doing so (in the absence of BNF) was no N limitation to NPP of the non N fixer across the climate gradient (Fig. S3).

Finally (in terms of temporal variability in simulated precipitation), we evaluated the joint effects of temporal variation and losses of DON both with and without BNF, recognizing that temporal variation in precipitation would drive losses of DON at lower rainfall than would be the case without such variation. The results of this analysis are shown in Fig. S4a,b; not surprisingly, without BNF N limitation to the NPP of the non N fixer occurred everywhere on the climate gradient (Fig. S4a), and (given the effect of precipitation variability on the responsiveness of the N fixer to N deficiency discussed above) the sum of NPP by the non N fixer plus 1.5 times the NPP of the N fixer was substantially less than the NPP of the non N fixer with added N (in the absence of BNF) at simulated precipitation levels above ~ 12 cm/month (Fig. S4a,b).

The final pathway of nitrogen losses that we evaluated was ecosystem-level disturbance. Large scale consumptive disturbances like fire and harvest remove nitrogen from ecosystems and provide no pathway by which biota within systems can prevent those losses. These disturbances, and even disturbances that don’t remove nitrogen from ecosystems directly, like windthrow, also disrupt the cycling of N within ecosystems and so drive losses of N in forms like nitrate that normally are available to organisms, with such losses occurring at times when N is superabundant to organisms. Once N is lost at such times, the lower quantity of N within ecosystems that results could reduce N supply and so drive N limitation at times when N is not superabundant.

We evaluated the consequences of disturbance by introducing one of three forms of disturbance (loosely corresponding to harvest, fire, and windthrow, as described above) after 45,000 months of simulation, and repeating that disturbance every 500 months thereafter. This analysis was carried out with losses of dissolved organic forms of nutrients set to 0, with no temporal variation in precipitation, and with ongoing losses via nitrate leaching and ammonia volatilization and with and without BNF, using the program newprodnofixcr, Results of these analyses are shown in Fig. 9 for the situation in which BNF was absent, and in Fig. 10 for the situation with BNF active. N limited the NPP of the non N fixer for every type of disturbance (Fig. 9a–c). However, with BNF active (Fig. 10), it was able to add enough N to offset N limitation to NPP (Fig. 10a, b). The effects of fire and harvest were identical in this analysis, because both were modeled as removing equal amount of N, and P does not limit NPP under these conditions, so we show only fire and windthrow in Fig. 10.

Fig. 9
figure 9

The effects of simulated ecosystem-level disturbance on potential N limitation to a non N fixer, in the absence of BNF, calculated using the program newprodnofixcr. The solid black line shows the NPP of the non N fixer with no nutrients added; the thick gray line shows the consequence of simulated additions of N. a Harvest, which was modeled by removing all of the plant biomass and the N and P it contains from each site. b Windthrow, which was modeled by adding all of the biomass and the N and P it contains to the soil C, N, and P pools. c Fire, which was modeled by removing all of the plant biomass and the N it contains from each site; the P in plant biomass in transferred to the PO4-P pool in the soil

Full size image
Fig. 10
figure 10

The consequences of disturbance for the NPP of a non N fixer and an N fixer (with BNF active) across the climate gradient, calculated using the program newprodnofixcr. NPP of a non N fixer (solid black line) and an N fixer (dashed black line), of an N-fertilized non N fixer without BNF (dotted line), and the sum of NPP by a non N fixer plus 1.5 times NPP of an N fixer. The factor of 1.5 represents the greater concentration of N in the N fixer, and so it puts the NPP of both plants on a comparable basis. This sum is greater than is NPP of the N-fertilized non N fixer in all sites with more than 3 cm/month of precipitation, showing that ecosystem-level disturbance does not constrain the responsiveness of BNF to N deficiency under these conditions. a Fire. b Windthrow. Harvest is not reported here because it yields results identical to those of fire, under these conditions

Full size image

N limitation caused by disturbance occurred in part because N supply can limit NPP in non-equilibrium situations (Fig. 5) and the disturbed ecosystems as modeled are perpetually recovering from disturbance, after the first one, and in part because there are enhanced losses of N from disturbed sites. In this respect, windthrow is more interesting than harvest or fire, in that we modeled no direct removal of nutrients by windthrow, rather the biomass killed by windthrow was added to the soil organic C, N, and P pools and any N losses that occur in this simulation are of potentially biologically available forms after disturbance. To explore controls of N losses with disturbance further, we used the program tcnewprodnofixcr at a simulated precipitation of 20 cm/month, and evaluated the first disturbance imposed upon the system after 45,000 months without disturbance, and a second disturbance that occurred after 55,000 months. This second disturbance was the 11th in the sequence of the disturbances. Results of this analysis are shown for windthrow in Fig. S5a and for fire in Fig. S6a. Losses of nitrate were much greater after the first than the 11th disturbance in the sequence, and the losses of nitrate that occurred following simulated windthrow at simulated precipitation levels greater than 15 cm/month of simulated precipitation (above which leaching losses of water and nitrate occur) were sufficient to induce N limitation in the regrowing ecosystem following windthrow (Fig. S5c). We also simulated NPP and N dynamics following the first disturbance at a simulated precipitation of 10 cm/month (Figs. S5b, S6b), a level at which leaching losses of water and nitrate did not occur. Here (as in wetter sites) simulated windthrow added plant material to the soil with C:N ratios well above the threshold for net N mineralization, so losses of nitrate in wetter sites were delayed following windthrow more than were losses following fire (Figs. S5a, S6a). This effect also contributed to N limitation in sites receiving 10 cm/month of simulated precipitation (Figs. S5b, S6b), and in any site drier than the leaching threshold of 15 cm/month of simulated precipitation. In contrast, nitrate losses following fire at a simulated precipitation of 20 cm/month were less than half of total N losses (including those caused directly by fire), and the ecosystems regrowing following fire were more profoundly limited by N supply than those regrowing following windthrow (Figs. S5c, S6c).

Our analysis demonstrates that there are at least three pathways of losses of N that cannot be prevented even where the growth of plants is limited by N supply - losses of biologically unavailable forms of DON, temporal variability in precipitation that causes losses of forms of N at times when it is in excess of plant’s demand for N, and ecosystem-level disturbances. The effect of each of these alone and in combination is summarized in Fig. S7. However, in our analyses, BNF responded to N losses caused by DON leaching or ecosystem-level disturbance and prevented profound and/or prolonged N limitation; in contrast, temporal variation in precipitation constrained the responsiveness of BNF to N deficiency as well as driving biologically uncontrollable losses of N.

What constraints to BNF could allow N to limit primary productivity where uncontrollable N losses are important? We evaluated this question by removing (one by one and together) the additional constraints to symbiotic N fixers that have been identified (in addition to the absolute priority that the non N fixer in the model enjoys for light, water, and P resources as long as N is available to it—a priority that isn’t sufficient to cause N limitation that is more than marginal and/or ephemeral). These additional constraints include differential growth responses to shade, differentially higher P requirements of symbiotic N fixers, and differentially intense grazing on symbiotic N fixers because of their N- and P-rich stoichiometry (Vitousek and Field 1999). We carried out this analysis for ecosystems influenced by a combination of DON losses and losses associated with temporal variation in precipitation. Although ecosystem-level disturbance is pervasive and represents a major source of uncontrollable N losses across much of Earth (Vitousek 2003) (Fig. S7), these disturbances are idiosyncratic, in that we need to know the history of each place in order to model them appropriately. In contrast losses of DON and losses driven by temporal variation in precipitation are ubiquitous and general.

We carried out our analyses using the complete model, including BNF and deep soils, using the program newprodall. Without any additional constraints to BNF, the sum of NPP by a non N fixer plus 1.5 times the NPP of an N fixer is less than the NPP of the N-fertilized non N fixer (Fig. S8a) indicating that BNF did not respond fully to N deficiency. As was demonstrated above, temporal variability in precipitation has this effect (Fig. 7), in addition to causing uncontrollable losses of N. Including differential effects of shading on the NPP of an N fixer (Fig. S8b) or the effect of a differentially P-rich stoichiometry for the N fixer (Fig. S8c) had no additional effect on this result (not surprisingly for P, because the NPP of an N fixer was not constrained by P supply under these conditions). However, differentially intense grazing on an N- and P-rich N fixer suppressed the growth of both N fixer and non N fixer substantially, and increased the extent of N limitation (Fig. S8d). The pattern for all of these additional constraints to N fixation combined was identical to that for the effects of differentially intense grazing on the N fixer (Fig. S8d, e), demonstrating that in this model under these conditions the dominant constraint to the responsiveness of BNF to N deficiency was differential grazing.

Conclusions

This analysis shows that temporal variation in precipitation like that predicted (and observed) to occur with anthropogenic climate change as well as anthropogenic or natural disturbance can be important vectors of biologically uncontrollable losses of N from ecosystems, and that those losses can in turn reduce the biological availability of N meaningfully. We also show that under the conditions simulated here, differential grazing on symbiotic biological N fixers (caused by their N and P-rich stoichiometry) represents the single most important constraint to the ability of N fixers to respond to N deficiency (Ritchie et al. 1998).