Toward Determining the Spatio-Temporal Variability of Upper-Ocean Ecosystem Stoichiometry From Satellite Remote Sensing

Introduction

Ever since Redfield first reported on the C:N:P ratio of particulate organic matter (POM) more than 85 years ago (Redfield, 1934), the ratio has been widely assumed to be stable. A fixed C:N:P ratio has long played a central role in ocean biogeochemistry because this ratio largely determines the strength of the biologically mediated ocean carbon cycle. However, recent studies show convincingly that the C:N:P stoichiometry of POM varies substantially on ocean-basin scales. For example, Martiny et al. (2013a) showed a globally coherent pattern, with C:N:P ratio of 195:28:1 in the subtropical gyres, 137:18:1 in the warm upwelling zones, and 78:13:1 in the nutrient-rich polar regions. An inverse model of ocean biogeochemistry also inferred a similar spatial pattern of the global C:P and N:P ratios (Teng et al., 2014; Wang et al., 2019).

As carbon export is inversely related to atmospheric CO₂ (Volk and Hoffert, 1985), carbon-enriched particulate organic matter in subtropical gyres could lead to lower atmospheric CO₂ and higher export production of carbon, thereby influencing climate (Galbraith and Martiny, 2015; Tanioka and Matsumoto, 2017; Matsumoto et al., 2020a; Ödalen et al., 2020). The ocean carbon modeling community is beginning to respond to this development. For example, the state of the art CMIP5/6 models developed by various climate modeling teams around the world represent phytoplankton stoichiometry with varying degree of flexibility, from no flexibility (i.e., fixed C:N:P ratio) to fully flexible (e.g., Bopp et al., 2013; Arora et al., 2020).

A major challenge to adopting fully flexible stoichiometry in biogeochemical models is our current inability to observationally constrain the temporal variability of the C:N:P in the global ocean. Although some progress has been made to explore a temporal shift in C:N:P using local time-series data (Hebel and Karl, 2001; Karl et al., 2001; Singh et al., 2015; Martiny et al., 2016; Talarmin et al., 2016), our holistic global view of the global C:N:P ratio variation is still unclear. In situ C:N:P measurements of POM inherently suffer from bias toward regions and periods of active oceanographic research, and large parts of the global ocean remain devoid of data. For example, there is a considerable paucity of POM sampling efforts in the South and Equatorial Atlantic regions (Sharoni and Halevy, 2020).

Satellite ocean-color sensors have the potential to provide a unique tool to constrain the temporal evolution of organic matter C:N:P ratio. Ocean color provides global, synoptic views of the spectral remote-sensing reflectance of the ocean that can be used to generate estimates of marine inherent optical properties (IOPs) at various timescales (Werdell et al., 2018). Satellite ocean color (i.e., remote-sensing reflectance) provides an unparalleled tool to capture climate-driven signals in the upper biological functions of the global ocean (Dierssen, 2010; Dutkiewicz et al., 2019), and has the potential to yield crucial information on the modes of C:N:P variability. Previous field studies have shown that C:N:P ratio is significantly influenced by interannual climate variabilities such as ENSO and Pacific Decadal Oscillation (Martiny et al., 2016; Fagan et al., 2019).

One possible approach to assess the spatio-temporal variability in the C:N:P of POM is to directly estimate the change in the total concentration of particulate organic carbon (POC), particulate organic nitrogen (PON), and particulate organic phosphorus (POP) using satellite ocean color data. Multiple methods of estimating total POC from satellite ocean color have been developed over the years, and the satellite estimates are extensively calibrated with in situ measurements (Evers-King et al., 2017; Rasse et al., 2017). More recently, Fumenia et al. (2020) have developed a method to link the backscattering coefficient (b_bp) at 700 nm with PON and POP concentrations in the oligotrophic Western Tropical South Pacific. However, the reliability of b_bp as a quantitative proxy of PON and POP still needs to be investigated in other oceanographic areas, including non-oligotrophic regions.

Another possible approach of deriving C:N:P of bulk POM is to predict phytoplankton’s elemental composition and use it as a proxy for the bulk composition, assuming phytoplankton makes up the largest proportion of POM. The study by Arteaga et al. (2014) showed a seasonally variable global C:N:P ratio of phytoplankton by using a combination of remote sensing data and a mechanistic growth-model of phytoplankton (Pahlow et al., 2013). More recently, Roy (2018) developed a method to estimate the macromolecular content of phytoplankton protein, carbohydrate, and lipid via satellite ocean color by using empirical relationships between the particulate backscattering coefficient, phytoplankton cell size, and cellular macromolecular concentrations. However, this method cannot derive phytoplankton C:P as there is no empirical link between cell size and P-rich macromolecules such as RNA and DNA (Finkel et al., 2016; Tanioka and Matsumoto, 2020b). Furthermore, a fundamental limitation in both of these studies is that the elemental composition of phytoplankton may not be able to explain the full dynamics of bulk POM because, in reality, phytoplankton biomass typically constitute only 30–50% of bulk particulate organic matter in the open ocean (Eppley et al., 1992; Durand et al., 2001; Gundersen et al., 2001; Behrenfeld et al., 2005).

