Projected expansion of Trichodesmium’s geographical distribution and increase in growth potential in response to climate change

Estimates of marine N2 fixation range from 52 to 73 Tg N/year, of which we calculate up to 84% is from Trichodesmium based on previous measurements of nifH gene abundance and our new model of Trichodesmium growth. Here, we assess the likely effects of four major climate change‐related abiotic factors on the spatiotemporal distribution and growth potential of Trichodesmium for the last glacial maximum (LGM), the present (2006–2015) and the end of this century (2100) by mapping our model of Trichodesmium growth onto inferred global surface ocean fields of pCO2, temperature, light and Fe. We conclude that growth rate was severely limited by low pCO2 at the LGM, that current pCO2 levels do not significantly limit Trichodesmium growth and thus, the potential for enhanced growth from future increases in CO2 is small. We also found that the area of the ocean where sea surface temperatures (SST) are within Trichodesmium's thermal niche increased by 32% from the LGM to present, but further increases in SST due to continued global warming will reduce this area by 9%. However, the range reduction at the equator is likely to be offset by enhanced growth associated with expansion of regions with optimal or near optimal Fe and light availability. Between now and 2100, the ocean area of optimal SST and irradiance is projected to increase by 7%, and the ocean area of optimal SST, irradiance and iron is projected to increase by 173%. Given the major contribution of this keystone species to annual N2 fixation and thus pelagic ecology, biogeochemistry and CO2 sequestration, the projected increase in the geographical range for optimal growth could provide a negative feedback to increasing atmospheric CO2 concentrations.


| INTRODUC TI ON
Marine phytoplankton account for about 45% of global net primary production (Field, Behrenfeld, Randerson, & Falkowski, 1998), and as such play an important role in the global carbon cycle (Arrigo, 2007).
Approximately 20% of the annual marine net primary production is exported from the surface to the deep via sinking particles (Buesseler & Boyd, 2009). This export production contributes to the draw-down of CO 2 from the atmosphere and its sequestration for hundreds or thousands of years in the deep ocean. Maintaining net primary production requires a source of fixed N (e.g. ammonium, nitrite, nitrate and organic N) which can be supplied to the euphotic region by mixing and upwelling of nitrate from the deep, deposition of nitrate and organic N from the atmosphere and N 2 fixation by diazotrophic cyanobacteria (Duce et al., 2008).
Nitrogen fixation accounts for more than half of the export of organic carbon from the surface ocean to the deep ocean in some parts of the oligotrophic tropical and subtropical oceans (Capone, 2005), and is likely to be two to three times more important than the atmospheric delivery of fixed N to the sea (Duce et al., 2008). Changes in export production could significantly affect the ocean's ability to sequester CO 2 from the atmosphere and store it in the deep ocean.
When operating over geological timescales, for example, between glacial and interglacial periods, even small changes in the balance between N 2 fixation and the loss of fixed N due to denitrification can significantly affect the amount of CO 2 that can be stored in the ocean (Falkowski & Raven, 1997;Kohfeld & Ridgwell, 2009).
One manifestation of global warming is the increase in sea surface temperature (SST), which enhances water stratification and leads to the expansion of hot tropical and warm subtropical regions (Doney, Yeager, Danabasoglu, Large, & Mcwilliams, 2007). Although the reduced flux of nitrate into the upper mixed layer associated with increased stratification will be detrimental to many phytoplankton groups, the expansion of these nitrogen-limited regions may give diazotrophic cyanobacteria like Trichodesmium a competitive advantage and accentuate competition for other limiting nutrients (e.g. Fe or P). Alongside global warming, increasing atmospheric CO 2 is driving increases in seawater CO 2 concentrations, lowering pH and changing the inorganic carbon (Rost, Zondervan, & Wolf-Gladrow, 2008) and iron chemistry (Millero, 2009;Shi, Sun, & Falkowski, 2007;Shi, Xu, Hopkinson, & Morel, 2010), all of which will have consequences for the structure and functioning of marine ecosystems.
In contrast, Monteiro, Dutkiewicz, and Follows (2011) concluded that Trichodesmium's niche will be constrained at higher latitudes due to the loss of oligotrophic conditions and competition for nutrient resources. The generation and transport of Fe-containing dust into the ocean (Mahowald & Luo, 2003;Tegen, Werner, Harrison, & Kohfeld, 2004) will also affect the growth and productivity of Trichodesmium and other diazotrophs (Moore et al., 2006). In addition, Jiang et al. (2018) have shown a temperature dependence of Trichodesmium's ability to assimilate Fe that in turn affects growth and nitrogen-fixation under Fe-limiting conditions.
To gain insight into changes in the distribution and growth potential of Trichodesmium from the last glacial maximum (LGM) to the present (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) and from the present to the end of this century, we combined information on how the growth of Trichodesmium responds to four key abiotic factors (temperature, CO 2 , irradiance and iron availability) with inferred global surface ocean maps of these variables for these time periods.

