Introduction

At a planetary scale, the eukaryotic phytoplankton group collectively known as diatoms contribute approximately one-fifth of all of Earth’s photosynthesis1 and up to 30–40% of global marine primary production each year2. They are particularly important in boreal systems where they dominate seasonal phytoplankton production and occupy key ecosystem roles3,4. Diatoms and the role of nutrients and ocean biogeochemistry on their biological patterns have been at the forefront of intense studies5,6,7. Excluding mathematical models8, least known, however, is the role of climate on in situ phytoplankton populations particularly at macro and multidecadal scales. How climate variability and diatom abundance interact at these scales is paramount to our understanding of the effects and prediction of the future impact of climate change on pelagic communities as well as having important future consequences to carbon sequestration (the biological carbon pump) and North Atlantic productivity.

Due to costs in maintaining consistent measurements over time and logistical complexities, multidecadal observations are rare. Here we used one of the few ecological time-series at this scale available derived from the Continuous Plankton Recorder (CPR) Survey to investigate spatial and temporal patterns of diatoms over 60 years in the NE Atlantic and adjacent seas. We relate these spatial and temporal patterns to the changing climatic conditions in the region, which include climate warming and wind patterns as well as an index that reflects Atlantic Multidecadal Variability (AMV).

Results

We first started to examine long-term changes in total diatom abundance from January 1958 to December 2017 (see methods9) in the North Sea where sampling by the CPR survey has been most dense. For reasons of statistical robustness, for the analysis of phytoplankton at the macroscale level we used diatom abundance for the whole NE Atlantic but for individual species trends we only used data from the North Sea region. For this analysis, we included a total of 49 species/taxa (Supplementary Table 1). Annual diatom abundance (Fig. 1) is dominated by the spring bloom with some further seasonal peaks occurring in autumn following the classic diatom pattern found in temperate and boreal marine environments3. Peak diatom abundance occurs between March and May followed by a decline during the summer months while the water-column becomes stratified and is then followed by a floral renaissance in production in the autumn when the water column becomes vertically mixed again. Long-term trends in total diatom abundance show a sharp decadal decline during the 1960s which continued to a very low point in the late 1970s associated with anomalously cold and low-salinity hydro-climatic conditions before increasing again in the mid 1980s10,11. After recording low values in the early 1990s, diatom abundance substantially increased again in the mid-1990s and has continued to remain high since this phase shift in abundance (Fig. 1, all trends statistically significant using Pearson correlation coefficient). While diatom phenology has remained fairly temporally static over time12, in particular the timing of the spring bloom, there have been changes to the amplitude in abundance over the seasonal cycle (Fig. 1). Particularly noticeable is the sharp decline in the autumnal bloom of diatoms that occurred during the 1960s. Over the last 20 years there appears to have been an increase in diatoms across the seasonal growth period with a return to increased autumnal growth and more abundance dispersed throughout the seasonal cycle including the summer period. In particular, the taxa Thalassiosira spp., Rhizosolenia imbricata and Pseudo-nitzschia seriata complex have increased their seasonal spread and autumnal bloom formation over the last few decades. The autumn bloom is less well studied than its vernal cousin but is a significant contributor to overall annual phytoplankton production13.

Fig. 1: Total monthly diatom abundance in the North Sea from 1958 to 2017.
figure 1

a Total diatom (49 taxa) counts per month. b Total diatom counts per month log transformed log10(x + 1).

Full size image