Here, we propose a new remote-sensing approach that uniquely combines established methodologies in order to quantify the spatio-temporal variability of the upper-ocean stoichiometry of phytoplankton and bulk POM (Figure 1). We demonstrate that this method specifically for the C:P ratio as C:P is more commonly compared to observations than C:N and N:P under different environmental conditions (e.g., Sterner et al., 2008). Furthermore, the C:P is a key variable for converting phosphorus-based fluxes to carbon-based fluxes in biogeochemical models, and therefore has important implications for the carbon cycle (e.g., Matsumoto et al., 2020b). The framework, however, can theoretically be expanded to include C:N and N:P ratios. In this approach, we first determine C:P of phytoplankton by combining satellite-derived estimates of growth rate, Chl:C ratio, and nutrient depletion temperatures (NDTs) with a newly developed mechanistic model of phytoplankton stoichiometry (Inomura et al., 2020). We then convert phytoplankton C:P ratio to the total POC:POP using remotely sensed concentrations of phytoplankton biomass and POC. This approach is unique in that all inputs are derived from satellite remote sensing and does not rely on in situ measurements, thereby enabling us to predict the “real-time” evolution of phytoplankton and bulk POM C:P on various temporal and spatial scales of interest.

The relative importance of the two main drivers of POC:POP variability are discussed following their (1) variability due to change in phytoplankton C:P that reflect changes in environmental condition such as nutrient supply (e.g., Martiny et al., 2013a; Garcia et al., 2018), and (2) variability due to change in community plankton composition (e.g., Weber and Deutsch, 2010; Talmy et al., 2016; Sharoni and Halevy, 2020). Finally, we discuss caveats, limitations, and future directions. Our ultimate goal in this paper is to demonstrate the method’s feasibility, given all the assumptions and limitations. We envision that future advances in satellite instrumentation will enhance the accuracy of satellite-derived input parameters and improve the overall estimate of C:N:P from space.

Materials and Methods

Satellite-Informed Modeling Framework

The flowchart shown in Figure 1 provides an overview of how we determine phytoplankton C:P and bulk POC:POP ratios from satellite products. In the sections below, we briefly describe the phytoplankton stoichiometry model and the method of estimating the bulk C:P of POM.

Phytoplankton Stoichiometry Model

In this study, we determined the C:P ratio for a single phytoplankton functional type using a recently developed phytoplankton stoichiometry model (Inomura et al., 2020). The phytoplankton stoichiometry model of Inomura et al. (2020) is conceptually simple but facilitates the accurate computation of phytoplankton C:P and C:N ratios under a variety of environmental conditions. The input variables required in calculating phytoplankton C:P are light intensity, growth rate, and the presence/absence of limiting nutrients. The model is based on four empirically supported lines of evidence: (1) a saturating relationship between light intensity and photosynthesis, (2) a linear relationship between RNA-to-Protein ratio and growth rate, (3) a linear relationship between biosynthetic proteins and growth rate, and (4) a constant macromolecular composition of the light-harvesting machinery. Also, it follows from these assumptions that chlorophyll-to-carbon ratio (Chl:C_phyto) and growth rate are directly linked for any given light intensity (Laws and Bannister, 1980). Inomura et al. (2020) calibrated their model parameters subject to constraints from published laboratory chemostat studies for several key prokaryotic and eukaryotic phytoplankton species. For this study, we used the model parameter set for the cyanobacteria Synechococcus linearis because the parameters for this species were most rigorously calibrated with laboratory data compared to the other two possible options (cf. a diatom, Skeletonema costatum, and a haptophyte, Pavlova lutheri). Also, picocyanobacteria such as Synechococcus and Prochlorococcus are the most abundant phytoplankton types in the global ocean (Flombaum et al., 2013; Berube et al., 2018). Thus, if we are choosing a single group of phytoplankton to represent the whole phytoplankton community, as we do in this study, Synechococcus would be a reasonable choice. However, as this particular species of Synechococcus is a freshwater species, further calibration efforts specific to the marine cyanobacteria species would be necessary. Although the ecology of marine and freshwater phytoplankton is fundamentally similar (Kilham and Hecky, 1988), a previous meta-analysis study has shown that the interactive effects of multiple environmental stressors (e.g., temperature, light, and nutrients) on C:N:P may be significantly different for marine and freshwater species (Villar-Argaiz et al., 2018). A complete description and evaluation of the phytoplankton stoichiometry model are provided in the original model description paper (Inomura et al., 2020).

In order to determine phytoplankton C:P, we made three minor modifications to the original stoichiometry model by Inomura et al. (2020). First, we drove the stoichiometry model directly with depth-integrated Chl:C_phyto in the mixed layer obtained from the satellite ocean color instead of calculating Chl:C_phyto as a function of photon-flux density. This way, we could circumvent the need to estimate depth-dependent irradiance, which is complicated by issues such as self-shading and particle scattering (Jamet et al., 2019). Second, we imposed a fixed maximum growth rate of 2 d^–1 in calculating C:P, which is equal to the maximum growth rate commonly imposed on the satellite-based estimates of growth rate (Westberry et al., 2008; Laws, 2013). Third, we imposed a constant C:P value of 102 under the P-replete condition regardless of the P supply. A C:P of 102 corresponds to the maximum capacity of cellular P ($Q_P^{\mathit{max}}$), and was used in the original stoichiometry model by Inomura et al. (2020). This assumption for Synechococcus is supported by culture experiments (e.g., Healey, 1985) and makes it possible to circumvent the need to calculate C:P based on the absolute values of external nutrient supply. Under P limitation, all cellular P tends to be allocated toward DNA, RNA, and phospholipids in the thylakoid membranes. Therefore, by default, the stoichiometry model does not require information on external nutrient concentration in calculating cellular C:P. With these three modifications, we were able to predict phytoplankton C:P using only satellite ocean color products as inputs.

Satellite-Derived Inputs

