Global declines in net primary production in the ocean color era

Domain: earth_science
Published: 2025-07-01
Editorial outcome: minor_revision

Abstract
The majority of heat associated with climate change has been absorbed in the sunlit surface ocean where phytoplankton carry out half of biospheric net primary production (NPP). The physical entrainment of nutrients from depth into the surface constrains NPP across most of the ocean, therefore it has been widely hypothesized that a warmer and more thermally stratified ocean will diminish NPP. As phytoplankton are the dominant driver of ocean color, the satellite remote sensing record is the best approach to assess global NPP trends. Here we show that statistically significant decreases in NPP have occurred in almost half of the ocean and these changes are dominated by declines in the tropical and subtropical stratified ocean. A deeper analysis confirms that strengthening nutrient limitation is largely driving declining NPP. Climate-mediated shifts in NPP represent a fundamental perturbation to biogeochemical cycles that can further weaken global fisheries.

Introduction
The ocean has absorbed more than 90% of the heat generated by anthropogenic climate change, resulting in a warmer and more thermally stratified environment across the vast majority of the ocean1,2. Most of this heat is confined to the upper few hundred meters of the ocean2, and it is in these environments where phytoplankton, a diverse assemblage of microscopic protists and cyanobacteria, carry out approximately half of biospheric photosynthesis3. Phytoplankton net primary production (NPP), the rate at which organic carbon is generated through photosynthesis, sets the foundation for marine food webs and fisheries4, is a key component of global biogeochemical cycles and a primary driver of biogenic carbon sequestration to the deep ocean5. The physical entrainment of vital nutrients from depth into the well-lit surface layer constrains NPP in most of the ocean6,7. Therefore, the prevailing conceptual view is that a warmer ocean with stronger thermal gradients will weaken nutrient entrainment and diminish NPP in much of the ocean8,9. Yet, how climate change is impacting NPP is still unclear. For example, multi-decadal direct measurements in the Atlantic and Pacific Oceans10,11 and earth system model projections show opposing trends12.

The key role of NPP in ecosystem services and the global carbon cycle has given rise to an ongoing pursuit of accurate, synoptic estimates at high resolution in both time and space13–15. Satellite remote sensing is the only platform that offers this capability, and past and present ocean color missions collectively provide an unbroken global time. Continual advances algorithms derived from remote sensing reflectance (Rrs), the primary metric of ocean color, have fundamentally improved the way we observe and estimate phytoplankton biomass, physiology, growth rates, and ultimately NPP from space16–18.

Here we present a 25-year global assessment across two satellite missions (SeaWiFS and MODIS-Aqua) to investigate how a warming ocean has influenced NPP in conjunction with underlying metrics of nutrient limitation. We use an absorption-based model that derives NPP as the product of photosynthetic active radiation (PAR), the fraction of PAR absorbed by phytoplankton (aϕ/a, where QPAR = PAR x aϕ/a), and the efficiency with which absorbed photosynthetic energy is converted into carbon biomass (ϕ μ)19. Direct field measurements20,21 and culture-based experiments22,23 have shown that ϕ μ has muted variance across taxonomic and nutrient gradients, consequently the vast majority of NPP variance is governed by ocean color measurements of QPAR. Physiological responses to increased nutrient limitation in phytoplankton are well understood and visible in the ocean color record24. As nutrient limitation increases, phytoplankton manufacture fewer light-harvesting pigments and associated photosynthetic machinery to match the lower energetic demands for growth25,26. This phenomenon is observable as concomitant decreases in the phytoplankton pigment absorption coefficient (aϕ) and phytoplankton growth rates (μ), where μ is derived by normalizing NPP to phytoplankton carbon biomass (Cphyto) as estimated from particulate backscattering27.

Monthly time series of NPP, growth rates, and all other ocean color parameters presented here (n=303, global 1/12° data) were deseasonalized using monthly climatologies (1997–2022). Spatially explicit trend analysis (Type 1 ordinary least-squares regression and uncertainty analysis accounting for autocorrelation, see Methods), was performed on deseasonalized data and is presented as absolute trends and percent change relative to the 25-year mean. Empirical orthogonal function (EOF) analysis was performed to identify dominant modes of variance in time and space and compared to a suite of climate indices. The specific NPP model variant chosen has no direct dependence on sea surface temperature (SST)19, allowing for an independent assessment of NPP trends and variance in a warming ocean.

