Hotspots for rockfishes, structural corals, and large-bodied sponges along the central coast of Pacific Canada

The Wuikinuxv, Kitasoo/Xai’xais, Heiltsuk and Nuxalk First Nations hold Indigenous rights to their territories, where all data were collected. Scientific staff who are members of these Nations or who work directly for them had direct approvals from Indigenous rights holders and were exempt from other research permit requirements. Collaborating DFO scientists worked in partnership with the First Nations to collect data in their territories..

Sampling targeted rocky reefs, the preferred habitat for most Sebastidae³⁸, which we located through local Indigenous knowledge or a bathymetric model⁴⁹. Data were collected by four fishery-independent methods—shallow diver transects, mid-depth video transects, deep video transects, and hook-and-line sampling—detailed in earlier publications^{32,33,34,35,50,51} and summarized in Table 1. Data had a spatial resolution of ≤ 130 m² and each sampling location (N = 2936 for Sebastidae, 2654 for sponges, 2321 for corals) was ascribed to a 1-km² planning unit within the standardized grid used to design the MPA network (N = 632 for Sebastidae, 525 for sponges, 529 for corals, 516 inclusive of surveys for all taxonomic groups).

Although sampling encompassed 11 years (2006–2007, 2013–2021: Table 1), 84% of 1-km² planning units were sampled during only one year (Appendix S2). Analyses, therefore, focus on spatial variability in species distributions and do not address temporal variability within planning units. When all years and methods are combined, 1-km² planning units had a median of 3 samples (range = 1 to 80, Q₁ = 2, Q₃ = 6) (i.e., sum of dive transects, video sub-transects, and hook-and-line sessions). Supplementary Data Set 1 reports sampling effort by 1-km² planning unit, survey type, and year (see Data Availability for link to these data).

For each 1-km² planning unit, u, we calculated hotspot indices for Sebastidae (B_SEB,u), structural corals (B_Cor,u), and large-bodied sponges (B_Sp,u). These indices did not consider cup corals, whip-like corals or encrusting corals or sponges.

As detailed below (Eqs. 1–4), each species of Sebastidae or genera of corals contributed to B_SEB,u or B_Cor,u, according to their abundance weighted by W_t: a conservation prioritization score based on taxon characteristics. For the 26 species of Sebastidae that we observed, W_t equaled the sum of scores for (1) fishery vulnerability, using intrinsic population growth rate, r, as a proxy variable^52,53, (2) depletion level, using the ratio of recent biomass to unfished biomass as a proxy variable, (3) ecological role, with trophic level as proxy, and (4) evolutionary distinctiveness¹⁴ (Table 2; Appendix S3). Because several rockfishes are very long-lived (i.e., have low values for r) and depleted, maximum potential scores were twice as large for fishery vulnerability and depletion level than for ecological role and evolutionary distinctiveness. Data for depletion level and evolutionary distinctiveness were unavailable for some species, and score calculations (detailed in Table 2) account for missing values (Appendix S3).

For the 6 genera of structural corals analyzed (Appendix S4), W_t depended on mean height (estimated from video transect images: Table 1), which correlates positively with vulnerability to physical damage from bottom-contact fishing gear (including longer time to recovery)^20,54,55 and with strength of ecological role (e.g., amount of biogenic habitat and carbon sequestration increases with height)^44,56 (Table 2, Appendix S4). W_t for corals did not include depletion level due to lack of data.

The hotspot index for large-bodied sponges, B_Sp,u did not differentiate between species characteristics (i.e., ({W}_{t}=1)) and we pooled the abundances of all observed species of Hexactinellidae (Aphrocallistes vastus, Farrea occa, Heterochone calyx, Rhabdocalyptus dawsoni, Staurocalyptus dowlingi) and Demospongiae (Mycale cf loveni). This approach is consistent with regional fishery bodies worldwide, which treat large-bodied sponges as a single functional group⁵⁷.

To derive hotspot indices for each taxonomic group (Sebastidae, structural corals, or large-bodied sponges), we first developed a set of candidate generalized linear mixed models (GLMM) to explain relative abundance data for rockfish, corals, and sponges. For each GLMM, we estimated ({lambda }_{t,i,l}), the expected counts (or expected percent cover) for taxa t obtained with survey method i at point location l. (Point locations are individual dive transects, video transect bins, or hook-and-line timed sessions: Table 1.) Specifically,

where g was the link function for the GLMM and f the distribution for the likelihood function modelling either the observed counts C (negative binomial) for Sebastidae and structural corals or a combination of counts (negative binomial) and percent cover D (beta distribution) for large-bodied sponges. We used multiple GLMMs to model large-bodied sponges because deep video transects recorded actual counts whereas dive or mid-depth video transects recorded percent cover categories (Table 1).

