Global oceanic diazotroph database version 2 and elevated estimate of global oceanic N₂ fixation

Marine diazotrophs convert dinitrogen (N₂) gas into bioavailable nitrogen (N), supporting life in the global ocean. In 2012, the first version of the global oceanic diazotroph database (version 1) was published. Here, we present an updated version of the database (version 2), significantly increasing the number of in situ diazotrophic measurements from 13 565 to 55 286. Data points for N₂ fixation rates, diazotrophic cell abundance, and nifH gene copy abundance have increased by 184 %, 86 %, and 809 %, respectively. Version 2 includes two new data sheets for the nifH gene copy abundance of non-cyanobacterial diazotrophs and cell-specific N₂ fixation rates. The measurements of N₂ fixation rates approximately follow a log-normal distribution in both version 1 and version 2. However, version 2 considerably extends both the left and right tails of the distribution. Consequently, when estimating global oceanic N₂ fixation rates using the geometric means of different ocean basins, version 1 and version 2 yield similar rates (43–57 versus 45–63 Tg N yr⁻¹; ranges based on one geometric standard error). In contrast, when using arithmetic means, version 2 suggests a significantly higher rate of 223±30 Tg N yr⁻¹ (mean ± standard error; same hereafter) compared to version 1 (74±7 Tg N yr⁻¹). Specifically, substantial rate increases are estimated for the South Pacific Ocean (88±23 versus 20±2 Tg N yr⁻¹), primarily driven by measurements in the southwestern subtropics, and for the North Atlantic Ocean (40±9 versus 10±2 Tg N yr⁻¹). Moreover, version 2 estimates the N₂ fixation rate in the Indian Ocean to be 35±14 Tg N yr⁻¹, which could not be estimated using version 1 due to limited data availability. Furthermore, a comparison of N₂ fixation rates obtained through different measurement methods at the same months, locations, and depths reveals that the conventional ¹⁵N₂ bubble method yields lower rates in 69 % cases compared to the new ¹⁵N₂ dissolution method. This updated version of the database can facilitate future studies in marine ecology and biogeochemistry. The database is stored at the Figshare repository (https://doi.org/10.6084/m9.figshare.21677687; Shao et al., 2022).

1 Introduction

Dinitrogen (N₂) fixation is a process carried out by select prokaryotes (diazotrophs) capable of converting N₂ gas, which is not usable by most organisms, into bioavailable nitrogen (N). In the sunlit surface ocean, where dissolved inorganic forms of N such as nitrate (NO${}_{\mathrm{3}}^{-}$) and ammonium (NH${}_{\mathrm{4}}^{+}$) are scarce, N₂ fixation plays an important role in providing N that can contribute to primary production, particularly in oligotrophic regions (Wang et al., 2019; Gruber, 2008). Globally, N₂ fixation serves to compensate, at least partially, for fixed N removed via denitrification and anammox (Deutsch et al., 2007; Gruber, 2019).

Marine diazotrophs include three main types of cyanobacteria (Zehr, 2011): (1) nonheterocystous filamentous cyanobacteria (e.g., Trichodesmium); (2) heterocystous cyanobacteria like Richelia, which may form diatom–diazotroph associations (DDAs); and (3) unicellular cyanobacteria (UCYNs). Non-cyanobacterial diazotrophs (NCDs) have also been widely detected in the ocean (Bombar et al., 2016; Delmont et al., 2021; Moisander et al., 2017). However, the contribution of NCDs to marine N₂ fixation has not been directly quantified, despite a few studies that have reported N₂ fixation by putative NCDs at the cellular level (Harding et al., 2022; Bentzon-Tilia et al., 2015a).

Diazotroph abundance has been estimated from nifH gene copies using qPCR assays (Church et al., 2005b) or droplet digital PCR (ddPCR; Gradoville et al., 2017). The abundance of some cyanobacterial diazotrophs can also be obtained by counting them directly using microscopy-based techniques and in some cases flow cytometry. A recent work combined an image recognition pipeline with molecular mapping of the nifH gene to quantify diazotrophs in the Tara Oceans dataset (Karlusich et al., 2021). Gene copies of nifH have been more frequently measured than microscopy-based cell counts and can be more useful when evaluating the abundance of different diazotrophic groups. Caution must be taken because there can be discrepancies between cell-count-based and nifH-based diazotrophic abundances (Luo et al., 2012), a finding largely attributed to large variations in the number of nifH copies per diazotroph cell, thus far observed particularly in Trichodesmium and heterocystous cyanobacteria (Sargent et al., 2016; White et al., 2018; Karlusich et al., 2021). However, a recent regional study spanning over 200 km of the North Pacific Subtropical Gyre has found a statistically significant linear correlation between the abundances of the nifH gene and cell counts in UCYN-B (i.e., Crocosphaera; linear slope = 1.82) and heterocystous cyanobacteria (Richelia and Calothrix; linear slope from 1.51–2.58) but not in Trichodesmium (Gradoville et al., 2022). A recent discussion highlighted the influence of the uncertainty in gene copy conversion to biomass and the need for further investigation of how to best take advantage of gene copy data for global diazotroph biogeography modeling purposes (Meiler et al., 2022; Zehr and Riemann, 2023); however, there is agreement that quantifying gene counts is a powerful tool for studying marine diazotroph distributions (Meiler et al., 2023; Zehr and Riemann, 2023). Meiler et al. (2023) proposed a number of topics of study for this field moving forward; Gradoville et al. (2022) concluded that “we hope that future studies report nifH : cell and explore the mechanisms controlling this ratio.” Both gene-based and microscopy cell counts have innate biases, which should be elucidated in future studies.

Given the importance of N₂ fixation to ocean ecology and biogeochemistry, it is imperative that a database of up-to-date N₂ fixation and diazotrophic abundance measurements be maintained. Currently, global estimates of marine fixed N inputs calculated via the N₂ fixation rate mostly range from 100 to 170 Tg N yr⁻¹ (see summary in Zhang et al., 2020). This value, together with other bioavailable N sources to the ocean including riverine input and atmospheric deposition, is considerably lower than estimates of N losses from the ocean such as denitrification, anammox, and sediment burial (Zhang et al., 2020; Gruber, 2008; Zehr and Capone, 2021). While the overestimation of the N losses cannot be ruled out, one of possible reasons for this imbalance is the inaccurate estimation of global marine N₂ fixation due to limited spatiotemporal coverage of rate measurements and the different methods employed in N₂ fixation assays (White et al., 2020). Another possible reason is the limited knowledge of ecological niches of N₂-fixing organisms. Over the last decade, the realm of marine N₂ fixation has been expanded to include numerous non-paradigmatic habitats. Coastal (Mulholland et al., 2012; Bentzon-Tilia et al., 2015b; Mulholland et al., 2019; Tang et al., 2020; Turk-Kubo et al., 2021), subpolar (Sato et al., 2021; Shiozaki et al., 2018a), and even polar ocean regions (Blais et al., 2012; Sipler et al., 2017; Harding et al., 2018; Shiozaki et al., 2020) have demonstrated N₂ fixation. Notably, N₂ fixation in aphotic waters remains debated (Bonnet et al., 2013; Farnelid et al., 2013; Selden et al., 2021b; Rahav et al., 2013a; Hamersley et al., 2011; Benavides et al., 2018a; Moisander et al., 2017). Other studies have also suggested that NCDs may be significant contributors to marine N₂ fixation (Shiozaki et al., 2014b; Turk-Kubo et al., 2022; Geisler et al., 2020; Delmont et al., 2021; Karlusich et al., 2021; Bombar et al., 2016; Moisander et al., 2017) and may occupy different niches than cyanobacterial diazotrophs (Shao and Luo, 2022).

Luo et al. (2012) compiled the first global oceanic diazotrophic database including in situ measurements of N₂ fixation rates and cell-count-based and nifH-based diazotrophic abundance. Several years later, two studies supplemented the database with a collection of some newly reported diazotrophic data (Tang and Cassar, 2019; Tang et al., 2019), although a substantial amount of additional data remained to be included. Here, we present an updated version of the global oceanic diazotrophic database with data not yet compiled. We describe the database information, a summary of the data updates, measurement methods, and data distribution. Furthermore, we conduct a first-order estimation of the global oceanic N₂ fixation rate using the updated version of the database. In light of the aforementioned concerns of nifH : cell and various N₂ fixation methods (see Sect. 2.3), we also discuss the significance of employing different methodological approaches to estimate N₂ fixation rates and abundance metrics. We use the data available in the database to analyze the discrepancies between N₂ fixation rates using ¹⁵N₂ bubble and dissolution methods, and compare the observed ranges of nifH gene copies and diazotrophic cell abundance.

2 Data and methods

2.1 Database summary

This study updated the original global oceanic diazotrophic database of Luo et al. (2012; version 1 hereafter) with new in situ measurements of N₂ fixation rates and abundances of diazotrophic cells and nifH gene copies. Together there were 55 286 diazotrophic data points in the updated database (version 2 hereafter; Tables 1–3), including 13 565 data points from version 1 (Luo et al., 2012), 6736 measured in 2012–2018 and compiled by two previous studies (Tang et al., 2019; Tang and Cassar, 2019), 26 597 data points measured in 1979–2023 and compiled by this study, and 8388 NCD data mostly from Turk-Kubo et al. (2022; see below). In version 2, some errors in the datasets of Tang et al. (2019), mostly caused by unit conversions, were also corrected.

Version 2 was composed of six main sub-databases: (1) 9231 volumetric N₂ fixation rates (5853 new data points; Tables 1 and 4); (2) 2590 depth-integrated N₂ fixation rates (1805 new data points; Tables 1 and 4); (3) 9040 volumetric cell abundances (4154 new data points; Tables 2 and 5); (4) 1784 depth-integrated cell abundances (859 new data points; Tables 2 and 5); (5) 29 655 volumetric nifH gene copy abundances (26 506 new data points; Tables 3 and 6); and (6) 2986 depth-integrated nifH gene copy abundances (2544 new data points; Tables 3 and 6). Please be aware that 2416 N₂ fixation rates were measured with incubation periods less than 24 h; they were listed in separate spreadsheets in the database for reasons discussed in Sect. 2.3. Additionally, we included a compiled NCD dataset (Turk-Kubo et al., 2022) in the database, which contained 7919 nifH gene copy abundances of primarily the most studied phylotype NCD Gamma A (Shao and Luo, 2022; Langlois et al., 2015), also referred to as 24774A11 (Moisander et al., 2012) and UMB (Bird et al., 2005), as well as other phylotypes, and updated the compilation with 469 additional nifH gene copy abundances of NCDs published more recently (Turk-Kubo et al., 2021; Sato et al., 2022; Moore et al., 2018; Reeder et al., 2022; Wen et al., 2022; Bonnet et al., 2023). We also collected 468 cell-specific in situ N₂ fixation rates and added them to version 2 (Table 7).

Depth-integrated data were either provided directly in published papers or calculated as part of this study for those vertical profiles with at least three volumetric data points in each profile. The measurements within a profile were first interpolated linearly with depth, with the shallowest datum representing the level between the sea surface and the depth of that datum. The profile was then integrated from the sea surface to the deepest recorded measurement. Most vertical profiles of N₂ fixation rates were measured within the euphotic zone, with a few studies extending measurements to several hundred meters or deeper. In these cases, we only integrated to the deepest data point above 200 m, taking into account the scarcity of aphotic N₂ fixation measurements in the global ocean and their controversial contribution to the global budget (Benavides et al., 2018a). As a result, it was possible that certain measurements below the euphotic zone but above 200 m were included in the integration. However, these measurements would typically have minimal impact on the depth-integrated N₂ fixation rates due to their low rates and limited vertical extent in this range.

N₂ fixation rates were measured for whole seawater samples, for different size fractions (>10 µm and <10 µm), or specifically for Trichodesmium and heterocystous cyanobacteria. When whole-water N₂ fixation rates were not reported, total N₂ fixation rates were calculated as the sum of the N₂ fixation rates of available groups.

The cyanobacterial diazotrophic abundance data in version 2 were grouped into three taxonomic categories: Trichodesmium, UCYN, and heterocystous cyanobacteria. The UCYN abundance data were further grouped into UCYN-A, UCYN-B, and UCYN-C. Four sublineages of UCYN-A, including UCYN-A1, UCYN-A2, UCYN-A3, and UCYN-A4, have been identified (Thompson et al., 2014; Farnelid et al., 2016). UCYN-A1 and UCYN-A2 have significant distinctions in the sizes and species of their symbiotic hosts, with the former living in relatively smaller hosts (Thompson et al., 2014; Martínez-Pérez et al., 2016; Cornejo-Castillo et al., 2016). Hence, in addition to recording the total nifH gene copy abundance of UCYN-A in our database, the nifH gene copy abundances of its sublineages were also included if reported. Heterocystous cyanobacterial abundance was grouped into Richelia intracellularis (het-1 and het-2, associated with Hemiaulus and Rhizosolenia, respectively) and Richelia rhizosolenia (het-3, previously named Calothrix sp., associated with Chaetoceros; Foster et al., 2022b).

Sampling information (latitude, longitude, depth, and time) was provided for each data point. Physical, chemical, and biological parameters, including temperature, salinity, and concentrations of nitrate, phosphate, iron, and chlorophyll a, were also included when available.

Table 1 Summary of the number of data points for N₂ fixation rates by category. Measurements with incubation periods of 24 h or less are summarized separately.

Table 2 Summary of the number of data points for diazotrophic cell abundances. UCYNs include UCYN-A, UCYN-B, and unclassified UCYNs. Heterocystous cyanobacteria include Het-1, Het-2, and Het-3.

Table 3 Summary of the number of data points for nifH gene copy abundances. UCYNs include UCYN-A1, UCYN-A2, UCYN-B, and UCYN-C. Heterocystous cyanobacteria include Het-1, Het-2, and Het-3.

Table 4 Summary of new data points of N₂ fixation rates added to version 2 of the database.

^a Data are reported by data providers as depth-integrated N₂ fixation rates (unlabeled data computed by integrating profiles of volumetric N₂ fixation rate data). ^b N₂ fixation rate incubation time for 24–30 h.

Table 5 Summary of new data points of cell-count-based abundances added to version 2 of the database. The data were measured using the microscopy-based method (method A), TSA/catalyzed reporter deposition fluorescence in situ hybridization (CARD-FISH; method B), flow cytometer (method C), or image recognition (method D). UCYNs include UCYN-A, UCYN-B, and unclassified UCYNs. Heterocystous cyanobacteria include Het-1, Het-2, and Het-3.

^∗ Data are reported by data providers as depth-integrated cell abundance (unlabeled depth-integrated abundances computed from volumetric data).

Table 6 Summary of new data points of nifH gene copy abundances added to version 2 of the database. UCYNs include UCYN-A1, UCYN-A2, UCYN-B, and UCYN-C. Heterocystous cyanobacteria include Het-1, Het-2, and Het-3.

^a Data are reported by data providers as depth-integrated nifH gene copy abundances (unlabeled depth-integrated abundances computed from volumetric data). ^b rnpB gene copies were determined.

Table 7 Summary of data points of cell-specific N₂ fixation rates added to version 2 of the database. The rates were measured either by using the combination of CARD-FISH and nanoscale secondary ion mass spectrometry (nanoSIMS; method A) or via the measurements of bulk N₂ fixation rates incubated with a known number of diazotrophic cells (method B; see Sect. 2.3). Note that all the data were reported as N₂ fixation rates per cell, except for Filella et al. (2022) in which biomass-normalized rates in units of d⁻¹ were reported.

2.2 Quality control

The data of N₂ fixation rates and diazotrophic abundance in the database spanned several orders of magnitude. Extremely high rates and abundance values of both usually occurred during algal blooms, and zero values indicated that diazotrophic activity was below detection or truly absent at the sampling time and stations. The positive-value data were first logarithmically transformed and then analyzed for outliers, considering that they were approximately log-normally distributed (Figs. S1–S5). For each parameter, we used Chauvenet's criterion to identify suspicious outliers whose probability of deviation from the means is lower than $\mathrm{1}/\mathrm{2}n$, where n is the number of data points (Glover et al., 2011). Because N₂ fixation rates and diazotroph abundances in the ocean can be extremely low, this filtering only applied to data on the high side. Although these outliers (labeled in the database) could be true values, we flagged them to caution users.

2.3 Nitrogen fixation rate data

The commonly used methods for marine N₂ fixation rates include ¹⁵N₂ tracer methods and the acetylene reduction assay (Mohr et al., 2010; Montoya et al., 1996; Capone, 1993). However, in the last decade, the community has turned largely to the use of ¹⁵N₂ tracer methods. The acetylene reduction assay estimates gross N₂ fixation rates indirectly from the reduction of acetylene to ethylene. Theoretical conversion factors of 3:1 and 4:1 have been used to convert acetylene reduction rates to N₂ fixation rates (Postgate, 1998; Capone, 1993; Wilson et al., 2012), although a wide range of conversion factors from 0.93 to 56 have been reported (e.g., Mague et al., 1974; Graham et al., 1980; Montoya et al., 1996; Capone et al., 2005; Mulholland et al., 2006; Wilson et al., 2012). When using the ¹⁵N₂ tracer method, samples are incubated in seawater with ¹⁵N₂ gas; the ${}^{\mathrm{15}}\mathrm{N}/{}^{\mathrm{14}}\mathrm{N}$ ratio of particulate nitrogen is measured at the beginning and the end of the incubation to calculate the N₂ fixation rate (Capone and Montoya, 2001). Most measurements using the ¹⁵N₂ tracer method only counted the fixed N in particulate forms and ignored the N that was fixed but then excreted by diazotrophs in the form of dissolved organic N (DON) during incubation, which could theoretically be counted by the acetylene reduction assays (Mulholland, 2007). In some studies using the ¹⁵N₂ tracer method, this missing N was counted by also measuring the ¹⁵N enrichment in DON (Berthelot et al., 2017; Benavides et al., 2013a; Berthelot et al., 2015; Benavides et al., 2013b).

Compared to the ¹⁵N₂ tracer method, the acetylene reduction assay requires less incubation time. However, in addition to the uncertainty in converting ethylene production to N₂ fixation, the purity of acetylene gas, trace ethylene contamination, and the Bunsen gas solubility coefficient of produced ethylene can also affect the accuracy of estimated N₂ fixation rates (Hyman and Arp, 1987; Breitbarth et al., 2004; Kitajima et al., 2009). Acetylene used in the assay can even impact the metabolic activities of diazotrophs (Giller, 1987; Hardy et al., 1973; Flett et al., 1976; Staal et al., 2001). Moreover, the acetylene reduction assay needs to preconcentrate cells for signal detection when diazotrophic biomass is low, which may lead to underestimated N₂ fixation rates by perturbing cells during concentration and filtration (e.g., Capone et al., 2005; Barthel et al., 1989; Staal et al., 2007). In recent years, the acetylene reduction assay has undergone significant advancement. The sensitivity of ethylene detection has been improved by utilizing a reduced gas analyzer (Wilson et al., 2012) and by using highly purified acetylene gas to minimize the ethylene background (Kitajima et al., 2009). However, preparing high-purity acetylene with a low level of ethylene contamination remains a challenge. More recently, a new method named Flow-through incubation Acetylene Reduction Assays by Cavity ring-down laser Absorption Spectroscopy (FARACAS) has been introduced for high-frequency measurements of aquatic N₂ fixation (Cassar et al., 2018). This method involves continuous flow-through incubations and spectral monitoring of acetylene reduction to ethylene. By employing short-duration flow-through incubations without cell preconcentration, potential artifacts are minimized. This approach also allows for near-real-time estimates, enabling adaptive sampling strategies.

The original ¹⁵N₂ tracer method involves the addition of a known volume of ¹⁵N₂-labeled bubbles to the incubation bottle (named original ¹⁵N₂ bubble method hereafter). However, this method was found to underestimate rates because N₂ gas solubility is low and tracer additions take a long time to equilibrate (Mohr et al., 2010; Großkopf et al., 2012; Jayakumar et al., 2017). To address this issue, the ¹⁵N₂ dissolution method has been employed, which involves pre-preparing ¹⁵N₂-enriched seawater to maintain constant ¹⁵N₂ atom % enrichment throughout the incubation (Mohr et al., 2010), similar to the method described in Glibert and Bronk (1994). However, the ¹⁵N₂ dissolution method does not always yield higher N₂ fixation rates than the original ¹⁵N₂ bubble method (Table S4 in Großkopf et al., 2012; Saulia et al., 2020); it is still not conclusive what controls the magnitude of the underestimation (if it exists) in the original ¹⁵N₂ bubble method. Compared to the original ¹⁵N₂ bubble method, the ¹⁵N₂ dissolution method is more susceptible to the introduction of contaminants (e.g., metals) during the preparation of the ¹⁵N₂ inoculum due to its more complex process, which can alter the diazotrophic activities and abundance, thereby impacting the accuracy of N₂ fixation measurements (Dabundo et al., 2014; Klawonn et al., 2015). For example, Needoba et al. (2007) reported that a low but detectable amount of Fe³⁺ contamination can be measured when protecting the needle of the gas-tight syringe with commercially available tubing. Additionally, pH and other chemical properties of the inoculum may be altered during its preparation, further affecting the measurements of N₂ fixation. Despite these limitations, the ¹⁵N₂ dissolution method remains the predominant assay for measuring N₂ fixation rate due to its ability to satisfy the fundamental assumption of constant ¹⁵N₂ atom % enrichment over the incubation period.

More recently, a modified ¹⁵N₂ bubble method, known as the ¹⁵N₂ bubble release method, has been proposed as an alternative to the ¹⁵N₂ dissolution method (Klawonn et al., 2015; Chang et al., 2019; Selden et al., 2019). This method involves adding ¹⁵N₂ gas to the incubation bottles and mixing for a brief period (∼15 min) to facilitate ¹⁵N₂ equilibration and then removing the gas bubble. Compared to the original ¹⁵N₂ bubble method, the ¹⁵N₂ bubble release method ensures uniform ¹⁵N₂ atom % enrichment throughout the incubation. Moreover, it causes less interference with the incubation matrix than the ¹⁵N₂ dissolution method. However, the mixing of incubation bottles required to stimulate gas dissolution has been suggested to negatively affect diazotrophs, although no robust studies have yet been performed to assess this critique (Wannicke et al., 2018; White et al., 2020). Moreover, the ¹⁵N₂ bubble release method requires a handling step, and additional costs for preparing tracers may be another challenge for researchers (White et al., 2020). Ultimately, White et al. (2020) “advise employing either the dissolution or bubble release method, whichever is best suited to the specific research objectives and logistical constraints” with additional recommendations on the need for determination of detection limits for all rate measurements.

We compared volumetric N₂ fixation rates in the upper 50 m and depth-integrated N₂ fixation rates in the database measured using acetylene reduction assays, the original ¹⁵N₂ bubble method, and the ¹⁵N₂ dissolution method and found that they span a similar range (Fig. 1). Meanwhile, in the analysis for volumetric N₂ fixation rates in the upper 50 m, the peak of the log-normal distributions of the measurements using the ¹⁵N₂ dissolution method was approximately double that of the original ¹⁵N₂ bubble method (Fig. 1a). The measurements using the ¹⁵N₂ bubble release method were limited to several study sites and their distribution was thus not presented in this study. A further analysis comparing the original ¹⁵N₂ bubble method and the ¹⁵N₂ dissolution method will be presented later (see Sect. 4.1).

Figure 1 Distribution of N₂ fixation rates measured using acetylene (C₂H₂) reduction assays, the original ¹⁵N₂ bubble method, and the ¹⁵N₂ dissolution method. (a) Volumetric data in the upper 50 m; (b) depth-integrated data. Only rates measured with incubation periods of 24 h are shown. Please note that the bars in the plots do not represent cumulative data.

The majority of N₂ fixation rates (9405) were measured with incubation periods of 24 h and were reported as daily rates. In contrast, 2416 samples were incubated for less than 24 h and hourly N₂ fixation rates were reported. Diel cycles of N₂ fixation vary among samples and/or diazotrophic groups, and substantial errors may be introduced when extrapolating N₂ fixation rates incubated for less than 24 h to daily rates (White et al., 2020). Therefore, the N₂ fixation rates measured with incubation periods of less than 24 h were collected into separated data sheets in our database and were not used in further analysis within this study. Please note that the incubation periods of whole diurnal cycles (e.g., 24, 48, or 72 h) were used in Konno et al. (2010). The incubation of samples in Yogev et al. (2011) lasted from 24 to 30 h. The reported daily N₂ fixation rates by these two studies were also included in the 24 h data sheets and were used in our estimation of the global marine N₂ fixation rate (see below).

Figure 2 Spatial distribution of the number of volumetric and depth-integrated data points in version 2 of the diazotrophic database, binned in 1^∘ latitude × 1^∘ longitude grids. (a, b) N₂ fixation rates, (c, d) cell abundance, and (e, f) nifH gene copy abundance. The data sources include the original version of this database (Luo et al., 2012; blue diamonds), two compiled datasets (Tang et al., 2019; Tang and Cassar, 2019; orange triangles), and this study (orange circles).

Figure 3 Spatial distribution of the number of volumetric and depth-integrated data points binned in 1^∘ latitude × 1^∘ longitude grids for (a, b) N₂ fixation rates, (c, d) cell abundance, and (e, f) nifH gene copy abundance. The size of the circles represents the number of data points in each bin. The blue diamonds mark the location of outliers identified using Chauvenet's criterion.

Cell-specific N₂ fixation rates of diazotrophs (or symbioses) were mostly measured using catalyzed reporter deposition fluorescence in situ hybridization (CARD-FISH) and nanoscale secondary ion mass spectrometry (nanoSIMS), in combination with ¹⁵N₂ addition experiments (Mills et al., 2020; Berthelot et al., 2019). Using specific oligonucleotide probes, CARD-FISH enables the visualization and location of the regions of interest in diazotrophs at a single-cell level using a epifluorescence microscope. This is subsequently prepared for the secondary electron image in nanoSIMS analysis. Importantly, the handling, fixation, and processing of the samples with CARD-FISH has been demonstrated to significantly impact the enrichment measured by nanoSIMS (see Musat et al., 2014; Woebken et al., 2015; Meyer et al., 2021). The nanoSIMS technique detects the enrichment of ¹⁵N atoms in the targeted regions, allowing for the calculation of the cell-specific rate. Additionally, in one study, handpicked Trichodesmium colonies or trichomes were incubated and the measured total N₂ fixation rates were normalized to number of cells (McCarthy and Carpenter, 1979).

2.4 Estimation of the global marine N₂ fixation rate

Using these data, we performed a first-order estimation of the global marine N₂ fixation rate. In a previous study (Luo et al., 2012), version 1 was utilized to estimate the global marine N₂ fixation rate, which included all the depth-integrated N₂ fixation rates. However, in this study, we employed more rigorous criteria to estimate the global rate using both version 1 and version 2, taking into account the reliability of different N₂ fixation rate data discussed in the preceding section. Specifically, we exclusively used depth-integrated N₂ fixation rates that met the following criteria: (1) measurements were taken from whole seawater samples, (2) incubation periods of 24 h were used, and (3) the three ¹⁵N₂-based methods were employed, although we acknowledged that the rates obtained using the original ¹⁵N₂ bubble method might be underestimated. N₂ fixation rates obtained through the acetylene reduction method were excluded from this estimate due to the significant uncertainties described above.

Figure 4 Latitudinal distribution of volumetric and depth-integrated (a, b) N₂ fixation rates, (c, d) cell abundance, and (e, f) nifH gene copy abundance, including the data from version 1 of the database (blue) and the new data added to version 2 of the database (orange).

Applying these criteria, we selected 309 and 1642 depth-integrated N₂ fixation rates from version 1 and version 2, respectively. The greater number of data in version 2 potentially provided more constraints on estimating global marine N₂ fixation. We applied Chauvenet's criterion to identify outliers, using the log-transformed values of the selected data (see Sect. 2.2). As a result, two high-value outliers were removed in version 1 (one in the North Pacific and one in the South Pacific) while no outliers were detected in version 2. This difference can be attributed to the larger number of data samples in version 2, which allowed for a more relaxed threshold in identifying outliers.

The estimation of the global marine N₂ fixation rate involved four steps. First, we calculated the arithmetic or geometric means of the depth-integrated N₂ fixation rates within each 3^∘ latitude × 3^∘ longitude bin. Second, these mean values were further averaged using either arithmetic or geometric methods to determine the mean N₂ fixation rates for different ocean basins, which included the North Atlantic, South Atlantic, North Pacific, South Pacific, Indian, Arctic, and Southern oceans, as well as the Mediterranean Sea. Third, we multiplied the arithmetic or geometric mean of each basin by its respective area to estimate the total N₂ fixation rate for that specific basin, except when there was insufficient spatial coverage available. Finally, we obtained the global marine N₂ fixation rate by summing up the individual rates calculated for each basin, with the errors associated with the basin rates propagated properly (Glover et al., 2011).

In the first two steps, the geometric means were derived from positive N₂ fixation rates (NF₊): if μ and SE represented the mean and standard error of ln (NF₊), respectively, the geometric mean was e^μ. The confidence interval for the geometric mean, based on the standard error, ranged between $e^{\mathit{\mu }}/e^{\mathrm{SE}}$ and e^μ⋅e^SE (Thomas, 1979). To address the issue of not including zero-value N₂ fixation rates, we adjusted the geometric means by multiplying them with the percentage of zero-value data within each 3^∘ latitude × 3^∘ longitude bin (in the first step) or within each basin (in the second step).

2.5 Diazotrophic abundance data

Diazotroph cell abundances were determined by using standard light microscopy, and in some cases by using epifluorescence microscopy. A recent study used machine learning techniques to detect and enumerate diazotrophs in a large dataset of microscopic images (Karlusich et al., 2021). In the original database, only the cell abundances of Trichodesmium and heterocystous cyanobacteria were recorded. Version 2 also included datasets of enumerated abundance of all UCYN groups detecting them by tyramide signal amplification–FISH (TSA-FISH) using a specific DNA probe UCYN-238 (Biegala and Raimbault, 2008; Le Moal and Biegala, 2009; Le Moal et al., 2011; Riou et al., 2016). This method is also called CARD-FISH and was used specifically to enumerate UCYN-A (Martínez-Pérez et al., 2016; Biegala and Raimbault, 2008; Le Moal et al., 2011; Table 5).

Cell abundance of Trichodesmium was recorded as the number of trichomes per volume of water in our database, although it was also reported in some studies as the number of cells or colonies per volume of water. In the latter cases, the data were converted to trichomes per volume of water by using a commonly used factor of 200 (132–241) trichomes colony⁻¹ (Letelier and Karl, 1996), similar to the conversion used in the original database (Luo et al., 2012).

Figure 5 Monthly distribution of volumetric and depth-integrated (a, b) N₂ fixation rates, (c, d) cell abundance, and (e, f) nifH gene copy abundance, including the data from version 1 of the database (blue) and the new data added to version 2 of the database (orange).

Figure 6 Seasonal distribution of volumetric and depth-integrated (a, b) N₂ fixation rates, (c, d) cell abundance, and (e, f) nifH gene copy abundance, including the data from version 1 of the database (blue) and the new data added to version 2 of the database (orange). Spring: March–May in the Northern Hemisphere and September–November in the Southern Hemisphere; summer: June–August in the Northern Hemisphere and December–February in the Southern Hemisphere; autumn: September–November in the Northern Hemisphere and March–May in the Southern Hemisphere; and winter: December–February in the Northern Hemisphere and June–August in the Southern Hemisphere.

The abundance of heterocystous cyanobacterial cells was also recorded in this database. Based on the number of DDAs was reported in several studies, we assumed that two (reported range: 1–2) and five (reported range: 1–5) Richelia spp. filaments were associated with each Hemiaulus and Rhizosolenia cell, respectively (Villareal et al., 2011; Caputo et al., 2019), and that five (reported range: 3–10) Richelia rhizosolenia filaments were associated with each Chaetoceros cell (Tuo et al., 2021; Caputo et al., 2019). Richelia have terminal heterocysts, and the number of vegetative cells varies depending on the host diatom. In Hemiaulus and Chaetoceros spp. diatoms, Richelia filaments are shorter (e.g., 3–4 vegetative cells) compared to in Rhizosolenia, where Richelia filaments are longer (e.g., 5–6 vegetative cells; Foster et al., 2022b).

In measurements of nifH gene copy abundances, different qPCR or ddPCR assays were designed to target specific diazotrophic groups (Church et al., 2005a; Foster et al., 2007; Gradoville et al., 2017; Benavides et al., 2016b), mainly including Trichodesmium, UCYN subgroups (A1, A2, B, and C) and heterocystous groups (het-1, het-2, and het-3; Table 6).

All the uncertainties reported in this paper reflect one standard error of the means unless specified.

3 Results

3.1 Data distribution

Version 2 of the database significantly expanded N₂ fixation rate measurements, filling spatial gaps, particularly in the Indian Ocean and the Southern Hemisphere (Table 1; Figs. 2a and b, 3a and b). The number of depth-integrated N₂ fixation rate measurements was tripled (Table 1; Figs. 2b and 3b). The largest fraction of new data derived from inclusion of nifH gene abundances, in particular data contributions from the Pacific and Atlantic oceans (Table 3; Figs. 2e and f, 3e and f). Compared to other parameters, the new database contained only a modest increase in new cell abundances, mostly from subtropical oceans (Table 2; Figs. 2c and d, 3c and d). Overall, there remained more limited data on N₂ fixation and diazotrophic abundance in the Arctic and Southern oceans, with a number of rate measurements reporting values below detection limits.

Version 2 added data at all latitudinal ranges (Fig. 4). In particular, version 2 extended the range of data from tropical and subtropical areas to include polar regions in the Arctic Ocean (Harding et al., 2018) and Antarctic coast (Shiozaki et al., 2020).

The data in version 2 reduce the difference in the number of data points across months, especially for nifH gene copies, in which substantially more samples were collected in January and February (Fig. 5). When considering seasons in both the Northern Hemisphere and the South Atlantic and Pacific, the data were distributed more evenly (Fig. 6).

Although most of the new data were measured in near-surface waters, numerous nifH gene copy abundance data were also sampled in deeper layers in the euphotic zone (Fig. 7). Additionally, active N₂ fixation and the existence of diazotrophs were found below the euphotic zone (e.g., depth >200 m; Benavides et al., 2016a, 2018b; Selden et al., 2019; Hamersley et al., 2011; Bonnet et al., 2013; Loescher et al., 2014; Benavides et al., 2015; Fig. 7).

Figure 7 Vertical distribution number of (a) N₂ fixation rates, (b) cell abundance, and (c) nifH gene copy abundance data, including the data from version 1 of the database (blue) and the new data added to version 2 of the database (orange).

3.2 N₂ fixation rates

The volumetric N₂ fixation rates in five vertical layers and the depth-integrated N₂ fixation rates were binned in 3^∘ latitude × 3^∘ longitude bins, and the arithmetic means in each bin are displayed (Fig. 8). The depth-integrated N₂ fixation rates ranged over orders of magnitude, from 10⁻⁴–10³ $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ (mostly from 1 to 10² $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$; Fig. 8a). Some high rates (i.e., 10²–10³ $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) were found in the western Pacific Ocean, the regions near the Hawaiian Islands, and the western tropical Atlantic Ocean. Approximately 10 % of the depth-integrated N₂ fixation rates were <1 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ and were mainly from the North Atlantic and Indian oceans. Within the water column, the N₂ fixation rates were highest in the upper 25 m (Fig. 8b and c), below which the rates rapidly decreased with depth (Fig. 8d–f). In the upper 25 m, volumetric N₂ fixation rates in the southwestern Pacific were higher than those in other areas, mostly ranging from 1 to 100 $\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{d}}^{-\mathrm{1}}$. Undetectable N₂ fixation rates were reported mostly in subpolar regions, as well as in certain tropical and subtropical regions (Fig. 8).