Results
Trend analysis
Across the satellite record, NPP is decreasing throughout vast regions of the subtropical and tropical ocean and increasing in polar and subpolar waters (Fig. 1a). Overall, NPP has significantly declined in 48% of the ocean at an average rate of −34 mg C m−2 d−1 decade−1, equivalent to a –7.1% normalized decadal rate. Conversely only 4% of the ocean has experienced a significant increase in NPP but at a higher average rate of 48 mg C m−2 d−1 decade−1, this is equivalent to a +20% normalized decadal rate reflecting the fact that significant increases in NPP are largely occurring in low productivity regions (Fig. 1b, c). Spatially contiguous trends broadly correspond to major ocean biomes and currents (Supplementary Fig. 1). NPP is declining in all subtropical gyres and along currents that advect warm water poleward (Brazil, Gulf Stream, Kuroshio, Norwegian). Conversely, NPP is increasing along the austral subtropical front and either increasing or unchanged along eastern boundary currents that advect colder water towards the equator (California, Greenland, Malvinas, Oyashio, and Peru). A clear latitudinal anomaly is increasing NPP centered around 15°S in the southeast Pacific Ocean (Peru current). This region, along with the North Atlantic where NPP is also increasing, is one of the few areas globally where SST has cooled across the ocean color record1. The strong correspondence of SST and ΔNPP is further shown by grouping trends into respective 1°C mean SST (1998–2022) bins and computing the area within each bin where NPP trends are significant (Fig. 1d, e). These data show that the 15°C isotherm, which has been used to delineate temperate seasonally mixing biomes from warmer permanently stratified biomes24, also represents the transition from increasing to decreasing ΔNPP (Fig. 1d). Above 15°C most declining trends exceed the 90% confidence interval, whereas below this temperature most of the increasing trends are not significant, and the largest increases in NPP are occurring in polar waters (Fig. 1e).

Trend analyses of the atmospheric (PAR), bio-optical (aϕ/a), and physiological (ϕ μ) components of NPP demonstrates that changes in the fraction of energy absorbed by phytoplankton (aϕ/a) is largely driving NPP trends (Fig. 2). Overall, aϕ/a has significantly decreased in 43% of the ocean at an average normalized decadal rate of −7.6% and increased in 6% of the ocean at an average normalized decadal rate of +18%. Significant trends in PAR and ϕ μ are present in only 16 and 22% of the ocean, respectively, and are comparatively muted, with average decadal rates less than 4% (both positive and negative for each parameter). There is strong spatial overlap in regions where both NPP and aϕ/a have significant trends. 86% (81%) of regions where NPP is significantly decreasing (increasing) has seen commensurate significant trends in aϕ/a. This analysis reveals that the dominant driver of NPP is embedded within the ocean color record, and the computationally simple term aϕ/a, derived directly from satellite data, is a robust indicator of shifts in marine productivity.

Declining NPP and aϕ/a across most of the stratified ocean is consistent with stronger and more pervasive nutrient limitation. However, these declines could also occur through nutrient-independent increases in the spectral absorption coefficient of dissolved organic matter (adg) or decreases in phytoplankton biomass (Cphyto) without commensurate changes in growth rates (μ). Therefore to better resolve the mechanisms causing changes in NPP embedded in the ocean color record, normalized trends, slopes divided by their respective mean, of aϕ and adg, as well as euphotic zone integrated Cphyto and μ were derived and directly compared. Significant decreases (increases) in aϕ are occurring in 36% (6%) of the ocean (Fig. 3), significant decreases (increases) in adg are occurring in 25% (3%) of the ocean, a direct comparison of normalized decadal trends confirm that declining aϕ is the globally dominant trend in 33% of the ocean (Supplementary Fig. 2). Significant decreases (increases) in μ are occurring in 44% (4%) of the ocean and euphotic zone integrated Cphyto also registers statistically significant decreases (increases) in 17% (6%) of the ocean. Moreover, the normalized decadal Cphyto rates (−2% and +4%) are comparatively muted relatively to the normalized decadal μ rates (−11%, +23%, Fig. 3). In this framework, decreases in aϕ, μ and NPP strongly suggest that the warming and more thermally stratified surface ocean is an increasingly nutrient limited environment, and that cellular reductions in photosynthetic pigments rather than changes in biomass is the dominant response observable in the ocean color record. This observation is mathematically consistent with statistically significant changes in Rrs that have emerged in more than half of the global ocean28.