Subsequently, we examined individually species trends over time of the most common 49 species/taxa recorded by the CPR Survey (Fig. 2). A Principal Component Analysis (PCA) was performed on the long-term trends of the diatom species/taxa at the scale of the North Sea. Where we determined the first (Fig. 2a) and second (Fig. 2b) principal components (PCs, 26.14% and 11.01% of the total variance, respectively) and their significance using a broken-stick test (see Methods). The first PC clearly shows a large decline in the 1960s with some notable low values in the late 1970s and early 1990s before a shift in the time-series after the mid-1990s. Examples of the highest species/taxa correlations with the first principal component were mainly neritic spring diatom species which included Paralia sulcata (the only tychopelagic species), Bellerochea malleus, Biddulphia alternans, Thalassionema nitzschioides, Asterionellopsis glacialis, Pseudo-nitzschia delicatissima complex, Pseudo-nitzschia seriata complex, Rhizosolenia imbricata and Rhizosolenia setigera (Supplementary Table 1). The second PC (Fig. 2b) shows a clear positive linear trend through time. Highest positive species/taxa correlations for this PC were with Gyrosigma spp. and Rhaphoneis amphiceros, and high negative correlations included Chaetoceros (Phaeoceros) spp., Chaetoceros (Hyalochaete) spp. and Odontella regia (Supplementary Table 1). Most of the species associated with a positive PC2 are autumnal bloom forming species that are found in shallow neritic waters of the North Sea.

Fig. 2: Long term trends of all diatom species by using a Principal Component Analysis in the North Sea from 1958–2017.
figure 2

a Principal Component (PC) 1 (bold line) and associated 49 species trends (light lines). Standardised abundance of all species related to the first PC (species with a correlation with PC1 higher than 0.7). b Principal Component 2 (bold line) and associated 49 species trends (light lines). Standardised abundance of all species related to the second PC (species with a correlation with PC2 higher than 0.7). Species were selected on the basis of the examination of normalised eigenvectors.

Full size image

We then examined the relationships between long-term trends in diatom abundance and AMV at the scale of the Northeast (NE) Atlantic by using the Atlantic Multidecadal Oscillation (AMO) index and Northern Hemisphere Temperature (NHT) anomalies (see Methods). Patterns of changes in total diatom abundance were similar to those observed in the North Sea (Fig. 2a versus Fig. 3b). The positive normalised eigenvector 1 (Fig. 3a) shows that the trend observed in the first PC (28.92% of the total variance) reflected changes that took place in most regions of the Northeast Atlantic (Fig. 3a, b). The multidecadal trend between the diatom PC1 and the AMO index is remarkably similar with positive correlations found throughout the NE Atlantic (r = 0.71, pACF < 0.05; Fig. 3b). Long-term trends in the AMO index and the diatom PC1 show a strong decline in the early part of the time-series, low values between the 1970s to the 1990s and a sharp shift in both time-series occurring in the mid-1990s (Fig. 3b). High values have remained since this phase shift up to the current time-period examined.

Fig. 3: Long-term (1958–2017) changes in diatom abundance as revealed by Principal Component Analysis (PCA) performed at the scale of the Northeast Atlantic and relationships with the Atlantic Multidecadal Oscillation (AMO) and Northern Hemisphere Temperatures (NHT).
figure 3

a First normalised eigenvector showing the spatial pattern of correlations between spatial changes in total diatom abundance and the first principal component (PC). b Long-term relationship between the first PC of diatom abundance (blue) and the AMO index (red). c Second normalised eigenvector showing the spatial pattern of correlations between spatial changes in total diatom abundance and NHT. d Long-term relationship between the second PC of diatom abundance (blue) and NHT (red).

Full size image

Long-term changes in the second PC (8.46% of the total variance) reflects an opposition between a negative trend in total diatom abundance in the southern part of the NE Atlantic and a positive trend in the northeastern part (Fig. 3c, d). Both time-series PC2 and NHT anomalies show a strong upward trend in values from the 1960s to the present (r = 0.93, pACF < 0.05). Because of the spatial patterns observed in the second normalised eigenvector, strong positive correlations between diatom abundance and the NHT are mainly spatially found in the shallow regions of the North Sea and off southern Norway (Fig. 3c). In contrast, negative correlations between diatom abundance and NHT are geographically found in the Bay of Biscay.