Figure 8 N₂ fixation rates in version 2 of the database. The panels show (a) depth-integrated data and volumetric data at (b) 0–5 m, (c) 5–25 m, (d) 25–100 m, (e) 100–200 m, and (f) below 200 m. For a clear demonstration, arithmetic mean N₂ fixation rates in 3^∘ latitude × 3^∘ longitude bins are shown. Zero-value data are denoted as black empty circles. Only rates measured with incubation periods of 24 h are included.

Cell-specific N₂ fixation rates span a range from 10⁻⁴ to 10³ fmol N cell⁻¹ d⁻¹, although mostly on the order of 10⁻² to 10² fmol N cell⁻¹ d⁻¹ (Fig. 9). The mean cell-specific N₂ fixation rates of Trichodesmium, UCYN-A2, and heterocystous cyanobacteria were 1 to 2 orders of magnitude higher than those of other diazotrophic groups (Fig. 9 and Table S1).

Figure 9 Violin plot of cell-specific N₂ fixation rates, including measurements for Trichodesmium, UCYN-A, UCYN-A1, UCYN-A2, UCYN-B, heterocystous cyanobacteria, unclassified cyanobacteria, and NCDs. The range of each box spans the 25th–75th percentile of data, the black line in each box is the median, and the red dot represents the arithmetic mean.

