Introduction

Phytoplankton is the major primary producer group in marine ecosystems and at the same time a highly suitable indicator for ecosystem change. Due to its high diversity and short generation time, phytoplankton responds quickly to changes in the environment (Reynolds 1998). Especially changes in temperature and nutrient supply are known to change phytoplankton composition and biodiversity, with higher temperatures and lower nutrients favouring small species (Litchman et al. 2007; Winder and Sommer 2012), whereas Si:N and N:P supply ratios change the relative importance of diatoms and other groups such as dino(flagellates) and cyanobacteria (Sommer et al. 2004; Makareviciute-Fichtner et al. 2020).

Consequently, phytoplankton is a significant component of marine environmental monitoring programs. Under the Annex V of the Water Framework Directive (WFD, 2000/60/EC), phytoplankton is defined as one of the “biological quality elements” to assess the ecological status in European coastal waters. Water quality status is assessed by measuring phytoplankton composition, abundance, biomass and bloom frequency and defined by five categories (high, good, moderate, poor bad). The Marine Strategy Framework Directive (MSFD, 2008/56/EC) requires European member states to assess Good Environmental Status (GES), for pelagic habitats by analysing a set of indicators that include information on phytoplankton and zooplankton communities and plankton biodiversity. For the assessment of pelagic habitats (PH) of the Northeast Atlantic, OSPAR (the Northeast Atlantic regional seas convention) has developed indicators based on the abundance of plankton “life forms”, which are a group of organisms that share functional traits (PH1) (OSPAR 2017a), on phytoplankton and zooplankton abundance/biomass (PH2) (OSPAR 2017b) and on species diversity (PH3) (OSPAR 2017c).

Whereas the positive biomass response to high nutrient and light availability is well known (Tilman 1982; Leibold 1999; Siegel et al. 2013), the biodiversity response is much more complex. First, biodiversity itself is a multifaceted construct, comprising the number of different taxa (e.g. species richness) and the equality of species contribution to the assemblage as aspects of local alpha diversity (evenness). But the change in richness is only the net effect of immigration and local extinction, while the turnover in species composition with time is an additional aspect of biodiversity. Extending to larger spatial scales, heterogeneity in space (beta diversity) and composition of the regional species pool (gamma diversity) are further biodiversity dimensions. Second, for each of these facets, all metrics are effort- and scale-dependent. Thus, in contrast to other aspects of ecosystem assessments, biodiversity metrics do not qualify to set absolute target values for good or bad status and effective ecosystem management. But, the relevance of the metrics arises from analysing their temporal trends over a long time period. Third, no single variable captures even the most important aspects of community composition and change. Therefore, the analysis of phytoplankton biodiversity requires a multi-metric approach. For the MSFD assessment requirements, Rombouts et al. (2019) recommend a multivariate approach, which we in general follow with some modifications to reflect recent findings on statistical performance of these metrics (Chase and Knight 2013; Hillebrand et al. 2018).

The complexity of biodiversity assessment has led to the development of other approaches addressing compositional change in plankton communities, most prominently the creation of functional groups. Functional groups are meant to comprise species that show similar responses to environmental drivers and share morphological and/or physiological characteristics. Functional group indicators have been shown to be relevant for describing community structure and biodiversity and are more comparable with other studies than species-based indicators (Mouillot et al. 2006). This approach is also frequently used in marine ecosystem models where varying numbers of functional groups are used (Van Leeuwen et al. 2023) such as the generic ecological model (GEM), which has been developed and improved by different Dutch institutes (Blauw et al. 2009). GEM takes into consideration four phytoplankton functional groups — diatoms, flagellates, dinoflagellates and Phaeocystis — and can simulate phytoplankton processes as primary production, mortality and respiration (Blauw et al. 2009). Here, we analysed the same four functional groups as in the GEMs.

In this study, we showcase how a multifaceted analysis of standing diversity, temporal turnover in species composition and the changes between functional groups allow for a holistic assessment of phytoplankton changes in the Wadden Sea. The Wadden Sea is the largest tidal flat system in the world and is home to a large number of different species. Because of its unique ecological and geological characteristics, it is listed as a UNESCO World Heritage site, but it is a habitat that has been heavily modified and is under intensive pressures by human activities. In this study, we combine two major phytoplankton monitoring time series from the Wadden Sea, which in total comprise 4177 unique phytoplankton samples from 13 stations. We ask whether joint temporal trends can be observed in different diversity metrics and at different locations.