To investigate the relationship between wind and temperature in detail we examined spatio-temporal correlations between diatom abundance and various components of wind and Sea Surface Temperature (SST) in the NE Atlantic. Figure 4 shows the correlation maps between total diatom abundance and SST (Fig. 4a), U wind (zonal or westerlies) (Fig. 4b), V wind (meridional or southerlies) (Fig. 4c) and W wind (wind intensity) (Fig. 4d) for the NE Atlantic over the 60 years time period (see Methods). There is a positive correlation between temperature and diatoms found in the northern regions of the NE Atlantic such as the northern North Sea and the Icelandic basin region (Fig. 4a). Negative correlations between total diatom abundance and SST are found in the Celtic Sea and Bay of Biscay (Fig. 4a). Negative correlations are found between U wind and diatoms in the oceanic region south of Iceland (Fig. 4b) Positive correlations are found between V wind and total diatom abundance in the central North Sea and in some central oceanic regions of the NE Atlantic (Fig. 4c). Finally, the overall wind intensity is negatively correlated with diatom abundance with this signal mainly found in the central North Sea. To examine these temporal correlations in more detail we focused on correlation coefficients between these various climatic indices and the total diatom abundance found in the North Sea between different seasons. Supplementary Table 2 shows the seasonal correlation coefficients between diatoms and the AMO, North Atlantic Oscillation (NAO) index, SST, W wind, U wind and V wind. Very low correlations (excluding the AMO) are found between all climate variables and spring diatom abundance. A significant positive correlation is found between summer diatoms and the AMO. During autumn the abundance of diatoms in the North Sea are significantly correlated with temperature (AMO and SST) and negatively correlated with U wind and wind intensity (Supplementary Table 2).

Fig. 4: Long-term spatial and temporal (1958–2017) correlations between total diatom abundance and climate variables in the Northeast Atlantic.
figure 4

a diatoms and Sea Surface Temperature (SST), b diatoms and U wind (westerlies), c diatoms and V wind (southerlies), d diatoms and W wind (intensity).

Full size image

Discussion

In this study, we have established a strong relationship between multidecadal variability in the Atlantic and the abundance of diatoms in the NE Atlantic and have also established a significant relationship between trends in neritic diatoms in the North Sea and Northern Hemisphere Temperature. In the temperate and boreal North Atlantic where CPR sampling is most prolific there have been a significant change in climate due to natural variability and global warming over the duration of this multi-decadal study. At the global level climate warming of pelagic surface waters is thought to increase the density gradient between the surface layer and the underlying nutrient-rich waters, in turn enhancing the stability of the water-column and hence restricting nutrient availability. Availability of the principal nutrient (nitrate) that limits phytoplankton growth has been found to be negatively correlated to global temperatures14,15. In another global analysis using satellite derived chlorophyll data it has been shown that a strong inverse relationship exists between SST and chlorophyll concentration16. In terms of climate feedback mechanisms, it is thought that phytoplankton driven biological carbon pumps will be less effective in warmer conditions due to changes mainly in the phytoplankton assemblage. These changes favour small flagellates8 to the disadvantage of diatoms, and less overall nutrient mixing due to enhanced stratification. It is also predicted that warmer conditions would shift the metabolic balance between production and respiration in the world’s oceans towards an increase in respiration thus reducing the capacity of the oceans to capture CO217. However, in the case of our study, the manifestation of climate variability in the NE Atlantic is less than homogenous with distinct spatial and temporal patterns emerging with diatom abundance increasing in some areas and declining in others. In the NE Atlantic, a prevalent multidecadal trend has been established where climate warming is increasing diatom populations in more northerly regions (e.g. Icelandic Basin), but deceasing populations in more southerly regions (e.g. Bay of Biscay).

Over the whole NE Atlantic there has been an increase in phytoplankton biomass during spring and autumn (where diatoms dominate) with increasing temperatures in cooler regions but a decrease in phytoplankton biomass in warmer regions18. This is possibly due to increased phytoplankton metabolic rates caused by warming temperatures in colder regions but conversely a decrease in nutrient supply in warmer regions. This pattern is probably only applicable to oceanic regions of the NE Atlantic such as the Bay of Biscay and the Icelandic Basin and not as important to diatoms in well-mixed neritic waters such as those found in regions of the North Sea (as nutrient replenishment is not cut off). Our study also shows that autumnal diatom abundance is positively correlated with SST and NHT. Regional climate warming in some areas of the North Sea has been linked to an increase in certain diatoms that are associated with Harmful Algal Blooms, in particular the taxa Pseudo-nitzschia seriata complex19,20. Diatom growth in such well mixed areas may be enhanced by temperature as these regions are not inhibited by stratification and hence nutrient availability (see Fig. 3b).