NPP anomaly variability in space and time
Previous analyses of global NPP anomalies (1997–2008) have shown that interannual variability is on the order of ±2 Petagram (Pg) C year−1 and is strongly influenced by SST (r=0.80), particularly in the permanently stratified surface ocean29. Past analyses should be interpreted with some caution, however, as they often employ NPP models with first-order dependencies on SST30 and, unlike the model used here, their behavior cannot be evaluated independently from changes in SST. Across the 25-year record, NPP averaged 56.8 Pg C year−1 with a standard deviation of 1.9 Pg C year−1 (Fig. 4a). Global NPP declined at an average rate of 1.8 Pg C decade−1, though interannual variability is pervasive. This change represents a net 7.4% decrease over the satellite record, which greatly exceeds the CMIP6 high emissions scenario models (RCP8.5) that project NPP to decrease by 2.99% by the end of the 21st century12.

To better elucidate the drivers and global patterns of NPP anomalies, we conducted EOF analysis to identify the dominant modes of NPP anomaly variance in space and time and then compared these modes to a suite of oceanic and atmospheric climate indices, including monthly global SST anomalies over the same time period (Fig. 4b). The first mode of global NPP anomaly variance (PC1) explains 13.3% of total variance, encapsulates the largest magnitude of NPP variance with a total range of 2.8 Pg C, and corresponds to global declines in NPP that become increasingly pronounced towards the equator. Global SST anomalies were the highest correlated climate index (r=-0.63), suggesting that a warming ocean is now the dominant driver of interannual variability in marine photosynthesis. The second mode of variance explains 5.3% of variability, is weakly positively correlated with SST anomalies (r=0.23) and encapsulates a comparatively smaller range of NPP (0.7 Pg C). This second mode has a strong latitudinal signal, with modest negative and positive NPP anomalies in the subtropical northern and southern hemispheres, respectively. The most striking feature of PC2 is the pseudo-seasonal time component that emerges from deseasonalized data. A more detailed inspection of this mode shows that regions of high PC2 variance occur near the 15°C isotherm in both hemispheres and these regions have undergone pronounced phenological shifts in NPP (Supplementary Fig. 3). In regions of high PC2 variance, NPP is increasing in the winter and decreasing in the summer in both hemispheres, resulting in diminished seasonality (Supplementary Fig. 3) that gives rise to pseudo-seasonal anomaly variance. That shifting phenology explains more of the global NPP anomaly variance than other established climate indices (e.g., El Niño) is notable. The third mode of variance explains 3.9% of NPP anomaly variability, encapsulates a total range of 0.9 Pg C, and is highly negatively correlated with the Pacific Decadal Oscillation (PDO, r=−0.63). The transition from a positive to strongly negative PDO in part explains the reversal in declining global rates from 2018 to 2021. Finally, the fourth mode of variance explains 2.8% of variability, encapsulates a total range of 0.5 Pg C, and is positively correlated with the El Niño 3.4 index (r=0.61).