| Growth rate model
The dependence of the steady state balanced growth rate of Trichodesmium on temperature, irradiance and iron was modelled as a multiplicative function: where ′ max is the hypothetical maximum growth rate (day −1 ) at the optimum temperature for growth before taking into account photoinhibition; T is the temperature (°C); T min and T max are the minimum and maximum temperature limits for growth (°C); θ is a shape-determining parameter which alters the skewness of the µ-T relationship; Φ is a shape-determining parameter which alters the kurtosis of the µ-T relationship; Feʹ is the sum of inorganic iron complexes (Iron hydroxides + Fe(II)) (pM); K m is the half saturation concentration for Fe (pM); E is the irradiance (mol photons m −2 day −1 ); E kα is the light saturation parameter (mol photons m −2 day −1 ); and E kβ is the photoinhibition parameter (mol photons m −2 day −1 ).
The parameterization of our Trichodesmium growth rate model was based on a series of long-term laboratory experiments Boatman, Oxborough, et al., 2018). In total, 184 treatments were cultured semi-continuously under well-defined growth conditions. Temperature (19°C-32°C), irradiance (10-1,400 µmol photons m −2 s −1 ) and iron growth response curves were collected in parallel using identical methods and equipment and the data were analysed using the same script-based code. All experiments consisted of a low (180 ppm), mid (380 ppm) and high (720 ppm) CO 2 treatment, as well as a low (40 µmol photons m −2 s −1 ) and high (400 µmol photons m −2 s −1 ) irradiance treatment (Figure 1).
A more detailed description of the laboratory methodology, culturing technique and analytical procedure is reported in our previous studies Boatman, Oxborough, et al., 2018).
As highlighted by Low-Décarie et al. (2017), each equation was objectively selected based on the shape and data resolution of the growth response curves. The temperature response was modelled (1) using a sine function , the light response using a two-phase exponential (initial slope and photoinhibition) function (Platt & Gallegos, 1980) and the Fe-response using a rate saturating function (Michaelis & Menten, 1913).

| Parameter optimization
Curve fitting to our multifactorial dataset allowed growth rate to be modelled by specifying eight parameter values (Table S1), seven of which were independent of the CO 2 , with one ( ′ max ) being CO 2dependent. Curve fitting was performed on the median growth rate using a weighted nonlinear least squares algorithm, where weights were the reciprocals of the standard errors associated with the growth rates. Initially, each CO 2 dataset (i.e. 180, 380 and 720 ppm) was modelled independently; producing three sets of parameterizations; each set consisting of eight parameters (Table S2). Then starting with T min , data were combined between CO 2 datasets, and the model re-optimized. The data groupings consisted of: (a) combined low/mid CO 2 , (b) combined low/high CO 2 , (c) combined mid/ high CO 2 and (d) all CO 2 data combined. Note, as previous studies have shown that Trichodesmium's maximum growth rate (µ max ) is significantly reduced at low CO 2 relative to mid CO 2 , and is consistently ~10% higher from mid to high CO 2 (Boatman, Davey, Lawson, & Geider, 2018;Boatman et al., 2017;Boatman, Mangan, Lawson, & Geider, 2018;Boatman, Oxborough, et al., 2018), ′ max was maintained independent between CO 2 conditions. Having processed all model iterations by optimizing all parameter groupings, maximum likelihood (L) values were calculated using each model's residual sum of squares and the number of data points (n = 184). The best model (Model 29 in Table S2) consisted of seven global constants (T min , T max , E kα , E kβ , θ, Φ and K m ) with one CO 2 -dependent parameter ( ′ max ; Table S1). Modelled against observed growth rates gave an r 2 value of .934 ( Figure S1).

