Investigating Pseudo-nitzschia australis introduction to the Gulf of Maine with observations and models.

In 2016, an unprecedented Pseudo-nitzschia australis bloom in the Gulf of Maine led to the first shellfishery closures due to domoic acid in the region's history. In this paper, potential introduction routes of P. australis are explored through observations, a hydrodynamic model, and a Lagrangian particle tracking model. Based on particle tracking experiments, the most likely source of P. australis to the Gulf of Maine was the Scotian Shelf. However, in 2016, connectivity between the Scotian Shelf and the bloom region was not significantly different from the other years between 2012 and 2019, nor were temperature conditions more favorable for P. australis growth. Observations indicated changes on the Scotian Shelf in 2016 preceded the introduction of P. australis: increased bottom salinity and decreased surface salinity. The increased bottom salinity on the shelf may be linked to anomalously saline water observed near the coast of Maine in 2016 via transport through Northeast Channel. The changes in upstream water mass properties may be related to the introduction of P. australis, and could be the result of either increased influence of the Labrador Current or increased outflow from the Gulf of St. Lawrence. The ultimate source of P. australis remains unknown, although the species has previously been observed in the eastern North Atlantic, and connectivity across the ocean is possible via a subpolar route. Continued and increased monitoring is warranted to track interannual Pseudo-nitzschia persistence in the Gulf of Maine, and sampling on the Scotian Shelf should be conducted to map upstream P. australis populations.

Introduction

The harmful algal genus Pseudo-nitzschia is a lightly silicified diatom of growing global presence and increasing concern (Bates et al., 2018; Trainer et al., 2012). Species of the genus are often described as “cosmopolitan” because they can persist in a wide range of temperature conditions (Hasle, 2002), and because their distributions span estuaries, coastal environments, and the open ocean (Bates et al., 2018). Some Pseudo-nitzschia species produce domoic acid (DA), a neurotoxin responsible for Amnesic Shellfish Poisoning. More than 52 Pseudo-nitzschia species have been identified worldwide, and at least half are confirmed DA producers (Bates et al., 2018).

In 2016 in the Gulf of Maine (GOM), DA concentrations exceeded the regulatory limit of 20 μg DA g−1 of shellfish tissue, leading to the first regional DA-induced shellfishery closures (Bates et al., 2018; Clark et al., 2019; Hubbard et al., 2017; Lewis et al., 2017). The bloom began in the Bay of Fundy in late September 2016, and progressed along the coast of Maine, continuing into the second week of October. Fourteen Pseudo-nitzschia species had been identified in the region prior to the DA event (Fernandes et al., 2014), but the record DA concentrations were caused primarily by the novel appearance of P. australis (Bates et al., 2018; Clark et al., 2019).

In observations from the 2016 DA event, neither cell concentrations nor relative species abundance correlated strongly with various environmental parameters (temperature, salinity, nitrate, ammonium, silicic acid, phosphate, or nutrient ratios). In addition, except for salinity, environmental parameters were not significantly different in 2016 compared to previous years (2012–2015) (Clark et al., 2019). Salinity was significantly higher in 2016 compared to the four previous years (0.4–1.5 PSU greater), but the salinity difference was not enough to result in improved growth conditions based on the literature for P. australis (Doucette et al., 2008; Thessen et al., 2005). Clark et al. argued that, because of the similarities in environmental parameters across the study period, it was unlikely that P. australis was able to bloom because of a sudden favorable change in conditions. Moreover, P. australis was not detected in the region prior to 2016 in either microscopy or genetic analysis.

An alternative explanation to changing environmental conditions for the novel bloom in 2016 is that P. australis was introduced to the region in a process linked to the anomalously saline conditions in the GOM (Clark et al., 2019). Of the primary source waters flowing into the GOM, Gulf Stream Ring water (GSRW) is saltier than Slope Water or Scotian Shelf Water (Townsend et al., 2015). Previous studies have shown that GSRW can contribute to or modify Northeast Channel inflows (Brooks, 1987), and measurements at the Northeast Regional Association of Coastal Ocean Observing Systems (NERACOOS) Buoy M suggested the presence of GSRW in the GOM interior in July 2013 (Townsend et al., 2015) and again in July 2016 (Clark et al., 2019). GSRW is not the only potential source for anomalously saline water, however. From 1990 to 2015, Brickman et al. (2018) observed near-bottom anomalies upstream on the Scotian Shelf that alternate between warm and saline or fresh and cool. With a model they showed that such anomalies form due to interactions between the Gulf Stream and Labrador Current off of the Grand Banks in the North Atlantic Ocean, propagate toward the southwest, and penetrate into the GOM via the Northeast Channel. Based on both observations and model results, Brickman et al. argued that warm/saline anomalies increased in frequency and magnitude from 2006 to 2015.

The correspondence of P. australis and an anomalous water mass is not without precedence. In 2015, a P. australis bloom on the West Coast of the United States led to record-breaking DA concentrations in Monterey Bay, CA (Ryan et al., 2017). The bloom occurred during the North Pacific Ocean Warm Anomaly, during which sea surface temperatures were more than 2.5 °C higher than the long-term average (McCabe et al., 2016). The anomalously warm water is thought to have contributed to the bloom, either as a primary factor (McCabe et al., 2016) or as a contributing factor along with intermittent upwelling of waters with an anomalously low silica-to-nitrate ratio (Ryan et al., 2017). Thus, the anomalies observed in recent years in the GOM region (increased salinity along the coast of Maine and in Jordan Basin (Clark et al., 2019), and warm, saline anomalies on the Scotian Shelf (Brickman et al., 2018)) may relate to the introduction of P. australis.

This paper uses field data, a hydrodynamic model, and a Lagrangian particle tracking model to address the following questions:

What are the most likely sources of P. australis to the Gulf of Maine?

Was the connectivity from potential source regions to the Gulf of Maine different in 2016 compared to previous and subsequent years?

Do model results indicate differences in hydrographic conditions in the Gulf of Maine that may have affected P. australis bloom timing, location, or distribution?

The following sections describe the study site, field data, hydrodynamic model, and particle tracking model and analysis. Simulations of the hydrodynamics and Lagrangian transport are compared for the years 2012–2019, with a focus on 2016. We evaluate hypotheses for P. australis introduction, effects of changing hydrodynamics on the 2016 bloom, influences of Scotian Shelf processes on the GOM in 2016, potential links to other P. australis populations in the North Atlantic, and implications of this research for future studies and monitoring practices.

Methods

Study site – the Gulf of Maine

The GOM lies between Cape Cod, MA, USA at 42°N and the Bay of Fundy, Canada at 44.5°N (Fig. 1A). Sea surface temperatures range from 6 °C in the winter to 22.5 °C in the summer, and salinity ranges from 29 PSU near the coast to 33.5 PSU offshore (Li and He, 2014). It is considered a gulf because its offshore boundaries, Georges Bank and Browns Bank, are shallower than 100 m, while its three major basins – Georges Basin, Jordan Basin, and Wilkinson Basin – are deeper than 200 m (Townsend et al., 2006). Gulf inflows occur via the Northeast Channel and south of Nova Scotia. The general circulation around the GOM is cyclonic, and outflows occur via the Great South Channel and the Northeast Channel (Brooks, 1985; Lynch et al., 1997; Pettigrew et al., 2005; Xue et al., 2000). Water properties in the GOM interior are set by the water mass properties at the inflow locations: Northeast Channel deep inflows are a mix of relatively warm and salty Warm Slope Water and slightly cooler and fresher Labrador Slope Water, while Nova Scotia surface inflows are made up of cool, fresh Scotian Shelf Water (Smith et al., 2012; Townsend et al., 2015).

