Harmful Algae and Oceanographic Conditions in the Strait of Georgia, Canada Based on Citizen Science Monitoring

Esenkulova, Svetlana; Suchy, Karyn D.; Pawlowicz, Rich; Costa, Maycira; Pearsall, Isobel A.

doi:10.3389/fmars.2021.725092

ORIGINAL RESEARCH article

Front. Mar. Sci., 09 September 2021

Sec. Marine Ecosystem Ecology

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

This article is part of the Research Topic Harmful Algal Blooms in Marine Environments: From Biologically Active Compounds to Species Diversity View all 6 articles

Harmful Algae and Oceanographic Conditions in the Strait of Georgia, Canada Based on Citizen Science Monitoring

$\r\nSvetlana Esenkulova*&#x;$ Svetlana Esenkulova^1*†

Karyn D. Suchy^1,2†

Rich Pawlowicz³

Maycira Costa²

Isobel A. Pearsall¹

¹Pacific Salmon Foundation, Vancouver, BC, Canada
²Department of Geography, University of Victoria, Victoria, BC, Canada
³Department of Earth, Ocean, and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada

In British Columbia (BC), harmful algal blooms (HABs) regularly cause severe economic losses through finfish mortalities and shellfish harvest closures due to toxin accumulation, gill damage, or hypoxia. As there is no routine governmental monitoring of HAB phenomena in BC, HAB variability, and its potential links to environmental drivers are not well understood. Here we present results from a well-managed citizen science program which collected an unprecedented 4 year, high-resolution (∼bi-monthly, ∼80 stations) dataset of harmful algae (HA) concentrations and corresponding physical and chemical properties of seawater throughout the Strait of Georgia (SoG), BC. Analysis of this dataset revealed statistically significant interannual and seasonal relationships between environmental drivers and the most common HA taxa: Rhizosolenia setigera, Dictyocha spp., Alexandrium spp., Heterosigma akashiwo, Chaetoceros convolutus, and C. concavicornis. HABs exhibited significant interannual variations; specifically, no HABs were found during the summer of 2015, blooms of Dictyocha occurred in 2016 and 2017, and dense blooms of Heterosigma and Noctiluca occurred in 2018. In addition, HA prevalence corresponded with negative effects observed in local aquaculture facilities where higher toxins concentrations (causing Paralytic and Diarrhetic Shellfish Poisonings) in shellfish flesh were detected during years with greater abundance of Alexandrium and Dinophysis. Furthermore, salmon mass mortality at fish farms corresponded to years with high concentrations of Heterosigma and Dictyocha. As such, these results highlight the need for long-term data to evaluate the potential role of HA as a stressor on the SoG ecosystem.

Introduction

The public health and economic impacts of harmful algal blooms (HABs) in the coastal ocean have increased in frequency, intensity, and distribution in recent years (Hallegraeff, 2004; Heisler et al., 2008), and it is now accepted that “global expansion” of HAB phenomena is real (e.g., Smayda, 1990; Anderson et al., 2012). In British Columbia (BC), HABs regularly cause severe economic losses through finfish/shellfish mortalities and shellfish harvest closures due to toxin accumulation, gill damage, or hypoxia (Horner et al., 1997; Whyte et al., 1997; Taylor and Harrison, 2002; Rensel and Whyte, 2004; Haigh and Esenkulova, 2011, 2014; Bates et al., 2020).

Multiple species are responsible for HABs in BC coastal waters, with the majority outlined below. Heterosigma akashiwo, the most prominent fish-killing algae in the world (Hallegraeff et al., 2004), has had devastating effects on BC aquacultured salmon since the 1960s (Taylor and Haigh, 1993). Currently, it is believed to be one of the primary causes of BC fish mortalities with estimated annual losses of millions of dollars (Taylor and Horner, 1994; Rensel and Whyte, 2004; Haigh and Esenkulova, 2014). Further, H. akashiwo blooms may adversely affect herring and wild salmon survival (Rensel, 2007; Rensel et al., 2010). Chaetoceros convolutus and C. concavicornis are known to kill fish as a result of gill damage (Albright et al., 1993; Taylor and Harrison, 2002; Hallegraeff et al., 2004) and have caused fish losses in coastal BC (Albright et al., 1993; Haigh and Esenkulova, 2014). Rhizosolenia setigera is a non-toxic diatom with long, needle-like spikes, and generally, is not considered harmful (Moestrup et al., 2008). However, in BC, R. setigera is known to cause fish mortalities due to gill clogging (Haigh and Esenkulova, 2011). Blooms of Dictyocha have previously been associated with fish kills in Europe (Henriksen et al., 1993) and BC (Rensel and Whyte, 2004; Haigh and Esenkulova, 2014; Haigh and Johnson, 2019). Several Alexandrium species produce saxitoxin causing Paralytic Shellfish Poisoning (PSP) (Hallegraeff et al., 2004). In BC, it results in numerous shellfish harvesting closures every year (Taylor and Horner, 1994) and is implicated in fish kills at fish farms (Haigh and Johnson, 2019). Some species of Dinophysis produce toxins that cause Diarrhetic Shellfish Poisoning (DSP) (Yasumoto et al., 1980) and have been reported to be common in BC since at least the 1990s (Taylor and Horner, 1994), with the first DSP outbreak confirmed by Taylor et al. (2013). Other species known to cause HABs in BC coastal waters include Chrysochromulina spp., Margalefidinium fulvescens, Noctiluca scintillans, Prymnesium spp., Pseudochattonella verruculosa, Pseudo-nitzschia spp. (Whyte et al., 1997, 2001; Haigh and Esenkulova, 2011; Haigh et al., 2014).

Managing the risk from HABs is therefore necessary. However, the most important problem facing managers of such risk in BC is the lack of information on the abundance, distribution, and population dynamics of harmful algae (HA) (Taylor and Horner, 1994). As there is no routine governmental monitoring and research of HABs in BC, HAB variability and its potential links to environmental drivers are not well understood. Knowledge about HAB ecology in BC is critically lacking. While a number of studies have addressed the distribution and ecology of some HA species (e.g., Haigh and Taylor, 1990; Taylor and Harrison, 2002; Pospelova et al., 2010; Rensel et al., 2010; Esenkulova et al., 2020), these studies have been limited to specific locations and were based on three or fewer years of data. Thus, there remains a pressing need for across-Strait analysis of multi-year time series to better understand HABs in BC (Haigh and Taylor, 1991; Taylor and Harrison, 2002).

Recognizing HAB problems and subsequently establishing routine monitoring programs can help minimize the effects of HABs, thus protecting public health and marine resources (Hallegraeff, 2004; Anderson et al., 2012; Wells et al., 2015). A citizen science program was started in 2015 as part of the Salish Sea Marine Survival Project (SSMSP), an initiative of the Pacific Salmon Foundation (Canada) and Long Live the Kings (United States). While one of the main goals of the SSMSP was to better understand bottom-up factors impacting salmon survival in the region, an additional concern of interest to the project is the possible direct impacts of harmful algae on wild salmon populations in the Salish Sea. The citizen science program was designed to collect data on HA concentrations and corresponding physical (temperature, salinity, stratification, and Secchi depths) and chemical (nitrates, phosphates, and silicates) properties of seawater, as well as monitor the abundance and types of HA, over the Strait of Georgia (SoG) in southwestern BC. Data were obtained from a coordinated sampling program comprising ∼80 stations, each visited approximately 20 times per year. The resulting dataset is the most detailed spatial and temporal in situ record of HA and associated oceanographic data in the SoG to date and is one of the most detailed regional records from the North Pacific’s coastal waters.

This study aims to analyze this dataset to advance our limited understanding of HAB ecology and provide insights useful to fisheries management in the entire Strait of Georgia. Specifically, we: (1) determine spatial and temporal (intra- and inter-annual) patterns in HAB taxa from 2015 to 2018, (2) examine links with environmental parameters for each major HAB taxa, and (3) provide findings that can be used in future retrospective studies evaluating HABs as a regional ecosystem stressor to both natural and aquaculture resources. This study also contributes to establishing a local baseline of HAB dynamics, which, with longer-term datasets, may allow for the development of predictive models.

Materials and Methods

Study Site