| Modelled oceanographic data
Mean monthly oceanographic data for sea surface chlorophyll (kg/m 3 ), temperature (°K), net downward radiation (W/m 2 ), dissolved iron (dFe; mol/m 3 ) and mixed layer depth (MLD; m) were obtained from the Institute Pierre Simon Laplace (IPSL), using the IPSL-CM5A-LR earth system model (Dufresne et al., 2013). The data resolution was 1.875° × 3.75° with 39 vertical levels for the atmosphere and about 2° (with a meridional increased resolution of 0.5° near the equator) with 31 vertical levels for the ocean. All data files were transformed from a tripolar grid to a 180° × 360° latitudinal/longitudinal rectangular grid format, and interpolated to a 1° resolution.
Data files for the aforementioned variables were collected using the lgm and rcp60 experiment, where the lgm experiment yielded data for the LGM while the rcp6.0 experiment yielded data for the present and future timescales. Mean monthly data files for each variable, at all three time periods, were averaged over a decade; present data were averaged between 2006 and 2015 and future data were averaged between 2091 and 2100. Based on the past, present and future climate emissions within the CM5A-LR lgm and rcp6.0 model, these time periods best match the low, mid and high CO 2 growth conditions of the laboratory experiments.

F I G U R E 1
The response of Trichodesmium erythraeum IMS101 growth rates to iron concentration (a, d), temperature (b, e) and irradiance (c). Growth rates (day −1 ) were calculated as a multiplicative function of temperature, iron and irradiance using Equation (1). Note, temperature and iron responses were measured at low (LL = 40 μmol photons m −2 s −1 ) and high (HL = 400 μmol photons m −2 s −1 ) light. Shipboard observations of Trichodesmium biomass (Luo et al., 2012) are plotted against the SST (f); where the dashed lines represent the minimum (T min ) and maximum (T max ) temperature limits for growth, and the red line is the temperature growth response curve under present conditions modelled using Equation (1) with the parameter values reported in Table S1 2.4 | Irradiance within the mixed layer Net downward radiation (W/m 2 ) was converted into a photosynthetic photon flux density (PPFD) by assuming half of the irradiance is photosynthetically active radiation (PAR) and that 1 W/m 2 of PAR is equivalent to a PPFD of 4.57 µmol photons m −2 s −1 (Langhans & Tibbitts, 1997). The mean PPFD within the mixed layer was calculated as follows (Helbling, Villafañe, & Holm-Hansen, 1994): where E z is the mean PPFD within the mixed layer (µmol photons m −2 s −1 ), E 0 is the surface PPFD (µmol photons m −2 s −1 ), z is the MLD (m) and K d (PAR) is the light attenuation coefficient (m −1 ) which was calculated using the Chl a data as follows (Dennison et al., 1993): The mean PPFD within the mixed layer was then converted into units of mol photons m −2 day −1 by calculating the photic period for each 1° longitudinal-latitudinal cell as follows: where D is the daylength (hr −1 ), p is the daylength coefficient (°), L is latitude (°), Φ is the Sun's declination angle, J is day of the year and θ is the revolution angle. This calculation is time accurate to within 1 min for latitudes between 40°N/S, increasing to 7 min up to latitudes between 60°N/S (Forsythe, Rykiel, Stahl, Wu, & Schoolfield, 1995).

| Iron within the mixed layer
Based on measurements from Rijkenberg et al. (2008), Fe(III) speciation in the tropical North Atlantic Ocean produces free iron (Feʹ) concentrations below what is required to support Trichodesmium growth.
In contrast, Fe(II) concentrations are several orders of magnitude greater than Fe(III)ʹ, and are therefore an important source of bioavailable iron for Trichodesmium. Previous studies report that Fe(II) accounts for 20% of surface dissolved Fe (dFe) concentrations in the Baltic (Breitbarth et al., 2009), 12%-14% in the Pacific (Hansard, Landing, Measures, & Voelker, 2009) and 5%-65% in the South Atlantic and Southern Ocean (Bowie, Achterberg, Sedwick, Ussher, & Worsfold, 2002;Sarthou et al., 2011). Based on these observations, we assumed a global estimate where 25% of dFe is present as Fe(II).

| Growth rate limitation maps
For each month, at each timescale, the degree to which temperature, irradiance and iron limits Trichodesmium IMS101 growth was calculated as follows: where T (°C), E (mol photons m −2 day −1 ) and Feʹ (pM) correspond to the spatial resolved oceanographic data values, and ′ max , T min , T max , E kα , E kβ and K m are the model parameterization values as reported in Table S1.
All limitation plots exhibit a scale ranging from 0 (no growth potential) to 1 (no limitation to maximal growth; Figures S2-S4).