Fig. 1.

(left) Climatological circulation of the GOM as described in Pettigrew et al. (2005) and adapted in Anderson et al. (2005).1 GOMCP stands for Gulf of Maine Coastal Plume, WMCC stands for Western Maine Coastal Current, and EMCC stands for Eastern Maine Coastal Current. Two-letter abbreviations are for Massachusetts (MA), New Hampshire (NH), and Nova Scotia (NS). (right) Map of the Gulf of Maine showing the ROMS domain (blue box), the Scotian Shelf portion of the ROMS domain (red shaded box), NERACOOS buoys (green circles), Halifax Station 2 (magenta circle), 2016 DFO survey locations (grey diamonds), transects for analysis (red lines), and particle analysis regions (solid blue boxes are sources and destinations, while dashed blue boxes are destinations only). The contour lines are drawn at 25m, 50m, 75m, every 100m to 1000 m, and 2000 m.

1The left panel was published in Deep-Sea Research II, Vol 52, Anderson, Donald M., Keafer, Bruce A., McGillicuddy, Dennis J., Mickelson, Michael J., Keay, Kenneth E., Libby, P. Scott, Manning, James P., Mayo, Charles A., Whittaker, David K., Hickey, J. Michael, He, Ruoying, Lynch, Daniel R., Smith, Keston W. “Initial observations of the 2005 Alexandrium fundyense bloom in southern New England: General patterns and mechanisms.” p. 2858, Copyright Elsevier (2005).

Field data

The Department of Fisheries and Oceans Canada (DFO) has conducted annual summer surveys of the Scotian Shelf since 1970. In 2016, the survey was conducted from June 28 to Aug 15, and a total of 250 stations were sampled with Conductivity-Temperature-Depth (CTD) profiles (Fig. 1b). At the surface, 50 m, 100 m, and the bottom, temperature and salinity values were interpolated onto a 0.2-by-0.2° grid via optimal interpolation (Hebert et al., 2018).

NERACOOS buoys have recorded temperature, salinity, conductivity, potential temperature, air temperature, air pressure, and wind speed at hourly intervals since 2001 (Morrison, 2019) (Fig. 1). Instruments were positioned at 1, 20, and 50+ m below the surface (site-depending), allowing for characterization of the vertical structure of the water column at high temporal resolution. Salinity and temperature data from buoys N, M, and I from 2001 to 2019 were used to characterize water mass characteristics and climatology.

Models

Hydrodynamic model

The Regional Ocean Modeling System (ROMS) is a free-surface, hydrostatic, primitive equation circulation model with split-explicit time-stepping for the baroclinic and barotropic modes (Shchepetkin and McWilliams, 2005). It has been applied to the GOM, where it is used annually to forecast A. catenella blooms (He et al., 2008). The GOM ROMS configuration includes 36 terrain-following sigma layers, 1-km resolution inshore, and 3-km resolution offshore. The model domain stretches from Georges Bank and Cape Cod in the Southwest to the Bay of Fundy and Halifax in the Northeast (Fig. 1), encompassing 38 to 47°N and 61 to 73°W. Other models have been developed for the region, including an operational model for the Gulf of Maine using the Finite Volume Coastal Ocean Model (FVCOM) (Beardsley et al., 2013; Chen et al., 2013), which has an unstructured grid that can allow for finer resolution of complex coastal geometry. ROMS was chosen for this study, however, because it is part of an existing framework that benefits from extensive model development and assessment in previous circulation and HAB studies in the region (e.g. He et al., 2008; Li et al., 2020, 2015, 2009; McGillicuddy et al., 2014, 2011). In addition, this study focuses on the large-scale low-frequency hydrodynamic variability, for which ROMS and FVCOM solutions are very similar. The choice to use ROMS also lays the foundation to develop future Pseudo-nitzschia models in the same manner as the existing A. catenella model.

Three-dimensional temperature, salinity, and velocity, and two-dimensional sea surface height from HYCOM (HYbrid Coordinate Ocean Model) (Wallcraft et al., 2003) were interpolated to the ROMS grid to create initial and boundary condition files. HYCOM experiment GOFS3.0, with 33 vertical layers and 1/12° horizontal resolution, was used for years 2012–2018, while GOFS3.1, with 41 vertical layers and 1/12° horizontal resolution, was used for year 2019, because GOFS3.0 was only available through November 2018. HYCOM utilizes hybrid vertical coordinates, with isopycnal vertical layers in the open stratified ocean, terrain-following sigma layers in the coastal ocean, and z-coordinates in unstratified areas. HYCOM also assimilates data from the Navy Coupled Ocean Data Assimilation (Cummings, 2005) and estimates surface fluxes with bulk parameterization and data from the Navy Operational Global Atmospheric Prediction System (NOGAPS).

Atmospheric forcing in the GOM ROMS was parameterized via bulk formulation with data from the North American Regional Reanalysis (NARR) (Mesinger et al., 2006), which assimilates data from various sources to calculate air temperature, air pressure, wind speed, cloud fraction, long- and shortwave radiation, relative humidity, and precipitation. Output resolution was 6 h and 1/6°. Five rivers are included in the GOM ROMS forcing files as volumetric transport (m3 s−1) as measured by the U.S. Geological Survey (Interior, 2021): (in order of decreasing transport) St. John River, Penobscot River, Kennebeck River, Androscoggin River, and Merrimack River. A volumetric adjustment was added during forcing file generation to account for drainage area downstream of the gauge.

For each year from 2012 to 2019, GOM ROMS was initialized in February and run through December. These years were chosen to compare four non-P. australis years (2012–2015) with the four years after P. australis’ first appearance (2016–2019). Each year was initialized with HYCOM output to reduce potential error caused by model drift over time. Because of strong tidal variability in the region, model results were saved hourly for use with the particle tracking model.

Particle tracking model

The Lagrangian TRANSport Model (LTRANS) is an offline, three-dimensional particle tracking model (Hunter, 2021; North et al, 2006, 2008; Schlag and North, 2012). It reads in the output from ROMS and calculates particle trajectories through a combination of advection (4th Order Runge-Kutta advection scheme), turbulence, and particle behavior. Turbulence was parameterized with a random walk model in the horizontal and a random displacement model in the vertical, with horizontal diffusivity set to 2 m2 s−1 and vertical diffusivity taken directly from ROMS. Particles in these simulations were passive and neutrally buoyant. LTRANS operates on the native ROMS grid, using bilinear interpolation in the horizontal and spline interpolation in the vertical to interpolate velocities, temperature, and salinity to sub-grid scales. Land-ocean boundaries are reflective, as is the surface boundary, and open ocean boundaries allow particles to escape the domain, but not reenter. For all experiments listed below, the internal time step for particle movement was 2 min, the external time step (from GOM ROMS) was 1 h, and particle tracking results were saved every 6 h.

Particle tracking experiments

A series of reverse and forward particle tracking experiments were designed to explore three related questions: (1) what are the most likely source regions for P. australis, (2) which of the source regions had the strongest connectivity to locations where contaminated shellfish and P. australis were observed, and (3) was the connectivity from source regions to the GOM different in 2016 compared to previous and subsequent years?