We drove the modified Inomura model with satellite-derived phytoplankton growth rates (μ), Chl:C_phyto (a measure of light intensity), and P limitation [as the difference between Sea Surface Temperature (SST) and cubic root-corrected phosphate depletion temperature (PDT3)] to estimate phytoplankton C:P (r_C:P) in the surface mixed layer (Eq. 1):

The required input data in Eq. 1 are monthly binned and averaged observations from the Aqua Moderate Resolution Imaging Spectroradiometer (MODIS-Aqua) acquired from January 2003 to December 2018 and re-gridded on a regular 1°-latitude by 1°-longitude grid. All satellite-derived input data and estimates of mixed-layer depth are available for download from the Oregon State Ocean Productivity Website¹.

The carbon-based specific growth rate μ (measured in d^–1) is estimated by dividing the depth-integrated net primary productivity (NPP, measured in mg C m^–2 d^–1) by the standing stock of phytoplankton carbon (C_phyto, measured in mg C m^–2):

Growth rate is influenced by light, temperature, and nutrient availability, and previous culture measurements (e.g., MacIntyre et al., 2002) provide robust and well-understood empirical relationships between μ and satellite products (Silsbe et al., 2016). There are multiple NPP data products available to date (Westberry and Behrenfeld, 2014; Bisson et al., 2018). In order to illustrate the robustness of our C:P determination to the choice of the NPP products, we used the following four NPP satellite data products: (1) the Carbon, Absorption and Fluorescence Euphotic-resolving model (CAFE) (Silsbe et al., 2016), (2) the Vertically Generalized Productivity Model (VGPM) (Behrenfeld and Falkowski, 1997), (3) the Eppley-VGPM Model (Eppley, 1972; Behrenfeld and Falkowski, 1997), and (4) the Carbon-based Productivity Model (CbPM) (Westberry et al., 2008). A previous study showed that CAFE compares best with in situ NPP measurements (Bisson et al., 2018). Because the growth rates from VGPM, Eppley-VGPM, and CbPM are similar quantitatively (Supplementary Figure 1), we only present results from VGPM as representing the three models in the main text. Throughout the text, we use the phrases “CAFE-informed phytoplankton C:P” and “VGPM-informed phytoplankton C:P” to refer to C:P calculated using μ from CAFE-based NPP and VGPM-based NPP, respectively.

For C_phyto, the satellite data product of Westberry et al. (2008) was used, who computed C_phyto as a linear function of the particulate backscatter coefficient at 443 nm, b_bp(443). We only considered a single algorithm of C_phyto in this study because the previous intercomparison study showed that no single algorithm outperforms any of the other algorithms when compared with in situ data (Martínez-Vicente et al., 2017). We excluded from our analyses the coastal regions with C_phyto exceeding 1,000 mg C m^–3 and we multiplied the monthly mean surface concentration of C_phyto with monthly mean mixed layer depth (MLD) from the Hybrid Coordinate Ocean Model (HYCOM) to get the depth-integrated C_phyto. Here, MLD is defined as the depth where the density of water is greater than that of water at a reference depth of 10 m by 0.125 kg m^–3 (Levitus, 1982). The growth rate calculated this way in Eq. 2 is representative of a well-mixed, photoacclimated community subject to the median PAR in the mixed layer. The satellite-derived seasonal variability in μ reflect changes in light and nutrient limitation, as well as phytoplankton community composition (Behrenfeld et al., 2005).

Supplementary Figure 1 shows satellite-derived estimates of μ during summer and winter. CAFE predicts a higher μ during summer months compared to winter months for the large parts of the ocean (Supplementary Figures 1A–C). VGPM and the other two NPP products (CbPM and Eppley-VGPM) show similar trends at high latitudes but show the opposite trend in the subtropics with lower μ during summer compared to winter. As a result, the range of estimated μ amongst NPP products are higher during the summer compared to winter and is most extensive in the subtropics. Throughout the rest of this paper, the “summer” average refers to average values during July–September in the Northern Hemisphere and January–March in the Southern Hemisphere. For the “winter,” the target months are reversed between two hemispheres.

The Chl:C_phyto ratio, a proxy for light limitation (Falkowski et al., 1985; MacIntyre et al., 2002), is computed here by dividing MODIS-derived Chl-a with C_phyto. Chl-a concentration is depth-integrated and therefore converted from mg Chl m^–3 to mg Chl m^–2 by multiplying the monthly mean surface concentration with monthly mean MLD. Like for growth rate, Chl:C_phyto is assumed to be vertically uniform in the mixed layer. Considering ocean-color measurements are typically representative of the first optical depth of the surface ocean (Volpe et al., 2012; Bellacicco et al., 2018), this assumption is widely used in satellite remote sensing models, including the four NPP models used here. However, this assumption overlooks depth-dependent changes in Chl:C_phyto, for example when the deep chlorophyll maximum (DCM) becomes much shallower than the mixed layer depth (Cullen, 1982). Supplementary Figures 2A–C shows estimates of Chl:C_phyto during summer and winter. In general, Chl:C_phyto is higher during winter than summer as the reduced incident irradiance causes phytoplankton to allocate more of the cellular component to the light-harvesting apparatus (Geider, 1987; MacIntyre et al., 2002; Arteaga et al., 2016). High Chl:C_phyto in the sunlit layer of the continental margins are known to be relatively inaccurate and biases due to interferences by the high and variable amounts of colored dissolved organic matter (CDOM) and detritus (Siegel et al., 2005; Morel and Gentili, 2009; Loisel et al., 2010). As we excluded coastal regions in the subsequence analyses, this issue should not affect our satellite-informed C:P estimates.