3.3 Diazotrophic abundance

The depth-integrated cell abundances and volumetric cell abundances in the upper 25 m are also shown as the arithmetic means in 3^∘ latitude × 3^∘ longitude bins (Fig. 10). Trichodesmium abundance generally decreased from the west to the east in the Atlantic Ocean (Fig. 10a and b). In the Pacific Ocean, Trichodesmium appeared more abundant in the west. The abundance data of heterocystous diazotrophs were still scarce (Fig. 10c and e). The volumetric cell-count-based abundance data are also displayed in three additional depth intervals (Fig. S6).

Figure 10 Depth-integrated cell abundances and volumetric cell abundances in the upper 25 m in version 2 of the database. The panels show (a, b) Trichodesmium, (c, d) het-1/2, and (e, f) het-3. For a clear demonstration, data are binned to 3^∘ latitude × 3^∘ longitude, and arithmetic means in each bin are shown. Zero-value data are denoted as open black circles.

Gene copies of nifH had better spatial coverage than the cell-count data (Fig. 11). Depth-integrated Trichodesmium nifH copies were also more abundant in the western Pacific and western Atlantic oceans (Fig. 11a). Some high depth-integrated nifH abundance of UCYN-A and UCYN-B were also reported in the northwestern and southwestern Pacific Ocean (Fig. 11c and e). High nifH abundances of Richelia were found in the southwestern Pacific Ocean and western Atlantic oceans (Fig. 11i). The nifH abundance data for UCYN-C and het-3 were sparse. The volumetric nifH abundance data are displayed in three depth intervals (Figs. 11 and S7). Almost all diazotrophs were more abundant in the upper 25 m than in deeper water.

