Controlling factors on the global distribution of a representative marine non-cyanobacterial diazotroph phylotype (Gamma A)

Non-cyanobacterial diazotrophs may be contributors to global marine N₂ fixation, although the factors controlling their distribution are unclear. Here, we explored what controls the distribution of the most sampled non-cyanobacterial diazotroph phylotype, Gamma A, in the global ocean. First, we represented Gamma A abundance by its nifH quantitative polymerase chain reaction (qPCR) copies reported in the literature and analyzed its relationship to climatological biological and environmental conditions. There was a positive correlation between the Gamma A abundance and local net primary production (NPP), and the maximal observed Gamma A abundance increased with NPP and became saturated when NPP reached ∼ 400 mg C m⁻² d⁻¹. Additionally, an analysis using a multivariate generalized additive model (GAM) revealed that the Gamma A abundance increased with light intensity but decreased with increasing iron concentration. The GAM also showed a weak but significant positive relationship between Gamma A abundance and silicate concentration, as well as a substantial elevation of Gamma A abundance when the nitrate concentration was very high (≳ 10 µM). Using the GAM, these climatological factors together explained 43 % of the variance in the Gamma A abundance. Second, in addition to the climatological background, we found that Gamma A abundance was elevated in mesoscale cyclonic eddies in high-productivity (climatological NPP > 400 mg m⁻² d⁻¹) regions, implying that Gamma A can respond to mesoscale features and benefit from nutrient inputs. Overall, our results suggest that Gamma A tends to inhabit ocean environments with high productivity and low iron concentrations and therefore provide insight into the niche differentiation of Gamma A from cyanobacterial diazotrophs, which are generally most active in oligotrophic ocean regions and need a sufficient iron supply, although both groups prefer well-lit surface waters. More sampling on Gamma A and other non-cyanobacterial diazotroph phylotypes is needed to reveal the controlling mechanisms of heterotrophic N₂ fixation in the ocean.

1 Introduction

Dinitrogen (N₂) fixation, mostly conducted by prokaryotic bacteria (termed “diazotrophs”), is an important bioavailable nitrogen (N) source to the ocean (Moore et al., 2018; Karl et al., 2002). Although autotrophic cyanobacteria have been recognized as important diazotrophs in the ocean (Zehr, 2011), non-cyanobacteria diazotrophs (NCDs) that are presumably heterotrophic (probably including photoheterotrophic) bacteria (Bombar et al., 2016) have been widely detected (e.g., Moisander et al., 2008; Langlois et al., 2008; Halm et al., 2012; Moisander et al., 2014; Shiozaki et al., 2014) and sometimes even found to dominate the diazotrophic gene pools in surface oceans (Farnelid et al., 2011). For example, NCD nifH (a gene-encoding subunit of the nitrogenase enzyme) amplicons had higher relative abundances than autotrophic diazotrophs at some sampling sites in the South Pacific Ocean (Halm et al., 2012; Moisander et al., 2014), Indian Ocean (Shiozaki et al., 2014; Wu et al., 2019) and South China Sea (Ding et al., 2021). Metagenomic studies also revealed the abundant presence of diverse N₂-fixing Proteobacteria in ocean genomic databases (Delmont et al., 2018, 2021). Additionally, nifH of NCDs was also detected in subphotic seawaters (Benavides et al., 2018) and oxygen-deficient zones (Jayakumar and Ward, 2020; Loescher et al., 2014) where nitrogen loss was considered significant (Lam and Kuypers, 2011). Although the N₂ fixed by NCDs has not been quantified, substantial N₂ fixation found in aphotic zones (Rahav et al., 2013; Bonnet et al., 2013) and in experiments with photosynthetic inhibitors (Rahav et al., 2015; Geisler et al., 2020), as well as recovered transcripts of the NCD nifH gene (Fernandez et al., 2011; Gradoville et al., 2017), provides a line of indirect evidence of heterotrophic N₂ fixation in the ocean.

The major known NCD classes include bacteria such as Alphaproteobacteria, Gammaproteobacteria, Epsilonproteobacteria, Betaproteobacteria and Firmicutes belonging to Cluster I of nifH clusters and some obligate anaerobic bacteria and archaea belonging to Cluster III of nifH clusters (Zehr et al., 2003; Riemann et al., 2010). Among them, Gamma A is the most sampled and studied phylotype. Gamma A represents a part of uncultured Gammaproteobacterial sequences isolated from the open ocean, and its cluster is distantly related to cultured NCD (Langlois et al., 2015). Gamma A nifH gene expression has been widely found in the global ocean, suggesting its important role in marine N₂ fixation (Bird et al., 2005; Moisander et al., 2014; Langlois et al., 2015; Shiozaki et al., 2017).