LTRANS was run “in reverse” by multiplying the velocities by −1 and reversing model fields (salinity, temperature, and velocity) with respect to time. The accuracy of this method was checked by plotting individual particle tracks without turbulence in the reverse and forward runs and ensuring that they matched. Particles were seeded at each of the 2016 ship survey sample locations in Clark et al. (2019) at 1 m, 10 m, and 20 m, in accordance with sample depths. 2000 particles total were released at each unique (x,y,z) location, according to the method outlined in Simons et al. (2013). The particles were distributed evenly with one release at each location every 6 h during the period of shipboard observations (Oct 5 00:00 to Oct 7 18:00) for a total of 420,000 particles. In addition to focusing on the timing of the 2016 ship survey, this release time period was sufficient for capturing the stochasticity of reverse connectivity, as the particle tracks fully sampled the known inflow pathways in the region (Section 3.2.1). The experiment was repeated in each year from 2012 to 2019 and run backward to the beginning of June (120 days).

The forward experiments were designed to complement, and were guided by, the reverse experiments. A preliminary analysis of the reverse experiment pointed to the inner Scotian Shelf and Northeast Channel as the two most likely source regions for P. australis. For each forward experiment, a latitude-longitude point was made the center of a square with side lengths that were chosen to encompass the region of interest (see Table 1). Particles were randomly distributed horizontally within the square at a density of 3 particles km−2, which was informed by the methods in Simons et al. (2013). At each (x,y) point, particles were released at the bottom, at mid-depth, and at 50 m, 25 m, and 1 m. Particles were released every 6 h for 3 days starting at midnight on July 15, August 1, and September 1 of each year and experiments were run through October 5. A 6-h release interval was chosen to avoid aliasing the tidal signal, and the three months were chosen to include mesoscale and other low-frequency variability. At the Northeast Channel (NEC) 9375 particles were released at one time, and on the inner Scotian Shelf (SS) 15,360 particles were released at one time. This totals 121,875 (NEC) or 199,680 (SS) particles per location-release, 365,625 (NEC) and 599,040 (SS) particles per location (July, August, and September combined), 964,665 particles per year (2 locations and 3 release times), and 7,717,320 particles in all.

Table 1

Locations of particle releases for the forward experiments.

Release Location	Center latitude/longitude	Region Side length (km)
Inner Scotian Shelf	43.4, −64.7	32
Northeast Channel	42.3, −66.0	25

Data analysis

Field data and ROMS

To see if the hydrographic conditions or transport pathways were distinct in 2016, both model output and observations were used to plot water mass properties, velocity, transport, and vertical water column structure. For all analyses, ROMS output was extracted from the nearest grid cell to the observation location, rather than interpolated.

On the Scotian Shelf, ROMS surface and bottom salinity and temperature were averaged over the region between 42.55°N and 44.64°N and 65.53°W and 62.01°W (Fig. 1A). Stratification was defined as the difference in potential density between the bottom and the surface.

ROMS results were further analyzed as transects. Four transects were chosen to capture various upstream locations from the 2016 sample site, where “upstream” was determined by knowledge of the climatological circulation and the results from the reverse particle tracking simulations. The four transects were on the Scotian Shelf (64.66°W, 44.00°N to 63.62°W, 42.9°N), south of Nova Scotia (66.57°W, 42.62°N to 66.02°W to 43.69°N), across the Northeast Channel (66.11°W, 42.02°N to 65.66°W, 42.37°N), and across the mouth of the Bay of Fundy (67.16°W, 44.68°N to 66.26°W, 44.11°N) (Fig. 1B). ROMS u and v velocities were projected in the alongshore direction according to the angle of the coast relative to East such that u > 0 was toward the GOM. In the Northeast Channel, velocities were projected parallel to the channel walls. Transport toward (away from) the GOM was calculated via where i is each grid cell where u > 0 (u < 0), u is the alongshore component of the velocity in grid cell i, A is the grid cell cross-sectional area, and N is the total number of grid cells where u > 0 (u < 0).

Salt transport was calculated via where S is the salinity in grid cell i. Transport estimates were filtered with a weekly running mean before plotting to reduce the tidal signals.

LTRANS

Particle density and particle probability.

To compare connectivity between one source region and multiple destination regions, particle density was calculated by dividing the sum of all particles in the destination region by the region area (l2, where l is the side length defined in Table 1). This method was chosen to enable comparisons between regions of different sizes. The destinations for the reverse run, the inner Scotian Shelf (43.4°N, 64.7°W, l = 32 km) and Northeast Channel (42.3°N, 66.0°W, l = 25 km), were the same as the forward release locations (Table 1).

For intra- and interannual comparisons in connectivity, particle probability was calculated according to where P is the total number of particles initially released from the source region. The destination regions were the Bay of Fundy, centered on 44.7°N, 66.4°W, with a side length of 50 km, and the eastern Maine coast, centered on 44.2°N, 67.9°W, with a side length of 45 km (Fig. 1B). The Bay of Fundy was chosen because it was where the 2016 bloom was first observed. The location along the eastern Maine coast was based on the 2016 ship survey location in Clark et al. (2019) and subsequent reports of shellfishery closures in eastern Maine.

Growth delivery potential.

An estimate of potential growth for P. australis as a function of temperature was calculated for those particles that were in the Bay of Fundy at any point between September 5 and September 19, the two weeks leading up to when shellfish toxicity exceeded the regulatory limit (Canadian Food Inspection Agency, personal communication). Growth rates were also estimated for particles in the Bay of Fundy between September 12 and 19, and results were similar to the two-week time period, so only the two-week time period is reported here.

Temperature-based laboratory growth rates were measured for a Gulf of Maine P. australis isolate obtained from the 2016 toxic bloom. The isolate was maintained in batch culture in sterile GOM seawater (salinity 33 PSU) amended with f/4 nutrients at 13°C under a 12:12-h light:dark cycle with cool-white light at ~150 μmol photons m−2 s−1. For acclimatized growth rates, cells were transferred to a thermal-gradient device (Blankley and Lewin, 1976) while keeping other growth conditions the same. Growth was monitored by measuring in vivo fluorescence daily with a Turner AU-10 Fluorometer (Turner Designs, San Jose, CA). After reaching steady state growth (i.e., for 3 generations) at 13°C, the culture was similarly acclimatized to additional temperatures (7, 9, 11, 13, and 15°C) using a stepwise increase or decrease of 2°C. If positive growth was achieved for a given temperature, the culture was transferred to a new temperature after three generations, even if not in steady state.

To calculate per-particle potential daily growth, each particle’s along-track temperature was extracted on days when it was entirely in the euphotic zone (defined as 100 m according to the data on NASA’s Ocean Color Website, Feldman, 2014) and averaged over the entire day. Euphotic zone depths of 50 and 75 m were also tested with no discernible changes to the results. Growth rate at each daily temperature was assigned via linear interpolation of the growth curve. On days when the particle exited the euphotic zone, the growth rate was set to zero.

Per-particle average potential daily growth rate was calculated over each particle’s trajectory by dividing the sum of daily potential growth rates by the number of days in the experiment. Interannual comparisons of the distributions of average potential daily growth rate were done with the Kolmogorov-Smirnov test. This pair-wise test compares the maximum difference between two empirical distributions to a theoretical maximum difference to determine the likelihood that the two samples are from the same distribution. It does not assume normal distribution.

Model validation

