Changes in Phytoplankton and Zooplankton Communities

Indicators based on plankton lifeforms have been used to assess community response to sewage pollution (Charvet et al 1998; Tett et al., 2008), anoxia (Rakocinski 2012), fishing (Bremner et al., 2004), eutrophication (HELCOM 2012), climate change (Beaugrand 2005; Bedford et al., 2020a; McQuatters-Gollop et al., 2019), and ocean acidification (Keys et al., 2018). Indicators based on functional groups have been proven relevant for the description of the community’s structure and biodiversity and are more easily inter-compared than other indicators based on taxonomy (Estrada et al., 2004; Gallego et al., 2012; Garmendia et al., 2012; Mouillot et al., 2006).

In practice, it is often preferable to aggregate species with similar traits into functional groups, such as lifeforms (Figure 1), rather than assessing the dynamics of individual species. Measures of species abundance are frequently subject to large interannual and regional variation, often due to natural physical dynamics and habitat preferences rather than anthropogenic stressors (de Jonge, 2007). Functional group abundance is often less variable because variability in the abundances of the group’s constituent species averages out. Cryptic speciation (species with near-identical appearance) within the plankton community, alongside the limitations of identifying plankton using routine light microscopy techniques, make it difficult to generate accurate counts at a species or genus level. Functional group abundance is more reliable as many plankton lifeforms are easily identified, making comparisons between different laboratories and institutes feasible. Both abundance and biomass data can be used to inform lifeform time series, depending on the lifeform in question and data availability from monitoring programmes.

In addition to studying change in individual lifeforms, these concepts can be extended to investigating changes in ecologically relevant pairs of lifeforms in the form of a lifeform pairs index. The precise combination of lifeforms composing the pairs will depend on the habitat and the objective of the indicator, e.g., as a measure for change in pelagic habitats, food webs, seafloor integrity or eutrophication. Change in the abundance of ecologically linked lifeforms over time can also provide an indication of changes in various aspects of ecosystem function. These include, for example: the transfer of energy from primary to secondary producers (changes in phytoplankton and zooplankton); the pathway of energy flow and top predators (changes in gelatinous zooplankton and fish larvae); benthic / pelagic coupling, i.e., changes in holoplankton (fully planktonic) and meroplankton (only part of the lifecycle is planktonic, the remainder is benthic) (see also Gowen et al., 2011; McQuatters-Gollop et al., 2019).

While studying lifeform pairs can be helpful for detecting changes in the annual cycle of ecologically linked lifeforms, assessing lifeforms separately is more suitable for evaluating gradual change over time. It is simpler to interpret results of an analysis of correlation between lifeform abundance and environmental pressures, than between an index of lifeform pairs and environmental pressures (Bedford et al., 2020a). A method has been developed to detect potential links between environmental pressures and change in lifeform abundance over time. The main features of the method are: (i) the grouping of planktonic taxa into functional types or lifeforms; (ii) spatially division of plankton samples to construct distinct time-series of lifeform abundance; (iii) using a robust nonparametric test to quantify long term changes in lifeform abundance; (iv) relating change in lifeform abundance to trends in environmental pressures and climate indices.

The OSPAR Intermediate Assessment 2017 (IA2017) stated that conclusions on the trends observed could not yet be made as additional work was needed to investigate environmental pressures and potential impacts of climate change and other human pressures on plankton communities. This QSR 2023 assessment addresses plausible links between environmental pressures and long term changes in plankton communities.

Previous assessment

The previous assessment (Intermediate Assessment 2017; IA2017) was based largely on state space theory (see below) and examined differences in the relative abundances of ecologically relevant lifeform pairs between a past comparison period (2004 to 2008) and a contemporary assessment period (2009 to 2014). While this method is suitable for detecting changes in the annual cycle of ecologically linked lifeforms, it is less suited to evaluating gradual change over time. Further, stronger statistical links can be drawn by examining relationships between the abundances of distinct planktonic lifeforms and environmental pressures, than between an index of lifeform pairs and environmental pressures (Bedford et al., 2020a). For the current assessment a simplified approach was developed, based on measuring long-term changes in the abundance of key lifeforms and evaluating how environmental pressures co-varied with changes in lifeform abundance. A long-term perspective helped improve confidence in establishing links to environmental pressures. Short-term changes within the assessment period were also evaluated to identify emerging trends. The PH1 indicator as originally defined in IA2017 was still evaluated for this assessment and is described under the Lifeform pairs indicator approach section of the Assessment Methods and the respective results are described in the Results (extended) section.