P limitation was assessed by utilizing nutrient depletion temperatures (NDTs), which are temperatures above which nutrients are no longer detectable by traditional wet-chemistry techniques (Zentara and Kamykowski, 1977; Kamykowski and Zentara, 1986). The method leverages an observed inverse empirical relationship between surface nutrient concentration and sea-surface temperature (SST). In this relationship, phytoplankton is considered nutrient-limited if the difference between SST and NDT is higher than 0 and vice versa if the difference is lower than 0. We used a global NDT mask of the percentile-based, cubic root-corrected phosphate depletion temperatures (PDT3) re-gridded to a 1-by-1° spatial resolution (Supplementary Figure 2F; Kamykowski et al., 2002). PDT3 was subtracted from MODIS-derived monthly mean SST to determine the absence/presence of P limitation in the surface ocean. P limitation as a result of SST exceeding phosphate depletion temperature is globally prevalent during summer (Supplementary Figure 2D). Phosphate depletion is alleviated during winter months at high latitudes and in some parts of the equatorial regions as the surface ocean cools in part because of enhanced vertical mixing (Supplementary Figure 2E). For the current work, we limited our study to latitudes ranging from 50°S to 70°N as the original data on PDT3 beyond this latitudinal range are sparse (Kamykowski et al., 2002). We further discuss the caveats and limitations of this approach in section “Caveats, Limitations, and Future Needs.”

MODIS-derived total monthly averaged POC (0.7 μm < D < 17 μm) was obtained from the NASA Ocean Color Product webpage². This total POC determination is based on an empirical relationship between POC and the blue-to-green band of spectral remote-sensing reflectance (Stramski et al., 2008). The algorithm employed here is widely implemented for producing maps of surface POC. The global mean C_phyto:POC is ∼30% (Supplementary Figures 2G,H), consistent with previous estimates (Behrenfeld et al., 2005). The C_phyto:POC is generally higher in the subtropical gyres than other regions reaching up to 50–70% during summer (Supplementary Figure 2G). C_phyto:POC ratio rarely exceeds a value of 1 except during episodic events in coastal regions, which we disregard in our analyses. Although C_phyto and POC are independently determined, the fact that C_phyto:POC ratio rarely exceeds a value of 1 increases our confidence in the predictability of C_phyto:POC.

Estimating C:P of Bulk POM

Globally, phytoplankton-derived organic matter represents on average ∼30% of bulk organic matter (Eppley et al., 1992; Durand et al., 2001; Gundersen et al., 2001; Behrenfeld et al., 2005), and the rest is due to contributions from zooplankton and non-living detrital materials. A previous satellite-based study showed that phytoplankton carbon biomass makes up 30–70% of the total POC pool in the tropical oligotrophic regions and ∼10–30% in higher latitudes and equatorial Pacific (Arteaga et al., 2016). The relative contribution from phytoplankton biomass is lower in nutrient-rich regions, possibly due to the top-down control on phytoplankton biomass by zooplankton (Ward et al., 2014; Talmy et al., 2016), resulting in higher zooplankton fraction and smaller phytoplankton fraction in the total POC pool.

In order to estimate C:P of bulk POM, we split the POC and particulate organic phosphorus (POP) into two components: (1) phytoplankton-derived organic matter with C:P ratio following the stoichiometry model in the previous section, and (2) non-algal component with fixed C:P of 117:1 following Anderson and Sarmiento (1994). Throughout the rest of this paper, the “community composition” refers to the relative balance between the algal and non-algal components of organic matter, not the community composition of different phytoplankton functional types.

The non-algal component of particulate organic matter with fixed C:P represents a combination of zooplankton and other non-living detrital materials such as fecal pellets and other organic matter left over from sloppy feeding (Martiny et al., 2013a, b; Talmy et al., 2016). Previous studies have shown that zooplankton generally has a C:P close to the Redfield ratio even under P-limited conditions (e.g., Copin-Montegut and Copin-Montegut, 1983; Sterner and Elser, 2002). Isopycnal analysis of export and remineralization stoichiometry of the deep ocean (>400 m) also indicates a relatively constant C:P of around ∼117 globally (Anderson and Sarmiento, 1994).

In calculating the C:P ratio of bulk POM, we solve for three unknowns: (1) the carbon content of non-algal POM (C_non), (2) the phosphorus content of non-algal POM (P_non), and (3) total POP. This is achieved with three equations:

The subscript “phyto” refers to the phytoplankton component, and “non” refers to the non-algal component of POM. All the quantities are in mol per unit volume. Eqs 3 and 4 describe the conservation of carbon and phosphorus, respectively, and the Eq. 5 describes the fixed C:P ratio of non-algal organic matter. Essentially, Eqs 3–5 constitute a simple two end-member mixing model of the algal and non-algal components. We can obtain C:P of the bulk organic matter as a function of the known quantities from section “Satellite-Informed Modeling Framework,” C_phyto, r_C:P, and total POC by rearranging Eqs 1, 3–5:

Eq. 6 shows that the bulk C:P ratio is a non-linear function of phytoplankton C:P (r_C:P) and the relative abundance of C_phyto over total POC (C_phyto/POC).

Model-Data Comparison of POC:POP