The Strait of Georgia is a major part of the semi-enclosed Salish Sea bounded by Vancouver Island and mainland British Columbia in Canada, and Washington State, United States, and is an economically important North American waterway. The SoG is a very biologically productive ecosystem (Parsons et al., 2011), supporting large commercial and recreational fisheries as well as aquaculture (Ketchen et al., 1983). Water circulation in the SoG is driven by estuarine exchange and modulated by tidal and wind mixing (Pawlowicz et al., 2007). In general, the water column in the SoG is highly stratified with the top 50 m characterized by strong seasonal variation with average sea surface temperatures (SST) of 13–16°C in August and 6–7°C in February (Thomson, 1981). In contrast, the layer below 50 m is more uniform in temperature and salinity (∼8 and 30.5°C, respectively). The Fraser River is the most significant source of fresh water in the SoG with an average annual flow of ∼140 km m³ (Mosher and Thomson, 2002), thus supplying about 70% of the dissolved and particulate organic carbon into the system (Johannessen et al., 2003). The Fraser River exhibits peak discharge in summer in response to sunshine-driven snowmelt in the interior mountains.

Algae species are typical for the coastal North Pacific and cold-temperate to temperate waters worldwide (Haigh and Taylor, 1990; Hallegraeff, 2004; Kudela et al., 2005; Taylor et al., 2008; Lewitus et al., 2012). The general characteristics of plankton dynamics have been well described (Stockner et al., 1979; Harrison et al., 1983; Masson and Peña, 2009; Peña et al., 2016). Chlorophyll-a concentrations typically range from <1 mg m³ in the winter to >15 mg m³ during blooms (Suchy et al., 2019); the average annual primary productivity is ∼280 g (Harrison et al., 1983; Sutton et al., 2013).

Field Sampling

Field sampling was carried out by citizen scientists in a program sponsored and managed by the Pacific Salmon Foundation (PSF) in partnership with Fisheries and Oceans Canada (DFO) and Ocean Networks Canada (ONC). Citizen science is a research practice that enlists members of the public to gather data and samples for scientific purposes, a practice which has been remarkably successful in providing a vast quantity of data (Bonney et al., 2009) that would otherwise be unobtainable at similar spatiotemporal scales. Citizen scientists of this project were members of local communities selected because they work on the water in some capacity (e.g., retired fisherman, biologists) and own their own vessels. They were organized into 7–10 patrols (depending on the year), each comprised of about 3 or 4 citizen scientists and a vessel. Fuel expenses were paid, and citizen scientists were provided with a small honorarium.

Citizen scientists were meticulously trained to gather the required data and collect samples/measurements. An initial training was carried out during a 2-day workshop, and each team then performed their first at sea sampling with a staff member with extensive experience in oceanographic sampling for quality control. Each team was visited by the program coordinator about three times a year, newsletters and or updates via group e-mails and social media were issued on a monthly basis, and an annual symposium was held to share citizen scientists experience and provide updates from the technicians and scientists who analyzed measurements and samples.

Sampling was carried out in a coordinated fashion on scheduled dates, approximately 1–3 times per month between February and October from 2015 to 2018. All sampling was conducted during the daylight hours of the maxima and minima periods of the fortnightly tidal cycle, with gaps for weekends and holidays. A total of 92 stations in 13 areas (Figure 1 and Table 1) were sampled between 2015 and 2018; 67 of these stations from nine of the areas were sampled in all 4 years (Supplementary Material 1). The Victoria area was sampled in 2015 only, whereas Galiano and Malaspina areas were only sampled from 2016 to 2018. The Ladysmith area was sampled opportunistically by the Stz’uminus First Nation.

FIGURE 1

Figure 1. Map of the sampling locations. Shading represents southern, central and northern parts of the Strait of Georgia.

TABLE 1

Table 1. Sampling areas, stations, and number of phytoplankton samples.

Sampling sites included both nearshore and open water locations. The proximity of the sampling stations to the shore ranged from 70 to 13,941 m, while station depth ranged between 17 and 418 m. Sampling included vertical profiles of seawater temperature, salinity, chlorophyll fluorescence, and dissolved oxygen content, obtained with standard internally recording conductivity-temperature-depth (CTD) profilers outfitted with additional sensors, as well as surface water samples for phytoplankton analysis, and Secchi disc (Cialdi and Secchi, 1865) measurements. Water samples for nutrient analysis at 0 and 20 m depth and additional samples from depths of 5, 10, and 20 m for phytoplankton analysis were collected at select stations (Figure 1). Extra samples were occasionally collected on non-scheduled dates if visible HAB events were detected. A more complete description of the program is provided in Pawlowicz et al. (2020) and in a Web atlas based on this program (sogdatacentre.ca/atlas), with quality-controlled data freely available from the Strait of Georgia Data Centre (sogdatacentre.ca). As of 2021, the program is ongoing.

Data and Sample Analysis

Conductivity, Temperature, and Depth

Hydrographic profiling for the data analyzed here was carried out using RBR Concerto CTD probes, which were slowly lowered into the water at the standard rate of 0.5 m s^–1 from the surface to 150 or ∼5 m above the seafloor, whichever was shallower. Each CTD cast recorded high resolution (6 Hz) vertical profile data which were subsequently transmitted to an Ocean Networks Canada (ONC) data management system via a tablet application developed by ONC. CTD data were processed using standard procedures (Pawlowicz et al., 2020) and binned at 1 m depth intervals, therefore surface water CTD data hereafter refers to the 1 m depth. Salinities were given on the TEOS-10 Reference Composition Salinity Scale with the salinity anomaly of zero (IOC et al., 2010). Further quality control was provided by the Ocean Dynamics Laboratory, UBC (Pawlowicz et al., 2020).

Based on the 1 m-binned CTD data, a stratification parameter △_t (kg m^–3) was calculated as the difference in density between the surface (average of densities for the top 10 m) and 20 m (Drinkwater and Jones, 1987). A depth of 20 m was chosen for stratification parameter calculations to ensure inclusion of shallow sampling sites (bottom depth <30 m) as previous work carried out in nearby Saanich Inlet, BC showed little discrepancy in results calculated using bottom depths of 20, 50, and 100 m (Suchy et al., 2016).

Nutrients

Water samples for nutrient analysis were collected at 40 of the 92 stations using a Niskin bottle (Figure 1). Approximately 50 mL of seawater was collected at the surface and 20 m depth and (after being filtered through a 0.7 um GFF from 2016 onward) was stored in polypropylene bottles previously rinsed three times with the sample water. Samples were immediately placed in a cooler at 4°C then frozen at −20°C within 3–6 h; no preservatives were added. Nutrient analysis was carried out at Fisheries and Oceans Canada (DFO) and the University of Victoria (UVic) in different years. Nutrient samples collected in 2015 were analyzed using a three channel SEAL Autoanalyzer with an AACE data acquisition program followed by the procedures outlined in Barwell-Clarke and Whitney (1996). Nutrient samples collected in 2016 and 2017 were analyzed on a Lachet Autoanalyzer following QuikChem Method 31-107-04-1-G (nitrate and nitrite; referred to as nitrates thereafter), 31-115- 01- 1H (phosphate), and modified determination of reactive silicate in seawater (Strickland and Parsons, 1972). In 2018, samples were analyzed using an Astoria Autoanalyzer following methods described in Barwell-Clarke and Whitney (1996). Further quality control was provided by the Ocean Dynamics Laboratory, UBC (Pawlowicz et al., 2020).

Phytoplankton

Surface seawater samples for phytoplankton analysis were collected at each station with a bucket. Additional samples were collected at 11 stations from 5, 10, and 20 m using Niskin bottles (Figure 1 and Table 1). The number of phytoplankton samples collected at each sampling station per year, along with station identification and coordinates are provided in Supplementary Material 1. Once collected, samples were immediately preserved in Lugol’s iodine with a final concentration in the sample of ∼1–2% (Throndsen, 1978). All water samples were resuspended, and a 1 mL aliquot was settled on a Sedgewick-Rafter slide (Throndsen, 1995) for at least 5 min prior to identification. The entire slide was examined with a compound light microscope (OMAX Lab Binocular) at 100X magnification.

Phytoplankton analysis followed a method developed by the Harmful Algae Monitoring Program (HAMP, Haigh et al., 2004). Species were identified to the lowest taxonomic level possible based on morphology (Hasle, 1978). The dominant species or group in each sample was enumerated (reported as cells mL^–1), as well as all species known or suspected to have a negative effect on finfish and shellfish in BC (Haigh and Esenkulova, 2011). Since this method does not generally permit reliable identification of nanoplankton, HA species <10 μm (e.g., Chrysochromulina, Prymnesium) were not identified. In addition, a biomass index estimation was employed, and percentages of total biomass in five main groups (diatoms, dinoflagellates, raphidophytes, other phytoplankton, and microzooplankton) were estimated (Haigh et al., 2004).