For each taxonomic group, we estimated a set of coefficients (beta) for the vector of X covariates that best estimated counts or percent cover. Our hypothesized covariates included the 1-km² planning unit (modelled as a random intercept to control for repeated measures within a given planning unit), survey method, depth (including both linear or a 2^nd order polynomial), and taxa. Each GLMM controlled for sample effort as an offset—effort was measured either as area covered by dive transects or video bins, or the duration of hook-and-line sessions. We also tested for possible covariate’s effects on the dispersion parameter (for the negative binomial GLMMs) and zero-inflation terms (for both the negative binomial and beta GLMMs). The best set of covariates to predict counts or percent cover were then chosen based on AIC model selection criteria. All models were fitted using ‘glmmTMB’⁵⁸ in R version 4.0.2⁵⁹, and simulated residuals and diagnostic tests performed for each best-fit model using the package ‘DHARMa’⁶⁰. For example, our best model for Sebastidae counts predicted 2% fewer zero counts than were observed.

We applied depth and survey method selectivity criteria to reduce excessive zeroes in the count data that may be biologically unjustified (Appendix S5). For all taxon, if i detected t, then the method was valid for that taxon. If i did not detect t and t is a Sebastidae, then the method was valid (i.e., count = 0) only if the overall 10th and 90th percentiles of depths sampled by that method encompassed the expected depth range of t (Appendix S5). If i did not detect t and t is a coral or sponge (which are rarer than Sebastidae), then the method is valid only if the depth of the sampling event exceeded or equaled the minimum expected depth of t. Also, hook-and-line gear cannot systematically sample sessile benthic organisms or planktivores and this method was valid only for non-planktivorous Sebastidae (Appendix S5).

Using the best-fit models from above, we calculated the expected count (or percent cover) per unit of effort, (mu), for taxa t observed with method i at each planning unit u:

where ({n}_{i,u}) was the total number of point locations sampled by that method within the planning unit and effort was either the cumulative area covered by dive or video surveys or the cumulative duration of hook-and-line sampling sessions within the planning unit. Because survey methods differed in their maximum values and potential biases (e.g., field of view is greater for divers than for video cameras; hook-and-line gear samples one fish at a time while visual methods can observe multiple fish simultaneously),({mu }_{t,i,u}) was rescaled as a min–max normalization,({mu }_{t,i,u}^{^{prime}}) (i.e., difference between the observed value and the minimum value across all u, divided by the range of values across all u).

The hotspot index for each of Sebastidae, structural corals, and large-bodied sponges (denoted as taxonomic group g) was then calculated for each planning unit as:

where W_t was the taxon-specific weighing factor (Table 2, Appendices S3, S4), ({n}_{s,g}) was the number of species in taxonomic group g, and ({n}_{m,g}) was the number of valid methods to sample group g.

For each 1-km² planning unit where all taxonomic groups were surveyed (N = 518), we then calculated the overall hotspot index:

where H is Shannon’s evenness index, with proportional abundance of each taxonomic group represented by B_SEB,u, B_Cor,u, and B_Sp,u.

Hotspot index values were normalized as the proportion of the maximum value and converted to decile ranks. Relationships between decile ranks and index values were nonlinear (Appendix S6).

For Sebastidae, large-bodied sponges, and the overall hotspot index, we defined hotspots as planning units containing decile ranks 9 or 10: criterion which we deemed appropriate for the small spatial scales of conservation planning being used for the central portion of the Northern Shelf Bioregion (16-km² planning units in Fig. 2). We are aware that other studies define hotspots based on a narrower range of values (e.g., top 10%²⁶; top 2.5%²⁸) but their context is generally one in which conservation planning is done at a much greater scale (e.g., ≈50,000-km² grid cells²⁶;1° latitude × 1° longitude grid cells²⁸). For structural corals, which had near-zero index values in all but the top-ranking planning units (Appendix S6), we defined hotspots as planning units containing decile rank 10.

Maximum depths sampled within planning units were deepest in the Mainland Fjord and shallowest in the Aristazabal Banks Upwelling Upper Ocean Subregion (Appendix S7). Accordingly, we used multiple logistic regression implemented with the ‘glm’ function in R to estimate the probabilities hotspot occurrence within 1-km² planning units in relation to maximum depth sampled (including a 2nd-order polynomial) and Upper Ocean Subregion. Competing models were compared with AIC model selection procedures.

Following the directive of Central Coast First Nations, decile rank distributions were mapped as 16-km² planning units, u₁₆ (N = 283 for Sebastidae, 264 for sponges, 263 for corals, 260 inclusive of surveys for all taxonomic groups), thereby protecting sensitive locations that would be revealed at smaller scales. To do so, we took the average between the maximum index value and the mean of the remainder of index values among the 1-km² planning units, u, contained within each u₁₆, and converted these values into decile ranks. This approach balances conservation prioritization among u₁₆ that may have good average index values for multiple u, and u₁₆ with a single high-ranking u among multiple low-scoring u. Relationships between decile ranks and hotspot index values also were nonlinear at this scale (Appendix S6). The same hotspot definitions developed for u apply to u₁₆.

Eighty one percent of 16-km² planning units were sampled during only one or two years (Appendix S2). When all years and methods are combined, 16-km² planning units had a median of 6 samples (range = 1 to 110, Q₁ = 3, Q₃ = 13). Supplementary Data Set 2 reports sampling effort by 16-km² planning unit, survey type, and year (see Data Availability for link to these data).

Source: Ecology - nature.com