We compared the satellite-informed bulk POC:POP with a recently compiled data set of 5573 in situ observations of suspended oceanic POC:POP ratios from cruises and other marine stations distributed globally (Martiny et al., 2014). The suspended POM samples were collected on 0.7 μm filters (GF/F), and their C:P ratios reflect contributions from phytoplankton, microzooplankton, detrital material, and mixed particle aggregates. We note that the nominal pore size of 0.7 μm cannot fully capture biomass of small bacterioplankton cells (Lee and Fuhrman, 1987; Gundersen et al., 2001), which typically have C:N:P close to the Redfield Ratio (Gundersen et al., 2002). Here, we only used samples from the upper 100 m of the water column, representative of an average mixed layer (Kara et al., 2003) and excluded samples with POP concentrations inferior to the reported detection limit of 5 nM. We also removed samples from coastal waters, which often include a substantial contribution of allochthonous POM (e.g., benthic, riverine) (Liénart et al., 2018).

When comparing the large-scale temporal variability of in situ C:P with satellite estimates, we binned the measured C:P data into 10°-latitude increments. At each sampling station, we calculated the mean C:P in the top 100 m. After this screening process, we were left with 185 observational points for summer and 111 observational points for winter (Figure 2). We compared the seasonally averaged, satellite-informed POC:POP with the C:P of suspended POM spanning from 50°S to 70°N.

To further evaluate the performance of our modeling framework, we compared our satellite-informed estimates of C:P to direct POC:POP measurements at the BATS and HOT sites. The time-series data of POC and POP measurements from these two stations are included in the global POM database. Here, we selected data in the top 100 m that were collected between 2003 and 2010 for the “point-to-point” comparison with the satellite estimates of C:P.

Results and Discussion

Large-Scale Seasonal Variability in Phytoplankton C:P

Combining the estimates of growth rate, Chl:C_phyto, and P limitation can help determine the seasonal variability in r_C:P (Figure 3). The satellite-informed r_C:P is highest in the stratified oligotrophic gyres and lowest in the higher-latitude, seasonally stratified seas and equatorial upwelling regions, consistent with existing field observations (Martiny et al., 2013a). Both the CAFE (Figures 3A–C) and VGPM-informed r_C:P (Figures 3D–F) show elevated r_C:P in the higher-latitude region during the summer months compared to the winter months as ocean warming enhances stratification and phytoplankton becomes P-limited. The increase in light availability during summer, shown by a decrease in Chl:C_phyto, also helps in increasing r_C:P at higher-latitude regions.

Although the spatio-temporal pattern of r_C:P is consistent across four satellite-informed cases for high-latitude regions and equatorial regions (Supplementary Figure 3), the range of the four satellite r_C:P estimates is large in the subtropics (Figures 3G,H). This larger range reveals a relatively large uncertainty in r_C:P in the subtropics. Here, we used range instead of arithmetic standard deviation as a measure of uncertainty because the ratios by their nature do not follow normal distributions (Isles, 2020) and the sample size is small (n = 4). We show the uncertainty of r_C:P expressed in terms of geometric standard deviation from the mean in Supplementary Figure 5. Note that the spatial pattern of geometric standard deviation (Supplementary Figures 5A,B) is very similar to the spatial pattern of range (Figures 3G,H). Considering that the oligotrophic gyres tend to be P-limited throughout the year and the change in Chl:C_phyto is small, large uncertainties in μ are predominantly responsible for this uncertainty in r_C:P in those regions of the global ocean. While the CAFE-informed r_C:P shows a noticeable decrease during summer by ∼100–200 molar units (Figure 3C), VGPM-informed r_C:P shows an increase during summer (Figure 3F).

In theory, r_C:P should decrease as growth rate increases, and the fractional change in r_C:P should be highest for low growth (Droop, 1974; Burmaster, 1979; Goldman et al., 1979; Morel, 1987). In other words, a small change in growth rate should lead to a large change in r_C:P when the growth rate is low. Multiple culture experiments support this prediction, where phytoplankton growing at a high rate is both P-rich and has reduced stoichiometric flexibility (e.g., Hillebrand et al., 2013). If we assume P-limited growth condition and we replace growth rate with PO₄ concentration, this pattern would also be true for PO₄ vs. r_C:P where phytoplankton growing under low P environment are frugal (high r_C:P) and more stoichiometrically flexible (Galbraith and Martiny, 2015; Tanioka and Matsumoto, 2017, 2020a). As subtropical regions are strongly P limited and the growth is suppressed (Wu et al., 2000; Martiny et al., 2019), this reasoning can explain the elevated r_C:P with large uncertainty and sensitivity.

Figure 4 illustrates how r_C:P varies under varying growth rates and Chl:C_phyto in specific regions. Contour lines (isopleths) representing the theoretical values of r_C:P are predicted by the Inomura phytoplankton stoichiometry model for different combinations of μ and Chl:C_phyto under the P limited scenario. In order to illustrate the regional variability, we superimposed monthly averaged, CAFE-informed r_C:P in four oceanographic regions. These four regions are: (1) the high latitude bloom-forming North Atlantic Ocean (NAT: 25°W – 35°W, 45°N – 50°N), (2) the North Atlantic subtropical gyre (NASG: 25°W–70°W, 25°–35°N), (3) the South Pacific subtropical gyre (SPSG: 90°W–150°W, 15°S–40°S), and (4) the Equatorial upwelling region (EQU: 5°S – 5°N), following Westberry et al. (2016). The size of the symbol indicates the extent of P limitation. “P-replete” symbolizes <20% of grid boxes in the region are P-limited, “Moderate” symbolizes 20–80%, and “Deplete” >80% based on the seasonally varying SST. The numbers represent the month of the year.

