Main

Generating roughly half the planetary primary production1, marine phytoplankton affect the abundance and diversity of marine organisms, drive marine ecosystem functioning, and set the upper limits to fishery yields2. Phytoplankton strongly influence climate processes3 and biogeochemical cycles4,5, particularly the carbon cycle. Despite this far-reaching importance, empirical estimates of long-term trends in phytoplankton abundance remain limited.

Estimated changes in marine phytoplankton using satellite remote sensing (1979–86 and 1997–present) have been variable6, with reported global decreases7 and increases8,9, and large interannual10 and decadal-scale variability11. Despite differences in scale and approach, it is clear that long-term estimates of phytoplankton abundance are a necessary, but elusive, prerequisite to understanding macroecological changes in the ocean10,11,12,13.

Phytoplankton biomass is commonly inferred from measures of total chlorophyll pigment concentration (‘Chl’). As Chl explains much of the variance in marine primary production14 and captures first-order changes in phytoplankton biomass, it is considered a reliable indicator of both phytoplankton production and biomass15. Shipboard measurements of upper ocean Chl have been made since the early 1900s, first using spectrophotometric and then fluorometric analyses of filtered seawater residues, and more recently through in vivo measurements of phytoplankton fluorescence16. Additionally, measurements of upper ocean transparency using the standardized Secchi disk are available from 1899 to present and can be related to surface Chl through empirically based optical equations17,18. Although the Secchi disk is one of the oldest and simplest oceanographic instruments, Chl concentrations derived from Secchi depth observations are closely comparable to those estimated from direct in situ optical measurements or satellite remote sensing18.

We compiled publicly available in situ Chl and ocean transparency measurements collected in the upper ocean over the past century (Fig. 1a–c; see Supplementary Information for data sources). Transparency measurements were converted to depth-averaged Chl concentrations using established models17. Systematic filtration algorithms were applied to remove erroneous and biologically unrealistic Chl measurements, and to exclude those in waters <25 m deep or <1 km from the coast, where terrigenous and re-suspended substances introduce optical errors. In situ and transparency derived Chl measurements (monthly averages for each year, 0.25° resolution) were strongly correlated (r = 0.52; P < 0.0001). After log-transforming these data to achieve normality and homoscedasticity, model II major axis regression analysis revealed linear scaling of transparency-derived and in situ-derived Chl (intercept, 0.18; slope, 1.08 ± 0.016; r2 = 0.60). Both this and additional analyses indicated that both data sources were statistically similar enough to combine (see Methods and Supplementary Figs 2, 3). The blended data consisted of 445,237 globally distributed Chl measurements collected between 1899 and 2008 (Fig. 1a). Data density was greatest in the North Atlantic and Pacific oceans and after 1930 (Fig. 1b, c), and broadly reproduced spatial patterns of phytoplankton biomass derived from remote sensing7 (Fig. 1d and Supplementary Fig. 3).

Figure 1: Data availability.
figure 1

a, Temporal availability of ocean transparency (red), and in situ Chl (blue) measurements. Bars represent the proportion of total observations collected in each year, coloured ticks on x axes represent years containing data. b, c, Spatial distribution of in situ Chl (b) and transparency data (c). Colours depict the number of measurements per 5° × 5° cell (ln-transformed). d, Averaged Chl concentration from blended transparency and in situ data per cell.

PowerPoint slide

Full size image

Chl trends were estimated using generalized additive models (GAMs)19. These models are extensions of generalized linear models that do not require prior knowledge of the shape of the response function. To ensure robustness, Chl trends were estimated at three different spatial scales—local, regional and global.

Local-scale phytoplankton trends

To estimate local Chl trends, blended data were binned onto a 10° × 10° global grid and GAMs of Chl as log-linear functions of covariates were fitted to data within each cell. Phytoplankton declines were observed in 59% (n = 214) of the cells containing sufficient data (Fig. 2a, b). Clusters of increasing cells were found across the eastern Pacific, and the northern and eastern Indian Ocean (Fig. 2b). High-latitude areas (>60°) showed the greatest proportion of declining cells (range: 78–80%).

Figure 2: Local-scale trends in phytoplankton.
figure 2

a, Baseline year and temporal span of Chl data used in local models. b, Mean instantaneous rates of Chl change in each 10° × 10° cell (n = 364). Yellow and red represent cells where Chl has increased, while blue represents a Chl decrease. Cells bordered in black denote statistically significant rates of change (P < 0.05) and white cells indicate insufficient data. c, Mean instantaneous rates of Chl change for each 10° × 10° cell, estimated as a function of distance from the nearest coastline (km) and baseline year of trend. Colour shading depicts the magnitude of change per year. All effects used to fit the trend surface were statistically significant (P < 0.05).

