Temporal and Spatial Patterns of Harmful Algae Affecting Scottish Shellfish Aquaculture

Gianella, Fatima; Burrows, Michael T.; Swan, Sarah C.; Turner, Andrew D.; Davidson, Keith

doi:10.3389/fmars.2021.785174

ORIGINAL RESEARCH article

Front. Mar. Sci., 22 December 2021

Sec. Marine Fisheries, Aquaculture and Living Resources

Volume 8 - 2021 | https://doi.org/10.3389/fmars.2021.785174

This article is part of the Research Topic Current Challenges in Providing Early Warning of Harmful Algal and Microbiological Risk to Aquaculture View all 13 articles

Temporal and Spatial Patterns of Harmful Algae Affecting Scottish Shellfish Aquaculture

$\r\nFatima Gianella*$ Fatima Gianella^1*

Michael T. Burrows¹

Sarah C. Swan¹

Andrew D. Turner² $Keith Davidson\r\n$ Keith Davidson¹

¹Scottish Association for Marine Science, Scottish Marine Institute, Oban, United Kingdom
²Centre for Environment, Fisheries and Aquaculture Science, Weymouth, United Kingdom

Consistent patterns of Harmful Algal Bloom (HAB) events are not evident across the scientific literature, suggesting that local or regional variability is likely to be important in modulating any overall trend. This study summarizes Scotland-wide temporal and spatial patterns in a robust 15-year high temporal frequency time series (2006–2020) of the incidence of HABs and shellfish biotoxins in blue Mussels (Mytilus edulis), collected as part of the Food Standards Scotland (FSS) regulatory monitoring program. The relationship between the countrywide annual incidence of HAB events and biotoxins with environmental variables was also explored. Temporal patterns exhibited interannual variability, with no year-on-year increase, nor any correlation between annual occurrences. Within years, there was a summer increase in bloom frequency, peaking in July for Dinophysis spp. and Pseudo-nitzschia spp., and a plateau from May to July for Alexandrium spp. Temporal-spatial patterns were analyzed with multivariate statistics on data from monitoring sites aggregated monthly into 50-km grid cells, using Principal Component Analysis (PCA) and cluster K-means analysis. PCA analyses showed correlation between areas with similar temporal dynamics, identifying seasonality as one of the main elements of HAB variability with temporal-spatial patterns being explained by the first and second principal components. Similar patterns among regions in timing and magnitude of blooms were evaluated using K-means clusters. The analysis confirmed that the highest risk from HABs generally occurred during summer, but demonstrated that areas that respond in a similar manner (high or low risk) are not always geographically close. For example, the occurrence of the most prevalent HAB genus, Dinophysis spp., is similar countrywide, but there is a regional trend in risk level with “very-high” and “high” clusters located primarily on the southwest coast, the islands of the central and northern west coast and the Shetland Islands. “Early” and “late” blooms were also associated with certain areas and level of risk. Overall, high risk areas mainly face in a southwest direction, whilst low risk locations face a south or southeast direction. We found relatively few countrywide relationships between environmental variables and HABs, confirming the need for regional analysis to support HAB early warning.

Introduction

Harmful algal blooms (HABs), associated with human and shellfish toxicity, are temporally and spatially variable. The occurrence, intensity and distribution of HABs is a problem worldwide (Hallegraeff, 1993, 2010; Van Dolah, 2000; Smayda, 2002; Glibert et al., 2005; Anderson et al., 2017; Gobler et al., 2017; Wells et al., 2019). However, trends of increasing HAB frequency and/or abundance are not evident in all studies (Moore et al., 2009; Díaz et al., 2016; Dees et al., 2017) suggesting that local or regional variability is likely in any overall trend. This is consistent with recent results from Hallegraeff et al. (2021) who used a meta-analysis of Harmful Algae Event Database and Ocean Biodiversity Information system data to demonstrate that there is no empirical support for a global increase in HAB events. Understanding the temporal-spatial variability of HABs and environmental drivers at a regional level is therefore of considerable importance to producers and consumers of shellfish.

The presence of toxic phytoplankton and synthesized biotoxins is a threat for shellfish aquaculture, since the consumption of shellfish that have concentrated these toxins through ingestion of harmful phytoplankton impacts negatively on human health (Smayda, 1990; Berdalet et al., 2016). The economic consequences of HABs for shellfish aquaculture have recently been shown to be significant (Mardones et al., 2020), with Martino et al. (2020) estimating that the incidence of biotoxins associated with Diarrheic Shellfish Poisoning (DSP) causes an annual average loss of 15% in mussel production in Scottish waters, equivalent to £1.37 m.

Scottish waters are impacted by a number of HAB genera, but concern over shellfish toxicity is primarily related to the genera Dinophysis, Alexandrium and Pseudo-nitzschia (Davidson and Bresnan, 2009). Dinophysis spp., responsible for the production of Diarrheic Shellfish Toxins (DSTs) by the synthesis of okadaic acid, dinophysistoxins and equivalents (Reguera et al., 2014) have caused prolonged shellfish harvesting closures in Scotland (Whyte et al., 2014; Swan et al., 2018; Bresnan et al., 2021). The toxic Dinophysis species most commonly associated with harmful events in Scotland are D. acuminata and D. acuta (Swan et al., 2018). Some Alexandrium spp. cause Paralytic Shellfish Poisoning (PSP) through the synthesis of saxitoxins (Brown et al., 2010; Touzet et al., 2010). The most common toxin producer in Scottish waters is thought to be Alexandrium catenella (Bresnan et al., 2021) with A. minutum and A. ostenfeldii also having been reported (Brown et al., 2010; Lewis et al., 2018). Non-toxic A. tamarense (Touzet et al., 2010) and A. tamutum (Brown et al., 2010) are also present. Pseudo-nitzschia spp. are linked to Amnesic Shellfish Poisoning (ASP) through the production of domoic acid (DA) (Fehling et al., 2006, 2004; Bresnan et al., 2017; Rowland-Pilgrim et al., 2019). In Scottish waters, toxic and non-toxic members are found in the Pseudo-nitzschia delicatissima group, with toxic species include P. calliantha, P. pseudodelicatissima and P. cf. pseudodelicatissima (Fehling et al., 2006; Smayda, 2006). The P. seriata group includes species that synthesize high levels of DA, P. australis, P. seriata, P. fraudulenta, P. multiseries (found occasionally), P. pungens and P.cf. subpacifica (Fehling et al., 2006; Smayda, 2006). Azaspiracids are infrequently recorded, and as their causative phytoplankton taxa (Azadinium/Amphidoma) are not enumerated within the regulatory monitoring program they have not been included in our analyses. The national monitoring program enumerates and reports the main harmful taxa to genus level only. While it is recognized that both toxic and non-toxic species may therefore be enumerated, this approach is used to (a) ensure rapid sample turn around and (b) on a precautionary basis to best protect human health. Hence, the bulk of the analysis and discussion is based on taxonomic resolution to genus level. However, where appropriate, based on other published studies, we comment on how these trends may relate to particular harmful species and their toxins.