Figure 11 Volumetric and depth-integrated nifH gene copy abundances in version 2 of the database. For volumetric abundances, only data in the upper 25 m are shown. The panels show gene copy abundances of (a, b) Trichodesmium, (c, d) UCYN-A, (e, f) UCYN-B, (g, h) UCYN-C, (i, j) het-1 + het-2, (k, l) het-3, and (m, n) Gamma A (an NCD phylotype). Depth-integrated data for Gamma A are not available. For a clear demonstration, data are binned to 3^∘ latitude × 3^∘ longitude and arithmetic means in each bin are shown. Zero-value data are denoted as open black circles.

3.4 First-order estimate of global oceanic N₂ fixation rate

Compared to version 1, the spatial coverage of data in version 2, in terms of the fraction of 3^∘ latitude × 3^∘ longitude bins, was greatly increased in all ocean basins (Table 8). The spatial data coverage was very low in the Southern and Arctic oceans (1 % and 2 % of total bins, respectively; Table 8), and we therefore did not estimate total N₂ fixation rates for these two basins. Please note that the inaccurate areas of the North and South Pacific oceans used in estimating the global oceanic N₂ fixation rate by Luo et al. (2012) was corrected in this study (Table 8).

Table 8 First-order estimates of N₂ fixation rates based on their arithmetic means in different ocean basins. Data are first binned to 3^∘ latitude × 3^∘ longitude grids before being used to calculate arithmetic means in each basin. The arithmetic means are multiplied by the basin areas to calculate the N₂ fixation rates of each basin. NQ: not quantified due to limited data points. ND: no data. The values in the parentheses are the percentages of 3^∘ × 3^∘ bins in each basin that have measurements. The reported uncertainties are one standard error of the mean.