ROMS has been used for many years to study circulation and A. catenella (formerly A. fundyense) blooms in the GOM (Fennel and Wilkin, 2009; He et al., 2008; Li et al, 2009, 2015; López et al., 2020; McGillicuddy et al., 2011). From a comparison between ROMS data and coastal CTD casts and NERACOOS moorings, Li et al. (2009) found that the model accurately captured surface and sub-surface temperature, salinity, and velocity to reproduce the general hydrodynamics along the coast on monthly and seasonal time scales, but that it missed processes on an event-by-event basis and at local scales. Our analysis agrees with this assessment.

ROMS accurately depicted the GOM’s climatic circulation, with cyclonic circulation in the gulf interior, anticyclonic circulation around Georges Bank, inflows via Nova Scotia and the Northeast Channel, and outflows via the Northeast Channel and Great South Channel (Fig. 2, compare with Fig. 1A). The model also depicted the episodic nature of Northeast Channel inflows: at depth (>100 m) near NERACOOS buoy N, both ROMS and the observations showed changes in salinity up to 2 PSU and changes in temperature up to 6°C in as little as a week (Fig. S-1 in the Supplementary Material). The model also captured seasonal and spatial patterns in sea surface temperature: time series of model output vs NERACOOS Buoy measurements and Halifax Station 2 were in good agreement (Li et al., 2009). On the Scotian Shelf, ROMS estimates of water mass properties were in general agreement with observations from the 2016 DFO Summer Survey (Fig. 3). The mean temperature biases in comparison with the July 2016 DFO Survey on the Scotian Shelf were −0.22°C (Surface), 0.47°C (50 m), −1.42°C (100 m), and −1.10°C (Bottom). The mean salinity biases from the DFO Survey in the same time period were 0.30 PSU (Surface), 0.09 PSU (50 m), −0.32 PSU (100 m), and −0.38 PSU (Bottom).

Fig. 2.

Year-long averaged surface (left) and bottom (right) circulation in the Gulf of Maine in 2016 as simulated by ROMS. Output grid was decimated by a factor of ten before plotting. The scale of the velocity arrows is given in the bottom right.

Fig. 3.

ROMS (red circles) and DFO (blue triangles) temperature vs. salinity on the Scotian Shelf at (clockwise from top left) the surface, 50 m, the bottom, and 100 m. ROMS data were extracted from the nearest grid point to the DFO sample location.

Although the mean biases in temperature and salinity over the Scotian Shelf were overall small, point comparisons between the observations and ROMS output revealed local inaccuracies in model predictions. This emphasizes the need to focus on long-term, low-frequency processes. At some locations, ROMS inaccurately simulated salinity: at Halifax Station 2, surface salinity in the model was more saline than observations by up to 1.8 PSU (Supplementary Material Fig. S-3), while salinity at depth was fresher than the DFO observations by up to 0.3 PSU (Fig. 4). The disagreement at depth could be partly due to the scarcity of nearshore observations in the survey, but inshore salinity values in ROMS were often fresher than observations. Despite these inaccuracies, ROMS represented the key spatial gradients in the salinity distribution: increasing salinity with depth and with distance from shore. At NERACOOS buoy I, 50 m modeled salinity was sometimes fresher than observations by up to 1.3 PSU (Supplementary Material Fig. S-4). Finally, ROMS did not capture the rapid increase in salinity (~0.5 PSU over 6 days) observed at 250 m at NERACOOS Buoy M in July 2016 (Clark et al., 2019). However, as mentioned above, ROMS did have rapid changes in salinity at the Northeast Channel and episodic salinity inflows, which are the processes thought to drive the rapid salinity increase at Buoy M.

Fig. 4.

(left) Scotian Shelf bottom salinity as measured during the DFO Summer Survey in July 2016. Black points indicate sample locations, and values are indicated by a color bar on the right. Spaces between sample locations were filled via optimal interpolation; (right) Bottom salinity as simulated by ROMS on July 15, 2016.

ROMS captures key features that are relevant to the research questions in this study. Because the model captures the alongshore flow along the coast of Maine, bloom progression and P. australis connectivity in the gulf interior can be addressed. With the model’s representation of Northeast Channel and Nova Scotia inflows, the P. australis gulf introduction hypothesis can be tested. Third, because the model captures both high salinity pulses at the Northeast Channel and the location and magnitude of the observed near-bottom salinity on the Scotian Shelf in 2016 (Fig. 4), potential origins of the high salinity signal near the coast of Maine in 2016 (Clark et al. 2019) can be explored.

For the purposes of this paper, ROMS is used to draw relative comparisons between years and locations, because large-scale, low-frequency flow is accurate in spite of local disagreements. In addition, large-scale patterns in salinity allow for analysis of nearshore vs. offshore inflow routes and water mass changes with depth. The model is also used to assess interannual variations in connectivity, inflows through the Northeast Channel, and water mass properties on the Scotian Shelf, to compare 2016 to the other years from 2012 to 2019, and to evaluate P. australis introduction hypotheses.

Results

Hydrodynamics from ROMS output and observations

Water mass properties

From the DFO survey in July 2016, bottom waters on the Scotian Shelf were anomalously warm (up to 4 °C) and saline (up to 0.7 PSU) relative to the 1981–2010 climatology (Fig. 5). Summer water mass properties are relevant to the fall bloom because of transport time scales (Section 3.2.1). Although ROMS did not capture the exact location of the anomalies, the model had relatively high salinity water (36 PSU) at a similar location to the observations (south of Halifax), as shown in Fig. 4. Along with the increased salinity at depth, modeled sea surface salinity (SSS) decreased on the Scotian Shelf in 2016 compared to previous years (Fig. 6): the annual average was lowest in 2017 and second-lowest in 2016 during the 8-year time period, with yearly minima in August 2016 and July 2017. Similarly, at Halifax Station 2, surface salinity (5–20 m) in 2016 decreased below any other year in the 8-year time period besides 2013, and mid-depth salinity (50–100 m) in 2016 and 2017 was less variable than in the other 6 years studied (Supplementary Material Fig. S-8). However, given that the temporal resolution of the 2016 observations was coarser than in the other years from 2012 to 2019, it is difficult to assess the significance of these changes.

Fig. 5.

Anomalies relative to the 1981–2010 climatology for (left) salinity and (right) temperature on the bottom of the Scotian Shelf in July 2016 as measured by the annual DFO Summer Survey.

Fig. 6.

Mean sea surface salinity (SSS) on the Scotian Shelf (as simulated by ROMS) vs time, from 2012 to 2019.

Transport

At the same time and location that positive salinity anomalies were observed on the Scotian Shelf in July 2016, modeled volume and salt transport toward the GOM increased as a result of increased velocities and increased salinity (Supplementary Material Fig. S-9). When the salt transport was calculated with just data from offshore (>80 km from the coast), the results were qualitatively similar to the full transect. This occurred simultaneously with the decreased sea surface salinity (Fig. 6), so the overall increase in salinity transport suggests that the offshore, deep salinity anomalies (Fig. 5) dominated the transport signal. (The potential relationship between increased salinity transport yet decreased nearshore salinity in July 2016 is discussed more in Section 4.3.) In addition, transport was calculated at additional transects downstream of the Scotian Shelf transect (not pictured in Fig. 1), and increased velocity and salt transport were consistent features along the shelf in mid-July 2016. This suggests an anomaly that propagated along the outer shelf.

LTRANS

The next section focuses on connectivity and potential growth estimates from the particle tracking experiments, with an emphasis on the August release. August was chosen as an illustrative release date, because the connectivity time scale from the October 2016 sample region to the inner Scotian Shelf and Northeast Channel in the reverse simulation was 30–60 days (see below), or early August.

