A full description of the methods used to develop these layers is available in Thompson et al. (in review).

Survey Data - The species distribution models are based on data from three fisheries-independent scientific surveys conducted within Canadian Pacific Waters: the Fisheries and Oceans Canada (DFO) Groundfish Synoptic Bottom Trawl Surveys (Trawl; Sinclair et al. 2003; Anderson et al. 2019), the DFO Groundfish Hardbottom Longline Surveys (HBLL; Lochead and Yamanaka 2006, 2007; Doherty et al. 2019), and the International Pacific Halibut Commission Fisheries Independent Setline Survey (IPHC; IPHC 2021). In all cases, we based our analyses on data from 2003 to 2020 and utilized presence-absence records to facilitate integrating across surveys and gear types. Our analysis included all species which had at least 150 observations (to exclude species that are rarely sampled by the surveys) across all the survey datasets, resulting in a final set of 66 species.

Environmental data - As environmental predictors in the models, we used log seafloor depth, bathymetric position index (BPI), substrate indices for muddiness and rockiness, log tidal speed (m s^(-1)), mean summer ocean current speed (m s^(-1)), mean summer salinity (practical salinity units - PSU), and mean summer salinity range (PSU). In our species distribution model (described below), we included all environmental predictors other than depth and BPI as linear fixed effects. We included depth and BPI as second-order-polynomials to allow species occurrences to peak at mid values. We also included year as a penalized spline (i.e., generalized additive model; Wood 2017) to account for year-to-year variation in species’ occurrences.

Species distribution model - We fit species distribution models using Generalized Linear Mixed Effects Models with the sdmTMB package (Anderson et al. 2021). This model included the environment variables described above as fixed effects as well spatiotemporal random effects to account for variation in species occurrences that is not associated with the environmental predictors. In models that integrated multiple surveys (see below) we included survey type (i.e., trawl, HBLL, IPHC) as a fixed effect to account for differences in catchability across surveys. See Thompson et al. in review for a full specification and details on the model.

Candidate models - For each species, we compared the predictive power of three classes of models. These were: 1) models fit on the individual surveys separately, 2) models that integrated surveys that shared the same gear type (integrated gear model), and 3) a model that integrated both gear types and all surveys (fully integrated model). We excluded surveys that included fewer than 10 occurrences for the focal species because including these surveys resulted in poor model convergence. Thus, for some species, the fully integrated model or one or both (i.e., trawl and longline) of the integrated gear models were not run. The Trawl survey was conducted in the Strait of Georgia in 2012 and 2015 and is typically excluded in analyses of the broader dataset (Anderson et al. 2019) because of the differences in gear and the limited temporal coverage. Thus, we considered the SoG Trawl survey as a separate survey in our single survey models.

All models included the fixed and spatiotemporal random effects described previously except for the SoG Trawl survey on its own where fixed effects other than depth were dropped because more complex models failed to converge due to the low number of observations in that dataset. We integrated the surveys using a categorical fixed effect (i.e., HBLL vs. IPHC vs. Trawl) to account for methodological differences between the surveys. Thus, species responses to the environmental conditions were shared across surveys in the integrated models, but the overall probability of occurrence could vary depending on the survey. Although there are some methodological differences in the trawl surveys in the SoG compared to the other regions, we assume that these differences would have minimal effect on presence-absence data. Furthermore, it would not be possible to distinguish these effects from the spatial and environmental differences that exist between the SoG and the outer shelf, which we assume have a larger impact on species occurrences. Thus, we elected to consider all trawl surveys as a single survey in the integrated models.

Model comparison - We compared models based on their predictive accuracy estimated using spatial block cross-validation. We used three-fold cross-validation with folds arranged sequentially over 24 blocks that divided the coastline equally from north to south. We assessed predictive accuracy using the area under the curve (AUC; Pearce and Ferrier 2000) of the withheld data in the cross-validation.

Selection of the final model for each species was made by applying a series of decisions:

1) Exclude SoG Trawl single survey models because of their spatial coverage constraints.

2) Exclude all models with a cross-validation AUC less than 0.7 because of low predictive accuracy.

3) Select the highest level of integration within the remaining models.

4) When a fully integrated model is selected, select the survey for making predictions (i.e., survey fixed effect) based on which survey has the higher proportion of occurrences for that species except when the proportion of occurrences is greater than 98% of all sets. This occurs in the IPHC survey for Pacific Halibut and North Pacific Spiny Dogfish, species which target the longline hooks, thus resulting in extremely high predicted probability of occurrence in a presence-absence model.

Spatial species occurrence predictions - We made predictions of occurrence using the final selected model for each species. These predictions were made at a 1 km resolution across the full British Columbia coast in waters for areas where seafloor depth was between 10 and 1400 meters and mean summer salinity was above 28 PSU. These values correspond to the approximate range of environmental conditions sampled in the surveys. Environmental values were obtained by averaging the 100 m resolution environmental data described previously within each 1 km grid cell. 1 km resolution for our selected models was chosen because it corresponds with the average distance between the start and end coordinates of the HBLL surveys and is in the same magnitude as the roughly 2 km average trawl. Finer scale predictions could be made if necessary (e.g., Nephin et al. 2020). However, this would require using the downscale environmental values and making the assumption that a model trained on trawl data, which integrates multiple kilometers of habitat, can make predictions at finer spatial resolutions. All species distribution predictions were made using models fit to the full data set, rather than using the cross-validation models. Predictions were made for 2012–2015 and then averaged across years in each grid cell. We elected to use this time span because the SoG Trawl survey was only conducted in 2012 and 2015 and so our predictions would not require extrapolation beyond the years in which we have data. Using the mean value across multiple years ensured that our predictions are not heavily influenced by year-to-year variation in predicted species occurrences. Predictions should be interpreted as the probability of catching a species in a given location if sampling using the survey identified as most appropriate in the model selection process.

We estimate model uncertainty by running 500 simulations based on draws from the joint precision matrix of our model, assuming multivariate normal parameter covariance (sensu Gelman and Hill 2007). We estimate the median standard error for each 1 km grid cell by taking the median standard error across all species. To identify species occurrence hotspots, we calculated the proportion of simulated values in each grid cell that exceeded a high occurrence threshold. This threshold was the lower of: 1) 0.8, or 2) the 80th percentile of the predicted probability of occurrence values across all grid cells that had a probability of occurrence greater than 0.05. Excluding grid cells with a probability of occurrence lower than 0.05 ensured that habitats that were unsuitable were not included. Using an upper bound of 0.8 for the hotspot threshold ensured that thresholds were not unreasonably high for species with high occurrence probability where it is predicted to be found (e.g. > 0.999 for Longspine Thornyhead).

Species richness predictions - We predicted species richness by summing the predicted species occurrences within each grid cell across all the selected final models (sensu Ovaskainen and Abrego 2020). This value should be interpreted as the number of species that would be expected to be caught if that location were surveyed using the most appropriate survey method for each species. This assumes, for example, that if one species has a probability of occurrence of 0.6 and another species has a probability of 0.8 in a given location, then a survey would catch an average of 1.4 species that location if sampled repeatedly.