We first compared the N₂ fixation rates estimated based on arithmetic means of version 1 and version 2 (Table 8). Using available data in version 2, the global N₂ fixation rate was determined to be 223±30 Tg N yr⁻¹, which was 3 times that obtained from version 1 (Table 8). The substantial increase was mostly driven by notable changes in the South Pacific, North Atlantic, and Indian oceans. In the South Pacific Ocean, numerous high N₂ fixation rates were observed in the western subtropical region over the past decade (Fig. 12), resulting in a substantial increase of 68±23 Tg N yr⁻¹ in the estimated N₂ fixation rate for this basin (Table 8). It is worth noting that these newly recorded measurements in the western subtropics of the South Pacific Ocean might even be underestimated since most of them were obtained using the original ¹⁵N₂ bubble method. In the North Atlantic Ocean, the estimated N₂ fixation rate also experienced an increase of 30±9 Tg N yr⁻¹ for (Table 8), without any discernible pattern regarding the locations of the new high N₂ fixation measurements (Fig. 13). Furthermore, in the Indian Ocean, the improved data coverage in version 2 (Fig. 8a) supported the estimation of an N₂ fixation rate of 35±14 Tg N yr⁻¹ for this basin (Table 8), which was not possible to calculate using version 1 due to insufficient data availability.