Trends observed in Scotland and other site specific studies suggest HABs are spatially and temporally variable, with blooms likely linked to large scale oceanic variations (Belgrano et al., 1999), meteorological or oceanographic anomalies (Moita et al., 2016), and specific environmental conditions, rather than a steady increase every year. Studies involving atmospheric dynamics in the Atlantic Ocean use the North Atlantic Oscillation (NAO) index as the difference between the low pressure system around Iceland and the high pressure of the Azores Islands in the mid North Atlantic Ocean. A positive index indicates a stronger wind intensity predominantly from the southwestern direction, whilst a negative NAO index is related to a smaller difference between the systems leading to lower wind speeds (Phillips et al., 2013). These atmospheric patterns influence the climate and consequently environmental conditions in the water column. Several studies have suggested the NAO index is a useful predictor for the development of HAB events, specifically Dinophysis (Belgrano et al., 1999; Báez et al., 2014; Ruiz-Villarreal et al., 2016). Changes in the prevalent wind direction have, through their impact on surface oceanographic currents, led to the development of HAB events in the Shetland Islands (Whyte et al., 2014). In other locations, changes in wind patterns are associated with upwelling/downwelling conditions that cause HAB events, as seen in the northwest coast of Spain (Fraga et al., 1988) and Bay of Bantry, Ireland (Cusack et al., 2016). An increase in sea surface temperature (SST) has also been associated with increased growth rate of phytoplankton, including the toxin-producing taxa Dinophysis and Alexandrium, found in Scotland (Peperzak, 2003; Gobler et al., 2017; Wells et al., 2019). In some locations, the increase of temperature above a certain limit is the main requirement for blooms of toxic species, such as Alexandrium minutum in the Bay of Brest (Chapelle et al., 2015).

Coastal Dinophysis and Pseudo-nitzschia events have been associated with advective oceanographic transport of blooms that initially develop offshore (Fehling et al., 2012; Whyte et al., 2014; Paterson et al., 2017). Cyst forming Alexandrium are more likely to be controlled by local hydrodynamics (Bresnan et al., 2021), with Alexandrium blooms potentially related to the germination of cysts introduced into the sediment in previous years (Anderson et al., 2005; Brown et al., 2010). Characterizing the spatial and temporal trends of different HAB genera over multiple years is therefore not straightforward. Many authors have attempted to understand HAB dynamics, using linear modeling in relation to single or multiple environmental drivers (Hinder et al., 2012; Ruiz-Villarreal et al., 2016; Dees et al., 2017), but such approaches are best suited to the evaluation of decadal trends. Aquaculture practitioners and coastal zone managers have a more pressing need to better understand the likely timing and location of HAB taxa and biotoxins to allow more informed farm management and regulatory decision-making.

High temporal resolution (typically weekly) regulatory monitoring of HABs and their associated shellfish biotoxins in Scottish waters has been undertaken by Food Standards Scotland (FSS) at 143 sites since 2006, generating a unique spatial and temporal dataset of 256 920 observations. Using these data as a regional case study, here we demonstrate the use of statistical techniques to describe, analyze and identify the main HAB and biotoxin patterns in Scottish waters. The relationship between the frequency of HAB events and environmental variables was also assessed. Understanding the dynamics of HABs is useful to support predictive models and other risk assessment approaches to provide early warning of shellfish toxicity (Davidson et al., 2016, 2021). This has the potential to improve measures for safeguarding shellfish consumption and minimize negative impacts in the aquaculture sector.

The statistical approaches used to analyze temporal patterns include the autocorrelation of the frequency of events in the time series. Autocorrelation (or serial correlation) is defined as the correlation of a variable with itself at different times (lag) to examine how current observations relate to those in the past (Cowpertwait and Metcalfe, 2008). Partial autocorrelation is the association with a particular lag, excluding the indirect effect explained by earlier lags. Autocorrelation has previously been used to assess the influence of environmental variables in a HAB time series, with Fischer et al. (2018) finding that cooling temperatures of previous years influenced Alexandrium cyst hatching in the Gulf of Maine. In contrast, the abundance of Pseudo-nitzschia and the autocorrelation within years presented a non-significant trend from 2008 to 2018 in Scotland (Rowland-Pilgrim et al., 2019), and a study evaluating the toxicity of PSP in Puget Sound did not find a significant correlation between toxicity and year from 1993 to 2007 (Moore et al., 2009).

The clustering analyses uses a descriptive approach to assess similarities between samples and identify patterns. K-means analysis is a method that groups a data set into a specified number of clusters (k) (Macqueen, 1967). Samples associated with a particular mean group together, hence the intra-cluster variation is reduced, differentiating these measurements from those in clusters. This approach has been used in the analysis of phytoplankton patterns and to identify the environmental conditions related to their assemblages and composition (Herrera and Escribano, 2006; Marchese et al., 2019; Barth et al., 2020). It has also been used for assessing clusters of toxin profiles from shellfish samples across geographical areas in Great Britain (Turner et al., 2014). Principal Component Analysis (PCA) is a tool for reducing high-dimensional data to fewer dimensions that explain variance and identify the key variables and their role in this variance (Lever et al., 2017). This method has been used to identify phytoplankton temporal patterns, assess the similarity between annual abundances in the time series, links with environmental variables, and similarity of conditions between sampled stations (Solic et al., 1997; Philippart et al., 2000; Kane, 2011; Fehling et al., 2012; Siemering et al., 2016).

This study therefore sought to examine and describe the temporal patterns of shellfish biotoxin-producing HABs and their biotoxins in Scotland. The spatial patterns and the role of seasonality were analyzed for the phytoplankton taxa only, as the suspension of biotoxin sampling and analysis during periods of shellfish farm closure resulted in missing data and inconsistent spatial patterns. The relationship between the annual HABs frequency, biotoxins and environmental variables was evaluated to identify whether there are countrywide relationships and if these can be used as predictors of a toxin event.

Materials and Methods

Data Access and Preparation

Phytoplankton and shellfish biotoxin data were obtained from the Food Standards Scotland (FSS) monitoring program¹. The FSS scheme is operated by the Centre for Environment, Fisheries and Aquaculture Science (CEFAS) for biotoxins, and the Scottish Association for Marine Science (SAMS) for harmful phytoplankton identification. Samples are collected from representative monitoring points (RMPs), primarily located on the west coast and islands (Figure 1) where aquaculture is active. Collection of seawater for phytoplankton analysis is carried out on a weekly basis from March to September and fortnightly in October at all sites, with monthly samples collected from November to February from a reduced number of sites. Harmful phytoplankton genera are identified and enumerated by light microscopy at SAMS. Further details of the program are found in the most recent annual program report (Coates et al., 2020). The following analysis comprises a 15-year time series database, from January 2006 to December 2020.

FIGURE 1

Figure 1. Names and location of the 25 aggregation boxes used (Geographical order North to South).

Phytoplankton samples were collected by officers operating on behalf of several contractors appointed by FSS. This follows the UK National Reference Laboratory (UKNRL) Standard Operating Procedure for the collection of water samples for toxic phytoplankton analysis. These are taken as close to the shellfish bed as possible and at the same location from where shellfish samples for tissue analysis are collected. The sampling method used depends on the depth of water at the site, water samples are usually collected with a 10 m “Lund tube” but occasionally a bucket is used in shallow waters. A well-mixed 500 mL sub-sample of this water is preserved using Lugol’s iodine prior to analysis by light microscopy (Uttermöhl method) at SAMS.

