Predicting the Effects of Climate Change on the Occurrence of the Toxic Dinoflagellate Alexandrium catenella Along Canada’s East Coast

Introduction

Harmful algal blooms (HABs), as defined by the UNESCO, designate proliferations of microalgae that are detrimental to public health, local marine fauna, fisheries and the aquaculture industry. The worldwide socio-economic threat posed by HABs is often linked to marine dinoflagellates, such as Alexandrium spp. (Anderson et al., 2012). Formation of long-lasting cysts at the end of the blooming season allows for recurrent Alexandrium spp. blooms. In this study we focus on Alexandrium catenella (previously known as Alexandrium tamarense), a widespread alga that produces paralytic shellfish toxins (PSTs). PSTs bioaccumulate in tissues of shellfish and zooplankton, making filter-feeding invertebrates highly neurotoxic to marine fish, mammals and birds when ingested (Shumway, 1990; Bricelj and Shumway, 1998; Cembella et al., 2002; Landsberg, 2009; Starr et al., 2017). Bioaccumulation of PSTs also prevents human consumption of contaminated shellfish, such as the economically valuable blue mussel Mytilus edulis (Etheridge, 2010). Bloom events often lead to mussel farm and aquaculture closures in affected areas (Blasco et al., 2003; Díaz et al., 2019). Aquaculture sites along Canada’s East Coast are affected by A. catenella when favorable environmental conditions occur (Bates et al., 2020; McKenzie et al., 2020).

It has been reported that A. catenella excystment, growth and bloom occurrence are largely controlled by oceanographic conditions (McGillicuddy et al., 2011; Condie et al., 2019). High sea surface temperature and low sea surface salinity have been identified as favorable to A. catenella (Weise et al., 2002; Fauchot et al., 2005; Ní Rathaille and Raine, 2011; Starr et al., 2017). Sea surface salinity decreases associated with high runoff were also linked to higher stratification and the presence of humic substances, which were found to stimulate A. catenella growth (Gagnon et al., 2005). High wind speeds > 20 km h^–1 associated with mixing of the water column inhibit A. catenella bloom events (Weise et al., 2002; Fauchot et al., 2005). Other variables that influence A. catenella life history are nutrients (Fauchot et al., 2005; Collos et al., 2007; McGillicuddy et al., 2011), pH and CO₂ (Raven et al., 2020). Changes in these environmental factors are expected under a warming climate (Long et al., 2016; Lavoie et al., 2017; IPCC, 2019).

Alexandrium catenella hotspots are widespread along the Canadian East Coast. Known hot spots include the Estuary and Gulf of St. Lawrence (EGSL), the Bay of Fundy (BoF) and the Gulf of Maine (GoM), where fisheries and aquaculture closures are often linked to A. catenella toxins (McGillicuddy et al., 2011; Bates et al., 2020; McKenzie et al., 2020). Furthermore, toxic bloom initiation is frequently associated with increased freshwater input in these regions (Weise et al., 2002; Keafer et al., 2005; Starr et al., 2017). On the Canadian East Coast, climate change is expected to alter sea surface temperature as well as the hydrological cycle, including an earlier spring freshet, and a more intense runoff overall (Long et al., 2016; IPCC, 2013). This new input of fresher, warmer water could have important repercussions on A. catenella occurrence along the coast. Effects of climate change on oceanographic conditions are expected to alter A. catenella bloom frequency, duration and amplitude (Hallegraeff, 2010; Townhill et al., 2018; Trainer et al., 2020) and are already causing range shifts in some phytoplankton species (Barton et al., 2016; Gobler et al., 2017). Furthermore, higher CO₂ concentrations as well as higher temperatures in temperate regions might increase harmful algae growth rates and toxicity (Brandenburg et al., 2019; Roggatz et al., 2019). In this perspective, predicting bloom events is becoming increasingly important in order to better understand climate change effects on A. catenella bloom frequency, spatial extent and temporal variability.

To date, various simple modeling approaches, such as regression-based models, have been used to predict HABs worldwide (Lane et al., 2009; Anderson et al., 2010; Raine et al., 2010; Cusack et al., 2016; Dabrowski et al., 2016). While valuable, simple models are not flexible enough to account for multiple non-linear predictors or for underlying random effects (Wood, 2017). Furthermore, very complex models, such as machine learning approaches (Valbi et al., 2019) and hydrodynamic/oceanographic models (Fauchot et al., 2008; Aleynik et al., 2016; Gillibrand et al., 2016) can require time-intensive and costly computer calculations that are less suitable for daily to weekly predictions on a large spatio-temporal scale. In the present study, we investigate a new approach to HAB modeling and prediction that combines flexibility, simplicity and computing cost-effectiveness. In order to understand how climate change could affect HABs along Canada’s East Coast, we developed a generalized additive mixed model (GAMM) that can predict A. catenella occurrence using an ecological (realized) niche approach (Yee and Mitchell, 1991; Guisan et al., 2002). We ran the GAMM with downscaled simulations from two Earth system models (CanESM2 and MPI-ESM-LR) under two climate scenarios (RCP 4.5 and 8.5) projected throughout the 21st century.

Materials and Methods

A. catenella Database and Environmental Conditions