Kendall trend test and statistic

Abundance trends for planktonic lifeforms can be evaluated over time by applying the Kendall trend test to annual mean abundance values (Bedford et al., 2020a; Desmit et al., 2020). The test was performed on annual log10 transformed mean abundance values, rather than monthly or seasonal values, to remove the seasonal variation typical of plankton time series data (Figure 1). This nonparametric test generates a statistic which is derived by comparing each value in a time-series with each of the values preceding it. If a latter value is greater than a previous one, the pairwise comparison is assigned a value of 1. If it is lower it is assigned a value of -1, with 0 assigned to cases when there is no difference between values. The sum of the pairwise comparisons for the time-series produces Kendall’s S-statistic. The variance in the S-statistic is used to derive a Z-score with an approximately normal distribution; thus, confidence in this statistic can be assessed with an associated p-value, with p ≤ 0,05 generally accepted as statistically significant change. The sign of the test statistic reveals the direction of the trend, with a positive statistic indicating an increasing trend and a negative statistic indicating a decreasing trend. The magnitude of the statistic is proportional to the strength of the trend. A great benefit of this nonparametric test is that it yields identical results irrespective of the data transformation method and is not sensitive either to gaps in data or to non-linear or irregular trends. It also has the advantage that abundance trends are comparable across datasets, lifeforms, and assessment units.

Time-series plots (Figure 1) indicate variation in lifeform abundances through time for a single assessment unit, with the blue lines displaying monthly variability (thinner line), and annual mean abundance (thicker line) values used to derive the Kendall statistic. Data correspond to the assessment unit representing the western entrance to the English Channel (‘Channel well mixed’ see Figure a) and contains the fixed-point station ‘L4’, consistently monitored by the Plymouth Marine Laboratory since 1988.

The ability to aggregate complex plankton data into ecologically relevant lifeforms is a useful tool for assessment because it can help identify where changes are occurring (McQuatters-Gollop et al., 2019). The ability to understand which taxa are driving changes in a lifeform can also be valuable. By calculating the Kendall statistic for each taxon, additional information can be gleaned to assist in interpreting lifeform time-series trends. For example, Figure a displays change in the abundance of diatoms for the ‘Channel well mixed’ COMP4 assessment unit described in Figure 1. The abundance of diatoms has been increasing since 1960 (z = 2,66, p ≤ 0,05). Through examining the annual abundance time-series for individual diatom taxa it becomes possible to also conclude that overall trend of increasing diatom abundance is driven by positive abundance trends in 11 out of 54 diatom taxa. Increasing trends and greater absolute abundances for Proboscia alata and Thalassiosira spp indicate that they are predominantly responsible for the increasing abundance trend in diatoms. Additionally, while the net trend is increasing, four diatom taxa (Guinardia striata, Skeletonema costatum, Thalassionema nitzschioides, and Odontella sinensis) are actually decreasing in abundance.

Since this assessment was based on data from multiple sources, the Kendall statistic was calculated independently for each combination of dataset, assessment unit, and lifeform. Since the datasets used also varied in duration, to avoid discarding data the full duration of each dataset was assessed up until the end of 2019. For a dataset to be included in this analysis, it needed to contain samples collected within the assessment period (2015-2019) and prior to the assessment period.

Spatial Scale:

Because plankton community composition, distribution, and dynamics are closely linked to their environment, the analysis was performed at the scale of the ‘COMP4 assessment units’ (COMP4 v8a; Figure b, Table a). Assessment units within the Greater North Sea and Celtic Seas (OSPAR Regions II and III, respectively) were initially developed by Deltares and partner institutes as part of the EU ‘Joint Monitoring Programme of the Eutrophication of the North Sea with Satellite data’ (JMP-EUNOSAT; Enserink et al., 2019) and further refined in the revision process of the eutrophication assessment by OSPAR expert groups ICG-EMO and TG-COMP. Assessment units with similar phytoplankton dynamics were derived from cluster analysis of satellite data for chlorophyll a and primary production. Boundaries between assessment units were derived by relating clustering results to the best-matching gradients in environmental variables obtained from the three-dimensional hydrodynamic Dutch Continental Shelf model version 6 (DCSMv6 FM). The variables which best matched the divisions highlighted by clustering were depth, salinity, and stratification regime. Additional geographic areas were added such as the Channel, Irish Sea, and Kattegat. These assessment units are a geographical representation of the conditions most likely to drive plankton distribution, dynamics, and community composition. The assessment units are being regularly updated and may not perfectly reflect the true state of the pelagic environment in their present form.