Shellfish were collected using UKNRL guidance and shipped to the CEFAS Weymouth laboratory for biotoxin analysis (Coates et al., 2020). These data were available from March 2007 to December 2020. Monitored shellfish species include common mussels (Mytilus edulis), common cockles (Cerastoderma edule), native oysters (Ostrea edulis), Pacific oysters (Magallana gigas) and razor clams (Ensis spp.). Common mussels constitute ∼62.2% of the total samples and ∼87.7% of samples with total toxin concentrations above maximum permitted level (MPL). As the diverse mechanisms of assimilation and depuration of biotoxins vary according to the shellfish species, our analysis exclusively selected common mussels as the study species due to the high proportion of samples that exhibited toxicity.

Biotoxins levels were quantified by the use of liquid chromatography with tandem mass spectrometry (LC-MS/MS) for lipophilic toxins (LT), including DSTs (Dhanji-Rapkova et al., 2018), and High Performance Liquid Chromatography (HPLC) with either ultra-violet (domoic acid-DA) (Rowland-Pilgrim et al., 2019) or fluorescence detection (Paralytic Shellfish Toxins—PST) (Turner et al., 2014) following official methods specified in EU regulations (Coates et al., 2020). Biological assays were previously used for the purpose of determining shellfish toxicity for PST and DST, with the year the methods changed to chemical instrumentation methods being specific to each biotoxin class. In the case of PST, HPLC was used from 2007 onward, depending on the shellfish species analyzed, with all PST quantitation performed by HPLC by 2011 (Turner et al., 2014). Quantitative determination of DSTs by LC-MS/MS started from July 2011 (Dhanji-Rapkova et al., 2018) replacing the original qualitative biological assay. The biotoxins laboratory quantifies the total toxin concentrations and uses measurement uncertainties of the methods to report “low,” “actual” and “high” values for LT and PST results. FSS uses the latter for risk assessment and as a precautionary measure; this value has been used in the analyses below.

Phytoplankton abundance and biotoxin levels are classified according to the safety threshold value set by the UKNRL following regulation (EC) 2019/627 (Table 1). These thresholds were used to categorize whether a HAB or biotoxin event occurred. The proportional frequency of occurrence of an event above the safety threshold was therefore calculated by dividing the frequency of these events by the total number of samples in a time period or at a specific location. Proportional frequencies of occurrence were averaged to evaluate yearly, monthly and the spatial distribution of HABs and toxic events in Scotland. The proportional frequency of these above regulatory threshold events is used in our analyses below.

TABLE 1

Table 1. Safety threshold limits delimited for phytoplankton species abundance in water samples and maximum permitted level (MPL) for biotoxin in shellfish tissues in Scotland, set by the UKNRL and EC No 853/2004, respectively.

Spatial Range Studied

The data included information from 130 RMPs. Areas located less than 5 km apart were grouped and averaged due to their close proximity and those with fewer than 10 samples over the time series were excluded. Samples were collected at shellfish farms, some of which moved between 2006 and 2019, in which case the median value of their location was calculated.

Temporal and Spatial Statistical Analysis

The statistical analysis was performed using the R program (R Core Team, 2021). Temporal analysis examined the yearly and monthly autocorrelation in the time series using the ACF function in R. This correlated the proportional frequency of events with previous frequencies, at different time lags (x-axis) in yearly and monthly time series. The partial autocorrelation function (PACF) was applied when the autocorrelation function showed significant lag values.

Sampling effort within and between RMPs varied throughout the time series. Hence, spatial analysis and statistics were applied using “aggregation boxes” (n = 25) of 50 km × 50 km to group values after calculating the median (Figure 1). Each box contained between one and 15 RMPs. Boxes recording eight or more missing monthly values in the time series were removed from the analysis. The coordinates, corresponding to the samples and estimated values, were rounded into the whole number to group them in the aggregation boxes.

PCA was applied to identify dominant temporal and spatial patterns of the incidence of HAB events throughout the time series, using monthly values of HABs and toxic events. This analysis also sought to evaluate the role of seasonality (temporal) and spatial patterns in the frequency of HAB events, shown by the first and second principal components (PC). These capture the highest proportion of variance and are the main axes of variation in the data set. The loadings (or eigenvectors) indicate the importance of each variable (month, year, and area) according to the PC. The sample scores refer to the correlation between these variables and the loadings for each PC. PCA cannot be performed on datasets with missing values, so these were estimated from generalized least squares models that included year, month and area averages of the time series.

The spatial risk of the blooms was assessed by clustering the aggregation boxes with similar bloom occurrence, timing and intensity. The K-means analysis identified clusters of high and low risk areas according to the bloom occurrence using the kmeans() function (R Core Team, 2021), with the “elbow” method first being used to determine the most appropriate number of clusters for a particular genus/biotoxin. This was assessed by plotting within clusters sum of squares (WSS) against the number of clusters (k) and evaluating the point at which there was a reduction in the slope (the elbow) using the fviz_nbclust function from the “factoextra” library in R (Kassambara and Mundt, 2020). The proportion of variation and K-means for each k number were calculated and compared for seasonal (monthly) and annual averages series for each species. The analysis sought to explain a high proportion of variance, maximizing similarity within the clusters and minimizing similarity between clusters, with care being taken not to overfit the values. The annual value was calculated using averaged monthly events. The clusters were categorized as “very-high,” “high,” “moderate,” and “low” risk.

The spatial risk of the taxa was evaluated using the monthly cluster risk classification, ranging from 1 to 4 (1 being “low” risk and 4 “very-high”) for quantifying the risk across the taxa and aggregation box. The average risk per aggregation box was calculated and compared to identify the highest and lowest combined risk. The association between HABs risk for all taxa and their location was evaluated using Spearman correlation rcorr() function, which does not assume the data follow a normal distribution.

Relationship With Environmental Variables

Our study explored the relationship between HABs occurrence and environmental variables, including SST, the North Atlantic Oscillation (NAO) index (Supplementary Figure 1) and wind direction and intensity. To assess the strength and the relationships between annual HABs frequency and the variables we used the correlation [rcorr()] function in R. The time series period comprised from 2006 to 2018 for all variables. The SST data was accessed from the Physical Oceanography Distributed Active Archive Center (PODAAC) website² carried out by NOAA National Centers for Environmental Information (PO.DAAC, 2020). Daily NAO index database was downloaded from the NOAA data and averaged (NOAA/National Weather Service, 2021). Wind data were taken from the MIDAS database “Open UK hourly weather observation data” (Met Office, 2019). A total of 39 stations located a maximum distance of 5 km from the coast were used. The average wind speed (knots) was calculated, whilst average wind direction (radians) used the directional component of the wind [sine() and cos()]. Linear regression was used to evaluate the relationship between wind direction and intensity, blooms and toxic events.

Results

Temporal Variability of Harmful Algal Blooms and Toxic Events

Temporal Trends

While all the HAB genera of interest exhibit a typical annual pattern of increased abundance in the spring, summer and autumn months, the frequency of events above the safety threshold on a countrywide scale is highly variable between years with no increasing pattern throughout the time series (Figure 2) for any of the taxa considered. Dinophysis spp. recorded the highest annual frequency of events above the safety threshold (≥100 cells/L) in 2013, reaching 27.53% of events above the safety threshold of the total samples. On a monthly basis, the maximum frequency of events above the safety threshold was recorded in July 2013 for Dinophysis spp. and Alexandrium spp., recording ∼69 and 62% of bloom events, respectively; whilst Pseudo-nitzschia spp. reached a maximum of ∼38% of monthly bloom events in August 2011 (Figure 2). The annual occurrence of blooms of Dinophysis spp. and Alexandrium spp. followed a similar pattern with a strong positive and statistically significant correlation (r = 0.61, p = 0.02). Correlation values between annual bloom events of these genera with Pseudo-nitzschia spp. were not significant.

