Biological nitrogen fixation (BNF) is one of the most fundamental processes supporting life on Earth. Despite decades of work quantifying BNF across a range of terrestrial systems, biome-scale BNF inputs and their primary drivers remain highly uncertain (Davies-Barnard & Friedlingstein, 2020). Quantifying rates and controls on BNF in natural systems is relevant to several key challenges in ecosystem science. For example, given evolving recognition that N constrains ecosystem responses to increasing atmospheric carbon dioxide (CO 2 ; Wieder et al., 2015), it is increasingly clear that accurate projections of feedbacks to climate change depend on incorporating BNF more accurately into Earth system models (Davies- Barnard et al., 2020;Wieder et al., 2015). Constraining natural BNF rates is also required to estimate the extent of human perturbation to the N cycle (Vitousek et al., 1997), the magnitude of which has been continually revised upwards as global estimates of natural BNF have generally decreased in recent years (Davies-Barnard & Friedlingstein, 2020;Staccone et al., 2020;Sullivan et al., 2014;Viousek et al., 2013).
Finally, via its influence on ecological processes such as competition and facilitation, BNF has important implications for understanding and forecasting ecosystem succession and recovery from disturbance (Batterman et al., 2013;Menge & Hedin, 2009).
Although measurement challenges remain a significant barrier to studying BNF, especially in natural ecosystems, the technical advantages and disadvantages of available methods are generally well understood and have been extensively reviewed (see Section 3.1).
However, measurement methods are only one component of study design. Less systematic attention has been given to the design of protocols that capture the significant temporal and spatial variation in BNF across different ecosystem niches and, if desired, how to scale from point measurements to larger scales using appropriate statistical and error propagation approaches. In surveying the contemporary BNF literature, is it apparent that approaches are not standardized and virtually all studies involve many unexamined assumptions.
This paper outlines a conceptual road map for the design, execution and scaling of terrestrial BNF studies based on field-collected samples, although much advice also applies to more controlled (pot or mesocosm) experiments (Figure 1). Some topics are specific to generating large scale annual field rates (e.g. kg N ha -1 year -1 ), but many are relevant to more limited spatial/temporal ranges or comparative studies focused on elucidating factors regulating BNF (e.g. manipulative experiments). We also make recommendations for data standardization and reporting to support synthesis efforts aimed at defining broad controls on BNF across systems. Despite progress in recent years, many aspects of estimating BNF still present significant logistical and technical challenges and unknowns. We highlight these not to dissuade, but rather to acknowledge that defining the limits of data interpretation can help researchers think critically about which questions and conclusions their data can and cannot address.
The ability to fix atmospheric N 2 is scattered throughout organisms in the bacterial and archaeal domains with a diverse range of metabolic strategies and habitat preferences (Raymond et al., 2004). In addition to symbiotic (rhizobial and actinorhizal) associations with vascular plants, BNF occurs in many niches (ecosystem compartments) even within a single biome including bulk soil, plant rhizospheres, decaying wood and roots, leaf litter, biological soil crusts, inside/on the surface of roots or above-ground plant tissues, in association with lichens and bryophytes, in the colonies and guts of insects and animals, and on snow and ice, and new examples emerge regularly (Cleveland et al., in review;Reed et al., 2011). While symbiotic inputs can be extremely high in certain cases, free-living niches can represent the dominant BNF input in ecosystems as diverse as tropical forests, boreal forests and pine savannas (DeLuca et al., 2002;Reis et al., 2020;Tierney et al., 2019). Spatially, the distribution of niches spans a gamut from highly discrete (symbiotic root nodules and lichens) to continuous (bulk soil), superimposed with vertical gradients such as soil depth or canopy position. For studies that sample continuous niches, it is thus important to clearly define boundaries on the system (such as soil horizon) prior to sampling.
The ubiquity and diversity of BNF means that the choice of clearly defined niche(s) and their relevant spatial extent are important considerations for study design (Figure 1). The choice of measurement niche for a given ecosystem often reflects historical precedence (especially where appreciable rates have been identified before) and methodological constraints, creating systematic biases in our understanding. These include coverage gaps for niches with very low mass-based rates (particularly where detection limits may lead to false negatives), even if those niches are relatively abundant.
BNF by canopy leaf epiphylls in broadleaf forests is one example;
although rarely comprehensively measured it may potentially represent an important BNF source due to its potentially large spatial extent (Freiberg, 1998;Reed et al., 2008Reed et al., , 2011)). Knowing where BNF is not occurring can be just as important as knowing where it is. However, publishing undetected rates is challenging, and this data gap has implications for accurate quantification/scaling of BNF and duplicated sampling efforts. It is still relatively common for the sum
asymbiotic, biological nitrogen fixation, free-living, nitrogen cycle, spatial variation, symbiotic, temporal variation F I G U R E 1 Conceptual road map for study design, using a hypothetical example scenario. Bold numbering corresponds to sections within this review of one or a small number of niche-specific rates (e.g. legume trees, plus litter and soil) to be treated as equivalent to whole-ecosystem rates, although this can substantially underestimate total BNF (Tierney et al., 2019). For studies that aim to estimate total ecosystem BNF, we suggest creating a comprehensive list of potential BNF niches in the study system a priori and using pilot sampling to rule out any negligible BNF sources after considering both the mass-specific fixation rates and the abundance of each potential niche. Many other arguments can be made for where measurements should be focused, however, depending on the goal. At the scale of the continental United States, for example, Staccone et al. (2020) found that although some widespread woody fixing genera have been proportionally under-sampled (e.g. Prosopis and Cercocarpus), improving continent-wide BNF estimates would be best served by further increasing precision for the best-studied but highest-fixing genus (Alnus).
Consider also that a niche's BNF rate is not always indicative of its ecological significance, which is context and question dependent.
Ecosystems can differ in their reliance on new N inputs versus internally recycled N (Cleveland et al., 2013), which in turn can depend on contributions of individual niches to overall supply. For example, the dominant forest N supply often shifts from symbiotic BNF in early-successional forests to free-living in late-successional forests (Menge & Hedin, 2009;Sullivan et al., 2014;Taylor et al., 2019).
Likewise, individual niches may differ considerably in their importance to other ecosystem processes. Although most fixed N is probably ultimately cycled through multiple niches, symbiotic or endophytic BNF initially directly supports plant growth and ecosystem C gains, whereas free-living BNF in leaf litter and soil more directly supplies N for decomposition that regulates ecosystem C loss.
BNF is highly spatially and temporally heterogeneous, both within and across niches and ecosystems, as a function of the patchy distribution and diversity of N-fixing organisms (Cleveland et al., in review;Staccone et al., 2020;Taylor et al., 2020), environmental heterogeneity, and the variable controls on the process itself (Reed et al., 2011;Zheng et al., 2019). Identifying the major potential drivers of variability for each niche in question serves two purposes.
First, understanding the drivers of BNF variability is itself a primary scientific goal of many BNF studies. Second, understanding the heterogeneity in BNF within a study system is a prerequisite for deciding where and when to replicate sampling (Figure 1). Identifying major sources of heterogeneity through pilot studies or relevant literature can aid with selection of methods, which integrate variation across different scales.
Many of the primary controls on BNF have both spatial and temporal axes of variation, which may operate across small (e.g. microsite, hourly) to large (e.g. landscape, decadal) scales (Reed et al., 2011;Smercina et al., 2019a; Figure 2). At the molecular to organismal level, direct controls include conditions that affect nitrogenase enzyme activity, metabolic state, the absolute and relative availability (stoichiometry) of resources, or source-sink dynamics (e.g. growth, herbivory). Controls at this level can act independently (e.g. light) but interactions between controls (such as moisture and oxygen) are also very common. Complex, higherlevel drivers typically integrate changes in multiple direct controls; disturbance, for example, is likely to simultaneously impact the availability of light, C, nutrients, and moisture. At larger scales, ecological controls shape the competitive ability or relative abundance of N-fixing organisms, which in turn can interact with local controls to shape realized BNF rates. Such multi-scale interactive controls on BNF are an important yet complex research frontier. Deliberate a priori consideration of these scale(s) of interaction can focus study design to address the major drivers of BNF variation, which can include:
• Temperature (direct control on nitrogenase activity and microbial growth; indirect control on plant species composition and growth): diurnal and seasonal variation; spatial variation with microclimate and elevation (Caputa et al., 2013;Houlton et al., 2008).
• Light (for photoautotrophs): diurnal and seasonal variation; spatial variation with microclimate and successional stage (Myster, 2006;Taylor & Menge, 2018).
Temporal and spatial scales of variation of some potential controls on terrestrial BNF in free-living and symbiotic niches. Some controls operate broadly across space and time (e.g. moisture), while others operate at much smaller spatial scales (e.g. O 2 concentration) and relevant scales of variation differ between niche types. For example, BNF may be strongly regulated by micro-scale variation in nutrient availability for soil microbes, but at scales larger than the individual root system for N-fixing trees and shrubs
• Moisture (direct control on microbial growth and nutrient diffusion, and on physiological activity for most organisms): shortterm, diurnal, seasonal and inter-annual variation with rainfall, evaporative demand and inundation; spatial variation with microclimate and soil type (Caputa et al., 2013;Hicks et al., 2016).
• Oxygen concentration (direct control for free-living niches as it inhibits nitrogenase activity): temporal and spatial variation with soil/substrate conditions including respiration rate and moisture; spatial variation with soil structure/texture and microtopography (Smercina et al., 2019b).
• Carbon availability (for heterotrophs, direct control on energy available for BNF and growth/nutrient demand): seasonal variation with root exudation, litterfall and root mortality; temporal variation with rainfall and degree of decomposition of litter/wood (Reed et al., 2011).
• Nitrogen availability (increased N availability commonly suppresses both free-living and symbiotic BNF activity; Dynarski & Houlton, 2017;Vitousek et al., 2013): seasonal variation with soil temperature and moisture effects on microbial cycling and plant phenology; spatial variation with microsite, plant species distribution, successional stage and atmospheric deposition (Dynarski & Houlton, 2017;Zheng et al., 2019).
• Availability of other nutrients (particularly phosphorus, as a limitation on growth and metabolic activity) and nitrogenase enzyme co-factors (molybdenum, vanadium, iron): spatial variation with plant species distribution, soil type/parent material and external inputs of dust or salt spray (Vitousek et al., 2002;Wurzburger et al., 2012).
• pH (indirect control on microbial community composition; reduced soil BNF observed at low pH): spatial variation across the rhizosphere and with soil type (Smercina et al., 2019b).
• Plant phenological state (indirect control influencing source-sink dynamics and C supply for symbiotic BNF): seasonal variation, often with reduced activity during periods of extreme temperature or low moisture (Gei, 2014;Myrold et al., 1999;Pearson & Vitousek, 2001).
• Degree of plant host-symbiont regulation: on a spectrum from obligate to facultative (Menge et al., 2015).
• Host community composition (e.g. for endophytes): spatial variation in abundance (Reed et al., 2008).
• Disturbance (indirect control that influences multiple direct physical and biological controls such as light availability, temperature, moisture, C and N availability): temporal and spatial variation in fire-affected areas, tree fall gaps etc. (Tierney et al., 2019).
• Successional stage (indirect control on long-term changes in light and nutrient availability; Zackrisson et al., 2004): temporal and spatial variation in species composition (Batterman et al., 2013).
It appears from the BNF literature that temporal variation is probably under-sampled in almost all cases. As with many processes, there is likely bias against sampling during periods with anticipated low activity or degree of inconvenience, such as night time, during periods of desiccation, or during winter in highly seasonal ecosystems (Heath et al., 1988). We recommend testing these assumptions explicitly whenever possible (e.g. taking at least some measurements during these periods) if they are necessary for downstream scaling to annual rates.
As with many biogeochemical processes, some BNF drivers may be characterized by legacy effects where rates depend not only on conditions at the time of measurement, but also those preceding.
Availability of C and N in soil, for example, is highly dependent on previous wetting duration and intensity as well as current moisture content (Fierer & Schimel, 2002). Legacy effects may apply to ecosystems such as tropical dry forests or deserts characterized by pulsed resource availability-for example, where an initial flush of C and nutrient release from necromass at the onset of the wet season may diminish with subsequent rainfall events (Waring & Powers, 2016)-and this temporal heterogeneity should be factored into sampling design.
Hotspots are areas that have high activity relative to the surrounding matrix and are a common phenomena in biogeochemical nutrient cycling (Kuzyakov & Blagodatskaya, 2015). Similarly, hot moments are transient periods where increased rates of a process are often associated with abrupt changes in physical conditions (such as the onset of rainfall), increased resource availability, or crossing of thresholds (e.g. from an oxic to an anoxic state). Many studies have reported hotspots of BNF in both symbiotic and free-living niches (Reed et al., 2011;Winbourne et al., 2018;Wong et al., 2019) and response curves for various controls suggest that hot moments almost certainly also occur (Caputa et al., 2013;Roley et al., 2019). Of the two, hotspots are generally easier to identify because it is far more common to collect multiple samples in space than to repeatedly resample the same point over time. Many questions remain as to the scale of free-living BNF hotspots in natural systems (do they operate at the scale of centimetres? Micrometres?) because studies generally do not measure BNF at multiple scales simultaneously. However, it seems reasonable to hypothesize that they correlate with the spatial scale of resource availability, microsite conditions, N-fixer distribution, or their interactions.
Hotspots can make up a large proportion of BNF activity and capturing them is an essential consideration in experimental design. In tropical forests, Wong et al. (2019) found that 1%-6% of leaf litter measurements accounted for between 57% and 71% of total free-living BNF and Winbourne et al. (2018) found that 7% of nodule-containing soil cores accounted for ~50% of symbiotic BNF (Figure 3). Because most measurement methods incur lag times of days to months between sampling and access to results, hotspots/moments generally cannot be detected 'on the fly'. Pilot studies can be one way to estimate the distribution of rates a priori and thus assess sampling intensity necessary to capture these dynamics (Figure 1, see Section 3.4).
Choosing the 'best' estimation method depends on the niche and ecosystem, but also reflects logistical and financial resources for both field sampling and sample analysis. If applicable, this choice should explicitly consider the availability of method-specific data needed to upscale rates (Tables 1 and 2, Section 4.1). Different methodological approaches integrate varying spatial and temporal scales. Instantaneous incubations measure nitrogenase enzyme activity either directly (based on the rate of incorporation of 15 Nlabelled N 2 into tissue; Unkovich et al., 2008), or indirectly (by substituting acetylene-which is reduced to ethylene by nitrogenase-and determining the rate of ethylene accumulation, i.e. the acetylene reduction assay, ARA; Hardy et al., 1968), both over a period of minutes to hours. At the individual plant scale, the xylem sap ureide method (applicable only to certain Fabaceae species) measures the transport of assimilated fixed N to above-ground tissues (Unkovich et al., 2008). Also at the individual plant scale, the 15 N natural abundance technique uses the isotopic composition of plant tissues as an indirect measure of the proportion of N derived from BNF versus soil sources (Boddey et al., 2000). Usually applied at the pot or plot scale in plants, 15 N dilution measures the dilution of an applied isotopic label by fixed N over the course of months to years (Unkovich et al., 2008). Mass balance approaches calculate BNF as the difference between system N inputs, outputs, and accrual in a bounded system, which can vary from a pot to a watershed (Cleve et al., 1971;Soper & Sparks, 2016). In addition, radioisotope ( 13 N 2 ) incubations have also been used where the goal is to simply distinguish the presence/absence of BNF activity (Moyes et al., 2016). The potential for BNF activity can also be identified by molecular methods that detect the presence, quantify the abundance, or measure the expression of the nif genes encoding the nitrogenase enzyme and can also be used to parse the taxonomic identity of the N-fixing organism (Zehr et al., 2003).
Several previous studies provide detailed instructions for these methods and discussion of their technical pros and cons, including Danso et al. (1992), Myrold et al. (1999), Unkovich et al. (2008) and Chalk et al. (2015). Here, we highlight two specific methodological issues. First, the 15 N natural abundance method is often used to estimate the percent of N derived from the atmosphere (%N dfa ) due to its relative ease of application. Yet, rarely have all of the original assumptions of the method (Shearer & Kohl, 1986) been tested or met in natural systems, including sufficient isotopic separation, comparable root phenology/morphology and comparable isotopic composition of assimilated soil N between fixers and reference plants (Boddey et al., 2000;Soper et al., 2015). Given this potential for low accuracy, we caution against using this approach to derive quantitative rates except in very specific circumstances (e.g. controlled agricultural settings). Second, ARA is a contentious method due to an inconsistent relationship between acetylene reduction and N 2 fixation and other well-recognized artefacts (Giller, 1987). This has led to a push to use the 15 N 2 incorporation method in its place (e.g. Smercina et al., 2019a), although the latter also comes with drawbacks (Chalk et al., 2017). However, ARA remains a useful tool, and many of its issues can be addressed with use of appropriate controls and 15 N 2 calibrations (see Section 3.3). This less resource-intensive method remains a technique of choice (especially for free-living BNF) where larger numbers of replicate samples are necessary to capture spatial and temporal heterogeneity, to increase measurement precision, or where relative differences rather than absolute rates are of interest. Recent improvements in instrumentation also mean that some ARA applications can detect changes in BNF over time scales (seconds to minutes) not possible with other approaches (Bytnerowicz et al., 2019;Cassar et al., 2012).
Optimization of all sampling and measurement approaches should consider whether it addresses accuracy (how close the measure is to the true value) or precision (how well constrained values are). For example, increasing effort expended on counting nodules increases precision, as does increasing the number of replicate samples analysed for ARA because these measures are often highly variable. On the other hand, carefully calibrating ARA incubations and using appropriate statistical treatments increases accuracy.
F I G U R E 3 Distributions of biological nitrogen fixation rates in leaf litter (ng N g -1 hr -1 ) and legume root nodules in individual soil cores (ng N cm -2 hr -1 ) in tropical forests shown on a linear or log scale with arithmetic (dashed line) and geometric (solid line) means. Data from Osborne et al. (2020) and Winbourne et al. (2018)
In many instances where the goal is to estimate niche-specific BNF rates, a survey approach is recommended that applies instantaneous (ARA orfoot_0 N 2 incorporation) measures to replicate samples in space and/or time to capture representative variation in major environmental drivers such as temperature and moisture. The advantage of surveys is that they capture BNF under realistic combinations of environmental variables, and in response to conditions that are difficult to manipulate experimentally. Downsides include the requirement for temporal interpolation between sampling points, the fact that multiple drivers may covary in ways that are difficult to disentangle, and the potential to miss or overrepresent hotspots and moments that contribute disproportionately to overall rates.
An alternative approach is to experimentally establish relationships between BNF rates and environmental drivers, which helps to determine mechanisms and can aid with upscaling.
Manipulative experiments generating BNF response curves (typically to light, temperature, moisture and nutrients) are relatively common (Caputa et al., 2013;Hicks et al., 2016;Wolf et al., 2017).
Where it is assumed that these small number of factors capture the majority of the variation in BNF rates, response curves can be matched to recorded environmental variation to scale BNF using actual experienced conditions. At the ecosystem scale, this is
TA B L E 1 Common methods for measuring biological nitrogen fixation in terrestrial ecosystems Method Niche Temporal scale Spatial scale Cost per sample Acetylene reduction assay (ARA) Free-living, symbiotic Short (mins-hrs) Small (cm-m) $ ($$ for 15 N 2 calibration) Sap ureides Certain Fabaceae only Short (mins-hrs) Medium (individuals) $ 15 N dilution Symbiotic Medium (weeks-years) Medium (individuals) $$$ 15 N natural abundance Some symbiotic Medium (months-years) Medium (individuals) $$ Mass balance Integrated (free-living and symbiotic) Long (years-decades) Medium-Large (m-km) $-$$$ TA B L E 2 Minimum reportable variables and ancillary or co-variate data for scaling or assessing common controls for BNF studies Minimum data for all studies Ancillary data Symbiotic Free-living Controlling variables Scaling variables Controlling variables Scaling variables Latitude and longitude Biome type and vegetation description MAT, MAP (available via global databases) Land use (unmanaged vs. managed, description) Niche and habit (e.g. ground, epiphyte) Taxon (species or genera and family, if known) Sampling time (year, month, season) Number of biological and technical replicates (reported with error terms) Where applicable R ratio for converting ARA data Isotopic value assigned to fixed N Soil taxonomy Incubation conditions: temperature, moisture (ambient or manipulated), light intensity and duration, incubation length Soil characteristics (pH, total N and P, inorganic N and P at the time of BNF sampling) Light availability Tissue C:N:P Timing and extent of major disturbance (e.g. fire) Land use history, age of secondary forests Spatial Above-ground biomass Below-ground biomass (measured or estimated allometrically) NPP N fixer % cover, stem abundance or basal area Tissue N content and NPP (for %N dfa ) Nodule mass and density Temporal Seasonal rate variation, or temporal variation as a function of temperature/ moisture Nodule biomass and density Substrate C:N:P Soil characteristics (pH, total N, P and C, extractable inorganic N and P at a depth relevant to BNF sampling) Ambient substrate moisture Wood decomposition stage Spatial Bulk density and profile depth (soil) Rhizosphere volume (root-associated) Mass per area (leaf litter, wood, canopy leaves, bryophytes, lichen) Cover area (biocrusts, lichens, bryophytes) Abundance per unit mass (endophytes and epiphytes) Temporal Diurnal rate variation Seasonal rate variation Seasonal abundance usually achievable only for free-living niches (e.g. decaying wood or biological soil crusts; Caputa et al., 2013;Hicks et al., 2016) although this approach has been applied to ecosystems dominated by a single climate-sensitive fixing tree (Mitchell & Ruess, 2015).
Controlled measurements increase the ability to identify and capture hot moments or threshold responses (those these may rely on specific combinations of multiple variables) that might be missed with a survey, and can be used to interpolate BNF rates between field sampling events. In addition, Earth system and dynamic vegetation models often represent BNF as a function of continuous environmental variables (evapotranspiration, net primary productivity or temperature; Wieder, Cleveland, Lawrence, et al., 2015) so that defining empirical relationships between BNF and these environmental controls also supports modelling efforts.
Symbiotic BNF is commonly estimated by excavating and quantifying root nodule biomass, measuring instantaneous activity rates in a subsample of nodules, and extrapolating across space and time. Because this approach assumes a strong correlation between nodule biomass and BNF rate, this relationship should be explicitly quantified across a broad sample of nodules (Sullivan et al., 2014;Winbourne et al., 2018).
Differences in phylogeny, nodule age/morphology and fixation strategy, for example, may lead to nonlinear relationships. Spatial sampling designs vary in their assumptions about the heterogeneity of BNF across a landscape and require different scaling data, and can be grouped into two categories: area-based and plant-stratified sampling. Both approaches can be applied to scaling studies, but for comparative studies (e.g. different levels of fertilization) the choice will probably depend on the spatial extent of the treatment imposed.
Either approach can adopt an adaptive cluster sampling (ACS) scheme (which assumes that nodules are rare across the landscape and found in clustered populations) where additional soil cores are taken adjacent to ones where nodules are known to be present (e.g. Sullivan et al., 2014). Regardless of the approach, repeated measurements of nodule biomass are often necessary to account for possible temporal (e.g. seasonal) variation (Gei, 2014). In mesic systems nodules are typically found to a depth of 10 cm; however, in dry or seasonally dry systems, nodulation may occur too deep in the profile (down to 1 m or more) to be readily accessible (Johnson & Mayeux, 1990) and alternative methods must be considered.
Approaches that make no assumptions about the underlying distribution of nodules include simple random sampling (Pearson & Vitousek, 2001) and systematic grid- (Taylor et al., 2019) or transect-based surveys (Sullivan et al., 2014;Winbourne et al., 2018). These systematic approaches ensure broad spatial coverage and sample areas away from conspicuous N fixers that plant-stratified sampling omit. They may be a good initial choice when knowledge of BNF dynamics is limited, when individual N fixers are abundant, or where individual rooting systems cover a large spatial extent. Random sampling more strictly adheres to assumptions of many statistical tests, but small sample sizes are prone to omitting large areas of a plot. Where ACS is applied (assuming that underlying assumptions are true), this approach can reduce estimation error by increasing certainty in the magnitude of hotspots detected as a result of more intensive sampling. While plot-based random or fixed area sampling approaches rely on fewer assumptions in the scaling process than plantstratified approaches, the number of cores required to obtain robust areabased estimates is strongly dependent on the nodule detection rate. This can make planning sampling efforts challenging; for ACS the total number of cores is not known a priori (Winbourne et al., 2018) but pilot data on nodule distributions and variability in BNF rates can help inform sampling design. Example code for simulating these distributions based on an ACS sampling scheme and scaling to plot-scale BNF rates is described in the Supporting Information and publicly accessible (Taylor, 2021, Part 1).
Sampling schemes that stratify by the area around N-fixing plants assume that the majority of symbiotic BNF occurs in their immediate proximity, and may be a better choice when N fixers are relatively rare. This can be done at random, at defined locations (such as along cardinal directions from a stem), and can also apply ACS (although we know of no published examples). Stratified sampling provides relatively robust BNF rate estimates at the individual plant scale, which are scaled spatially by multiplying by fixer prevalence (Batterman et al., 2013), but require greater sampling intensity when fixer diversity is high. If nodules are clustered around plants, this approach may be preferable because it allows for greater representation per unit sampling effort, and can reduce estimation error if variability in rates of BNF per core is greater across stratum than within stratum.
We recommend verifying the spatial distribution of nodules around individual N fixers prior to adopting this approach to avoid omitting nodules that may be positioned relatively far beyond the canopy.
15 N dilution uses trace levels of enriched 15 N-labelled compounds to increase the isotopic ratio of the soil N pool. By comparing tissue 15 N of fixers to non-fixing reference plants, the fraction of N derived from fixation (%N dfa ) can be calculated and then scaled to BNF rates using biomass N accumulation. Collection of these scaling data must also be factored into sampling effort (see Section 3.4). This method integrates BNF rates over time and the spatial extent of 15 N application (typically plots of up to tens of square metres, over months to years) with sensitivity that is independent of BNF rate, making it the benchmark for estimating symbiotic rates (Unkovich et al., 2008).
However, the high cost of applying 15 N to plots large enough to encompass plant lateral rooting extent limits its practical use to herbs, shrubs or juvenile trees (Yelenik et al., 2013); application to mature trees is possible but has not yet occurred. Continued declines in the cost of 15 N-enriched materials and analysis may lead to increased use of this approach, particularly in ecosystems where estimating nodule abundance and seasonal activity is difficult.
Several considerations must be adopted into study designs applying 15 N dilution. Selecting application rates requires careful estimation of the quantity of enriched 15 N needed to reliably increase soil available 15 N values above natural background variability (and ideally uses multiple 15 N applications over time), while at the same time maintaining low enough overall additions to avoid suppressing BNF via a fertilization effect. This requires a priori estimation of soil and vegetation N pools. The availability of suitable nearby reference species of similar growth habit and rooting characteristics is a prerequisite for site selection. The method is best suited to settings where surface-applied 15 N can be available for root uptake, such as in coarse-textured soils or mesic-to-wet ecosystems (Busse, 2000;Yelenik et al., 2013). Example code for simulating values derived from 15 N 2 dilution experiments and converting to %N dfa values is provided in Taylor (2021, Part 2) and described in the Supporting Information.
Free-living BNF is typically measured via instantaneous approaches and subsequent scaling of rates is based on the mass and/or areal cover of the niche across the study area. If the goal is to estimate total free-living BNF at a landscape scale, or to assess how environmental variables affect free-living BNF broadly, care should be taken to ensure that all potentially important niches are sampled (see Section 2.1). Sampling efforts are generally divided between the BNF assay itself and estimating niche distribution for scaling.
Estimating the relative variation in each term is essential to plan sampling intensity (Figure 1). Estimating cover area/biomass/volume is relatively straightforward for niches such as leaf litter, bulk soil or biological soil crusts, but may require time-consuming data collection for less accessible or discrete niches such as canopy endophytes, epiphytes or rhizosphere soil.
Because free-living BNF is heterogeneous across small spatial scales, sampling should encompass the breadth of substrate conditions experienced by a given niche in a given environment. As with symbiotic niches, spatially random (e.g. using randomized GPS points) or systematic (grid or transect) sampling schemes can be used to minimize bias, but must be balanced with the potential to under-sample spatially rare niches and account for high variance due to hotspots.
Because free-living BNF is highly sensitive to abiotic conditions (Reed et al., 2011), seasonality should be explicitly incorporated into sampling designs and diurnal variation may also be important, especially for photoautotrophs (e.g. cyanobacteria). Note that accurate 15 N 2 incubation estimates also require assessment of background 15 N levels; the greater the background variation between individual samples, the higher the replication needed for precise quantification of BNF.
As an alternative to summing inputs from individual niches, N mass balance methods infer total BNF as the mass difference between N inputs (wet and dry deposition), N outputs (gaseous emissions, leaching and hydrologic loss of inorganic and particulate N, biomass removal, fire) and N accrual (in soil and biomass pools; Soper & Sparks, 2016;Soper et al., 2017;Tierney et al., 2019). Calculating N balance is relatively straightforward in pot or microcosm studies where inputs can be tightly controlled. For ecosystems, however, establishing system boundaries where data are available or sampling is feasible (such as a watershed or defined forest area) is essential. This method is best suited where BNF rates are likely to be high, particularly relative to measurement errors in other components. If an ecosystem is assumed to be at steady state (e.g. a mature forest; Cleveland et al., 2010) only inputs and outputs are measured, while accrual must also be accounted for during primary or secondary succession. For the latter, at least two sampling sites of different known age or time since disturbance are required (typically separated by years to decades); additional chronosequence sites add the ability to estimate changes in BNF inputs over time, as this function is not necessarily linear (Soper & Sparks, 2016).
Certain ecosystem N fluxes are known to be highly episodic (especially gaseous and hydrologic losses; Taylor et al., 2015) and obtaining accurate estimates can require frequent or continuous sampling (Barton et al., 2015). It can also be difficult to measure small changes in large N pools, such as bulk soil, even with large sample sizes (Chalk, 2020;Turner et al., 2019). The impact of unmeasured fluxes (such as N 2 loss) on N balance can be explored using sensitivity analyses (e.g. Soper et al., 2015Soper et al., , 2017)). Redistribution of N from unmeasured pools (such as deep soil or rock weathering) into measured soil and biomass is often unaccounted for (Turner et al., 2019) but could be similarly constrained. Appropriate error propagation is important for mass balance, where uncertainty can be high for multiple parameters.
Incubation-derived BNF rates are highly sensitive to measurement conditions such as moisture, temperature, C availability, disturbance (e.g. detaching nodules from roots), incubation length, and headspace gas concentrations (e.g. of O 2 ; Vessey, 1994;Chalk et al., 2017;Smercina et al., 2019a). These conditions can be manipulated by researchers (e.g. via additions of water or glucose) to measure 'potential' BNF rates. Because these do not accurately represent expected field BNF rates, potential rates have limited utility for scaling.
However, manipulations can deliver insight into the response of BNF to conditions like rain events that might not be captured in field sampling (and 'wet' rates are sometimes scaled assumpting the fraction of time samples would be wet in situ), to standardize conditions (to isolate the response of BNF to a specific field or laboratory treatment), to bound the upper limits of BNF potential, or to increase BNF rates above instrument detection limits.
Several common pitfalls of ARA or 15 N 2 incorporation can be avoided with careful methodological design. Typically, free-living BNF rates per unit mass are much lower than symbiotic and require longer incubation times for detection (typically 8-24 hr for soil or leaf litter).
Very short incubations may preclude measurable ethylene production (for ARA using gas chromatography) or sufficient isotopic enrichment (using 15 N). By contrast, long ARA incubations can impact rates and reduce accuracy by causing microbial community turnover (Fulweiler et al., 2015), altering microbial metabolism by blocking processes such as methanogenesis and nitrification (Hynes & Knowles, 1982), and derepressing nitrogenase synthesis (Silvester et al., 2011). For both ARA and 15 N 2 incorporation, overly long incubations can also result in shifts in the CO 2 /O 2 partial pressure of the headspace, although steady concentrations can be maintained in recirculating setups (Chalk et al., 2017).
Pilot studies can help determine the shortest incubation time that allows for reliable detection. For all niches (especially those that may have high rates of background ethylene production, such as leaf litter, fresh leaves or forest soils; Nohrstedt & Muller, 1983) blanks (sample-only or acetylene-only) and paired blanks (pre-sampling ethylene immediately prior to ARA) can account for non-BNF-derived ethylene and thus distinguish between true high rates and artefacts (Reed et al., 2008).
Although not yet widely adopted, cavity ring-down spectroscopy-based ethylene analysers have significantly lowered the ethylene detection limit for ARA (to ~0.2 ppb) and are an option for niches with very low rates, precluding the need for manipulations such as wetting (Cassar et al., 2012). This technique can reduce acetylene concentrations required from ~10 to ~2% and thus reduce artefacts arising from physiological disruption by acetylene (Bytnerowicz et al., 2019).
For ARA incubations where data will be converted from measured units (acetylene converted to ethylene) to units of N 2 fixed, a conversion factor ('R ratio') must be applied. Although many studies rely on the theoretical R ratio of 3:1 or 4:1 (Hardy et al., 1968), this ratio varies widely across niches and samples (Soper et al., 2021) and may be a significant source of inaccuracy. R ratios should thus be calculated directly using 15 N 2 calibrations performed at the same time and under comparable conditions as ARA. Studies concerned with the response of BNF to a treatment often bypass calibration and compare units of ethylene production directly. However, 15 N 2 calibration is still necessary for comparative studies when samples are measured under different physical conditions (such as degree of water saturation) that could impact the ratio. Calibration should also be used in systems where alternative nitrogenase (V-or Fe-only based isoenzyme) activity may be high (e.g. Mo-poor soils), because isoenzymes vary in their relative BNF efficiency (Bellenger et al., 2020). Appropriate controls are necessary at each step because some commercial 15 N 2 gas stocks have been shown to be contaminated with reactive N compounds, which can lead to overestimated BNF rates (Dabundo et al., 2014).
Finally, commercial (even 'high purity') acetylene often has high background ethylene concentrations; this contamination is substantially lower when fresh acetylene is produced directly using calcium carbide (Bytnerowicz et al., 2019).
Measurements of BNF frequently yield both undetectable and episodic high values (e.g. Roley et al., 2019), leading to an often non-normal, right-skewed distribution (Figure 3) and necessitating adequate replication especially for scaled-up rates. As such, it is highly beneficial to determine appropriate sample sizes a priori.
Power analyses can address two issues: the sampling effort needed to detect differences between experimental treatments or across environmental gradients, and to gain robust BNF rate estimates. For symbiotic studies, Winbourne et al. (2018) showed that the parameters that most strongly determine the sampling effort (e.g. number of soil cores) needed to gain accurate BNF estimates are 1) the fraction of sampled cores that contain nodules and 2) the variation in BNF rate per core. In theory, these parameters could be applied to studying other BNF niches by defining them as the proportion of samples where BNF is >0 and the coefficient of variation in BNF among non-zero samples. Both Winbourne et al. (2018) and Barton et al. (2015) offer good demonstrations of the effect of sampling frequency on capturing hotspots/hot moments.
A complementary approach is to conduct a sensitivity analysis of the various levels of data that are incorporated into the final rate calculation. For example, implicit in estimates of annual free-living BNF in leaf litter measured via ARA are variation and error in BNF per unit of litter (which incorporates error in background ethylene production and C 2 H 4 :N 2 conversion factor), variation and error in mass of litter per unit ground area, and temporal heterogeneity in these measurements throughout the year. Using sensitivity analyses to understand which sampling component has the largest influence on the overall BNF estimate can help determine priorities for sampling effort (to increase precision) among these components (Figure 1). As a first pass, assessing the ratio of variance to the mean for each data component will indicate which variables contribute disproportionately to the propagated error. Both power and sensitivity analyses rely on reasonable a priori estimates of means and variation for each data source. These estimates may be available in the literature, but the substantial regional variation often necessitates site-specific pilot data in order to run informative analyses.
Many empirical measurement methods for BNF produce rate data that require conversion, for example from 'mass of N fixed per unit sample mass/area per incubation time' to commonly reported units of 'mass of N fixed per unit ground area per unit time'. This presents temporal and spatial scaling challenges. How to estimate the biomass and/or spatial coverage of each niche? Does it vary temporally? Should linear interpolation be used to infill between measurement points or should seasonal averages be applied? Are necessary data available at a scale compatible with BNF measurements? The categories of data generally required to scale BNF from different niches are listed in Table 2, but the answers to these questions are casespecific. Although not available for every site, data for vascular plant species abundance and cover for certain niches (e.g. lichens, mosses and decomposing wood) are increasingly accessible via public repositories and national forest inventories.
Although decisions and assumptions made during scaling of point data strongly influence final extrapolated rates, surprisingly, descriptions of scaling calculations are often brief. Where data are scaled, we advocate explicitly identifying:
1. Source of scaling data and its sampling frequency, temporal and spatial extent, and appropriate descriptive statistics (mean, sample size and variation or distribution).
2. Spatial assumptions: relevant system boundaries, abundance and extrapolations (e.g. extrapolation of rates measured at a given soil depth to another depth).
3. Temporal assumptions and extrapolations (e.g. presence/absence of BNF at night, extrapolation of summer-measured rates to the entire growing season).
4. Values assigned to R ratios for ARA incubations and isotopic end members, their sources and associated error, and parameters for any response functions (e.g. between BNF and temperature/ moisture). 5. Methods used for propagating error from the various data sources.
The inherent methodological issues discussed above (capturing heterogeneity, sampling layout) make the appropriate statistical analysis critical for obtaining accurate, precise estimates of ecosystem N inputs or robust BNF rate comparisons.
BNF rates often follow a lognormal distribution (Figure 3) and data are often further complicated by being strongly zero-inflated, with several important implications. First, because hotspots and moments of activity make up the right-hand tail of the lognormal distribution (Kuzyakov & Blagodatskaya, 2015), BNF datapoints that look like statistical outliers are often not only real, but the most important points in the dataset. Anomalously high points are critical to include in analysis, although there remains uncertainty in how to weight such data when computing total BNF. Second, BNF data distributions often violate the assumptions of normality required for many traditional statistical tests. The issue of zero-inflated data can be addressed by aggregating samples (i.e. even if most individual samples exhibit no detectable BNF rates, there will be many fewer zeros when aggregating all samples within a subplot or plot), but this may reduce the effective sample size and statistical power such that it is only viable for large datasets. While some zeros are true (i.e. absence of nodules), for free-living BNF it is often not possible to distinguish between a true zero (no activity) and activity that falls below the analytical detection limit. One option is to assign these ambiguous zeros a value of 1/2 or 1/3 of the detection limit (Barron et al., 2008) or to use maximum likelihood to estimate the distribution of values below detection (Menge & Hedin, 2009). The lognormal nature can be addressed either by transforming data prior to analyses (typically a log transformation) or by analysing using model structures that do not assume normally-distributed data (e.g. Taylor et al., 2019). The geometric mean (mean of the lognormal distribution) should probably also be used in place of arithmetic means to reduce the influence of very high values (Taylor et al., 2019). Other statistical methods appropriate for zero-inflated continuous data, such as bootstrapping (Paneru et al., 2018), may provide better estimates of distributions but have not yet been widely adopted by the biogeochemical community.
Propagating error or uncertainty from multiple sources is critical for reporting BNF estimates transparently. The most straightforward approach is to propagate the proportional variance of each term (example code in Taylor, 2021, Part 1). A second is to apply parametric or nonparametric bootstrapping, the latter of which has the advantage of not assuming any particular distribution. A good example of error propagation for isotopic %N dfa calculations can be found in Menge et al. (2015). Where applicable, we recommend moving away from applying discrete values for certain scaling parameters (e.g. the isotopic value assigned to fixed N, isotope end members and R ratios) and instead applying a distribution that reflects the uncertainty in their measurement and generates a probability distribution of BNF rates (Taylor & Menge, 2018, example code in Taylor, 2021, Part 2).
Comparison or meta-analysis of broad rates and controls on BNF is often hampered by a paucity of accompanying information on site and sample characteristics, collection methods, and analysis conditions. At a minimum, we recommend a set of basic descriptive variables that should accompany any BNF study (Table 2). We also suggest reporting where feasible additional site and sample characteristics that can greatly increase the utility of a dataset beyond the immediate study to identify BNF patterns and controls, aid scaling for regional and global rate estimates, and improve, parametrize, and test model representations. For example, despite ample evidence from individual studies that the availability and stoichiometry of C and nutrients (N, P, Mo) control BNF rates (Dynarski & Houlton, 2017;Reed et al., 2011), identifying whether these relationships hold at larger scales has been hampered by a lack of data.
As recently highlighted by Davies-Barnard and Friedlingstein (2020), we strongly advocate that publications include null results when no BNF is detected because the absence of this information biases extrapolations at larger scales. In addition, publication of either raw BNF rate data (best case) or means, error terms and sample numbers as supplementary tables (in addition to figures) makes secondary data usage significantly more streamlined and accurate. A standardized template for reporting BNF and associated metadata is provided in Table S1.
The final step is to consider data in light of the original goal. Did high N losses from a system suggest there is (or was) high BNF to find? If the goal was to explicitly evaluate BNF rates in the context of other known N fluxes of the system, caution should be applied in comparing measurements reflecting different temporal or spatial scales. In almost all cases, it is probably preferable to treat any estimate of BNF as a range, and to consider sensitivity analyses when trying to balance a N budget. The numerous sources of uncertainty involved in estimating BNF mean that plot-scale rates can have percent uncertainties of ±100% or more, and logistical constraints on sampling effort often prohibit sampling exhaustively to increase precision (e.g. when very large sample sizes are required; Winbourne et al., 2018).
As for using traditional statistical tests, the inherent variability in BNF rates often masks strong statistical significance even when BNF differs substantially across a gradient or between treatments.
Relating BNF rates to other measured ecosystem properties and processes using techniques such as model selection with the Akaike information criterion (Akaike, 1974), the Bayesian information criterion (Neath & Cavanaugh, 2011) or a classification and regression tree approach (De'ath & Fabricius, 2000) can provide useful information about ecological significance but are not yet widely applied.
Given that different methods rely on inherently different principles to arrive at BNF rates-many associated with large potential sources of measurement error-cross-calibration between methods to determine systematic biases would seem to be a priority. In fact, this has rarely been performed in natural ecosystems (Boddey et al., 2000), likely because often only one or two methods are suitable for a given site, the resource cost of replicating methods is high, and the scaling uncertainty associated with varying temporal and spatial extents is considerable (Winbourne et al., 2018). Despite this, a few conclusions stand out and can be supplemented by findings from managed ecosystems. First, agreement between methods tends to be greater at higher BNF rates and/or dependence on BNF as an N source (Binkley et al., 1992(Binkley et al., , 1994;;Busse et al., 2007;Chalk et al., 2015).
Second, natural abundance isotopic comparisons vary greatly with local site characteristics (availability of reference species, degree of isotopic separation and physiological diversity) and should be applied with great caution, including explicit tests of relevant assumptions (Soper et al., 2015). Finally, woody species are a particularly important avenue for further methodological development because most comparisons focus on juvenile seedlings which differ in their physiology and source/sink dynamics compared with mature trees (Thomas & Winner, 2002).
In woody legumes, the 15 N dilution and natural abundance methods vary in their agreement and have been shown to vary by as little as 10% in Robinia (Marron et al., 2018) or up to a factor of five in Acacia (Bouillet et al., 2008;Hamilton et al., 1993). Seasonal or interannual variability (for both approaches) and lack of adequate isotopic separation between fixing and reference plants (natural abundance) are identified as a likely source of error or disagreement (Hamilton et al., 1993;Soper et al., 2015). Comparisons are more common for herbaceous and agricultural legumes, where greater spatial homogeneity may reduce measurement error. A meta-analysis by Chalk et al. (2015) regressed the estimates derived from 15 N natural abundance and dilution for predominantly agricultural herbaceous legumes and some woody seedlings. Despite a positive linear relationship, only around half of values fell within 10% of each other and the degree of agreement was related to the proportion of N derived from symbionts versus soil. For context, however, a 10% variation may be considered as relatively precise in a natural ecosystem setting.
Mass balance and ARA/ 15 N 2 incubations arguably offer the most independent cross-comparisons between methods because they do not rely on overlapping measurements and scaling data, as comparisons of isotopic approaches usually do. Several studies contrast these approaches for chronosequences and regenerating secondary forests. When scaled ARA rates do not match measured N accumulation, the difference has generally been attributed to poorly constrained fluxes (such as soil organic matter turnover or gaseous/hydrologic losses), rather than a failure of ARA to accurately capture rates (Anderson et al., 2004;Pearson & Vitousek, 2001;Uliassi & Ruess, 2002). In managed systems, in particular North American alder Alnus rubra plantations, relatively high BNF rates (75-130 kg N ha -1 year -1 ) and greater species/substrate homogeneity of single-aged stands probably increases the accuracy of both methods. In three cases where these techniques were performed side-by-side, alder BNF rates fell within error ranges of each other for a maximum difference of ~50% (Binkley et al., 1992(Binkley et al., , 1994)).
Finally, some studies have contrasted mass balance approaches with other techniques, with variable agreement. For the United States, Staccone et al. (2020) compared estimates of woody symbiotic BNF obtained via N accretion (inferred from plot-level biomass increases) and N demand (scaled from growth rates, stoichiometry and %N dfa ) and found that estimates were twice as high for the accretion method (3.4 vs. 1.4 kg N ha -1 year -1 ). The authors argue that accretion was likely to overestimate BNF via bias towards studying younger stands where accounting is easier but BNF rates are higher, despite the potential for unaccounted N losses to underestimate BNF. Conversely, N demand likely underestimated rates due to probable overestimation of the wood C:N ratio, as well as lack of accounting for N sinks such as bark, fruit and herbivory. Using a similar approach combined with N accretion estimates, data from a well-characterized sub-tropical Prosopis woodland suggested that BNF meets around 40%-45% of total plant N demand (Soper, unpublished data;Soper & Sparks, 2016). At the same site, however, natural abundance %N dfa generated significantly greater estimates (65%-80%) but was highly sensitive to inter-annual variability associated with climate (Soper et al., 2015).
Measurement and scaling remain among the most significant challenges to understanding and contextualizing BNF in terrestrial ecosystems, but careful study design and execution can alleviate many of the associated concerns and ensure that we continue to advance a robust, quantitative and predictive understanding of BNF. Clear opportunities for improvement include:
• More holistic sampling of both well-recognized and comparatively cryptic niches, including simultaneous consideration of both symbiotic and free-living sources (e.g. Tierney et al., 2019).
• Sampling/scaling that explicitly explores sources of heterogeneity and mechanisms of hotspots/moments of BNF, and improved statistical techniques for dealing with the zero-inflated, non-normal data typical of BNF.
• More widespread adoption of tools such as cavity ring-down spectroscopy that reduce artefacts for ARA (Cassar et al., 2012), allow for discrimination between low and zero BNF rates, and facilitate easier measurement of continuous BNF response functions for abiotic variables.
• More standardized data reporting to facilitate interpretation of patterns and controls at a broader scale.
• Additional explicit inter-comparison among methods, to both constrain rates and improve decision-making for future methods selection. For systems that are adequately bounded, basic biogeochemical modelling approaches based on N demand, stoichiometry, isotopic mass balance and fixer abundances could be used to 'reality-check' estimates derived from other methods (Cleveland et al., 2010).
Together, these advances can be used to address some of the significant remaining unknowns in ecological BNF research, such as the drivers of interspecific variation in BNF potential, the persistence of ecosystem N limitation, 'missing' ecosystem N inputs (Turner et al., 2019), the distribution of BNF strategies (Menge et al., 2015), and the importance of N-fixer diversity for BNF rates (Taylor et al., 2020).
This paper evolved from the 'Novel multi-scale synthesis of nitrogen fixation rates and drivers across the terrestrial biosphere' working group supported by the USGS John Wesley Powell Center for Analysis and Synthesis. FMS, CCC and SR acknowledge support from a National Science Foundation RCN grant (DEB-1754126). We thank Sian Kou-Giesbrecht, Caroline Dallstream and Mia Marcellus for helpful comments. The authors report no conflict of interest. Any use of trade, firm or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
N 2 incorporation Free-living, symbiotic Short (mins-hrs) Small (cm-m) $$$