Although it is difficult to tease out absolute mechanisms between climate and phytoplankton due to the nature of macroscale correlative studies such as this analysis, our study does shows that the AMO has a very strong influence on diatoms in the NE Atlantic at the multidecadal scale. During a negative phase of the AMO, total diatom abundance is reduced despite the fact that temperature decreases and that wind increases. We observe the opposite effect during a warm phase of the AMO and show that the influence of the AMO mainly manifests itself during the end of summer and in autumn for the North Sea regions. Higher temperatures and less wind mixing during the autumn months could enhance diatom growth in these well mixed regions. It also appears that since the positive phase shift in the AMO beginning the mid-1990s to the present time has substantially tightened the relationship between diatom abundance and climate warming (AMO and NHT). These findings are important because mechanisms will differ between different regional areas of the world with varying hydro-climatic regimes. This study illustrates the complexity of the responses of phytoplankton and their functional types to climate variability as regional nuances will not always be represented well by homogenised global models.

Methods

Biological data

The CPR Survey is a long-term, sub-surface marine plankton monitoring programme consisting of a network of CPR transects towed monthly across the major geographical regions of the North Atlantic. It has been operating in the North Sea since 1931 with some standard routes existing with a virtually unbroken monthly coverage back to 194621,22. The CPR survey is recognised as the longest sustained and geographically most extensive marine biological survey in the world23. The dataset comprises a uniquely large record of marine biodiversity covering ~1000 taxa over multidecadal periods. The survey determines the abundance and distribution of phytoplankton and zooplankton (including fish larvae) in our oceans and shelf seas. Using ships of opportunity from ~30 different shipping companies, it obtains samples at monthly intervals on ~50 trans-ocean routes. In this way the survey autonomously collects biological and physical data from ships covering ~20,000 km of the ocean per month, ranging from the Arctic to the Southern Ocean24.

The CPR is a high-speed plankton recorder that is towed behind ‘ships of opportunity’ through the surface layer of the ocean (~10 m depth)9. Water passes through the recorder and plankton are filtered by a slow-moving silk (mesh size 270 µm). A second layer of silk covers the first and both are reeled into a tank containing 4% formaldehyde. Upon returning to the laboratory, the silk is unwound and cut into sections corresponding to ten nautical miles and ~3 m3 of filtered sea water25. There are four separate stages of analysis carried out on each CPR sample, with each focusing on a different aspect of the plankton: (1) overall chlorophyll (the phytoplankton colour index; PCI); (2) larger phytoplankton cells (phytoplankton) in the case of this study we have just analysed diatoms which are all on the large size compared with smaller flagellates and coccolithophores that are also recorded.; (3) smaller zooplankton (zooplankton ‘traverse’); and (4) larger zooplankton (zooplankton ‘eyecount’). The collection and analysis of CPR samples have been carried out using a consistent methodological approach, coupled with strict protocols and Quality Assurance procedures since 1958, making the CPR survey the longest continuous dataset of its kind in the world24. Of the more than 200 phytoplankton species recorded by the Survey we analysed a sub-set of data consisting of over 49 diatom species/taxonomic entities (Supplementary Table 1). The selection of taxa was based on their frequency of occurrence, consistency of recorded taxa across the whole time-series (1958–2017). Rare species were removed during the selection to make the dataset statistically more robust. Practically, we chose diatoms that had an annual abundance >0 for at least 50 years.

Potential biases and limitations