FIGURE 2

Figure 2. Time series of the proportional frequency of monthly HAB events for genera: Dinophysis spp. DSTs -related (red line), Alexandrium spp. PSTs-related (blue line), Pseudo-nitzschia spp. DA -related (black line) from 2006 to 2020 in Scotland. Total number of samples 16,987.

Yearly frequencies of HABs and biotoxin events for each taxon and toxin group were not correlated with frequencies in preceding years. The only exception is observed for a lag of 1 year for Pseudo-nitzschia spp.

The biotoxins also showed a high variability of monthly occurrence of events above the MPL throughout the time series in a country scale (Figure 3), but in contrast to the phytoplankton, not all the biotoxin groups followed a clear seasonal trend, likely caused by the slow depuration of some toxins. The biotoxins sampling was variable across the country and throughout the time series. 8 545 samples analyzed by mouse bioassay (MBA) for DSTs and 760 for PSTs were included in our analysis. Chemical methods superseded the MBA for the analysis of DST s and PSTs during 2011. For these LC-MS/MS analyses we used 16 679 for DSTs, 19 130 for PSTs and 13 694 for DA. There is a clear difference in the frequency of biotoxin events when determined by MBA (2007–2011) and LC-MS/MS (2011–2018), with the chemical analyses displaying a higher frequency of events. DSTs had the highest frequency of events above the safety threshold, with the highest values during summer months, reaching a maximum monthly average in August 2013 (∼54%). PSTs and DA exhibited sporadic events above the MPL throughout the time series. PSTs exceeded regulatory threshold in ∼10.95% of samples in May 2018. DA only exceeded the maximum permitted limit in seven events, in 4 months, reaching a maximum occurrence (∼1.6%) in August 2007.

FIGURE 3

Figure 3. Time series of the proportional frequency of biotoxin monthly events above safety limits from 2007 to 2020 in Scotland, for DSTs-MBA (dashed red line), PSTs-MBA (dashed blue line) (biological assay); DSTs (red line), PSTs (blue line), and DA (black line) (using instrumental method). Vertical dotted line in July 2011 marks the start of the use of the LC-MS/MS method for DSTs quantification in mussels. In previous years, MBA were mainly used. See Table 1 for definition of abbreviations. Total number of tests performed 58,808.

A comparison between the time series of above regulatory threshold events for monitored HAB taxa and the biotoxins they synthesize does not show a monthly pattern. There is a significant and strong correlation between Dinophysis spp. and DSTs (r = 0.75, p < 0.05), Alexandrium spp. and PSTs (r = 0.60, p < 0.05), whilst Pseudo-nitzschia spp. and DA were not significantly correlated (r = 0.04, p > 0.05).

Seasonal Trends

Countrywide, HAB events above the safety threshold were observed mainly from March to October. HAB occurrence reached maximum monthly median values in July for Dinophysis and Pseudo-nitzschia; Alexandrium reached similar median values in May and June (Figure 4). However, the frequency of events and their timing varied, with both early and late year blooms of all organisms. Early season (March) blooms of both Alexandrium and Pseudo-nitzschia are evident, but in contrast, the average frequency of these early events is nearly zero for Dinophysis although a few exceptional years are observed. For Dinophysis, the median number of events per month increased as the year progressed and reached its highest value, and the highest amongst all taxa, in July before declining toward the end of the year. Alexandrium bloom frequency increased until June, with the frequency remaining high in July before a subsequent decline. Pseudo-nitzschia spp. blooms increased in April, before decreasing in May. Following the July peak, a further increase is evident in September and subsequently the frequency of events above the safety threshold decreases.

FIGURE 4

Figure 4. Boxplots of the monthly incidence of HAB taxa: Dinophysis spp., Alexandrium spp., Pseudo-nitzschia spp., and HAB toxins: DSTs, PSTs, DA. Boxplots display the median (bold horizontal line), upper and lower quartile (hinge), maximum and minimum range (whiskers) and outliers (dots above whiskers) of the yearly average of events by month.

Elevated shellfish toxins are observed mainly in 8 months of the year (Figure 4), with DSTs present at a noticeably higher occurrence in comparison with the other groups. The earliest occurrence of toxicity for DSTs is mainly in June, although there were exceptional toxic events in previous months (March to May). Toxin events above the permitted level reached a peak in August (median above ∼ 13.25%), following a decreasing trend with a low median in October (∼0.88%). PSTs showed a peak of median values in June (∼3.52%); the maximum occurrence decreased in July and presented exceptional cases in August. The occurrence of DA above the safety threshold were mainly zero, with the exception of events in May, July, August, and September.

Relationship Between Monthly Blooms

The apparent seasonal pattern of blooms evident from Figure 4 was further investigated using autocorrelation analysis in the time series. All of the phytoplankton taxa presented a significant seasonal and cyclical monthly correlation at all monthly lags (Supplementary Figure 2). The positive autocorrelation at shorter lags (1–3) suggests the frequency of HABs is linked to their occurrence in the previous month. The negative correlation at lags from 4 to 8 months indicate the different seasonal response. This confirms that seasonality is an important driver for the occurrence of HAB events. While this general trend was confirmed, PACF analysis indicated the seasonal pattern was specific for each HAB taxa. Dinoflagellates Dinophysis spp. and Alexandrium spp. showed five lags partial negatively autocorrelated (up to 6 months); whilst two and one positive lag, respectively. Pseudo-nitzschia spp. showed a low number of autocorrelated lags (three negative and three positive lags), suggesting seasonality has less of an impact in influencing the occurrence of events than for the other two genera.

The biotoxins autocorrelation pattern was group specific (Supplementary Figure 3). DSTs and PSTs showed a cyclical pattern by ACF, the former displaying both positive and negative correlations at all lags values while the latter presented only positive autocorrelated lags at 11, 12, and 13 months. DA only showed two positive lags at 1 and 13 months. DSTs showed a higher number of lags autocorrelation by PACF, at low lags and high lags months. PSTs showed one negative partial autocorrelated lag and two positive (10, 11 months). DA only presented one lag autocorrelated (13 months) by PACF.

Temporal-Spatial Analysis

Temporal and Spatial Dynamics

