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

Alexandrium catenella produces paralytic shellfish toxins that affect marine fisheries and aquaculture as well as ecosystem and human health worldwide. This harmful algal species is extremely sensitive to environmental conditions and potentially to future climate change. Using a generalized additive mixed model (GAMM) we studied the potential effects of changing salinity and temperatures on A. catenella bloom (≥1000 cells L–1) occurrence along Canada’s East Coast throughout the 21st century. Our GAMM was applied to two high greenhouse gas emissions scenarios (RCP 8.5) and one mitigation scenario (RCP 4.5). Under present-day conditions, our model successfully predicted A. catenella’s spatio-temporal distribution in Eastern Canada. Under future conditions, all scenarios predict increases in bloom frequency and spatial extent as well as changes in bloom seasonality. Under one RCP 8.5 scenario, A. catenella bloom occurrences increased at up to 3.5 days per decade throughout the 21st century, with amplified year-to-year variability. Blooms expanded into the Gulf of St. Lawrence and onto the Scotian Shelf. These conditions could trigger unprecedented bloom events in the future throughout our study region. In all climate scenarios, the bloom season intensified earlier (May–June) and ended later (October). In some areas of the Gulf of St. Lawrence, the thermal habitat of A. catenella was exceeded, thereby locally reducing bloom risk during the summer months. We conclude that an increase in A. catenella’s environmental bloom window could further threaten marine fauna including endangered species as well as fisheries and aquaculture industries on Canada’s East Coast. Similar impacts could be felt in other coastal regions of the globe where warming and freshening of waters are intensifying.


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 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 . 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 Collos et al., 2007;McGillicuddy et al., 2011), pH and CO 2 (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 2 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 timeintensive 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.

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 2 − ), nitrate (NO 3 − ), silicate (SiO 4− ), phosphate (PO 3 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 1 . Extracted FIGURE 1 | Map of the Atlantic coast of Canada. Names of harmful algae sampling stations are in bold and their locations are indicated by black dots. Extent of the regional boxes are defined by colored polygons. Hundred meter bathymetry is indicated by gray contour lines. Main rivers (with average discharge > 100 m 3 s −1 ) and their average discharge between 1999 and 2008 are indicated by diamonds. variables were precipitation summed over 1, 3, and 7 days (P 1 , 3 , 7 ), wind speed averaged over 1, 3, and 7 days (S 1 , 3 , 7 ) and wind direction averaged over 1, 3 and 7 days (D 1 , 3 , 7 ). 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).
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 3 − , and P 3 were significant. P 3 and NO 3 − 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 highfrequency downscaling ensemble for the region. In terms of future atmospheric CO 2 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 2 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 2 .

Days =
(Grid cell area × number of risk days )

Region area
( 2) 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 longterm 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.  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 riskday 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.      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 "riskdays") 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 2 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 2 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 spatiotemporal 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 . 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 Collos et al., 2007;  2 | Changes in mean annual risk days predicted over the period 2001-2099 under MPI-ESM-LR RCP 8.5 within each region defined in Figure 1 (LSLE, Lower St. Lawrence Estuary): linear model slope ± Standard deviation (SD) and p-value, average sum of risk days ± SD of the contemporary period (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020), year of emergence of climate signal and climate signal emergence time. Year of emergence is defined by the year when the linear model slope has a higher or equal value to the mean of the contemporary period + 1 SD. † † † Emergence time is defined as the difference between the year of emergence and the year corresponding to the middle of the contemporary period (2010.5).
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 3 2− and NO 3 − ) such as suggested by Condie et al. (2019). In fact, SSS was a more significant predictor variable when combined to SST, compared to NO 3 − or PO 3 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 .
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 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 presentday 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 . 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(Doucette et al., , 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. catenellafavorable 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.

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.