HA that were identified and enumerated in this study were Alexandrium spp., Chaetoceros convolutus, and C. concavicornis (combined), Dictyocha spp. [Octactis speculum, which prior to 2017 was described as Dictyocha speculum (Chang et al., 2017) was included in the counts], Dinophysis spp. (Phalacroma spp. were included in the counts as this genus is closely related to Dinophysis and share the same negative effects), Heterosigma akashiwo, Margalefidinium fulvescens [prior to 2017 described as Cochlodinium fulvescens (Gómez et al., 2017)], and Rhizosolenia setigera. Pseudo-nitzschia spp. was only enumerated if it was one of the dominant taxa in the sample or if its concentration was significant (i.e., > 50 cells mL^–1) as this genus was not considered harmful in the SoG [during decades of monitoring, ASP toxins almost never exceeded the regulatory limit (McIntyre, 2012)].

Additional Data

Weather data (wind speed, rainfall, and cloud cover) were obtained from the Vancouver weather station (Environment and Climate Change Canada, n.d.a) and Fraser River flow from the archived hydrometric data with daily discharge measured at Hope, BC (Environment and Climate Change Canada, n.d.b).

For fisheries management purposes, HA data were linked to corresponding DFO fishery management areas (Supplementary Materials 1, 2), and HA concentrations were converted to index levels (low, moderate, and high) based on observed effects of toxic algae on net pen-reared Atlantic salmon (Salmo salar) in the SoG and west Coast of Vancouver Island (Haigh et al., 2019). At high levels, algae are likely to cause harm to farmed salmon in net pens. Toxin concentrations in shellfish flesh were obtained from the Canadian Food Inspection Agency (CFIA) through the Access to Information Act.

Statistical Analysis

Pearson Product-Moment Correlation (r) was used to assess the relationships between the dominant HA taxa and physico-chemical variables. Intra-annual relationships were examined using mean monthly data from March to September (2015–2018) averaged over the entire SoG. Due to the possible likelihood of seasonal multicollinearity between physico-chemical variables in our analyses, we used the intra-annual relationships as an indicator of the timing and conditions across seasons wherein specific HABs taxa might be expected to occur. Furthermore, to address the limitations of interpreting these relationships when multicollinearity is present, we removed the seasonal signal to examine inter-annual relationships between HA taxa and environmental variables by using average values for the individual stations during each of the summer months only (June, July, and August), as these were the months when HABs were typically most abundant. HA data from surface samples were used for analysis. All variables were log-transformed to normalize the data in order to meet the assumptions of the statistical analyses, and a standard value of 0.005 was added to all zero data to allow for transformation (O’Brian et al., 2013).

In addition, Redundancy Analysis (RDA), which combines multiple regression with principal component analysis (Abdi and Williams, 2010; Borcard et al., 2018), was used to explore the spatial relationships between HABs taxa and explanatory variables during only the summer months. Mean monthly concentrations of the HABs taxa were averaged across all summer months (i.e., June, July, and August combined) for each station and Hellinger-transformed prior to analysis. All statistical analyses were performed using Sigmaplot version 13.0 and R version 3.5.1 (R Core Team, 2018). Due to their sporadic occurrence (present in < 3% of the samples), Dinophysis spp., M. fulvescens, and N. scintillans were not included in the statistical analysis.

Results

CTD, Secchi Depth, and Additional Environmental Parameters

The general seasonal cycle of surface water temperature, salinity, stratification, and Secchi depth showed similar patterns during all 4 years of this study. Seasonal variations in the surface water temperatures were consistent, with the coldest temperatures generally observed in February (∼8.4°C) and warmest in July (∼17.5°C); however, spring of 2017 was noticeably cooler than the other years (March, April, and May temperatures were ∼7.2, 9.3, and 12.4°C, respectively, in 2017 and ∼8.5, 10.4, and 14.0°C in other years) (Figure 2A). Seasonal variations in surface salinities were also similar across years, with the lowest salinities usually observed in June (∼23.3 g/kg) and the highest in October (∼27.4 g/kg); however, surface salinities in late spring and summer of 2016–2018 were noticeably lower (mean May–August values were ∼24.2 g/kg) compared to that of 2015 (∼27.6 g/kg) (Figure 2B). Stratification, as reflected by the stratification parameter, was generally low in early spring (<1 kg m^–3), gradually increasing from mid-spring (∼1.3 kg m^–3) to the highest values in summer (∼3.4 kg m^–3) before decreasing again. Stratification was considerably higher in the summer of 2018 (∼3.7 kg m^–3) than in 2015–2017 (<3.0 kg m^–3; Figure 2C). Secchi depths were high in February (∼11 m) and low in August (∼5.8 m); Secchi depths were very high (up to 23 m) in April of 2018 (Figure 2D).

FIGURE 2

Figure 2. Water temperature (A) and salinity (B) at 1 m, stratification parameter (C) and Secchi depth (D) in the Strait of Georgia.

In addition to seasonal variations, regional variations in CTD and Secchi depths were apparent. For example, on average between February and October, southern stations adjacent to Juan de Fuca Strait (Victoria stations) had colder surface waters (∼10.4°C), whereas northern areas were slightly warmer (∼14.2°C) (Figure 2A). Spatial variability in salinities was also pronounced, with the southernmost area (Victoria) showing higher salinity (∼30.89 g/kg) compared to the rest of the SoG (25 g/kg). Within the SoG, central and southern surface salinities (∼24.2 g/kg) were significantly lower than northern SoG waters (∼26.1 g/kg) (Figure 2B). The highest stratification parameter values (3.5 kg m^–3) were observed in the eastern SoG (i.e., Steveston, Irvine-Sechelt, and Malaspina areas) owing to outflow from the Fraser River, and the lowest (0.2 kg m^–3) were observed in the Victoria area (Figure 2A). Secchi depths were generally greatest (∼8.0 m) in the northern SoG (Campbell River, Lund, and Powell River areas), with central and southern SoG having shallower values (∼6.7 m). The shallowest Secchi depths (∼3.4 m) were consistently observed in the Steveston area at sites close to the mouth of the Fraser River.

There were substantial interannual variations in weather parameters and Fraser River discharge (Figure 3). Specifically, discharge was high in March and April of 2015 but remained somewhat low during that summer. In 2016, flows reached an early summer peak (in April), and peak flows were only half of those observed in the other 3 years. In 2017 and especially in 2018, high flows occurred in late May and early June.

FIGURE 3

Figure 3. Daily Fraser River flow (A), average monthly precipitation (B), cloud cover (C), and winds (D) in 2015–2018. Data source: Environment and Climate Change Canada.

Nutrients

All surface (0 m) nutrients (nitrates, phosphates, and silicates) decreased in spring of each year, remained low over the summer, and subsequently returned to high concentrations (∼21.2, 1.6, and 48.2 μM, respectively) at the end of each sampling year (Figure 4 and Supplementary Material 3). Surface nitrate concentrations, in particular, were nearly zero over the summer. Phosphate values in summer were low but non-zero and relatively constant (∼0.2–0.5 μM depending on the year), and silicate values were somewhat lower in summer (∼10–30 μM) relative to winter values but were also highly variable. Nutrient concentrations at 20 m had a similar annual cycle with values decreasing during late spring and summer; however, nitrate concentrations did not approach zero as they did in the surface waters and remained ∼8.0 μM during summer months. Overall, mean nutrient concentrations (except for silicate) were higher in 2015 and 2018 compared to 2016 and 2017 (Figure 4).

FIGURE 4

Figure 4. Concentrations (μM) of nitrates (A) and phosphates (B) at the surface; nitrates (C), and phosphates (D) at 20 m depth.

Surface nutrients had very noticeable regional patterns. High mean nitrate and phosphate concentrations occurred at the southern (Victoria) and northern (Campbell River) ends of the SoG (∼20.1 and 12.6 μM, respectively, for nitrates, 1.74 and 1.26 μM for phosphates for the south and north, respectively). Concentrations at stations in the central SoG were less than half (∼4.7 μM for nitrates, 0.62 μM for phosphates) those found in the northern and southern regions. Higher mean silicate values were observed in areas directly affected by river flow (e.g., ∼36.17 μM at Steveston, 35.23 μM at Cowichan Bay, and 35.23 at Campbell River) compared to the rest of SoG (∼24.09 μM).

Nutrients at 20 m also showed strong regional patterns, however, these patterns differed from those observed for surface waters. Notably, high mean values of all nutrients at 20 m occurred in samples from the southern part of the SoG (Victoria area; ∼ 22.53 for nitrates, 1.82 for phosphates, and 35.27 for silicates) compared to those in the rest of the SoG (∼14.08, 1.31, and 27.2, respectively).

Phytoplankton

Interannual Variations