The PCA analysis identified common patterns that emerged from the similarities (correlation) of temporal dynamics and their association with particular aggregation boxes across the coastline. These are portrayed in a low number of dimensions, principal components (PC), that identify the main variables that play a role in regulating HAB variance. The first principal component (PC1) described ∼56% of Dinophysis spp. monthly variance, the highest value amongst the phytoplankton groups. PC1 also described ∼42% of the variance for Alexandrium spp. and the lowest ∼27% for Pseudo-nitzschia spp. Dinophysis, Alexandrium, and Pseudo-nitzschia showed high and positive PC1 sample scores in years and months of high occurrence of HAB events (Figures 5A–F). Thus, PC1 could be related to the common coastal temporal patterns and seasonality marking HABs occurrence. Spatial patterns were assessed using monthly factor loadings (correlations between frequencies and PC1 scores) for the time series in each aggregation box. Dinophysis spp. had positive PC1 loadings for all aggregation boxes and a noticeable difference between areas in the southwest coast, the northwest coast, east coast, Isle of Lewis (boxes 13,14), and the Northern Isles (Orkney and Shetland Islands) (Figure 5G). This suggests blooms follow a regional pattern, with an important difference in seasonality between southwest and northwest areas. For Alexandrium, locations with similar loadings were dispersed throughout the coastline (Figure 5H). The highest Alexandrium PC1 loadings (0.21–0.3) were in two southwest areas, the central west coast, the northwest coast and the Northern Isles. Lower Alexandrium PC1 loadings were found for boxes in sheltered inner lochs of the west coast, Isle of Lewis, east coast and southern areas of the Shetland Islands. Areas showing similar correlations with PC1 for Pseudo-nitzschia spp. were widely dispersed with high positive values across the southwest, central coast, northwest, Isle of Lewis and the Shetland Islands (Figure 5I). In summary, PC1 for each species captured the temporal pattern of blooms and hence the role of seasonality. The lower proportion of variation explained by PC1 for Pseudo-nitzschia spp. suggests local conditions or other variables might have a greater influence on the incidence of HABs for this genus than for Dinophysis spp. and Alexandrium spp.

FIGURE 5

Figure 5. First principal component (PC1) of Dinophysis spp., Alexandrium spp. and Pseudo-nitzschia spp. (A–C) Monthly sample scores, (D–F) annual sample scores, (G–I) monthly factor loadings per aggregation box.

The second principal component (PC2) explained a lower proportion of the variance of HAB incidence than PC1: ∼6% for Dinophysis, ∼7% for Alexandrium and ∼8% for Pseudo-nitzschia but was instructive in demonstrating the timing and location of blooms around the country. The PC2 sample scores portray positive or negative scores according to the timing in the year (Figures 6A–C). Annual sample scores display an up-down trend of between 1 and 3 years, with no particular association with high incidence years, which could indicate an interannual shifting of bloom location (Figures 6D–F). Dinophysis showed high and positive monthly PC2 loadings (0.2–0.4) in the southwest and one area in the northeast coast (Figure 6G). Negative Dinophysis PC2 loadings were rare, the lowest (−0.4–−0.2) observed in the Isle of Harris and Northern Isles (Shetland and Orkney Islands). The pattern suggests there is a temporal shift of bloom location with the high incidence of blooms that peak in July in the Isle of Lewis and the Northern Isles, with “early” (May) and “late” (September) blooms occurring on the southwest coast (Firth of Clyde). Alexandrium PC2 loadings were primarily positive in west coast areas, whilst negative values were focused on the Isle of Lewis (aggregation box 13, 14), the Northern Isles (aggregation box 1–5), and two areas in the central west coast (aggregation box 18–20) (Figure 6H). This suggests “early” blooms are focused on the south, central and north coast and “summer” blooms are associated with the Isle of Lewis and the Northern Isles (Figure 6H). For Pseudo-nitzschia, “early” and “late” season blooms (March, April, and September, October, respectively) were associated with negative PC2 sample scores; whilst “summer” blooms during May to August mainly showed positive scores (Figure 6I). Pseudo-nitzschia PC2 loadings were spatially variable with little regional pattern. For example, negative PC2 were evident in the west coast and near the island of Mull, but the adjacent Loch Fyne (aggregation boxes 21, 22) exhibited high values. “Early” (March-April) and “late” (September-October) blooms were associated with the central west coast (aggregation box 19, 20), east coast (aggregation boxes 6–7) and the south of the Shetlands (aggregation box 4), while “summer” blooms were related to the Isle of Lewis, Orkney Island, Shetland Islands (except the previous mentioned aggregation box) and the southwest coast.

FIGURE 6

Figure 6. Second principal component (PC2) of Dinophysis spp., Alexandrium spp., and Pseudo-nitzschia spp. (A–C) Monthly sample scores. (D–F) Annual sample scores. (G–I) Monthly factor loadings per aggregation box.

Harmful Algal Bloom Occurrence Risk

The PCA allowed the identification of areas that shared temporal dynamics, and the role of seasonality in regulating the different HAB genera. However, PCA results are independent of the magnitude of blooms, and hence to evaluate spatial trends in abundance we used the K-means test to cluster aggregation boxes that showed similar bloom size and timing, hence identifying similar levels of risk on a monthly and annual scale. In all cases, the elbow plot suggested using four clusters for analyses of monthly data. The K-means clustering explained almost ∼70% of variance in monthly Dinophysis and Alexandrium HAB frequencies (between_SS/total_SS = 71.9%, between_SS/total_SS = 68.2%, respectively) and Pseudo-nitzschia ∼50% (between_SS/total_SS = 52.1%). The monthly clusters of Dinophysis followed a similar trend, peaking in July (Figure 7A). However, bloom occurrence across the year varies between clusters with the “very-high” cluster exhibiting elevated risk both “early” (in April) and “late” (in October) in the year. In contrast to Dinophysis where clusters rarely overlapped, Alexandrium (Figure 7B) and Pseudo-nitzschia (Figure 7C) exhibited similar frequency of “early” and “late” blooms for all clusters. The “high” cluster exhibited “early” blooms by the increase of K-means from March to a peak in May for Alexandrium and April for Pseudo-nitzschia. The “very-high” cluster rapidly increased later in the year than the “high cluster” but then markedly exceeded it and all other clusters from June to September. Alexandrium “moderate” and “low” clusters showed similar seasonal patterns, but at different magnitudes. An increase in Pseudo-nitzschia risk was noted for the “very-high” and “high” clusters in spring (March-April), with the “moderate” cluster also being somewhat elevated in April. All clusters decreased in May, before a further summer increase during which bloom patterns differed amongst all clusters from June to October. For all genera, there was a distinct “low” cluster that exhibited consistently lower K-means peaking in July for Dinophysis, May and June for Alexandrium and June for Pseudo-nitzschia, remaining somewhat elevated until September for the latter genus. Thus Dinophysis blooms followed a similar temporal occurrence across clusters while Alexandrium and Pseudo-nitzschia showed a shift between “early” and “summer” blooms from low to high risk clusters.

FIGURE 7

Figure 7. Dinophysis spp., Alexandrium spp., and Pseudo-nitzschia spp. (A–C) Monthly K-means for each cluster. (D–F) Monthly K-means according to the HAB risk per aggregation box distributed in Scotland.

Members of “very-high” and “high” Dinophysis HAB-risk clusters were in the southwest, central west, northwest coast and south Shetland Islands (aggregation boxes: 3, 4, 9, 10, 12, 13, 15–17, 19–23, 25) (Figure 7D). The “low” risk cluster is scattered around the country, primarily on the Isle of Harris, Orkney Island, one area in the Shetland Islands, but also including the only east coast area and one mainland west coast location (aggregation boxes 2, 5, 6, 14, 18). For Alexandrium, the “very-high” HAB-risk cluster was only in one area in the west of Lewis and another in the south of the Shetland Islands; whilst the “high” cluster areas were more widespread, located in the southwest coast, central west coast, northwest coast, northeast and north of the Shetlands (Figure 7E). Areas assigned to the “low” Alexandrium HAB-risk cluster were in the southwest, southeast, Orkney and one in the Shetlands. The “very-high” risk cluster for Pseudo-nitzschia is located uniquely in two areas in the south of the Shetland Islands (Figure 7F). A high proportion of the areas correspond to the “low” cluster, the majority being in the southwest coast.

