Quantifying the Impact of Climate Change on Marine Diazotrophy: Insights From Earth System Models

Nitrogen fixation is a major source of new nitrogen to the ocean, supporting biological productivity in the large nitrogen-limited tropical oceans. In Earth System Models, the response of nitrogen fixation to climate change acts in concert with projected changes to physical nitrogen supply to regulate the response of primary productivity in nitrogen-limited regions. We examine the response of diazotrophy from nine Earth System Models and find large variability in the magnitude and spatial pattern of nitrogen fixation in both contemporary periods and future projections. Although Earth System Models tend to agree that nitrogen fixation will decrease over the next century, strong regional variations exist, especially in the tropical Pacific which may counteract the response of the Atlantic and Indian oceans. As the climate driven trend of nitrogen fixation emerges by mid-century in the RCP8.5 scenario, on regional scales it may modulate the broad climate trends in productivity that emerge later in the century. The generally poor skill and lack of agreement amongst Earth System Models indicates that the climate response of nitrogen fixation is a key uncertainty in projections of future ocean primary production in the tropical oceans. Overall, we find that the future evolution of nitrogen fixation plays an important role in shaping future trends in net primary production in the tropics, but the poor skill of models highlights significant uncertainty, especially considering the role of multiple concurrent drivers.


INTRODUCTION
Future warming is projected to enhance ocean stratification, reducing the vertical supply of nutrient-rich subsurface waters to the nutrient-depleted surface waters of the euphotic zone in the tropical ocean in particular (Steinacher et al., 2010;Laufkötter et al., 2015). Two resources of particular importance in the tropics are phosphorus (P) and nitrogen (N), which are required for the production of DNA, ATP, proteins, pigments, and enzymes. A reduction in the vertical supply of N and P to the surface ocean due to ocean warming will cause phytoplankton to become nutrient limited, thereby restricting their growth and net primary production (NPP), with impacts on the biological carbon pump Moore C.M. et al., 2013;Moore et al., 2018). A common way to illustrate the impact of climate change on ocean properties is by using the concept of "time of emergence, " which assesses the time when the climate change signal in a given property emerges from the background variability. Such analyses illustrate that changes linked to NPP can emerge much later than those associated with sea-surface temperature (SST) and pH (Keller et al., 2014;Frölicher et al., 2016;Henson et al., 2017).
As dinitrogen (N 2 ) fixation is a globally significant N supply mechanism (68-164 Tg N yr −1 , Gruber and Sarmiento, 1997;Jickells et al., 2017;Luo et al., 2014;Luo et al., 2012;Tang et al., 2019a;Wang et al., 2019), future changes in N 2 fixation by diazotrophs may modulate any N limitation arising from enhanced stratification. This means that future changes to diazotrophy, and the supply of new fixed N, will be an important component of the response of tropical ocean productivity to climate change in Earth System Model projections. The process of marine N 2 fixation involves the reduction of N 2 gas to ammonia (NH 3 ) by marine organisms known as diazotrophs (Sohm et al., 2011). Fixing N 2 gas requires the iron-rich enzyme nitrogenase, which requires 5 to 100-fold more iron (Fe) than non-diazotrophs (Raven, 1988;Kustka et al., 2003). N 2 fixation also has an associated high energetic cost (Großkopf and Laroche, 2012), and like all phytoplankton, diazotrophs are susceptible to P limitation (Sañudo-Wilhelmy et al., 2001;Sohm et al., 2008). Overall, the need for light, Fe, and P, tends to confine diazotrophy to regions relatively rich in Fe and high light availability (such as the subtropical North Atlantic), and regions with ample P (such as the Western Pacific and subtropical South Atlantic, Falkowski, 1997;Mills et al., 2004;Goebel et al., 2008;Moore et al., 2009;Sohm et al., 2011). In N limited regions of the ocean, diazotrophs are only able to outcompete and proliferate provided that sufficient Fe and P are available (Moore et al., 2009;Dutkiewicz et al., 2012).
The most well studied diazotroph is the colony forming cyanobacterium, Trichodesmium, which occurs predominantly in the subtropical and tropical ocean (Sohm et al., 2011). Trichodesmium is one of a few diazotrophs that is able to be cultured in the laboratory allowing for detailed analysis of its physiology (Hutchins et al., 2013) and, it is one of the major contributors to global marine nitrogen fixation supplying ∼60-80 Tg N yr −1 (Mahaffey et al., 2005). For these reasons, Trichodesmium has been used as a model organism to parameterize diazotrophy in several global ocean models. These include the Earth System Models used for climate projections, which tend to restrict rates of diazotrophy based on light, Fe, P, and N availability. Recent laboratory work has demonstrated that diazotrophs are likely to be responding to a whole suite of interacting drivers such as temperature, CO 2 and nutrient co-limitation in the future (Hutchins et al., 2007(Hutchins et al., , 2013Walworth et al., 2016). Moreover, the contributions of unicellular diazotrophs such as Crocosphaera and the unicellular cyanobacterial group A (UCYN-A) which have high abundances and broad ecological ranges (Moisander et al., 2010), have not been incorporated into Earth System Models. UCYN-A for example has been shown to occupy a larger ecological range than Trichodesmium and may be a more important player in marine N 2 fixation away from the tropics (Martinez-Perez et al., 2016). It has recently been identified that heterotrophic diazotrophs are widespread and highly diverse (Bombar et al., 2016;Moisander et al., 2017), with the ability to fix N 2 in environments that are traditionally believed to exclude diazotrophy (Rahav et al., 2013;Benavides et al., 2016). Heterotrophic diazotrophs are also not currently incorporated into Earth System Models, which may indicate that these models are underestimating the areal extent of N 2 fixation by neglecting the niche occupied by heterotrophic diazotrophs.
As part of the IPCC 5th assessment report (IPCC-AR5), 9 of the Earth System Models used to project future changes in ocean biogeochemistry included marine N 2 fixation (Ciais et al., 2014). These models represent diazotrophy with a range of complexity (Supplementary Table 1), at the simplest only including the effect of temperature and N inhibition, while others incorporate singular limitation by Fe and P. More recently, climate projections have also been made using models with variable non-diazotroph cellular stoichiometry (IPSL-Quota, Kwiatkowski et al., 2018). In this study, we compare 8 IPCC-AR5 models and the more recent variable cellular stoichiometry model IPSL-Quota to examine how ocean diazotrophy is projected to respond to climate. While previous studies have conducted intermodel comparisons for N 2 fixation, they used arbitrary scaling factors for some models (Riche and Christian, 2018;Tang et al., 2019a) and did not assess the timescales and regionality of change. In this study, we update the reported global trends using the correct units for all models, address the regional distribution of N 2 fixation and explore how the time of emergence of N 2 fixation varies on global and regional scales both within and across models.

METHODS
Model outputs of depth integrated N 2 fixation from the coupled model intercomparison project 5 (CMIP5) that compiled model outputs were obtained from 8 models, alongside a more recent quota-based formulation of the IPSL-CM5A-LR model from Kwiatkowski et al. (2018). For each model, data from the historical (1850-2005) and future (RCP 8.5, 2006(RCP 8.5, -2100 scenarios was re-gridded to fit a standard 1 • × 1 • global grid. To calculate differences in N 2 fixation rates, both spatially and globally integrated, we used 20-year averaged time slices from the historical (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) and RCP 8.5 (2081-2100) simulations for all 9 models. We further investigated regional patterns of depth integrated N 2 fixation rates using an established biogeochemical province approach (Longhurst, 2007) and derived the average trend of depth integrated N 2 fixation for each province. We then compared the percentage of models that agreed with the average trend to interpret model agreement. We note here that prior work used an arbitrary scaling of 0.17 for the two IPSL model in the CMIP5 archive (Riche and Christian, 2018;Tang et al., 2019a). This is due to IPSL output being stored in the CMIP5 archive in carbon, rather than nitrogen, units. We here report results from those models using the correct adjustment of 1/7.625 (or 0.131), which reflects the model N/C ratio (Aumont and Bopp, 2006). Models vary markedly in the way in which they represent N 2 fixation. Some represent an explicit diazotroph functional group, which obtains fixed N from N 2 fixation. Others include N 2 fixation implicitly as a source of NH 3 under certain conditions. While others simply apply a N flux that restores the N:P ratio back to Redfield values of 16:1 (Supplementary Table 1).
Time of emergence (ToE) analysis was performed to detect when the climate driven N 2 fixation trend emerges from the natural variability within the models (e.g., Keller et al., 2014). The ToE rests on the concept of comparing the temporal climate change signal to the background noise, which is determined using different thresholds. Similar to Frölicher et al. (2016), we assume that the N 2 fixation signal had emerged when it was greater or less than the mean of the preindustrial control simulation ± 1 standard deviation consistently for 4 years in the direction of the overall trend (to avoid the false diagnosis of shorter oscillations). The sensitivity of the time of emergence to our noise threshold was assessed by also using 2 standard deviations from the preindustrial control simulation. The ToE at each model grid cell was calculated by comparing the annual mean global N 2 fixation time series from the RCP8.5 scenario  to the parallel preindustrial control simulations which had been run for at least 300 years. We averaged the ToE across each biogeochemical province to facilitate comparisons.

Contemporary Nitrogen Fixation: Spatial Patterns and Global Trends
As in prior work (Riche and Christian, 2018;Tang et al., 2019a), there are large variations between models regarding the magnitude and spatial pattern of diazotrophy in the contemporary ocean (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005). Across the range of models, the updated magnitude of total global N 2 fixation using the correct units for all models, ranged 2.5-fold from 76.5 to 185.6 Tg N yr −1 . Despite this range, the models broadly fell into two groups based upon their spatial patterns. The first group of models placed emphasis on higher rates of N 2 fixation in the upwelling regions of the ocean, while the second group of models placed the majority of N 2 fixation in the low nutrient subtropical gyres (Figure 1). Some models (CanESM2, CESM1-BEC, IPSL-CM5A-LR, IPSL-CM5A-MR, and IPSL-Quota) imposed a minimum temperature restriction, which restricts diazotrophs to low latitudes, but for those models without thermal limits, diazotrophy can occur over the entire global ocean (albeit at very low rates at high latitudes, Figure 1) this may be an imperfect solution to account for diazotrophy in the high latitudes. These models without thermal limits on diazotrophy tended to fall into the upwelling group which on the whole had larger global N 2 fixation values than gyre models. One exception to this general rule is the N only model CanESM2, which does not include the role of secondary limiting factors, such as Fe and P which exert more control on diazotrophy than N. The large disagreement between the models suggests that our understanding of spatial distribution of N 2 fixation is limited. Interestingly, the recent re-evaluation of global ocean N 2 fixation tends to support the concept of the greatest rates of N 2 fixation occurring in low nutrient subtropical gyres rather than in upwelling regions (Wang et al., 2019).
The difference in the spatial pattern between gyre models and upwelling models may be reflected in differences in the rate of  total net primary production. For instance, when the ratio of N 2 fixation to net primary productivity is plotted across the models we find stronger spatial agreement than for N 2 fixation alone (Supplementary Figure 1). N 2 fixation is a small fraction of net primary production in the upwelling regions for all models, increasing to a larger fraction in the gyres. This suggests that the higher rates of N 2 fixation associated with upwelling regions in the GFDL and MPI models is a result of a much larger rate of upwelling zone NPP in these models, which is allied to a stronger contrast in net primary production between the upwelling and gyre regions.

Future Changes to Marine Diazotrophy: Spatial Patterns and Global Trends
Across the models evaluated, 8 predict a decrease or very little change in global N 2 fixation between 2006 and 2100 (between +1 and −31%), with only one model (CanESM2) projecting an increase in diazotrophy. This is because CanESM2 has no limitation by Fe and P, and in this model N 2 fixation responds positively to future ocean warming (Figure 2A). The N 2 fixation trend is not consistent between models when assessed over the Atlantic, Pacific and Indian ocean basins, such that the global trend can mask important inter-basin variability. For example, stronger trends were predicted for the Atlantic and Indian ocean basins and were relatively consistent with the direction of the global trend (Figures 2B,D). However, the Pacific Ocean exhibited the most variability, with the sign of the trend being opposite to the global trend for the majority of models ( Figure 2C). This suggests that diazotrophy in the Atlantic and Indian oceans will be negatively affected by climate change, but the Pacific Ocean may counteract the decrease occurring in the other basins in some cases. It is, however, difficult to assess how realistic these trends are as currently there are very few observations of N 2 fixation in the Indian Ocean basin (Luo et al., 2012;Tang et al., 2019b).
The changes in the total global N 2 fixation rate also exhibited large variability once updated from prior work to use the correct units for all models, ranging from −50.1 Tg yr −1 to +58.0 Tg yr −1 ( Table 1). When models are grouped based on their contemporary pattern of diazotrophy, upwelling models had very little coherency for the spatial pattern of diazotrophy (comparing the 2 MPI and 2 GFDL models, Supplementary  Figure 2). For gyre models, however, moderate coherence existed between the 3 IPSL models and the CESM1-BEC model (CanESM2 was overlooked hereafter as an outlier). This shows that although upwelling models display a similar spatial pattern for contemporary diazotrophy, there is a little agreement on how the spatial pattern of diazotrophy will change, suggesting that mechanisms controlling diazotrophy in these models may actually be very different. For example, the 2 GFDL models  1986:2005)], time of emergence of global N 2 fixation when the natural variability is calculated using 1 standard deviation (σ) and 2 standard deviations (2σ) of the control simulation, model skill represents the correlation between model depth integrated N 2 fixation (1986:2005) and observed depth integrated N 2 fixation (Tang et al., 2019b).

Model
Total N 2 fixation anomaly (Tg N yr −1 ) impose Fe limitation on diazotrophy but the MPI models do not. On the other hand, the relatively good agreement between the gyre based models matches their common underlying drivers (Supplementary Table 1).

Time of Emergence
Taking a province-based view to synthesize the temporal and spatial model comparison revealed regions where the average N 2 fixation trend for the 9 models was projected to be strongly modified. The strongest decreases in diazotrophy occurred in the equatorial Atlantic, western equatorial Pacific and northern Indian. Provinces with the strongest decreases in diazotrophy coincided with provinces with the highest intermodel agreement (Figures 3A,B). The strongest overall increase in diazotrophy was found in the Pacific upwelling region, however, this province coincided with relatively low intermodel agreement (Figures 3A,B). In general, the provinces with lowest agreement corresponded to regions of low average trend, suggesting that these lower average trends were a consequence of divergence between models.

Future Changes to Marine Diazotrophy: Time of Emergence
The N 2 fixation trend due to climate change emerges from the natural variability across the models between 2018 and 2040 when the natural variability is calculated using one standard deviation of the control simulation (Table 1). If we allow greater natural variability (by using 2 standard deviations), then the emergence is delayed by 5-20 years ( Table 1 and Supplementary Figure 3). There is no clear link between the magnitude of the projected change in future N 2 fixation and the time of emergence, however, when considering both sensitivities (1 and 2 standard deviations) the gyre based models have better inter model agreement of emergence (2018-2042) than the upwelling based models, which have a wider range of emergence times (2020-2054). Regionally, the N 2 fixation trend emerges earliest (∼2030-2050) in the lower latitudes and much later (∼2080->2100) in the high latitude Southern Hemisphere (Figure 3C), suggesting that N 2 fixation in the low latitudes is most susceptible to change. The high latitude emergence of the N 2 fixation signal may be difficult to assess as not all models allow high latitude N 2 fixation in their mean state. On the whole N 2 fixation decreases (e.g., in the tropical Atlantic) or shows very little change (e.g., in the subtropical gyres) in the tropical ocean, with a marked exception being the equatorial Pacific, which shows a significant increase in N 2 fixation. This is likely due to climate driven changes in upwelling, which reduces N supply therein. Any change to the magnitude of "new" N being supplied by N 2 fixation in the tropics will serve to amplify or dampen any stratification driven reduction in N supply. N 2 fixation trends emerge later than trends in both sea surface temperature and pH (2016)(2017)(2018)(2019)(2020)(2021)(2022)(2023)(2024)(2025)(2026)(2027)(2028)(2029)(2030)(2031)(2032)(2033)(2034)(2035), but earlier than trends in both oxygen and NPP (>2095, Frölicher et al., 2016;Henson et al., 2017;Keller et al., 2014;Rodgers et al., 2015).
Although climate change will directly impact NPP, it is important to also investigate the processes that indirectly affect NPP, such as N 2 fixation and how it may modulate the NPP signal. For instance, in regions where N 2 fixation declines (western and subtropical N. Pacific, Indian Ocean and the subtropical N. Atlantic) there is a strong contribution (between 4 and 73%) from the decrease in N 2 fixation to the overall decrease of the upper ocean nitrate (NO 3 ) inventory (Figures 3A,D). On the other hand, areas of increased N 2 fixation act against the loss of NO 3 from the upper ocean due to enhanced stratification (Figures 3A,D). In Figure 3D, a positive ratio indicates that the N 2 fixation trend is in the same direction as the overall upper NO 3 inventory trend, whereas a negative ratio indicates opposing trends ( Figure 3A shows the overall directional trend in N 2 fixation). As NPP in the topics is largely responding to changes to upper ocean N, the impact of climate driven changes in stratification may be amplified or dampened by the more rapidly emerging trends in rates of N 2 fixation. Future work should also focus on evaluating the processes underpinning this potentially important feedback (e.g., rates of nitrification, zooplankton recycling and the sinking and remineralisation of diazotroph organic material) within the nitrogen cycle in Earth System Models.

Mechanistic Limitations of Current Models
Our analysis has shown that decadal trends in N 2 fixation in response to climate change are likely to contribute to the changes in NPP simulated in Earth System Models. But we have found that a range of complexity exists for how N 2 fixation is parameterized within the current models and large disagreement occurs in the spatial pattern and trends of N 2 fixation projected by the different models. However, despite these differences, all models show poor skill compared to compiled N 2 fixation observations ( Table 1). We took the observations from the nitrogen fixation database of Tang et al. (2019b), regridded them onto the model horizontal grid and compared them to FIGURE 3 | Regional comparison of the 9 global nitrogen fixation models using Longhurst provinces: (A) average trend of nitrogen fixation across all 9 models in each province, (B) percentage of models that agree with the average trend for each province, (C) average time of emergence across all 9 models for each province using 1 standard deviation of the control and (D) the contribution of the nitrogen fixation anomaly to the NO 3 anomaly. the model fields averaged over 1985:2005. Overall, model skill, taken as the correlation coefficient, ranges from as low as 0.06 for GFDL-ESM2M or up to 0.35 for MPI-ESM-LR (Table 1). This suggests that the change in N 2 fixation in response to climate change and the associated impact on nitrogen supply to primary producers in the oligotrophic ocean remains a very poorly constrained component of Earth System Models. The overall trend in N 2 fixation across the majority of the ESMs is at odds with other perturbation experiments predicting the opposite (Boyd and Doney, 2002;Dutkiewicz et al., 2014;Jiang et al., 2018;Tang et al., 2019a), highlighting the need for more consideration of the competition between diazotrophs and non-diazotrophs within ESMs in the future.
There are a range of potentially significant drivers not yet included in the parameterization of ocean N 2 fixation in Earth System Models, which could affect the response to climate change and contribute to the general lack of model skill. Firstly, increasing CO 2 has been shown to produce a "nutrient-like" response for growth and N 2 fixation rates for both Trichodesmium and Crocosphaera spp. (Hutchins et al., 2013;Hutchins and Fu, 2017), with CO 2 potentially defining an upper limit for maximum N 2 fixation rates (Hutchins et al., 2013). However, increasing CO 2 will increase ocean acidification, which has been shown to negatively impact diazotrophy by reducing the efficiency of the nitrogenase enzyme (Hong et al., 2017). Secondly, warming may push tropical diazotrophs past their thermal optima and/or facilitate the expansion of recently observed cold water diazotrophs occurring below 20 • C (Needoba et al., 2007) and even in polar waters (Harding et al., 2018;Shiozaki et al., 2018). Breitbarth et al. (2007) found that the maximum rate of growth and N 2 fixation in Trichodesmium occurs between 24 and 30 • C and both Trichodesmium and Crocosphaera growth has been shown to rapidly decline at temperatures above 28 and 30 • C, respectively . Most models do not include an upper thermal threshold and so depending on how growth is parameterized, rates of diazotrophy could continue to increase past observed thermal thresholds. If the ecological niche for diazotrophy moves poleward due to excessive warming, light limitation may become a more important factor than previously (Boyd and Doney, 2002;Boatman et al., 2017). Thirdly, laboratory studies show that when Fe and P are both present at low concentrations, Trichodesmium cell sizes decreased, allowing uptake rates to be maintained and promoting both growth and diazotrophy, suggesting that Fe-P co-limitation may be as important in shaping N 2 fixation rates (Walworth et al., 2016(Walworth et al., , 2018. However, dust deposition is predicted to increase due to enhanced desertification (Mahowald et al., 2017). This may promote diazotrophy by supplying more Fe, but would enhance P limitation (Mills et al., 2004) and negatively impact the effect of Fe:P co-limitation.
Previously it was believed that the contribution of coastal regions to global marine N 2 fixation was very small, but over the last decade the importance of coastal diazotrophy has been highlighted, with it being estimated to contribute 16.7 ± 14.3 Tg N yr −1 (Mulholland et al., 2012;Tang et al., 2019b). Due to the coarse resolution of climate models, coastal environments are poorly represented and will need to be improved in the future to gain a more complete understanding of global marine N 2 fixation. The majority of models have an explicit parameterization of diazotrophy use only one diazotroph for the entire ocean, but actual diazotroph biodiversity is far more complex (Tang et al., 2019b). Similarly, the reliance of models on the colonial cyanobacterium Trichodesmium as a model organism neglects other important groups such as the UCYN-A, which may be predominant globally (Martinez-Perez et al., 2016), and heterotrophic diazotrophs which occupy niches previously believed to be unfavorable to diazotrophs (Rahav et al., 2013;Benavides et al., 2016). Model studies tend to focus upon the impact of bottom-up control on diazotrophs and generally overlook top-down controls. For instance, Wang et al. (2019) showed that when the zooplankton maximum grazing rate on diazotrophs was decreased to that of diatoms, global N 2 fixation increased by 62%. Diazotroph growth rates are also not well constrained, and models are hence biased toward relying on observed growth rates of Trichodesmium, which will restrict the modeled diazotroph to their niche (Dutkiewicz et al., 2009). In the future, a new generation of N 2 fixation models that include these multiple interacting effects and account for diazotroph biodiversity would help increase confidence in our projections of how climate change will affect the ocean N cycle (Tang and Cassar, 2019) and net primary production in the tropics.

CONCLUSION
Nitrogen fixation will be crucial to supporting productivity in the future tropical ocean that is dustier and more stratified. Current models project different spatial patterns and magnitudes of diazotrophy in the contemporary ocean, models, however, broadly fall into two groups based upon their spatial distribution of N 2 fixation. Generally models project that globally diazotrophy will decrease in the future, at the regional scale the models reveal more complexity in the response to climate change with the Pacific counteracting the decline projected for the Atlantic and Indian Oceans. Gyre models exhibited more coherency in the response of N 2 fixation at the regional scale than upwelling models. For the models as a whole, the strongest intermodel agreement coincided with regions of strong decreases in diazotrophy in the Atlantic and Indian Oceans, whilst low model agreement occurred with strong increases in diazotrophy in the Pacific. The change in new nitrogen supply from diazotrophy may alleviate any changes to N limitation in the oligotrophic ocean due to increased stratification. The rapid emergence of the climate driven trend of diazotrophy from background variability by around 2040, implies that changes in new nitrogen supply from diazotrophy is likely to contribute to the modulation of the later emerging trends in net primary productivity on a regional scale. Overall, there is a need for a new generation of N 2 fixation Earth System Models that account for the multiple interacting factors which drive diazotrophy including both topdown and bottom-up control. Understanding how climate will impact diazotrophy is essential to gain insight into how patterns of nutrient limitation and net primary productivity are likely to respond to future climate perturbations.

DATA AVAILABILITY STATEMENT
The model data that was analyzed in this study was generated as part of CMIP5 and is freely available from the Earth System Grid Federation (https://esgf-node.llnl.gov/projects/cmip5/).

AUTHOR CONTRIBUTIONS
LW and AT designed the study and wrote the manuscript. LW performed the data analysis. Both authors contributed to the article and approved the submitted version.

FUNDING
This work contained in this publication was conducted during a Ph.D. study supported by the Natural Environment Research Council (NERC) EAO Doctoral Training Partnership and is fully-funded by NERC (NE/L002469/1).