PowerPoint slide

Full size image

Owing to sparse observations in early years, local trends were also estimated using post-1950 data only. This yielded almost identical results, although the magnitude of change was amplified in some cells (see Supplementary Fig. 7).

Local models further suggested that Chl has declined more rapidly with increasing distance from land (Fig. 2c). This agrees with results derived from satellite data, documenting declining phytoplankton in the open oceans8,20,21, and expansion of oligotrophic gyres, probably due to intensifying vertical stratification and ocean warming10,22. These trends are noteworthy, because most (75%) aquatic primary production occurs in these waters23. In shelf regions, Chl trends switched from negative to positive in more recent years (since ∼1980), consistent with reported Chl increases due to intensifying coastal eutrophication and land runoff8.

Regional and global phytoplankton trends

To estimate regional Chl trends, we divided the global ocean into ten regions, in which similar variability in phytoplankton biomass was observed in response to seasonality and climate forcing24 (Fig. 3a). To capture the range of potential Chl trajectories, regional trends were estimated from GAMs as linear functions of time on a log scale in three different ways: (1) continuous (linear trend), (2) discrete (mean year-by-year estimates), and (3) smooth functions of time (non-monotonic trend). This approach is comprehensive; it allows both the quantitative (magnitude) and qualitative nature (shape) of trends to be estimated (see Methods Summary and Supplementary Information for full details).

Figure 3: Regional and global trends in phytoplankton.
figure 3

a, Ocean regions (n = 10) used to estimate regional trends in Chl. N., North; Eq., Equatorial; S., South. b, Mean instantaneous rates of Chl change estimated for each region, with 95% confidence limits. Diamonds indicate the global meta-analytic mean rate of Chl change, with 95% confidence intervals. Trends were estimated using all available data (red symbols) and data since 1950 only (blue symbols). Individual estimates are displayed as tick-marks on the x axis. All estimates were statistically significant (P < 0.05), except for the North Indian region (P = 0.27).

PowerPoint slide

Full size image

Estimation of Chl trends as continuous log-linear functions of time revealed phytoplankton declines in eight out of the ten regions. The largest rates of decline were observed in the South Atlantic (−0.018 ± 0.0015 mg m−3 yr−1), Southern (−0.015 ± 0.0016 mg m−3 yr−1), and Equatorial Atlantic (−0.013 ± 0.0012 mg m−3 yr−1) regions (P < 0.0001 for all trends; Fig. 3b). Increases were observed in the North Indian (0.0018 ± 0.0015 mg m−3 yr−1; P = 0.268) and South Indian regions (0.02 ± 0.0011 mg m−3 yr−1; P < 0.0001). The global meta-analytic mean rate of Chl change derived from individual regional model estimates was −0.006 ± 0.0017 mg m−3 yr−1 (P < 0.0001; Fig. 3b), representing an annual rate of decline of ∼1% relative to the global median chlorophyll concentration (∼0.56 mg m−3).

Regional trends were also estimated using data since 1950 only, but the direction of all trends remained unchanged and the magnitude of changes was minimal (Fig. 3b). Post-1950 trends were amplified in some regions, resulting in a greater but more variable global rate of decline (−0.008 ± 0.0068 mg m−3 yr−1; P < 0.0001). Estimating regional trends separately for each data source yielded similar results (see Supplementary Fig. 4).

Modelling Chl trends as both discrete and smooth functions of time revealed pronounced interannual to decadal fluctuations superimposed on long-term trends (Fig. 4a). We observed greater Chl fluctuations in the Southern Hemisphere regions and greater uncertainty about estimates before 1950; both issues probably reflect limitations in data availability for these regions and time periods. In the polar and Atlantic regions, Chl increased until ∼1950, before undergoing prolonged declines (about 1950–95). After ∼1995, sharp increases were observed in the South Indian and Southern regions (Fig. 4a).

Figure 4: Temporal variability in phytoplankton trends.
figure 4

a, GAM estimates of Chl as a discrete (points) or smooth function (lines) of yearly variability in each region (n = 10). For each initial year, Chl is the arithmetic, rather than model-estimated, mean. Temporal data availability is displayed as tick-marks on the x axis. b, Seasonal patterns of Chl as a smooth function of day of year in each northern (blue), equatorial (red) and southern (green) ocean. Shaded areas represent approximate 95% Bayesian credible limits around each estimate.