Figure 12 Depth-integrated N₂ fixation rates in the South Pacific Ocean ($\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). The shown data are arithmetic mean rates in 3^∘ latitude × 3^∘ longitude bins. Empty diamonds and filled circles denote the existing data in version 1 of the database and the new data added to version 2, respectively.

Figure 13 Depth-integrated N₂ fixation rates in the North Atlantic Ocean ($\mathrm{\mu }\mathrm{mol}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$). The shown data are arithmetic mean rates in 3^∘ latitude × 3^∘ longitude bins. Empty diamonds and filled circles denote the existing data in version 1 of the database and the new data added to version 2, respectively.

However, when estimating the global marine N₂ fixation rate using geometric means, both version 1 and version 2 yielded similar rates of approximately 50 Tg N yr⁻¹ (Table 9). The N₂ fixation rates in each basin tended to follow a log-normal distribution (Fig. 14), with the geometric mean aligning near the peak of the distribution. In the South Pacific Ocean, as discussed earlier, version 2 included a substantial number of newly observed high N₂ fixation rates, but it also incorporated a significant number of rates that were much lower than those in version 1 (Fig. 14c). This could be partially attributed to enhanced detection limits in measurements. Consequently, while version 2 yielded a much higher arithmetic mean N₂ fixation rate compared to version 1 for the South Pacific Ocean (Table 8), their geometric means remained quite similar (Table 9). In the North Pacific Ocean, for the same reasons, the arithmetic mean N₂ fixation rates obtained from both versions were very close, while the geometric mean of version 1 was even higher than that of version 2 (Tables 8 and 9; Fig. 14a). These analyses reveal that, despite the similarity in geometric means of N₂ fixation rates obtained from both versions of the database, the higher arithmetic means in version 2 were not coincidental. Instead, they were the direct outcome of the improved measurement methods and the expanded spatial and temporal coverage of marine N₂ fixation over the past decade. Consequently, previous assessments of the global marine N₂ fixation rate were likely underestimated due to the absence of these new measurements.