| Global primary production
To estimate Trichodesmium's contribution to current ocean production, we applied our growth rate model (Equation 1) to a map of Trichodesmium biomass; inferred from a map of nifH gene abundance (Tang & Cassar, 2019). The modelled nifH gene abundance was converted into Trichodesmium biomass (mg C/m 2 ) by multiplying by the same assumed 0.3 conversion factor (mg C/10 6 nifH copies) as reported in Luo et al. (2012).
Global growth rates were re-calculated using the same mean monthly oceanographic data used by Tang and Cassar (2019) to generate their map of nifH gene abundance, which included SST (°C) from SeaWIFS, PAR in the MLD (mol photons m −2 day −1 ) from MODIS and dissolved iron (dFe; nM) from Community Earth System Model. Note, we used the same global assumption that 25% of dFe is present as Fe(II). Data from the supporting information of Tang and Cassar (2019) were downloaded (https://doi.org/10.1594/PANGA EA.905108) and interpolated to a 1° resolution. Trichodesmium's mean annual primary production (g C m −2 year −1 ) was calculated as follows: where B(t) is the mean biomass (g C/m 2 ) in each of the 12 months, µ(t) is the mean growth rate (day −1 ) and N(t) is the number of days in the month. Annual primary production was converted to annual nitrogen fixation using a C:N ratio of 7.8 (mol:mol; Boatman, Mangan, et al., 2018), which corresponds to the present day (mid CO 2 ) condition. (2)

| The role of temperature and pCO 2 in Trichodesmium biogeography
Trichodesmium's fundamental niche is principally set by its thermal niche width (w), which is the range between the minimum (T min ) and maximum (T max ) temperature tolerance limits for growth. Our growth rate model estimates T min as 19.50°C and T max as 30.94°C, with both parameters being independent of CO 2 , irradiance and iron concentration (Table S1). Although values for T max as high as 36.5°C have been reported for this species (Breitbarth et al. (2007), data on the in situ distribution of Trichodesmium are consistent with T max close to 31°C (Figure 1).
Modelling Trichodesmium growth as a function of SST and CO 2 reveals a single equatorial species distribution belt during the LGM, with the niche constrained to lower latitudes during the LGM than at present day ( Figure 2). In contrast, Trichodesmium's niche will continue to expand into higher latitudes by the end of this century, and will exhibit a niche reduction in certain equatorial regions due to mean SST's exceeding T max ; this in turn creates a dual equatorial species distribution belt ( Figure 2). Overall, the geographical area corresponding to Trichodesmium's fundamental niche has increased by ~32% from the LGM to present day, but will decrease by ~9% from present day to the end of this century (Table 1). This estimate relates to the changes in total ocean area only, and does not account for how the area of optimal growth conditions will vary.
Our Trichodesmium growth model shows that the maximal growth rate that can be achieved under low, mid and high CO 2 conditions associated with the LGM, present day and end of this century are 0.204, 0.324 and 0.357 day −1 respectively. It also estimates that the optimal temperature for Trichodesmium growth (T opt ) is 25.85°C.
At this optimal temperature, and under the low CO 2 condition associated with the LGM, regions where SST was sufficient to support a growth rate >0.15 day −1 (a rate equal to 72% of μ max under these conditions) were 1.7 times lower than at present, but are projected to decrease by 10% from the present day to the end of this century (Table S4). Our results also show that at the optimal temperature and under present day and future CO 2 conditions, the area where SST allows Trichodesmium to grow at >0.25 day −1 (>70% and 75% of maximum growth rate under optimal temperature at mid and high CO 2 levels) will not significantly change by the end of this century (Table S4).

| The role of temperature, irradiance and pCO 2 in Trichodesmium biogeography
Incorporating the effect of irradiance into calculations based on our SST and CO 2 model did not change the spatial distribution of Trichodesmium at the LGM, present day and end of this century, but did reduce areas associated with high growth rates (Figure 3). Our Trichodesmium growth model estimates that the optimal irradiance for Trichodesmium growth (E opt ) is 320 µmol photons m −2 s −1 , and is independent of temperature, CO 2 and Fe. As such, under the low, mid and high CO 2 conditions associated with the LGM, present day and end of this century, maximum in situ growth rates were still found to be 0.204, 0.324 and 0.357 day −1 respectively. Regions where both SST and irradiance were sufficient to support a growth rate >0.15 day −1 were 1.9 times lower at the LGM than at present, but are projected to decrease by 10% from the present to the end of this century (Table S5). In contrast, our results also show that under current and projected future CO 2 conditions, the area where SST and irradiance allow Trichodesmium to grow at >0.25 day −1 (>70% and 75% of maximum growth rate under optimal temperature at mid and high CO 2 levels) is projected to increase by 7% by the end of this century (Table S5).

F I G U R E 2
The distribution of Trichodesmium erythraeum IMS101 growth rate (day −1 ) calculated for the mixed layer as a function of sea surface temperature and CO 2 , assuming optimal irradiance and iron-replete conditions. Maps were generated for February (a-c) and August (d-f) during the last glacial maximum (LGM; a, d), for the present (b, e) and projected for the future (est. 2100; c, f). Note, a delta growth rate map for (Present-LGM) and (Future-Present) is presented in Figure S5 3.3 | The role of temperature, irradiance, iron and pCO 2 in Trichodesmium biogeography Incorporating the effect of iron into calculations based on our SST, irradiance and CO 2 model did not change the spatial distribution of Trichodesmium at the LGM, present day and end of this century, but did cause an even greater reduction of areas associated with high growth rates (Figure 4). Our model estimates that the half saturation concentration of iron for Trichodesmium growth (K m ) is constant at 185.5 pM, and is independent of CO 2 , and irradiance at a temperature of 26°C. We did not assess the interaction of Fe limitation with temperature in our experiments, a point that we discuss below. Modelling growth rates onto maps of SST, irradiance and iron, showed that the highest achieved in situ growth rates were lower than the maximum Fe-replete rates (i.e. 0.172 day −1 vs. 0.204 day −1 at low CO 2 , 0.269 day −1 vs. 0.342 day −1 at present CO 2 and 0.305 day −1 vs. 0.357 day −1 at high CO 2 conditions). Regions where SST, irradiance and iron were sufficient to support a growth rate >0.15 day −1 were almost an order of magnitude lower at the LGM than at present, but are not projected to change substantially (decrease by 1%) by the end of this century (Table S6). Our results also show that under present and projected future CO 2 conditions, the area where SST, irradiance and Fe are sufficient to allow Trichodesmium to grow at >0.25 day −1 is projected to increase by 173% by the end of this century (Table S6).

| Temperature dependence of iron saturation
Our Trichodesmium growth model (Equation 1) assumes a multiplicative interaction of the effects of temperature, CO 2 , irradiance and Fe on growth rate. The experimental data used for parameter estimation included three CO 2 levels at both a light-saturating and light-limiting irradiance at 21 temperatures in the range 19°C to 32°C giving us high confidence that a multiplicative interaction of these variables is correct. However, the effect of Fe limitation on growth was only determined at the optimum temperature for growth of 26°C.
Although our data at this temperature indicate a multiplicative interaction of the effects of CO 2 , irradiance and Fe on growth rate, we did not measure the Fe dependence of growth at other temperatures.

TA B L E 1
The geographical areas (M km 2 ) of Trichodesmium's fundamental niche (µ > 0 day −1 ) and regions of optimal growth conditions (µ > 0.25 day −1 ) at the last glacial maximum (LGM), the present and projected for the future (est. 2100). Ocean area associated with growth rates (µ > 0.25 day −1 ) were calculated as a function of: (i) sea surface temperature (SST) assuming both irradiance and Fe are not limiting, (ii) SST and irradiance assuming Fe is not limiting and (iii) SST, irradiance and iron concentration. Note, the maximum growth rate at the LGM was 0.204 day −1 ; therefore, the geographical areas were always zero. The projected ocean area for each month, calculated at varying growth rate thresholds can be found in Tables S3-S6. Values in parenthesis are calculated using a temperature-dependent K m (

F I G U R E 3
The distribution of Trichodesmium erythraeum IMS101 growth rate (day −1 ) calculated for the mixed layer as a function of sea surface temperature, CO 2 and irradiance, assuming iron-replete conditions. Maps were generated for February (a-c) and August (d-f) during the last glacial maximum (LGM; a, d), for the present (b, e) and projected for the future (est. 2100; c, f). Note, a delta growth rate map for (Present-LGM) and (Future-Present) is presented in Figure S6 In our model, we assume that Trichodesmium's K m for Fe-limited growth is independent of temperature, which is consistent to the observation that the relationship between the relative degree of iron limitation of growth rate and the concentration of biologically available dFe varied little with temperature in the diatom Thalassiosira pseudonana (Sunda & Huntsman, 2011). In contrast to this assumption, Jiang et al. (2018) reported that Trichodesmium exhibits a temperature dependency of K m ; with values decreasing by up to 84% from the value at the optimal growth temperature at sub-and supra-optimal temperatures. Jiang et al. (2018) concluded that this temperature dependence of K m could result in a large increase in growth and N 2 fixation under the SST projected for the end of this century, with the global marine N 2 fixation rates increasing by ~22%.
We assessed the effect that using a temperature dependence of K m similar to that reported by Jiang et al. (2018; Figure S8) has in our growth rate model. When we did this, ocean area with optimal SST, optimal irradiance and iron (µ > 0.15 day −1 ) were 8% lower in the present and 37% greater in the future when a temperature-dependent K m instead of a constant K m (see Table 1). Furthermore, regions where SST, irradiance and iron were sufficient to support a growth rate >0.25 day −1 were 90% lower in the present and 29% greater in the future when a temperature-dependent K m was used instead of a constant K m (see Table 1). Thus, regions where Trichodesmium can achieve high growth rates are projected to increase in area by 787% from the present to the end of this century when a temperature dependent K m is used (Table S7), much greater than the 173% increase that we calculated using a constant temperature invariant K m (Table S6).

| D ISCUSS I ON
Our model of the dependence of Trichodesmium growth rates on four key abiotic variables (temperature, irradiance, pCO 2 and iron; Figure 1a-e) was based on previous experiments Boatman, Oxborough, et al., 2018), where the difference in growth rate between experiments at identical growth conditions was <5%. The growth rate data were obtained for cultures that were in balanced growth. Other studies report somewhat different thermal tolerance limits (Boyd et al., 2013;Breitbarth et al., 2007), and CO 2 dependencies of growth rate (Eichner, Kranz, & Rost, 2014;Shi, Kranz, Kim, & Morel, 2012); however, these were obtained from experiments with lower data resolution at the temperature limits for growth, and short culture acclimation times. Overall, we have confidence in our previous studies and Equation (1), which is based on those previous results.