PowerPoint slide

Full size image

GAMs also accounted for mean seasonal variation in Chl (Fig. 4b) and closely reproduced known patterns24,25, providing a measure of confidence in our approach. Strong seasonality in polar regions reflects pronounced variability in mixing, irradiance and ice cover26 whereas weak seasonality in equatorial regions is a function of near-constant solar irradiance. Complex seasonality in the Indian Ocean relates to the effects of monsoon dynamics and freshwater inputs on nutrient delivery27. Temperate regions are affected by seasonally changing solar irradiance and trade winds, and their effects on upper ocean nutrient delivery28.

Climate effects on phytoplankton

Regional phytoplankton trends display both short-term variation and longer-term trends. We tested the hypothesis that the short-term (interannual to decadal) component in Chl variation may be explained by the effects of leading climate oscillators, such as the El Niño/Southern Oscillation (ENSO) or the North Atlantic Oscillation (NAO). After de-trending and removing seasonal variation, yearly Chl anomalies were strongly negatively correlated with the bivariate ENSO index in the Equatorial Pacific (r = −0.45; P < 0.0001; Fig. 5a). Positive ENSO phases are associated with warming sea surface temperatures (SSTs), increased stratification, and a deeper nutricline, leading to negative Chl anomalies in the Equatorial Pacific10,11. Negative correlations were also found between the NAO index and Chl in the North Atlantic (r = −0.31; P = 0.0002; Fig. 5b) and Equatorial Atlantic (r = −0.44; P = 0.001) regions, in accordance with results from Continuous Plankton Recorder surveys29. Positive NAO phases are associated with intensifying westerly winds and warmer SST in Europe and the central North Atlantic30. Possibly, the observed effects relate to increased westerly wind intensity during the winter months, when annual phytoplankton productivity is limited by light availability associated with deep mixed layer depths (MLD)29. We put forward a hypothesis: that an observed coupling of NAO and wind intensity to regional zooplankton abundances29,30 represents a trophic response to the observed phytoplankton fluctuations.

Figure 5: Effects of climate variability on phytoplankton.
figure 5

a–d, Linear relationships between normalized and de-trended yearly Chl anomalies (red) and smoothed climate indices (black) for the Equatorial Pacific (a), the North Atlantic (b), the Arctic (c) and the Southern Ocean (d) regions. Pearson correlation coefficients and P-values are shown. Climate indices and correlation coefficients have been inverted in order to better visualize correlations.

PowerPoint slide

Full size image

No significant relationship was found between the Indian Ocean Dipole index and Chl in the North Indian region (r = −0.23; P = 0.18). The Atlantic Multidecadal Oscillation was positively correlated with Chl in all Atlantic regions (range: r = 0.31–0.43; P < 0.05 for all). Chl anomalies in the Arctic region were negatively correlated with the Arctic Oscillation index (r = −0.31; P = 0.01; Fig. 5c). Chl anomalies in the Southern region were negatively correlated with the Antarctic Oscillation index (r = −0.48; P = 0.029; Fig. 5d), again, possibly owing to intensifying westerly winds and deep MLDs. The strength of all relationships increased after 1950, indicating that phytoplankton may be increasingly driven by climate variability or, alternatively, that model accuracy increased because of increased data availability.

Physical drivers of phytoplankton trends

Long-term trends in phytoplankton could be linked to changes in vertical stratification and upwelling10,11,22, aerosol deposition31, ice, wind and cloud formation8,32, coastal runoff20, ocean circulation33 or trophic effects34. For parsimony, we focus on three variables that may reflect the coupling between physical climate variability and the Chl concentration in the upper ocean: ocean MLD (1955–2009), wind intensity at 10 m (1958–2009) and SST (1899–2009). These physical variables (monthly averages, 1° resolution) were matched by time (year, month) and location (1° cell) with Chl data in order to estimate their effects on Chl within a single model framework (see Supplementary Information for details). SST was the strongest single predictor of Chl. Rising SSTs over most of the global ocean (Fig. 6a) were associated with declining Chl in eight out of the ten regions (range: −0.21 to −0.019 mg m−3 °C−1; P < 0.0001 for all). Positive relationships between SST and Chl were found in the Arctic (0.067 mg m−3 °C−1; P < 0.0001) and Southern regions (0.002 mg m−3 °C−1; P = 0.11). Likewise, inclusion of SST as a covariate in our local models revealed negative SST effects on Chl in 76% (n = 118) of 10° × 10°cells (Fig. 6b). Negative effects prevailed at low latitudes and strong positive effects at high latitudes, particularly in the Southern Ocean (P < 0.05 for all; Fig. 6b and c).