Table 9 Same as Table 8 but based on the geometric means of N₂ fixation rates. The numbers in parentheses are estimated ranges based on one standard error of log-transformed N₂ fixation rates (see Sect. 2.4).

Figure 14 Comparison of the distribution of log-transformed N₂ fixation rates between the two versions of the database. Note that the zero-value data are not included because of the log transformation. The comparison is performed for data in (a) North Pacific, (b) North Atlantic, (c) South Pacific, and (d) South Atlantic oceans.

We must emphasize that this calculation simply used the average N₂ fixation rates in different ocean basins; therefore, our calculation can only be considered a first-order estimate. Furthermore, limited measurements have shown a large range of N₂ fixation rates in the Southern Ocean (Fig. 8). Considering its vast area, future measurements expanding coverage of N₂ fixation rates in the Southern Ocean (see White et al., 2022) may help to better constrain the contribution of N₂ fixation to the N budget of the global ocean. The new database presented here also expands opportunities for improved statistical estimates of N₂ fixation patterns and global rates based on the modeling of environmental controls (Luo et al., 2014).

Figure 15 Comparison of measured N₂ fixation rates using the original ¹⁵N₂ bubble method and the ¹⁵N₂ dissolution method. The pink dots are measurements. The fitted results of the two methods by the generalized additive model (GAM) and confidence intervals are represented by the red solid line and the dashed black lines, respectively. Only the N₂ fixation rates measured with incubation periods of 24 h were included in this analysis.

4 Discussion

4.1 Comparison of N₂ fixation measured using ¹⁵N₂ bubble and dissolution methods

To date, the origin of the discrepancy in the N₂ fixation rates estimated using different ¹⁵N₂ tracer methods remains unclear. As shown above, the volumetric N₂ fixation rates obtained by the original ¹⁵N₂ bubble method and the ¹⁵N₂ dissolution method spanned a similar range (Fig. 1), while the average rates using the former method were significantly lower than that measured using the latter method (one-tailed Wilcoxon test, p<0.001, n=2460 and 1128). With substantial data accumulated over the past decade, we further compared N₂ fixation rates measured using the two methods at close locations and sampling time, although the samples were not identical. We first binned data collected from the same months, horizontal locations (3^∘ latitude × 3^∘ longitude), and depth intervals (0–5, 5–25, 25–100, and 100–200 m) and calculated the average rates for each method in each bin. The results showed that the original ¹⁵N₂ bubble method produced lower rates than the ¹⁵N₂ dissolution method in 69 % of the cases (Fig. 13). Furthermore, our analysis employing the generalized additive model (GAM) revealed that the relationship between the rates measured using the original ¹⁵N₂ bubble method and those obtained through the ¹⁵N₂ dissolution method closely adhered to the 1:1 line, albeit with slightly lower values in the former (Fig. 15). Please note that these slightly lower values can still result in significant underestimation in measured N₂ fixation rates because the GAM model was applied in a logarithmic space. It is crucial to reiterate that the rates being compared were derived from different samples, emphasizing the necessity for future investigations that directly compare the two methods using the same samples with controlled parameters such as temperature, volume of injected ¹⁵N₂, and incubation volume. Despite this limitation, our analysis suggests that the extensive body of historical marine N₂ fixation rate data obtained through the original ¹⁵N₂ bubble method is still valuable, particularly in the examination of spatial and temporal variations in N₂ fixation.

We also used the same procedure to compare the N₂ fixation rates measured using acetylene reduction assays and the ¹⁵N₂ tracer methods. However, there were insufficient pairs of data available for reliable comparisons (n=16 for acetylene reduction versus the ¹⁵N₂ dissolution method; n=6 for acetylene reduction versus original ¹⁵N₂ bubble method).

4.2 Comparison between diazotrophic cell counts and nifH copies

Whether or not nifH copies can be used to infer diazotrophic abundance and to study diazotrophic biogeography, some still challenges remain in the conversion of gene counts to biomass, as a large range in the number of nifH copies per diazotrophic cell has been reported (Table S2). In version 2, we first converted Trichodesmium trichome abundance to cell abundance using the same conversion factor of 100 cells trichome⁻¹ as that used in Luo et al. (2012). This conversion resulted in the mean and variance of log₁₀-transformed Trichodesmium cell abundance (10^6.5±1.3 cells L⁻¹) very similar to that those Trichodesmium nifH gene copies (10^6.6±1.5 copies L⁻¹; Fig. 16a). More recently, however, a much lower conversion factor of 13.2±2.3 cells trichome⁻¹ was suggested for Trichodesmium based on larger sample sizes, although a very large range of 1.2–685 cells trichome⁻¹ was reported (White et al., 2018). Hence, when a conversion factor of 10 cells trichome⁻¹ was applied, the Trichodesmium nifH gene copy abundance was 1 order of magnitude higher than its cell abundance (Fig. 16a). This result was within the reported mean nifH : cell ratios for Trichodesmium, albeit based on sparse samples, on the order of 10–100 (Table S2). It is worth noting that there have been suggestions that the observed nifH : cell ratio for Trichodesmium may be overestimated due to methodological limitations (Gradoville et al., 2022). Our analyses underscore the importance of enumerating Trichodesmium cells, rather than solely focusing on trichomes, in correctly evaluating Trichodesmium abundance, which has been suggested for future studies by White et al. (2018). While counting all Trichodesmium cells may be impractical, it would be valuable to report the number of cells in random samples of Trichodesmium trichomes.

Figure 16 Comparison of all cell-count and nifH gene copy abundance data in the database. The box plots show the median (central line), 25th and 75th percentiles (upper and lower edges of the boxes), 5th and 95th percentiles (error lines), and outliers (red crosses) of the log₁₀-transformed data. The comparisons are conducted for (a) Trichodesmium, (b) het-1/2, and (c) het-3. Note that the two conversion factors of 10 and 100 cells trichome⁻¹ are used for Trichodesmium.

The same analyses for heterocystous cyanobacteria showed that the nifH gene copy abundances were approximately 2 orders of magnitude greater than the cell abundances in terms of both mean and distribution (Fig. 16b and c). It must be noted that this simple analysis used all the data in our database. The limited in situ measurements for identical samples resulted in a mean nifH : cell ratio of 76 for heterocystous cyanobacteria (Table S2), consistent with our simple analysis.

In contrast, much lower nifH : cell ratios (1.51–2.58) were derived from regression analysis for heterocystous cyanobacteria and UCYN-B collected in the subtropical North Pacific (Gradoville et al., 2022). Considering these overall scarce measurements and the outcomes of our analysis, it is plausible that there is substantial variability in nifH : cell ratios. We expect that future studies, focusing on constraining these ratios and identifying mechanisms underlying variability in these ratios, will contribute to a more comprehensive understanding of the connection between nifH gene counts and diazotrophic cell abundance.