Data of A. catenella abundances for biological model fitting were obtained from the DFO/IML harmful algae monitoring program (DFO-HAB) spanning the period, 1994–2017, at 12 recurrent coastal monitoring stations located across the EGSL on the Canadian East Coast (Figure 1). The stations were regularly sampled (usually on a weekly basis – notably for the 1999–2008 period) from May to October. Sampling included A. catenella abundance as well as sea surface temperature (SST), sea surface salinity (SSS), sea surface concentrations of nitrite (NO₂^–), nitrate (NO₃^–), silicate (SiO^4–), phosphate (PO₃^2–), and Secchi disk depth. A total of 5373 samples were retrieved throughout the 23 y monitoring period. Environmental and phytoplankton data were collected as described in Starr et al. (2017).

To complete this dataset, we also extracted continuous atmospheric measurements from Environment and Climate Change Canada’s historical climate data archives¹. Extracted variables were precipitation summed over 1, 3, and 7 days (PΣ₁,₃,₇), wind speed averaged over 1, 3, and 7 days ($\overline{S}$₁,₃,₇) and wind direction averaged over 1, 3 and 7 days ($\overline{D}$₁,₃,₇). We chose meteorological stations that best represented conditions at the DFO-HAB monitoring sites. The criteria for their selection included station height above sea level (<100 m), distance from sampling site (<50 km) and general near-coast location. When data had not been collected at the same meteorological station for all 23 y covered by the HAB monitoring program, then to complete the data set, we used data from new meteorological stations located near the previous stations.

GAMM Fitting for A. catenella Bloom Occurrences

All GAMMs and associated statistics were computed using R (R Core Team, 2018), with a 95% confidence interval (α = 0.05). The algae count data were zero-inflated (>80% of our count data were zeros) and abundance modeling with a Poisson distribution was deemed inappropriate (Zuur, 2012). Instead, we fitted an occurrence model using a binomial distribution to illustrate A. catenella presence (1) or absence (0). In the present paper, samples containing ≥ 1000 A. catenella cells L^–1 were considered as presence days at bloom level, and were expressed in “risk-days.” We selected this occurrence threshold by taking into account A. catenella associated risks, such as the toxin bioaccumulation in mussels at the level at which they become unsafe for human consumption and cause the closure of shellfish harvesting activities (Blasco et al., 2003). We performed our modeling using a logit link function implemented in package mgcv (Wood, 2017). Smoothing functions for the predictors were based on thin plate regression splines (s) as described by Wood (2003).

Week, month, and year of sampling, SST, SSS, NO₂^–, NO₃^–, SiO^4–, PO₃^2–, Secchi disk depth, PΣ₁,₃,₇, $\overline{S}$₁,₃,₇, and $\overline{D}$₁,₃,₇ were used as initial predictor variables for the model. These variables were chosen by considering their potential effects on A. catenella excystment, growth and bloom occurrence (Parkhill and Cembella, 1999; Weise et al., 2002; Fauchot et al., 2005; McGillicuddy et al., 2011; Starr et al., 2017). To assess model collinearity, we used Pearson’s correlations.

Multiple variable combinations were considered through the process of model selection. Only combinations of the statistically significant variables were tested for GAMM fitting. In our full model, time markers (week, month, and year), SSS, SST, NO₃^–, and PΣ₃ were significant. PΣ₃ and NO₃^– were excluded from our model since they were less significant compared to SSS. We used combinations of remaining variables to fit an ensemble of GAMMs that we evaluated through cross-validation as described in Albouy-Boyer et al. (2016). Random selection was used to generate a training set comprised of 70% of our full data set. The remainder (30%) was used for validation. The Akaike Information Criterion (AIC) was used as a first marker of model quality. We chose the final model in regards to its predictive power using the True Skill Statistic (TSS; Allouche et al., 2006) and the Area Under the Curve statistic (AUC; Fielding and Bell, 1997). Our best model for predictions of A. catenella presence at ≥1000 cells L^–1 was based on SST and SSS. This selected model was appropriate for climatic predictions (decadal scale) on a large spatial scale and was named the HARM (Harmful Algae Risk Model- Equation 1).

To take into account spatial patterns potentially present in the data, sampled stations (Station) were added to the model as a random effect (bs = “re”). Temporal autocorrelation was calculated by means of full and partial autocorrelation function plots following model fitting. Since a lag of approximately 20–30 days was detected in our model, a second order autoregressive moving average process nested within each year sampled was included. We transformed predictions to presence/absence data with a threshold determined by an approach that maximizes (sensibility + sensitivity)/2 (Freeman and Moisen, 2008).

A. catenella-Ocean Model and Validation

The HARM was used in conjunction with a Regional Ocean Model, named CANOPA, developed by Brickman and Drozdowski (2012) for the EGSL, the Scotian Shelf, and the BoF/GoM to model and study the occurrence of A. catenella blooms under current and future atmospheric and oceanographic forcing. This regional ice-ocean model is based on the code from the Nucleus for European Modeling of the Ocean (NEMO) and includes ice cover, tides, oceanic surface momentum, heat and salt fluxes. The grid features a horizontal resolution of 1/12° (∼6.5 km) and a vertical resolution of 46 vertical layers with vertical thickness resolutions varying from 6 m at the surface to 250 m at 5000 m depths. The model includes runoff from 78 main rivers to account for the freshwater fluxes into the domain (main rivers and their average discharge are shown in Figure 1). A simulation for the 1999–2008 periods using the updated atmospheric conditions obtained from the National Centers for Environmental Predictions (NCEP) was successfully compared to field data including temperature and salinity in the EGSL, Scotian Shelf, and BoF/GoM (Kalnay et al., 1996; Brickman and Drozdowski, 2012; Lavoie et al., 2017). This hindcast of daily temperature and salinity was first used to test predictive capabilities of the HARM to predict A. catenella risk-days from May through October 1999–2008 in the EGSL. We compared hindcasted predictions for A. catenella to field observations coming from the DFO-HAB monitoring program to look for spatial and temporal coherencies. Since monitoring stations were not sampled every day, field results needed to be adjusted for visual comparisons with our daily predictions. To achieve this, sums of bloom occurrences for every month and every year, sampled at the monitoring stations, were multiplied by the 1999–2008 average number of days between samples (namely by 7.8 days).