| A representative model for the Trichodesmium genus
A notable consideration of our Trichodesmium growth model (Equation 1) is that it is based on data for Trichodesmium erythraeum, but the Trichodesmium genus includes six species assigned to four clades (Lundgren, Janson, Jonasson, Singer, & Bergman, 2005); although the lower and upper temperature limits that we measured for T. erythraeum IMS101 are nearly identical to the temperature limits currently observed for Trichodesmium spp. in nature from shipboard samples (Figure 1f). Recent studies indicate that other species including T. thiebautii (Rouco, Warren, Mcgillicuddy, Waterbury, & Dyhrman, 2014) and T. tenue (Chappell & Webb, 2010) are more abundant in nature than T. erythraeum. Unfortunately, these other species have not been subjected to the same rigorous laboratory investigations of their growth requirements. In addition, our model does not account for the potential evolution of Trichodesmium to elevated temperature and/or CO 2 as the ocean warms. For example, previous research has indicated that prolonged (~6.5 years) exposure to future elevated CO 2 concentrations causes Trichodesmium to F I G U R E 4 The distribution of Trichodesmium erythraeum IMS101 growth rate (day −1 ) calculated for the mixed layer as a function of sea surface temperature, CO 2 , irradiance and iron concentration. Maps were generated for February (a-c) and August (d-f) during the last glacial maximum (LGM; a, d), in the present (b, e) and projected for the future (est. 2100; c, f). Note, a delta growth rate map for (Present-LGM) and (Future-Present) is presented in Figure S7 significantly, and irreversibly increase N 2 fixation and growth rates .