The CPR survey currently records ~1000 plankton entities (many to species level) in routine taxonomic analysis with strict Quality Assurance protocols in operation for its sampling and plankton counting procedures since 1958. Due to the mesh size of CPR silks, many plankton species are only semi-quantitatively sampled owing to the small size of the organisms. In the case of phytoplankton there is thus a bias towards recording larger armoured flagellates and chain-forming diatoms and that smaller species abundance estimates from cell counts will probably be underestimated in relation to other water sampling methods24. However, the proportion of the population that is retained by the CPR silk reflects the major changes in abundance, distribution and specific composition, i.e. the percentage retention is roughly constant within each species even with very small-celled species down to 10 µm19,24. In certain circumstances clogging due to increased loading associated with mucilage and high densities of plankton may contribute to the number of small plankton captured. Despite these sporadic occurrences the CPR has still been shown to capture a consistent fraction of the in situ abundance of the species assayed25. A similar potential under estimation of zooplankton abundances has recently been thoroughly statistically explored and found that while the CPR survey does underestimate abundance in some cases it provides a realistic quantification of both temporal (i.e. seasonal and diel scales) and spatial (i.e. regional to basin-scale) changes in zooplankton taxa (In this case for the species Calanus finmarchicus)24,26. The study also showed that while the CPR sampling is restricted to surface waters ~10 m in depth the seasonal and diel patterns of abundance of C. finmarchicus were positively correlated to patterns of abundance to a depth of 100 m26.

While all sampling and analytical tools for measuring plankton have their own individual strengths and weaknesses, the CPR Survey has had perhaps the most well documented and transparent examination of its consistency and compatibility of any ecological time series over time21,23. Perhaps the most obvious limitation of CPR sampling is its underestimation of components of the plankton, e.g. large plankton like fish larvae and delicate gelatinous plankton. This has been well documented and users of the data are advised of the CPR’s potential biases24,25. It is also widely recognised that all plankton sampling systems have their own limitations and nuances and all underestimate abundance to some degree and that the varying mechanisms are not always directly comparable24,27. A detailed study has been conducted on flow rate and ship speed on CPR sampling28 given that the speed of the ships has, in some circumstances, increased since the 1960s, which may impact sample efficiencies. However, no significant correlation was found between the long-term changes in the speed of the ships and two commonly used indicators of plankton variability: Phytoplankton Colour and Total Copepods indices24. This absence of relationship may indicate that the effect found is small in comparison with the influence of hydroclimatic forcing24,28. For further details on the technical background, methods, consistency, and comparability of CPR sampling, see ref. 21.

Hydroclimatic indices