Future Predictions of A. catenella Bloom Occurrences

Future A. catenella bloom occurrences were obtained by the HARM using daily SST and SSS predictions obtained from CANOPA under future atmospheric, hydrologic and boundary oceanic conditions. The future atmospheric conditions were obtained from three downscaled climate scenarios adapted to our study region using the Canadian Regional Climate Model (CRCM), following Long et al. (2016). Three oceanic simulations with high frequency output were performed separately with the CRCM downscaled forcing, as driven by the coarser resolution global simulations, by the Earth System Model (ESM) MPI-ESM-LR, under the representative concentration pathway (RCP) scenarios 8.5 and 4.5 and, as a comparison, the Canadian Earth System Model (CanESM2) under RCP 8.5 (Chylek et al., 2011; Giorgetta et al., 2013; IPCC, 2013). MPI-ESM-LR was first selected as a model because of its global performance in reproducing the historical conditions in the Northwest Atlantic region (Lavoie et al., 2013); one of the most dynamically complex regions of the globe with cold and low salinity waters flowing southward from the Labrador Sea and warm and saline waters flowing northward with the Gulf Stream. We are confident that these key oceanographic processes, which can influence our bordering model conditions (and per se our predictions), are well represented by the MPI-ESM-LR RCP 4.5 and RCP 8.5 scenarios. The second model, CanESM2, that is structurally different to MPI-ESM-LR, was arbitrarily selected among the other CMIP5 multi-model ensemble prior to have in hand a complete high-frequency downscaling ensemble for the region. In terms of future atmospheric CO₂ emissions, RCP 4.5 is a scenario in which mitigation occurs and greenhouse gas emissions are slowed by 2050. RCP 8.5 is considered a “business as usual” scenario in which anthropogenic CO₂ emissions continue to increase until 2100. It is important to note here that predictions of the future evolution of regional oceanic conditions include changes in runoff of the 78 rivers included in the CANOPA model, such as estimated from Lambert et al. (2013)’s hydrological model, which uses precipitation and temperature outputs from the CRCM downscaling of the ESMs’ atmospheric conditions.

From the HARM predictions, we extracted averages of the number of risk-days for each model cell from May through October for decades 1999–2008, 2050–2059, and 2090–2099. Monthly numbers of risk-days were also computed for these decades. The differences in risk-days (ΔDays) were also computed between years 2050–2059 and 1999–2008 (Δ2050–2059) as well as between 2090–2099 and 1999–2008 (Δ2090–2099) in order to highlight the long-term trends and mitigate the biases of the model (overestimations/underestimations). In addition, we present the full time series (2001–2099) of the predicted average number of risk-days from May through October using the MPI-ESM-LR RCP 8.5 predictions. These results are presented for eight “regions” in order to look into spatio-temporal variations of A. catenella risk-days throughout our study area (Figure 1). For each region, we computed simple linear models to investigate how the average number of risk-days varied as a function of time during the next decades. Numbers of risk-days y^–1 averaged over each region were obtained using Equation 2, where the grid cell area is in km².

Additionally, we computed a simple linear model in each grid cell to determine how the number of risk-days varied as a function of time throughout our domain.

Year of emergence of the climate signal (Hawkins and Sutton, 2011) was also estimated for A. catenella bloom occurrence patterns and defined as the year when the long-term trend line is above the mean of the contemporary period + 1 SD. Here, we used the 2001–2020 period as a contemporary reference. In addition, we computed the climate signal emergence time defined as the difference between the year of emergence of the climate signal and the year corresponding to the middle of the contemporary period (Hawkins and Sutton, 2011).

Results

The HARM has two significant predictor variables, SST (p < 2^∗10^–16) and SSS (p < 2^∗10^–16), that amount to a TSS of 0.62 and an AUC of 0.88, which are considered as “useful” and “good” modeling capabilities respectively, according to Zhang et al. (2015). Differences between internal and external validation statistics are very small (TSS SD = 0.03), showing that the model can produce accurate predictions on a wide range of data points in the EGSL (Table 1). Partial residuals showed the positive influence of predictors on A. catenella presence when SST was between 9 and 19°C (optimal value: 14°C) and SSS was between 15 and 28 (optimal value: 23) (Table 1 and Figure 2).

A. catenella Ocean Model Validation

The hindcasted predictions of A. catenella obtained with the HARM using daily SST and SSS from the CANOPA ocean model with NCEP past atmospheric forcing fields were compared to field observations. Figures 3, 4 show that the predictions reproduced relatively well the observed distribution patterns of maxima and minima, as well as spatial extent and seasonal variability of A. catenella risk-days. Overall, the model of risk-days over May–October statistically matched A. catenella observations (AUC = 0.68). Temporal predictions of A. catenella risk-days were visually similar to observations made from May through October 1999–2008, and the large spatial extent simulated in July successfully matched observations (Figures 3, 4). Indeed, the hindcast’s number of risk-days throughout the months of June, July, and August has a good fit (AUC = 0.64, 0.60, and 0.69, respectively). Although AUCs were much lower for May and October (Figure 4), the model is skillful enough in reproducing the very low or non-existent observations for these months. Close examination of predictions and observations reveals nevertheless that our model somewhat overestimates the frequency of risk-days at certain stations – notably along the southern shore of the St. Lawrence Estuary and in Chaleur Bay. For example, through the May–October average, the maximum number of risk-days is 101 for hindcasted data whereas the adjusted observed data show a maximum of 37 risk-day y^–1 (Figure 3). Monthly patterns show a similar overestimation, with hindcasts showing maxima of 31 risk-day along the southern shore of the LSLE (Lower St. Lawrence Estuary) compared with a maximum of 16 risk-day A. catenella adjusted observations in July (Figure 4). Overestimations near the Gaspé Peninsula and in Chaleur Bay are notable during August, September and October. In spite of these overestimations, the predictions of the HARM reproduced well the bloom recurrence in the EGSL.