| Trichodesmium's fundamental niche under nutrient saturation
Based on SST alone, Breitbarth et al. (2007) projected an 11% increase in the spatial distribution of Trichodesmium by the end of this century, but with a 16% decrease in areas they defined as optimal for growth. Our results do not support these findings, instead showing a 9% decrease in the spatial distribution by the end of this century, with no significant difference (−0.2%) in areas with optimal SST for growth. The causes of these contradictory findings are: (a) a difference in the predicted global maps of oceanographic SST, (b) a difference in Trichodesmium's maximum temperature limit for growth while we opted to use a growth rate threshold (µ > 0.25 day −1 ).
Incorporating irradiance into our SST and CO 2 model led to a projected 7% increase in areas defined as optimal for growth from the present to the end of this century. This increase is due to a projected increase in water stratification and a shallower mixed layer ( Figure S9). This in turn reduces the supply of nitrate from deep water into the surface waters, which in turn leads to lower phytoplankton biomass, as indicated by the modelled reduction in surface chlorophyll a ( Figure S10) and dissolved organic carbon concentrations ( Figure S11). The decrease in biomass causes a decrease in light attenuation ( Figure S12), resulting in higher mean irradiances within the mixed layer, and a 7% increase in Trichodesmum's growth potential.

| The role of nutrient limitation in Trichodesmium biogeography
Although temperature and irradiance are the principal factors that determine the potential geographical range of Trichodesmium, iron is a major factor that determines growth rates. This is because Trichodesmium has a high cellular requirement for iron due to an extensive suite of metalloenzymes, including nitrogenase, proteins in the photosynthetic reaction centres and ferredoxin (Kustka et al., 2003), and to the low solubility and short residence time of Fe in seawater (Liu & Millero, 2002).