Material and methods

Monitoring stations and data harmonization

The phytoplankton and environmental datasets used in this study are part of the extensive monitoring programs in the Wadden Sea conducted by Rijkswaterstaat, in the Netherlands, and the Lower Saxony Water Management, Coastal Defence and Nature Conservation Agency (NLWKN), in Germany.

We analysed long time series data from 9 Dutch and four German monitoring stations in the Wadden Sea (Fig. 1, Table 1). We used the original datasets but a series of harmonization steps were undertaken. In the phytoplankton data, we first removed all species characterized as purely heterotrophic according to Olenina (2006). Second, we harmonized and updated the species nomenclature between the two datasets (Dutch and German) by using the World Register of Marine Species (WoRMS) open access database (www.marinespecies.org). Third, species-specific biomass was estimated from biovolume using the same C-conversion equations described by Menden-Deuer and Lessard (2000). Irrespective of which method the original programme used, we calculated biomass in pg C cell−1 using one equation for diatoms with biovolumes > 3000 μm3, where C = 0.288 × volume0.811, while for smaller diatoms and other groups, C = 0.216 × volume0.939.

Fig. 1
figure 1

Map of the study area, including the phytoplankton monitoring stations in the Wadden Sea. Data is available for four German stations (triangles) and 9 Dutch stations (circles). Wadden Sea area is represented in blue

Full size image
Table 1 Monitoring stations in the Dutch and German monitoring scheme
Full size table

Sampling frequency

The sampling effort of phytoplankton varied considerably across stations and years in the monitoring datasets (Fig. S1). For most of the Dutch stations, the monitoring frequency was monthly in winter and twice a month in summer (Baretta-Bekker et al. 2009). For the German stations, sampling occurred year-round at two stations (Nney_W_2 and JaBu_W_1) but with different sampling frequency, while the other two stations were sampled in the vegetation period only (March–September/October) (NLWKN 2013). In Bork_W_1, for example, phytoplankton samples were taken over the four seasons from 2007 to 2010, afterwards only in spring, summer and autumn (Fig. S1). WeMu_W_1 was never sampled in winter. In some stations, the sampling frequency increased over time (ex. In DOOVBWT), while in other stations it decreased over time (ex. GROOTGND and HUIBGOT). Also, the number of sampled months per season varied across stations and years. In general, the Dutch stations have longer time series data than the German stations (Table 1).

A more detailed description of the stations and sampling methods is given in Hanslik et al. (1998) for the German stations and Prins et al. (2012) for the Dutch stations.

Annual averaging and biodiversity metrics

A major part of phytoplankton dynamics is between seasons, and these seasonal dynamics can mask more subtle long-term trends, which were the focus of our study. Therefore, we remedied seasonal differences by collapsing the biotic monitoring data to annual medians. Thus, for each species i at station j and year k, we derived the median C-biomass Bijk, for which absences were ignored. The use of a median is superior over means, as the latter are highly influenced by outliers and non-normal distribution of the variables. Moreover, the median can be considered more robust in cases of different sampling regimes between stations. However, in the Supporting Online Material, we document an analysis by which medians were not taken across the entire year but for different seasons (ESM). As these season-specific long-term trends were highly comparable to the outcome of the trends based on annual medians, we focus on the latter within the manuscript.

Based on the annual medians from each species, we first calculated two metrics of phytoplankton standing diversity. We used richness (S) as the raw number of species observed in a year. The incidence-based S treats rare and dominant species equally, which makes S very susceptible to sampling effort as it increases the chance to detect rare species. By contrast, the effective number of species (ENS) is a dominance weighted measure with high statistical robustness (Chase and Knight 2013). It equals the number of species you would encounter in an assemblage having the same diversity but if all species were equally abundant. It can be envisioned as the number of species effectively taking part in the community.

Congruently, we also analysed temporal turnover in two different ways, taking only absence versus the presence or dominance into account. For incidence-based turnover, we used the Jaccard dissimilarity or richness-based species exchange ratio (SERr) (Hillebrand et al. 2018). For abundance-based turnover (SERa), we used Wishart’s dissimilarity as it equally weighs dominance as ENS does (Hillebrand et al. 2018). Both measures range from 0 to 1, with 1 = complete exchange of community and 0 = no change in community composition with time. The analysis of turnover, for both SERr and SERa, can comprise two different time scales: immediate and cumulative. Immediate turnover is calculated between consecutive years and thus reflects whether turnover from 1 year to the next becomes faster or slower with time. Cumulative turnover is calculated between all samples, and by comparing turnover against temporal distance between years, it allows to test whether changes in composition continue (linear relationship between cumulative turnover and distance) or whether previous assemblages are found again, e.g. when assemblages recover from perturbations.