Future Predictions of A. catenella Occurrences Under the MPI-ESM-LR RCP 4.5 Scenario

Predicted data for decade 1999–2008 under MPI-ESM-LR RCP 4.5 are relatively comparable with our hindcast using NCEP past atmospheric forcings, showing a maximum of 97 risk-day y^–1 on the South shore of the St. Lawrence Estuary, along the Gaspé peninsula (Figure 5). Risk-days along the Gaspé peninsula remain high in 2050–2059 and 2090–2099, with a maxima of 111 and 108 risk-day y^–1, respectively. Zones near rivers along the Lower North Shore show up to 20 risk-day y^–1 in 2090–2099, where no risk-days were reported during the 1999–2008 reference period. In 2090–2099, up to 30 risk-day y^–1 are also predicted along the northern coastline of Prince Edward Island (PEI) and risk-days extend to the Scotian Shelf. For Δ2090-2099, an increase of up to 30 risk-day y^–1 is recorded in the BoF, near Penobscot Bay. Maximum increases of 42 and 57 risk-day y^–1 in the EGSL are forecasted for Δ2050–2059 and Δ2090–2099, respectively. Differences in the monthly number of risk-days are highest in June, September and October, with increases of up to 20 risk-day y^–1 for Δ2090–2099 (Figure 6). In July, August, and September, increases in A. catenella risk-days are aggregated in the St. Lawrence Estuary and in the Northwestern Gulf (NW Gulf), whereas increases in May, June, and October are recorded near PEI and Chaleur Bay. Large monthly decreases of 9 risk-day y^–1 for Δ2050–2059 and 20 risk-day y^–1 for Δ2090–2099 are obtained in July, August, and September mostly around PEI and Chaleur Bay.

Future Predictions of A. catenella Occurrences Under the MPI-ESM-LR RCP 8.5 Scenario

Under the MPI-ESM-LR RCP 8.5 scenario, the annual A. catenella occurrences increase by 45–60 risk-day y⁻¹ (Figure 7). For the 2090–2099 periods, a small number of risk-days are visible on the Scotian Shelf and on Newfoundland’s West Coast. For the same decade, predictions indicate higher risk-days annually compared to MPI-ESM-LR RCP 4.5 along the Lower North Shore, with up to 40 risk-day y^–1. Up to 20 risk-day y^–1 are visible near PEI’s northern point and up to 50 risk-day y^–1 are recorded in the NW Gulf for the 2090–2099 periods. The number of risk-days during the months of May, June, September and October reach 20 days all throughout our study region for the 2090–2099 period (Figure 8). An increase in number of risk-days compared to MPI-ESM-LR RCP 4.5 is observed in the BoF, with maximum increases of 45 and 86 risk-day y^–1 for Δ2050–2059 and Δ2090–2099, respectively. Maximum decreases reach 15 risk-day y^–1 for Δ2050–2059 near Anticosti Island and 35 risk-day y^–1 for Δ2090–2099 near Chaleur Bay compared to the reference period. July and August show decreases of 20 risk-day y^–1 along the southern shore of the LSLE in Δ2090–2099. In October, increases are visible mostly near PEI for Δ2050–2059, and concentrated in the LSLE for Δ2090–2099.

Future Predictions of A. catenella Occurrences Under the CanESM2 RCP 8.5 Scenario

Forecasted A. catenella risk-days under CanESM2 RCP 8.5 are most frequent in the LSLE compared to our MPI-ESM-LR RCP 8.5 scenario (Figure 9). In this scenario, future conditions on the Lower North Shore are predicted to become more favorable to A. catenella’s presence, where up to 80 risk-day y^–1 are predicted in 2090–2099. Only rare occurrences are forecasted on the Scotian Shelf and in the BoF/GoM in 2090–2099, with maxima of 10–20 risk-day y^–1. Compared to MPI-ESM-LR RCP 8.5, CanESM2 RCP 8.5 predicts a lower number of risk-days for PEI, with maxima of 10 risk day y^–1. Increases mostly occur in the LSLE for Δ2090–2099, reaching up to 122 risk-day y^–1. In CanESM2 RCP 8.5, maximum decreases are more visible than in MPI-ESM-LR RCP 8.5 and reach up to 34 and 46 risk-day y^–1 for Δ2050–2059 and Δ2090–2099, respectively. As with other models presented herein, CanESM2 RCP 8.5 predicts increases throughout the months of May, June, and October in the St. Lawrence Estuary throughout the 21st century (Figure 10). Similar to MPI-ESM-LR RCP 8.5, decreasing risk-days in CanESM2 RCP 8.5 are most visible in July, August and September for Δ2090–2099 and take place in the Gaspé Current and in Chaleur Bay.

Evolution of A. catenella Occurrences Throughout the 21st Century Under the MPI-ESM-LR RCP 8.5 Scenario