There are two key features in this plot. The first is that different oceanographic regions occupy a unique space. For example, North Atlantic (NAT) experiences large seasonal variability in growth rate, P limitation, and r_C:P, while EQU experiences small seasonal changes. The second important feature is that the contours representing r_C:P become increasingly close together as the growth rate decreases. This reiterates the fact that a small change in satellite-derived growth rate can lead to a large change in r_C:P in chronically nutrient-deplete subtropical gyres (NASG and SPSG).

Light availability also affects r_C:P as light modulates the cellular allocation between light-harvesting apparatus, biosynthetic apparatus, and energy storage reserves (Falkowski and LaRoche, 1991; Moreno and Martiny, 2018). The Inomura phytoplankton stoichiometry model predicts that increased light limitation increases cellular allocation toward photosynthetic proteins and decreases allocation toward C-rich biosynthetic proteins. Therefore, an increase in Chl:C_phyto (i.e., increased light limitation) will lead to a decrease in r_C:P at a constant growth rate (Figure 4).

As expected, satellite-derived Chl:C_phyto shows maxima during winter months (January-March in Northern Hemisphere and July–September in Southern Hemisphere) due to decreased exposure to sunlight (Figure 4). As shown in previous modeling studies, the effect of light on r_C:P is disproportionally large when the growth rate is low, and an increase in Chl:C_phyto can effectively reduce r_C:P during winter months (Arteaga et al., 2014; Talmy et al., 2014). Compared to the growth rate, however, the effect of light limitation on r_C:P is weak, as shown by the vertically steep contour lines in Figure 4. Indeed, a meta-analysis on published laboratory studies has shown that the effects of light on r_C:P are significantly weaker than that of macronutrients and temperature (Tanioka and Matsumoto, 2020a).

Large-Scale Seasonal Variability in Bulk POC:POP

By combining the satellite-informed phytoplankton C:P and the community composition measured by C_phyto:POC, we can determine POC:POP of the bulk POM (Figure 5). Similar to r_C:P, bulk POC:POP ratios are highest in the gyres compared to the equatorial upwelling and high-latitude regions. Globally, satellite POC:POP is higher during the summer compared to the winter. This seasonal trend can be explained by the higher C_phyto:POC during summer than winter (Supplementary Figure 2I). This makes intuitive sense because the phytoplankton biomass concentration is kept low in the mixed layer during winter months due to the deepening of MLD, strong light limitation, and zooplankton grazing (Behrenfeld and Boss, 2018). As C:P of P-limited phytoplankton is higher than C:P of non-algal organic matter, increase in C_phyto:POC during summer leads to an increase in POC:POP. The most noticeable increase is visible in the South Pacific Subtropic Gyre, where summertime POC:POP is higher than the winter value by ∼200 as C_phyto:POC increases by ∼50% during summer compared to winter. The range (uncertainty) in satellite-informed POC:POP (Figures 5G,H) is much smaller compared to that of phytoplankton C:P (Figures 3G,H), and all the four satellite-informed estimates agree on a general increase in POC:POP during summer compared to winter (Figure 5I and Supplementary Figure 4).

Figure 6A illustrates how the bulk POC:POP is non-linearly related to the community composition (measured by C_phyto:POC) for a given change in r_C:P. We observe from the satellite-derived data of C_phyto and POC that C_phyto:POC is, on average, ∼30% and rarely exceeds 50% of the total POC pool. The increase in POC:POP with respect to increase in r_C:P reaches a plateau quickly when C_phyto:POC < 30%. In other words, the community dominance of non-algal POM over algal POM can effectively put a cap on the increase in bulk POC:POP, even when phytoplankton C:P is very high (e.g., NASG and SPSG). This top-down control on POC:POP due to community composition also explains the low uncertainty in the estimates of satellite-informed POC:POP despite the large uncertainty in satellite-informed r_C:P.

Figure 6B is an alternative way of illustrating this top-down control on bulk POC:POP by community composition. Contour lines representing POC:POP based on our simple two end-member algal/non-algal mixing model are widely separated when C_phyto:POC is low, indicating that POC:POP is relatively stable when C_phyto:POC is relatively low. On the other hand, when C_phyto:POC is high, contour lines become closer together, and bulk POC:POP quickly approaches r_C:P. If we plot monthly averaged estimates of satellite-derived bulk POC:POP under different regions, two distinct clusters become apparent. Subtropical gyres (NASG and SPSG) are characterized by high r_C:P and C_phyto:POC resulting in sizeable seasonal variability in bulk POC:POP. On the other hand, NAT and EQU experience a smaller seasonal change in POC:POP as the C_phyto:POC remains relatively constant around 15%. The take-home message from Figure 6 is that a community composition can exert a strong top-down control on POC:POP even when phytoplankton C:P is much higher than the Redfield ratio. Indeed, multiple studies emphasize this point, including recent studies on C:N (e.g., Talmy et al., 2016), N:P (e.g., Sharoni and Halevy, 2020), as well as the original study by Redfield et al. (1963).

Model-Data Comparison

To assess our model predictions, we first compare our seasonally resolved zonally averaged satellite POC:POP estimates with measurements of sampled POC:POP (Figure 7). Globally, both the satellite estimates and the in situ observations show higher POC:POP in summer (Figure 7A) than in winter (Figure 7B). This increase during summer is likely to be driven by a change in community composition, with an increased C_phyto:POC during summer. At high-latitudes, an increase in phytoplankton C:P also drives an increase in POC:POP during summer. Therefore, the combination of the change in community composition and phytoplankton C:P is responsible for the increased bulk POC:POP during summer.