Clustering the annual data into four clusters explained a similar proportion of the variation to the clustering of monthly data: ∼60% of the variance in Dinophysis spp. (between_SS/total_SS = 63.5%), ∼50% for Alexandrium (between_SS/total_SS = 50.8%) and ∼60% for Pseudo-nitzschia (between_SS/total_SS = 61.1%). Risk varied from year to year, Dinophysis clusters displayed similar annual K-means patterns across the years (Figure 8A), but the relative level or risk remains consistent (i.e., “Very-high” risk areas had consistently higher risk than “high” risk areas). A noticeable difference was observed for the annual trends of Alexandrium (Figure 8B) and Pseudo-nitzschia (Figure 8C), with K-means showing considerable interannual variability for all clusters. Despite the “very-high” cluster scoring the highest risk in the majority of years, the “high” cluster exhibited similar risk levels from 2011 onward for Alexandrium, and exceeded “very high” cluster risk levels in 2015 for Pseudo-nitzschia. The distribution of clusters, in general, follows the monthly Dinophysis (Figure 8D) and Pseudo-nitzschia (Figure 8F) K-means analysis. For Dinophysis, the “very-high” and “high” clusters were mostly on the southwest coast, the islands of the central and northern west coast and the Shetland Islands, with “low” risk clusters on the east coast and Orkney Islands. For Pseudo-nitzschia, the “very-high” cluster is only in the south of Shetland Islands, with the “low risk cluster mainly on the southwest coast. The location of the Alexandrium clusters differed between monthly and annual values, “very-high and “high” risk is scattered in the west coast, west and Northern Isles (Figure 8E).

FIGURE 8

Figure 8. Dinophysis spp., Alexandrium spp., and Pseudo-nitzschia spp. (A–C) Annual K-means for each cluster. (D–F) Annual K-means according to the HAB risk per aggregation box distributed in Scotland.

The spatial risk of the aggregation boxes used the monthly K-means classification, this was averaged and classified in a cluster rank (Supplementary Table 1). The total average risk was 2.36, with a maximum value of 3.67 in the south of the Shetland Islands and a minimum of 1.00 in the Orkney Islands and the southeast coast. The highest combined risk locations (average above 3) were scored by the south Shetlands Island (box 4) and the South Isle of Skye (box 16). High scores (average of 3) included seven locations: the southwest of the Shetland Islands (box 3), the northwest coast (box 9), the north of the Isle of Lewis (box 13), southeast of Skye (box 15), Loch Ailort (box 17), the Isle of Mull (box 19) and the Isle of Arran (box 25). Five locations were identified as the lowest combined risk (average below 1.50), located in the south of the Kintyre Peninsula (box 24), central west inner loch-Loch Linnhe (box 18), the south east coast (box 6), the Orkney Islands (box 5) and the northwest Shetland Islands (box 2). High and low risk locations were significantly associated for Dinophysis spp. and Alexandrium spp. (correlation among risk scores, r = 0.44, p < 0.05), and Alexandrium spp. with Pseudo-nitzschia spp. (0.48, p < 0.05), but risk scores for Dinophysis spp. showed only a weak non-significant correlation with those for Pseudo-nitzschia spp. (0.29, p > 0.05).

Harmful Algal Bloom Relationships With Environmental Variables

Correlation was applied to evaluate any possible association between HABs annual frequency and annual mean of environmental variables (Table 2). This was conducted at a countrywide level as environmental information is not collected with sufficient spatial resolution for a regional analysis. The only positive significant correlation was found between the NAO index and the genus Alexandrium and corresponding PSTs toxins. Dinophysis spp. and SST were correlated negatively significant on an annual basis. There were no significant relationships between other taxa or biotoxins with SST or wind (intensity and direction).

TABLE 2

Table 2. Analyses of relationships of annual incidence of HABs and toxic events above the safety threshold with annual averages of environmental variables from 2006 to 2018.

Discussion

Interannual Variability of Harmful Algal Blooms

The analysis of temporal-spatial patterns is crucial for determining the risk of phytoplankton species of major concern to shellfish aquaculture in Scotland. Over the 15 years of the time series, we observed high interannual variability between years with no statistically significant trend in harmful bloom or shellfish biotoxin frequency for all taxa and toxins. This is consistent with other site specific case-studies (Moore et al., 2009; Díaz et al., 2016) and agrees with the recent analysis of global data series (Hallegraeff et al., 2021), which attributed the perceived increase in HAB events to intensified monitoring.

Furthermore, the frequencies of annual blooms and toxic events were found not to be correlated across years, suggesting annual HAB occurrences are independent and are not influenced by the previous year. The only exception was observed at a lag of 1 year for Pseudo-nitzschia spp. Given that this taxon is not thought to form resting stages (Azanza et al., 2018) this observation could potentially be related to a common biological response in consecutive years to similar environmental conditions. An earlier study (Rowland-Pilgrim et al., 2019) did not find a correlation between Pseudo-nitzschia spp. cell density and DA concentration in the same time series from 2008 to 2018.

The annual frequency of blooms of the dinoflagellates Dinophysis spp. and Alexandrium spp. were similar and significantly correlated. Hence, although the temporal and spatial dynamics of blooms of the two taxa within and between years differ, the data suggest that on a countrywide basis “good” and “bad” years exist for harmful dinoflagellate blooms. However, no correlation with readily available environmental drivers (SST, wind) were found. Dinophysis is a holoplanktonic genus with blooms having been described to develop offshore with advective transport to the coast (Whyte et al., 2014). In contrast, Alexandrium spp. forms cysts that are deposited in the benthos and hatch under the right conditions. Hence, the observed dynamics could be related to the group trait or succession pattern described by Smayda and Reynolds (2001), overcoming life cycle differences. A deeper understanding of the role of environmental variables is needed to determine the triggers for blooms of these harmful taxa.

Temporal and Spatial Bloom Dynamics

HAB occurrence in general followed the well-known seasonal pattern, with most blooms occurring from March to October. Dinophysis spp. and Pseudo-nitzschia spp. blooms reached a peak in July, and Alexandrium spp. exhibited a plateau from May to July. Pseudo-nitzschia, similarly to Alexandrium, exhibited two peaks throughout the year in spring and summer associated with specific locations. This confirms, countrywide, the characteristic temporal profile for this taxon demonstrated from local monitoring by Fehling et al. (2006). The crucial role of seasonality in regulating bloom dynamics was confirmed by the temporal patterns and significant autocorrelation between monthly occurrences for all genera. Seasonality is considered as a key driver of the phytoplankton composition, including dinoflagellate and diatom species in the west and east coast of Scotland (Bresnan et al., 2015). The environmental conditions during winter hinder the development of blooms, since Scottish waters are exposed to frequent storms leading to intense mixing of the water column; with light intensity and water temperature being below the optimal conditions of growth. During spring, increasing irradiance intensity and duration, and an increase of water temperature and general decrease of storms result in ideal conditions for phytoplankton cell growth and higher incidence of blooms. Initially these are primarily diatom dominated, and then as nutrient concentrations decrease and stratification intensifies dinoflagellates become more prevalent (Davidson et al., 2011). These patterns are associated with Scottish transitional waters being regulated by the seasonal heating and cooling cycle (Edwards et al., 2013). However, despite the critical role of seasonality, its influence varied for the different HAB taxa. Greatest control was observed for Dinophysis spp. with seasonality, reflected in PC1, explaining over 50% of variance, perhaps because of its heterotrophic nature and the requirement for ciliate prey that are most prevalent in summer months (Reguera et al., 2012). In contrast, the lower proportion of variance explained and number of monthly lags autocorrelated for Pseudo-nitzschia spp. suggest seasonality is not as crucial factor for the development of blooms of this genus, with other factors such as local environmental conditions making dissolved silica available or stochastic effects being more important. As noted above, the differences in the seasonal peaks of the seriata and delicatissima groups of Pseudo-nitzschia are also responsible for the relative lack of seasonality of the genus as a whole.