From 2001 to 2099, our model based on the downscaled MPI-ESM-LR RCP 8.5 scenario predicts significant increases of A. catenella risk-days in all regions previously defined in Figure 1, except for Chaleur Bay and the Scotian Shelf (Figure 11). The LSLE, NW Gulf, and PEI show the most dramatic increases in risk-days (SD), with increases of 3.97 (0.35), 3.53 (0.53), and 1.11 (0.29) risk-day per decade (Table 2). The Magdalen Islands, the Lower North Shore, the Scotian Shelf, and the BoF exhibit increases in risk-days (SD) on a smaller scale, with increases of 0.70 (0.18), 0.20 (0.05), 0.03 (0.02), and 0.01 (0.00) risk-day per decade respectively. The LSLE has the earliest climate signal emergence year (2031) and the shortest climate signal emergence time (20.85 years). The Lower North Shore is the second fastest to emerge in 2032 with an emergence time of 21.14 years. The NW Gulf, PEI, the Magdalen Islands and the BoF regions emerge in 2038, 2043, 2057, and 2038, respectively. Their times of emergence vary between 27.78 and 46.81 years. Increases in the interannual variability of risk-days are predicted until 2099 in most areas of the model domain.

Figure 12 confirms that increases in A. catenella risk-days from 2001 to 2099 are most intense in the LSLE and in the NW Gulf, with a maximum increase of 9.44 risk-day per decade. Increases near river mouths are particularly visible, notably along the Lower North Shore as well as in the BoF and the GoM. Small increases are also forecasted throughout the EGSL, around PEI and the Magdalen Islands, on the Scotian Shelf and on Newfoundland’s West Coast (<2 risk-day per decade). The standard error of the slope estimate is highest in Chaleur Bay and in the St. Lawrence Estuary (±1.1 risk-day per decade). Linear trends are significant mostly in the LSLE, the NW Gulf and along the Lower North Shore, as well as around PEI and the Magdalen Islands (p < 0.05). Significant risk-day increases are also predicted in the BoF, near Penobscot Bay, Saint-John River and Casco Bay.

Discussion

In this study, we developed a model that can predict A. catenella potential occurrences ≥ 1000 cells L^–1 (referred herein as “risk-days”) using an ecological niche approach. The model was run using two ESMs (MPI-ESM-LR andCanESM2) under two climate scenarios (RCP 4.5 and 8.5). RCP 4.5 is a scenario in which greenhouse gas emissions increase steadily until 2040, before slowing down by 2050. In this mitigation scenario, global temperatures are expected to rise by 1.1–2.6°C by 2081–2100 compared to 1986–2005. In RCP 8.5, anthropogenic CO₂ emissions increase rapidly until 2100. In this scenario termed “business as usual,” global temperatures are expected to rise by 2.6–4.8°C by 2081–2100 relative to 1986–2005. Changes in global temperatures and CO₂ levels are expected to have major effects on the hydrological cycle as a whole, including water temperature and salinity. Projection results were used to evaluate future A. catenella risk days throughout the 21st century (IPCC, 2013). All climate change scenarios explored predict increases in frequency, spatial expansion, and changes in seasonality of A. catenella blooms with expected changes in water temperature and salinity along Canada’s East Coast.

A. catenella GAMM Description

The HARM succeeded at modeling average A. catenella spatio-temporal extent over a 23 years monitoring program that took place in the EGSL. The model fit our samples accurately (AUC = 0.88). Our GAMM predicts A. catenella bloom occurrences to be most likely when SSS is 15–28 and SST is 9–19°C (optima at 23 and 14°C, respectively), which is consistent with several laboratory and field studies showing that optimal growth rates and/or high abundances of this species are measured within this environmental window (Etheridge and Roesler, 2005; Condie et al., 2019; Paredes-Mella et al., 2020). Our results compare to past observations made by Weise et al. (2002) in the EGSL near Sept-Îles, who showed that the highest A. catenella concentrations occurred when SSS was 20–26 and SST was ≥12°C. In July 1998, during a red tide event in the LSLE, A. catenella in situ growth rates were highest when SSS was ≤24.5 and SST was ≥10°C (Fauchot et al., 2005). Etheridge and Roesler (2005) have shown experimentally that optimal growth rates were about 15°C and 25 salinity units for two A. catenella strains (previously known as A. fundyense) isolated from the GoM and BoF.

Temperatures and salinity are already known to have direct effects on Alexandrium’s growth (Bill et al., 2016; Brosnahan et al., 2020). Salinity, however, could fulfill multiple roles in our model as it is linked to a variety of different factors which are known to affect A. catenella bloom occurrence. Low salinity values are a marker for increased freshwater input in the EGSL, either from river runoff and/or from precipitation, and have been identified as a trigger for bloom occurrence (Weise et al., 2002; Starr et al., 2017). A fresher surface layer increases stratification strength, which in turn favors active vertically migrating dinoflagellates such as A. catenella (in contrast to less mobile phytoplankton) to accumulate in the euphotic zone and grow with an access to nutrients at depth if depleted in surface (Margalef, 1978; Weise et al., 2002; Hallegraeff, 2010; Glibert, 2016; Condie et al., 2019). In the EGSL’s coastal areas, salinity is often identified as the main driver of vertical stratification (Koutitonsky and Bugden, 1991). This typical association between A. catenella and freshwater plumes has also been previously attributed to the beneficial effects of riverine input of nutrients that can serve as growth stimulants (Fauchot et al., 2005; Collos et al., 2007; McGillicuddy et al., 2011). However, model fitting during our study suggested that SSS as a predictor variable is likely mostly linked to physical factors, rather than nutrients (notably PO₃^2– and NO₃^–) such as suggested by Condie et al. (2019). In fact, SSS was a more significant predictor variable when combined to SST, compared to NO₃^– or PO₃^2– combined with SST (data not shown). We do not, however, exclude that freshwater runoff often causes humic substance injections into the EGSL, which have been shown to favor A. catenella growth (Gagnon et al., 2005).