It is unclear what controls the growth and distribution of NCDs as most of them, including Gamma A, are uncultivated (Bombar et al., 2016). NCDs can be different from their autotrophic counterparts in whether they are heterotrophic or photoheterotrophic, depending on organic matter as their carbon and energy source, which can be supported by experimental evidence that N₂ fixation is stimulated by adding dissolved organic matter (DOM) (Rahav et al., 2016, 2015; Bonnet et al., 2013; Bentzon-Tilia et al., 2015). However, DOM addition has sometimes not stimulated nifH expression of Gamma A even when its DNA copies have been ambient (Benavides et al., 2018), implying that DOM may not always stimulate the activity of Gamma A. Due to sensitivity to O₂ and high energy requirements of N₂ fixation (Bombar et al., 2016), abundant NCDs were found to associate with particles that supposedly provide diazotrophs with a microenvironment with depleted oxygen and rich organic matter (Riemann et al., 2010; Farnelid et al., 2010; Scavotto et al., 2015; Pedersen et al., 2018; Geisler et al., 2019). NCDs were also detected in diatom mats (Martínez et al., 1983), implying another novel habitat for NCDs. An isolated strain of diazotrophic Alphaproteobacteria from the Baltic Sea found to contain photosynthetic genes (Bentzon-Tilia et al., 2015) may complicate this issue, casting doubt on whether NCDs can be mixotrophic and also depend on light.

Although dissolved inorganic nitrogen (DIN) is generally considered to inhibit marine N₂ fixation (Karl et al., 2002; Zehr and Kudela, 2011), substantial presence of NCDs is found in DIN-replete environments such as estuaries (Geisler et al., 2020), coastal zones (Li et al., 2020), upwelling regions (Geisler et al., 2020; Moreira-Coello et al., 2017; Dekaezemacker et al., 2013) and other eutrophic seas (Bird and Wyman, 2013). Culture experiments showed that the DIN inhibition effect on NCDs can be strain specific (Bentzon-Tilia et al., 2015; Martínez-Pérez et al., 2018). Temperature could be another factor controlling NCDs, which may prefer warm oligotrophic surface oceans (Langlois et al., 2015; Shiozaki et al., 2017), the same region where the majority of autotrophic N₂ fixation occurs (Wang et al., 2019; Luo et al., 2014). Similarly to cyanobacterial diazotrophs, phosphate can also limit the growth of NCDs in oligotrophic environments (Rahav et al., 2015). Regarding other important factors that control autotrophic diazotrophs, iron (Fe) may potentially impact NCDs if they also depend on the high-Fe-containing nitrogenases to fix N₂ (Bombar et al., 2016), although, as discussed above, the N₂ fixation by NCDs is still not quantified. Strong stratification may also benefit NCDs by causing organic matter to accumulate in the upper water column (Langlois et al., 2015). However, to our knowledge, no studies have analyzed the effects of Fe or stratification on NCDs.

Mesoscale eddies may also impact NCD abundance. Although anticyclonic eddies have generally been considered to benefit autotrophic diazotrophs by inhibiting vertical DIN input from deep waters (Liu et al., 2020; Fong et al., 2008; Church et al., 2009), a class of NCDs, Gammaproteobacteria, has been found to dominate diazotrophic communities inside cyclonic eddies in the South China Sea (Zhang et al., 2011; Liu et al., 2020). Different types of mesoscale eddies may have discrepancies in impacting the ecophysiology of NCDs and their autotrophic counterparts (Benavides and Robidart, 2020).

Langlois et al. (2015) analyzed the distribution of the Gamma A phylotype in the Pacific and Atlantic oceans and suggested that Gamma A prefers warm and oligotrophic surface oceans. With more data becoming available in recent years, we collected, to the best of our knowledge, all the reported in situ measurements of Gamma A nifH copies using quantitative polymerase chain reaction (qPCR) assays, compiling a dataset with 80 % more data than those used in Langlois et al. (2015). We then analyzed the relationship between this nifH-based Gamma A abundance and the long-term background of ecological and environmental factors by using their climatological monthly averages. In addition to temperature and concentrations of nitrate, phosphate and silicate that were used in Langlois et al. (2015), we included five more variables (primary production, Fe, DOC concentrations, solar radiation and mixed-layer depth) to more thoroughly analyze potential controlling factors on Gamma A. We further explored the influence of mesoscale eddies on Gamma A abundance. Our analyses suggested that local primary productivity, temperature, dissolved Fe concentration and the occurrence of cyclonic eddies can be the main factors impacting the distribution of Gamma A in the global ocean.