The application of qPCR assays for nifH-based abundance (DNA) and expression (RNA) has emerged as a critical step forward in our understanding of the distribution, abundance, and physiology (e.g., expression of nifH) of diazotrophs (Short and Zehr, 2005; Zehr and Riemann, 2023). Previously, estimating the abundances of diazotrophs was limited to those that could be identified by microscopy, e.g., Trichodesmium, heterocystous cyanobacteria (e.g., Richelia, Calothrix, Anabaena, Nodularia, Aphanizomenon), and some unicellulars (e.g., Cyanothece, later Crocosphaera). Thus, qPCR enabled the study of diazotrophic targets (and their activity) without the need for microscopy to identify them, which came later as some diazotrophs did (and still do) require application of FISH techniques for identification (Biegala and Raimbault, 2008). Additionally, qPCR allowed the study of in situ activity (gene expression) by diazotrophs without the need for cultivation. Although beyond the scope of the work presented here, important considerations should be taken into account when applying microscopy and qPCR datasets (Table S3), for example, to biogeochemical models (Meiler et al., 2023).

4.3 Biomass conversion factor

For possible further usage of cell-counted abundance data, here we suggest carbon biomass conversion factors for different diazotrophic groups (Tables 10 and S4). Most biomass conversion factors suggested here are the same as those used in Luo et al. (2012), excluding UCYN-A and heterocystous cyanobacteria, where new information has become available or additional consideration is necessary. A recent study has discovered a new symbiosis association between the unicellular diazotroph (UCYN-C) and diatom Epithemia strains (Schvarcz et al., 2022). However, the conversion factor of UCYN-C could not be updated in this study due to insufficient information on the biovolumes of host cell.

The conversion factor for UCYN-A was updated because it has been found to live symbiotically with haptophyte Braarudosphaera bigelowii and relatives (Thompson et al., 2012; Hagino et al., 2013). Because the host and UCYN-A should function together, the host biomass is allocated to UCYN-A. It has been reported that each haptophyte cell hosts one UCYN-A1 cell (Cornejo-Castillo et al., 2019) or one UCYN-A2 cell (Suzuki et al., 2021). We used the empirically derived equation (Verity et al., 1992)

$$\begin{array}{}\text{(1)}& C=\mathrm{0.433}\times V^{\mathrm{0.863}}\end{array}$$

to estimate the biomass of UCYN-A and their hosts. The biomasses of a UCYN-A1 cell with a diameter of 1 µm and a UCYN-A2 cell with a diameter of 1.6–3.3 µm (Cornejo-Castillo et al., 2019; Martínez-Pérez et al., 2016) are 0.2 and 0.8–5.5 pg C, respectively. The biomasses of the host cell for UCYN-A1 or UCYN-A2 is 1.5–2.2 or 6.8–43 pg C according to their reported cell diameters (2–2.3 or 3.6–7.3 µm), respectively (Martínez-Pérez et al., 2016; Cornejo-Castillo et al., 2019). Hence, the biomasses of the UCYN-A1 and the UCYN-A2 symbioses are 1.7–2.4 and 7.6–48 pg C, respectively. After normalizing the symbiotic biomass to the number of UCYN cells in each symbiosis (1 for both UCYN-A1 and UCYN-A2), the biomass conversion factors are 1.7–2.4 pg C (UCYN-A1 cell)⁻¹ and 7.6–48 pg C (UCYN-A2 cell)⁻¹.

Because heterocystous cyanobacteria and their host diatoms form DDAs, similar to UCYN-A, we also suggest allocating the biomass of host diatoms to each associated diazotrophic cell (Table S4). The biomasses of heterocystous cells and vegetative cells in Richelia filaments were updated according to the cell dimension data reported in Caputo et al. (2019) using the same empirical equation above. The carbon biomass of host diatom cells was calculated using an empirical equation (Menden-Deuer and Lessard, 2000):

$$\begin{array}{}\text{(2)}& C=\mathrm{0.117}\times V^{\mathrm{0.881}},\end{array}$$

where C is the diatom cell carbon biomass (pg C cell⁻¹) and V is the average cell biovolume (µm³) of each diatom genus, for which values from a database (Harrison et al., 2015) were used in this study (Table S4). Each host diatom associates with multiple heterocysts. The numbers of Richelia heterocysts associated with Hemiaulus, Rhizosolenia, and Chaetoceros were observed to be within the range of 1–2, 1–5, and 3–10 respectively (Villareal et al., 2011; Yeung et al., 2012; Caputo et al., 2019); we selected both the maximum and minimum to do the estimation. The number of vegetative cells in each heterocyst was also updated according to Caputo et al. (2019). Conversion factors for DDAs were estimated by dividing the total biomass of each DDA by the number of associated heterocysts. Changes in the number of Richelia in Rhizosolenia (1 or 5) would make a large variation in its conversion factor, possibly due to the large host biomass; therefore, we keep them both to let users take caution when using this conversion factor. The resulting biomass conversion factors of Richelia–Hemiaulus and Richelia–Chaetoceros associations were estimated to be 280 pg C heterocyst⁻¹ (range: 150–1250) and 430 pg C heterocyst⁻¹ (range: 10–1900), respectively (Table S4), as the number of filaments did not have a large impact on the conversion factors.

It is important to reiterate that these biomass conversion factors are only applicable to cell-count data. Attempting to convert nifH gene copies to biomass is not recommended due to significant uncertainties associated with nifH : cell ratios, as previously discussed.

Table 10 Recommended carbon biomass conversion factors and their likely ranges for diazotrophic groups.

5 Data availability

The database is available in a data repository (https://doi.org/10.6084/m9.figshare.21677687, Shao et al., 2022).

6 Conclusions

In this study, we updated the global oceanic diazotrophic database by Luo et al. (2012) by adding new measurements reported in the past decade. Although the spatial coverage of the data was greatly expanded by this effort, the data distribution is still uneven, with most measurements reported from the Pacific and Atlantic oceans. Using the updated database, the estimation of global oceanic N₂ fixation based on arithmetic rates in ocean basins was increased from 74±7 to 223±30 Tg N yr⁻¹. This change is largely attributable to a new estimate for the Indian Ocean and a much elevated estimate for the South Pacific Ocean, the latter of which would account for ∼40 % of global N₂ fixation. This high estimation for the South Pacific Ocean is in line with its qualification as a hotspot for diazotrophy (Messer et al., 2016; Bonnet et al., 2017), partly due to iron fertilization processes in this region (Bonnet et al., 2023). Due to data sparsity, our updated estimation did not include N₂ fixation in the Southern and Arctic oceans. Furthermore, data were more concentrated in surface seawater, and a significant amount of data were measured with incubation periods shorter than a daily cycle (24 h), limiting reliable evaluations of depth-integrated N₂ fixation rates. Although this result suggests more balanced N inputs and losses in the global ocean than the previous estimate suggested, large uncertainties still exist. We also compared the N₂ fixation rates measured using the addition of a bubble of labeled gas or the addition of dissolving ¹⁵N₂ gases reported at the same location and month (not necessarily in identical samples). The results indicated that the original ¹⁵N₂ bubble method produces lower rates than the ¹⁵N₂ dissolution method in 69 % of the cases. These results reveal that, despite decades of effort, the ocean is still undersampled in terms of the distribution of diazotrophs and N₂ fixation rate measurements. Our analyses suggest that prioritizing N₂ fixation measurements in the South Pacific Ocean, Indian Ocean, and high northern latitudes can significantly reduce the current uncertainty of N₂ fixation rates in the global ocean. Nevertheless, we believe that this updated diazotrophic database, supplemented with enhanced data from the past decade, is timely and can be helpful to scientists studying the marine biogeochemical cycle of N.