Connectivity

The 2016 LTRANS reverse simulation indicates connectivity from the 2016 sample region to the Northeast Channel and the inner Scotian Shelf (Fig. 7). Connectivity was quickest to the Northeast Channel, rising at about 20 days post release in early/mid-September, while connectivity to the inner Scotian Shelf increased about 30 days post release, or late August/early September (Fig. 7). Coastal connectivity to the inner Scotian Shelf was roughly consistent once established. Meanwhile, connectivity to the Northeast Channel was more episodic, with peaks near 20, 30, 45, and 65 days, possibly due to the highly variable nature of Northeast Channel inflows. Reverse connectivity to the Northeast Channel was greater than to the inner Scotian Shelf via the coast, but the majority of particles that initially passed through the Northeast Channel ultimately continued upstream onto the outer Scotian Shelf (Fig. 7). Therefore, regardless of inflow pathway, the most likely source region for P. australis was the Scotian Shelf.

Fig. 7.

(Top) Particle tracks from the 2016 reverse simulation on a planar map of the Gulf of Maine. Tracks are color-coded by the first origin region the particle passes through: (blue) inner Scotian Shelf or (yellow) Northeast Channel. The thick brown contour indicates the shoreline, the thick black contours indicate 100 m and 500 m, and the thin black contours indicate 25 m, 50 m, 75 m, every 100 m–1000 m, 2000 m, and 3000 m. (bottom) Particle density for each of the potential transit regions (blue =inner Scotian Shelf; yellow =Northeast Channel) vs day since release from the 2016 reverse simulation. Particle densities are shown, rather than probabilities, because the Northeast Channel area is smaller than the inner Scotian Shelf.

The forward simulations complement the reverse simulations by illustrating connectivity between each source region and the Bay of Fundy, where toxicity was first observed in 2016, and the eastern Maine coast, where P. australis was observed in water samples (Bates et al., 2018; Clark et al., 2019). Regardless of release time, particles were more likely to pass through the Bay of Fundy if released from the inner Scotian Shelf than the Northeast Channel. In 2016 this difference in connectivity was up to an order of magnitude (Fig. 8), and for all 8 years it ranged from a factor of 2 to a factor of 8 (Fig. 9). In contrast, the likelihood of reaching the eastern Maine coast was within the same envelope of variability for particles released on the inner Scotian Shelf and in the Northeast Channel in all years (Figs. 8 and 9). On average, about 70 % of the particles released on the inner Scotian Shelf that passed through the Bay of Fundy also passed along the eastern Maine coast. Connectivity varied interannually (Fig. 9), but from neither the inner Scotian Shelf nor the Northeast Channel did particles in 2016 have an increased probability of reaching either the Bay of Fundy or the eastern Maine coast relative to other years from 2012 to 2019. In addition, weekly variability in connectivity was as large as interannual variability: in a series of tests for 2016, connectivity differences between release dates August 1st, 7th, 14th, and 21st were as large as connectivity differences among the August 1st releases in all 8 years (not shown). Interannual differences are therefore considered insignificant in this context.

Fig. 8.

Probability of a particle being in the Bay of Fundy (left) or E. Maine (right) vs. days since release as a function of release location (blue = inner Scotian Shelf; yellow = Northeast Channel) and release time in 2016 (o = July 15; * = August 01; x = September 01).

Fig. 9.

Probability of being in the (top) Bay of Fundy or (bottom) E. Maine coast vs. days since release for particles released on August 1 on the inner Scotian Shelf (left) and in the Northeast Channel (right). Each year is represented by a different color, which is the same in all four plots. The vertical dashed line in the top plots indicates September 19, the time when toxic shellfish were first observed in the Bay of Fundy.

Because of their proximity to the coast, particles released from the inner Scotian Shelf were associated with warmer and fresher water than particles released in the Northeast Channel (Fig. 10). On average, the particles released in August on the inner Scotian Shelf that passed along the eastern Maine coast were 0.8 °C warmer, 1.2 PSU fresher, and 7 m shallower than those released in the Northeast Channel that passed along the eastern Maine coast (Table 2). In 2016, particles released on the inner Scotian Shelf that passed along the eastern Maine Coast were typically 2 °C warmer, 0.6 PSU fresher, and 14 m shallower than particles released in the Northeast Channel that passed along the eastern Maine coast (Table 2).

Fig. 10.

Planar view of particle tracks in 2016 for particles released on August 1 from (left) the Scotian Shelf and (right) the Northeast Channel. Particle tracks are color-coded by (top) Salinity, or (bottom) Temperature, and colors are defined by the color bars on the right.

Table 2

Summary of average temperature, salinity, and depth (2016 mean and interannual mean) for particles released from either the inner Scotian Shelf or the Northeast Channel that passed through the 2016 ship survey region.

	Temperature Mean		Salinity Mean		Depth Mean
Release Location	Overall	2016	Overall	2016	Overall	2016
Inner Scotian Shelf	11.3 °C	11.5 °C	32.1	31.8	33 m	30 m
Northeast Channel	10.1 °C	9.5 °C	33.3	32.4	40 m	44 m

Potential growth

Growth rates were measured in laboratory experiments for 7, 9, 11, 13, and 15 °C (Fig. 11). As evidenced by high standard error, cells grown at 7 and 9 °C did not achieve steady growth. Nevertheless, they maintained positive growth at both temperatures for multiple transfers and thus mean – rather than acclimated – growth rates were calculated. Growth rates were extrapolated beyond the 7–15 °C range using the rate of change of growth with temperature (0.02 d−1°C−1 if T < 12 and −0.03 d−1°C−1 if T > 12). These rates were derived from the growth response of P. australis maintained at 13°C following short-term exposure to temperatures <7°C and >15 °C. To account for the variability across temperatures as well as the uncertainty associated with testing a single P. australis isolate, a growth curve was fit to the growth rates at 4, 7, 9, 11, 13, 15, and 18°C according to the function given in Thomas et al. (2012): T is the temperature in degrees Celsius, w is the thermal niche width, and z is the location of the maximum of the quadratic portion of the function, or the temperature at which the growth curve is tangent to the Eppley Curve (Eppley, 1972). Parameter a is the growth rate estimate at 0°C, and b is the Eppley coefficient, or the slope of the growth curve where growth rate rises as a function of temperature. Parameters a, b, z, and w were fit within defined bounds. The bounds for a and b were [0,1], which was informed by the values given in (Norberg, 2004). The bounds for z were [7,15] and the bounds for w were [10,20], which were based on the range of temperatures at which growth rates were measured. The best fit by nonlinear least-squares was achieved when a = 0.198, b = 0.0746, z = 7, and w = 15.75. These a and b values are the same order of magnitude as in Norberg (2004), and these z and w values are reasonable given the temperature values used in laboratory experiments. The sum of squared errors from this fit was 0.0026. The full growth curve is given in Fig. 11 and the measured growth values are given with corresponding standard errors in the Supplementary Material, Table S-1.

Fig. 11.

Laboratory-based growth rates as a function of temperature as measured for a GOM P. australis isolate. The blue points indicate values that were measured, and the error bars show standard error. The orange diamonds were estimated using a rate of change of growth rate with temperature as measured during short-term exposure experiments. The black line shows the theoretical curve as estimated from the function given in Thomas et al. (2012).

