Three decades of nearshore surveys reveal long-term patterns in gray whale habitat use, distribution, and abundance in the Northern California Current

The nearshore waters of the Northern California Current support an important seasonal foraging ground for Pacific Coast Feeding Group (PCFG) gray whales. We examine gray whale distribution, habitat use, and abundance over 31 years (1992–2022) using standardized nearshore (< 5 km from shore) surveys spanning a large swath of the PCFG foraging range. Specifically, we generated density surface models, which incorporate detection probability into generalized additive models to assess environmental correlates of gray whale distribution and predict abundance over time. We illustrate the importance of coastal upwelling dynamics, whereby increased upwelling only yields higher gray whale density if interspersed with relaxation events, likely because this combination optimizes influx and retention of nutrients to support recruitment and aggregation of gray whale prey. Several habitat features influence gray whale distribution, including substrate, shelf width, prominent capes, and river estuaries. However, the influence of these features differs between regions, revealing heterogeneity in habitat preferences throughout the PCFG foraging range. Predicted gray whale abundance fluctuated throughout our study period, but without clear directional trends, unlike previous abundance estimates based on mark-recapture models. This study highlights the value of long-term monitoring, shedding light on the impacts of variable environmental conditions on an iconic nearshore marine predator.


Figure S1
. Bioregions of the California Current Large Marine Ecosystem (Schroeder et al. 2022).The survey area for this study is shown in the blue polygon, with study area latitudinal boundaries denoted by the dashed blue lines.

Figure S2
. Schematic of the sampling scheme within a primary sampling unit (PSU).The shoreline is denoted by the dark gray.Survey segments are illustrated as the red lines, in both the inshore and offshore subunits.S1.      (hazard-rate, halfnormal).Models were compared using the difference in Akaike's information criterion (ΔAIC, lower value indicates better relative performance), and the realism of the predicted effective strip width (ESW) under different covariate conditions.In this case, the models that used a hazard-rate key produced unrealistically small ESW values despite having a lower AIC value, necessitating revised selection based on ecological knowledge of the system, as can be the case in cetacean distance sampling studies (Williams and Thomas 2007).Therefore, we elected to proceed with the best performing model that used a half-normal key, which included BSS as a covariate.Corresponding plots of the detection function fit to the distance data are shown for all candidate detection functions in Fig. S5.

Figure S3 .
Figure S3.Static environmental predictors included in the gray whale density surface models: A) bathymetric depth, B) distance to the coast, C) shelf width, measured as the distance between the coast and the 200 m isobath, illustrated by the white line, D) distance to the nearest prominent cape, with negative values indicating distances north of the nearest cape, positive values indicating distances south of the nearest cape, and capes shown by the white triangles, E) distance to the nearest river estuary, with river estuaries show by the white crosses, and F) benthic substrate type.

Figure S4 .
Figure S4.A) Map of the study area, with the region 1 and region 2 boundaries shown by the dashed lines, and major placenames denoted.B) Heatmap illustrating the survey effort throughout the study period, where the y-axis corresponds to the latitude on the map, the x-axis corresponds to the study year, and the fill color indicates the survey effort in km.

Figure S5 .
Figure S5.Candidate detection models fit with either a hazard rate (left column) or half-normal (left column) key, including no covariates (first row), Beaufort sea state (BSS, second row), sightability (third row), or both BSS and sightability (categorical classifier ranging from very bad to excellent, fourth row).The truncation distance has been set to the 98 th percentile of the data, which is 1,000 m.Vertical bars show a histogram of the perpendicular distances between the gray whale observations and the survey trackline, and the black line shows the candidate detection function fit to the data, with detection probability shown on the y-axis.The selected detection function for the DSMs was fit with a half-normal key and BSS as a covariate.Corresponding comparisons of the candidate detection function performance are shown in TableS1.

Figure S6 .
Figure S6.Diagnostic plots for GAM model outputs: normal Q-Q plots and histograms of residuals.The left column is for region 1, and the right column is for the region 2 model.

Figure S7 .
Figure S7.Daily predicted abundance estimates across the study period each year, for region 1 (top panel) and region 2 (bottom panel).Gray lines each represent one year, and the black line represents the mean of all years (1992-2022).

Figure S8 .
Figure S8.Timeseries cross-correlation results between predicted annual gray whale abundance and either the Pacific Decadal Oscillation (PDO) or Multivariate ENSO Index (MEI), in the spring (left column) or summer (right column), for each region.

Figure S9 .
Figure S9.Relationship between predicted annual gray whale abundance and the variability in PDO and MEI, calculated over the preceding 1, 2, 3, and 5 years.

Table S1 .
Comparison of candidate detection function models, fitted with different covariates (none, Beaufort sea state, overall sightability) and different key functions

Table S2 .
Results of linear regression models evaluating the relationship between annual predicted gray whale abundance and variability in ocean basin-scale indices.Variability was calculated as the standard deviation of PDO and MEI, computed at four different timeframes: 1, 2, 3, and 5 years prior to the year of interest.Relationships were assessed separately for each region.Bold p-values with an asterisk (*) indicate a statistically significant linear relationship.