Functional groups’ biomass over time

In this study, we considered four functional groups, namely diatoms, dinoflagellates, flagellates and Phaeocystis. Taxa that did not fit into any of these groups were considered 'other'. We calculated the annual biomass (carbon content in μg L−1) of each functional group as the sum of carbon biomass per sample and functional group and then the median per year.

Statistical analysis

For analysis of statistical trends of standing diversity and turnover over time, we used a mixed effect linear model (LMM) of the form “Response ~ Year + 1|StationID”. Thereby, the LMM allows for additional variation between sample locations through independent intercepts, but still tests for a common linear slope across stations. We constructed separate LMMs for each of the response variables (S, ENS, immediate SERr, immediate SERa and functional group biomass). We analysed the cumulative change in composition by a LMM of the same form, with temporal distance replacing year. All analyses were conducted using R (R Development Core Team 2021) and RStudio (Rstudio Team, 2022).

Although our aim is to analyse temporal trends of phytoplankton in the Wadden Sea, we are dealing with two different monitoring programmes with an unbalanced number of stations and sampling frequency. For this reason, we also present separate models by country in the ESM.

Results

Phytoplankton biodiversity over time

Both metrics of standing diversity, richness (S) and effective number of species (ENS), significantly decline with time across stations (Fig. 2, Table 2), but this negative trend is only associated to the Dutch stations, which in general report more species, but lower ENS (ESM, Table S1). The slope year−1 corresponds to a reduction of S by 17.1 species per decade and 3.02 effective species, which is around 20% of the initial standing S or ENS, respectively (ESM, Table S1). In the German stations, however, there is an increasing trend in richness, corresponding to a gain of 13.5 species in a decade (ESM, Table S1).

Fig. 2
figure 2

Temporal trend of the phytoplankton standing diversity: richness (left) and effective number of species (ENS, right) at the Wadden Sea coastal stations. Overall predicted time effects from the LMM with their confidence interval (grey shaded area) as well as separate line trends for German (green dashed lines) and Dutch stations (orange solid lines). Data input: annual median

Full size image
Table 2 Results of the linear mixed effect model (LMM), analysing the change in phytoplankton diversity and turnover between years. Station ID is included as random effect. σ2 represents the variance component associated with the random effect of Station ID; τ00 measures the variability among the intercepts of different stations in the model, capturing the differences between stations that are accounted for by the random effect; the intraclass correlation coefficient (ICC) indicates the extent to which the change in phytoplankton diversity between years is influenced by the differences between stations; N represents the total number of stations included in the analysis
Full size table

Immediate turnover between neighbouring years is variable over time but does neither speed up nor slow down (Fig. 3a, b, Table 2). However, when separating trends by country, we find a significant increase in SERr in Dutch stations and a decrease in German stations (ESM, Table S2). Cumulative turnover significantly increases with time across all stations (Fig. 3c, d). Richness-based turnover is larger in German stations and abundance-based turnover larger in Dutch stations (Table S3), but both increase consistently with increasing temporal distance. The consistency of the pattern indicates that turnover does involve both shifts in species identity and dominance. It is also noteworthy that the down right corner of the diagram is void of data (Fig. 3c, d), indicating that there is no “return” to a previous assemblage over long time scales, indicating a strongly directional compositional drift over time (Table 3).

Fig. 3
figure 3

Temporal trend of phytoplankton turnover at the Wadden Sea coastal stations. a, b Turnover measured between years. Separate line trends for German (green dashed lines) and Dutch stations (orange solid lines); c, d turnover accumulated over years. Green triangles represent German stations and orange circles represent Dutch stations. Overall predicted time effects from the LMM with their confidence interval (grey shaded area). Data input: annual median

Full size image
Table 3 Results of the linear mixed effect model (LMM), analysing the change in phytoplankton cumulative turnover over years (dist). Station ID is included as random effect. σ2 represents the variance component associated with the random effect of Station ID; τ00 measures the variability among the intercepts of different stations in the model, capturing the differences between stations that are accounted for by the random effect; the intraclass correlation coefficient (ICC) indicates the extent to which the change in cumulative turnover over years is influenced by the differences between stations; N represents the total number of stations included in the analysis
Full size table