Because particles released from the inner Scotian Shelf were in warmer water, they had higher average potential growth rates than those released from the Northeast Channel (Fig. 12). The largest potential growth rates generally occurred nearshore and in the Bay of Fundy, with occasional elevated growth rates farther offshore (not shown). The exception to this is the year 2012, the year of a marine heat wave (Pershing et al., 2015), when potential growth rates were higher at all locations and there was no cross-shore gradient (not shown).

Fig. 12.

Histograms of particle-averaged potential growth rates in 2016 for particles released on the inner Scotian Shelf (top) and in the Northeast Channel (bottom) that were in the Bay of Fundy between September 5 and September 19, 2016. Note that the y-axis scales are different.

The Kolmogorov-Smirnov test indicated that potential growth rates from the inner Scotian Shelf were significantly higher than from the Northeast Channel, whether the growth rates were from all years combined or from 2016 individually. However, potential daily growth rates were neither significantly higher nor significantly lower in 2016 compared to other years from 2012 to 2019, regardless of release location.

Discussion

The particle tracking model, physical model, and in situ data provide clues regarding P. australis inflow route and timing, as well as changes to physical processes driving the 2016 bloom. The following section uses connectivity results from LTRANS and physical parameters from ROMS to determine the most likely inflow route. It also explores changes upstream that could have caused the 2016 event and hypothesizes the risk of regional P. australis HABs in the future.

Likely P. australis introduction pathway

The LTRANS reverse results indicate a P. australis source region on the Scotian Shelf, and inflow routes via either the Northeast Channel or around the southern tip of Nova Scotia. The current along the latter inflow route is not well defined in the literature, having been referred to as both “an extension of the Nova Scotia Current” (Townsend et al., 2015) and the beginning of the Gulf of Maine Coastal Current (Pettigrew et al., 2005). This route will therefore be called the coastal route for simplicity. In 2016, particles released in the Northeast Channel were as likely to reach the eastern Maine coast as those released on the inner Scotian Shelf, the majority of which were transported along the coastal route (Fig. 8). From this alone, both the Northeast Channel and the coastal route were equally likely inflow routes for P. australis cells in 2016. However, particles released from the Northeast Channel had a low probability of reaching the Bay of Fundy compared to those released from the inner Scotian Shelf (Fig. 8). An additional reverse experiment was run with particles released just below the surface in the Bay of Fundy from September 12 to September 19, which was the period of shellfish closures in that region. It confirmed that particles in the Bay of Fundy were connected to the inner Scotian Shelf via the coastal route, but were not connected to the Northeast Channel (Supplementary Material Fig. S-12).

In addition to connectivity, along-path temperature and salinity provide important context. Particles that entered the gulf via the Northeast Channel were generally in a higher salinity range (Fig. 10) than those that entered via the coastal route, although the salinity average was not as high as the coastal observations in Clark et al. (2019). However, particles that entered via the coastal route were more likely to be between 7 and 15°C, the range in which the GOM strain of P. australis can sustain growth, than particles that entered via the Northeast Channel. The highest potential growth values occurred near the coast and in the Bay of Fundy, where water is shallower and warmer. As a result, potential growth rates were higher for particles coming in via the coastal route than for those coming in via the Northeast Channel.

Thus, there are arguments for and against both inflow routes, and particles entering both via the coastal route and the Northeast Channel connected to the eastern Maine coast. Northeast Channel particles were in a salinity range closer to the 2016 observations, but they lacked connectivity to the Bay of Fundy and were in relatively cool water. Meanwhile, coastal route particles had the strongest connectivity to the Bay of Fundy and were in a temperature range that promoted growth, but were in relatively fresh water. Despite the fact that the modeled salinity did not agree with the 2016 coastal observations, the warmer temperatures and connectivity to the Bay of Fundy mean that inflows via the coastal route were more likely to introduce P. australis than inflows via the Northeast Channel. In order for the cells to have been introduced to the Bay of Fundy before the first observation of toxic shellfish in mid-September 2016, they would have had to come from the Scotian Shelf sometime in or before early August 2016 and entered the Gulf of Maine near the coast of Nova Scotia.

2016 in an interannual context

A distinct aspect of the 2016 P. australis bloom compared to previous Pseudo-nitzschia blooms of other species in the Gulf of Maine is that it occurred in late September to October, while other observed blooms associated with DA had occurred from July to early September (Clark et al., 2019; Fernandes et al., 2014). The September bloom could be related to a “marine heat wave” that took place in 2016 across the GOM and Scotian Shelf (Pershing et al., 2018). From ship observations in 2016 (Clark et al., 2019, Fig. 10B), sea surface temperatures along the coast of Maine in early October 2016 were warmer than the 2001–2016 October mean as measured at NERACOOS Buoy I, and as many as half of the measurements were more than one standard deviation warmer than the 2001–2016 October mean. The average water temperature from the ship observations was 13.5 °C, which is near the optimal temperature for the tested strain of P. australis (Section 3.2.2, Fig. 11). Increased temperatures have also been associated with increased DA production on the U. S. West Coast (Trainer et al., 2020). Temperature has also been found to contribute to regime shifts between species of Pseudo-nitzschia near Denmark (Lundholm et al., 2010), and between zooplankton communities along the United States East Coast (Morse et al., 2017), although the shifts occurred over decades, not years. Therefore, increased temperatures may have been important for both growth and DA production in the GOM in 2016. The 2016 marine heat wave might also have selected against other organisms, whether competitive or predatory, to allow P. australis a greater chance to survive and grow. Water temperature effects on inter-community competition and long-term trends in the GOM are speculative and warrant further research.

Warmer temperatures may have supported a later bloom and more particles with relatively high growth potential in 2016, but the connectivity and potential growth rates that year were not significantly different from the other years between 2012 and 2019. Although the warm temperatures could have provided favorable growth conditions at other locations within the Gulf of Maine, they did not correspond with improved growth delivery potential along the inflow routes. It should be noted that the growth model used here incorporated only temperature-dependent growth rates. Although growth rates as a function of salinity have been reported for other strains of P. australis (Ayache et al., 2020), they are not currently available for this strain. To improve our understanding of bloom dynamics in the region, it is important to explore the growth rate estimates as a function of salinity, irradiance, and nutrients, as these can vary considerably depending on when and where cells occur. Expanding the model, however, is beyond the scope of this study.

Changes upstream

On the Scotian Shelf in 2016, two changes occurred simultaneously: an increase in bottom salinity and a decrease in surface salinity. The increased bottom salinity, and its connection to elevated nearshore salinity in 2016 as explained in Clark et al. (2019), led to the “initial introduction” of P. australis hypothesis that motivated this study. The following section will explore this hypothesis more carefully, evaluate each salinity change separately, and investigate how the salinity anomalies relate to both each other and the 2016 DA event.

Increased bottom salinity on the Scotian Shelf

A notable change on the Scotian Shelf in 2016 was the increased bottom salinity in the summer, which was captured both by ROMS and by the DFO summer survey. Brickman et al. (2018) observed and modeled these anomalies from 1990 to 2015 and argued that they occur as a result of warm, saline Gulf Stream water cutting off cool, fresh Labrador Current water near the Grand Banks off the coast of Newfoundland, Canada. Positive anomalies increased in frequency from 2006 to 2015, leading to more positive temperature and salinity anomalies at depth on the Scotian Shelf, but this might be the result of interdecadal variability in the region, rather than an indicator of a long-term trend (Brickman et al., 2018).