Deepening of NPP
Satellite ocean color measurements are restricted to the upper few meters of the surface mixed layer. In most of the ocean, the euphotic zone (depths where NPP > 0) extends beyond the surface mixed layer, and phytoplankton residing at depths beyond the reach of ocean color sensors often form so-called deep chlorophyll maxima (DCM). A key feature of DCMs is that they often represent a photoacclimative response to low-light, whereby cellular pigmentation increases with no change in biomass31. The NPP model employed here includes a photoacclimative response beneath the surface mixed layer that performs well compared to direct measurements19 and, when applied globally, constitutes more than 25% of global annual NPP (Fig. 5). Declining absorption in the surface deepens euphotic depths, thus permitting greater NPP at depth (Supplementary Fig. 4). The fraction of NPP within the DCM is sensitive to the surface mixed layer depth (MLD) which is regulated by wind stress, buoyancy input, and other physical processes. The MLD is deepening in some regions of the ocean and becoming shallower in other regions (Supplementary Fig. 4), although in general, the density contrast across the base of the surface mixed layer is strengthening32. On aggregate, based on a photoacclimative response alone, NPP beneath the surface mixed layer has increased by almost 3 Pg C across the satellite record and now constitutes more than 30% of global NPP (Supplementary Fig. 4).

Methods
Data sources
Absorption-based NPP is calculated using the CAFE model19. Level 3, monthly, 1/12° ocean color data (2022 reprocessing version) were downloaded from the NASA ocean color webpage for the SeaWiFS (09/1997–12/2007) and MODIS Aqua (07/2002–12/2022) missions. Inherent optical properties are determined via inversion of remote-sensing reflectance using the generalized inherent optical property model18. The CAFE model is spectrally explicit and evaluated at 10 nm increments in the visible spectrum (400–700 nm). Monthly MLD were retrieved from the Oregon state ocean productivity website (https://sites.science.oregonstate.edu/ocean.productivity/). The MLD is calculated as the depth where the density of water is 0.125 kg m−3 greater than the density of near-surface water (z=10 meters), where density is derived from HYCOM salinity and temperature data. Monthly 1/12° SST data are from the optimum interpolation sea surface temperature advanced very high resolution radiometer version 2.145.

Satellite blending
To account for persistent biases between SeaWiFS and MODIS Aqua missions, we compared NPP from collocated monthly 1/12° data for the period 2002–2007. Overall, NPP showed remarkable correspondence between satellite missions (Supplementary Fig. 5); Type-1 linear regression returned an r value and slope of 0.96 and 0.97, respectively, and a mean absolute error of 1.09. That said, MODIS Aqua generally had higher NPP (Supplementary Fig. 5), therefore to accommodate for persistent spatial biases we followed the general convention of inter-satellite comparisons17,46: SeaWiFS data were adjusted (SeaWiFS*) by the monthly climatology differences derived from overlapping periods between the two missions. This adjustment resulted in an r value and slope of 0.98 and 0.99, respectively, and mean absolute error of 1.04. Satellite data (e.g., NPP, IOPs) were averaged between SeaWiFS* and MODIS-Aqua for the 2002–2007 overlap.

Statistical and spatial analysis
For all time series, monthly anomalies were first calculated by subtracting monthly data from monthly climatology (n=303). All statistics were performed using the Python package SciPy47. Spatially explicit trends and residuals were computed using ordinary least-squares regression. Serial autocorrelation in anomaly residuals was quantified using a first-order model (i.e., monthly) that reduced the effective sample size of the time series. The standard error of the slope was recomputed using the reduced effective sample size. Uncertainties are expressed as 90% two-tailed confidence intervals derived from the effective sample size and recomputed standard error of the slope. EOF were derived after downsampling the data to an equal area 50×50 km grid. The full list of climate and oceanographic indices and their sources that were compared to the principal components are given in Supplementary Table 1. Gap filling of data using monthly climatological data was only performed to compute global NPP totals (Fig. 4a). On average less than 2% of ocean surface area where PAR > 1 mol photons m−2 day−1 required gap filling. All spatial analyses account for areal changes in pixel size and apply a global land mask created for satellite ocean color data48.

Data availability
Monthly global CAFE model output is available at https://doi.org/10.5281/zenodo.12586956. Remotely sensed data are available at https://oceancolor.gsfc.nasa.gov/l3/. Monthly global MLD data are available at http://orca.science.oregonstate.edu/2160.by.4320.monthly.hdf.mld030.hycom.php. Data generated in Figs. 1 through 5 are available at https://doi.org/10.5281/zenodo.15497141. All other data that support the findings of this study are available from the corresponding author upon reasonable request.

Code availability
All major codes are available at doi.org/10.5281/zenodo.15497980. This repository includes Phyton and C code to compute NPP and associated data, and Python scripts for the trend and EOF analysis presented here.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Using long-term spatiotemporal time series data of remote sensing ocean color, this study aims to show a long-term decline of NPP in the ocean, with clarification of spatial variation. The first reading gives the impression that this manuscript is very similar to Behrenfeld, et al. (1996), albeit with updated data. However, careful investigations of their analyses and results lead me to conclude that this work is interesting, particularly in providing a better mechanistic understanding of the global decline of NPP. My main concerns are about the presentation: 1. The novelty of this work is unclear due to the way of writing; as such, the manuscript reads like an update of Behrenfeld, et al. (1996) for none- specialists of ocean color remote sensing; 2. The methods are unclear. Below, I list some suggestions. My comments aim to make this manuscript accessible to a broad audience of Nat. Comm.\nMain concerns:\n1. The Abstract should include a better mechanistic explanation in addition to stating only the observed patterns. Otherwise, the findings seem similar to Behrenfeld et al. (1996), except that the time series is longer.\n2. I feel the Introduction needs substantial re-writing. The current version reads too much similar to Behrenfeld et al. (1996). The critical \"knowledge gap\" (the key innovation of this work) only starts to be unveiled from Line 61. Unfortunately, Lines 61 to 66 are too technical. Thus, the current writing does not successfully bring out the novelty of this work. I suggest some sentences in Lines 130-140 could be moved to the Introduction with some modifications. Importantly, readers expect an improved mechanistic understanding of the research topic.\nI, by no means, suggest the authors criticize Behrenfeld et al. (1996). Nevertheless, the authors should make the sentences clear about what achievement is beyond Behrenfeld et al. (1996) in this work.\n3. The sentences in the Discussion read like implications of findings instead of a real discussion of the results.\n4. What I feel is missing in the Discussion is more detailed physiological components of remote sensing modeling. See a nice review paper: Westberry et al. (2023), actually by the same group of authors of this manuscript. I wish to see the authors comment on how the remote sensing modeling of NPP can be moved forward, for better mechanistic understanding.\n5. Fig. 2a. shows a very interesting reverse pattern from 2018 to 2021. Please discuss this.\n6. Extended Data Figures 3 and 4 in SI are more interesting than Figures 2 and 3 in the main text because Extended Data Figures 3 and 4 provide mechanistic explanations of decreasing NPP.\nThe spatial pattern of NPP (PC 3-4) is, to a large extent, affected by changes in circulation and SST that are driven by climate. Thus, the correlation between the PC3-4 and the climate index is not surprising. I feel detailed explanations of Figure 2 are not the most novel findings of this work. Instead, I wish to see more detailed explanations and discussion for Extended Data Figures 3 and 4, which provide “mechanistic understanding from the viewpoint of pigment absorption.”\n7. Monthly anomalies were calculated by subtracting monthly data from monthly climatology. According to Fig 2b, PC2, the seasonality remains very visible. This issue is probably not critical to the key findings of this work. However, please discuss this, to avoid confusion.\n8. Methods: How to calculate and obtain Extend Data Figures 2 to 5 need detailed explanations in the Method section. Codes to reproduce the results should be provided along with the paper. The majority of readers are not experts in remote sensing.\n9. In many places, the wording is too technical, as indicated in the editorial suggestions below. Please clarify the biological meanings of those. I just made some examples; the authors should carefully revise the manuscript.\nEditorial suggestions:\nLine 28: Please clarify \"phytoplankton pigmentation\"?\nLine 33: Please explain the biological meaning of \"declining pigment absorption\". I do not think the broad audience of Nat. Comm. understand this term.\nFinally\nI look forward to using this as teaching material in my Biological Oceanography class.\nChih-hao Hsieh", "reviewer_2": "Review of Silsbe et al. “Global declines in net phytoplankton production in the ocean color era”\nThis manuscript presents an analysis of a merged 25 year dataset of net primary production derived from satellite data. The results suggest a decline in NPP in low latitude regions, and furthermore explores the photo-physiological drivers underlying the trends. The manuscript presents some interesting and relevant results, but I would urge the authors to assess the uncertainties introduced by their choice of satellite time series data and how the merged time series is created.\nMethods: It would be wise to check the merged SeaWiFS-MODIS-Aqua time series the authors have created for any discontinuities introduced when different combinations of satellites are included in the time series. Change point detection is one option. Existing packages are available in matlab as findchangepts or R as changepoint (and I’m sure Python has something equivalent). If the algorithms pick out dates associated with when an additional data source is introduced/removed, then artificial trends could be created in the time series.\nI also wonder why the authors didn’t use the OC-CCI dataset, as it incorporates multiple ocean colour satellites (not just the 2 used here) and has made a huge effort in producing a consistent, climate record-quality timeseries. That would overcome my misgivings about using the 2-satellite merged time series here. There has also been reported deterioration in the MODIS-Aqua dataset for the last couple of years due to sensor degradation. Has this been taken into account here?\nGiven the lack of data at high latitudes in winter, due to both cloud/ice cover and low sun angle, the gap-filling routine applied here must be effectively “making up” a significant fraction of the time series at high latitudes. Some discussion of the implications of this for your results needs to be added.\nLines 26-27: the authors use “net phytoplankton production (NPP)” here rather than the more common “net primary production” (also commonly abbreviated as NPP). Is that intentional?\nLine 104: “a 7.4% over” – this is incomplete.\nLines 115-119: the second mode of variance appears to represent the opposed seasonal cycles of the Southern and Northern Hemispheres, and indeed I think I counted 1 peak/trough per year of the time series too. The authors hoped to remove the seasonal cycle by taking anomalies from the climatological year, but there seems to be a large amount of residual seasonality in the time series. The authors should consider whether there is a better way to remove the seasonal cycle so it isn’t so dominant in the EOF analysis. Could they try applying an EOF to the full time series (before seasonality is removed), and the 1st mode will probably be the seasonal signal, then subtract that from the full time series, and run the EOF again??\nLines 120-124: the aseasonal pattern of the Southern Ocean is emphasised here, but the EOF spatial pattern is as strong in the Northern Hemisphere as it is in the SOcean. The discussion needs to be more balanced here.\nLine 135: define phi_mu again here to help the reader\nLine 145-147: is the result reported in this sentence shown in this manuscript (if so, direct the reader to a figure) or is it from literature and so needs a reference?\nLine 171: how is the SML defined in this study?\nLines 188: specify that the NPP models used here are satellite data-based\nLine 206: I find use of the word ‘evolution’ here potentially misleading as the authors aren’t really talking about an evolutionary process (I think? Given the 25 year time scale of the study)\nFigure 1a: mark the non-significant pixels on the map with hatching (or something)\nFigure 1c: it’s very hard to distinguish between non-significant and small trends with this colour scheme. As for figure 1a, add hatching to identify non-significant pixels.\nFigure 1d: specify that the inset is the SST and what timeframe it covers\nFigure 1e: I didn’t really understand this figure! What are the ‘top’ and ‘bottom’ sections of the plot. Why is the shaded area both positive and negative?\nLines 402-403: please acknowledge here the lack of data in winter at high latitudes, and that a lot of gap-filling will need to have taken place to make a complete timeseries.\nExtended Figure 2c: I didn’t follow what “spatially averaged latitude of austral waters” meant\nExtended Figure 3d-f: are these panels needed? They’re not very informative.\nExtended Figure 3g-i: add to caption that the legend also includes the proportion of the ocean in each category."}, "author_response": "The abstract has been revised with a better mechanistic explanation (Lines 31:34, please note that Lines refer to the manuscript where track changes are hidden).\nThe introduction has been substantially revised. The critical knowledge gap is now more explicitly stated with less technical jargon (Line 60-74). Following your suggestion, lines from the Results have been incorporated into this new text. The critical knowledge gap still begins in the third paragraph, as we feel it is still important to introduce NPP (Paragraph 1) and remote sensing (Paragraph 2) to a general audience before arriving at this text.\nThe Discussion has been revised (see next comment).\nWe have added several sentences summarizing some of the key concepts how remote sensing modeling can be moved forward, as also discussed in Westberry et al. (2023) (Lines 225:238).\nThank you for this comment. We cannot fully explain the complete reversal pattern, however some of this reversal is due to the transition from a strongly positive to strongly negative PDO (Lines 173-174).\nConsistent with changes to the abstract and introduction, we have brought forth the most pertinent data from the Extended Data Figures 3 and 4 to form a new figure (Figure 3) in the main text and changed the text accordingly, with commensurate changes in Extended Data Figures.\nThank you for this comment, given a similar comment from Reviewer #2 we have strived to provide a clearer and more balanced analysis of PC2. Here we briefly summarize this change, please see Lines 161-169 and the newer Extended Data Figure (now Extended Data Figure 4) that provides a clearer explanation of PC2. While the PC2 time-series appears seasonal, it in fact represents regions of the ocean where summer NPP is decreasing but winter productivity is increasing, leading to an overall shift towards more muted seasonality. This phenological shift is simultaneously occurring in both hemispheres, such that declining summer NPP in one hemisphere is offset by increases in winter NPP in the other hemisphere and vice-versa, a phenomenon that gives rise to the pseudo-seasonal time series and strong latitudinal gradients shown in Figure 4. That this phenomenon constitutes 5.3% of the global NPP anomaly variance is notable as it exceeds the variance imparted by both the PDO and ENSO.\nTrend analyses are identical for NPP and all other satellite data. We have attempted to make this more clear in the Results (Lines 125-127) as well as the Methods. We have created a GitHub repository with all relevant code (Line 444-446).\nThis has been changed simply to ‘light-harvesting pigment capacity’ in our revised abstract (Line 33). We hope the inclusion of light-harvesting capacity is enough for a general audience to understand the physiological role of pigments.\nSee previous comment.\nThank you! We strive for effective science visualization.\nWe conducted a change point analysis of global NPP anomalies: Autoregression in the anomalies was removed as described in the methods, then multiple change point detection was determined using the Python package ‘ruptures’ (C. Truong et al. 2020) parameterized with a least square deviation cost function and a minimum length between change points of 3 years (results were robust across different cost functions and lengths). This analysis identified 3 change points (Jun-2004, Jul-2011, Jun-2019) that are different than the introduction of MODIS Aqua (Aug-2002) and our discontinuation of SeaWiFS (Dec-2008).\nWe considered the OC-CCI dataset but decided upon exclusively SeaWiFS and Aqua for four reasons. 1: The OC-CCI does not have a PAR product that is essential to all NPP models. 2: The MODIS Aqua record, that dominates the time series (20 years), has greater spectral resolution than OC-CCI. 3: The spectral inversion algorithm employed by the OC-CCI (Lee et al. 2009) has higher uncertainty (Brewin et al. 2015, http://dx.doi.org/10.1016/j.rse.2013.09.016) than the NASA operational GIOP spectral inversion model (Werdell et al. 2013) employed here. 4: At time of submission, the OC-CCI was using an older reprocessing version (R2018.0) of NASA ocean color data. More to the point, our methods essentially follow the bias correction method adopted by the OC-CCI. We thank you for the comment and to allay future misgivings, we have made this more clear in our Methods (Line 413-414, please note that Lines refer to the manuscript where track changes are hidden).\nAs we were aware of the deterioration on MODIS-Aqua, we purposefully ended our time series in 2022 to minimize instrument artefacts and to align with the most recent MODIS Aqua reprocessing at time of submission (R2022.0, https://oceancolor.gsfc.nasa.gov/data/reprocessing/r2022/).\nWe have clarified this in our methods to more clearly state that gap filling was only performed to compute the global NPP Trend (Figure 4A), and that the missing area constituted less than 2% of the ocean surface where PAR > 1 mol photons m-2 day-1. In other words once you factor in areas subjected to polar night when NPP is negligible, gap filling is in fact minor. (Lines 429-433).\nWe have changed NPP to net primary production, and when introducing NPP the first time we state ‘phytoplankton Net Primary Production (NPP)’. We originally used the term phytoplankton, and add the distinction above, as non-phytoplankton autotrophs (i.e. kelps, sea grasses) contribute a very small fraction of marine primary production.\nThe sentence now states ‘a 7.4% reduction over the satellite record’ Line 147-148.\nThank you for this comment, given a similar comment from Reviewer #12 we have strived to provide a clearer and more balanced analysis of PC2. Here we briefly summarize this change, please see Lines 161-167 and the newer Extended Data Figure (now Extended Data Figure 4) that provides a clearer explanation of PC2. While the PC2 time-series appears seasonal, it in fact represents regions of the ocean where summer NPP is decreasing but winter productivity is increasing, leading to an overall shift towards more muted seasonality. This phenological shift is simultaneously occurring in both hemispheres, such that declining summer NPP in one hemisphere is offset by increases in winter NPP in the other hemisphere and vice-versa, a phenomenon that gives rise to the pseudo-seasonal time series and strong latitudinal gradients shown in Figure 4. That this phenomenon constitutes 5.3% of the global NPP anomaly variance is notable as it exceeds the variance imparted by both the PDO and ENSO.\nFollowing the suggestion of reviewer 1, this sentence has been moved to the introduction and immediately follows phi_mu being defined (Line 65-66).\nOriginally, this sentence was missing a reference to an extended data figure. Following comments by Reviewer 1, this figure is now in the main section of the manuscript and the text has been modified. (Lines 123-129).\nThe definition of the surface mixed layer is provided in the Methods (Lines 405-406).\nThe line now reads ‘and satellite-based NPP models’ (Line 192).\n‘oxygen evolution’ has been changed to ‘oxygen production’ (Line 213).\nHatching has been added to all relevant figures.\nHatching has been added to all relevant figures.\nWe have added an inset label.\nIn the original version, the y-axis was mirrored to show how much of the ocean, binned by average SST, had significant increasing NPP (top half), or decreasing NPP (bottom half). We have simplified this figure to show areas of significant NPP trends regardless of direction and changed the figure caption accordingly.\nThis section has changed following both reviewers suggestion, and the sentence no longer appears.\nThis figure has been changed, and following comments from Reviewer #1, is now in Figure 3.\nExtended Data Figure 2 Caption now reads ‘Frequency map showing the largest statistically significant normalized trend (increasing or decreasing) for inherent optical properties.’"}, {"round": 2, "reviewer_comments": {"reviewer_1": "My comments have been well addressed. The revision is satisfactory.", "reviewer_2": "I find the manuscript improved in terms of clarity for the reader and outlining the implications of the analysis. I just have a couple of minor corrections:\n- Title still includes \"net phytoplankton production\" instead of \"net primary production\", although it has been changed in the rest of the manuscript.\n- Line 186 should refer to Figure 5 (currently says Figure 4).\n- Line 196: anthropogenic heat is not predominantly being absorbed by ecosystems, but rather by the ocean water itself. Of course there will be effects on the ecosystem, but it is wrong to say that \"marine ecosystems are the predominant sink of anthropogenic heat\"."}, "author_response": "Reviewer 1: No revisions given.\nR2: I just have a couple of minor corrections:\n- Title still includes \"net phytoplankton production\" instead of \"net primary production\", although it has been changed in the rest of the manuscript.\nTitle Changed.\n- Line 186 should refer to Figure 5 (currently says Figure 4).\nCorrection made (now Line 187).\n- Line 196: anthropogenic heat is not predominantly being absorbed by ecosystems, but rather by the ocean water itself. Of course there will be effects on the ecosystem, but it is wrong to say that \"marine ecosystems are the predominant sink of anthropogenic heat\".\n‘Marine Ecosystems’ replace with ‘Oceanic waters’ (now Line 197)."}]

Decision letter