Environmentally Driven Seasonal Forecasts of Pacific Hake Distribution

Changing ecosystem conditions present a challenge for the monitoring and management of living marine resources, where decisions often require lead-times of weeks to months. Consistent improvement in the skill of regional ocean models to predict physical ocean states at seasonal time scales provides opportunities to forecast biological responses to changing ecosystem conditions that impact fishery management practices. In this study, we used 8-month lead-time predictions of temperature at 250 m depth from the J-SCOPE regional ocean model, along with stationary habitat conditions (e.g., distance to shelf break), to forecast Pacific hake (Merluccius productus) distribution in the northern California Current Ecosystem (CCE). Using retrospective skill assessments, we found strong agreement between hake distribution forecasts and historical observations. The top performing models [based on out-of-sample skill assessments using the area-under-the-curve (AUC) skill metric] were a generalized additive model (GAM) that included shelf-break distance (i.e., distance to the 200 m isobath) (AUC = 0.813) and a boosted regression tree (BRT) that included temperature at 250 m depth and shelf-break distance (AUC = 0.830). An ensemble forecast of the top performing GAM and BRT models only improved out-of-sample forecast skill slightly (AUC = 0.838) due to strongly correlated forecast errors between models (r = 0.88). Collectively, our results demonstrate that seasonal lead-time ocean predictions have predictive skill for important ecological processes in the northern CCE and can be used to provide early detection of impending distribution shifts of ecologically and economically important marine species.