The North Atlantic experiences climate variability on multidecadal scales often referred to as AMV29. The AMV has been operating over the last few centuries as an observed phenomenon and has also been identified at millennial timescales using paleoclimatic proxy data and models. However, identifying the mechanisms behind the AMV, be they externally forced or driven by internal dynamics or a combination of both is still widely debated29,30,31. Here we use the AMO index to reflect this low frequency variability in the North Atlantic. The AMO is a large-scale oceanic phenomenon and a source of a natural variability in the range of 0.4 °C in the oceanic regions of the North Atlantic32. The AMO is also detected in paleoclimatic reconstructions with a cycle ranging from 60 to 100 years33. The index is calculated from the Kaplan SST dataset which is updated monthly and represents an index of the North Atlantic temperatures; we used the unsmoothed version of the index (http://www.esrl.noaa.gov/psd/data/timeseries/AMO/). The AMO has been related to various biological shifts in the North Atlantic sector34.

NAO index is a basin scale alternation of atmospheric mass over the North Atlantic between the high pressures centred on the subtropical Atlantic and low pressures around Iceland and is a prominent mode of variability with important impacts from the polar to subtropical Atlantic and surrounding landmasses. The index is calculated on the surface sea-level pressure difference between the Subtropical (Azores) High and the Subpolar Low (https://www.ncdc.noaa.gov/teleconnections/nao/).

We also used Northern Hemisphere temperature anomalies (HadCRUT4). This climatic index was developed by the Climatic Research Unit (University of East Anglia) and the Hadley Centre (UK Met Office)35. CRUTEM4 and HadSST3 are the land and ocean components of this dataset, respectively (https://crudata.uea.ac.uk/cru/data/temperature/). Global temperature anomalies have been frequently correlated with long-term changes in biological and ecological systems36.

Climatic data (wind and sea surface temperature)

We used wind intensity (W wind), the zonal (U wind, from west to east) and meridional (V wind, from south to north) components of the wind to investigate potential relationships between wind intensity and direction and diatom abundance. Wind intensity and direction influence the structure of the water column and are known to be important environmental parameters for diatoms37. Wind intensity and direction (monthly basis and spatial resolution of 2.5° latitude × 2.5° latitude) was obtained from the National Centers for Environmental Prediction and the National Center for Atmospheric Research, who have run a reanalysis project using a state-of-the-art analysis/forecast system to perform data assimilation using past data from 1948 to the present38. For the analysis we constructed a spatially gridded (2.5° latitude × 2.5° longitude) dataset for the northeast Atlantic region (40.5°N to 65.5°N and 30.5°W to 10.5°E) and the North Sea (52°N to 60°N and 2°W to 10°E), corresponding to the period 1958–2017.

SST was extracted from the Extended Reconstructed Sea Surface Temperature (ERSST)_V5 data-set, which is derived from a reanalysis of the most recently available International Comprehensive Ocean-Atmosphere Data Set SST data. Improved statistical methods were used to produce a stable monthly reconstruction from relatively sparse data39. For this analysis we constructed a spatially gridded (2° latitude × 2° longitude) annual average SST dataset for the Northern Hemisphere.

Spatialised standardised principal component analysis on total diatoms in the northeast Atlantic

We calculated the sum of the abundance of all diatoms recorded semi-quantitatively by the CPR survey (264,837 samples). Because data from the CPR survey are not regularly spaced in time and space, we interpolated the data for total diatoms on a 2-monthly basis between 1958 and 2017 using an inverse squared distance method with a search radius of 150 km and a number of values for calculating the weighted mean ranging between 3 and 2040,41. The boundary of the studied region ranged from 40.5°N to 65.5°N and from 30.5°W to 9.5°E for every 1° latitude × 1° longitude; inside the area, we kept geographical cells that had at least 50 non-zero estimates of total diatom abundance. Annual means were subsequently calculated for every geographical cell. We kept cells that had at least 30 years of estimates of total diatom abundance.

A spatialised standardised PCA was performed on the matrix 1066 geographical cells (26 latitudes × 41 longitudes) × 60 years. The matrix of standardised eigenvectors 1066 × 1066 geographical cells represented the linear correlations with the principal components that had a dimension of 60 years41. We applied a broken-stick model to test the significance of the different axes of the PCA using 10,000 simulations42.

Standardised principal component analysis on species of diatoms in the North Sea

We calculated the abundance of all diatom species (including total diatom abundance) occurring in the North Sea (52°N to 60°N and 2°W to 10°E) on a monthly basis for the period 1958–2017. Long-term monthly changes in total diatoms were examined first in relation to the hydro-climatic indices (AMO, NAO and NHT). Data were then annually averaged and a standardised PCA was performed on the matrix 60 years × 49 species/taxa. This analysis therefore focussed on all species or taxa of diatoms, which was not investigated on the first PCA that jointly considers spatial and temporal changes in total diatom abundance. Here also, we applied a broken-stick model to test the significance of the different axes of the PCA using 10,000 simulations42.

Relationships between changes occurring in the northeast Atlantic and the North Sea

We correlated the long-term changes in the first principal components (PCs) originating from the PCA performed at the scale of the North Atlantic with changes in diatom abundance at the species or taxonomic level. For this analysis, we represented changes in diatom species in addition to a PC when the PC-species correlation was higher than 0.7 (i.e. percentage of explanation >49%).

Maps of correlations

We calculated maps of correlation between W, U, V wind and SST and the first PC originating from the PCA performed on total diatoms at the northeast Atlantic scale. For this analysis, we applied an order-2 simple moving average to alleviate the influence of short-term higher frequency variability. This analysis was conducted for the time period 1958–2017.

Correlation analyses between variables

Correlations were calculated to assess the strength of the relationships between total diatom abundance and climate in the North Sea for different seasons. The probability of significance of the coefficient of correlation, called pACF, was adjusted to correct for temporal autocorrelation.