Although it is promising that our predictions are mostly consistent with observations, there are two distinct regions where the satellite POC:POP and the observations do not agree. The first is the equatorial region during summer (Figure 7A), where satellite-informed POC:POP is around 150 but observed POC:POP is close to the Redfield ratio of 106. This discrepancy stems from the fact that our method likely overestimates the degree to which the equatorial regions are P-limited, which in turn leads to an overestimation of the phytoplankton C:P to as high as ∼200. In addition, our phytoplankton C:P model is tuned to data for Synechococcus. In reality, fast-growing opportunistic eukaryotic plankton such as diatoms and other eukaryotes with lower C:P are more predominant in the equatorial region (Arrigo, 2005; Martiny et al., 2013a; Kostadinov et al., 2016). The second region where we observed a noticeable difference between satellite estimates and in situ observation is around 20°S during winter (Figure 7B). Given the paucity of observations in this region, however, (n = 12 and 8 for summer and winter, respectively), it is challenging to determine the exact cause for the increase in POC:POP during the winter.

To further assess our model capability, we compare time-series data of suspended POC:POP in the top 100 m from BATS and HOT with satellite estimates in the seasonally mixed-layer depth from 2003 to 2010 (Figure 8). It is important to note that suspended POC:POP is a “point” value reflecting elemental composition at a particular location and at a particular time, whereas the satellite-informed POC:POP is a monthly and area-averaged value for a 3-by-3-pixel area around the BATS and HOT stations. We use the median satellite-informed phytoplankton C:P and POC:POP values from four satellite products (CAFE, VGPM, Eppley-VGPM, and CbPM) for comparison with the data. Here, we use the median rather than the mean as the measure of average C:P across satellite products because the ratios by their nature do not follow normal distributions (Isles, 2020).

In general, measured POC:POP ratios lie between our satellite estimates of phytoplankton C:P and bulk POC:POP ratios at both BATS and HOT (Figures 8A,B). Measured POC:POP, on average, is closer to the satellite-informed POC:POP than to satellite-informed r_C:P (Figures 8C,D). This makes intuitive sense because both satellite estimates (Figure 8E) and in situ observations (Casey et al., 2013) show that the biomass of picocyanobacteria (Prochlorococcus, Synechococcus) only contributes to <∼40% of the POC pool in the gyres. Qualitatively, our satellite estimates of bulk POC:POP seem to capture the general seasonal variability, with POC:POP being lowest in the winter and highest in the summer and the fall. Also, both the satellite-informed and the observed C:P are lowest in the late winter as a result of deep mixing and increased supply of nutrients, which cause phytoplankton C:P to decrease (Singh et al., 2015).

The amplitude of POC:POP variability is greater at BATS than at HOT, reflecting larger temporal variability in C_phyto:POC (Figure 8E). Phytoplankton carbon biomass can reach as high as 60% of the total POC pool at BATS during the growing season while the peak C_phyto:POC at HOT does not exceed ∼40%. This difference between the two sites may reflect the fact that the North Atlantic is a bottom-up ecosystem, while the North Pacific is a top-down ecosystem (Steele and Henderson, 1992; Siegel et al., 2016). Satellite-informed bulk POC:POP, however, underestimates the observed POC:POP by ∼50 on average at BATS (Figure 8C) and ∼20 at HOT (Figure 8D), and this may reflect the fact that non-algal organic matter has a higher ratio than Redfield of 117. Also, our satellite-informed estimate may not be fully capturing episodic temporal changes in C:P, for example, during the spring bloom when phytoplankton C:P is expected to increase rapidly (Polimene et al., 2015).

The satellite-informed estimates of phytoplankton C:P and POC:POP presented here are still preliminary and, therefore, should not be treated as accurate estimates. Nevertheless, even with this simple two-end-member mixing model approach, we can make a testable hypothesis regarding the underlying mechanisms causing the observed temporal change in suspended POC:POP. First, to model temporal shifts in POC:POP, we need to consider the contribution that non-algal organic matter makes to POM and the change in phytoplankton C:P. Our results indicate that phytoplankton C:P alone leads to a considerable overestimation of bulk POC:POP, regionally and globally. Second, our satellite-informed bulk POC:POP can capture the seasonal trend in POC:POP, which shows elevated values during summer compared to winter. We are optimistic that with more sophisticated parameter calibration of the phytoplankton stoichiometry model and non-algal C:P, it will be possible to predict the temporal variability of POC:POP accurately in future studies.

Caveats, Limitations, and Future Needs

Satellite estimates of phytoplankton and bulk C:P have considerable uncertainty in the subtropical gyres during summer. This mainly stems from the fact that satellite-derived growth-rate estimates are considerably different depending on which NPP product is used. In the future, we also need to conduct careful sensitivity analyses of how different satellite-based algorithms of Chl, C_phyto, and POC as well their depth-dependent changes would affect satellite-informed estimates of ecosystem stoichiometry. For example, a recent study suggests a substantial shift in POC:POP from the surface to DCM in the Indian Ocean during summer (Sahoo et al., 2020). It is inherently challenging to characterize C:P accurately in subtropics with phytoplankton stoichiometry models (Garcia et al., 2020) as phytoplankton turnover happens quickly on a time scale of days (Malone et al., 1993). For a complete understanding of the temporal variability of phytoplankton and bulk C:P, measurements of phytoplankton-specific C:P using high-throughput flow cytometry (Graff et al., 2015; Kirchman, 2016) or single-probe mass spectrometry (Sun et al., 2018) would be necessary. Linking metagenomics data with the phytoplankton stoichiometry model and remote sensing may also help improve C:P estimates in the subtropics (Garcia et al., 2020).