No other environmental variables examined improved the model significantly, although these variables had initially been chosen for their potential influence on A. catenella. For example, high wind speeds are known to inhibit A. catenella occurrences (Weise et al., 2002; Fauchot et al., 2005; Starr et al., 2017). However, the addition of wind intensity as a predictor was found to be insignificant when GAMM fitting was performed. In addition, the close examination of partial response curves indicated that the effect of wind intensity on A. catenella occurrence was weak. Wind speed is known to generate vertical mixing of the water column that can eventually erode stratification. As this latter is directly linked to SSS, this may explain weaker power of prediction of winds. Together, these results highlight that SST and SSS are the most important environmental factors controlling A. catenella bloom occurrence in the EGSL, as well as in other regions of the globe (Fauchot et al., 2005; Li et al., 2009; McGillicuddy et al., 2011; Bill et al., 2016; Condie et al., 2019). Globally, these environmental drivers are well recognized as important for the structuring of ocean plankton communities (Glibert, 2016).

Observations and Hindcasts

When driven by the 1999–2008 hindcast, the HARM reproduced A. catenella occurrences fairly well (AUC = 0.68) and succeeded in hindcasting past general patterns of A. catenella spatial extent and bloom seasonality in the EGSL. Occurrence maxima happened in July and August and occurrence minima took place in May, September and October. Maximum A. catenella concentrations were found in the LSLE, on the northern shore of the St. Lawrence Estuary, and around the Gaspé Peninsula. Blasco et al. (2003) observed similar spatio-temporal patterns during the 1984–1994 periods, with A. catenella maxima in the LSLE and around the Gaspé Coast, through the months of June, July, and August. Fauchot et al. (2008) modeled similar seasonal and spatial patterns when looking at A. catenella abundance in the LSLE. The model correctly reproduced the bloom initiation in June and abundance maxima in the northern LSLE.

Occurrence overestimations are present on the southern shore of the St. Lawrence Estuary along the Gaspé Current, and in Chaleur Bay. In the Gaspé Current, our model predicts a high number of A. catenella risk days. In real conditions, however, A. catenella cells along the Gaspé Current are exposed to strong currents and are likely to be advected away from this area (El-Sabh and Silverberg, 1990). Our model does not account for cell advection, which explains why higher risk days are predicted in this section of the EGSL. A similar mechanism could also explain the differences observed in Chaleur Bay. Although relatively simple, our GAMM is deemed appropriate for A. catenella occurrence prediction in our study region.

The BoF and GoM usually show high A. catenella occurrences during the bloom season (Anderson et al., 2014; Martin et al., 2014; Bates et al., 2020; McKenzie et al., 2020), which is not completely captured by our model. The BoF/GoM is an oceanographically complex region where A. catenella cells are advected from a “donor” to a “receiver” region (Keafer et al., 2005). Note that the analysis of Lavoie et al. (2017) also shows that the oceanographic predictions obtained in the GoM have a greater uncertainty than those obtained in the EGSL. Lavoie et al. (2017) attribute these biases to weaknesses in the predictions obtained from the ESMs along the open boundaries. In spite of these biases, the predicted trends remain of interest.

Our model does not account for A. catenella acclimation and adaptation to future climate conditions. However, HAB range shifts induced by warming temperatures are expected and have already been observed, as well as lengthening of the bloom season window in temperate regions (Wells et al., 2015; Barton et al., 2016; Gobler et al., 2017; Brandenburg et al., 2019). Considering these findings, we determined that adding acclimation and adaptation to our model by shifting or widening A. catenella’s niche was not necessary. Instead, we have kept our model flexible by avoiding temporal variables in our GAMM (such as weekly or monthly sampling) to allow for future A. catenella occurrences to happen outside of their current bloom window.

Projected Changes in Occurrence of the Toxic Algae Alexandrium catenella Along Canada’s East Coast

Our results show that A. catenella risk days will increase and show more year to year variability throughout the 21st century along Canada’s East Coast. Future climate conditions will favor a warmer and fresher surface layer (Supplementary Figures S1–S18), that will expand across our study region. The effect of increased runoff and expansion of freshwater plumes on A. catenella risk days are evident along the Lower North Shore, in the LSLE and in the NW Gulf under all climate scenarios investigated. Under MPI-ESM-LR RCP 4.5, A. catenella risk days increase at the lowest rate, with maximum increases > 50 risk-day y^–1 recorded for 2090–2099. Changes are most intense under CanESM2 and MPI-ESM-LR RCP 8.5, the “business as usual” scenarios, and reach > 80 risk-day y^–1 increases in 2090–2099. In these scenarios, rates of increase are projected to be more pronounced in the LSLE and in the NW Gulf. Under MPI-ESM-LR RCP 8.5, the climate signal emergence years vary between 2031 and 2057. Regions with the earliest emergence years, such as the LSLE, the NW Gulf and PEI, could be the first to see associated effects on socio-economically valuable marine fauna.