2 Methods

2.1 Data summary and quality control of Gamma A abundance

A total of 1861 in situ measurements of nifH copies of Gamma A in the Pacific, Atlantic and Indian oceans were collected from 21 published papers (Table 1) and are available in a data repository (https://doi.org/10.6084/m9.figshare.17284517) (Shao and Luo, 2021). Gamma A was sometimes also named γ-24774A11 in the collected papers (Moisander et al., 2008). All these data were measured using qPCR. Note that the primer of Gamma A used by Langlois et al. (2015) in the North Atlantic was slightly different from in other studies. Most samples (88 %) were collected in the upper 100 m of the water column. In the following analyses, we represented Gamma A abundance using its nifH copies, although we note that variations in nifH copies in different cyanobacterial diazotrophic cells have been reported (White et al., 2018; Sargent et al., 2016) and nifH copy numbers in Gamma A genome remain unknown.

The non-zero nifH-based abundance data of Gamma A were approximately log-normally distributed (Fig. S1 in the Supplement). There were 748 data points reporting zero nifH copies, which theoretically could indicate that either Gamma A in the samples was truly absent or its abundance was below the detection limit. As the reported detection limit of qPCR usually ranges from 10¹ to 10² copies L⁻¹, the number of the Gamma A nifH data that could be below detection in our dataset, according to the log-normal distribution of observed non-zero data, was very likely less than 72, even assuming a high detection limit of 10² copies L⁻¹ (Fig. S1). The fact that there were far more zero-value data (748) in our dataset indicated that a high fraction of the zero-value data could represent true absence of Gamma A.

The zero-value abundance data of Gamma A were not included in our further analyses, mainly for two reasons. First, the fact that Gamma A was absent in many samples, as well as the spatially mixed distribution of the zero-value and non-zero Gamma A abundance data (see “Results and discussion”), indicated the patchy distribution of Gamma A, which was also widely found for other diazotrophs as a consequence of lateral transport and mixing of water masses (Robidart et al., 2014).

The patchiness of Gamma A implicated that it could be either present or absent even when the environmental conditions were suitable for its growth. It can also indicate our limited understanding of environmental conditions: the currently available environmental data do not include all the controlling factors of Gamma A. Nevertheless, if the zero-value data were included, similar environmental conditions could possibly be associated with both high abundance and zero abundance of Gamma A (Fig. S2), which would bias the response function of our statistical analyses, particularly as the fraction of the zero-abundance data was large (ca. one-third). Second, it is difficult to identify whether the zero-value data represented true absence or below-detection abundance of Gamma A, considering that the accuracy of qPCR was highly sensitive to sample preservation, extraction protocol and the reliance of the standard curve (Smith and Osborn, 2009).

Chauvenet's criterion was used to identify outliers in the non-zero Gamma A abundance data by first log-transforming all the data (Glover et al., 2011). Two outliers (0.22 and 0.33 copies L⁻¹) were removed because their probability of deviation from the mean was smaller than $\mathrm{1}/\left(\mathrm{2}n\right)$, where n is the number of data. Even though they can be reliable, we excluded them from the analyses to avoid possible biases.

Table 1 Data source of complied Gamma A nifH samples.

2.2 Environmental and ecological parameters

Monthly climatological environmental and ecological parameters were used as predictors for Gamma A abundance (Table 2). Temperature and concentrations of nitrate, phosphate and silicate were the products of the World Ocean Atlas (WOA) 2018 (https://www.nodc.noaa.gov, last access: 11 June 2022) (Boyer et al., 2018), and excess phosphate (P^∗) was derived from concentrations of phosphate and nitrate based on the Redfield ratio (P^∗ = [phosphate] − [nitrate]$/$16). Dissolved iron (Fe) concentrations were obtained from the Community Earth System Model (CESM) – biogeochemistry module (Misumi et al., 2014). Dissolved organic carbon concentration used a product estimated by an artificial neural network (Roshan and DeVries, 2017). Mixed-layer depth (MLD) was downloaded from Ifremer (https://www.ifremer.fr/, last access: 11 June 2022) using the criterion that the potential density of water parcels at the depth was 0.03 kg m⁻³ higher than that at the surface (De Boyer Montégut et al., 2004). Net primary production used satellite data based on the VGPM algorithm (Behrenfeld and Falkowski, 1997) (http://sites.science.oregonstate.edu/ocean.productivity/, last access: 11 June 2022). Surface photosynthetically active radiation (PAR) was downloaded from the MODIS Aqua program (https://oceancolor.gsfc.nasa.gov, last access: 11 June 2022). To estimate the vertical profile of PAR, we first obtained the estimated euphotic zone depth Z_e (https://oceancolor.gsfc.nasa.gov/, last access: 11 June 2022) at 1 % surface PAR based on an inherent optical property (IOP)-centered approach (Lee et al., 2005) and used it to estimate the attenuation coefficient:

$$\begin{array}{}\text{(1)}& k_d={\displaystyle \frac{\ln \left(\mathrm{0.01}\right)}{Z_{\mathrm{e}}}}.\end{array}$$

The PAR at depth z can be calculated assuming organisms in the mixed layer were exposed to PAR homogenously:

$$\begin{array}{}\text{(2)}& \mathrm{PAR}\left(z\right)=\left\{\begin{array}{ll}{\mathrm{PAR}}_{\mathrm{0}}{e}^{-k_{d}z}& \left(z>\mathrm{MLD}\right)\\ \frac{\mathrm{1}}{\mathrm{MLD}}{\int }_{\mathrm{0}}^{\mathrm{MLD}}{\mathrm{PAR}}_{\mathrm{0}}{e}^{-k_{d}z}\mathrm{d}z& \left(z<\mathrm{MLD}\right)\end{array},\right.\end{array}$$

where PAR₀ is the surface PAR.

To identify whether the Gamma A abundance was sampled in cyclonic or anticyclonic eddies, we extracted the satellite-merged daily sea level anomaly (SLA) from the AVISO program (https://www.aviso.altimetry.fr, last access: 11 June 2022) for the sampling days of the Gamma A data. The cores of mesoscale eddies were identified by the outermost closed contour lines of the SLA field. Only those sampling points located in the cores of cyclonic (negative SLA) or anticyclonic (positive SLA) eddies were recorded. Otherwise, data points were recorded as “outside eddy”.

All the variables used in the analyses are available in a data repository (https://doi.org/10.6084/m9.figshare.17284517) (Shao and Luo, 2021).

Table 2 Environmental and ecological parameters.

^* E represents einsteins.

2.3 Statistical analyses

For Gamma A data points sampled in the same months and the same depth bins (defined in WOA), they were binned to 2^∘ × 2^∘ grids to help eliminate possible biases caused by concentrated samplings in specific regions, resulting in 939 binned means of log₁₀-based Gamma A nifH abundance. The corresponding environmental and ecological parameters were also averaged to the same bins when necessary.

We used a generalized additive model (GAM) using R package “mgcv” (Wood, 2017) to demonstrate nonlinear relationships between the multiple predictors and the Gamma A abundance:

$$\begin{array}{}\text{(3)}& y=\mathit{\alpha }+{\sum }_{i=\mathrm{1}}^{n}s\left(x_i\right)+\mathit{\epsilon },\end{array}$$

where y is the response variable (Gamma A abundance), x_i is the ith predictor (i.e., the environmental or ecological variable), α is the intercept, s(x_i) is a linear combination of smooth functions of predictor x_i, n is the number of predictors and ε is the standard error. To avoid over-fitting to noise, the restricted maximum likelihood (REML) method was selected for the GAM smoothing parameters of every predictor with the basis function number (k) set to 9 (Wood et al., 2016). In the model selection of GAM, a double-penalization approach was used to identify and remove those insignificant predictors with large smoothing parameters and set them to zero functions (Marra and Wood, 2011).

The scientific color maps are used in several figures to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies (Crameri et al., 2020).

3 Results and discussion

3.1 Global distribution of Gamma A nifH abundance

The nifH gene abundance ranges from 1 to 10⁷ copies L⁻¹ in the global ocean and shows an approximately log-normal distribution (Fig. S1). High abundance of Gamma A nifH abundance of over 10⁵ copies L⁻¹ is prevalent in the subpolar North Pacific, tropical Atlantic and Bay of Bengal (Indian Ocean) (Fig. 1). Most Gamma A abundance data were sampled above 100 m, particularly in the upper 25 m. The deepest datum with detectable Gamma A nifH was sampled at 885 m in the Southern California Bight (Hamersley et al., 2011). Available data showed that nifH abundance decreased with depth in the southwestern Pacific Ocean, the Indian Ocean and the South China Sea but did not have an apparent trend from the surface down to 200 m in the tropical Atlantic Ocean (Fig. S3). More data particularly in deep waters are needed to better and more reliably reveal the vertical pattern of Gamma A abundance.

Although high Gamma A abundance of over 10⁶ nifH copies L⁻¹ was observed in the surface North Pacific Ocean, the number of zero-value data was also massive (215 of a total of 608 data points) and even located close to those high-abundance data (Cheung et al., 2020) (Fig. 1), indicating the patchy distribution of Gamma A. As discussed already (Sect. 2.1), zero-abundance data were not included in the further analyses due to the patchiness of Gamma A and the limitations of the qPCR method in detecting the true absence of Gamma A.

Figure 1 Gamma A abundance (nifH copies L⁻¹). The panels show data in depth ranges of (a) 0–25 m, (b) 25–100 m and (c) 100–200 m and (d) below 200 m. For clear demonstration, data are binned to 2^∘ × 2^∘ and geometric means in each bin are shown. Zero-value data are denoted as black triangles.

3.2 Relationship between primary production and Gamma A abundance

If the presumption that Gamma A is heterotrophic or photoheterotrophic bacteria (Bombar et al., 2016; Zehr and Capone, 2020) is true, a positive relationship between the Gamma A abundance and net primary production (NPP) can be expected because its energetically intensive N₂ fixation can benefit from a sufficient supply of organic matter from primary producers. The significant positive correlation between the logarithm of Gamma A nifH abundance and the logarithm of NPP in our data (correlation of 0.21, p<0.01) (Fig. 2) was consistent with this presumption. However, this positive correlation could just reveal that Gamma A and primary producers share certain common controlling factors. For example, even if Gamma A were autotrophic or mixotrophic and were able to harvest energy from solar radiation, it could also positively correlate with NPP as both of them would be supported by high light intensity. Although the capability of Gamma A to fix N₂ has not been quantified, it could also be possible that the fixed N by Gamma A, if it occurred, could in turn support NPP. It is also interesting that the upper bound of the observed Gamma A abundance increased with the NPP [log₁₀(Gamma A) = 3.7 log₁₀NPP − 3.2] when NPP was less than 10^2.6 (≈ 400) mg C m⁻² d⁻¹, above which the upper bound of Gamma A abundance was saturated at ∼ 10⁷ nifH copies L⁻¹ (Fig. 2). This pattern indicated that the maximal potential of Gamma A abundance was positively associated with NPP, although it is unclear why the maximal Gamma A abundance ceased to increase when NPP increased above 400 mg C m⁻² d⁻¹. The observed reduction of Gamma A from its maximal potential could indicate that Gamma A was limited by other environmental factors and/or simply reflect the vertical distribution of Gamma A in the water column (see more analysis below), noting that the NPP used here is a depth-integrated rate of the euphotic zone.

Nevertheless, our results show that Gamma A tends to inhabit high-productivity waters. Our finding contradicts the hypothesis mentioned above that Gamma A prefers oligotrophic waters based on samples mainly in the tropical and subtropical Pacific and Atlantic oceans, in which Gamma A reached 8×10⁴ nifH copies L⁻¹ (Shiozaki et al., 2018a; Langlois et al., 2015). However, the new dataset (Cheung et al., 2020) included in the present study showed even higher (over 10⁵ nifH copies L⁻¹) Gamma A abundance in the subarctic North Pacific (Fig. 1), where nutrient concentrations and NPP are generally high.

Figure 2 The relationship between Gamma A abundance and net primary production. Both Gamma A abundance and net primary production (NPP) are log₁₀-transformed. The data with NPP of 10^2.0 to 10^2.6 mg C m⁻² d⁻¹ are divided into 12 groups with equal log-NPP intervals of 0.05, and the highest Gamma A abundance in each group is identified (red dots). The upper bound of Gamma A abundance (red line) is estimated by linearly fitting the red dots in the low NPP range and is saturated at ∼ 10^7.0 nifH copies L⁻¹ for NPP > 10^2.6 mg C m⁻² d⁻¹.

3.3 Multivariate nonlinear relationships between environmental factors and Gamma A abundance using the GAM

We then analyzed which environmental factors can be the predictors for Gamma A abundance. Multivariate GAM analysis was used to obtain a nonlinear relationship between Gamma A abundance and environmental parameters. Note that phosphate was not included in the GAM because its variance can be partly represented by P^∗.

The GAM results also showed that Gamma A abundance increased with NPP at low or intermediate levels of NPP (Fig. 3a), while the Fe concentration and PAR even contributed more than NPP to the variance in Gamma A abundance (Fig. 3b and c). The prediction of Gamma A abundance is also contributed by nitrate (only when its concentration is high) and silicate (Fig. 3d and e). P^∗ and temperature, however, did not show a clear pattern with Gamma A, particularly in the ranges where most data were located, although Gamma A abundance was reduced at temperature lower than 15 ^∘C (Fig. 3f and g). The GAM did not identify a significant relationship between Gamma A abundance and DOC concentration or MLD (Fig. 3h and i). In the following, we discuss the relationship of Gamma A abundance to the substantial predictors, including Fe, PAR, nitrate (with reference to P^∗) and silicate.

Figure 3 Partial effects of environmental variables on Gamma A abundance using GAM multivariate analysis. The analyses show the anomaly of Gamma A abundance contributed by the smooth function (red line) and its 95 % confidence interval (dashed blue line) of each environmental variable. Data (circles) are shown as partial residuals after all other partial effects have been considered. The numbers in the parentheses of y-axis labels are the degrees of freedom of the smooth functions. A degree of freedom smaller than 1 is equivalent to a linear line, and higher degrees of freedom represent more wiggly curves. The black ticks on the x axis and y axis indicate the location of the data. In panel (c), the unit E represents einsteins.

3.3.1 Iron

The Gamma A abundance generally showed a decreasing trend with increasing dissolved Fe (Fig. 3b). Our dataset showed that a high abundance of Gamma A was prevalently observed in the North Pacific Ocean (Fig. 1a), where Fe was considered the dominant limiting factor for N₂ fixation (Sohm et al., 2011). Other Gammaproteobacterial phylotypes, such as Gamma 3 and Gamma ETSP2, were also found to dominate the diazotrophic community in the eastern South Pacific (Turk-Kubo et al., 2014; Halm et al., 2012) where Fe heavily limited primary production (Knapp et al., 2016; Bonnet et al., 2008). It has also been suggested that Gamma A and unicellular cyanobacterial diazotroph UCYN-B share niches in the Fe-depleted western and southern Pacific Ocean (Moisander et al., 2014; M. M. Chen et al., 2019), possibly to avoid competing with other Fe-demanding diazotrophs. Gammaproteobacterial diazotrophs may be equipped with siderophore-releasing genes, such as those already reported in another versatile phylotype, Gamma 4 (Cheung et al., 2021), and the released siderophores are an efficient tool for scavenging low-level Fe in the ocean (Boyd and Ellwood, 2010). Although more studies are certainly needed to further explore the ecological and physiological mechanisms and their evolutionary reasons, the good survival of Gamma A in a low-Fe environment is an intriguing finding that may expand our recognized space of active N₂ fixation in the ocean.

3.3.2 Light

In the range that most PAR data spanned (≳ 1 E m⁻² d⁻¹; note that E represents einsteins), Gamma A abundance increased with PAR (Fig. 3c). This result represented a general pattern in which the Gamma A abundance decreased with depth (Fig. S3a, c and d) (Moisander et al., 2008; Langlois et al., 2015; T. Y. Chen et al., 2019; Shiozaki et al., 2014; Wu et al., 2019). This could be because Gamma A is photoheterotrophic, as hypothesized by previous studies (Moisander et al., 2014), but could also be because it took advantage of using photosynthetic products in the surface ocean (Langlois et al., 2015). The limited number of data showing the substantial presence of Gamma A under very low PAR (Fig. 3c), as well as the nearly constant Gamma A abundance with depth in the tropical Atlantic Ocean (Fig. S3b), raises the question of whether the active transport of organic matter from the surface can support the growth of Gamma A in the dark deeper ocean.

3.3.3 Nitrate and P^∗

Neither nitrate nor P^∗ generally had a substantial effect on Gamma A abundance, except that Gamma A abundance increased substantially when the nitrate concentration was higher than 10 µM (Fig. 3d). This is consistent with a previous review showing that nitrate did not show an immediate inhibition of Gamma A (Moisander et al., 2017) and the findings of abundant Gamma A in oceans with high nitrate concentrations (Bird and Wyman, 2013) or a shallow nitracline (Shiozaki et al., 2014). How and to what extent Gamma A, including the NCDs in general, is inhibited by nitrate remain unknown (Bombar et al., 2016). The hypothesis that low-nitrate and high-P^∗ environments favor autotrophic diazotrophs is based on the competition of inorganic nutrients between diazotrophs and other phytoplankton (Karl and Letelier, 2008; Deutsch et al., 2007), while our results tentatively suggest that competition may not occur strongly between NCDs and phytoplankton, although it is still unclear whether NCDs use inorganic or organic P sources. Nevertheless, our results do not suggest that high inorganic nutrients can be an inhibiting factor for Gamma A.

3.3.4 Silicate

Most Gamma A abundance data were associated with silicate concentrations higher than 10^0.5 (∼ 3.2) µM, and a positive relationship between them was revealed by the GAM (Fig. 3f), suggesting a possible association between Gamma A and diatoms. NCDs have been found on the surface of diatoms or on the diatom mats (Martínez et al., 1983), as discussed above. Diatom-dominant ecosystems tend to produce abundant large particles from either dead diatoms and their aggregates or the fecal pellets generated by zooplankton (Tréguer et al., 2018). Large particles can be a good habitat for NCDs, as already discussed. Our results then provide indirect evidence for the association between Gamma A and diatoms.

3.3.5 Predictions based on the GAM

Overall, the multivariate GAM model explained 43 % of the variance in observed Gamma A abundance (Fig. 4). The predicted Gamma A abundance generally followed the observed values, although it tended to underestimate the observed high Gamma A abundance (>10⁵ copies L⁻¹) or overestimate the low Gamma A abundance (<10² copies L⁻¹) (Fig. 4). The moderate explained variance indicated that although the tested environmental factors can substantially influence Gamma A abundance, there must be other untested factors and unknown mechanisms that can also substantially impact the Gamma A distribution.

Figure 4 Predictivity of the GAM. Comparison of predicted versus observed Gamma A nifH abundance. The red line is the 1 : 1 ratio of prediction to observation.

Although the overall R² was at a moderate level of 43 %, we applied this model to give a first-order estimate of Gamma A abundance in the surface ocean (Fig. 5a) from climatological NPP and environmental factors (Fig. S4), admitting that this demonstration did not fully cover the observed spatial variance in Gamma A abundance. The results suggest that Gamma A was most abundant in the upwelling region of the eastern tropical South Pacific (Fig. 5a). High Gamma A abundance was also predicted in certain regions of the Southern Ocean (Fig. 5a), which was mostly caused by its high nitrate concentration (Fig. S4g–h). However, the largest uncertainties for the predictions also exist in the Southern Ocean (Fig. 5b) as there were no Gamma A samples in this high-nitrate area (Fig. 1). Future sampling in the Southern Ocean can test our predictions and reduce the uncertainties. The lowest Gamma A abundance predicted was in the subtropical North Atlantic and in the southwestern Pacific near Australia (Fig. 5a), which was mostly caused by high Fe concentrations (Fig. S4c–d).

Figure 5 Prediction of Gamma A abundance. (a) Predicted annual mean surface (0–25 m) Gamma A abundance (nifH copies L⁻¹) and (b) the standard errors estimated by the GAM.

3.4 Impact of mesoscale eddies on Gamma A

The root-mean-square error (RMSE) of 0.84 and an R² of 43 % in the prediction model (Fig. 4) indicate that there was still substantial unexplained variance in Gamma A abundance. One possible reason for this was that we used the climatological monthly means for the environmental factors, but the in situ conditions can differ greatly from the climatological values. For example, oceanic mesoscale eddies can influence biogeochemical processes not only by advective transport but also by variations in the biological and chemical environments (McGillicuddy, 2016). Particularly, as discussed above, some regional studies have suggested that mesoscale eddies may influence the distribution of autotrophic diazotrophs and/or NCDs. We then explored whether the occurrence of mesoscale eddies can impact Gamma A abundance. We identified 59 data points of Gamma A abundance that were sampled in cyclonic eddies, while more (268) were sampled in anticyclonic eddies. This is consistent with the fact that eddies are more likely anticyclonic in the Northern Hemisphere, where most (74 %) of our sampling points were located (Chelton et al., 2011).

The results showed that the average Gamma A abundance within cyclonic eddies was substantially higher than that in anticyclonic eddies or outside eddies (t test, p<0.01) (Fig. 6a). We further separated the data into a low-NPP group and a high-NPP group at an NPP of 400 mg m⁻² d⁻¹, a value selected because our data showed that Gamma A abundance was saturated when NPP was higher than this level (Fig. 2). Interestingly, the phenomenon that Gamma A was more abundant in cyclonic eddies was mostly caused by the data in the high-NPP group (Fig. 6c) but not by the low-NPP group (Fig. 6b). To avoid the possible biases caused by different climatological conditions at the locations of cyclonic eddies, we then checked the residuals of the predicted Gamma A abundance using climatological factors (i.e., Fig. 4) and still found that the Gamma A abundance in cyclonic eddies in the high-NPP group was significantly higher (one-tailed t test, p<0.05) than the climatology-based predictions by on average nearly an order of magnitude, while this was not the case for samples in anticyclonic eddies and outside eddies in the high-NPP group or for all the samples in the low-NPP group (Fig. 6d–f).

Figure 6 Influence of mesoscale eddies on observed Gamma A abundance. (a–c) Gamma A abundance and (d–f) residuals of predicted Gamma A abundance using climatological NPP and environmental factors in Fig. 4, grouped according to the data sampled in cyclonic eddies (CE), anticyclonic eddies (AE) or outside eddies. The box plots show the median (central line), 25th and 75th percentile (upper and lower edges of the box), 5th and 95th percentile (error lines), and outliers (red crosses). The error bars within boxes show the mean value (purple dots) and its standard error. Values above black lines are p values of a two-tailed t test if the means of observed Gamma A abundance are equal (a–c) or a one-tailed t test if the residuals are greater than zero (d–f). The analyses are conducted for the whole dataset (a, d), for the data with local climatological NPP lower than 400 mg m⁻² d⁻¹ (b, e) and for the data with local climatological NPP higher than 400 mg m⁻² d⁻¹ (c, f).

These results indicated that cyclonic eddies could stimulate Gamma A abundance but only in the high-productivity oceans. This finding is contrary to a previous hypothesis on autotrophic diazotrophs that anticyclonic eddies form a nitrate-depleted and well-lit environment favorable to N₂ fixation (Davis and McGillicuddy, 2006; Fong et al., 2008; Church et al., 2009; Liu et al., 2020). However, our results suggested that Gamma A can benefit from the vertical pumping of nutrition-rich water driven by cyclonic eddies (McGillicuddy et al., 1998). Nevertheless, the biogeochemical consequences of mesoscale eddies can be complex (Gaube et al., 2014; McGillicuddy, 2016). For example, in addition to vertical pumping, the eddy stirring and trapping generated by mesoscale eddies can also have spatial effects on phytoplankton (Abraham, 1998; Wiebe and Joyce, 1992). Further sampling and studies are still needed to improve our mechanistic understanding of the effects of mesoscale eddies on both autotrophic cyanobacterial diazotrophs and NCDs.

3.5 Reliability of Gamma A nifH data

It is questionable whether the nifH copies measured using qPCR and collected in this study can reliably represent the abundance of Gamma A or even NCDs in general. When metadata are used, the reliability of comparison among absolute quantifications of nifH copies can be affected by methodological factors of qPCR assays. For example, even highly reproducible standard curves may result in significant variations in quantities of the same template in separated qPCR assays due to the log nature of the curve (Smith et al., 2006). The extraction method of nucleic acids, sample preparation, variations in the efficiencies of qPCR and differences in the qPCR platform can also impact the quantitative results (Smith and Osborn, 2009). In addition, the copy numbers of the nifH gene in Gamma A's genome remain unknown. There exists a large uncertainty regarding the extent to which nifH gene copies can represent Gamma A abundance, especially in contrast to its autotrophic counterparts. All these problems will need better technology to be resolved in the future.

4 Summary and outlook

With more measurements becoming available, we explored in this study what factors controlled the distribution of a representative phylotype of non-cyanobacterial diazotrophs, Gamma A, in the global ocean. The results of our study did not fully agree with the conclusion of a previous study that Gamma A preferred warm oligotrophic oceans (Langlois et al., 2015). Instead, our study suggests that Gamma A tends to inhabit high-productivity waters. In addition, our analyses also found that Gamma A was more abundant in Fe-depleted areas, possibly to avoid competition with autotrophic diazotrophs in high-Fe environments. Overall, our study suggests that productivity and Fe can be factors differentiating niches between non-cyanobacterial and cyanobacterial diazotrophs in the ocean, with the former favoring a high-productivity and low-Fe niche, while the latter occupies the opposite. Our study is also consistent with previous findings that Gamma A prefers well-lit surface oceans. The phenomenon of elevated Gamma A abundance in cyclonic eddies found in this study is worth further exploration for mechanistic understanding.

However, the moderate explanatory power of our prediction model indicates that there must be other unknown factors and mechanisms that also impact non-cyanobacterial diazotrophs. For instance, non-cyanobacterial diazotrophs found in the guts of copepods (Scavotto et al., 2015) imply that they are subject to top-down controls, which has also been suggested for marine autotrophic diazotrophs (Landolfi et al., 2021; Wang et al., 2019; Wang and Luo, 2022). Future studies should consider qPCR primers and probe sets targeting other NCDs such as Alphaproteobacteria and Cluster III phylotypes, which can also be important diazotrophs, particularly in previously unrecognized regions for marine N₂ fixation (Wu et al., 2019; Langlois et al., 2008; Martínez-Pérez et al., 2018; T. Y. Chen et al., 2019). The combination of PCR amplification and metagenomic data can identify a broader NCD community (Delmont et al., 2018) and may help us design a better universal primer that targets major NCDs. Last, the uneven spatial samplings of Gamma A, particularly the relatively scarce samples in the Southern Hemisphere, may also introduce biases into our analyses. More samples and studies are needed in the future to improve our understanding of the controlling factors, niches and distributions for non-cyanobacterial diazotrophs so that their contribution to global marine N₂ fixation can be better evaluated.