In this study, we used parameters for Synechococcus, a cosmopolitan phytoplankton species with a broad biogeographic distribution that extends from tropics to subpolar regions (Flombaum et al., 2013; Berube et al., 2018). This parameterization should be representative of another picocyanobacterium, Prochlorococcus. Together, Prochlorococcus and Synechococcus are responsible for roughly a quarter of the total ocean net primary productivity (Flombaum et al., 2013). Given that the current satellite-derived products cannot easily resolve size-partitioned phytoplankton physiologies such as growth rate and Chl:C_phyto, it seems reasonable to tune the phytoplankton stoichiometry model to these most common phytoplankton types. With new advances in satellite instrumentation, such as the development of reliable hyperspectral ocean color measurements (Werdell et al., 2018; Schollaert Uz et al., 2019), we may be able to better resolve the size-specific C:P of different phytoplankton functional types. This would enable us to fully capture the spatio-temporal variability of community phytoplankton C:P, particularly in nutrient-rich upwelling and coastal regions where nano- and micro-phytoplankton are more dominant than picophytoplankton (Kostadinov et al., 2016).

We inferred the P limitation of phytoplankton by comparing satellite-based SST and the previously compiled mask of nutrient depletion temperature. Although our method can provide a first-order pattern of P limitation, this method cannot resolve the degree to which phytoplankton are P-stressed. In other words, we cannot determine whether phosphate is the primary or secondary limiting nutrient for phytoplankton growth (Moore et al., 2013). Also, a recent study suggests that we cannot determine for sure that phytoplankton is P-limited even when the observed phosphate concentration is below the detection limit (Martiny et al., 2019). Accurate determination of nutrient concentration from space is inherently challenging (Goes et al., 2000; Steinhoff et al., 2010; Arteaga et al., 2015), and this is one of the major bottlenecks for accurately probing phytoplankton nutrient limitation from space. Although there are no standard protocols or algorithms currently available, we may be able to accurately retrieve surface nutrient concentrations by using advanced statistical and machine learning techniques applied to satellite-retrieved sea surface salinity, temperature, and remote-sensing reflectance (e.g., Wang et al., 2018).

A fundamental assumption made when predicting bulk POC:POP is that C:P of non-algal organic matter is constant with a Redfield Ratio of 117. There is a consensus from previous marine and freshwater studies that C:P of heterotrophs is generally lower and more homeostatic (relatively constant) than that of phytoplankton (Elser and Urabe, 1999; Persson et al., 2010). The bulk POM, however, can be modified due to decomposition (Schneider et al., 2003), viral shunt (Jover et al., 2014), preferential remineralization (Shaffer et al., 1999), as well as the interplay between the dissolved and particulate pools. Measuring the elemental composition of separate constituents of organic matter should better help us constrain the most appropriate end-member C:P for non-algal organic matter. Alternatively, we can mechanistically predict C:P of bulk POM by coupling the phytoplankton stoichiometry model with models of prey-predator interaction and decomposition (e.g., Anderson et al., 2005; Butenschön et al., 2016; Tanioka and Matsumoto, 2018).

Conclusion

We showed that it is possible to determine spatially and temporally coherent patterns of the C:P ratios of a single class of marine phytoplankton and bulk POM using only remotely sensed information. In principle, the same method can be extended to the C:N ratio, which together with the C:P ratio can yield the N:P ratio. The results shown here should not be treated as accurate estimates of upper-ocean C:P but rather as a feasibility study that can benefit from more accurate remotely sensed estimates of growth rate and from a better understanding of the links between growth rate and stoichiometry for different size classes of marine phytoplankton. The number of measurements of the C:P ratio of individual POM components (i.e., algal and non-algal components) is currently insufficient to validate our estimates. However, our main conclusion highlighting the importance of community composition in controlling bulk POC:POP does not depend on the accuracy of stoichiometry estimates. This conclusion has important implications for estimating carbon and phosphorus fluxes to the deep ocean and for the trophic transfer to higher organisms. Indeed, if the POC:POP of exported POM is controlled by community composition rather than phytoplankton C:P, we would not expect large “stoichiometric buffering” of carbon export under climate-change scenarios as proposed by previous studies (Teng et al., 2014; Galbraith and Martiny, 2015; Tanioka and Matsumoto, 2017; Matsumoto et al., 2020a). However, the effects of change in phytoplankton C:P will become more critical for carbon export if the total % of phytoplankton in organic matter increases, or if the C:P of the non-algal component increases. We hope the questions raised here will foster collaborative work combining satellite remote sensing, field sampling, and numerical modeling specialists to improve our ability to predict organic-matter dynamics and reduce uncertainties in our projections of the future carbon cycle.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.5281/zenodo.4136353.

Author Contributions

TT, CF, and KM designed the study. TT gathered and analyzed data. All the authors wrote the manuscript.

Funding

This research was supported by US National Science Foundation (OCE-1827948, KM). The funding supports the project “A power law model of dynamic marine phytoplankton stoichiometry,” aimed at developing a new approach to representing stoichiometric diversity in ocean models and by quantifying the global impacts of that diversity under different climate conditions.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.604893/full#supplementary-material

Footnotes

1. ^ http://sites.science.oregonstate.edu/ocean.productivity/index.php (last access: June 22, 2020)
2. ^ http://oceancolor.gsfc.nasa.gov (last access: June 22, 2020)

References