Biomass and functional groups over time

Overall phytoplankton biomass increases over time, but this trend is only associated with the Dutch stations (Fig. 4a). In most of these stations, we observe an increase in biomass of the functional groups, but especially of diatoms and of the ‘other’ group (Fig. 4, Table 4 and ESM Table S4). Diatoms have the highest biomass over all functional groups in both countries, whereas dinoflagellates show no overall trend as their biomass varies across stations (Table 4). Although the biomass at the German stations does not change significantly over the period analysed, there is a decrease in biomass of flagellates and of the ‘other’ group (ESM Table S5).

Fig. 4
figure 4

Annual median of phytoplankton biomass (a); at the Wadden Sea stations, also shown separately for diatoms (b); dinoflagellates (c), flagellates (d), Phaeocystis (e), and ‘other’ (f). Overall predicted time effects from the LMM with their confidence interval (grey shaded area) as well as separate line trends for German (green dashed lines) and Dutch stations (orange solid lines). Data input: annual median

Full size image
Table 4 Results of the linear mixed effect model (LMM), analysing the change in phytoplankton biomass (LN μgC L−1) and biomass of each functional groups (LN μgC L−1) between years. Station ID is included as random effect. σ2 represents the variance component associated with the random effect of Station ID; τ00 measures the variability among the intercepts of different stations in the model, capturing the differences between stations that are accounted for by the random effect; the intraclass correlation coefficient (ICC) indicates the extent to which the change in phytoplankton biomass between years is influenced by the differences between stations; N represents the total number of stations included in the analysis. A separate model for each country is included in the ESM (Table S4 and S5)
Full size table

Discussion

The investigation of almost 20 years of phytoplankton monitoring data from the Wadden Sea, including data from Germany and the Netherlands, provided an opportunity to better understand temporal changes in the diversity of primary producers and in the biomass of their functional groups. Our main results show that despite the proximity and connectedness of Dutch and German Wadden Sea stations, trends in biomass, diversity and turnover differ, with only Dutch stations experiencing an overall decline in both richness and ENS but an increase in carbon biomass. Immediate richness-based turnover also increases in the Dutch stations, but it decreases in the German ones. But accumulated turnover over time (SERr and SERa) increases in both countries.

We are not able to fully resolve the question whether the different trends are real or partly based on systematic differences in the monitoring protocols. Although the Dutch and German data are within themselves highly curated, it is important to note that the monitoring programs themselves are not harmonized. Consequently, technical differences and different taxonomic schools may contribute to some of the observed differences, which, considering the regular tidal flooding of this intertidal system, may appear significant. We are analysing long time series of data where probably different analysts have counted the phytoplankton samples. It is well known that the personal component has an influence on how well species are identified. Even the same person counting the samples can evolve over time, gaining more experience, which can also influence the results (Löder et al. 2012; Nohe et al. 2018). A recommendation for further improvement of the assessment is a case wise cross-validation of some samples by the monitoring agencies involved. Another important issue to be considered are the differences in sampling frequency between countries and also between stations (see ESM Fig. S1). In highly dynamic and turbid waters such as the Wadden Sea, low frequency and infrequent sampling is not sufficient to properly capture temporal signals (Fettweis et al. 2023).

The decrease in richness and ENS observed in the Dutch stations, together with the increase in biomass, indicates that fewer species are becoming more dominant in the system. In addition, one of our main conclusions is the continuous shift in species composition over time. This was — in contrast to some of the other results — consistent across stations and time, indicating that the dynamics of compositional change were characterized by a directional drift. Most functional groups (except dinoflagellates with a weaker slope) showed similar increases in biomass over time, the compositional change thus mainly occurred at lower taxonomic levels and only at the Dutch stations. Addressing the size distribution of phytoplankton in the German Wadden Sea, Hillebrand et al. (2022) observed an increase in smaller cells both by smaller species becoming more important and by cells within species declining in cell size. They relate this to higher water temperature and lower nutrient concentrations, an observation confirmed at larger spatial scales (Marañón et al. 2012; Peter and Sommer 2013; Mousing et al. 2018).