The timing of the increased salinity on the Scotian Shelf (Fig. 4) and the connectivity from the Northeast Channel to the eastern Maine coast support the hypothesis that a salinity anomaly entered via the Northeast Channel. Brickman et al. (2018) estimated that salinity anomalies on the Scotian Shelf propagate southwest at about 150 km month−1. The distance from the Scotian Shelf transect (where the increased salt transport was observed) to the Northeast Channel is about 230 km, which means that the observed pulse would have taken about 1.5 months (45 days) to travel from the transect to the Northeast Channel, arriving at the Northeast Channel in mid-to-late August. From the model output, salt transport toward the Gulf of Maine across the Northeast Channel in mid-August increased more than one standard deviation above the interannual mean (as it did on the Scotian Shelf) and remained elevated for about 5 days. This salt transport was not larger than all other years from 2012 to 2019, and some years besides 2016 had summer salt transports of similar magnitude, but 2016 was the only year with elevated summer transport on both the Scotian Shelf and the Northeast Channel. Meanwhile, connectivity between the coast and the Northeast Channel in the reverse experiment in 2016 had a peak near 45 days post-release, or mid-August. This suggests that a water mass that was present on the Scotian Shelf in July and entered via the Northeast Channel in August could have arrived at the coast of Maine by early October, in time to be observed during the 2016 ship survey (Clark et al., 2019).

This discussion has focused on the timing of a particular anomaly observed in mid-July 2016, but as described by Brickman et al. (2018), multiple anomalies can form on the Scotian Shelf in a summer. Perfect timing is therefore not necessary for this introduction route to be plausible. Transport time varies with depth as a function of velocity shear: transport time scales from the Northeast Channel to the eastern Maine coast were calculated as a function of particles’ initial depth and found to range from 45 days (surface to 50 m release) to 90 days (bottom release). From the DFO survey data, anomalous salinity signals were present from 50 to 125 m. Considering the range of possible depths and associated transport time scales, an anomaly passing through the Northeast Channel any time between early June and mid-August 2016 could explain the elevated salinity values at the coast.

Brickman et al. also hypothesized that repeated salinity anomalies would, through an integrative effect, result in deeper on-shelf penetration via channels such as the Northeast Channel. Considering that 2016 followed several years of increased bottom salinity anomalies on the Scotian Shelf (Brickman et al., 2018), it is not unreasonable that an anomaly would have reached the coast of Maine in 2016, causing the high salinity values observed in Clark et al. (2019). Given these lines of evidence, the saline water mass that was observed on the coast of Maine in 2016 likely came from the Gulf Stream, propagated along the outer Scotian Shelf, and entered the GOM through the Northeast Channel.

Decreased sea surface salinity on the Scotian Shelf

As the bottom salinity on the Scotian Shelf increased in 2016, the surface salinity decreased by up to 0.6 PSU. Although the increased bottom salinity and decreased surface salinity occurred simultaneously in summer 2016 (Figs. 5 and 6), bottom salinity and surface salinity on the Scotian Shelf were not consistently correlated across the 8 years studied. For each year, linear regressions were calculated between surface nearshore salinity and bottom offshore salinity at the ROMS eastern boundary (on the Scotian Shelf), and R2 values ranged from 0.01 to 0.67, with an average of 0.16. Therefore, the surface and bottom salinity anomalies do not appear to be caused by the same process.

The modeled decrease in sea surface salinity on the Scotian Shelf in 2016 and 2017 (Fig. 6) might suggest increased outflow from the St. Lawrence River via the HYCOM boundary conditions, because upper layer salinity (30–50 m) on the Scotian Shelf is influenced by St. Lawrence River discharge (Dever et al., 2016). However, St. Lawrence River discharge was at an 8-year minimum in the 12 months prior to summer 2016 (Interior, 2021), and it has been shown that the salinity signal from the St. Lawrence River decreases to ± 0.1 PSU by the time it reaches the Halifax Line (Ohashi and Sheng, 2013). It is therefore unlikely that changes in the St. Lawrence River caused the sea surface salinity changes on the Scotian Shelf. Townsend et al. (2015) argued that there had been an increase in the frequency of strong Scotian Shelf Water inflows via Nova Scotia in the early 2010s, and presented evidence of increasing freshwater fluxes from the Labrador Sea via the Labrador Current. It is possible that the surface freshening modeled on the Scotian Shelf in 2016 and 2017 is indicative of these episodic changes. This is in contrast to the findings of Brickman et al. (2018), which outlined how the Gulf Stream effectively cuts off the Labrador Current near the Grand Banks, but this difference could be because the processes Townsend et al. described occur at the surface, while those described by Brickman et al. occur at depth.

Intermediate layer water on the Scotian Shelf (50–100 m, also known as the Cold Intermediate Layer, or CIL) is composed of either the inshore branch of the Labrador Current or CIL water formed in the Gulf of St Lawrence, which is distinct from St. Lawrence River outflow (Dever et al., 2016). In the intermediate layer at Halifax Station 2, observed salinity and temperature were both lower and less variable in 2016 than in the 7 other years in this study, with a water mass signature that was a mixture of Labrador Slope Water, Cabot Strait–Cold Intermediate Layer water, and Cabot Strait Subsurface water, as described in Dever et al. (2016) (Fig. 13). This suggests that changes on the Scotian Shelf are driven by changes in the Labrador Current or Gulf of St. Lawrence, either of which could be a link to the upstream source of P. australis. More work is necessary – either via an expanded model domain or ship surveys – to determine the relative likelihood of each source.

Fig. 13.

Year-round temperature vs salinity as measured at Halifax Station 2 at 85 m. 2016 is highlighted in green, while 2012–2018 are plotted in black. For comparison, the endmembers of CBS–CIL (Cabot Straight–Cold Intermediate Layer), inLC (Inshore Labrador Current), CBSS (Cabot Straight Subsurface Water), and LSW (Labrador Slope Water), as defined in Dever et al. (2016), are plotted. Contours show lines of constant density and are labeled accordingly.

To summarize, P. australis likely entered the Gulf of Maine via the inner Scotian Shelf and coastal route around southern Nova Scotia, was concurrent with decreased sea surface salinity on the Scotian Shelf, and had an upstream source in either the Labrador Current or Gulf of St. Lawrence. In contrast, the observed positive salinity anomalies on the Maine Coast likely originated in the Gulf Stream and entered the GOM via the Northeast Channel. Particle tracking experiments suggest that these two phenomena met along the upper coast of Maine and were mixed through strong tidal mixing, which is not unreasonable according to the data. If the positive salinity anomaly at the bottom (mean = 35) and negative salinity anomaly at the surface (mean = 32.5) were to mix, the required mixing ratio to achieve the salinity observed in Clark et al. (2019) (mean = 33.25) would be 30/70 (bottom/surface). The mixed water could have then flowed southwest where it was observed during the 2016 ship survey. For a more thorough analysis of environmental salinity and its correlations with different Pseudo-nitzschia species in the GOM during this period, the reader is referred to Clark et al. (2019).

It should be noted that several studies have found salinity to influence Pseudo-nitzschia growth (Ayache et al., 2020; Doucette et al., 2008; Thessen et al., 2005) and DA production (Ayache et al., 2018; Thessen and Stoecker, 2008). In each of these studies, Pseudo-nitzschia growth and DA production increased with salinity, suggesting that the elevated salinity in 2016 might have played a role in the DA event. However, the year-to-year salinity differences discussed in this paper (order 1 PSU) are much smaller than those examined in the literature (order 5–10 PSU), so salinity is not thought to be the leading causative factor of the 2016 P. australis bloom or the associated DA event. This is discussed in more detail in Clark et al. (2019).