Figure 6: Physical drivers of phytoplankton trends.
figure 6

a, Estimated SST change at 1° resolution from 1899 to 2009. Blue represents cells where SST has declined while yellow and red represent increases. b, Effects of SST changes on Chl estimated for each 10° × 10° cell with >10 yr of data (n = 205). Size of circles represents the magnitude and colours depict the sign of the standardized SST effect on Chl in each cell. c, d, Effects of SST (c) and MLD (d) on Chl as a function of latitude (data points). Relationships were best approximated as quadratic functions of latitude (fitted lines and test statistics). Blue shading represents 95% confidence limits.

PowerPoint slide

Full size image

The effects of SST on Chl are probably explained by its influence on water column stability and MLD10,22. Increasing SST leads to a shallower mixed layer, which further limits nutrient supply to phytoplankton in already stratified tropical waters, but may benefit phytoplankton at higher latitudes where growth is constrained by light availability and deep mixing35. Indeed, in our local models MLD was a significant, but weaker, predictor of Chl concentrations compared with SST, possibly owing to the reduced time series span (1955–2009). Latitudinal gradients in MLD effects were also observed, with predicted positive effects between 20° N and 20° S and negative effects in polar areas (r2 = 0.1; P = 0.018; Fig. 6d). Cumulatively, these findings suggest that warming SST and reduced MLD may be responsible for phytoplankton declines at low latitudes. This mechanism, however, does not explain observed phytoplankton declines in polar areas, where ocean warming would be predicted to enhance Chl (Fig. 6c). This may partially be explained by concurrent increases in MLD and wind intensity there (see Supplementary Fig. 9). Further work is needed to understand the complex oceanographic drivers of phytoplankton trends in polar waters.

Conclusions

Our analysis suggests that global Chl concentration has declined since the beginning of oceanographic measurements in the late 1800s. Multiple lines of evidence suggest that these changes are generally related to climatic and oceanographic variability and particularly to increasing SST over the past century (Fig. 6). The negative effects of SST on Chl trends are particularly pronounced in tropical and subtropical oceans, where increasing stratification limits nutrient supply. Regional climate variability can induce variation around these long-term trends (Fig. 4), and coastal processes such as land runoff may modify Chl trends in nearshore waters. The long-term global declines observed here are, however, unequivocal. These results provide a larger context for recently observed declines in remotely sensed Chl7,10,22, and are consistent with the hypothesis that increasing ocean warming is contributing to a restructuring of marine ecosystems36,37, with implications for biogeochemical cycling15, fishery yields38 and ocean circulation3. Such consequences provide incentive for an enhanced in situ and space-borne observational basis to reduce uncertainties in future projections.

Methods Summary

Data

Available upper ocean (<20 m) in situ Chl data were extracted from the National Oceanographic Data Center (NODC; http://www.nodc.noaa.gov/) and the Worldwide Ocean Optics Database (WOOD; http://wood.jhuapl.edu/wood/). After removing duplicate observations, mean in situ Chl over the upper 20 m was calculated for each cast. Ocean transparency data were extracted from NODC, WOOD, and the Marine Information Research Center (MIRC). Chl (mg m−3) was estimated from transparency measurements as

where D is Secchi depth in metres (ref. 17). As data may be affected by sampling and data entry errors, we filtered erroneous or biologically implausible measurements.

Analysis

Chl trends were estimated for each 10° × 10° cell containing adequate data (‘local’ models, Fig. 2b) and for each regional area (‘regional’ models, Fig. 3b). GAMs were fitted to the blended data to estimate Chl trends as follows:

where η is the monotonic link function of the expected mean Chl concentration µi , B0 is the model intercept, Bi and fi are respectively parametric and non-parametric effects estimated from the data, and εi is an error term. A Γ-distributed error structure and a log link were used. The global mean rate of Chl change (Fig. 3b) was estimated by calculating an inverse variance-weighted random-effects meta-analytic mean39 from the ten regional estimates (see Supplementary Information for full details).

SST changes (Fig. 6a) were estimated by fitting linear models to data in each 1° × 1° cell and area-weighted additive models to data in each of the ten regions. To examine the effects of physical drivers (SST, MLD, wind), Chl and physical data sets were merged by location (1° cell) and time (year, month), and GAMs were fitted with an added effect for the physical driver in question.