Increasing bloom frequency will be accompanied by spatial expansion of A. catenella risk days across the EGSL and on the Scotian Shelf. Spatial expansion and range shifts of phytoplankton species have already been documented and predicted (Hallegraeff, 2010 and references therein), especially poleward shifts, related to higher sea temperatures in the North Atlantic and Arctic Oceans (Neukermans et al., 2017; Renaut et al., 2018). Our downscaled predictions clearly show increased SST in the future as well as decreased SSS over most of the Northwest Atlantic coast (Supplementary Figures S1–S18). Consequently, stratification exhibits an increase in future predictions (Long et al., 2016). This process will contribute to increasing A. catenella’s spatial extent, by limiting mixing with deeper waters. Moreover, the estuarine circulation shows seasonal changes with a decrease in the estuarine ratio in winter and an increase in summer (Long et al., 2016). Since most A. catenella occurrences happen during summer, the amplification of the estuarine circulation during that period will likely contribute to a spatial expansion of HABs. MPI-ESM-LR RCP 4.5 and 8.5 both show increases in risk days near PEI, the Magdalen Islands and across the Scotian Shelf. Lower SSS in these regions (Supplementary Figures S1–S9) enables A. catenella to expand into the Gulf. Conversely, in CanESM2 RCP 8.5, the spatial expansion is not as intense as in MPI-ESM-LR scenarios; the expansion is limited by high temperatures (≥19°C) that reach PEI, the Magdalen Islands and the Scotian Shelf. In this instance, our climate “worst case scenario” might dampen the effects of A. catenella occurrences on PEI, the Scotian Shelf and the Magdalen Islands. However, mitigation efforts will most likely not prevent A. catenella from widening its reach across the EGSL.

Changes in bloom seasonality of A. catenella were observed across all climate change scenarios investigated. A. catenella is projected to bloom earlier in the season (May) and to remain active until September-October due to longer favorable conditions in the EGSL. Earlier blooming in phytoplankton has already been measured and verified in the North Atlantic (Barton et al., 2016; Gobler et al., 2017) and is predicted to continue during the 21st century (Brosnahan et al., 2020). In the mitigation scenario MPI-ESM-LR RCP 4.5, more risk days are predicted to take place through the entire May–October period. In MPI-ESM-LR RCP 8.5 and CanESM2 RCP 8.5 estimates, a decrease of risk days is nevertheless predicted mostly during the months of July and August in some regions of the EGSL. This period corresponds to the usual peak of A. catenella blooms observed for our reference period 1999–2008, previously reported by Blasco et al. (2003), Fauchot et al. (2005), and Starr et al. (2017) in the EGSL. In both RCP 8.5 scenarios, increasingly high temperatures during summer are projected in 2050–2059 and 2090–2099 (Supplementary Figures S10–S18). In some regions of the EGSL, these conditions exceed A. catenella’s thermal niche and consequently limit the number of risk days during this period of the year. Decreases are most noticeable along the Gaspé Peninsula, in Chaleur Bay and around PEI. Similar declines of Alexandrium midseason growth rates were reported for present-day conditions in the North Atlantic by Gobler et al. (2017), when SST values became higher than the species’ optimal yield. As for spatial expansion, changes in bloom seasonality could affect mussel farms as well as other open ocean aquaculture sites in these areas.

In all climate scenarios studied, increases in risk days are predicted at the mouth of the St. Croix River, which is a known cyst bed and nursery for A. catenella cells in the BoF. Moreover, increases are also predicted in Penobscot Bay and Casco Bay, regions of the GoM to which A. catenella cells are usually advected, from the St. Croix River (Keafer et al., 2005). These observations indicate that the BoF and the GoM could become more favorable to A. catenella by the end of the 21st century, as it is the case in the EGSL.

Impacts of the Projected Changes

Alexandrium catenella blooms represent an immediate threat to marine vertebrates. Toxins produced by A. catenella are efficiently transferred to higher trophic levels via bioaccumulation in invertebrate filter-feeders, who remain relatively tolerant to PSTs (Doucette et al., 2006a). In 2008, a mass mortality event of marine species was recorded during an intense A. catenella bloom in the LSLE (Starr et al., 2017). During this event, mortalities of beluga whales (Delphinapterus leucas), seals, porpoises, marine birds and fish were recorded. At-risk species vulnerable to PSTs, such as beluga whales, are important tourist attractions and hold historical value in the EGSL (Williams et al., 2017). Projected increases in HAB occurrence may further endanger these species and the socio-economic context surrounding them. In addition, the expansion of A. catenella risk days could pose a threat to the mussel industry on PEI’s coastlines, which is already exposed to domoic acid produced by Pseudo-nitzschia spp. (Bates et al., 1989). In 2016, PEI produced > 80% of Canada’s mussel production, which is valued at CAN$23.8 million (Bates et al., 2020). A. catenella occurrence in this region of the EGSL could trigger more frequent mussel farm closures. In addition, the endangered North Atlantic right whale is also potentially threatened by A. catenella’s expansion (Doucette et al., 2006b, 2012).

Perspectives

As emphasized by Wells et al. (2015) and by the IPCC’s Special Report on the Ocean and Cryosphere in a Changing Climate (2019), our results further support the importance of maintaining sustained HAB monitoring programs like the one presented herein. The GAMM developed in this study is eventually intended to be used as an operational prediction tool for governments and the industry. Predictions produced by our model could be verified a posteriori using short term, high resolution hindcasts, as well as various occurrence thresholds (500 cells L^–1, 5000 cells L^–1). Our model could also be refitted annually to account for any changes in A. catenella’s realized niche and to increase model accuracy. Such short-term predictions have been effectively applied in Europe through the ASIMUTH program (Maguire et al., 2016). Short term predictions would be beneficial to mussel farms and aquaculture sites as an “early warning” tool. This tool would enable industries to make prompt management decisions when A. catenella-favorable conditions are predicted in their region.

