https://orcid.org/0000-0003-4984-4143
Figures
Abstract
Euphausiids are a keystone species in coastal food webs due to their high lipid content and seasonally high biomass. Understanding the habitat and environmental drivers that lead to areas of high biomass, or ‘hotspots’, and their seasonal persistence, will support the identification of important foraging regions for mid- and upper- trophic level predators. We quantify the distribution of hotspots of the two dominant species of euphausiid in the north-east Pacific Ocean: Euphausia pacifica and Thysanoessa spinifera, as well as euphausiid larvae (mixed species). The Canadian coast encompasses the northern California Current Ecosystem and the transition zone to the Alaska current, and is a highly productive region for fisheries, marine mammals, and seabirds. We used spatiotemporal modelling to predict the distribution of these three euphausiid groups in relation to geomorphic and environmental variables during the important spring-summer months (April through September) when euphausiid biomass is highest. We quantified the area, intensity, and persistence of biomass hotspots across months according to specific oceanographic ecosections developed for marine spatial planning purposes. Persistent hotspots of both adult species were predicted to occur along the 200 m depth contour of the continental slope; however, differences were predicted on the shallower Dixon shelf, which was a key area for T. spinifera, and within the Juan de Fuca Eddy system where E. pacifica hotspots occurred. The continental slope along the west coast of Vancouver Island was the only persistent hotspot region common between both adult species and euphausiid larvae. Larval distribution was more correlated with T. spinifera than E. pacifica biomass. Hotspots of adults were more persistent across months than hotspots of euphausiid larvae, which were seasonally patchy. The persistence of biomass hotspots of forage species through periods of low overall biomass could maintain trophic connectivity through perturbation events and increase ecosystem resilience to climate change.
Citation: Evans R, English PA, Anderson SC, Gauthier S, Robinson CLK (2021) Factors affecting the seasonal distribution and biomass of E. pacifica and T. spinifera along the Pacific coast of Canada: A spatiotemporal modelling approach. PLoS ONE 16(5): e0249818. https://doi.org/10.1371/journal.pone.0249818
Editor: Andrea Belgrano, Swedish University of Agricultural Sciences and Swedish Institute for the Marine Environment, University of Gothenburg, SWEDEN
Received: January 22, 2021; Accepted: March 26, 2021; Published: May 14, 2021
Copyright: © 2021 Evans et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Zooplankton biomass data are available via Open Portal from https://open.canada.ca/data/en/dataset/9447ecf8-a7f7-4904-8ab0-3c597c534c4b The bathymetry dataset from which all geomorphic variable were calculated can be found at: https://bcmca.ca/datafeatures/eco_physical_bathymetry NOAA chlorophyll and SST data are open access and are referenced in the text.
Funding: S.C.G and C.LK.R received funding. This research was funded by Fisheries and Oceans Canada Aquatic Climate Change Adaptation Services Program (ACCASP).
Competing interests: The authors have declared that no competing interests exist.
Introduction
Euphausia pacifica and Thysanoessa spinifera are the most abundant euphausiid species in the north east Pacific Ocean [1–3]. Both species comprise a large proportion of the diet of numerous forage and commercially important fish species in this area, such as sardine (Sardinops sagax), anchovy (Engraulis mordax), Pacific herring (Clupea pallasii), Pacific hake (Merluccius productus) and salmon (Oncorhynchus spp.) [4], and a large proportion of seabird [1, 5] and whale diets [6, 7]. Euphausiids are an energy rich food source due to a high lipid content [8], and are seasonally abundant forming large aggregations or ‘swarms’. Zooplankton in general are considered indicator species of their environment as they show high sensitivity to changes in physical parameters due to short life spans. Large variability in the spatial and temporal distribution of E. pacifica and T. spinifera has been linked to seasonal changes in the timing and intensity of the spring transition, which is reflected by ocean temperature and productivity cycles [2, 9–11].
Euphausiid research in the California Current ecosystem (CCE) has focused on the physical mechanisms that lead to regions of elevated biomass, or ‘hotspots’, of euphausiids e.g. [12–18]. The distribution, intensity and temporal persistence of these areas has been quantitatively linked with local euphausiid predators [17, 19–24]. Physical processes that lead to hotspot formation tend to be localised and are linked to the timing and local intensity of upwelling and the complexity of seafloor topography [16]. Bathymetrical edges, such as along continental margins and around the edges of canyon systems, are regions where hotspots of euphausiids often occur [7, 15, 16, 25] and are therefore also often important foraging regions for their predators [6, 26]. Both E. pacifica and T. spinifera have been recorded to occur in order-of-magnitude larger densities along steep canyon edges than in surrounding water [7, 27]. E. pacifica is approximately 100 times more abundant than T. spinifera [10, 28], and they exhibit differences in life history, reproductive timing, organism size and habitat associations [e.g. 3, 10, 12, 29], which control their distributions.
The Canadian Pacific coast encompasses the transition zone between the California Upwelling Domain and the northerly Alaska Downwelling Domain, although a short upwelling season still occurs along the central and northern British Columbia (BC) coast [30]. The Juan de Fuca canyon off south west Vancouver Island is part of the upwelling domain and is an important region for euphausiids [25, 31, 32]. Whether hotspots of euphausiids are spatially and temporally persistent is unknown. The distribution of adults and larvae is spatially patchy and spatial segregation by age classes occurs perpendicular to the coast [3, 31, 32], with less-motile larvae more vulnerable to changes in current patterns [3]. As in the California Current ecosystem, E. pacifica is thought to be more numerically dominant overall, with T. spinifera more abundant over shallower areas on the shelf, and E. pacifica more abundant in deeper waters over the continental shelf and slope [12, 29, 33].
Monitoring of euphausiid dynamics for the Canadian coast has been spatially and temporally patchy, with a focus on the west coast of Vancouver Island [3, 34, 35]. Species-specific studies of the spatial structure of euphausiid hotspots in the northeast Pacific Ocean are few [but see 14]. It is currently unknown whether structural differences occur between aggregations of the two species, whether the distribution of these hotspots varies between species, and how temporally persistent the hotspots are after the spring transition when biomass is highest. The spatial structure and temporal patchiness of euphausiid larvae hotspots is even less well studied. Therefore, the main aim of this work was to identify areas of persistently high biomass of the two dominant species of euphausiids, and euphausiid larvae, using a long-term euphausiid net-capture dataset for the whole BC coast. Spatiotemporal modelling using geostatistical methods has been shown to generate more precise and accurate biomass and abundance predictions than design-based methods, or other modelling approaches used for species distribution modelling [36, 37]. Through these models, we identify 1) the physical habitat characteristics shaping euphausiid distribution; 2) where biologically important areas occur for the two species and euphausiid larvae; and 3) whether these hotspot areas are persistent through the spring and summer months. Given the known effects of complex bathymetry for balancing advective and retentive currents [16], we hypothesize that euphausiid hotspots will occur along depth isobars in areas such as the Juan de Fuca Eddy and canyon system and the continental slope. We compare the distribution of predicted hotspots across oceanographic ecosections for the British Columbia (BC) coast, which were developed for marine spatial planning purposes [38].
Materials and methods
Euphausiid biomass data
The Canadian government has collected net-based zooplankton samples to estimate abundance and biomass for the British Columbian coast for over 30 years [~1986-present, 32], as part of a variety of sampling programs. Prior to 1993 the spatial and temporal resolution of sampling was low; therefore, we truncated these years from the dataset prior to analysis. Data were collected by both vertical hauls and oblique hauls, which were considered comparable, however this will add some bias to results. We analyzed data from daylight hauls, as these were more numerous than night hauls. A correction factor [33], was applied to daylight counts to account for net avoidance by euphausiids.
In general, nets were lowered to near the bottom (~+/-5m) or to a maximum of 250 meters. Within the dataset there were some points with large disparity between bottom depth and sampling depth, and these were excluded from further analysis: at a bottom depth of less than 100m, sampling points were excluded if they were over 30m from the bottom, while at bottom depths between 100m and 250m, sampling points were excluded if they were over 50m from the bottom. We assigned sampling sites to an oceanographic ecosection (Fig 1), based on similar biological and physical characteristics, developed as part of the British Columbia Marine Ecological Classification System (BCMEC) [39]. Due to known differences in currents along the slope, the continental slope ecosection was further divided into north and south bioregions using boundaries used by the Canadian government for marine spatial planning purposes [38].
Coloured sections show the boundaries of the BCMEC oceanographic ecosections. Grey lines indicate the 200 and 1000 m isobars. The solid black line shows the division of north and south bioregion sections [38]. The dashed black lines indicate regions of more complex bathymetry including the Juan de Fuca canyon and Eddy system on the Vancouver Island shelf and three sea valley and canyons systems on the west edge of Queen Charlotte Sound.
For all stations sampled, formalin-preserved zooplankton samples were sub-sampled to contain ~300 individuals and identified to species level using a dissecting microscope. Euphausiids were separated from the rest of the plankton sample and measured. Larval stages were assigned a “stage” group comprised of eggs: euphausiid eggs and nauplii <3mm, s1: euphausiid zoea <5mm, s2: euphausiid juveniles ≥5mm<10mm, and S4: euphausiid ≥10mm, which were considered adults and separated into males and females. For volume calculations, a flow meter was used when available (+90% of samples), and sample volumes were compared with the expected value calculated from the amount of wire put out during sampling as a quality assurance check. Multiplicative species and stage-specific size coefficients (estimated mg dry weight per individual) were used to convert euphausiid abundance counts to biomass estimates [see Ref. 40 for coefficients]. Biomass values for each stage and sex were standardized via the volume to mg/m3. We later converted the euphausiid biomass from mg/m3 to vertically integrated dry weight biomass mg/m2 based on sampling depth, which is considered more robust to variability in sampling depth than cubic biomass estimates [41].
The majority of voyages occurred over the spring and summer months between April to September. This coincides with the important transition period for the upwelling domain of the coast [42] and the peak euphausiid biomass period [12, 29]. We focused on these months for further analysis. A map showing the distribution of sampling points across months (n = 1609) as modelled can be viewed in the (S2 Fig). The contribution of adult E. pacifica and T. spinifera to overall biomass was calculated as a ratio of adult euphausiid biomass (E. pacifica:T. spinifera).
Spatiotemporal modelling of seasonal changes in biomass
We estimated E. pacifica and T. spinifera biomass using spatiotemporal models to explore differences in the environment and habitat preference of these two dominant euphausiid species. Larval stages of euphausiids are difficult to differentiate; therefore, a large proportion of larval stages in the database were not identified to the species level. Because of this, we were unable to separately model the distribution of species-specific larval stages; instead, we modelled the biomass of all larval stages combined. For all models, we chose biomass over abundance, as it is more ecologically relevant in terms of predator consumption.
The contribution of many sampling programs meant that the euphausiid biomass dataset was unbalanced and irregular both spatially and across seasons. We opted to account for this irregular sampling design using a model-based approach rather than a design-based approach, by fitting a geostatistical model that makes use of a predictive-process stochastic partial differential equation (SPDE) mesh. This ‘mesh’ approximates a Gaussian Markov random field (GMRF) through a set of representative locations (‘knots’) and bilinear interpolation between those knots [43]. We fit the models via the R package ‘sdmTMB’ [44, 45], which interfaces the integrated nested Laplace approximation [INLA; 46] with Template Model Builder [TMB; 47] to find the maximum marginal likelihood while integrating over spatiotemporal random effects. This method is increasingly popular in ecology [37, 48–50]. Using this approach, biomass is modelled as a function of both ‘fixed’ effects, as a result of explicit habitat variables such as depth, and of ‘random’ effects as a product of unobserved or ‘latent’ spatiotemporal effects using Gaussian Markov random fields. There is evidence that such model-based standardization of survey biomass index trends better encapsulates in situ spatial correlation and processes than design-based estimates [30, 41]. A thorough outline of the model framework and underlying statistical structure is provided [44]. We made some modifications to this framework, which are discussed in more detail below.
Our models are GAMMs (generalized additive mixed-effect models) with spatiotemporal Gaussian Markov random fields. The random fields account for spatial and temporal autocorrelation between sampling events, as well as estimating unmeasured components of habitat suitability, allowing that suitability to change through time [36]. The models estimate a spatiotemporal random field that controls for remaining correlated spatial correlation processes each month that are not accounted for by the fixed effects. This random field follows a stationary autoregressive (AR1) process with a first-order correlation (see Supporting information). To allow for a smooth relationship between some predictors and the response variable, we used thin plate regression splines [51] with fixed basis dimensions as calculated via the ‘mgcv’ R package [44, 45]. In all cases, sampling year was included as a factor (estimating a separate mean per year), and numeric month as a spline, to account for seasonal fluctuations in biomass. Our biomass data contained zero and continuous positive values, therefore we used a Tweedie observation model [52, 53] in all models with a log link (Supporting information).
We chose covariates that represent certain features of the coastline and the in situ environment and are known to influence euphausiid distribution [16, 25, 32]. Depth and slope of the seafloor highlight the complex bathymetry of the region, including the three sea valleys of Queen Charlotte Sound and the Juan de Fuca Canyon system, and shallower regions on the Dixon Shelf (Fig 1). Distance to the coast and distance to the shelf edge (1000m isobar) combine to describe the location of land and the continental slope, representing the boundary between shelf waters and the open Pacific Ocean, as well as various sea canyons along the outer edge of the shelf (see S3A–S3D Fig). Distance to the 200m isobar and rugosity were also tested but omitted due to collinearity issues. These datasets were produced by the Marine Spatial Ecology Section, Fisheries and Oceans Canada, Nanaimo, BC [54]. Mean monthly sea surface temperature (SST, °C) [55] and mean monthly surface Chlorophyll A (mg/m3) [56] were obtained from online NOAA databases for April-September (averaged across 2012 to 2019) and were included to represent changes in the physical environment that occur through the spring transition and summer period when upwelling is strongest (see S3E and S3F Fig).
Global models, including all covariates (see Table 1) were used to highlight species-specific responses to environmental variables. For the larval model only, we included the local biomass of E. pacifica and T. spinifera adults as predictors as log(x +1). We standardized covariates by subtracting their mean and scaling by their standard deviation. This helped avoid computational issues when fitting the models and made the magnitude of the coefficients comparable. Depth and Chlorophyll A were log-transformed prior to this process to decrease the importance of large and small biomass values. Each variable was fit in the model as a spline. The number of basis functions (k) for each spline was selected based on the 95% confidence interval of the coefficient estimates that did not overlap zero and inspection of the model residual plots. Model fit was checked via residual plots. Covariates were included as fixed effects in all models, while separate spatial random fields were estimated for each month to allow seasonal (spring and summer) changes to be approximated. Final model formulas can be found in Table 1 and observed vs. predicted biomass plots can be viewed in S5 Fig. The shape of the relationship between each covariate and euphausiid biomass was further assessed to confirm that all relationships appeared biologically reasonable via conditional effects plots. We calculated conditional effect predictions and standard errors for each covariate while fixing the other covariates at the average value of conditions within the occupied depth range of each species (an average weighted by biomass found at each sampling location rather than the overall mean of all sample locations that is represented by zero for scaled variables). This made the illustrated conditional effects representative of abundance under realistic environmental conditions for each species.
Identification and seasonality of biomass hotspots.
To investigate seasonal variability in the spatial distribution of E. pacifica, T. spinifera and euphausiid larvae on the BC coast, we generated maps of mean predicted biomass (mg/m2) for each of the six months. We did not predict for April in the northern bioregion (see Fig 1) as we had no observed data points in this region. Predictions of the biomass of E. pacifica, T. spinifera, and euphausiid larvae were based on all fixed and random effects for each month along a 3km2 grid. We also calculated the spatial ratio of E:T biomass (E. pacifica / T. spinifera), to highlight regions where one species dominated and regions that were preferred habitat for both species.
A major aim of this work was to assess if hotspots of the three euphausiid groups modelled were persistent through the important spring and summer period. To accomplish this, we first determined the spatial coherence/patchiness of predicted euphausiid biomass by employing the Getis-Ord statistic [Gi; 57]. Gi is a statistical Z-score of local clustering (or spatial intensity) relative to the background spatial mean and standard deviation, and therefore is a quantitative measure of hotspot distribution and intensity [see Refs 13, 58, 59 for practical application of the method]. Moran’s I test for spatial autocorrelation was used to inform the spatial neighborhoods (d, distance matrix defining the distance between cells to be considered as ‘neighbours’), which was set at 10km. Separate Gi tests were run for each of the six months for each of the three groups of modelled euphausiids. Z-scores were calculated for the same 3km2 grid as model predictions and mapped. Moran’s I tests and Gi tests were run using ‘spdep’ [60] and ‘usdm’ [61] R packages, respectively. We also combined biomass of the two adult species (‘total adults’) to map combined hotspots of E. pacifica and T. spinifera, which make up the majority of biomass on the BC coast.
To characterize the distribution and intensity of biomass hotspots for each euphausiid group, we mapped the 90th and 95th (top 10% and top 5%) of Z-scores. Percentiles were calculated across the 6-month period for each species and larvae to allow comparisons between spring and summer periods (i.e., to highlight temporal peaks in hotspots for each group). To categorise the persistence of hotspots, we counted the number of times a cell fell in the 95th percentile across months, giving each hotspot cell a score of 1–6. Cells that were hotspots through each month had a perfect score of 6 and were considered persistent. Finally, we calculated the area (km2) of hotspots, and persistent hotspots, for each of the oceanographic ecosections in Fig 1.
Results
Model fit and selection
The inclusion of depth as a predictor reduced the standard deviation of the spatiotemporal random fields (σϵ; see S1 File for equations) in all models (Table 2). For all three groups, conditional-effect plots indicated that the shape of this depth-biomass relationship was roughly quadratic (Fig 2A). The addition of more geomorphic and dynamic covariates reduced the magnitude of this spatiotemporal random field standard deviation σϵ in two out of three models (full models—Table 2), with the largest reduction for euphausiid larvae. σϵ compared with the depth only model for T. spinifera (Table 2); however, the AIC and residual plots (not shown) indicated that the full model was a more parsimonious model fit. The larval model showed a substantial decrease in AIC between the depth-only model and the full covariate model (ΔAIC 130, Table 2). For larvae, the additional inclusion of adult distribution predictors led to an improvement in capturing spatiotemporal variability and therefore we included these in the final model. Based on AIC and inspection of residuals, the rest of the results will focus on analysing the full covariate models for all euphausiids.
Conditional effects of each of the fixed effects (x-axis) included in the full models, illustrating differences in peak biomass density (y-axis) between species (top row) and for total larval biomass (bottom row). The height of each curve represents the predicted density in the month with the highest estimated biomass densities for each group, at the weighted mean occupied conditions of all other covariates, and at the mean of spatiotemporal modelled effects. The shapes of these curves were not allowed to interact with other covariates. Lines represent means and shaded ribbons represent 50% CIs, with the most lightly shaded area representing 95% CIs.
Conditional-effect plot responses were not scaled to represent relative influence of predictors for each model, and therefore represent actual predicted biomass of the euphausiid groups. Fig 2 indicated a tendency for E. pacifica biomass to be highest in areas where the edge of the continental slope described here as the 1000m isobar coincided with the coast (see Fig 1 and S3C and S3D Fig). This environment is found in the north of the study region, to the west of Haida Gwaii, where the continental shelf is narrow. E. pacifica biomass was also highest in areas of increased seafloor slope, and during months with high chlorophyll. Conditional-effect plots indicated the relationship between T. spinifera biomass and the majority of covariates tested fit a quadratic shape. Relationships were weak but indicated the highest biomass of this species occurred nearer to the coast, further from the shelf, and in flatter areas than E. pacifica. Euphausiid larvae biomass was positively correlated with both adult biomass; however, this relationship was much stronger for T. spinifera than E. pacifica biomass (see S4 Fig for conditional effects plot). Larval biomass was also significantly linked to monthly means of the dynamic variables tested; there was a positive relationship between larval biomass and both monthly chlorophyll and SST (Fig 2).
Predicted spatiotemporal biomass variability and persistence of euphausiid hotspots
Seasonal biomass fluctuations were well represented by each model for the different groups. Mean predicted biomass of adults was similar, and low throughout spring, until August and September when T. spinifera exhibited a large peak in biomass (Fig 3). T. spinifera biomass increased across all regions during these months, but was particularly apparent on the Dixon Shelf, the North Coast Fjords, and the Continental slope (Fig 3). Predicated E. pacifica biomass exhibited a smaller and slightly later peak in September, which was mostly confined to the North Coast Fjords and along the Continental slope. This spatial pattern is also apparent in Fig 4A and 4B, with consistently higher biomass of T. spinifera across all regions in the study area, whereas higher E. pacifica biomass was only concentrated along the Continental slope. Overall, euphausiid larvae biomass was smaller than the two adults (Fig 4C) and rose through the spring to peak in June before decreasing slightly (Fig 3). A comparison of the observed vs. predicted biomass from each model can be found in the (S5 Fig).
Mean predicted biomass from each model for a. E. pacifica, b. T spinifera, c. euphausiid larvae averaged across spring and summer (April-September). d. The spatial distribution of the predicted ratio of E. pacifica:T. spinifera biomass. All predictions were made across a 3km2 prediction grid to 1000m depth. Mean calculations did not include predictions for April in the northern bioregion due to gaps in the observed data.
Mapped Z-scores used to calculate hotspots for each group can be found in the, S6A–S6C Fig. All three groups exhibited some seasonal variability in the distribution of hotspots within the 90th percentile across months (Fig 5A–5C); however, larval hotspots exhibited the least persistence (Fig 6A–6C). For larvae, the only persistent hotspot area (classified as present across all 6 months) was found along the west coast of Vancouver Island (or the southern part of the Continental slope ecosection); therefore, the only common persistent hotspot between all three groups was this area. (Fig 6A–6C). For adults, there was high commonality in predicted hotspot distribution, with persistent hotspots along the Continental slope and the fjords of the Northern Coast. Differences between the two species occurred on the Dixon shelf, where T. spinifera exhibited a large and persistent hotspot, whereas E. pacifica did not, and within the Juan de Fuca Eddy system, which was a hotspot region for E. pacifica but never for T. spinifera.
Predicted distribution of Gi Z-scores in the 90th (orange) and 95th (red) percentile for (a) E. pacifica, (b) T. spinifera, and (c) euphausiid larvae. For May to September, gaps off the west coast of Haida Gwaii indicate gaps in predictor surfaces (NAs), we did not predict the distribution of hotshots in the northern bioregion for April due to gaps in the observed data (see S6 Fig).
1 = grid cell counted as a hotspot only in one month, 6 = grid cell persistently ‘hot’ in all 6 months, for (a) E. pacifica, (b) T. spinifera, (c) euphausiid larvae and, (d) total euphausiid adults (E. pacifica + T. spinifera).
Overall adult hotspots displayed higher persistence than larvae, particularly along the along the Continental slope (Fig 6A and 6B and Table 3), with E. pacifica hotspots exhibiting the highest persistence by area. Hotspots in the 95th percentile for E. pacifica and euphausiid larvae encompassed a similar area for all three groups; however, for larvae this area was contained along the west coast of Vancouver Island, with few other hotspot areas in other regions around the coast (Table 3; Fig 6). Despite the overall predicted biomass of larvae being lower than for the two adults (Fig 4), hotspots of this group covered the largest area (Fig 6C and Table 3), indicating that larval biomass was patchier, with higher levels of spatial disparity than adult biomass. The intensity of hotspots (Z-score), or clustered biomass within hotspots, was also highest for larvae (Table 3).
The area (km2) of predicted hotspots showed high spatial variability for all three groups modelled. As well as exhibiting the most persistent hotspots, the Continental slope was also the region with the largest hotspots for all three groups (Table 3). However, there were differences between the northern and the southern slope regions (see Fig 1 for bioregion divisions). Hotspots of both adult species were larger and more persistent along the northern slope, off the coast of Haida Gwaii, and along the edge of Queen Charlotte Sound, while larval hotspots were larger and more persistent in the south, along the continental slope off the west coast of Vancouver Island (Table 3, Fig 6C).
We hypothesized that the Juan de Fuca Eddy system would be a hotspot region for euphausiids (see Fig 1 for location). However, this region was not predicted to be a persistent hotspot area for any of the euphausiids. Hotspots here only occurred for E. pacifica and larvae, and were temporally patchy (Fig 6A and 6C). We also hypothesized that the three sea valleys of Hecate Strait would be important hotspot areas due to increased seafloor complexity. Persistent hotspots occurred for both E. pacifica and T. spinifera along the outer edges of these sea valleys, where the valleys deepened to canyons along the edge of the continental shelf. E. pacifica and T. spinifera hotspots were associated with different valleys however (see ratio of E.T biomass, Fig 4D). Larvae also exhibited hotspots along these sea valleys, but they were unpersistent and occurred in shallower areas on the shelf than adult hotspots. Larval hotspots only occurred in these areas in May, August and September (Fig 5C).
Discussion
We used spatiotemporal modeling of more than 27 years of net-capture data to identify regions of high euphausiid biomass, or ‘hotspots’, for the first time along the Pacific coast of Canada. In addition, few studies have characterized species-specific hotspots of euphausiids. This work is the first step towards quantifying the physical and climate variables that lead to the development and persistence of euphausiid hotspots, with the ultimate goal of delineating important foraging regions for euphausiid vertebrate predators along the BC coast. This type of information is critical, but often lacking, for the development of conservation management strategies such as marine spatial planning [62–64]. However, as with any modelling exercise, some caution should be exercised when interpreting the results. For instance, the modelling framework explained a considerable proportion of variation through spatiotemporal random effects. In addition, as our focus was on seasonal changes, phenological differences between years such as differences in the timing of development of larval biomass, were not considered. Future work should focus on quantifying the effect of these interannual phenological changes, and on improving links between the spatial and temporal scale of response variables with environmental predictors. Both of these factors would likely improve model fit, and may result in tighter coupling between environmental variables such as temperature and the response, which almost certainly shaped euphausiid distribution, but for which we only found a strong correlation with larvae.
Spatiotemporal variability in biomass peaks
In both the observed and modelled data, biomass of E. pacifica and T. spinifera rose throughout the spring to peak in September and August respectively (see S7 Fig). We can infer spawning of both species occurred around the spring transition in April/May, as both E. pacifica and T. spinifera are thought to reach lengths of >10 mm (considered here to be adults) 2–4 months after spawning [2, 12, 29]. In both the northern and central CCE (Oregon and California), reproduction in E. pacifica has been linked to the onset of upwelling (usually April-June) [2, 28], which would also coincide with our observed peaks in euphausiid larvae along the coast in June.
Tanasichuk [12, 29] conducted an in-depth study of the growth and reproduction of E. pacifica and T. spinifera in Barkley Sound, an inshore area along the west coast of Vancouver Island. During 1991–1994, abundances of larger size classes mainly occurred May to September for T. spinifera, and May to July for E. pacifica; however, there was large interannual variability (Tanasichuk 1998a, b [12, 29]). We found that the peaks in biomass that lead to hotspots of both species were patchily distributed, both spatially and across months. These differences increase the potential for species-specific hotspots to occur, effecting the availability of these two dominant species for euphausiid predators.
In many cases, knowledge of the alongshore distribution of euphausiid hotspots comes from acoustic studies, which are unable to differentiate between species. T. spinifera are larger and have a higher lipid content than E. pacifica [65], and some predators preferentially select this species during chick rearing. Cassin’s auklets in the CCE have been found to feed on E. pacifica during pre-breeding and the early part of the breeding season, and then to switch to more energy rich T. spinifera during chick-rearing [1]. In addition, chick growth rates were positively correlated with the proportion of T. spinifera in auklet diets [11]. This suggests that species-specific hotspots of euphausiids may be important for breeding success in seabird predators. Whales can also be selective foragers; a recent study in southern California found that blue whales primarily selected T. spinifera, even when other euphausiid species were present [66].
The spatiotemporal distribution, intensity and persistence of hotspots
In agreement with studies further south in the CCE, our models indicate that euphausiid hotspots were orientated alongshore in areas associated with the shelf break and specific depth contours [14–16]. This study adds to the large body of literature that suggests that a combination of depth, and the edges that occur along topographic features such as the continental slope and seafloor valleys and canyons, are unequivocally important to the development of areas of higher euphausiid biomass [15, 25, 67–69]. We found similarities between the overall distribution of adults, and hotspots of both species were highest over the continental slope, particularly off the west coast of Haida Gwaii. Larger differences may have been observed between the distribution of these two species in waters deeper than 1000m [14], however, unfortunately sampling points in deeper waters were too sparse to be included in our models.
We identified a persistent hotspot region for both adult species and larvae along the shelf break on the west coast of Vancouver Island, which was the only common persistent region between groups through the spring and summer months. This region includes the highly productive La Perouse bank region (south-west Vancouver Island) and is well studied, being one of the most productive fishing grounds off western North America [42]. The southern coast of Vancouver Island is part of the upwelling California Current Ecosystem (CCE), and is highly productive following the spring transition, particularly within the Juan de Fuca Eddy system. Previous studies found that the eddy system was important in generating the biomass of euphausiids required by key euphausiid predators during the summer, particularly Pacific Hake [70, 71]. The northern end of Vancouver Island is part of the down-welling dominated Alaska current, and is considered less productive than the southern section [72]. However, this work highlights the importance of this northern region to euphausiid hotspot development, with persistent hotspots for both adults found here (Fig 6B and 6C).
Mackas et al. [31] suggested a model for euphausiid transport in this region related to local current patterns, euphausiid diel migration, and upwelling. They hypothesized that euphausiid biomass peaks at the shelf break occurred just below, or just inside, upwelling currents converging upwards along the continental slope, optimizing their food supply. Euphausiids avoid advection during the day and night using the surface California Current and subsurface California Current Undercurrent. They do this by positioning at the inshore margin of the shelf break current where equatorward currents are weak. Therefore, at night, they minimize surface (and offshore) advection and are transported onshore and slightly poleward during the day. In this way they maintain the same general position in the water column. This was supported by the results of Lu et al. [3] in the same area, and by Swartzman et al. [73], further south in the CCE. Due to these differences in surface and deep current directions, and the breadth of the continental slope, high production on the shelf and at the shelf-break from local productivity and during large scale upwelling events is extended over a large latitudinal range. We hypothesize that this mechanism may be partly responsible for the high persistence of euphausiid hotspots along the shelf.
A similar process may be responsible for hotspots on the shelf of Queen Charlotte Sound. The shelf encompasses three sea-troughs of ~300 m depth, and multiple sub-marine canyons at the shelf edge [74]. Sea valleys and canyons disrupt alongshore current flow, leading to areas of retention at daytime distributional depths of euphausiids [15, 16]. In addition, these troughs can act as conduits for oceanic waters during upwelling [75], so are often areas of high productivity [27, 76]. They present an ideal habitat for euphausiid larval development [15, 25], and were predicted areas of larval hotspots in this study. In upwelling systems, an optimum window for euphausiid hotspot development has been proposed that balances advective upwelling processes with retention [16, 77]. Canyons are considered regions of high importance to the development of euphausiid hotspots further south in the CCE [15], and were persistent hotspot regions for both adults.
Larval hotspots covered more area, and exhibited higher intensity than adult hotspots, despite overall biomass being lower, indicating tight spatial clustering of biomass. The persistence of larval hotspots was generally low apart from the west coast of Vancouver Island. This is consistent with previous observed patterns for this region [3] and in other regions [78]. Aggregations of euphausiid larval would be expected to be more dynamic than adult hotspots; larvae are less motile than adults and therefore, do not perform extended vertical diel migration. During summer upwelling on the west coast of Vancouver Island, larvae are advected seaward and equatorward in the top 10-30m of the water column as they develop, until the depth range of their diel vertical migration becomes deep enough that they are less exposed to surface seaward-currents [31]. Larval biomass is therefore likely to be more patchily distributed across the coast than adult biomass, perhaps with the majority of biomass deposited at the seaward edge of upwelling cells. However, there is large variability associated with upwelling current patterns, causing differences in the distribution of larval biomass between years, and perhaps leading to the predicted larger spatial area of hotspots averaged across these years. Lu et al. [3], compared larval distribution between two years along the west coast of Vancouver Island, and found differences in the distribution of larvae with respect to the coast between years due to an El Niño event. Further south off Oregon, the same El Niño event caused an increase in the coastal abundance of larvae due to strong onshore advection, ensuring larvae were retained in production centers near to the coast [78].
The Juan de Fuca Eddy and canyon system (see Fig 1), is an area of seasonally heightened primary productivity compared with surrounding water, and was therefore expected to be an area of importance for euphausiids [15, 32]. Hotspots in the 95th percentile occurred in the eddy for E. pacifica and larvae from June to August, however, and never for T. spinifera. Eddy systems enhance mixing and productivity in surface waters, and entrain less motile larvae and other plankton within the eddy structure, leading to higher larval and adult survival [79]. Predation pressure can account for the temporally patchy occurrence of euphausiid hotspots in the eddy during mid-summer. Robinson et al. [80] attributed the rapid decline in euphausiid biomass over the mid-summer period to predation by Pacific Hake. This commercial fish predator migrates onto the southern BC continental shelf in June and are heavy consumers of euphausiids [81]. However, the lack of high biomass of T. spinifera requires further investigation.
Areas of topographic complexity, which are important to euphausiids, are also hotspots for euphausiid predators [6, 21, 26, 82, 83]. The persistence of euphausiid hotspots will aid in the predictability of important foraging regions for predators, which is critical information for management strategies such as marine spatial planning. In an ecological sense, if hotspots of prey are persistent through high and low prey availability, as hotspots along the Continental slope are predicted to be, these regions will be important for ecosystem resilience through extreme events such as El Niño and marine heatwaves. Future work could focus on characterising the interannual persistence of these hotspots, and other productive euphausiid regions on the BC coast such as Dixon shelf and the west coast of Haida Gwaii. Given use of these areas for foraging by species of conservation concern—including Eulachon, Pacific Hake, Pacific Herring, salmon species, and seabirds [21, 31, 84–86]—this would contribute important information for ecosystem-based management of these coastal regions.
Supporting information
S1 Fig. SPDE mesh used in euphausiid modelling for the British Columbia coast, shelf and slope.
Red dots represent knot locations (n = 200) and open black circles represent the locations of euphausiid net-hauls. Triangles represent the SPDE mesh.
https://doi.org/10.1371/journal.pone.0249818.s002
(PDF)
S2 Fig. Distribution of euphausiid sampling stations as modelled.
https://doi.org/10.1371/journal.pone.0249818.s003
(PDF)
S3 Fig.
a-f. Model covariates included in all full models. a-d show static, geomorphic covariates (Depth, Slope, Distance to coast (km), Distance to shelf or 1000m isobar (km). e and f show the dynamic (time-varying) variables: monthly mean of chlorophyll and sea surface temperature (SST).
https://doi.org/10.1371/journal.pone.0249818.s004
(TIF)
S4 Fig. Conditional effects of adult distribution (fixed effects, x-axis) on larval biomass (y-axis).
See Fig 2 for full caption.
https://doi.org/10.1371/journal.pone.0249818.s005
(TIF)
S5 Fig.
Observed vs. predicted plots for (a) E. pacifica, (b) T. spinifera and, (c) euphausiid larvae from the full covariate model specified in Table 1. Vertical lines of dots represent values of 0 at an arbitrary location on the left edge of the y-axis.
https://doi.org/10.1371/journal.pone.0249818.s006
(TIF)
S6 Fig.
a-c. Spatiotemporal variability in Z-score from Getis Ord hotspot analysis for April to September for (a) E pacifica, (b) T. spinifera, and (c) Euphausiid larvae. Predictions were not made for April due to gaps in the observed data.
https://doi.org/10.1371/journal.pone.0249818.s007
(TIF)
S7 Fig.
Seasonal comparison in mean observed vs. mean model estimates of biomass for each of the three models for E. pacifica, T. spinifera, and euphausiid larvae (left to right).
https://doi.org/10.1371/journal.pone.0249818.s008
(TIF)
Acknowledgments
We would like to acknowledge all personnel from the Fisheries and Oceans Canada and the Canadian Coast Guard who collected this data from all over the BC coast through multiple sampling programs, aboard many ships including the CCGS John P. Tully and the CCGS W.E. Ricker. We especially thank Moira Galbraith and David L. Mackas who have collected much of this data, processed and identified the samples, maintained the database, and published much of the earlier zooplankton research from this region.
References
- 1. Abraham CL, Sydeman WJ. Prey-switching by Cassin’s auklet Ptychoramphus aleuticus reveals seasonal climate-related cycles of Euphausia pacifica and Thysanoessa spinifera. Marine Ecology Progress Series. 2006;313:271–83.
- 2. Brinton E. Population biology of Euphausia pacifica off southern California. Fishery Bulletin. 1976;74(4):733–62.
- 3. Lu B, Mackas DL, Moore DF. Cross-shore separation of adult and juvenile euphausiids in a shelf-break alongshore current. Progress in Oceanography. 2003;57(3–4):381–404.
- 4. Emmett RL, Brodeur RD, Miller TW, Pool SS, Bentley PJ, Krutzikowsky GK, et al. Pacific sardine (Sardinops sagax) abundance, distribution, and ecological relationships in the Pacific Northwest. California Cooperative Oceanic Fisheries Investigations Report. 2005;46:122.
- 5. Hipfner J. Euphausiids in the diet of a North Pacific seabird: annual and seasonal variation and the role of ocean climate. Mar Ecol Prog Ser. 2009 Sep 18;390:277–89.
- 6. Nickels CF, Sala LM, Ohman MD. The euphausiid prey field for blue whales around a steep bathymetric feature in the southern California current system. Limnology and Oceanography. 2019;64(1):390–405.
- 7. Schoenherr JR. Blue whales feeding on high concentrations of euphausiids around Monterey Submarine Canyon. Canadian Journal of Zoology. 1991;69(3):583–94.
- 8. Hagen W, Van Vleet E, Kattner G. Seasonal lipid storage as overwintering strategy of Antarctic krill. Mar Ecol Prog Ser. 1996;134:85–9.
- 9. Smith SE, Adams PB. Daytime surface swarms of Thysanoessa spinifera (Euphausiacea) in the Gulf of the Farallones, California. Bulletin of Marine Science. 1988;42(1):76–84.
- 10. Brinton E, Townsend A. Decadal variability in abundances of the dominant euphausiid species in southern sectors of the California Current. Deep Sea Research Part II: Topical Studies in Oceanography. 2003;50(14–16):2449–72.
- 11. Abraham CL, Sydeman WJ. Ocean climate, euphausiids and auklet nesting: inter-annual trends and variation in phenology, diet and growth of a planktivorous seabird, Ptychoramphus aleuticus. Marine Ecology Progress Series. 2004;274:235–50.
- 12. Tanasichuk R. Interannual variations in the population biology and productivity of Thysanoessa spinifera in Barkley Sound, Canada, with special reference to the 1992 and 1993 warm ocean years. Mar Ecol Prog Ser. 1998;173:181–95.
- 13. Dorman JG, Sydeman WJ, García-Reyes M, Zeno RA, Santora JA. Modeling krill aggregations in the central-northern California Current. Marine Ecology Progress Series. 2015;528:87–99.
- 14. Cimino MA, Santora JA, Schroeder I, Sydeman W, Jacox MG, Hazen EL, et al. Essential krill species habitat resolved by seasonal upwelling and ocean circulation models within the large marine ecosystem of the California Current System. Ecography [Internet]. 2020 [cited 2020 Sep 24];n/a(n/a). Available from: https://onlinelibrary.wiley.com/doi/abs/10.1111/ecog.05204 pmid:33304029
- 15. Santora JA, Zeno R, Dorman JG, Sydeman WJ. Submarine canyons represent an essential habitat network for krill hotspots in a Large Marine Ecosystem. Scientific Reports. 2018 May 15;8(1):7579. pmid:29765085
- 16. Santora JA, Sydeman WJ, Schroeder ID, Wells BK, Field JC. Mesoscale structure and oceanographic determinants of krill hotspots in the California Current: Implications for trophic transfer and conservation. Progress in Oceanography. 2011;91(4):397–409.
- 17. Santora JA, Schroeder ID, Field JC, Wells BK, Sydeman WJ. Spatio‐temporal dynamics of ocean conditions and forage taxa reveal regional structuring of seabird–prey relationships. Ecological Applications. 2014;24(7):1730–47. pmid:29210234
- 18. Fiechter J, Santora JA, Chavez F, Northcott D, Messié M. Krill hotspot formation and phenology in the California Current Ecosystem. Geophysical research letters. 2020;47(13):e2020GL088039. pmid:32728303
- 19. Manugian S, Elliott ML, Bradley R, Howar J, Karnovsky N, Saenz B, et al. Spatial Distribution and Temporal Patterns of Cassin’s Auklet Foraging and Their Euphausiid Prey in a Variable Ocean Environment. PLOS ONE. 2015 Dec 2;10(12):e0144232. pmid:26629818
- 20. Sydeman WJ, Thompson SA, Santora JA, Koslow JA, Goericke R, Ohman MD. Climate–ecosystem change off southern California: time-dependent seabird predator–prey numerical responses. Deep Sea Research Part II: Topical Studies in Oceanography. 2015;112:158–70.
- 21. Santora JA, Field JC, Schroeder ID, Sakuma KM, Wells BK, Sydeman WJ. Spatial ecology of krill, micronekton and top predators in the central California Current: Implications for defining ecologically important areas. Progress in Oceanography. 2012 Nov;106:154–74.
- 22. Santora JA, Ralston S, Sydeman WJ. Spatial organization of krill and seabirds in the central California Current. ICES Journal of Marine Science. 2011;68(7):1391–402.
- 23. Rockwood RC, Elliott ML, Saenz B, Nur N, Jahncke J. Modeling predator and prey hotspots: Management implications of baleen whale co-occurrence with krill in Central California. PloS one. 2020;15(7):e0235603. pmid:32634142
- 24. Benoit-Bird KJ, Battaile BC, Heppell SA, Hoover B, Irons D, Jones N, et al. Prey Patch Patterns Predict Habitat Use by Top Marine Predators with Diverse Foraging Strategies. PLOS ONE. 2013 Jan 3;8(1):e53348. pmid:23301063
- 25. Allen SE, Vindeirinho C, Thomson RE, Foreman MG, Mackas DL. Physical and biological processes over a submarine canyon during an upwelling event. Can J Fish Aquat Sci. 2001 Apr 1;58(4):671–84.
- 26. Burger AE. Effects of the Juan de Fuca Eddy and upwelling on densities and distributions of seabirds off southwest Vancouver Island, British Columbia. Marine Ornithology. 2003;31(2):113–22.
- 27. Croll D, Marinovic B, Benson S, Chavez F, Black N, Ternullo R, et al. From wind to whales: trophic links in a coastal upwelling system. Mar Ecol Prog Ser. 2005;289:117–30.
- 28. Feinberg LR, Peterson WT. Variability in duration and intensity of euphausiid spawning off central Oregon, 1996–2001. Progress in Oceanography. 2003;57(3–4):363–79.
- 29. Tanasichuk R. Interannual variations in the population biology and productivity of Euphausia pacifica in Barkley Sound, Canada, with special reference to the 1992 and 1993 warm ocean years. Mar Ecol Prog Ser. 1998;173:163–80.
- 30. Hsieh WW, Ware DM, Thomson RE. Wind-induced upwelling along the west coast of North America, 1899–1988. Canadian Journal of Fisheries and Aquatic Sciences. 1995;52(2):325–34.
- 31. Mackas DL, Kieser R, Saunders M, Yelland DR, Brown RM, Moore DF. Aggregation of euphausiids and Pacific hake (Merluccius productus) along the outer continental shelf off Vancouver Island. Can J Fish Aquat Sci. 1997 Sep 1;54(9):2080–96.
- 32. Simard Y, Mackas DL. Mesoscale Aggregations of Euphausiid Sound Scattering Layers on the Continental Shelf of Vancouver Island. Can J Fish Aquat Sci. 1989 Jul 1;46(7):1238–49.
- 33. Mackas DL, Thomson RE, Galbraith M. Changes in the zooplankton community of the British Columbia continental margin, 1985–1999, and their covariation with oceanographic conditions. Canadian Journal of Fisheries and Aquatic Sciences. 2001;58(4):685–702.
- 34. Burger AE, Hitchcock CL, Davoren GK. Spatial aggregations of seabirds and their prey on the continental shelf off SW Vancouver Island. Marine Ecology Progress Series. 2004;283:279–92.
- 35. Mackas DL, Peterson WT, Zamon JE. Comparisons of interannual biomass anomalies of zooplankton communities along the continental margins of British Columbia and Oregon. Deep Sea Research Part II: Topical Studies in Oceanography. 2004;51(6–9):875–96.
- 36. Shelton AO, Thorson JT, Ward EJ, Feist BE. Spatial semiparametric models improve estimates of species abundance and distribution. Can J Fish Aquat Sci. 2014 Jul 8;71(11):1655–66.
- 37. Thorson JT, Shelton AO, Ward EJ, Skaug HJ. Geostatistical delta-generalized linear mixed models improve precision for estimated abundance indices for West Coast groundfishes. ICES Journal of Marine Science. 2015;72(5):1297–310.
- 38. Rubidge E, Jeffery S, Gregr EJ, Gale KS, Frid A. Assessment of nearshore features in the Northern Shelf Bioregion against criteria for determining Ecologically and Biologically Significant Areas (EBSAs). DFO Can. Sci. Advis. Sec. Res. Doc; 2020.
- 39. Harper JR, Christian J, Cross WE, Frith R, Searing GF, Thompson D. A classification of the marine regions of Canada. Environment Canada, Ottawa. 1993;
- 40. Mackas DL. Interannual variability of the zooplankton community off southern Vancouver Island. Canadian Special Publication of Fisheries and Aquatic Sciences. 1995;603–15.
- 41. Ohman MD, Lavaniegos BE. Comparative zooplankton sampling efficiency of a ring net and bongo net with comments on pooling of subsamples. California Cooperative Oceanic Fisheries Investigations Report. 2002;162–73.
- 42. McFarlane GA, Ware DM, Thomson RE, Mackas DL, Robinson CLK. Physical, biological and fisheries oceanography of a large ecosystem (west coast of Vancouver Island) and implications for management. Oceanologica Acta. 1997;20(1):191–200.
- 43. Lindgren F, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2011;73(4):423–98.
- 44. Anderson SC, Keppel EA, Edwards A. A reproducible data synopsis for over 100 species of British Columbia groundfish. DFO Canadian Science Advisory Secretariat Research Document. 2019;41(2019):1–330.
- 45.
Anderson SC, Ward EJ, Barnett LAK, English PA. sdmTMB: Spatiotemporal Species Distribution GLMMs with “TMB” [Internet]. Available from: https://pbs-assess.github.io/sdmTMB/index.html
- 46. Rue H, Martino S, Lindgren F, Simpson D, Riebler A, Krainski ET. INLA: functions which allow to perform a full Bayesian analysis of structured additive models using Integrated Nested Laplace Approximation. R package version 00. 2009.
- 47. Kristensen K, Nielsen A, Berg CW, Skaug H, Bell BM. TMB: Automatic Differentiation and Laplace Approximation. Journal of Statistical Software. 2016 Apr 4;70(1):1–21.
- 48. Godefroid M, Boldt JL, Thorson JT, Forrest R, Gauthier S, Flostrand L, et al. Spatio-temporal models provide new insights on the biotic and abiotic drivers shaping Pacific Herring (Clupea pallasi) distribution. Progress in Oceanography. 2019;178:102198.
- 49. Marshall KN, Duffy-Anderson JT, Ward EJ, Anderson SC, Hunsicker ME, Williams BC. Long-term trends in ichthyoplankton assemblage structure, biodiversity, and synchrony in the Gulf of Alaska and their relationships to climate. Progress in Oceanography. 2019 Jan;170:134–45.
- 50. Quiroz ZC, Prates MO. Bayesian spatio-temporal modelling of anchovy abundance through the SPDE Approach. Spatial Statistics. 2018 Dec 1;28:236–56.
- 51. Wood SN. Thin plate regression splines. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2003;65(1):95–114.
- 52. Dunn PK, Smyth GK. Series evaluation of Tweedie exponential dispersion model densities. Statistics and Computing. 2005;15(4):267–80.
- 53.
Tweedie MC. An index which distinguishes between some important exponential families. In: Statistics: Applications and new directions: Proc Indian statistical institute golden Jubilee International conference. 1984. p. 579–604.
- 54. Fields C, Nephin J. Government of Canada; Fisheries and Oceans Canada; MSEA Section, 2020. Environmental Layers—Northern & Southern Shelf Bioregions (100 m). 2020.
- 55. JPL MUR MEaSUREs P. GHRSST Level 4 MUR Global Foundation Sea Surface Temperature Analysis [Internet]. 2015 [cited 2020 Dec 15]. Available from: https://doi.org/10.5067/GHGMR-4FJ04
- 56. NASA. NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group. Visible and Infrared Imager/Radiometer Suite (VIIRS) [Internet]. 2021 [cited 2020 Dec 15]. Available from:
- 57. Getis A, Ord JK. The Analysis of Spatial Association by Use of Distance Statistics. Geographical Analysis. 1992;24(3):189–206.
- 58. Kuletz KJ, Ferguson MC, Hurley B, Gall AE, Labunski EA, Morgan TC. Seasonal spatial patterns in seabird and marine mammal distribution in the eastern Chukchi and western Beaufort seas: Identifying biologically important pelagic areas. Progress in Oceanography. 2015;136:175–200.
- 59. Santora JA, Reiss CS, Loeb VJ, Veit RR. Spatial association between hotspots of baleen whales and demographic patterns of Antarctic krill Euphausia superba suggests size-dependent predation. Marine Ecology Progress Series. 2010;405:255–69.
- 60. Bivand R, Wong DWS. Comparing implementations of global and local indicators of spatial association. TEST. 2018;27(3):716–48.
- 61. Naimi B, Hamm N a s, Groen TA, Skidmore AK, Toxopeus AG. Where is positional uncertainty a problem for species distribution modelling. Ecography. 2014;37:191–203.
- 62. Bertram DF, Mackas DL, Welch DW, Boyd WS, Ryder JL, Galbraith M, et al. Variation in zooplankton prey distribution determines marine foraging distributions of breeding Cassin’s Auklet. Deep Sea Research Part I: Oceanographic Research Papers. 2017 Nov 1;129:32–40.
- 63. Montevecchi WA, Hedd A, Tranquilla LM, Fifield DA, Burke CM, Regular PM, et al. Tracking seabirds to identify ecologically important and high risk marine areas in the western North Atlantic. Biological Conservation. 2012;156:62–71.
- 64. Sherley RB, Botha P, Underhill LG, Ryan PG, van Zyl D, Cockcroft AC, et al. Defining ecologically relevant scales for spatial protection with long‐term data on an endangered seabird and local prey availability. Conservation Biology. 2017;31(6):1312–21. pmid:28248436
- 65. Fisher JL, Menkel J, Copeman L, Shaw CT, Feinberg LR, Peterson WT. Comparison of condition metrics and lipid content between Euphausia pacifica and Thysanoessa spinifera in the northern California Current, USA. Progress in Oceanography. 2020;102417.
- 66. Nickels CF, Sala LM, Ohman MD. The morphology of euphausiid mandibles used to assess selective predation by blue whales in the southern sector of the California Current System. Journal of Crustacean Biology. 2018;38(5):563–73.
- 67. Macquart-Moulin C, Patriti G. Accumulation of migratory micronekton crustaceans over the upper slope and submarine canyons of the northwestern Mediterranean. Deep Sea Research Part I: Oceanographic Research Papers. 1996;43(5):579–601.
- 68. Levine R. Temporal And Spatial Variability In Euphausiid Abundance, Biomass, And Species Composition At The Northwest Atlantic Shelf Break And Its Canyons. 2014;
- 69. Plourde S, Lehoux C, McQuinn IH, Lesage V. Describing krill distribution in the western North Atlantic using statistical habitat models. Canadian Science Advisory Secretariat; 2017.
- 70. Robinson CL, Ware DM. Modelling pelagic fish and plankton trophodynamics off southwestern Vancouver Island, British Columbia. Canadian Journal of Fisheries and Aquatic Sciences. 1994;51(8):1737–51.
- 71. Ware DM, McFarlane GA. Climateinduced changes in Pacific hake (Merluccius productus) abundance in the Vancouver Island upwelling system. Climate change and northern fish populations Canadian Special Publication of Fisheries and Aquat Sci. 1995;121:509–21.
- 72. Mackas D. Zooplankton on the west coast of Vancouver Island: distribution and availability to marine birds. The ecology, status, and conservation of marine and shoreline birds on the west coast of Vancouver Island. 1992.
- 73. Swartzman G, Hickey B, Kosro PM, Wilson C. Poleward and equatorward currents in the Pacific Eastern Boundary Current in summer 1995 and 1998 and their relationship to the distribution of euphausiids. Deep Sea Research Part II: Topical Studies in Oceanography. 2005;52(1–2):73–88.
- 74. Shaw J, Barrie JV, Conway KW, Lintern DG, Kung R. Glaciation of the northern British Columbia continental shelf: the geomorphic evidence derived from multibeam bathymetric data. Boreas. 2020;49(1):17–37.
- 75.
Dosser HV, Jackson J, Waterman S, Hunt B, Hannah CG. Connecting Small-Scale Processes to Water-Mass Transformation and Nutrient Availability in the Coastal Pacific. In: 2018 Ocean Sciences Meeting. AGU; 2018.
- 76. Rennie SJ, Pattiaratchi CB, McCauley RD. Numerical simulation of the circulation within the Perth Submarine Canyon, Western Australia. Continental Shelf Research. 2009 Sep 15;29(16):2020–36.
- 77. Cury P, Roy C. Optimal environmental window and pelagic fish recruitment success in upwelling areas. Canadian Journal of Fisheries and Aquatic Sciences. 1989;46(4):670–80.
- 78.
Peterson WT, Feinberg L, Kiester J. Ecological Zonation of euphausiids off central Oregon. In: Report of the 1999 Monitor and Rex Workshops and the 2000 Model workshop on lower trophic level modeling PICES, Sidney. Citeseer; 2000. p. 125–8.
- 79. Goldstein ED, Pirtle JL, Duffy‐Anderson JT, Stockhausen WT, Zimmermann M, Wilson MT, et al. Eddy retention and seafloor terrain facilitate cross‐shelf transport and delivery of fish larvae to suitable nursery habitats. Limnology and Oceanography. 2020;
- 80. Robinson CLK. The influence of ocean climate on coastal plankton and fish production. Fisheries Oceanography. 1994;3(3):159–71.
- 81. Tanasichuk RW, Ware DM, Shaw W, McFarlane GA. Variations in diet, daily ration, and feeding periodicity of Pacific hake (Merluccius productus) and spiny dogfish (Squalus acanthias) off the lower west coast of Vancouver Island. Canadian Journal of Fisheries and Aquatic Sciences. 1991;48(11):2118–28.
- 82. Flinn RD, Trites AW, Gregr EJ, Perry RI. Diets of fin, sei, and sperm whales in British Columbia: an analysis of commercial whaling records, 1963–1967. Marine Mammal Science. 2002;18(3):663–79.
- 83. Santora JA, Reiss CS, Cossio AM, Veit RR. Interannual spatial variability of krill (Euphausia superba) influences seabird foraging behavior near Elephant Island, Antarctica. Fisheries Oceanography. 2009;18(1):20–35.
- 84. Alton MS, Nelson MO. Food of Pacific hake, Merluccius productus, in Washington and northern Oregon coastal waters. US Fish and Wildlife Service Circular. 1970;332:35–42.
- 85. Vermeer K, Fulton JD, Sealy SG. Differential use of zooplankton prey by Ancient murrelets and Cassin’s auklets in the Queen Charlotte Islands. J Plankton Res. 1985;7(4):443–59.
- 86. Wells BK, Santora JA, Field JC, MacFarlane RB, Marinovic BB, Sydeman WJ. Population dynamics of Chinook salmon Oncorhynchus tshawytscha relative to prey availability in the central California coastal region. Marine Ecology Progress Series. 2012;457:125–37.