The temporal decline in cell size observed in the German stations (Hillebrand et al. 2022) could give an explanation for the different trends in carbon biomass observed between the Dutch and German stations. In the Dutch monitoring programme, cell size is not measured per sample. Instead, values are obtained from the literature and no temporal variability in cell size is recorded. Therefore, one hypothesis for the linear increase in carbon biomass (biomass per litre was measured as abundance × cell size × carbon factor, see “Methods”) in the Dutch stations and absent trend in the German stations could be the decrease in cell size observed in German samples (see Hillebrand et al. 2022), which counterbalance the biomass increase. In other words, if cells are getting smaller, smaller organisms become predominant (Hillebrand et al. 2022), the linear increase in biomass in the Dutch stations would have been offset or amenable if the decrease in cell size had been taken into account. As we cannot confirm this hypothesis, we strongly emphasise the importance of measuring cell size from the samples.

But why is the biomass as carbon in the Wadden Sea not decreasing, even though nutrient levels have been decreasing since the 1990s? The highest nutrient concentrations in the Wadden Sea were recorded between 1980 and 1990, which forced policy measures to be taken to reduce eutrophication. Since then, nutrient loads have continuously decreased (Philippart and Cadée 2000; van Beusekom et al. 2001, 2009; Jung et al. 2017), with phosphorus (P) showing a stronger decline than nitrogen (N) (van Beusekom et al. 2019). This suggests that P has become the main limiting nutrient in the region (Philippart et al. 2007; Kuipers and van Noort 2008; Ly et al. 2014; Leote et al. 2016). It would be expected that a reduction in nutrient loads would result in reductions in phytoplankton biomass, but our results do not confirm this expectation. De Jonge (1997), analysing Wadden Sea data up to 1995, suggested that the remaining high primary production in the most western Dutch Station (MARSDND) could be explained by other sources of nutrient inputs, such as from the English Channel. Other studies suggest that P may remain in sediments for long periods of time, indicating that even if inputs of P from external sources have decreased, it is still present at high levels in the sediment and eventually gets released into the water (de Jonge et al. 1993; Leote and Epping 2015; Leote et al. 2016; Jung et al. 2017). Therefore, P released from the sediment could play a crucial role in nutrient availability in the Wadden Sea, and it has been suggested to be the main source of PO4, at least in the western part of the Wadden Sea (Leote et al. 2016). Despite the decreased concentrations of P and N, silicate (Si), which plays a critical role in diatom growth, demonstrated minimal changes over time (Rick et al. 2023) or even showed an increase during the phytoplankton growing season (Prins et al. 2012). These findings could potentially explain the sustained phytoplankton biomass despite the overall decline in nutrient levels.

In addition to nutrient availability and nutrient input from various rivers (Rhine, Ems, Weser and Elbe), phytoplankton growth in the Wadden Sea is strongly influenced by the light availability, as it is a very turbid system, by zooplankton and benthic grazing. Also, the mixing of fresh water with sea water in the estuaries results in strong physical and chemical gradients that influence phototrophic production, composition and diversity (McLusky 1993; Muylaert and Sabbe 1999). Van Beusekom et al. (2019) suggested that the increase in phytoplankton biomass at the DANTZGT station was due to the improvement in light conditions observed since the 1980s. Light conditions and differences in local hydrography were found to be the most limiting factors for phytoplankton on the North Sea coast in a study including monitoring stations in Belgium, the Netherlands and Germany (Loebl et al. 2009). Therefore, it is reasonable to hypothesize that the observed increase in biomass could be attributed, at least in part, to the improved light availability in the Wadden Sea over the past few decades.

Conclusion

Our analysis of biodiversity of Wadden Sea phytoplankton strengthen previous conclusions that a multivariate approach is needed to capture even the most important aspects of community compositional change (Rombouts et al. 2019). We suggest that our 2 × 2 factorial assessment, which on the one hand uses incidence- and dominance-based measures and on the other hand trends in alpha diversity as well as turnover could prove useful for and beyond the Wadden Sea ecosystems as it allows the detection of detailing changes in composition and overall biodiversity. We have also shown that temporal trends in phytoplankton dynamics can be different or even opposite in stations that are very close to each other but monitored by different countries. These differences may be artefacts due to technical reasons or real differences due to regional environmental conditions. As we cannot confirm either of these hypotheses, we strongly recommend that monitoring programs go through a far-reaching process of harmonization to make results more comparable between countries. This would allow us to draw more congruent conclusions on the physical and biological status of an ecosystem, especially a transboundary ecosystem such as the Wadden Sea, and to better detect potential changes related to different drivers.