Conclusion

Long-term A. catenella predictions produced by the HARM based on three climate change scenarios show increasing frequency, spatial extent and lengthening of the bloom season along Canada’s East Coast by the end of the 21st century. The mitigation scenario, which follows RCP 4.5, still shows increasing A. catenella occurrences across our study region, albeit on a smaller scale than under both RCP 8.5 scenarios. An aggregation of A. catenella occurrences in the LSLE and spatial expansion of blooms into the Gulf and on the Scotian shelf could threaten endangered species, such as the beluga whale and North Atlantic right whale. Changes in bloom timing and increased interannual variability could become problematic to mussel farms and aquaculture sites in the EGSL. The model presented herein is proposed to become an early warning tool for governments and industries on the Canadian East Coast. Our modeling approach could eventually be applied to appropriately downscaled climate models in other coastal regions of the globe where A. catenella poses a threat, particularly in mid- and high latitudes where warming and freshening of waters are intensifying.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MS and JC designed the project. All the authors helped to design the experiment and edit the manuscript. AB-R analyzed the data and wrote the initial manuscript.

Funding

This work was supported by the Aquatic Climate Change Adaptation Services Program (ACCASP) within the Department of Fisheries and Oceans Canada.

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.

Acknowledgments

We thank all DFO laboratory technicians that analyzed HAB monitoring samples, and more specifically Jean-Yves Couture and Sylvie Lessard. We also thank Alain Caron (UQAR) and Caroline Lehoux (DFO) for the help with GAMM fitting and cross-validation. Finally, we are also very grateful to Nicolas Lambert (DFO) who prepared the atmospheric files for running the CANOPA model.

Supplementary Material

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

Supplementary Figure 1 | Average sea surface salinity from May through October along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) 1999–2008, (B) 2050–2059, and (C) 2090–2099 and difference in average sea surface salinity between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 2 | Average sea surface salinity along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 3 | Average sea surface salinity along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) May 1999–2008, (B) difference in average sea surface salinity between May 1999–2008 and 2050–2059 and (C) difference in average sea surface salinity between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 4 | Average sea surface salinity from May through October along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) 1999–2008, (B) 2050–2059, and (C) 2090–2099 and difference in average sea surface salinity between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 5 | Average sea surface salinity along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 6 | Average sea surface salinity along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) May 1999–2008, (B) difference in average sea surface salinity between May 1999–2008 and 2050–2059 and (C) difference in average sea surface salinity between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 7 | Average sea surface salinity from May through October along Canada’s East Coast under CanESM2 RCP 8.5 for (A) 1999–2008, (B) 2050–2059 and (C) 2090–2099 and difference in average sea surface salinity between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 8 | Average sea surface salinity along Canada’s East Coast under CanESM2 RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 9 | Average sea surface salinity along Canada’s East Coast under CanESM2 RCP 8.5 for (A) May 1999–2008, (B) difference in average sea surface salinity between May 1999–2008 and 2050–2059 and (C) difference in average sea surface salinity between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 10 | Average sea surface temperature (°C) from May through October along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) 1999–2008, (B) 2050–2059, and (C) 2090–2099 and difference in average sea surface temperature between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 11 | Average sea surface temperature (°C) along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 12 | Average sea surface temperature (°C) along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) May 1999–2008, (B) difference in average sea surface temperature between May 1999–2008 and 2050–2059 and (C) difference in average sea surface temperature between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 13 | Average sea surface temperature (°C) from May through October along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) 1999–2008, (B) 2050–2059, and (C) 2090–2099 and difference in average sea surface temperature between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 14 | Average sea surface temperature (°C) along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 15 | Average sea surface temperature (°C) along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) May 1999–2008, (B) difference in average sea surface temperature between May 1999–2008 and 2050–2059 and (C) difference in average sea surface temperature between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 16 | Average sea surface temperature (°C) from May through October along Canada’s East Coast under CanESM2 RCP 8.5 for (A) 1999–2008, (B) 2050–2059, and (C) 2090–2099 and difference in average sea surface temperature between (D) 1999–2008 and 2050–2059 and (E) 1999–2008 and 2090–2099.

Supplementary Figure 17 | Average sea surface temperature (°C) along Canada’s East Coast under CanESM2 RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 18 | Average sea surface temperature (°C) along Canada’s East Coast under CanESM2 RCP 8.5 for (A) May 1999–2008, (B) difference in average sea surface temperature between May 1999–2008 and 2050–2059 and (C) difference in average sea surface temperature between May 1999–2008 and 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October.

Supplementary Figure 19 | Average number of risk days with Alexandrium catenella ≥ 1000 cells L^–1 along Canada’s East Coast under MPI-ESM-LR RCP 4.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October estimated using the HARM.

Supplementary Figure 20 | Average number of risk days with Alexandrium catenella ≥ 1000 cells L^–1 along Canada’s East Coast under MPI-ESM-LR RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October using the HARM.

Supplementary Figure 21 | Average number of risk days with Alexandrium catenella ≥ 1000 cells L^–1 along Canada’s East Coast under CanESM2 RCP 8.5 for (A) May 1999–2008, (B) May 2050–2059, and (C) May 2090–2099. Same pattern applies for (D–F) June, (G–I) July, (J–L) August, (M–O) September, and (P–R) October using the HARM.

Footnotes

1. ^ https://climat.meteo.gc.ca/

References