Dinophysis presented the greatest bloom frequency. The countrywide interannual patterns observed are consistent with abundance observations at two targeted locations on the west coast of Scotland (Swan et al., 2018) and offshore patterns found in the Northeast Atlantic from Continuous Plankton Recorder data (Dees et al., 2017), and do not confirm the model based predictions of Gobler et al. (2017) of a climate driven increase in the abundance of Dinophysis acuminata in the region. The monthly and annual Dinophysis analysis using K-means categories, determined that the blooms respond similarly between clusters at a regional scale, but with varying levels of risk. The “very-high” cluster was predominantly located in the southwest coast, particularly in the Firth of Clyde, and was associated with “early” and “late” season blooms, as well as the highest mid-summer bloom frequency. The spatial distribution of the two highest K-means clusters was located predominantly on the west coast and the south of the Shetland Islands for both monthly and annual data series. Whereas the central west, northwest, Isle of Lewis and Shetland Islands are likely to respond to seasonal conditions in a shorter season, according to the principal components. These areas could be susceptible to the advective transportation of cells, since this has been attributed as the main reason of Dinophysis blooms rather than the in situ growth of cells (Smayda and Reynolds, 2001). Conditions such as local seasonal stratification or the development of frontal regions have previously been shown to regulate Dinophysis blooms at the species level with late summer Dinophysis acuta blooms occurring only in some years in two hydrologically different locations (Firth of Clyde and Loch Ewe) (Swan et al., 2018).

For Alexandrium, the spatial statistical tests indicated a shift of bloom location between “early” and “late” (summer) occurrence. “Early” blooms corresponding to the “high” risk cluster are geographically dispersed, predominantly on the west coast (south, central, north), while summer blooms portrayed by the “very-high” risk cluster are predominantly in the northeast of the Isle of Lewis and south of the Shetland Islands. The use of different cluster numbers (k = 5, 4, 3) consistently identified these areas as high-risk areas. These results confirm earlier observational studies that identified Orkney and Shetland Islands as Alexandrium “hotspots” in June and July (Bresnan et al., 2006). A decade-long time series (1996–2005) showed high densities in the Western Isles, Orkney and Shetland Islands; the location reaching the highest average annual cell densities and PSP levels amongst all areas (Brown et al., 2010). The “low” risk clusters located in two southwest locations, the east coast and the Orkney Islands differed from another study that found low maximum cell densities in lochs in the central west coast (Brown et al., 2010). Furthermore, the lowest combined risk for all taxa, including Alexandrium, is located in the south of the Kintyre Peninsula; however, the nearby area of Campbeltown recorded the highest PSTs value within the 2007–2020 data series (27,822 μg/kg) in April 2015.

Pseudo-nitzschia bloom occurrence showed, as did Alexandrium, a temporal-spatial shift between “early” and “summer” blooms. The monthly occurrence of all clusters increased during spring, commonly associated with the increase of nutrient availability, light intensity and duration. “Very-high” and “high” risk clusters exhibited similar K-means levels early in the year. A decrease of bloom events was observed after the spring bloom (March-April), showing low levels in May. This decreasing trend has been previously associated with the depletion of nutrients in the water column (Fehling et al., 2006). The spatial statistics identified the shift of bloom dynamics, with a higher occurrence during summer, especially for the “very-high” cluster, focusing on the south of the Shetlands (aggregation box 4). The increase in bloom occurrence during summer in our results is consistent with another study using the FSS data series, that identified an increasing trend of the taxon’s cell abundance during summer and rare events of domoic acid above the permitted levels in shellfish tissue (Rowland-Pilgrim et al., 2019).

Spatially, high and low-risk areas are mainly differentiated by their orientation rather than geographical proximity, but also varying according to the taxa. The “high-risk” locations mainly face in a southwest direction. On the Scottish west shelf, water flows in a predominantly northwards direction governed by the Scottish Coastal Current (Hill and Simpson, 1988). The advection of cells from offshore waters, driven by this current has previously been identified as a means of transportation of Alexandrium tamarense (Hill et al., 2008) and Karenia mikimotoi (Hill et al., 2008; Gillibrand et al., 2016). The advection from offshore waters has also been associated with the transport of Pseudo-nitzschia spp. (Fehling et al., 2012) and Dinophysis spp. cells (Paterson et al., 2017). This potentially explains the “very-high” risk of Dinophysis in southwest facing sea lochs in the South of the country.

Phytoplankton and Biotoxin Relationship

Contamination of shellfish with DST is related to the presence of Dinophysis cells in the water column, with the shellfish toxin concentration gradually increasing until it potentially surpasses the MPL. The high incidence and seasonal trend observed for Dinophysis and DSTs occurrence above the safety threshold and the significant correlation of such annual frequencies demonstrated the close link between this phytoplankter and shellfish toxicity. Dinophysis cells exhibit diverse levels of toxicity through the presence of different species with different toxin profiles (Swan et al., 2018). Dinophysis acuminata (linked with the production of okadaic acid and dinophysistoxin-1) has been found to be the predominant species in late spring and summer (June-August) in Scottish waters (Stern et al., 2014; Swan et al., 2018), whilst D. acuta (producer of dinophysistoxin-2) is found occasionally from July onward, but can be associated with extended toxic events (Swan et al., 2018). Our study found a prolonged toxicity of DSTs in shellfish events in late autumn (September-October) when the occurrence of Dinophysis above the safety threshold is low, likely demonstrating the slow depuration of these toxins from shellfish flesh (Svensson, 2003). Similar temporal dynamics of D. acuminata and D. acuta have also been found in Iberia (Moita et al., 2016).

Previous studies have demonstrated that a number of toxic and non-toxic species of Alexandrium exist in Scottish waters. PST shellfish toxicity in Scotland is thought to be primarily related to toxic Alexandrium catenella (formerly identified as Gonyaulax excavata, Alexandrium fundyense, and/or Group I ribotype/North American Alexandrium tamarense). However, monitoring identifies cells to genus only and the morphologically indistinguishable non-toxic A. tamarense can co-occur (Touzet et al., 2010) and hence be enumerated in the monitoring program. The presence of non-toxic strains of A. minutum in Scottish waters has been also detected (Lewis et al., 2018). In an ideal scenario, regulatory monitoring would enumerate only toxic Alexandrium. However, difficulty in discriminating between species on a rapid routine basis within a monitoring program has led the regulator to take a precautionary approach for the monitoring of Alexandrium (and our other genera of interest). Observations of high concentration of cells but low or absent toxicity exist (Higman et al., 2001). Moreover, temperature (and other environmental drivers) may impact excystment (Collins et al., 2009), competition (Davidson et al., 1999; Eckford-Soper et al., 2016) and toxicity (Flynn et al., 1994). However, our study found a correlation between the annual frequency of Alexandrium cells above the safety threshold and toxicity levels above the MPL. This concurs with other studies of cell abundance and PST concentrations (Bresnan et al., 2006; Stobo et al., 2008; Brown et al., 2010) and hence confirms the suitability of identifying Alexandrium to genus level as a means of ensuring shellfish safety.