Low iron concentrations limit
Trichodesmium's productivity and growth in many ocean regions  (Table S6). In addition, many regions within the geographic range that are close to the optimal temperature for nutrient-replete growth, exhibited low or no growth when iron was included into the model (Figure 4). Regions where growth remained relatively high included waters off the west coasts of Africa and Australia; most likely due to the contribution of dust deposition of iron into the surface waters. Aeolian dust has been shown to benefit Trichodesmium by providing a source of bioavailable iron (Basu, Gledhill, De Beer, Matondkar, & Shaked, 2019;Basu & Shaked, 2018).
We applied the same constants for the half saturation concentration for Fe and fraction of bioavailable dFe across all seasons and timescales (i.e. LGM, Present, Future). Thus, changes in growth potential over time, notably the 173% increase in the area associated with >0.25 day −1 by the end of this century, are solely driven by changes in the dissolved iron concentration (dFe; Figure S13). The IPSL-CM5A-LR earth system model uses Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES; Aumont & Bopp, 2006) to simulate all major nutrients including iron. In terms of inputs, atmospheric deposition is estimated from Integrated Nitrogen Catchment model (Aumont, Bopp, & Schulz, 2008); river discharge of carbon and nutrients is taken from Ludwig, Amiotte Suchet, and Probst (1996); and iron input from sediment mobilizsation as parameterized in Aumont and Bopp (2006). These sources are explicitly included but do not vary in time apart from a climatological seasonal cycle for the atmospheric input. Thus, the increase in dFe modelled from present day to the end of the century is not due to iron inputs, rather the relationship between pools. For example, a higher dFe could be due to less complexed iron, reduced demand by phytoplankton growth, increased zooplankton exudation or reduced loss of particulate iron from sinking sediment whether by higher remineralization, reduced scavenging/aggregation or lower bacterial uptake.
It is worth noting that the relationship between dFe and complexed iron within PISCES uses a basic description of iron-ligand interactions, and does not account for how future conditions (i.e. ocean acidification) will alter the organic chelation of iron. We acknowledge that this is a shortcoming in our projections of the future growth potential of Trichodesmium given the importance of dFe data to our iron-integrated modelling outcomes, and highlight this is an area requiring further work in order to improve projections of future ocean trace metal chemistry. We did not consider the potential role of phosphorus limitation in our model. Diazotrophs are reported to be iron limited in the Pacific and Indian Ocean, and phosphorus limited (Misumi et al., 2014) or phosphorus-iron co-stressed (Held et al., 2020;Mills, Ridame, Davey, La Roche, & Geider, 2004) in the Atlantic Ocean. A reduction in dust deposition could have major implications for the phosphorus pool as well as iron. While phosphorus (P) was not integrated into our growth rate model, Trichodesmium may be able to alleviate P limitation by utilizing organic sources of phosphorus (Dyhrman et al., 2006;Sohm & Capone, 2006). Alternatively, Trichodesmium colonies larger than 1 mm may alleviate P limitation by vertically migrating below the phosphocline to assimilate phosphate and rising back to the surface of the euphotic zone (White, Spitz, Karl, & Letelier, 2006). A recent study by Garcia, Fu, Sedwick, and Hutchins (2015) showed that under P deficiency, Trichodesmium cells grew and fixed nitrogen faster with concurrent iron limitation than when iron was replete. This could have significant implications for Trichodesmium's growth potential, particularly in the Atlantic Ocean, with more emphasis on SST and irradiance in determining the fundamental niche.