Changing upstream water masses could also have affected nutrient concentrations in the GOM, and, in turn, the 2016 bloom. The surface and mid-depth salinity changes suggest increasing influence from the Labrador Current, which typically has a silica-to-nitrate ratio greater than 1 (Townsend et al., 2010). The bottom salinity anomalies, meanwhile, are thought to have originated in the Gulf Stream, which has a silica-to-nitrogen ratio less than 1 (Townsend et al., 2010). It is reasonable from a nutrient perspective that P. australis were carried in with the surface Scotian Shelf inflows, in that an increased silica-to-nitrogen ratio could have supported a diatom population. In the 2016 ship survey, however, the silica-to-nitrogen ratio was significantly lower compared to two of the previous years’ ship survey samples, and this decrease was thought to exacerbate DA production (Clark et al., 2019). The low silica-to-nitrate ratio in the 2016 observations could be the result of Gulf Stream water mixing with Scotian Shelf water at the mouth of the Bay of Fundy, as hypothesized above. Thus, the two-pronged introduction hypothesis (P. australis via the coastal route and a positive salinity anomaly via the Northeast Channel), in combination with the nutrient characteristics of the source water masses, suggest that P. australis was advected into the Bay of Fundy and eastern Maine, where it met saltier waters with low silica, enhancing DA production.

P. australis distribution in the North Atlantic

Evidence suggests that P. australis originated outside the GOM, which necessitates a look at P. australis distributions in the North Atlantic. P. australis has been observed near the United Kingdom and Ireland (Bresnan et al, 2015, 2017; Cusack et al., 2002; Fehling et al., 2006; Hasle, 2002; Thorel et al., 2014), France (Ayache et al., 2020; Husson et al., 2016; Klein et al., 2010; Lema et al., 2017; Thorel et al., 2017), the Iberian Peninsula (Churro et al., 2009; Palma et al., 2010; Zapata et al., 2011), and Morocco (Ennaffah et al., 2012). The question remains how the cells could have been transported from these regions in the eastern Atlantic to the Scotian Shelf and Gulf of Maine.

The two likely transport routes are via the subtropical gyre and the Gulf Stream, or via the subpolar gyre and the Labrador Current, which meet near the Grand Banks of Newfoundland (Townsend et al., 2006). Introduction via the subtropical gyre would indicate a more southern origin of the species, possibly from Morocco or the Iberian Peninsula. To have been introduced along this route, the cells would have been carried west across the subtropical Atlantic and then northward by the Gulf Stream. However, growth rates of the 2016 GOM strain of P. australis declined rapidly at temperatures above 18 °C (Fig. 11), which is regularly exceeded along this route. In addition, we have argued that the source of P. australis was likely either the Labrador Current or Gulf of St. Lawrence. The strong connectivity from the inner Scotian Shelf to the Bay of Fundy and the decreased surface salinity on the Scotian Shelf in 2016 align with introduction via the subpolar gyre and Labrador Current. This suggests a cooler origin, which better aligns with the temperature growth curve for the GOM P. australis isolate. For comparison, P. australis from the English Channel grew best at 13.5–18.6°C, with growth rates from 0.47 to 0.83 day−1 (Thorel et al., 2014). This temperature range overlaps with the range of peak growth for the GOM strain of P. australis, but it is slightly warmer and the growth rates are slightly higher. This could be due to population variability in the GOM that was not represented by examining a single GOM strain under the chosen laboratory conditions. Alternatively, selection may have occurred en route to the GOM. A subpolar introduction route seems more likely, but this hypothesis can be neither proven nor refuted without additional data. The best way to test this hypothesis is via targeted sampling, further culturing or modeling, and a genetic comparison between the GOM P. australis strain and other strains across the Atlantic, which could be the subject of future research.

Looking to the future: persistence vs. persistent introduction

An important question for scientists, managers, and state officials is whether the P. australis bloom in 2016 was anomalous or if it indicates a growing DA concern for the region. The 2016 event was likely due to the import of P. australis from the Scotian Shelf, but future blooms could be seeded by cells retained in the GOM, recurring introductions of cells from outside the GOM, or some combination of both. In other words, do we expect introduction plus persistence, or persistent introduction?

P. australis cells and DA-related shellfish harvest area closures have occurred in Maine every year since 2016 (Maine Department of Marine Resources, personal communication). P. australis cells are classified as a large-sized Pseudo-nitzschia species, and concentrations of large Pseudo-nitzschia cells at the Bar Harbor monitoring station exceeded 100,000 cells L−1 in 2017, 2019, and 2020, with species-specific data being generated as part of an ongoing study. Low cell concentrations in 2018 could be the result of a spatially heterogeneous bloom in which cell concentrations were higher elsewhere, because data post-2016 are from a single-point time series.

However, despite observations of P. australis since 2016, the model output in this study did not point to any shifts that could explain why there were no DA events pre-2016 but annual DA events post-2016. Surface salinity was low in 2016 and 2017 and this corresponded with the original P. australis bloom in 2016 (Section 4.3.2), but the signal did not persist past 2017. In addition, model output did not indicate lasting changes in connectivity strength, connectivity timing, or growth potential. Evidence from the literature does point to long-term trends of increased inflows of Scotian Shelf Water (Townsend et al., 2015) and increasingly warm/saline near-bottom water on the Scotian Shelf (Brickman et al., 2018). This study did not identify any regime shift in water mass characteristics associated with the 2016 bloom, but other studies suggest changes in GOM inflows.

The evidence post-2016 indicates that P. australis blooms have continued in the GOM, but the available data are insufficient to determine whether the blooms are due to persistence or persistent introduction. Continued and expanded time series monitoring is recommended to map future blooms as they occur and to better understand interannual variability in the timing, extent, and severity of DA events in the region. Ideally this expanded effort would include sample collection on the Scotian Shelf and regions further upstream to better understand connectivity with established P. australis populations.

Conclusions

This paper used a physical model, particle tracking model, and observations to explore introduction and connectivity pathways to the Gulf of Maine in 2016, the year of a regionally historic DA event. Results from the particle tracking simulations revealed that P. australis cells likely originated on the Scotian Shelf and were carried into the GOM via inflows south of Nova Scotia, because only this pathway explains the highly toxigenic P. australis cells and observed toxicity in the Bay of Fundy. Anomalously saline water originated near the bottom on the Scotian Shelf, entered the GOM via the Northeast Channel, and mixed with lower salinity surface water near the mouth of the Bay of Fundy. Regional connectivity and potential growth rates in 2016 were not larger than in the other modeled years between 2012 and 2019, but surface salinity on the Scotian Shelf was up to 0.6 PSU lower than the 7 other years in this study, and bottom salinity on the Scotian Shelf was 0.7 PSU greater than the 1981–2010 climatology. This suggests that the bloom was controlled by large-scale processes, not local-scale ones, and we propose that P. australis arrived at the Scotian Shelf and GOM via either the Gulf of St. Lawrence or the Labrador Current. Future studies should consider adding nutrient-, light-, or salinity-dependency to the growth model, exploring mechanisms underlying DA production in the GOM, expanding the hydrodynamic model domain to include offshore dynamics, and running the model continuously over more years to test the question of interannual persistence. Ocean sampling is recommended on and beyond the Scotian Shelf to determine if there is an upstream source of P. australis. Finally, P. australis blooms are likely to remain a concern in the GOM, and continued monitoring is warranted.