When viewed as a time series, the annual frequency of Pseudo-nitzschia events above the safety threshold and domoic acid did not correlate, consistent with studies using cell densities and biotoxin levels (Hinder et al., 2011; Bresnan et al., 2017; Rowland-Pilgrim et al., 2019). Conversely, another study limited to two locations in Scotland for a period of 2 years (2001–2003) found a significant correlation using a 1.5 week temporal lag, associated with the fast uptake and depuration of the biotoxin by M. edulis (Bresnan et al., 2017). The non-correlation between Pseudo-nitzschia cells densities and domoic acid has been attributed to the variation in toxin content according to the group, species within the group, and size of the cells (Smayda, 2006). The rare occurrence of domoic acid events above the MPL during spring Pseudo-nitzschia spp. blooms could be explained by the presence of P. delicatissima group species containing low levels of DA per cell, or non-toxic members of the P. delicatissima group, predominant during spring (Fehling et al., 2006). Toxic events in late summer could be explain by the presence of the P. seriata group which contains species known to synthesize high levels of domoic acid (Smayda, 2006) and which are mostly found during the summer months (Fehling et al., 2006) and are infrequently observed during spring in Scottish waters (Paterson et al., 2017). A further confounding factor is that the decrease of the frequency of biotoxins above the maximum permitted level, specifically for PSTs, might be due to the high number of site closures for DSTs predominantly in July and August. Sampling is suspended, even if Alexandrium was still present in the water samples, hence the presence of a single biotoxin above the regulatory limit will cause under-sampling for the others.

Harmful Algal Blooms and Environmental Relationship

A range of conceptual models have been published to explain the relationship between different HAB genera and environmental conditions (see reviews of Reguera et al., 2012; Anderson, 2015; Bates et al., 2018). However, such models typically relate blooms to local environmental conditions that are not surveyed with the same fine temporal and spatial resolution as phytoplankton and biotoxins. Hence, to further evaluate the relationship between HAB events and those “large scale” environmental variables that are available in accessible databases, we used the correlation function to identify the direction and strength between these at a countrywide scale. At a climatological level, we found a significant positive correlation between the NAO index and the frequency of PSTs and producing genus Alexandrium. Positive NAO index leads to windier and turbulent conditions in the column water. Reports have suggested dinoflagellate blooms are enhanced under calm conditions (Smayda and Reynolds, 2001; Hinder et al., 2012). In contrast, the positive correlation between Alexandrium and PSTs could suggest turbulent conditions are strongly related to coastal bloom development. This again suggests HAB occurrences are not limited to regular seasonal environmental conditions; but that local conditions and atmospheric anomalies are likely to play a more important role in regulating these. The differences in correlation between the NAO index and the other environmental variables suggest the effect of the climatological system is not straightforward in that the NAO index is related to more climatological effects than just wind variation, and other variables are likely to be important in HAB development. The significant negative correlation between SST and Dinophysis blooms is consistent with another study that assessed the relationship between several dinoflagellate genera and this environmental variable (Hinder et al., 2012). The weak relationship between SST, wind and HAB events for most taxa and biotoxins might also be explained by the large geographical scale used and high variability of the values in the time series. While genus based monitoring is thought to be suitable to safeguard human health within regulatory shellfish safety monitoring, the results strongly suggest cell enumeration to the group or species level in combination with higher resolution environmental monitoring would enhance the potential for phytoplankton counts to provide the early warning of shellfish toxicity that the shellfish industry require to minimize the economic losses associated with HABs.

Conclusion

The frequency of phytoplankton and biotoxin events above the safety threshold sampled by the Scottish shellfish safety monitoring programme did not show a significant increasing or decreasing trend in the last 15 years. HABs exhibited considerable variability, with non-correlation between annual frequencies of event indicating that blooms are not interannually related. A traditional seasonal pattern is observed for all taxa, reaching a monthly peak in July for Dinophysis spp. and Pseudo-nitzschia spp., with a plateau of blooms from May to July for Alexandrium spp. DSTs surpassed the MPL every year in the time series, PSTs in fewer years. The annual frequency of bloom and toxic events were significantly correlated with Dinophysis and Alexandrium and their respecting toxins. DA in contrast showed rare toxic events throughout the time series and correlation with Pseudo-nitzschia annual frequencies were non-significant. This has been attributed to the presence of non-toxic organisms, different factors regulating cell growth and biotoxin production and mussels’ biotoxin assimilation rates. Dinophysis spp. and DSTs showed the highest frequency of events surpassing safety thresholds, representing the highest risk to the shellfish aquaculture in Scotland. The temporal analysis and spatial statistics (PCA and K-means) evidenced the crucial role seasonality plays on regulating timing, duration, intensity and location of HABs. But even with this consistent trend, known to be regulated by the cyclical Scottish waters, HAB dynamics vary according to the taxa and across the country. Relatively few countrywide relationships between HABs and environmental factors were found. Only the NAO index was significantly related to the HAB events of Pseudo-nitzschia spp. and PSTs. Despite this variable’s influence on environmental conditions in the water column, the non-correlation between blooms and SST and wind (direction, intensity) demonstrates the relationship between the NAO index and blooms is not straightforward and HAB events are regionally controlled.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Phytoplankton and biotoxin data can be accessed through the Food Standards Scotland (FSS).

Author Contributions

FG, KD, and MB contributed to the design, statistical tests used, analysis and writing. SS contributed to the data analysis, interpretation and writing. AT contributed to the writing, especially the toxin analysis. All authors reviewed the results and approved the final version of the manuscript.

Funding

FG was funded by an Ocean Risk Ph.D. Scholarship from AXA XL. Additional support (KD) was received from the RCUK projects CAMPUS (NE/R00675X/1) and Off-Aqua (BB/s004246/1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

All phytoplankton and biotoxin data were accessed through the Food Standards Scotland funded regulatory marine biotoxin monitoring program.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.785174/full#supplementary-material

Footnotes

References

Anderson, D. (2015). HABs in a changing world: a perspective on harmful algal blooms, their impacts, and research and management in a dynamic era of climatic and environmental change. Harmful Algae 2012, 3–17. doi: 10.1038/nature13478

PubMed Abstract | CrossRef Full Text | Google Scholar

Anderson, D. M., Boerlage, S. F. E., and Dixon, M. B. (2017). in Harmful Algal Blooms (HABs) and Desalination: A Guide to Impacts, Monitoring and Management, Vol. 78, eds D. M. Anderson, S. F. E. Boerlage, and M. B. Dixon (Paris: Intergovernmental Oceanographic Commission of UNESCO), 539.

Google Scholar

Anderson, D. M., Stock, C. A., Keafer, B. A., Bronzino Nelson, A., McGillicuddy, D. J., Keller, M., et al. (2005). Alexandrium fundyense cyst dynamics in the Gulf of Maine. Deep. Res. Part II 52, 2522–2542. doi: 10.1016/j.dsr2.2005.06.014