| Trichodesmium production
Calculating Trichodesmium's global primary production requires information on the oceanographic distributions of both biomass and growth rate. Recently, a data-driven map of the global Trichodesmium nifH abundance was provided by Tang and Cassar (2019). In this paper, a large but still geographically limited database of Trichodesmium nifH gene abundance was extrapolated to the global ocean. Here, we estimate Trichodesmium's contribution to present day ocean production by converting Tang and Cassar (2019) modelled nifH abundance to carbon biomass, and then multiplying by our map of mean growth rate ( Figure 5).
We estimate that Trichodesmium fixes 0.347 Pg C/year, which represents a contribution of 0.58%-0.96% of the total ocean primary production of 36-60 Pg C/year (Carr et al., 2006). Assuming  (Luo et al., 2012), Trichodesmium would account for 84% of global diazotrophy. These estimates were made using a constant K m value for iron-limited growth; however, using a temperature-dependent K m value increases primary production to 0.380 Pg C/year and N 2 fixation to 56.8 Tg N/year, which represents 0.64%-1.05% and 91.6% of the total annual mean C and N productivity respectively (Table S8).
Although such an approach was employed for the present data, we did not apply it to the LGM, or to project into the future as the current relationship between Trichodesmium abundance and oceanographic variables will most likely breakdown for other climate regimes. This is because the prevailing climate regime will affect not only ocean circulation but also competition of Trichodesmium with other phytoplankton, including other diazotrophs, as well as trophic interactions with grazers and pathogens. It is also worth noting that the growth rates used in these productivity calculations were dependent on temperature, light and iron only, and do not include rate limiting factors such as competition and trophic interactions. As such, these estimates should be considered as the upper limit for Trichodesmium's relative contribution to the C and N cycle.

| CON CLUS ION
Our analysis indicates that the increase in SST from the LGM to present has allowed Trichodesmium's range to expand to higher latitudes. Future increases in SST from the present to the end of this century are projected to cause a range expansion at high latitudes and a range contraction in the tropics (Figure 2). The former is driven by the poleward shift of the 20°C isotherm, while the latter will be due to the SST exceeding the maximum thermal tolerance limit for growth (31°C). We also found that although the increase in pCO 2 from the LGM to the present has allowed growth rates to nearly double under nutrient-replete conditions (Figure 2), future increases in CO 2 will have little direct effect on the growth rate of Trichodesmium.

F I G U R E 5
The annual mean growth rate (day −1 ) (a), modelled nifH gene abundance (log 10 copies/m 2 ) (b), biomass (mg C/m 2 ) (c) and primary production (g C m −2 year −1 ) (d) of Trichodesmium erythraeum IMS101. Growth rates and biomass values were averages across all months and were multiplied to generate the map of primary productivity. Note, global monthly mean primary production values are reported in Table S8 Increased water stratification and a shallower mixed layer in the future ocean will limit the supply of nitrate to the surface waters. Since these stratified waters are already N-limited, this should benefit Trichodesmium and other diazotrophs, as competition with N-limited phytoplankton for Fe and P may decrease. It is also likely that the surface Chl a concentration will decline as more of the ocean becomes increasingly N-limited (Zehr & Kudela, 2011); this should in turn decrease light attenuation leading to a higher mean irradiance within the shallower mixed layer. Higher irradiance will increase the rate of iron photoreduction. This coupled to the increase in Fe(II) caused by ligand dissociation under more acidified conditions (Shi et al., 2010) should increase the bioavailablity of Fe to Trichodesmium.
Such future changes could significantly increase Trichodesmium's global productivity (Hutchins et al., 2007), which, given the relative contribution to global annual mean N production, suggests a potentially significant negative feedback to the increasing atmospheric CO 2 concentrations that have been and continue to be caused by fossil fuel burning and deforestation.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available in the supplementary material of this article.