Overall, significant interannual variations were evident in phytoplankton dynamics. Specifically, the peak of the spring bloom in 2015 was recorded several weeks earlier (early March vs. late March/April) than in 2016–2018 (the 2018 bloom was only partially captured by the Citizen Science sampling but apparently occurred mid-April (Chandler et al., 2018), and its concentrations were higher than those observed from 2016 to 2018 (up to > 10,000 diatom cells mL^–1 or an average of ∼2,000 cells mL^–1 compared to ∼6,500 diatom cells mL^–1 or an average of ∼500 cells mL^–1, respectively). Furthermore, the duration of the 2015 spring bloom was shorter compared to the blooms observed in 2016–2018, with average diatom cell counts > 100 cells mL^–1 recorded during three sampling periods in 2015 (early/mid-March and early April), whereas similar cell counts in 2016–2018 were recorded across six sampling periods (from late-March to mid-May).

In terms of community composition, the spring bloom of 2015 was mainly comprised of Skeletonema, whereas the communities in 2016–2018 were comprised of Thalassiosira, Skeletonema, and Chaetoceros. Another notable difference in phytoplankton dynamics between years was that no significant blooms were observed in late summer and early fall of 2015, however, several blooms were observed throughout this same time period in 2016–2018. In addition, the biomass of the majority (∼85%) of the 2015 samples was dominated by diatoms. In contrast, 2016–2018 had a lower number of samples dominated by diatoms (∼65%) and a higher number of samples dominated by dinoflagellates, silicoflagellates, and raphidophytes.

Numerous potentially harmful algal taxa were identified (Table 2). Overall, the most abundant HA from 2015 to 2018 were R. setigera, Dictyocha spp., Alexandrium spp., and H. akashiwo, with each of these taxa occurring in >6% of the surface samples. The least abundant taxa was M. fulvescens occurring only in <1% of the samples, which did not allow for a reliable intra- and interannual or regional comparison. Harmful species of known concern in BC coastal waters (Haigh and Esenkulova, 2011) that were not encountered included Chattonella cf. marina, Gymnodinium mikimotoi, and Pseudochattonella cf. verruculosa. HA abundances and concentrations decreased with depth with the exception of C. convolutus, C. concavicornis, and Dinophysis spp. (Supplementary Material 4). The maximum abundances and concentrations of C. convolutus and C. concavicornis were recorded at 10 and 5 m and for Dinophysis spp. at 5 and 10 m.

TABLE 2

Table 2. Frequency of occurrence (percent of samples containing HA) and maximum concentration (cells mL^–1) recorded in surface samples in different years.

Interannual variations in HA occurrences were substantial (Table 2 and Figure 5). While in 2015 C. convolutus and C. concavicornis were at least twice as common as in 2016–2018 (10% compared to 1–4.5%); all other taxa were more prevalent in 2016–2018. In 2016–2018, R. setigera and Alexandrium spp. were more prevalent compared to 2015 (R. setigera 17.7–20.2% compared to 6.9% in 2015 and Alexandrium spp. 15.7–18.1% compared to 10.7%). Also, Dictyocha spp. was about five times more prevalent in 2016 and 2017 (25.1 and 26.7%) compared to 2015 and 2018, with frequency of occurrence of 3.6 and 5.4%, respectively. Furthermore, in 2018 H. akashiwo was much more prevalent (12.5%) compared to other years (1.1–7.4%), and its maximum concentration (11,000 cells mL^–1) was substantially higher (by two logarithmic scales) compared to 2015–2017 (150 cells mL^–1). Both Dinophysis spp., and N. scintillans were more common in 2018 (5.1 and 2.9%, respectively) than in 2015–2017 (Dinophysis spp. 0.7–1.7% and N. scintillans 0.3–0.7%). Thus, annually, 2015 was the least “eventful” year in terms of HABs occurrences (with the exception of C. convolutus and C. concavicornis), 2016 and 2017 were somewhat comparable, and 2018 exhibited the highest occurrences of HAB taxa (particularly H. akashiwo, Dinophysis spp., and N. scintillans).

FIGURE 5

Figure 5. Concentrations (cell mL^{– 1}) of Rhizosolenia setigera (A), Dictyocha spp. (B), Alexandrium spp. (C), and Heterosigma akashiwo (D) at the surface.

Seasonal Variations

In terms of seasonal variation, the highest abundances of HA in the surface samples were usually recorded from the end of spring to the beginning of fall, with the exception of C. convolutus and C. concavicornis, which were most abundant in March (13.2%) and October (8.8%) (Table 3). Cells of R. setigera, Dictyocha spp., and Alexandrium spp. were mostly observed during May through September, with the highest abundances recorded in August (51.4, 39.7, and 25.0%, respectively). While H. akashiwo was also typically seen in samples from May through September, its highest abundance (17.6%) was recorded in June. Dinophysis spp. typically occurred in May to August samples with its highest abundance (4.3%) observed in June. Cells of N. scintillans were the most common in April and May (2.1% in both months) as well as August and September (1.0 and 1.5%, respectively).

TABLE 3

Table 3. Percent (%) of surface samples containing harmful algae.

Regional Variability

Noticeable regional patterns in HA distribution were evident (Figure 6 and Supplementary Materials 5, 6). R. setigera was generally the most abundant throughout the SoG reaching high concentrations in central SoG (Baynes Sound, Powell River, Lund, Cowichan Bay, and Malaspina Inlet); Dictyocha spp. was the most abundant in the central east regions (Irvine’s Sechelt and Malaspina); Alexandrium spp. in shallow areas of the central west SoG (Cowichan Bay and Ladysmith); H. akashiwo in the central Strait in relatively close proximity to the Fraser and Cowichan rivers (Steveston, Galiano, Malaspina, Irvine’s Sechelt, and Cowichan); C. convolutus and C. concavicornis in northern SoG (Powell River, Campbell River, and Lund); Dinophysis spp. in northern SoG (Malaspina Inlet, Campbell River, and Lund; and N. scintillans in central SoG (Cowichan Bay, Malaspina Inlet, and Irvine’s Sechelt). Overall, most HA were frequently observed in the region between Texada Island and the Sunshine Coast, Mainland BC (specifically, in Malaspina Inlet and Irvine’s Sechelt). In contrast, HABs were the least common in the Juan de Fuca region (Victoria); although this area was sampled for only 1 year, the low frequency of HA occurrence compared to other regions was noticeable (Supplementary Material 6).

FIGURE 6

Figure 6. Average 2015–2018 frequency of occurrence (%) of Rhizosolenia setigera (A), Dictyocha spp. (B), Alexandrium spp. (C), and Heterosigma akashiwo (D) in surface samples. Stations with at least 40 measurements within 4 years (10 per each year) were used.

The most common HA species were associated with variable ranges of salinity, stratification index, and nutrient concentrations (Figure 7). On average, R. setigera was present at high water temperatures and low nutrient concentrations; Dictyocha spp. at low salinities but relatively high nutrient concentrations; Alexandrium spp., at intermediate values; H. akashiwo at the lowest salinities and highest stratification as well as at low nitrate and phosphate levels but very high silicate levels; and C. convolutus and C. concavicornis at low temperatures, high salinity, low stratification, high nitrates and phosphates but low silicates. Although the relatively low occurrence of Dinophysis spp. does not allow for definitive elucidation of environmental preference, this taxon was seen in relatively highly stratified waters (mean stratification parameter 3.19, n = 98) and relatively high nutrient concentrations (mean N = 5.24, P = 0.67, Si = 27.6; n = 166 for each nutrient). Hence most of the HAB taxa had distinct environmental niches.

FIGURE 7

Figure 7. Boxplots for temperature (A), salinity (B), stratification parameter (C), nitrates (D), phosphates (E), and silicates (F) for Rhizosolenia setigera (A), Dictyocha spp. (B), Alexandrium spp. (C), Heterosigma akashiwo (D), Chaetoceros convolutus and C. concavicornis (E). Within each box, the horizontal red line denotes median values and the circle represents mean; boxes extend from the 25th to the 75th percentile; vertical lines denote adjacent values within 1.5 interquartile range of the 25h-75 percentile of each group; red dots denote observations outside the range of adjacent values; the number above each box represents quantity of data points.

In addition to routinely monitored HA, several unusual events were noted. Scheduled samples collected in Cowichan Bay on July 7, 2017 had brown coagulated particles visible to the naked eye, which appeared upon microscopy examination to be a decomposing cyanobacteria bloom. An opportunistic (unpreserved and immediately analyzed) sample collected in bright, turquoise waters in the central SoG (49.52 N, 124.63 W) on August 14, 2017 revealed over 40,000 coccolithophores cell mL^–1. Another opportunistic sample collected directly from a thick, orange bloom on April 17, 2018 at one of the northern SoG stations (LD-2) revealed ∼3,000 cell mL^–1 of N. scintillans. Samples collected at the same location on the scheduled date of April 19, 2018 had only 1 cell mL^–1 of N. scintillans.

Statistical Analyses

Interannual Variability

Pearson Product-Moment Correlation based on summer-averaged values only revealed a significant, positive relationship between H. akashiwo and stratification (r = 0.237) and Si (r = 0.227), and negative relationships with both temperature (r = −0.131) and salinity (r = −0.178) (Table 4). Alexandrium spp. was significantly related to all environmental variables except for Secchi depth, whereas Dictyocha spp. abundance was correlated with both P (r = −0.165) and Si (−0.146), but none of the environmental variables (Table 4). Significant negative relationships were found between Dictyocha spp. and salinity and Secchi depth. R. setigera was only significantly correlated with salinity and Si (r = 0.226 and r = −0.287, respectively). No significant relationships were observed for C. convolutus and C. concavicornis with any of the environmental drivers or nutrients.

TABLE 4

Table 4. Pearson Product-Moment Correlations (r) to examine inter-annual relationships between mean HABs taxa concentrations and environmental drivers and nutrients during summer (June, July, August) using average values for each station in the Strait of Georgia from 2015 to 2018 (n = 259).

The interactive effects of sampling site, physico-chemical variables, and HA concentration is shown in the RDA ordination (Adj r² = 0.11; Figure 8). Together, these variables accounted for only 9% (RDA axis 1) and 3% (RDA axis 2) of the overall variability in HA concentration. The strongest axis, Axis 1, was most correlated with H. akashiwo (r² = 0.124), which was characterized by high Si (r² = 0.369) and associated mainly with sites nearest the mouth of the Fraser River (Steveston, Irvine-Sechelt, Nanaimo Qualicum), as well as a few sites in Campbell River. Additionally, Axis 1 was significantly correlated with R. setigera (r² = −0.194) and associated with sites exhibiting high salinity (r² = 0.369), such as those in the Baynes Sounds and Lund regions. RDA ordination grouped by year clearly shows that the majority of Heterosigma blooms occurred in 2018.

FIGURE 8

Figure 8. RDA triplots showing the relationships between HABs taxa and physico-chemical parameters averaged across all summers from 2015 to 2017 in the Strait of Georgia grouped by station (A) and year (B). RDA plots display 9% (RDA axis 1) and 3% (RDA axis 2) of HABs abundance variance. Eigenvalues are 0.045 and 0.016 for the first and second axes, respectively. Vectors show the relationships between physico-chemical variables and the two ordination axes. Vectors pointing in the same direction indicate variables that are highly correlated, whereas vectors pointing in opposing directions indicate the variables are negatively correlated. (Note: only the axes with the two highest variances are shown).

Intra-Annual Variability

Pearson Product-Moment Correlation revealed several seasonal links to physical and chemical variables, thus providing key information regarding the time of year, and resulting environmental conditions, most strongly related to the occurrence of certain HABs taxa (Table 5). R. setigera, Dictyocha spp., Alexandrium spp., and H. akashiwo had significant positive relationships with temperature (range of r-values from 0.498 to 0.753) and stratification (r = 0.443 to r = 0.709) as well as significant negative relationships with Secchi depths (r = −0.645 to r = −0.406) and nitrates (r = −0.651 to r = −0.408). In contrast, C. convolutus and C. concavicornis showed significant negative relationships with temperature and stratification (r = −0.509 and r = −0.537, respectively) and positive relationships with Secchi and nitrates (r = 0.429 and r = 0.609, respectively). R. setigera, Alexandrium spp., and H. akashiwo were negatively correlated with rainfall (r = −0.417, −0.505, and −0.466, respectively). Both Alexandrium spp., and H. akashiwo were negatively correlated with cloud cover (r = −0.653, r = −0.379) and phosphates (r = −0.557, r = −0.383). Dictyocha spp. was the only taxa that had a significant relationship with salinity (r = −0.441). No significant associations between HA and mean monthly silicate concentrations, wind speed, and Fraser River discharge were found.

TABLE 5

Table 5. Examination of intra-annual relationships using Pearson Product Moment Correlations (r) between mean monthly concentrations of HABs taxa and various physical and chemical variables from March to September 2015–2018 averaged over the entire SoG (n = 28).

Harmful Algae and Fisheries Management

The only HA species that reached high levels in terms of potential toxic effects on salmonids during the study period were H. akashiwo and Dictyocha spp. Heterosigma concentrations were high (>500 cell mL^–1) across the entire SoG, except in area 17 where its concentration was moderate (10–499 cell mL^–1) at the beginning of June of 2018 (Table 6). Heterosigma concentrations were high in area 16 at the end of May, June, and July 2018, end of May at area 29, and mid-August of 2018 at area 17. Therefore, the summer of 2018 had the most negative conditions during the study period with respect to prevalence of HABs toxic to finfish. High levels of Dictyocha spp. (>300 cell mL^–1) were occasionally detected in two areas: area 16 at the end of June and beginning of July 2016 and areas 16 and 18 in August of 2017 (Supplementary Material 7). Both taxa, H. akashiwo and Dictyocha spp. did not reach high or moderate levels anywhere in the SoG during 2015.

TABLE 6

Table 6. Heterosigma akashiwo levels (level 1 ≤ 10 cell mL^–1; level 2 = 10–499 cell mL^–1, level 3 ≥ 500 cell mL^–1) based on maximum cell concentration encountered in DFO management areas for each year 2015–2018.

The occurrence of HA coincided with negative effects in SoG aquaculture. Higher toxin concentrations of PSP (PSP-total) and DSP (TOX-DSP-LC) in shellfish flesh were detected by CFIA in years when Alexandrium spp. and Dinophysis spp. were more prevalent in the Citizen Science samples (Table 7). Mass mortalities of aquacultured salmon were reported in years with high concentrations of Heterosigma and Dictyocha. Finfish aquaculture sites exist in the northern SoG DFO management area 16 (DFO, 2019a). Four reports of salmon mortalities in 2016 (June 14, July 1 and 27, August 29) overlapped with high Dictyocha levels (Supplementary Material 7) and five mortality reports in 2018 (at 2 sites on June 2, June 6, June 12, and August 19) overlapped with high Heterosigma levels (Table 6); no mass mortalities were reported in 2015 and 2017 (DFO, 2019b) and no algal species harmful to salmon reached high levels in these years.

TABLE 7

Table 7. Alexandrium and Dinophysis frequency of occurrence (%) and maximum toxin concentrations (data source: Canadian Food Inspection Agency) in February–October, Strait of Georgia.

Discussion

This study presents the first multiyear, high-resolution observations of harmful algal dynamics in the Strait of Georgia. Our results illustrate clear spatial and temporal patterns in intensity and distribution of harmful algal blooms and provide some of the first information on both physical and chemical environmental preferences for the most commonly recorded taxa.

Strait of Georgia Environment

The environmental parameters and HAB dynamics measured during this study captured the transition from conditions affected by a very strong El Niño event (same class as the strongest El Niños of twentieth century; Santoso et al., 2017) in 2015 and early 2016 to a weak La Niña from the end of 2016 throughout 2017, and then transitioning to a weak El Niño in 2018 (NOAA, n.d.). As a result, surface water temperatures in the SoG, overall, in 2015 and 2016 were slightly above climatological means, whereas surface temperatures in 2017 and 2018 were below climatological means over this period (Pawlowicz et al., 2020). In addition, Fraser River flow is indirectly affected by El Niño as air temperature affects snow accumulation and subsequent snowmelt (Thorne and Woo, 2011). Given that the surface salinity of the SoG, and thus our stratification parameter, is greatly influenced by freshwater inflow on monthly scales (Pawlowicz et al., 2007), Fraser River inflow is the most important factor affecting surface salinity in the Strait. Along with its effects on stratification, the Fraser River plume is also quite turbid, especially between April and June (Phillips and Costa, 2017), thus contributing to the relatively shallow Secchi depths seen in the southern Strait. In contrast, Secchi depths in the northern Strait tend to reflect the magnitude of primary production in this region, with clear water in winter and shallower Secchi depths in summer. As a consequence of this variation in discharge, sea surface salinity, stratification and nutrient dynamics in the Strait were affected. The most pronounced summer signals included saltier and not completely nitrate depleted waters in 2015 (especially in the northern Strait) and highly stratified waters in 2018, with almost a complete absence of nitrate in the surface waters.

Phytoplankton Dynamics

Annual variations in phytoplankton dynamics were associated with differences in physical and chemical water properties. Results of the present study show that the spring bloom of 2015 occurred earlier than in 2016–2018, in agreement with studies based on in situ chlorophyll, HPLC data and satellite imagery (Chandler et al., 2018; Gower and King, 2018; Suchy et al., 2019; Del Bel Belluz et al., 2021), as well as detailed observations available from monitoring systems installed on ferries that cross the area (Wang et al., 2019). The early spring bloom of 2015 was associated with the anomalous conditions caused by a strong El Niño, including warm water temperatures, high February water stratification, weak winds (Environment and Climate Change Canada, 2018), and an early start of the Fraser River freshet (Figure 3). Although spring 2016 conditions were still influenced by El Niño, waters in the SoG were slightly cooler in February and March, less stratified, and affected by stronger winds (Environment and Climate Change Canada, 2018; Suchy et al., 2019); the latter of which is known to disrupt spring bloom development in the SoG (Yin et al., 1997; Collins et al., 2009). Altogether, these conditions led to a later start for the spring bloom in 2016 than was observed in 2015. Spring bloom timing in 2017 and 2018 (i.e., years not affected by strong ENSO events) were comparable to 2016, and altogether could be considered “typical” bloom years. Although the present study reports strong signals in bloom timing variability, these results are based on samples collected at 2-week intervals. Given the known importance of the spring bloom timing for the SoG ecosystem (El-Sabaawi et al., 2009; Schweigert et al., 2013; Riche et al., 2014), a finer time-scale for data (e.g., continuous chl data, satellite imagery) collection may be more appropriate for comprehensive spring bloom characterization.

Another notable observation from this study is that the spring bloom of 2015 was comprised mainly of Skeletonema as opposed to the combination of diatom taxa (Thalassiosira, Skeletonema, and Chaetoceros) observed during 2016–2018, which is more commonly reported for the SoG (Parsons et al., 1969; Huntley and Hobson, 1978; Harrison et al., 1983; Taylor et al., 1994; Hobson and McQuoid, 1997; McQuoid and Hobson, 1997). This 2015 spring bloom composition was unusual as, to our knowledge, a Skeletonema-dominated spring bloom in SoG was only previously reported in Howe Sound in 1974 (Stockner et al., 1979). The Skeletonema bloom in Howe Sound was initially linked to low irradiance (Stockner et al., 1979) as well as effluent from pulp mills (Stockner and Cliff, 1979), while Harrison et al. (1983) later suggested a possible causal link with both irradiance and temperature. In our study, the 2015 spring bloom developed in mid-February, wherein the number of hours of sunlight was about 1 h shorter compared to the later (mid-to-late March) blooms in 2016–2018 (Environment Canada, 2019a). Haigh et al. (1992) noted that in Sechelt Inlet, in 1989 and 1990, peaks of Skeletonema during mixed blooms were related to high N:P ratio and suggested that taxa dominance can be a response to optimal nutrient conditions. In this study, during February 2015, high N and P concentrations were observed (21.66 and 1.75 μM, respectively), thereby suggesting that the N:P ratio may also be an important contributing factor for a Skeletonema-dominated bloom. Future studies are needed to determine if Skeletonema-dominated spring blooms are rare for the SoG, and to elucidate the conditions driving these blooms and their ecological implications.

In addition to the variability observed for spring blooms during our study, late summer algae blooms had strong interannual variations, although 2015 was uneventful compared to 2016–2018. However, the reasons for low phytoplankton cell concentrations at the surface in the summer of 2015 remain unclear. For example, summer 2015 had stronger winds compared to 2016 and 2017 (Figure 3), and unlike spring winds, summer wind events generally have a positive effect on bloom development in the SoG because they bring growth-limiting nitrogen from deeper waters to depleted surface waters (Harrison et al., 1983; Harrison and Mackas, 2014). Indeed, nitrogen was replete at the surface in summer 2015, but no significant blooms were observed, suggesting that limiting factors other than nitrogen, or perhaps biological top-down controls (e.g., zooplankton grazing), were in place. In contrast, the summers of 2016–2018 had weaker winds compared to 2015, much higher freshwater input from heavier rains (Environment Canada, 2019b) and subsequent Fraser River discharge (Figure 3). Therefore, these conditions were more favorable for bloom development, and a number of blooms, including HABs, were observed. A detailed description of harmful algal taxa and environmental parameters are discussed below.

Harmful Algae Variability in Relation to Physico-Chemical Drivers

Seasonal and Interannual Variability in Harmful Algae

The most common HA taxa found during this study were R. setigera, Dictyocha spp., Alexandrium spp., H. akashiwo, and C. convolutus with C. concavicornis. All taxa, with the exception of Chaetoceros, were prevalent during summer months (when temperatures and stratification are high and Secchi depths as well as nitrates are low) while C. convolutus with C. concavicornis were more prevalent during spring/fall. The pennate diatom R. setigera was common and often abundant at the end of summer and early fall, consistent with previous studies (Shim, 1976; Haigh and Taylor, 1991; Hobson and McQuoid, 1997; McQuoid and Hobson, 1997; Cassis et al., 2011). Similar to the previously reported low abundance of R. setigera in SoG inlets during the El Niño summers of 1986 and 1987 (Sancetta, 1989), our results showed that this species was less abundant during the 2015 El Niño year, which exhibited saltier waters in summer. The observed negative relationship with silica was not unexpected as diatom blooms deplete silicates (Cupp, 1943). Seasonally, R. setigera had a negative relationship to rainfall, indicating that high abundances occurred during months with low precipitation levels (e.g., July and August).

The silicoflagellate Dictyocha spp. was abundant throughout the entire sampling season, and predominately formed blooms of the naked form during the summer, which is a pattern previously reported in the SoG (Haigh et al., 1992, 2019). Notably, Dictyocha was the only taxa whose abundance was related to nutrients, exhibiting a negative relationship with both phosphates and silicates, with prominent blooms developing in the summers of 2016 and 2017 when nutrients were especially low. These results provide the first evidence of significant differences in the environmental drivers of Dictyocha and Heterosigma, whose environmental niches were previously thought to be similar in BC waters (Haigh et al., 1992, 2019; Horner, 2002). Consistent with previous studies (Valkenburg and Norris, 1970; Haigh et al., 1992; Henriksen et al., 1993), Dictyocha was negatively related to salinity.

Dinoflagellates from the genus Alexandrium were abundant but did not result in heavy blooms (maximum concentration was only 18 cells mL^–1), a pattern common for some coastal waters (Anderson, 1997, 1998). Interannually, while Alexandrium in the SoG was more likely to thrive in summer months characterized by warm stratified waters as has been observed in other northwest Pacific regions (Anderson, 1997; Moore et al., 2009; Tobin et al., 2019), its abundance was negatively affected by the 2015 El Niño conditions and positively by the 2017 La Niña conditions.

The most significant fish-killing HA in the SoG, H. akashiwo (Whyte et al., 1999; Haigh and Esenkulova, 2014), displayed considerable interannual variability. Despite blooming (>100 cell mL^–1) in the SoG almost every year since the 1990s (Haigh et al., 2004; Brown et al., 2018), there were no Heterosigma blooms in 2015 and 2017, with particularly low cell counts occurring in 2015. In contrast, strong Heterosigma blooms did occur in 2018, and were positively correlated with both stratification and silica. Even though 2017 and 2018 had comparable water temperatures and high Fraser River spring discharge rates, the freshet during 2018 occurred several weeks earlier, and the waters were highly stratified. These results support the notion that the timing of the Fraser River freshet, and resulting stratification, are key factors for Heterosigma bloom development in the SoG (e.g., Haigh and Taylor, 1990; Taylor and Harrison, 2002). Seasonally, Heterosigma was negatively linked to phosphates, rainfall, and cloud cover, suggesting higher abundances may occur in nutrient-poor, drier summer months. To our knowledge, the relationship between Heterosigma abundance and cloud cover has not previously been suggested for BC coastal waters. Further investigation of Heterosigma and light in SoG is required given that this alga can grow at low light intensities (Smayda, 1998) due to strong light-harvesting capabilities (Zhang et al., 2019). Thus, it is important to determine if our results reveal a true relationship between Heterosigma and cloud cover, or if this finding is a result of multicollinearity amongst the environmental variables. In addition, longer time series including additional variables that are known to be important for Heterosigma (e.g., nutrients, minerals, vitamins, bottom water temperature, weather events (Smayda, 1998; Imai and Itakura, 1999; Granéli, 2006) would be crucial in future studies of the physico-chemical drivers of this species.

Spatial Variability of Harmful Algae

This study demonstrated several key spatial differences in HA distribution, thus reflecting the diverse environmental niches known to occur in the SoG ecosystem, which often results in a heterogeneous and patchy phytoplankton community (Haigh and Taylor, 1991). For example, Alexandrium was common at near-shore sites, whereas R. setigera was more abundant in deeper, saltier waters—a finding consistent with previous HA studies in the region (Sancetta, 1989; Haigh and Taylor, 1990; Horner et al., 2011). Alexandrium frequently occurred in shallow, restricted embayments and inlets in central SoG (Cowichan Bay, Ladysmith, Irvine’s Sechelt, and Malaspina). Furthermore, most (>65%) of the highest concentrations (5–18 cell mL^–1) of Alexandrium were seen in the western central SoG (Cowichan Bay, Ladysmith), where sampling sites were relatively close to the shore and the depth was under 100 m. In many cases, Alexandrium blooms near the coast are “point-source” events with direct coupling between cysts (stage of Alexandrium lifecycle) in the sediments and conditions of overlying waters (Anderson, 1997). Cowichan Bay, in particular, appeared to have a prominent “seedbed” of cysts. Future studies on cyst distribution and abundance in sediments (e.g., Pospelova et al., 2010; Horner et al., 2011; Tobin and Horner, 2011) at the Alexandrium hot-spots identified here can further improve our understanding of local bloom dynamics.

Another notable spatial difference was that the most southerly and northerly parts of the Strait (Victoria and Campbell River areas, respectively), which generally showed nitrate-rich, turbulent waters (as evidenced by low stratification) with low phytoplankton cell concentrations. These areas are known to have intense vertical circulation and fast surface currents (Thomson, 1981; LeBlond, 1983) and thus exhibit low phytoplankton standing stock despite high nutrient concentrations (Lewis, 1978; Masson and Peña, 2009). Although the HA concentrations at the southern end of the SoG were always low, this was not the case for the northern area, which consistently had the highest levels of C. convolutus and C. concavicornis, taxa adapted for survival in high nutrient, salty, turbulent waters (Albright et al., 1992; Anderson et al., 1996). In contrast, the densest blooms of Dictyocha and H. akashiwo occurred in the central SoG, close to river inputs (Fraser and Cowichan rivers), where waters were silicate-rich and highly stratified, conditions that are known to be associated with silicoflagellates and raphidophytes (Henriksen et al., 1993; Taylor and Horner, 1994; Horner, 2002). Lastly, Malaspina Inlet in the SoG historically had very high primary production (Stockner et al., 1979), conditions prone to HABs development (Anderson et al., 2002; Heisler et al., 2008). Our results indicated that this area tended to have heterogeneously moderate to highly stratified waters with high phytoplankton cell concentrations and frequent HABs.

Other Algae

During 2015–2018, several HA appeared too infrequently to establish statistically significant relationships with environmental conditions. Among routinely monitored HA, Dinophysis spp., N. scintillans, and M. fulvescens were seen in less than 2.5% of the samples, and several events (e.g., Pseudo-nitzschia and coccolithaphore blooms) were recorded only once during these 4 years. Despite their sporadic nature, these observations are important and emphasize the need for long-term data series in order to capture these events. For instance, similar to other studies, Dinophysis spp. were fairly common but not abundant in the SoG (Haigh et al., 1992; Taylor and Haigh, 1996; Taylor and Harrison, 2002), occurring most frequently at stations in the northern SoG (i.e., Malaspina and Lund) from mid-summer onwards as has been previously reported (Taylor and Haigh, 1996; Esenkulova and Haigh, 2012; Haigh et al., 2014; Pospelova et al., 2018). A notable result not previously reported from the SoG was the high interannual variability in Dinophysis, which was at least three times more abundant in 2018 compared to 2015–2017. This taxon appeared in moderately stratified waters with high nutrient concentrations; however, similar conditions frequently occurred where Dinophysis was not abundant. Therefore, this observation suggests that other factors, perhaps biological, may have played an important role in the dramatic variability in Dinophysis abundance.

In addition, widespread, vivid-orange blooms of N. scintillans occurred in April and May of 2018 (Esenkulova and Pearsall, 2019). While relatively high concentrations (up to 14 cell mL^–1) were recorded in routine samples in 2018, opportunistic sampling conducted directly from the visible bloom patch 2 days prior to scheduled sampling revealed extremely high concentrations (3,000 cell mL^–1) of N. scintillans. This observation reinforces the fleeting nature of blooms and the importance of high-frequency sampling. Relatively low Noctiluca occurrence and corresponding environmental parameters do not allow for an environmental niche characterization, however, an important observation was that thick blooms of Noctiluca occurred in 2018, the same year wherein Dinophysis abundance was high. Given that both taxa are phagotrophic, prey availability could have been a critical factor for their high abundances in 2018.

Another HA, Pseudo-nitzschia, is generally not considered harmful in the SoG (McIntyre et al., 2013) and was enumerated only when it was dominant in a sample or when its concentration exceeded 50 cell mL^–1. Only one thick bloom (1,000–4,500 cell mL^–1) of Pseudo-nitzschia was observed, occurring in the northern SoG (Lund and Powell River) at the end of August and beginning of September 2018. Simultaneously, R. setigera concentrations up to 300 cell L^–1 were observed in the same samples, suggesting conditions were favorable for bloom development of pennate diatoms. In striking contrast to the extreme, months-long Pseudo-nitzschia blooms occurring along the majority of the West Coast of North America during the spring and summer of 2015 (Du et al., 2016; McCabe et al., 2016), no blooms of Pseudo-nitzschia were observed in the SoG during the same period. In January and February 2015, Pseudo-nitzschia was dominant in southern SoG (Victoria) with low concentrations of < 10 cell mL^–1, though later spring samples were dominated by Skeletonema. It is likely that regional factors in the SoG in 2015, including the semi-enclosed nature of the SoG and the high nutrient conditions known to be important in Skeletonema bloom development (Smayda, 1957; Bates and Trainer, 2006; Granéli, 2006), allowed Skeletonema to outcompete Pseudo-nitzschia.

Lastly, most areas of the SoG turned unusually bright turquoise in July and August of 2016 (Chandler et al., 2017), and unpreserved samples revealed high concentrations (40,000 cell mL^–1) of coccolithophores. Coccolithophore blooms are typical for the west coast of Vancouver Island but are uncommon within the SoG (Haigh et al., 2015). Future studies are thus needed to determine the cause of this unusual event. While blooms of the dinoflagellate M. fulvescens do occur in BC (Whyte et al., 2001), none were observed during the study period, and only a few rare cells (<3 cell mL^–1) were detected in samples taken in summer from the northern and central SoG (Steveston, Malaspina Inlet, Campbell River, and Cowichan Bay). Longer datasets will likely reveal additional uncommon HA events and have the potential to advance our understanding of the linkages between HABs and changing oceanographic conditions in the SoG (Wells et al., 2015; Trainer et al., 2019).

Harmful Algae and Fisheries Management in the Strait of Georgia

The Marine Biotoxins program in Canada includes biotoxins monitoring in shellfish flesh by the CFIA and shellfish harvest openings/closures reinforced by DFO.¹ Phytoplankton monitoring is not part of the national program, although its need is recognized (McIntyre et al., 2013). This study revealed that higher toxin levels accumulated in shellfish flesh in years when toxic algae were more abundant. An in-depth study that spatially corresponds phytoplankton and shellfish sampling sites would help to establish direct, local links between algae and toxin concentrations. An ongoing HA monitoring program could also provide an early warning of biotoxin contamination in SoG.

High levels of toxic algae observed during this study corresponded with industry-reported mortalities from HABs at nearby salmon farms. Four mortality events were reported in 2016 when Dictyocha reached high levels for two consecutive sampling periods, and five events were reported in 2018 when Heterosigma levels reached high and moderate levels for five consecutive sampling periods. Fish-killing events of 2018 were especially severe, with economic losses of ∼ USD 3 million (Editors, 2018; Robinson, 2018). Current observations of fish-killing events associated with HABs in the SoG are consistent with studies from previous years (e.g., Taylor and Horner, 1994; Whyte et al., 1997; Rensel and Whyte, 2004; Haigh and Esenkulova, 2014). The fact that HABs severely impact local finfish aquaculture, coupled with a suggestion that HABs could affect wild SoG salmon (Rensel et al., 2010; Chittenden et al., 2018), make future studies using high resolution, long-term HAB time series, like the one presented in this study, especially important for evaluating the potential effect of HABs on wild salmon stocks.

Assessing the impacts of HABs on marine ecosystem health is a complex issue (Hoagland et al., 2002; Stauffer et al., 2019; Brown et al., 2020; McKenzie et al., 2020). HA recorded during this study included not only toxic species (Alexandrium spp., Dictyocha spp., Dinophysis spp., and H. akashiwo), which arguably cause the most direct damage and are the easiest to connect to negative outcomes, but also species that are non-toxic but cause other effects. For example, C. convolutus, C. concavicornis, and R. setigera are mechanically harmful to finfish, N. scintillans disrupts classic energy transfer, and any algae that produces high-density blooms decrease visibility and reduce dissolved oxygen. This study revealed frequent and sometimes highly concentrated toxic and non-toxic HA in the SoG and points to a likelihood of various negative impacts to wild and aquacultured fisheries, marine ecosystem health, and food safety. Without appropriate studies, these impacts are likely to remain unknown and overlooked. Because DFO’s Sustainable Fisheries Framework advocates an implementation of a precautionary approach (DFO, 2009), HAB impacts must be studied and assessed to evaluate the potential role of HA as a stressor on the SoG ecosystem to reduce uncertainty and ensure an ecosystem approach to local fisheries management.

Study Limitations

The described patterns of HA spatial and temporal abundance and their correspondence to specific environmental niches contribute to the current understanding of how phytoplankton dynamics vary in different regions of the Strait (e.g., Masson and Peña, 2009; Suchy et al., 2019), and can have significant implications for higher trophic levels, including commercially important salmon species. However, the short length of our time series (only a few years) is one of the major limitations in defining ecologically meaningful statistical relationships. While this study shows clear correlations between physico-chemical variables and HABs taxa, strong collinearity amongst variables may be present. Thus, we interpret these results with caution and suggest that they represent only a first estimate of the “ideal conditions” under which certain HAB taxa could persist. Our separate analyses for intra-annual (all months considered) and interannual (summer months only) help to further elucidate true relationships between HA taxa and their environmental drivers. However, despite this and various other challenges related to predicting HABs (e.g., Anderson et al., 1996; Smayda and Reynolds, 2001; Glibert et al., 2010), newer models, combined with sufficient data, have led to recent improvements (McGillicuddy et al., 2011; Franks, 2018; Glibert et al., 2018), so we believe it is worthwhile to present these results now. Long-term, high-resolution citizen science data, integrated with other available datasets, can clearly contribute more to future SoG HAB models and ecosystem-based research.

Citizen Science Programs

This work would not have been possible without the multi-year, high resolution observations of HA and associated environmental conditions in the SoG by the PSF Citizen Science Program. Although monthly sampling occurred for several years during the late 1960s at a small number of stations in the SoG, comprehensive monitoring of the whole Strait has been carried out since then, at best, only 4 times a year. Our dataset, which involves many thousands of CTD profiles and water samples, provides information at a temporal and spatial scale never accomplished before. This program is nationally unique, but, similar to some other citizen science projects, has been successful for advancing scientific knowledge and providing vast quantities of data (Bonney et al., 2009; Silvertown, 2009). The amount of effort involved in the HAB research reported here should not be underestimated. Upwards of 30 “citizens” each contributed at least 20 full days a year to collect samples, backed by several full-time technicians and a project manager who provided real-time support and periodically visited all teams. Over a dozen lab analysts contributed nutrient and chlorophyll analysis, plankton enumeration, CTD data processing and archiving, and quality control. Funding provided by PSF also covered the capital cost of outfitting boats with wire outrigger winches, as well as the purchase and maintenance of CTDs, bottles, Secchi disks, nets, and other equipment to outfit each patrol, at a cost roughly equivalent to 1 year of operating costs. At times, the heavily used equipment required expensive repairs, or even replacement (e.g., CTDs are now being replaced after 6 years). Crucial to the success of such programs is the ongoing involvement of interested scientists (as in the authors of this paper). Their professional input and oversight ensure that the data gathered is more than an archive, and of a high enough standard to be useful in addressing important research questions.

Conclusion

This is the first work producing multi-year, high-resolution records of HA distribution and corresponding environmental parameters in the SoG, made possible as a result of a well-organized citizen science program. Pronounced interannual variations of phytoplankton dynamics and especially late summer blooms were documented. The spring bloom of 2015 occurred very early, most likely due to strong El Niño conditions, and was unusually comprised of Skeletonema. Summer phytoplankton cell concentrations were low and there were no significant HABs despite the fact that surface waters were not nutrient-depleted. The spring blooms in 2016–2018 occurred later and lasted longer than in 2015 and were comprised of a mixture of diatoms (Thalassiosira, Skeletonema, and Chaetoceros). Summer cell concentrations in 2016–2018 were higher, with substantially higher abundances of dinoflagellates, silicoflagellates, and raphidophytes compared to summer 2015, and blooms of several HAB taxa were observed. Dictyocha spp. blooms occurred in 2016 and 2017, and heavy H. akashiwo and N. scintillans blooms were documented in 2018.

Our results provide new information about the spatio-temporal distribution, ecological niches, and statistically significant relationships between HA and the environmental parameters. Alexandrium was more common in shallow, near-shore regions, whereas R. setigera was more common in deeper, saltier waters. In contrast, Chaetoceros convolutus and concavicornis consistently had high abundances in the high nutrient, and turbulent, northern SoG. Although Dictyocha and H. akashiwo both tended to occur in the central SoG near freshwater inputs, Dictyocha was the only taxa whose summer abundance was related to nutrients, with prominent blooms occurring during years when nutrient concentrations were especially low. H. akashiwo blooms, on the other hand, were influenced by stratification patterns.

The interannual variability in phytoplankton dynamics observed in this study further highlights the need for long-term data collection to better characterize HAB events and establish which environmental factors, for specific taxa, are the most important for bloom development. Data collected during this study has the potential to be combined with higher trophic level data (e.g., zooplankton and fish) to further analyze interactions among physical oceanographic parameters, food web dynamics, and potential outcomes for fish stocks. Together, this information can aid in our understanding of the relationships between SoG conditions and the survival and growth of the juvenile fish that are of cultural and economic importance in the region.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Nutrients in the Strait of Georgia: https://soggy2.zoology.ubc.ca/geonetwork/srv/eng/catalog.search#/metadata/dc72987d-2eb0-4b36-a6bf-62d04 272bf29. Phytoplankton in the Strait of Georgia: https://soggy2. zoology.ubc.ca/geonetwork/srv/eng/catalog.search#/metadata/73 6c4721-db8e-4116-a5f8-d9084126b135. Secchi Depths in the Strait of Georgia: https://soggy2.zoology.ubc.ca/geonetwork/srv/eng/catalog.search#/metadata/7c9e31de-f5a3-439e-bf31-aacab1fadf11.

Author Contributions

IP: project design. SE, RP, and KS: data analysis. KS: statistical analysis. SE and KS: writing. SE, KS, RP, MC, and IP: manuscript revision and approval. All authors contributed to the article and approved the submitted version.

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.

Funding

The Citizen Science Program was funded through PSF’s Salish Sea Marine Survival Project, which was supported through the Pacific Salmon Commission, DFO, Ocean Networks Canada, and private donations. This is publication number 60 of the Salish Sea Marine Survival Project (marinesurvivalproject.org).

Acknowledgments

We thank the many individuals who contributed to the inception, development, and management of the Citizen Science program of the Strait of Georgia. We hereby thank all citizen science crews for the numerous days spent at sea and their continued enthusiasm for the Program. Ben Skinner, Trent Suzuki, and Brianne Boufford are thanked for plotting the environmental and HA time series and maps.

Supplementary Material

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

Supplementary Material 1 | Sampling areas and stations with coordinates, number of phytoplankton samples per year.

Supplementary Material 2 | DFO fishery management areas in the Strait of Georgia. Map adapted from the official website of the Government of Canada.

Supplementary Material 3 | Concentrations (μM) of silicates at the surface (A) and at 20 m depths (B).

Supplementary Material 4 | Frequency of occurrence (% of samples where algae were observed) and maximum cell concentration (cell mL^–1) at various depth. Note: only complete sets with data for all 4 depths at a given date at a given location were used.

Supplementary Material 5 | Average 2015–2018 frequency of C. convolutus and C. concavicornis occurrence (%) in the surface samples. Stations with at least 40 measurements within 4 years (10 per each year) were used.

Supplementary Material 6 | Average annual frequency of occurrence (%) of Rhizosolenia setigera (A), Dictyocha spp. (B), Alexandrium spp. (C), Heterosigma akashiwo (D), and C. convolutus and C. concavicornis (E) in the surface samples in 2015 (A–D1), 2016 (A–D2), 2017 (A–D3), and 2018 (A–D4). Stations with at least 10 measurements per each years were used. Note that break levels for charts (E) are different than breaks for (A–D).

Supplementary Material 7 | Dictyocha spp. levels (level 1 ≤ 10 cell mL^–1; level 2 = 10–299 cell mL^–1, level 3 ≥ 300 cell mL^–1) based on maximum cell concentration encountered in DFO management areas. Note: Samples collected within 4 days of the scheduled sampling date were used. Area 19 is not included as it was monitored for 1 year and harmful levels never exceeded low levels.

Footnotes

^ https://www.pac.dfo-mpo.gc.ca

References

Abdi, H., and Williams, L. J. (2010). Principal component analysis. Wiley Interdiscip. Rev. Comput. Stat. 2, 433–459. doi: 10.1002/wics.101