INTRODUCTION
Anticipating ecological change is an important component of living marine resource management where decisions often require lead-times of weeks to months. Yet, the lack of advanced warnings about the response of marine taxa to ecosystem shifts limits the ability of management systems to respond to rapidly changing ecosystem conditions (Clark et al., 2001;Dietze et al., 2018). Increasingly, seasonal ecological forecasts provide a means to reduce related uncertainties and play a key role in supporting management of living marine resources into the future (Hobday et al., 2016;Payne et al., 2017;Tommasi et al., 2017). Indeed, over the past decade seasonal ecological forecasts have been developed for a wide range of marine taxa including American lobster in the Gulf of Maine (Mills et al., 2017), sardines in the California Current (Zwolinski et al., 2011;Kaplan et al., 2016), and southern bluefin tuna in eastern Australia (Hobday et al., 2011a(Hobday et al., , 2016. Increases in the predictive skill of physical ocean states has partially driven the increased availability of seasonal ecological forecasts and has resulted in the availability of skillful ocean forecasts with seasonal lead-times for many of the world's large marine ecosystems (Stock et al., 2015;Tommasi et al., 2017;Jacox et al., 2020). In the northern California Current Ecosystem (CCE), the J-SCOPE (JISAO's Seasonal Coastal Ocean Prediction of the Ecosystem) model provides forecasts of physical, chemical, and biological ocean states with seasonal lead times (e.g., 6-9 months) . Skill assessments have shown the J-SCOPE model has considerable predictive skill at seasonal lead times for several ecologically relevant variables including subsurface temperature . In turn, J-SCOPE seasonal forecasts of ocean conditions can then be used to drive ecological forecasts, such as sardine distribution in the CCE .
Pacific hake (Merluccius productus, hereafter just hake) is an important mid-trophic-level species in the CCE that supports one of the largest United States groundfish fisheries outside of Alaska (Ressler et al., 2007;Berger et al., 2019). Hake are distributed from about 25 • to 55 • N and at depths typically between 100 and 1000 m. The dynamics of the hake stock are characterized by episodic recruitment events with a few large ageclasses dominating the stock (Hamel et al., 2015;Berger et al., 2019). Age-structure of the stock, in turn, influences distribution since older and larger hake tend to be distributed further north than smaller and younger conspecifics (Berger et al., 2019). Hake growth is variable across years and is at least partly influenced by ocean conditions (e.g., El Niño events) and availability of prey resources (Ressler et al., 2007;Hamel et al., 2015).
Pacific hake are seasonally migratory, with a northward spring migration from southern spawning grounds off the United States west coast, terminating as far north as southeast Alaska. This migration pattern results in hake being a trans-boundary resource fished commercially in the United States and Canada (Bailey et al., 1982). The fraction of the population that migrates into Canadian waters, however, can vary greatly across years, creating challenges for monitoring and management planning (Dorn, 1995). For instance, monitoring of the hake stock is conducted jointly by a United States/Canada summer acoustic-trawl survey that provides an index of hake biomass that is used for stock assessment and management planning (Berger et al., 2019). The ability of the monitoring survey to sample the full spatial extent of the stock partially determines the magnitude of uncertainty associated with the biomass index.
Environmental conditions influence the summer distribution of hake along the west coast of North America (Benson et al., 2002;Ressler et al., 2007;Agostini et al., 2008). Thermal conditions, in particular, have been positively associated with the fraction of the hake stock in Canadian waters, suggesting warmer ocean conditions drive a more northern distribution of hake (Dorn, 1995;Ware and McFarlane, 1995). More recently, evidence has suggested that thermal conditions have a spatially variable effect on hake distribution with strong positive associations with hake biomass north of Vancouver Island, British Columbia (BC) and strong negative associations offshore of Vancouver Island, BC and Washington, United States (Malick et al., 2020). This suggests that ocean temperatures could be a useful predictor of hake distribution in the northern CCE.
Skillful forecasts of hake distribution could help inform management and survey planning decisions in three important aspects. First, early warnings of changes in hake distribution can inform planning of fisheries independent surveys used to monitor the hake stock (Payne et al., 2017). For example, survey planning decisions, such as allocating survey effort between northern and southern areas, are made several months prior to the start of the survey. Thus, forecasts could inform decisions about allocating limited survey effort by predicting areas where hake are unlikely to be present in a given year. If vessel breakdowns or weather forced a reduction in survey effort, transect density could be reduced in regions predicted to have low probability of hake occurrence. Second, skillful forecasts provide information on the projected trans-boundary distribution of hake, and thus could help reduce uncertainties in the availability of the hake stock to fishers in Canada and the United States (Hobday et al., 2011b;Mills et al., 2017). Third, skillful forecasts provide early warnings of potential ecosystem shifts that can inform ecosystem-based management (Levin et al., 2009;Malick et al., 2017). For instance, Pacific hake are an important predator of fish and shellfish populations and are prey for larger fish and marine mammals in the CCE, thus advanced warnings of shifts in hake distribution could aid detection of consequential ecological shifts in the CCE (Bailey et al., 1982;Francis, 1983).
In this study, we examined whether seasonal forecasts of physical oceanographic conditions can be used to accurately predict hake distribution in the northern CCE. In particular, we developed and tested 8-month lead-time forecasts of summer hake distribution with the goal of providing forecasts to support management and survey planning decisions. We used 7 years of acoustic-trawl survey data to characterize hake distribution. We then used the J-SCOPE regional ocean model to develop 8month lead-time forecasts of subsurface temperatures that were used to force environmentally driven species distribution models for hake. We further evaluated whether multi-model ensembles improved forecast skill of hake distribution by comparing ensemble forecasts to single-model forecasts. This process of using oceanographic forecasts to predict hake distributional shifts in the CCE explicitly addresses fisheries management ecosystem-linkage goals and provides necessary context for short-term oceanographic variability within the scope of longerterm perturbations (e.g., climate change).

MATERIALS AND METHODS
The primary motivation for developing seasonal forecasts of hake distribution was to provide an early warning of changes in hake distribution to support management decisions and Pacific hake acoustic-trawl survey planning. As a result, forecasts were co-developed with survey planners, while stakeholder involvement occurred via meetings associated with the Pacific hake treaty process.

Pacific Hake Data
The hake survey aims to sample the full range of the hake distribution in summer, and survey extent and the number of transects are often adjusted in response to the presence or absence of hake following survey design guidelines. Therefore, we focus on forecasting the probability of hake occurrence, rather than density, because the acoustic-trawl survey is better informed by early warnings in the expansion or contraction of hake distribution than forecasts of hake density in a given location.
We used 7 years of spatially explicit biennial hake occurrence data collected via joint United States/Canada acoustic-trawl surveys from 2009 to 2019, with an additional 2012 survey ( Table 1). Surveys started in southern California and moved northward along the United States and Canada west coasts until hake were no longer observed (typically around 54.5 • N). The spatial extent of data analyzed here, however, was limited to the region 43-50 • N -the latitudinal domain of the J-SCOPE model used to generate forecasts of ocean conditions (see below). The number of annual survey transects within the study domain ranged from 34 in 2015 to 49 in 2011 ( Table 1). Survey timing was fairly consistent across years, with the southern third of the study domain typically sampled during the second half of July and the northern two-thirds typically sampled during August.
Acoustic backscatter measurements attributable to hake were converted to hake biomass using the procedures outlined in Fleischer et al. (2008) and Malick et al. (2020). We aggregated the N transects gives the annual number of survey transects. N bins gives the annual number of 10 km bins across the study region. % absent gives the annual percentage of 10 km bins with no Pacific hake. hake data into 10 km bins to reduce spatial autocorrelation in the data and coded bins with non-zero biomass as hake occurrences. We also tested smaller (e.g., 5 km) and larger (e.g., 20 km) bin sizes and found results were robust across different sized bins.

Ocean Forecasts
We used 8-month lead-time forecasts of temperature at 250 m depth from J-SCOPE to forecast hake distribution (Supplementary Figure S1). J-SCOPE is a Regional Ocean Modeling System (ROMS) (Haidvogel et al., 2008) simulation of seasonal ocean conditions spanning 43-50 • N on the outer coast of Washington, Oregon, and southern BC . The J-SCOPE model has a 1.5 km horizontal resolution with 40 vertical levels and includes both rivers and tides. The large scale oceanic and atmospheric forcing comes from NOAA's global Climate Forecast System (CFS). In this study, we focus on retrospective ocean forecasts, i.e., reforecasts, which are true forecasts for a historical period using a free-running model unconstrained by observations after initialization. The aim in using these reforecasts was to test the models skill for jointly predicting future ocean conditions and hake distribution 8months ahead.
We chose temperature at 250 m depth as our primary ocean variable because (1) previous research has shown strong correlations between temperature at depth and hake distribution (Malick et al., 2020), and (2) 250 m represents depths commonly occupied by hake (Ressler et al., 2007). In areas where bottom depth was less than 250 m, we used bottom temperature instead. July and August temperature forecasts were generated for each survey year using a January initialization period. For 2019, three model runs from CFS were used to quantify the uncertainty related to those forcing variables. The model runs were chosen from the beginning (January 5), middle (January 15), and end (January 25) of the forecast initialization month. The initial conditions for J-SCOPE ROMS consist of the average conditions from CFS-reanalysis for the initialization month of the forecast. As is typical in the oceanographic literature, we focus on anomalies -i.e., differences from the climatology or time-averaged field detailing the seasonal cycle. In this case, the J-SCOPE reforecast climatology was based on 2009-2017, building on Siedlecki et al. (2016).
In addition to the dynamic temperature variable, we also explored a static index of cross-shelf location as a predictor of hake distribution. In particular, we used the distance to the 200 m shelf break, where the distance was defined as the minimum euclidean distance between a hake observation and the 200 m isobath. Positive values of the shelf distance variable indicated the hake observation was offshore of the 200 m isobath and negative values indicated the observation was inshore.

Statistical Forecasting Models
We used both generalized additive models (GAM) and boosted regression trees (BRT) to model species distribution, because previous studies have shown the potential for differences in explanatory power and predictive skill across model types (Abrahms et al., 2019;Brodie et al., 2019).
We used binomial GAMs with a logit link to predict the probability of hake occurrence, where y is the logit transformed probability of occurrence for year t at location i, α is the intercept, s is a univariate smooth function of shelf break distance (S), f is a smooth function that describes the effect of temperature anomaly at 250 m depth (T; i.e., deviations from the long-term mean), and t,i are model residuals.
We considered two alternative formulations for the f temperature term ( Table 2). The first formulation assumed a spatially stationary temperature effect modeled as a univariate smooth function of temperature, s T t,i . The second formulation allowed for a spatially variable temperature effect by modeling the temperature effect as the product of T t,i and a bivariate smooth function of longitude (lon) and latitude (lat), i.e., g lon, lat · T t,i . Non-parametric thin plate regression splines were used for the univariate (s) and bivariate (g) smooth functions in the GAMs (Wood, 2003).
Two simpler GAMs that included static covariates were considered as alternative null forecast models ( Table 2). The first simpler model included a univariate smooth of shelf break distance. The second simpler model included a bivariate smooth of longitude and latitude. In total, we evaluated four alternative GAM forecast models ( Table 2).
The BRT used a Bernoulli distribution. Since BRT models can handle colinearity among predictors, we included four covariates: temperature anomaly at 250 m depth, distance to the shelf break, longitude, and latitude (Elith et al., 2008). BRT models are composed of a large number of decision trees constructed via recursive binary splits of the data with non-linear responses produced by evaluating these splits across many trees (Elith et al., 2008). The BRT was estimated using a maximum of three interactions among covariates, a learning rate of 0.02, and bag fraction of 0.6, which resulted in models with at least 1000 trees (Elith et al., 2008). Standard errors of predicted species distribution were calculated across 100 BRT fits to provide model error estimates (Hazen et al., 2018;Brodie et al., 2019).
Four ensemble model forecasts were generated by averaging each of the GAMs and the BRT, where each model in the ensemble was given equal weight (Clemen, 1989;Araujo and New, 2007). We also tested the sensitivity of our results to the inclusion of temperature in the BRT by re-running the analysis with temperature anomaly at 250 m excluded from the BRT model. All analyses were conducted using R v3.6.0 and the mgcv and dismo packages (Elith et al., 2008;Wood, 2017;R Core Team, 2019).

Forecast Evaluation
Previous work has shown considerable skill for the J-SCOPE temperature forecast model. Siedlecki et al. (2016) evaluated J-SCOPE's predictability for temperature, and found that predictive skill increased with depth. In addition, annual evaluation of J-SCOPE forecasts against observations are available on the NANOOS IOOS portal 1 . Supplementing earlier evaluations of J-SCOPE performance and skill, we further quantified performance of J-SCOPE predictions of bottom temperature by comparing against temperature data from the Northwest Fisheries Science Center's West Coast Groundfish Bottom Trawl Survey (Keller et al., 2017). That survey samples the United States West Coast slope and shelf (55-1280 m) annually from May-October, targeting bottom-dwelling species of commercial importance.
We evaluated the GAM and BRT model performance using a combination of in-sample and out-of-sample metrics including the area-under-the-curve (AUC) and mean squared error (MSE) (Fielding and Bell, 1997). The AUC measures how well a model can discriminate a presence from an absence. The AUC ranges between 0 and 1 where a value of 0.5 indicates a random classifier and values closer to 1 indicate higher forecast skill. The MSE measures the accuracy of the forecast model where lower values indicate a more accurate forecast.
We used leave-one-year-out cross validation to evaluate how the models performed at forecasting hake distribution for unobserved years (Fielding and Bell, 1997). In this procedure, a single year was left-out of the data, each model was re-fit using the remaining years of data, and forecasts were produced for the 1 http://www.nanoos.org/products/j-scope/home.php To further test the performance of the forecast models, we generated "true" out-of-sample forecasts for 2019 prior to the hake acoustic-trawl survey. We fit the hake forecast models using data through 2017 and then tested those models using true forecasts of the physical environment from the J-SCOPE January initialized temperature forecasts. True forecasts are produced every January and released on the web in February prior to the conditions being observed, thus referred to as a "true" forecast. In this case, J-SCOPE ocean forecasts were used to generate forecasts of hake distribution for August 2019. To better characterize uncertainty in the ocean forecasts, we generated three forecasts from the J-SCOPE model for 2019 using different initialization dates in January (January 5, January 15, January 25). The spread across the three forecasts was used as a measure of uncertainty in the J-SCOPE temperature forecasts. For each hake forecasting model, we generated separate forecasts for each of the three temperature forecasts and used an average across the three forecasts as our ensemble mean 2019 hake forecast for each model. Performance of the 2019 hake forecasts was evaluated using the AUC and MSE metrics.

RESULTS
J-SCOPE performed well when bottom temperatures observed in situ were compared with the simulated bottom temperatures from the same locations (R 2 = 0.88; Supplementary Figure S2). The predicted bottom temperatures were biased warm (RMSE = 0.48), which is not uncommon in ROMS applications and is addressed here by focusing on anomalies rather than raw temperature values (Giddings et al., 2014). This strong agreement between observed and predicted temperatures supports the use of numerical ocean model forecasts of sub-surface temperatures to predict suitable hake habitat.
Hake occurred across the latitudinal extent of the study region with the exception of 2011, when no fish were observed off  (Figure 1). Across all years, the cross-shelf distribution of hake was concentrated around the 200 m isobath with the majority of hake occurrences (62%) occurring within 10 km of the shelf break (Figure 1). The percentage of hake absences across years was consistent ranging from 64% in 2019 to 77% in 2013, although the latitudinal distribution of absences varied across years ( Table 1). For example, in 2015 hake were absent from most locations off-shore of Washington and Oregon, whereas in 2011 hake occurrences tended to be concentrated in this region (Figure 1).
The shelf-break term in the GAM and BRT models confirmed the strong preference of hake to be present slightly offshore of the 200 m shelf-break (Figures 2A, 3A-C). The shelfbreak preference was also the dominant pattern in the bivariate smooth of longitude and latitude in model GAM1 (Figure 2B).
The stationary temperature terms in the GAM3 and BRT models indicated that hake occurrence tended to have a positive association with temperature anomaly (Figures 2C, 3D). In contrast, the non-stationary temperature term in model GAM4 showed negative associations between temperature and hake occurrence off the Washington and Oregon coasts, but positive associations in more northern and southern areas ( Figure 2D).
All forecasting models had considerable forecast skill (both insample and out-of-sample) with AUC values greater than 0.79 and MSE values lower than 0.17 ( Table 2). The BRT tended to fit the data the best with the highest in-sample AUC (0.93) and lowest in-sample MSE (0.09). For out-of-sample, however, an ensemble model (ENS2) performed best with an overall AUC of 0.84 and MSE of 0.14 ( Table 2). Among the four GAMs, the longitude-latitude model (GAM1) tended to fit the data the best  (i.e., had the highest AUC and lowest MSE), whereas the shelfonly model (GAM2) had the best out-of-sample forecast skill ( Table 2). Between the two GAMs that included a temperature effect (GAM3 and GAM4), there was support for a spatially varying temperature effect with the GAM4 model having better in-sample and out-of-sample performance than GAM2.
The 2019 temperature anomaly forecasts from the J-SCOPE model indicated above average temperatures at depth across the study region ( Figure 4A) with the warmest forecast ( Figure 4C) having an average temperature anomaly of 0.78 • C and the coolest forecast ( Figure 4B) having an average anomaly of 0.36 • C. The three individual temperature forecasts displayed moderate variability in temperature anomalies with grid cell specific standard deviations across the three temperature forecasts ranging from 0.03 to 1.36 (Supplementary Figure S3).
All models had considerable skill in forecasting 2019 hake occurrence with the ENS1 model having the best 2019 forecast skill with an AUC of 0.85 and MSE of 0.16 (Table 2 and Figure 5). The 2019 forecasts showed higher probabilities of hake occurrence near the 200 m isobath, which is consistent with their historical distribution within the study region. The three models that included temperature anomaly (GAM3, GAM4, and BRT) showed more spatial variability in predicted hake occurrence and also tended to have higher standard errors of prediction compared to models that lacked temperature (Figures 5, 6).
When temperature was removed from the BRT, model fit declined compared to the original BRT model that included temperature (i.e., lower in-sample skill; Table 2 and Supplementary Table S1). In addition, the model with the highest out-of-sample skill changed to an ensemble of the BRT without temperature and GAM4, which includes a non-stationary temperature effect, suggesting temperature contributes to out-of-sample forecast skill.

DISCUSSION
Our objective was to develop and test environmentally driven seasonal forecasts of hake distribution to support management and survey planning decisions. The forecast models we tested showed appreciable out-of-sample forecast skill at 8-month time horizons. In addition, we found that: (1) the J-SCOPE model had considerable predictive skill of subsurface temperatures throughout the study domain, (2) distance to the 200 m shelf break was a strong predictor of historical hake occurrence and temperature at depth had a spatially varying effect on the probability of occurrence; and (3) the BRT model had moderately higher forecast skill than the GAMs and a multi-model ensemble forecast had slightly better out-ofsample forecast skill compared to the individual GAM and BRT models. Together, our results suggest that comparatively simple models can forecast hake distribution using seasonal projections of subsurface ocean temperature and distance to the 200 m shelf break. Most of the forecast skill was derived from the 200 m shelfbreak covariate. The strong affinity for hake occurrences to be concentrated in areas just-offshore of the 200 m shelf break is consistent with several previous studies that have shown areas near the 200 m isobath and areas with steeply sloping bathymetry provide good hake habitat (Dorn, 1995;Mackas et al., 1997;Swartzman, 1997Swartzman, , 2001. One possible explanation for this is that food availability may be high in these areas. In particular, euphausiids -an important prey item for haketend to concentrate in areas of steeply sloping bathymetry and submarine canyons (Buckley and Livingston, 1997;Mackas et al., 1997;Santora et al., 2018). In addition, areas just offshore of the 200 m isobath may also provide good physical ocean conditions. For instance, the California Undercurrent is strongest offshore of the 200 m isobath; the northward flowing undercurrent may act as a migration corridor for hake that could facilitate northward migration of Pacific hake and aggregate prey resources (Bakun, 1996;Agostini et al., 2006).
Our results indicated a moderate subsurface temperature effect on hake occurrence. The best GAMs included a spatially variable temperature effect and the BRT indicated the temperature term accounted for ∼17% of the variability in the response. This result broadly agrees with several previous papers that have shown associations between hake and temperature at depth (Ressler et al., 2007;Hamel et al., 2015;Malick et al., 2020). Temperature most likely acts as a proxy for other processes that have a more direct impact on hake distribution because the temperature ranges analyzed here are comparable to previously observed in situ temperature preferences of hake (Bailey et al., 1982;Ressler et al., 2007). Although using variables that have a more direct impact on hake habitat preferences (e.g., food availability) may provide better forecasts, skillful forecast of lower-trophic-level processes relevant for hake (e.g., euphausiid distribution) are currently not available. In contrast, temperature provides an ecologically relevant variable for which there is forecast skill at the lead-times important for decision makers Siedlecki et al., 2016;Jacox et al., 2017).
The nine forecasting models evaluated here (five individual models + four ensemble models) performed similarly across years when forecasting out-of-sample hake distribution, e.g., most forecasts for 2011 and 2017 had relatively low skill, whereas forecasts for 2009 and 2012 had relatively high skill (Supplementary Figures S4, S5). Two factors likely contributed to lower forecast skill in some years. First, gaps in the latitudinal distribution of hake reduced skill, which occurred in 2011 when hake stock size was lower and few were present off the west coast of Vancouver Island. Second, variability in the cross-shelf distribution also appears to reduce skill; in 2017, hake occurrences were concentrated just inshore of the 200 m isobath, but just offshore in all other years (Figure 1 and Supplementary Figure S6). This suggests that the consistency in which hake are present just offshore of the 200 m isobath across the latitudinal range of the study area drives differences in forecast skill among years. A priority for future research would be to examine additional covariates to better capture inter-annual deviations in distribution from the shelf-break, e.g., the California Undercurrent or subsurface oxygen concentrations.
Combining multiple forecasts into an ensemble forecast has been widely shown to produce increased forecast precision compared to individual model forecasts, given that the individual forecasts provide some independent information (Bates and Granger, 1969;Clemen, 1989;Abrahms et al., 2019). We found that the ensemble hake forecasts had only slightly better out-ofsample skill compared to the individual model forecasts ( Table 2). The weak increase in predictive performance for the ensemble models compared to the individual GAM and BRT models is likely due to high correlations among model prediction errors. Multi-model ensemble models tend to outperform individual model predictions when weakly or negatively correlated model predictions are combined due to cancelation of random errors (Clemen, 1989). In this study, however, the individual model forecast error (i.e., observed hake occurrence -forecasted probability of occurrence), were strongly correlated, e.g., the average pairwise correlation of forecast errors from the GAMs and BRT was 0.91 (Supplementary Figure S7). This indicates that the individual forecast models produce similar forecast errors, which reduces the effectiveness of multi-model averaging (Araujo and New, 2007).
The results presented here provide a critical first step in developing an early warning of hake distributional shifts. Yet, we believe future work on three areas could further improve the usefulness of hake forecasts for management and survey planning. First, extending the northern range of this work to include waters through SE Alaska could help inform survey planning by providing additional information on the projected northern extent of Pacific hake, which is a critical uncertainty during survey planning. Second, developing a forecast of hake density could improve how the forecasts inform management decisions by helping to reduce uncertainties regarding the proportion of the population expected to migrate into Canadian waters. Third, if the spatial extent of the study area is extended northward beyond 50 • N, the maximum latitudinal domain of this study, accounting for impacts of age-structure on Pacific hake distribution may be important. Exploratory analysis did not identify strong age-based differences in Pacific hake occurrences within the current spatial extent; however, evidence suggests that older and larger hake tend to migrate further north than smaller hake with few age-2 hake observed north of Vancouver Island (Ressler et al., 2007;Malick et al., 2020).
Collectively, our results provide evidence that hake distribution can be skillfully forecast at lead-times of 8-months in the northern CCE. Our results also illustrate the broader utility of using seasonal lead-time ocean predictions in an ecological context to provide early warnings of distribution shifts of ecologically and economically important marine species. Marine ecosystems are changing rapidly and experiencing extreme events more frequently. Thus, skillful ecological forecasts provide new tools to inform the management process by reducing uncertainties regarding future states of nature that management decisions are often dependent upon (Dietze et al., 2018).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
MM was the principal author of the manuscript and led the modeling analysis. MM, MH, SP-S, SS, IK, MH, KM, and AB designed the study and contributed to interpreting the results. SP-S and SG performed the data collection and processing. SS, EM, AH, and NB ran the J-SCOPE model forecasts and conducted the skill analysis. All authors contributed to writing and editing of the manuscript.

FUNDING
Funding was provided by NOAA's Fisheries and the Environment (FATE) program (grant #16-08). This study was partially supported by NOAA's Climate Program Office's Modeling, Analysis, Predictions, and Projections Program, Grant #NA17OAR4310112. This publication is partially funded by the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under NOAA Cooperative Agreement NA15OAR4320063, Contribution No. 2018-0183. The J-SCOPE forecasts were performed on the University of Washington Hyak supercomputer system, supported in part by the University of Washington eScience Institute.