Because the Bay of Biscay and Iberian Coast (OSPAR Region IV) extended beyond the boundaries of the DCSMv6 FM, assessment units within this Region were developed using a different methodology, based on phytoplankton dynamics (Spain) and salinity dynamics (Portugal). To delineate assessment units for the Spanish coast, a polygon was created to extend from the coast to the exclusive economic zone (EEZ) boundary. Daily NASA Aqua/MODIS Level-2 satellite images (https://doi.org/10.5067/AQUA/MODIS/L2/OC/2018) were used to calculate climatological mean values of chlorophyll a for each pixel. K-means clustering was then used to group pixels with similar dynamics, resulting in six distinct groupings within the main Spanish polygon. Portugal’s three Water Framework Directive assessment units were extended to the boundaries of the Portuguese EEZ. These assessment units were further divided longitudinally to separate pelagic waters from coastal waters more subject to eutrophication from river influence by applying a salinity threshold of 30,0–34,5 for coastal waters and >34,5 for offshore waters, followed by a bathymetry threshold of 100 m.

COMP4 assessment units and fixed-point stations were categorised according to the OSPAR Regions they intersected as well as the habitat type they represented (Table a, Table b, Figure c). The four habitat types considered were: variable salinity (representing areas of freshwater influence where estuarine plumes extend beyond waters designated as Transitional Waters under Directive 2000/60/EC), coastal (representing areas with mean salinity < 34,5 psu), shelf (representing areas with mean salinity > 34,5 and mean depth < 200 m) and oceanic / beyond shelf (representing areas with mean salinity > 34,5 and mean depth > 200 m).

Plankton Data:

The assessment has been carried out using 24 phytoplankton and zooplankton datasets from 14 sources (Table c, Figure d). Other datasets were provided but they were out of the scope of this assessment (e.g., time-series with comparison period duration shorter than assessment period duration, time-series ending prior to the assessment period, time-series commencing during assessment period, distributed sampling locations restricted to transitional waters only). These additional datasets were submitted by: Vlaams Instituut voor de Zee (BE), Bundesamt für Seeschifffahrt und Hydrographie (DE), Centre for Environment, Fisheries and Aquaculture Science (UK), Newcastle University (UK), and Natural Resources Wales (UK).

The data submitted by AU (DK), PML (UK), MSS (UK), NLWKN (DE) and IEO (ES) were from discrete fixed-point stations which were evaluated independently. The data submitted by SAMS (UK) were from multiple stations within a small area of freshwater influence. These data were aggregated and treated as a single fixed-point time series.

Data from the Continuous Plankton Recorder (CPR) survey, collected by the Marine Biological Association (MBA, UK), consisted of spatially distributed data collected along transects. The CPR survey, provides offshore open ocean data at a broad spatial scale using ships-of-opportunity (Reid et al., 2003). Due to the distributed nature of CPR data, CPR samples are typically aggregated across a grid or a set of polygons at monthly temporal resolution (Bedford et al., 2020a). The remaining fixed-point datasets, including VLIZ (BE), LLUR (DE), RWS (NL), IPMA (PT), SMHI (SE), Cefas (UK), and EA (UK), were also treated in this manner since samples were collected at several distributed stations within close proximity to one another. The dataset from IPMA (PT) only records abundances for the genera Pseudo-nitzschia (diatom) and Dinophysis (dinoflagellate), thus results from the Portuguese data should only be considered as a proxy for diatom and dinoflagellate lifeform abundances.

Data from the different providers were not combined for analysis due to differences in sampling, plankton analysis and enumeration methods. Instead, the datasets were analysed separately. Each dataset has internal QA/QC procedures to ensure consistency and accuracy of the data. Before total lifeform abundance values were log10 transformed, a nominal value equivalent to half the minimum non-zero observed value for each time-series was added to each sample. All spatially distributed data (e.g., excluding AU, PML, MSS, IEO, NLWKN and SAMS) were averaged per month within each assessment unit. In a few cases, this resulted in localised groupings of samples being extrapolated across much larger assessment units (see distribution of CPR samples in Figure d (i)), as was the case for CPR data along the west of Scotland (Intermittently Stratified 1) and Ireland (Atlantic Seasonally Stratified). For months when there were no samples within an assessment unit, gaps in the time series were filled by extracting a mean value from an inverse distance weighted interpolated surface generated from samples adjacent to the assessment unit (<250 nmi) from the same dataset. Years containing less than eight months of averaged and spatially interpolated values were discarded. Remaining gaps of three months or less were filled via linear interpolation. For the reporting of trends and linking to pressures, in cases where multiple suitable datasets contained samples from within the same assessment unit the trends displayed were based on the dataset with the greatest number of unique months (e.g., June 1978, February 2000, May 2014).

Confidence scoring:

A confidence scoring methodology, based on an approach developed by ICG-Eut to validate the output from their COMPEAT Tool, was applied to evaluate the robustness of reported trends for each plankton dataset for each assessment unit it intersected. For each assessment unit or fixed-point station temporal confidence was evaluated using two metrics, which are referred to here as specific temporal confidence (STC) and general temporal confidence (GTC).

STC was calculated on an annual basis as the proportion of months when samples were reported out of a total of 12 possible months. For example, if a dataset reported samples within an assessment unit only for June, July and August in a particular year, the STC for that year would be 0,25. GTC was calculated as the proportion of unique months with recorded samples out of the total number of potential months across the time-series. For example, a dataset commencing in January 2015 and ending in December 2019 contains a total of 60 potential months. If samples were reported for 30 of those months, the GTC would be 0,5. An overall temporal confidence score was then determined by calculating both the mean STC across all years and subsequently the mean of this value with the GTC.

Spatial confidence was also evaluated for distributed datasets such as the CPR. However, for fixed-point stations spatial confidence was not considered since there was no spatial element to the data. For each COMP4 assessment unit spatial confidence was calculated using two metrics, which similarly are referred to here as specific spatial confidence (SSC) and general spatial confidence (GSC).

SSC was calculated on an annual basis based on the proportion of longitudinal and latitudinal ranges of an assessment unit covered by the spatial extent of samples. For example, if an assessment unit spanned longitude from 1 to 5° E, and the distribution of samples collected within that assessment unit within a particular year ranged from 2 to 4° E, the proportion of longitudinal range covered by samples would be 0.5. The same procedure is also applied to the latitudinal range. The SSC for a particular year would be the mean of the longitudinal and latitudinal proportions covered by the spatial extent of samples.

GSC was derived using the same approach as SSC, but across all samples in a dataset that intersect with an assessment unit. An overall spatial confidence score was then determined by calculating the mean SSC across all years and subsequently calculating the mean of this value with the GSC.

Finally, an overall confidence score was calculated as the mean of temporal and spatial confidence. For fixed-point stations, spatial confidence was not considered, therefore only temporal confidence was used to represent the overall confidence of fixed-point station data.

Lifeform construction:

The eight plankton lifeforms highlighted in this assessment were chosen due to their ecological relevance and owing to the high confidence in their classification (Table d). All lifeforms were defined based on common ecological traits. The rationale for selecting the lifeforms and additional criteria containing supplementary information on lifeforms is listed in Table e.

The database master species list (Ostle et al., 2021) was built by assigning functional traits to each species for a dataset and then adding additional new datasets to expand the master species list (Figure e). The species list for each new dataset was assigned a unique Aphia ID via the World Register of Marine Species (WoRMS). The new dataset’s species list was then compared with the plankton database’s master species list via Aphia IDs and any new species were identified. This process ensures that each species is only entered in the database once. The new species were then manually assigned functional traits by searching the literature; fields were left blank where functional traits for species were unknown. Once traits were assigned, the new species were added to the master species list. Queries were constructed (Table e) to build lifeforms from the functional trait information (McQuatters-Gollop et al., 2019).

A simple method of confidence assessment was applied for each lifeform (Table f). Using expert opinion, each lifeform was evaluated on two characteristics: the ability to identify and speciate organisms in that lifeform using light microscopy and the understanding of the accuracy of determining traits assigned to the lifeform. For example, medium confidence is assigned to the lifeforms ‘autotrophic dinoflagellates’ and ‘mixotrophic dinoflagellates’ because the mode of nutrition of many dinoflagellate species is currently uncertain (Flynn et al., 2019). Thus, the accuracy of assigning the lifeform category is medium. Likewise, the lifeform ‘non-carnivorous zooplankton’ has a medium confidence designation since the feeding habits of many abundant and common zooplankton species remain not strictly defined. Further investigation must also be conducted to decide whether both harmful-bloom-forming algae and potentially toxin-producing lifeforms should be considered in future assessments of this indicator. Work is ongoing to increase confidence in lifeforms, however, this work is resource dependent.

Assessing links to pressures:

The Boruta algorithm (Kursa and Rudnicki, 2010) was applied to evaluate which environmental variables were the best predictors of lifeform abundance. The Boruta algorithm is a tree-based permutational variable importance method which uses a wrapper around random forest to evaluate the predictive performance of each variable being tested against a shuffled set of the same predictor variables. If a variable achieves better predictive accuracy than the highest-performing shuffled variable, it is determined to be important and is assigned a score based on the mean decrease in predictive accuracy when the variable is excluded from the model.

For each unique combination of dataset, assessment unit, and lifeform whenever a significant Kendall statistic (i.e., a significant change in lifeform abundance over the time-period assessed) was observed, a separate permutational variable importance process based on random forest (Boruta algorithm) was conducted to evaluate relationships between lifeform abundance and a set of 16 environmental variables described in Table g. Random forest is robust to detecting separate effects from collinear predictors, since it uses an ensemble approach based on agreement from many trees which each test a different subset of the predictor variables. Based on the results of the previous assessment, a table of implicating factors was developed, linking lifeforms with a set of environmental pressures with plausible links to variation in lifeform abundance (Table h). Only variables indicated in Table h were tested for each lifeform.

Table g: Environmental variables and their descriptions and sources evaluated in the random forest permutational variable importance process.

Prior to any manipulation of data, lifeform time-series were divided into separate training and testing sets. The training data were used for variable selection and for generating the random forest models, while the testing data were used to validate model predictive accuracy. For the purposes of this analysis, the training data were limited to all months prior to the assessment period (e.g., 1960 to 2014 for CPR data) and the testing data were restricted to the assessment period itself (e.g., 2015 to 2019).

This analysis covered the full temporal extent of each plankton dataset (1960 to 2019 in the case of CPR data). While gridded data from NOAA were available from 1960 to 2019, the longest duration modelled dataset for nutrients only spanned from 1993 to 2019, with in-situ nutrient datasets commencing as early as 1980. To evaluate long-term links to pressures and avoid excluding the first several decades of many plankton time-series from our analysis the ‘missRanger’ package for R (Mayer and Mayer, 2019) was used to perform multiple imputation by chained random forests. This method uses multiple random forest regressions to impute missing values based on collinearities among observed values in the predictors. For each variable containing missing values the algorithm generates a separate regression model based on all the other predictors. To improve imputation performance, a numeric variable representing ‘month’ was included in this step to better predict the consistent seasonal patterns in some variables. Missing values in the predictors were imputed separately for the training and testing datasets. It is important to note that no imputation of lifeform abundance time-series was performed to avoid introducing greater uncertainty into the process of linking to pressures.

All variables including lifeform abundance were smoothed to remove seasonality and uncover long-term trends by calculating the mean value across a 12 month moving window (Figure f). This step required the first and last six months of data from training and testing sets to be excluded from the analysis.

Values for each variable were calculated as the mean of monthly mean gridded values (modelled and remotely sensed; Table g) intersecting each COMP4 assessment unit. For fixed-point stations, mean values were calculated from all measurements within a 5-nautical mile radius of the station. Where in-situ data were available (total nitrogen, nitrate, phosphate, total phosphorous, silicate) they were evaluated instead of the modelled environmental variables. For Atlantic Multidecadal Oscillation (AMO) and North Atlantic Oscillation (NAO), monthly values were applied identically across all assessment units since these variables have basin-scale influence likely to cover most of the assessment region.

Environmental variables were grouped into ten categories to simplify visualisation in map format (Table h). Categorisations of the most important predictor variable for each model, or the predictor variable which resulted in the greatest mean decrease in accuracy when it was removed, were reported across a map. Instances when the validation step with the testing data indicated significant positive correlation (p ≤ 0,05) between predicted and observed lifeform abundance values were also indicated in this map to provide greater confidence in the variable importance results and to avoid reporting spurious relationships.

It is important to note that observing high importance of any environmental variable is simply an indication that it co-varied with lifeform abundance. Our approach cannot imply causation since it is impossible to hold all other potential influences on lifeform abundance constant, as would be the case in a controlled experimental setting.

Integration of indicator results:

A primary objective of this indicator assessment was to integrate results to facilitate an understanding of changes occurring across pelagic habitat types within OSPAR Regions II, III and IV. This required indicator results for each OSPAR Region to be integrated according to the following pelagic habitat categories: variable salinity, coastal, shelf, and oceanic / beyond shelf. This categorisation of COMP4 assessment units and fixed-point stations is described in Table a and Table b. To meet this objective, the focus was on the primary direction of change detected across assessment units and fixed-point stations within each pelagic habitat category for each of the 8 high confidence plankton lifeforms highlighted in this assessment. The mean confidence, spatial representativeness were then reported, and most likely links to environmental pressures.

As an example, changes in dinoflagellate abundance were assessed across 4 COMP4 assessment units and 4 fixed-point stations representing variable salinity habitats within OSPAR Region II. If 1 decreasing trend, 4 increasing trends, and 3 instances of no trend were detected across these locations, an increasing net trend would be reported and the proportion of assessment units studied where this trend was detected, in this case 0,5.

To report the level of confidence in this result, the mean confidence score for locations where dinoflagellate abundance was increasing, was calculated. COMP4 assessment units and fixed-point stations as equivalent were considered for this integration.

To report the spatial representativeness of the result, the proportion of the total number of COMP4 assessment units considered in the analysis was calculated, in this case 8, out of the total number of possible COMP4 assessment units representing variable salinity habitats within the OSPAR Region, in this case 16. Therefore, the spatial representativeness of the result would be 0,5. Note that fixed-point station datasets do not contribute to this score.

Finally, to report links to environmental pressures which can drive changes in lifeform abundance for the net trend, environmental variables for each location based on their relative variable importance were ranked, with 1 assigned to the variable with highest importance, 2 to the variable with second highest importance and so on. For locations where the net trend was increasing, the mean rank of each environmental variable was calculated and the variable with the lowest mean rank reported.

Assessing the status of pelagic habitats:

To assign a designation of assessment status to pelagic habitats based on the integration of indicator results, a semi-quantitative methodology described in McQuatters-Gollop et al., (2022) was applied, which was developed from the lessons gained from the previous OSPAR assessment (IA 2017). In this case, the status of a habitat type can be designated as either “Good”, “Unknown”, “Not good”, or “Not assessed” (Table i). Following the criteria outlined in this study, if a pelagic habitat has been assessed, it should by default be considered as either “Unknown” or “Not Good”. At this stage it is not realistic to assign “Good” status to pelagic habitats, since it is difficult to develop meaningful assessment thresholds for plankton and generally not possible to determine whether a particular state is desirable or undesirable, except under specific circumstances such as eutrophication. Following this logic, the status of pelagic habitats should be considered “Unknown” by default. In cases when change has been detected and that change can be confidently linked to the impact of an anthropogenic pressure, the status of this habitat is “Not good”.

Table i: Biodiversity status categories and colours used for the interpretation, by expert judgement, of indicator biodiversity state (McQuatters-Gollop et al., 2022).

The status at the level of each of the 8 high confidence lifeforms within each pelagic habitat within OSPAR Regions was considered and the results integrated for multiple lifeforms to assign a single quality status designation for the pelagic habitat type. For the status of a lifeform to be shifted from “Unknown” to “Not good” the results of the integration had to meet certain criteria:

• The net trend should be either increasing or decreasing and must be present in at least 0,5 of the locations assessed.
• The mean confidence for locations considered for the determination of the net trend must be at least 0,5.
• The spatial representativeness of the assessed locations out of the total number of possible locations for that habitat type must be at least 0,5.
• The environmental pressure with the greatest mean rank for locations expressing the net trend must represent either sea surface temperature, pH, or nutrients.
• The mean rank of the most important environmental pressure must be ≤ 3, indicating that across all locations the variable ranks in the top 3 most important for predicting the abundance of the lifeform.

If all the above criteria are met, the lifeform is assigned a status of “Not good”. If 25% or more of assessed lifeforms within a pelagic habitat type are assigned a status of “Not good” then the whole habitat type is also assigned a status of “Not good”. Otherwise, the habitat type is assigned a status of “Unknown”.

Lifeform pairs indicator approach:

Although this assessment focused primarily on long-term trends in lifeform abundance due to their stronger links with environmental pressures at large spatial scales, the lifeform pairs indicator was also evaluated for the current assessment (which constitutes the PH1 indicator as originally defined in IA2017).

Tett et al., (2008) proposed to track changes in the state of the plankton community by means of plots in a state space and calculating a Plankton Community Index, referred to here as a Plankton Index (PI). The conceptual framework is that ecosystems can be viewed as systems with an instantaneous state defined by values of a set of system state variables which are attributes of the system that change with time in response to each other and external conditions. Building on this approach and plotting plankton lifeform abundance in a multi-dimensional state space provides a means of monitoring changes in the structure of plankton communities. A state can be defined as a single point in state space, with coordinates provided by the values of the set of state variables, in this case two lifeform abundances, which together are used as a pair.

The plankton community index is calculated by initially mapping relative plankton lifeform abundances from an assessment period in a state space. The distribution of lifeform abundances within the state space is then used to define an envelope, representing the prevailing conditions for the assessment period. Finally, a set of abundance values from a comparison period are mapped onto the state space to evaluate the plankton community index, which represents the proportional similarity of the two periods.  

Mapping lifeforms in a state space

In the example illustrated in Figure g for diatoms and dinoflagellates, the axes of the two-dimensional space are the abundance values of the two lifeforms. Each point represents the state of the ecosystem in terms of the two lifeforms at the time of sample collection. Subsequent samples yield additional pairs of abundance values that can be mapped onto the lifeform state space.

Point A is the ecosystem state at the instant a water sample was taken and is characterised by the abundance values of two lifeforms. Another sample, taken in the same location, yields abundance values that map to a different point in the diatom-dinoflagellate state space (point B).

The path between the two states is called a trajectory, and the plankton condition is defined by the trajectory drawn in the state space by a set of points. Such trajectories reflect: (i) cyclic and medium-term variability (the higher order consistencies in the plankton that result from seasonal cycles, species succession and interannual variability); and (ii) long-term variability that might result from environmental pressure. The seasonal nature of plankton production and the succession of species in seasonally stratifying seas, result in this trajectory tending in a certain direction (as plankton growth increases in spring and declines during autumn), such that the trajectory tends towards its starting point. Given roughly constant external pressures, the data collected from a particular location over a period of years form a cloud of points in state space that can be referred to as a regime. Long-term variability may show a persistent trend of movement away from a starting point in state space.

Approach for defining the envelope

To define a regime, an envelope can be drawn about this group of points using a convex hull method. Because of theoretical arguments that the envelope should be doughnut-shaped with a central hole (Tett et al., 2008), bounding curves can be fitted outside and inside the cloud of points (Figure h). The data are from the CPR dataset for the ‘Channel well mixed’ COMP4 assessment unit. The colour of each point corresponds to the season within which it was sampled; blue = winter, green = spring, yellow = summer, red = autumn.

The size and shape of the envelope are sensitive to sampling frequency and the total number of samples. Envelopes are made larger by including extreme outer or inner points, and the larger the envelope, the less sensitive it will be to change in the distribution of points in state space and therefore to detect a change in condition. Conversely, if too many points are excluded the envelope will be small and even minor changes will result in a statistically significant difference. It is therefore desirable to exclude a proportion (p) of points, to eliminate these extremes, and so the 90-percentile was used. Envelopes are therefore drawn around the cloud of points to include a proportion (p=0,9) of the points, excluding the 5% of points that were most distant from the cloud's centre and the 5% of points that were closest to the centre.

For a Plankton Index to be calculated, it is necessary to establish a set of conditions as the basis for making comparisons. For this assessment the Plankton Index was calculated by using the assessment period (2015-2019) to define the envelope. For each dataset, all measurements collected prior to the assessment period were used for comparison. For the CPR data this meant the envelope was defined with samples from 2015 to 2019 and samples from 1960 to 2014 were used for comparison. This made it possible to maintain the same range of years to define the envelope across multiple datasets with different temporal coverage. Like the Kendall statistic, the PI was calculated independently for each combination of dataset, assessment unit and lifeform pair.

In the example shown in Figure i, CPR data collected from the ‘Channel well mixed’ COMP4 assessment unit during the assessment period (2015-2019), were used to create an envelope. The envelope, thus drawn (Figure i) defines a domain in state space that contains a set of trajectories of the diatom-dinoflagellate component of the pelagic ecosystem.

Calculating the plankton community index

The next step is to map a new set of data onto the same state space for comparison (Figure i). The value of the PI is the proportion of new points that fall between the inner and outer envelopes. In this example, 30% or 198 of the 660 new points lie outside, and the PI is 0,7. A value of 1,0 would indicate no change, and a value of zero would show complete change. The envelope was made by excluding 10% of points, so 66 (10% of 660) points are expected to fall outside. The exact probability of getting 198 by chance alone when only 66 are expected, can be calculated using a chi-square calculation (with 1 df and a 1-tail test). The value of 0,7 is significantly less than the expected value of 0,9, and so the difference between the two periods is statistically significant (p < 0,01).

Statistically significant continuation of long-term trends was only apparent in a few assessment units when the Kendall statistic was calculated for the assessment period only (2015 to 2019; Figure j), owing at least partially to the reduction in statistical power associated with detecting trends across such short time periods. Longer assessment periods increase the detection of many of the important changes which are occurring in plankton lifeforms, due to high inter-annual variability in lifeform abundance and gradual rates of change over time. Increasing the temporal extent also increases the detection of change because data from more historic periods are included (Bedford et al., 2020a; Bedford et al., 2020b)

The lifeform pairs indicator (PI), which indicates the degree of similarity in the co-abundance of lifeform pairs (Figure k) between the assessment period (2015 to 2019) and the comparison period (1960 to 2014) reveals that changes in diatoms and dinoflagellates were greatest in coastal waters surrounding the United Kingdom and within the Norwegian Trench, driven by reduced summer and autumn abundances of dinoflagellates. Diatoms and dinoflagellates showed the most stability within the Channel well mixed, Meuse and Elbe plume assessment units. Fixed-point stations, particularly those around Scotland, did not always reflect patterns in offshore CPR data. For gelatinous zooplankton and fish larvae/eggs there were insufficent records of gelatinous zooplankton to form any conclusions, however there was evidence of important changes in the Norwegian Trench and Kattegat Coastal assessment units. For holoplankton and meroplankton the greatest changes were along the Atlantic on the northern coast of Spain (Noratlantic Area NOR-NorO1) and in the Norwegian Trench. Widespread changes in the Dogger Bank, Eastern and Southern North Sea assessment units were caused by increasing trends in holoplankton winter abundance. Large and small copepods demonstrated a similar pattern to that of holoplankton and meroplankton. The PI also showed poor overlap with patterns in the Kendall statistic for individual lifeform abundances, which indicates there were important changes to seasonal trajectories which did not contribute to detectable trends in mean annual abundance.

The integration of the indicator results across pelagic habitat categories (Figure c) indicate that meroplankton in coastal and shelf areas of the Greater North Sea (Table j) are increasing in abundance and the increasing trends are closely linked to rising sea surface temperatures. Simultaneous declines in the abundance of holopankton in shelf regions have also been linked to temperature increase.

Table j: Integration of the indicator results for OSPAR Region II – Greater North Sea. Column names are described as follows: ↓: the number of COMP4 assessment units and fixed-point stations where decreasing trends have been detected, -: the number of COMP4 assessment units and fixed-point stations where no trends have been detected, ↑: the number of COMP4 assessment units and fixed-point stations where increasing trends have been detected, Dir: the net direction of change in the lifeform (↓: decreasing, ↑: increasing, -: stable), Trend: the percentage of assessment units exhibiting the respective trend, Conf: the mean confidence of datasets considered in the assessment, Change: a logical variable (TRUE/FALSE) to report whether a net trend is likely given the proportion of locations expressing the trend and the confidence and spatial represenativeness scores, Press1: the environmental pressure with the greatest mean rank for the respective trend, Rank1: the mean rank of the environmental pressure indicated under Pres1, nStn: the total number of fixed-point stations considered, totAssess: The total number of COMP4 assessment units and fixed-point stations considered, totCOMP4: The total number of potential COMP4 assessment units for the habitat category, spatialRep: the spatial representativeness score of the analysis. The status of the individual lifeforms are indicated by the colours in the Lifeform column, according to the legend in Table i. The overall status of the habitat category is indicated by the colour of the first column, which also identifies pelagic habitat types and follows the same colour legend.

In the Celtic Seas (Table k), declining abundance of diatoms in variable salinity habitats have been observed and linked to the increase in the N:P ratio. Declining abundance trends in dinoflagellates, as well as large and small copepods in coastal habitats were linked to temperature increases. Declining abundance of holoplankton in coastal areas may also be linked to decreasing pH. Changes in pH can be linked to changes in primary production, since phytoplankton ingest Dissolved Inorganic Carbon (DIC) to fuel growth and reproduction, but also to the concentration of atmospheric carbon dioxide, the main driver of global warming. In shelf areas a similar increase in meroplankton abundance has been observed in response to temperature, while decreasing abundance of holoplankton and dinoflagellates may be linked to reductions in pH and  increases in the N:P ratio, respectively. 

Table k: Integration of the indicator results for OSPAR Region III – Celtic Seas. An explanation of column names can be found in the caption for Table j.

In the Bay of Biscay and Iberian Coast region (Table l), holoplankton, large and small copepods in shelf areas are also exhibiting decreasing abundance trends linked to temperature increase, while large and small copepods in oceanic / beyond shelf habitats are experiencing declines linked to increasing temperatures. Trends for variable salinity habitats were not assessed due to the limited data available for this region.

Table l: Integration of the indicator results for OSPAR Region IV – Bay of Biscay and Iberian Coast. An explanation of column names can be found in the caption for Table j.