The Importance of Environmental Exposure History in Forecasting Dungeness Crab Megalopae Occurrence Using J-SCOPE, a High-Resolution Model for the US Pacific Northwest

The Dungeness crab (Metacarcinus magister) fishery is one of the highest value fisheries in the US Pacific Northwest, but its catch size fluctuates widely across years. Although the underlying causes of this wide variability are not well understood, the abundance of M. magister megalopae has been linked to recruitment into the adult fishery 4 years later. These pelagic megalopae are exposed to a range of ocean conditions during their dispersal period, which may drive their occurrence patterns. Environmental exposure history has been found to be important for some pelagic organisms, so we hypothesized that inclusion of recent environmental exposure history would improve our ability to predict inter-annual variability in M. magister megalopae occurrence patterns compared to using “in situ” conditions alone. We combined 8 years of local observations of M. magister megalopae and regional simulations of ocean conditions to model megalopae occurrence using a generalized linear model (GLM) framework. The modeled ocean conditions were extracted from JISAO’s Seasonal Coastal Ocean Prediction of the Ecosystem (J-SCOPE), a high-resolution coupled physical-biogeochemical model. The analysis included variables from J-SCOPE identified in the literature as important for larval crab occurrence: temperature, salinity, dissolved oxygen concentration, nitrate concentration, phytoplankton concentration, pH, aragonite, and calcite saturation state. GLMs were developed with either in situ ocean conditions or environmental exposure histories generated using particle tracking experiments. We found that inclusion of exposure history improved the ability of the GLMs to predict megalopae occurrence 98% of the time. Of the six swimming behaviors used to simulate megalopae dispersal, five behaviors generated GLMs with superior fits to the observations, so a biological ensemble of these models was constructed. When the biological ensemble was used for forecasting, the model showed skill in predicting megalopae occurrence (AUC = 0.94). Our results highlight the importance of including exposure history in larval occurrence modeling and help provide a method for predicting pelagic megalopae occurrence. This work is a step toward developing a forecast product to support management of the fishery.

The Dungeness crab (Metacarcinus magister) fishery is one of the highest value fisheries in the US Pacific Northwest, but its catch size fluctuates widely across years. Although the underlying causes of this wide variability are not well understood, the abundance of M. magister megalopae has been linked to recruitment into the adult fishery 4 years later. These pelagic megalopae are exposed to a range of ocean conditions during their dispersal period, which may drive their occurrence patterns. Environmental exposure history has been found to be important for some pelagic organisms, so we hypothesized that inclusion of recent environmental exposure history would improve our ability to predict inter-annual variability in M. magister megalopae occurrence patterns compared to using "in situ" conditions alone. We combined 8 years of local observations of M. magister megalopae and regional simulations of ocean conditions to model megalopae occurrence using a generalized linear model (GLM) framework. The modeled ocean conditions were extracted from JISAO's Seasonal Coastal Ocean Prediction of the Ecosystem (J-SCOPE), a high-resolution coupled physical-biogeochemical model. The analysis included variables from J-SCOPE identified in the literature as important for larval crab occurrence: temperature, salinity, dissolved oxygen concentration, nitrate concentration, phytoplankton concentration, pH, aragonite, and calcite saturation state. GLMs were developed with either in situ ocean conditions or environmental exposure histories generated using particle tracking experiments. We found that inclusion of exposure history improved the ability of the GLMs to predict megalopae occurrence 98% of the time. Of the six swimming behaviors used to simulate megalopae dispersal, five behaviors generated GLMs with superior fits to the observations, so a biological ensemble of these models was constructed. When the biological ensemble was used for

INTRODUCTION
The Dungeness crab fishery is one of the most economically important fisheries on the US West Coast, totaling over $200M in 2017(Pacific States Marine Fisheries Commission, 2019. However, this fishery experiences wide inter-annual fluctuations in catch size. For example, in Washington and Oregon, the lowest commercial crab catches were reported in the early 1980s (<1300 metric tons in Washington; <2300 metric tons in Oregon), and record high catches, nearly 10-fold higher, were reported in the 2004-2005 season (>11,300 metric tons in Washington; >15,000 metric tons in Oregon), with moderate variability observed even across consecutive years. 1 , 2 Variable catch rates have been accompanied by large swings in ex-vessel landing values, e.g. from $33.9M in the 2013-2014 season to $74.2M in the 2017-2018 season in Oregon. These large fluctuations have the potential to affect management strategies, fishermen's livelihoods, and local economies (Botsford et al., 1983;Methot, 1986). Due to consistently high fishing effort, it is thought that variable annual catch rates reflect changes in adult Dungeness crab population sizes. The precise causes of this variability are not completely understood but have long been a subject of research (Methot and Botsford, 1982;Botsford and Hobbs, 1995;Higgins et al., 1997).
One factor that has been linked to variability in the Metacarcinus magister fishery is the abundance of the final pelagic larval stage, the megalopal stage, 4 years prior (Shanks and Roegner, 2007;Shanks et al., 2010;Shanks, 2013). Shanks (2013) reported a significant, parabolic relationship between megalopae abundance in coastal habitats and recruitment into the adult M. magister fishery four years later: recruitment into the fishery was maximized at intermediate megalopae abundances, and otherwise was limited either by population levels (low abundance) or density-dependent effects (high abundance).
Every summer the coastline of the Pacific Northwest experiences a shift in wind direction that promotes upwelling of "corrosive" deep water, which is low in oxygen, pH, and calcium carbonate saturation states, to habitat on the continental shelf (Huyer et al., 1979;Feely et al., 2008Feely et al., , 2016Hickey and Banas, 2008). Though these winds vary in intensity and duration, hypoxia (O 2 < 1.4 ml l −1 ; 61 µmol kg −1 ; 62 mmol m −3 ) has increasingly developed over portions of the continental shelf in recent years, with occasional severe hypoxia (O 2 ∼ 0.5 ml l −1 ; 22 µmol kg −1 ; 22 µmol kg −1 ) occurring in M. magister habitat (Grantham et al., 2004;Chan et al., 2008;Connolly et al., 2010). Low pH conditions, sometimes as low as 7.6, and other carbonate chemistry parameters (e.g. delta calcite, pCO 2 , aragonite and calcite saturation states), accompany this low oxygen water, providing an additional stress to organisms (Feely et al., , 2016Harris et al., 2013;Busch and McElhany, 2016;Hodgson et al., 2016;Miller et al., 2016;Siedlecki et al., 2016;Bednaršek et al., 2017Bednaršek et al., , 2020. Summer upwelling conditions have been simulated by an experimental ocean model, JISAO's Seasonal Coastal Ocean Prediction of the Ecosystem (J-SCOPE; Siedlecki et al., 2016). J-SCOPE is a high-resolution (1.5 km horizontal resolution; 40 vertical layers), Regional Ocean Modeling System (ROMS)based, biogeochemical model for Washington and Oregon shelf waters. 3 J-SCOPE has demonstrated measurable skill for ocean conditions on seasonal timescales , and environmental variables from this model have been used to predict habitat for sardines (Sardinops sagax; Kaplan et al., 2016) and hake (Merluccius productus; Malick et al., in prep). Additionally, by pairing the J-SCOPE model system with a particle tracking model and simulated behaviors, environmental exposure history has been shown to influence pteropod survival (Bednaršek et al., 2017). When run in hindcast-mode, J-SCOPE simulates realistic historical ocean conditions and benefits from physical forcing (e.g. boundary conditions, atmospheric forcing, rivers, and tides) that is data-assimilated. These hindcasts, spanning years 2009-2017, output variables specifically relevant to M. magister megalopae, and are used in this study.
With the right combination of prognostic ocean conditions, M. magister megalopae occurrence and habitat could be forecast over a large spatial scale and variable conditions, supplementing field surveys of megalopae (Shanks et al., 2010) and improving forecasting for management applications. Current management relies on pre-season monitoring of crab conditions and realtime catch rates. By incorporating the megalopal stage into management, Dungeness crab fisheries managers in Washington and Oregon would benefit from an ocean model-based tool that would forecast catch with a lead time of >4 years, a timescale that is useful for coordinating state, tribal, and federal managers to develop realistic long-term management strategies, and for fishermen to anticipate changes in the resource (Hobday et al., 2016).
Our aim is to develop a statistical model driven by modeled ocean conditions to predict megalopae occurrence. We tested the hypothesis that occurrence of M. magister megalopae is affected by exposure to both concurrent and recent environmental conditions that are predictable on seasonal timescales. We propose that exposure history may influence megalopae occurrence patterns because prior exposure to lethal or sub-optimal environmental conditions would either increase megalopae mortality or potentially spur avoidance behaviors (which have been observed during settlement, e.g. Sobota and Dinnel, 2000), ultimately resulting in a decreased probability of occurrence. Thus, we investigated whether inclusion of recent exposure history would improve our ability to model megalopae occurrence over using only environmental conditions concurrent with sampling. Hence, simulated ocean conditions over an 8 year period were used to model the occurrence of megalopae within generalized linear models (GLMs) that included two distinct suites of predictor variables: (1) "in situ" GLMs were developed with physical and biogeochemical variables extracted from the ocean model at the times and locations concurrent with megalopae sampling and (2) "exposure history" GLMs included physical and biogeochemical variables extracted along simulated megalopae trajectories from six distinct dispersal experiments. We found that exposure history did improve our ability to model megalopae occurrence, and we assembled a "biological ensemble" of GLMs to generate a superior forecast for megalopae occurrence. An ancillary outcome of this study was the identification of ocean conditions that may influence spatial heterogeneity of M. magister megalopae, identified as the environmental variables that were included as predictors in the occurrence GLMs. The framework developed in this study could be applied to other pelagic species to assess the influence of exposure history on their habitat.

MATERIALS AND METHODS
Our methods rely on a range of interdisciplinary tools and procedures. We provide a flow chart to clarify the order of operations and linkages therein for the methods and results in this paper (Figure 1).

Metacarcinus magister Megalopae Collection
Metacarcinus magister larvae were collected on surveys conducted by NOAA Northwest Fisheries Science Center (NOAA/NWFSC) as part of a larger study of juvenile salmonids and associated nekton (Morgan et al., 2019) at 37 unique stations off of the Washington and Oregon coasts from 2009 to 2017 (Table 1 and Supplementary Table 1). Surveys were conducted during the daytime in late May and/or June for approximately week-long periods. Larvae were collected using a 0.6 m bongo net with 335 µm mesh size. Plankton nets were towed obliquely by letting out 60 m of cable and immediately retrieving it at a rate of 30 m/min while being towed at two knots. Thus, nets were fished from a maximum depth of 20-30 m to the surface, spanning a large portion of the expected depth habitat of megalopae (see details below). Samples were immediately preserved in 5% buffered formalin/seawater solution and returned to the laboratory for analysis. In the laboratory, plankton samples were rinsed and then sorted based on larval developmental stage and enumerated. For more details on field and laboratory methods, see Morgan et al. (2005). For this study, we modeled megalopae occurrence, characterized as presence or absence of megalopae at each sampling station.

J-SCOPE Historical Ocean Simulations
Historical ocean simulations (i.e. hindcasts) were used in lieu of empirical measurements because they provide spatially and temporally continuous ocean conditions. We have conducted extensive model evaluation with ocean observations for all the variables considered here (see Siedlecki et al., 2016, and results below). Modeled environmental variables, obtained from the J-SCOPE ocean model , were used to develop the GLMs to predict megalopae occurrence. Environmental conditions were extracted from historical ocean simulations for 2009-2016 corresponding to either the megalopae collection times and locations (in situ GLM) or along backtracking particle trajectories (exposure history GLMs; more details below). Historical ocean simulations for 2017 were reserved for GLM performance testing.

J-SCOPE Variable Skill Assessment
Validation of J-SCOPE's historical ocean simulations was performed for a wider range of environmental variables than in previous work (e.g. Siedlecki et al., 2015Siedlecki et al., , 2016 Year in Review webpages at http://www.nanoos.org/products/j-scope/). Additionally, assessments were performed for distinct depth intervals ( Table 2) to investigate the skill of particular ocean conditions as they were experienced by in situ or exposure history particles, since the skill of these predictor variables may be relevant to the subsequent performance of the megalopae occurrence GLMs. To accomplish this, empirical observations were matched temporally and spatially to modeled variables.

Ocean Observational Data
Observational data were compiled from regional moorings and surveys conducted from 2009 to 2017 (Figure 2 and FIGURE 1 | Flowchart summarizing the main methods and results of this paper, including validation of predictor variables skill (gray); development of in situ (blue) and exposure history (red) generalized linear models (GLMs); comparison of in situ and exposure history GLMs (purple); assembly of the biological ensemble (yellow); and forecasting of megalopae occurrence and habitat (green). Start and end nodes are in ovals, process steps are in boxes, and decision points are in diamonds. Table SX refers to Supplementary Tables (1-9). AICc = Akaike information criterion corrected for small sample size; AUC = area under the receiver operating characteristic (ROC) curve; EH = exposure history; GLM = generalized linear model; J-SCOPE = JISAO's Seasonal Coastal Ocean Prediction of the Ecosystem.   (DIC), and total alkalinity (TA) were measured at surface and sub-surface locations on Pacific Coast Ocean Observing System (PaCOOS) annual surveys in 2009 and 2010. Temperature, oxygen, salinity, nutrients, phytoplankton, DIC, and TA were also measured at surface and sub-surface locations on West Coast Ocean Acidification (WCOA) surveys conducted in 2011, 2012, 2013, and 2016. For all surveys with temperature, salinity, phosphate, silicate, DIC, and TA measurements, values of pH (at in situ temperature, salinity, and pressure, on the total scale) and aragonite ( ar ) and calcite saturation states ( ca ) were calculated in CO2SYS (Pelletier et al., 2007), using carbonate dissociation constants from Lueker et al. (2000), salinity to boron ratios from Uppström (1974), bisulfate equilibrium constants from Dickson (1990), and the "seacarb" option for fluoride (i.e. Perez and Fraga, 1987 when 33 > T > 10 [ • C] and 40 > S > 10 [PSU], otherwise Dickson and Riley, 1979). Additionally, empirical relationships using the proxy variables oxygen and temperature were developed to estimate carbonate chemistry variables for the region encompassed by the model domain based on calibration data sets collected on all WCOA surveys following the methods described in Alin et al. (2012;Eqs. 1-3):

Supplementary
where T r = 8.1903 and O r = 144.441 are reference values based on the calibration data set, and T = temperature ( • C) and O = oxygen (µmol kg −1 ) are modeled or observed values. These equations were applied at sites where both temperature and oxygen data were available. For phytoplankton observations (measured as either fluorescence or chlorophyll), units were converted following the methods of Davis et al. (2014) prior to comparison with the modeled phytoplankton variable (mmol m −3 ).

Statistical Skill Calculations
Observations and modeled data were paired within a given depth interval (Table 2), and two statistics were calculated to assess model variable skill: (1) the Pearson correlation coefficient (r), ranging from −1 (negative correlation) to 1 (positive correlation), indicates the degree of linear correlation between the observed and modeled variables; and (2) normalized root-mean-square error (NRMSE) estimates the magnitude of the difference between the observed and modeled variables, with the sign indicating the direction of model bias (positive sign indicates model overestimate; negative sign indicates model underestimate; see Supplementary Equations 1, 2).

In situ Variables
To assemble a suite of potential predictor variables for developing the GLM, we selected physical and biogeochemical variables identified in the biological literature as important for the development and survival of M. magister megalopae or a closely related organism, and that exist or can be derived from the J-SCOPE historical ocean simulations. We omitted synthetic or summary variables, such as the PDO or upwelling indices, despite their reported correlations with megalopae abundance Shanks and Roegner, 2007;Shanks, 2013) because we strive for a mechanistic understanding of how fundamental ocean conditions characterize megalopae habitate.g. is it the cooler water temperatures, lower oxygen levels, or more acidified conditions of the upwelled waters that influence habitat? Thus, a total of eight potential predictor variables were  Table 2 for additional information on sampling dates, depths, and measured parameters.
considered for inclusion in our statistical models: temperature, salinity, dissolved oxygen concentration, nitrate concentration, phytoplankton concentration, pH, aragonite saturation state ( ar ), and calcite saturation state ( ca ). Temperature was chosen because it is thought to influence M. magister larval development duration (Moloney et al., 1994; and to be an indicator of advective processes (Ferrari and Ferreira, 2011). Due to the high abundance of M. magister megalopae observed within Columbia River plume fronts, salinity was selected as there is substantial salinity variability within the plume (Morgan et al., 2005), and M. magister larvae are reportedly sensitive to changes in salinity (Pauley et al., 1989).
Dissolved oxygen was selected because M. magister megalopae and juveniles experience negative metabolic effects and high mortality rates, respectively, in hypoxic conditions (Bancroft, 2015;Gossner, 2018). In addition to temperature, nitrate concentration is a strong indicator of upwelling, indicating both nutrient sources for primary productivity and offshore Ekman transport (Hales et al., 2005;Palacios et al., 2013). Phytoplankton concentration was selected as an indicator of cross-shelf and alongshore currents and retentive features, and as a proxy for food availability (Largier et al., 2006;Kudela et al., 2008). pH was included because Miller et al. (2016) found that exposure of M. magister larvae to low pH conditions increased mortality and slowed development rates. We selected ca as a potential predictor because the M. magister megalopae exoskeleton contains calcite (Bednaršek et al., 2020;Boßelmann et al., 2007;Neues et al., 2007). Although ar may not directly affect the condition of the megalopal exoskeleton, we chose to include this variable as a potential predictor because of the impacts it has on the pelagic food web that may influence megalopae development and survival (Riebesell et al., 2000;Fabry et al., 2008).
Five variables (temperature, salinity, dissolved oxygen concentration, nitrate concentration, and phytoplankton concentration) were extracted directly from the J-SCOPE historical ocean simulations at the stations where megalopae were sampled, and then averaged over the upper 30 m of the water column, the sampling depth of the megalopae collections. The remaining three variables (pH, ar , ca ) were calculated with empirically derived formulae that used dissolved oxygen concentration and temperature at the station locations from the J-SCOPE historical ocean simulations (see equations above); these were then averaged over the surface 30 m.

Exposure History Variables Simulated megalopae dispersal
To test our hypothesis that recent environmental exposure history is important in determining megalopae occurrence, we used particle tracking to simulate virtual megalopae dispersal. Since megalopal stage duration is approximately 30 days for temperatures in this region (Poole, 1966;Ebert et al., 1983) and physical uncertainty in the particle trajectories due to unresolved vertical and horizontal advection and diffusion becomes large over approximately the same time span (sensitivity analyses not shown), particles were tracked backward in time for a period limited to 30 days. Future studies may incorporate earlier larval developmental stages (e.g. zoeae) or run full lifecycle individualbased models, but the current study focuses on whether the environmental exposure over the course of the 30 days prior to megalopae collection can be used to improve megalopae occurrence modeling. Thus, we simulated megalopae dispersal with an offline particle tracking model, the Larval TRANSport Lagrangian model (LTRANSv2b; North et al., 2008North et al., , 2011Schlag and North, 2012), driven by external physical forcing, random displacement, and directed swimming behaviors as prescribed by the user. To account for stochasticity, 100 particles were initialized at each of the 37 sampling stations on the last day of each survey ( Table 1). Ocean velocities from J-SCOPE hourly historical ocean simulations were "reversed" (i.e. negated) to force particle advection backward in time, and particle behavior was imposed.
Although researchers have used observations of the horizontal and vertical distribution of larvae to infer swimming behavior, these behaviors have not been precisely characterized. Thus, six particle back-tracking experiments were run to simulate the depth habitats occupied and the range of behaviors potentially exhibited by M. magister megalopae in the 30 days prior to collection ( Table 2). In the first two simulations, particles exhibited diel vertical migration (DVM) behavior to daytime depths of either 30 m ("EH-DVM30") or 60 m ("EH-DVM60"), per depth variations reported in the literature . These behaviors simulated active swimming down to the maximum depth at dawn, maintaining that depth throughout the day (∼16 h), swimming up to the water's surface at dusk, and maintaining the near surface habitat for the remainder of the night (∼8 h; swimming speed = 10 cm/s; Fernandez et al., 1994;Rasmuson, 2013;Rasmuson and Shanks, 2014). In the third simulation, particles sustained a near-surface habitat by constantly swimming upward ("EH-S1") to mimic anecdotal reports of surface aggregations of megalopae in swarms or attached to flotsam (Lough, 1976;Shenker, 1988). Although DVM and surface-following behaviors are most commonly reported in the biological literature for megalopae, we ran two additional simulations that allowed particles to disperse passively (i.e. no swimming behavior) after being initialized either at the surface ("EH-P1") or at 30 m depth ("EH-P30"), which represent the vertical limits of the plankton collection tows. Passive behavior was simulated for two reasons: (1) because the collected organisms were identified to developmental stage (i.e. megalopae) but their precise age was unknown, organisms that had recently molted from the last zoeal stage would have exhibited reduced swimming abilities commensurate with the earlier life stage (Jacoby, 1982) for potentially a large portion of the back-tracking period, which would be more closely approximated by passive dispersal; and (2) to investigate the effects of behavior on environmental exposure history, passive dispersal served as a null model for comparison with the active behaviors. The passive dispersal simulations included a backward random walk in the vertical direction, whose magnitude was calculated at each time step based on the stored local value of vertical diffusion calculated by the J-SCOPE historical ocean simulation. Finally, to account for megalopae who may have molted from the final zoeal stage mid-way through the 30-day period prior to collection, we simulated a behavior for which the particle exhibited DVM (0-30 m depth) for the first 15 days of the backtracking simulation, followed by passive dispersal for the second half of the simulation ("EH-D15P"). For all behaviors, simulated dispersal trajectories were updated every 60 s, and particle locations and ambient environmental conditions were recorded hourly.
Prior to calculating exposure history for the particle tracks, we removed records of particles that had exited the model boundaries, which included particles located in the Columbia River and Salish Sea, due to a lack of confidence in the biogeochemical modeling in those areas (i.e. particles located at longitude > −123.9 • E or < −126.5 • E or latitude > 49.5 • N or < 43.5 • N; 12.2% of exposure history records). Additionally, we removed any particle exposure data that had unrealistic negative values due to extrapolation errors when a particle was located either at the ocean surface or just above the seafloor (0.29% of exposure history records).

Exposure history variable calculations
Two types of exposure history statistics were calculated for each particle and then averaged across all particles (N = 100) initialized at each station. (1) Average environmental conditions for all variables (temperature, salinity, dissolved oxygen concentration, nitrate concentration, phytoplankton concentration, pH, ar , and ca ) were calculated by extracting the ambient environmental field from the J-SCOPE historical ocean simulations along the particle trajectory and then averaging over the entire 30-day simulation.
(2) Severity indices ("SI"), which are a combined metric of the intensity and duration of exposure to stressful conditions, were calculated for a subset of environmental variables (oxygen, pH, ar , and ca ) by multiplying the duration of time (in days) and magnitude beyond an environmental threshold that a particle was exposed to a stressful condition, and then summed over the 30-day backtracking simulation (Hauri et al., 2013;Bednaršek et al., 2017; Supplementary Equations 3, 4). Environmental thresholds for the severity indices were defined as oxygen < 1.4 ml l −1 (O 2 < 62 mmol m −3 ; 61 µmol kg −1 ; i.e. hypoxia), pH < 7.75 (Feely et al., , 2016Hodgson et al., 2016;Miller et al., 2016), ar < 1 (i.e. physical threshold for aragonite dissolution), and ca < 1 (i.e. physical threshold for calcite dissolution).

GLM Development
A GLM was used to identify the modeled environmental variables that best explain the temporal and spatial heterogeneity in M. magister megalopae occurrence in coastal Washington and Oregon waters. Although life stage processes and population dynamics are assumed to be non-linear, linear models are often used to approximate these processes in fisheries science (e.g. Ricker, 1973;Austin and Cunningham, 1981;Guisan et al., 2002;Venables and Dichmont, 2004), and GLMs have been used successfully to predict probability of occurrence for a wide range of species (Brotons et al., 2004;MacLeod et al., 2008;Krigsman et al., 2012;Froehlich et al., 2015). A GLM is a statistical model that relates a combination of predictor variables (i.e., modeled ocean variables) to a response variable (i.e. megalopae occurrence characterized as presence or absence at each sampling station). Because our response variable had a binomial distribution (i.e. megalopae were either "present" or "absent"), we used a logit link function (Eq. 4;Fisher, 1954), where X b is a linear combination of predictor variables (Eq. 5).
To examine the hypothesis that exposure history is an important driver of megalopae occurrence, a suite of potential GLMs were developed using environmental variables from one of two types of experiments ( Table 2): (1) "in situ" ocean conditions extracted from the model at the times and locations where megalopae were sampled; and (2) exposure history statistics for ocean conditions were extracted along particle trajectories from six distinct particle behavior simulations. We included modeled variables from all 12 surveys (2009-2016; i.e. the "calibration survey period") during GLM development. The automated GLM forward and backward stepwise function in Matlab (R2018b; stepwiseglm) was used to evaluate all possible combinations of the suite of predictor variables (eight variables for the in situ models; 12 variables for the exposure history models), and to add and/or remove variables until the best GLM with the lowest Akaike information criterion corrected (AICc) score for small sample size was derived (Burnham and Anderson, 2002). Nonsignificant and/or potentially collinear predictors were retained in GLMs if their removal resulted in an increased AICc score, indicating that they contributed to model fit despite their lack of statistical significance and/or independence from other predictors in the model. AICc scores attempt to balance the inclusion of additional predictor variables that improve model fit but result in increased model complexity by imposing a penalty for variable inclusion to prevent over-fitting and subsequent declines in model prediction performance (Wilks, 1995). AICc scores were used to compare GLMs developed using in situ and exposure history experiments.
To investigate the impacts of year effects due to inter-annual variation in environmental conditions in the northeast Pacific Ocean (e.g. El Nino-Southern Oscillation and an anomalous "marine heatwave" in 2014-2016; Bond et al., 2015;Fisher et al., 2015;Thompson et al., 2018) and variability in J-SCOPE model skill (e.g. Year in Review pages at http://www.nanoos.org/ products/j-scope/) on GLM model development, we conducted a modified sliding window analysis to develop additional GLMs. Here, we modified the calibration survey period to include only 10 of the 12 surveys available between 2009 and 2016 during GLM development, i.e. all combinations of 10 out of 12 surveys [mathematically, C(12,10)] were used to develop an additional 66 GLMs for the in situ and exposure history experiments, following the methods described above.

Biological Ensemble Assembly
Since the primary aim of this study was to generate the best model for forecasting inter-annual megalopae occurrence patterns, we relied on a model performance metric, the in-sample AUC value, to identify high-performance GLMs across experiments (i.e. in situ or exposure history behaviors). The in-sample AUC value ["area under the receiver operating characteristic (ROC) curve, " see Fielding and Bell, 1997] measures model performance on the data used to develop the model. AUC is calculated by comparing the GLM output probability to the observed presence/absence for megalopae at each station. AUC ranges from 0 to 1, with AUC < 0.5 indicating no model skill, and 0.5 < AUC ≤ 1 indicating skill above random chance.
Due to the high performance of several individual GLMs across exposure history experiments, we selected multiple GLMs to form a "biological ensemble" to represent a range of simulated megalopae behaviors for forecasting applications (see below). We used a criterion of in-sample AUC ≥ 0.64 to identify a cluster of top-performing GLMs to include in the biological ensemble. Each member GLM was weighted equally. To evaluate potential collinearity of variables in the final biological ensemble, correlation coefficients were calculated for predictor variables in each member GLM.

Biological Ensemble Performance Evaluation
To quantify the biological ensemble's ability to forecast megalopae occurrence, the ensemble was used to predict megalopae occurrence for the out-of-sample 2017 survey. Each member GLM of the biological ensemble was individually evaluated for 2017, and then the forecasted megalopae occurrence probabilities were averaged across the member GLMs to obtain the biological ensemble forecast for 2017. Specifically, particle simulations using J-SCOPE modeled ocean conditions for 2017 were conducted for each behavior represented by the GLMs in the biological ensemble. Exposure histories for each behavior were incorporated into the relevant member GLM to generate five independent forecasts of megalopae occurrence probabilities for each sampling site. These five sets of probabilities were then averaged (with equal weighting of each member GLM) to forecast megalopae occurrence probabilities for the biological ensemble as a whole. Ultimately, an AUC value was calculated by comparing the biological ensemble's forecasted megalopae probabilities to megalopae occurrence observed on the 2017 survey.

Habitat Forecasting
Finally, megalopae habitat throughout the J-SCOPE model domain was forecast for 2017 using the biological ensemble. Habitat prediction simulations were run for each behavior represented by the GLMs in the biological ensemble, using particles initialized over a grid throughout the x/y domain of the J-SCOPE ocean model. As described above, the particle exposure histories were incorporated into the individual GLMs comprising the ensemble, and then those output probabilities of megalopae occurrence were averaged across the member GLMs to generate a spatially comprehensive forecast of megalopae habitat for the biological ensemble for 2017. Table 2 were constructed using modeled ocean conditions from either in situ megalopae sampling locations or exposure history simulations based on particle tracking experiments with unique behavior and initialization depth. Comparing these GLMs allowed us to

Generalized linear models predicting megalopae occurrence based on the experiments outlined in
Environmental variables were validated in May and June from 2009 to 2017 and binned by megalopae depth habitat, corresponding to distinct exposure history particle dispersal simulations or "in situ" habitat (0-30 m).  Table 2 for details of in situ and exposure history behaviors.
investigate the hypothesis that exposure history is important for characterizing megalopae occurrence. Because these GLMs are built from modeled ocean variables, we begin by reporting skill assessments for these variables to understand their influence on GLM fit and performance.

J-SCOPE Variable Skill Assessment
Environmental variables from J-SCOPE historical ocean simulations performed well (as indicated by r and NRMSE; Table 3) within the depth habitats dictated by the simulated megalopae behaviors ( Table 2) and during the time period when particle dispersal was simulated (May-June). All but one of the variables (phytoplankton) had a significant correlation with the observations at all depth habitats (r > 0.5), and variable performance generally improved with depth. Phytoplankton did not perform as well as the other variables, except near the ocean surface (0-5 m). Year-round model validation showed similar patterns as seen in the May-June validations (Supplementary Table 3).

Environmental Exposure of Megalopae
The exposure histories generated by the particle tracking experiments were driven predominantly by the particle depth habitats, which were determined by particle behavior and initialization depth (Figures 3-5; ANOVA results shown in Supplementary Table 4). On average, particles exhibiting passive dispersal initialized at 30 m depth (EH-P30) inhabited significantly deeper waters (mean particle depth = 47.1 ± 13.3 m (mean ± std); Figure 3 and Supplementary Table 4) and originated farther offshore (mean particle ending isobath = 572 ± 455 m; Figure 4) than particles in any other experiment. The properties of the simulated ocean conditions that these particles were exposed to were ultimately different because of their unique depth habitat (Figure 5 and Supplementary Table 4). Exposure histories for these passive particles were characterized by significantly lower temperature, dissolved oxygen concentration, phytoplankton concentration, pH, ar , and ca , and significantly higher salinity and nitrate concentrations than particles in other exposure history simulations. Consequently, particles in this experiment were exposed to the most severe hypoxic stress and corrosive waters [severity index (SI) for oxygen = 2.14 ± 3.23 hypoxia-days; SI pH = 1.88 ± 1.45 days with pH < 7.75; SI ar = 2.90 ± 2.37 undersaturation-days; SI ca = 1.69 ± 2.13 undersaturationdays). In contrast, the surface-following particles (EH-S1) had a tightly constrained depth habitat within ∼5 m of the ocean surface (mean particle depth = 2.39 ± 0.04 m), which was significantly shallower than any of the other experiments (Figure 3). Their exposure histories were characterized by significantly higher temperature, oxygen, phytoplankton, pH, ar , and ca , and significantly lower salinity and nitrate concentrations (Figure 5 and Supplementary Table 4). Thus, these surface-following particles experienced minimal exposure to hypoxic and corrosive waters (SI Oxygen = 0.02 ± 0.06 hypoxia-days; SI pH = 0.07 ± 0.15 days below 7.75 pH; SI ar = 0.10 ± 0.20 undersaturation-days; SI ca = 0.002 ± 0.004 undersaturation-days). Particles in other experiments were exposed to intermediate environmental conditions at intermediate depth habitats. Passive particles initialized at the surface (EH-P1; mean particle depth = 20.4 ± 20.7 m) had exposure histories most similar to particles in the 30 m DVM experiment (EH-DVM30; mean particle depth = 20.3 ± 1.2 m) FIGURE 3 | Average depth and variability over time was driven by simulated particle behavior and initialization depth. Particles exhibiting DVM behaviors were tightly constrained between 30 m (EH-DVM30; A) or 60 m depths (EH-DVM60; B) and the surface. (C) Passive particles initialized at 1 m (EH-P1) and at 30 m depths (EH-P30) originated from deeper depths and traveled up in the water column over time, while surface-following particles (EH-S1) experienced a tightly constrained depth habitat near the surface. (D) These particles spent the first 15 days of the simulation exhibiting DVM between the surface and 30 m depth, and then they transitioned to passive dispersal (EH-D15P). and the experiment where particles transitioned from DVM to passive dispersal (EH-D15P; mean particle depth = 20.7 ± 7.3 m). Passive dispersal particles initialized at 30 m (EH-P30; mean particle depth = 47.1 ± 13.3 m) experienced conditions most similar to the 60 m DVM particles (EH-DVM60; mean particle depth = 30.8 ± 2.4 m). In situ extracted conditions were most similar to those experienced by passively dispersing particles initialized at the surface (EH-P1).

GLM Comparisons: In situ Versus Exposure History
Inclusion of environmental exposure history during GLM development improved our ability to predict megalopae occurrence. The in situ GLM had the worst model fit (i.e. highest AICc score) and worst in-sample model performance (i.e. lowest AUC) compared to the exposure history models ( Table 4; see additional statistics in Supplementary Table 5). Sliding window analyses that were used to evaluate the influence of individual surveys on GLM fit and performance (Supplementary Table 6) provided further support that, independent of the calibration survey period used to develop the GLM, the in situ model was out-performed by the exposure history models in 65 out of 66 cases (98% of cases).
Among the exposure history GLMs, model fit and insample performance were affected by the type of simulated behavior and the depth at which the particles were initialized ( Table 4 and Supplementary Table 6). Simulations that included passive dispersal (EH-P1, EH-D15P, and EH-P30) had the best model fit (i.e. lowest AICc scores), followed by DVM behaviors (EH-DVM30 and EH-DVM60), and then the surfacefollowing behavior (EH-S1). These rankings of model fit based on simulated behavior were supported by the sliding window analyses as well (Supplementary Table 6). In-sample model performance, however, showed a different pattern of GLM rankings. When all surveys from 2009 to 2016 were used for GLM development, in-sample AUC indicated highest model performance for the EH-P1 GLM, followed by EH-D15P, EH-DVM60, EH-S1, EH-DVM30, and finally EH-P30 GLMs ( Table 4). The sliding window analysis showed similar results, such that the EH-P1, EH-D15P, and EH-DVM60 GLMs generally had the highest performance, but on average, the EH-DVM30 and EH-P30 GLMs performed better than the EH-S1 GLMs (Supplementary Tables 6, 7).  Table 2 for details about exposure history experiments). Particles were initialized at the same 37 sampling stations. Back-tracked origination locations were averaged for all 100 particles initialized at each station, which sometimes resulted in the average origination location being on land; these particles were moved to the nearest shoreline. Initialization locations are uniquely colored for improved resolution of differing dispersal patterns across the spatial domain. Land is shaded gray, and the 200 m isobath is shown for reference.
FIGURE 5 | In situ conditions and particle tracking exposure histories differed among years and among simulations with different particle behaviors and initialization depths. These box and whisker plots show the interquartile range (bounds of the box) with a horizontal line at the median, ends of the vertical lines at the fifth and 95th percentiles, and outliers represented by dots. Averages were calculated across all particles within an in situ extraction or particle tracking simulation, and statistics were calculated across survey averages to show variability over time (2009-2017; see Supplementary Table 4 for ANOVA results).
The GLMs contained different significant predictors depending on whether or not exposure history was included, and which particle tracking behavior was simulated ( Table 4). The predictors included in the GLMs calibrated with all surveys from 2009 to 2016 were oxygen (three occurrences), salinity (two), nitrate (two), SI for ca (one), phytoplankton (one), SI for ar (one), temperature (one), and pH (one). Although the sliding window analysis indicated that the calibration survey period used to develop the GLM had some influence on which predictors were included in the GLMs (Supplementary Tables 6, 8), the relative frequency of the most common predictors was similar to that observed in the GLMs developed using all surveys: oxygen (194 occurrences), salinity (161)

Biological Ensemble Formation
Due to the overall high performance of several GLMs and slight variations in their relative performance when different calibration survey periods were used (Supplementary Table 6), we decided to select multiple models to generate a biological ensemble that we expect to be more robust to interannual variability than any single model ( Table 5). The biological ensemble is comprised of five GLMs with the highest overall in-sample model performance (i.e. AUC ≥ 0.64) from the following exposure history behaviors: 30 m DVM (EH-DVM30), 60 m DVM (EH-DVM60), surface-following (EH-S1), passive dispersal initialized at 1 m depth (EH-P1), and DVM transitioning to passive dispersal (EH-D15P). For the biological ensemble, megalopae abundance was positively correlated with salinity, oxygen, and pH, and negatively correlated with temperature, nitrate, and the SI for ca . Correlation coefficients calculated for pairwise comparisons of predictor variables within a single GLM indicated low dependency among variables (Supplementary Table 9).

Habitat Model Performance and Predictions
When model performance was tested for the biological ensemble using the out-of-sample 2017 survey, the ensemble performed better than random (AUC > 0.5), indicating skill in predicting megalopae occurrence ( Table 5). Each member of the biological ensemble performed better than random (AUC > 0.5; Supplementary Figure 1), and the ensemble as a whole  Table 4 for all models considered).

Experiment
Equation ( Out-of-sample model performance was evaluated for the 2017 survey for each model individually, and for the biological ensemble as a whole when probabilities from all five models were averaged to predict megalopae occurrence. O = oxygen (mmol m −3 ); S = salinity; T = temperature ( • C); N = nitrate (mmol m −3 ); and SI ca = severity index for calcite saturation state (undersaturation-days). See Table 2 for details about exposure history experiments.
performed better than any individual member (AUC = 0.94). An AUC value of 0.94 means that for 94% of the stations where megalopae were found to be present, the model predicted a higher probability of megalopae occurrence than for randomly sampled stations where megalopae were observed to be absent.
Finally, when we simulated dispersal of megalopae initialized throughout the J-SCOPE model domain, and applied their environmental exposure histories within the biological ensemble, we were able to generate a spatially explicit habitat model for Washington and Oregon (Figure 6). Over this larger domain, the biological ensemble predicted relatively high probabilities of megalopae occurrence seaward of the continental shelf break. Additionally, megalopae occurrence probabilities tended to increase with latitude. Markedly low probabilities of megalopae occurrence were predicted near the mouth of the Columbia River and near the Strait of Juan de Fuca.

Importance of Exposure History
Model fit and performance (indicated by AICc and AUC, respectively) improved when megalopae exposure history was used to develop the GLM compared to the in situ GLM, regardless of the type of particle behavior used to simulate megalopae dispersal. This suggests that prior environmental exposure is important to include in addition to in situ conditions when predicting megalopae occurrence. While behavior affected the exposure history of the particles, the type of behavior did not influence the GLM performance as much as the decision to include exposure history itself. A biological ensemble of five top performing exposure history GLMs was created to capture the range of behaviors that best predicted megalopae occurrence. Within the biological ensemble, megalopae occurrence was positively correlated with dissolved oxygen concentration, salinity, and pH, and negatively correlated with temperature, nitrate concentration, and the SI for ca . These predictor variables suggest that megalopae are less common in nutrientrich environments, potentially generated by upwelling of deep waters that are corrosive and hypoxic, or inflow from terrestrial  Table 5). The outline color and orientation of the triangle indicates whether megalopae were observed to be present (orange, upward-pointing triangles) or absent (blue, downward-pointing triangles) at that station. (B) Probability of megalopae occurrence throughout the J-SCOPE x, y model domain forecasted by the biological ensemble (see color bar). Land is shaded gray, and the 200 m isobath is shown.
sources, such as the Columbia River or the Strait of Juan de Fuca, characterized by warmer temperatures and low salinities (Fiedler and Laurs, 1990;Davis et al., 2014).
Below, we discuss in more detail (1) why a biological ensemble was assembled to encompass the range of potential behaviors exhibited by M. magister megalopae, (2) how the predictor variables in the biological ensemble generate a picture of the preferred habitat for M. magister megalopae, and (3) what the limitations of this work are and how they can be addressed in future studies.

Multiple Behaviors in the Biological Ensemble
The biological ensemble consists of the top performing GLMs which represent the breadth of possible behaviors that forecast megalopae occurrence most skillfully ( Table 5). Due to the uncertainty of the precise age of the megalopae at collection, potential life history changes over the 30 days prior to collection may explain why passive, DVM, and surface-following behavior models produced GLMs with strong model fit and performance (Tables 4, 5). M. magister larvae may spend as few as ∼7 days or as many as ∼33 days in the megalopal stage (Poole, 1966;Ebert et al., 1983;, so our 30-day backtracking experiments may have spanned a period when larvae had reduced swimming abilities as zoeae (Jacoby, 1982), or our DVM behaviors may have overestimated their swimming speeds . Including a GLM with a combination of passive and DVM behavior (EH-D15P) only increased the performance of the ensemble as a whole. Additionally, the surface-following model (EH-S1) may have simulated reported phenomena of megalopae attaching to flotsam or forming swarms near the water's surface during the daytime (Lough, 1976;Shenker, 1988;Roegner et al., 2003). Thus, we hypothesize that realistic megalopae behaviors may be more accurately represented through a combination of passive, DVM, and surface-following simulations compared to any one behavior alone.

Modeled Megalopae Habitat Conditions
The biological ensemble, when used to forecast megalopae habitat for the entire J-SCOPE domain and compared to the out-ofsample 2017 survey, performed very well (Table 5 and Figure 6). This spatially explicit model showed increasing probabilities of megalopae occurrence seaward of the continental shelf break. This modeled habitat pattern aligns with reports from the literature that M. magister larvae disperse offshore during early development before traveling back to the continental shelf to settle (Johnson et al., 1986;Pauley et al., 1989;Morgan and Fisher, 2010). On regional scales, low probabilities for megalopae occurrence were predicted near the mouth of the Columbia River, in the Strait of Juan de Fuca, and for nearshore areas in Oregon. Finally, high variability in megalopae occurrence was predicted on kilometer scales. The predictor variables identified in the biological ensemble may provide insight into the environmental conditions that create suitable (and unsuitable) habitat for M. magister megalopae.
The biological ensemble member GLMs can be divided into two depth habitat groups characterized by unique predictor variables. Four of the five models in the biological ensemble were defined by megalopae swimming behaviors that resulted in intermediate depth habitats, in which a small portion of dispersal time was spent at the ocean surface and the majority of time was spent at depth (EH-DVM30, EH-DVM60, EH-P1, and EH-D15P; Tables 2, 5 and Figure 3). In these GLMs, megalopae occurrence was positively correlated with oxygen concentration and/or salinity, or pH. These predictor variables may indicate that megalopae in mid-water depth habitats are sensitive to low oxygen or low pH conditions characteristic of upwelled waters and/or low salinities indicative of terrestrial sources, such as the Columbia River plume or the Strait of Juan de Fuca waters. Preferences for environments characterized by high oxygen, pH, and salinity are generally consistent with the habitat requirements for M. magister described in the literature (e.g. Reed, 1969;Sulkin and McKeen, 1989;Brown and Terwilliger, 1999;Curtis and McGaw, 2012;Descoteaux, 2014;Miller et al., 2016;Gossner, 2018). For example, negative impacts, such as increased mortality, decreased growth, and increased respiration rates, have been demonstrated when M. magister megalopae and juveniles are exposed to hypoxic conditions (Bancroft, 2015;Gossner, 2018). Prior laboratory experiments have also shown increased mortality of M. magister zoeae and megalopae when exposed to low salinity conditions (Reed, 1969;Brown and Terwilliger, 1999), and avoidance of adult crabs to low salinity conditions, except when starved (Curtis and McGaw, 2012). The positive correlation between pH and megalopae occurrence is consistent with reports that M. magister larvae are negatively impacted by exposure to low pH (Descoteaux, 2014;Miller et al., 2016).
In contrast to the intermediate depth habitat models discussed above, the fourth member GLM of the biological ensemble corresponds to megalopae occurrence in near-surface habitats (i.e. the surface-following behavior model, EH-S1). In surface waters, where oxygen concentrations and salinity are relatively consistent across temporal and spatial domains (Figure 5), alternative predictors were identified to characterize preferred habitat, such as minimal exposure to calcite-undersaturated conditions, relatively cool temperatures, and low nutrient concentrations ( Table 5). These habitat preferences again may signal avoidance of upwelled waters, characterized by low ca and rich in nutrients, or terrestrial inputs, with warm temperatures and high nutrient concentrations. The negative correlation between SI for ca and megalopae occurrence is supported by recent work by Bednaršek et al. (2020) which showed that M. magister megalopae may experience external carapace dissolution due to prolonged exposure to more severe calcite saturation state gradients. Several studies have also highlighted the importance of temperature on larval development and survival in M. magister (Sulkin and McKeen, 1989;Brown and Terwilliger, 1999), and our results indicate a preference for relatively cool temperatures in shallow habitats where thermal stress may be more common than at depth (Figure 5). If megalopae occurrence is linked to exposure history via a mortality mechanism, then our results may suggest that megalopae experience lethal temperatures in shallow habitats over the 30day particle tracking simulations. To our knowledge, no studies have investigated the direct effects of nitrate concentrations on megalopae survival or development. We propose that low nutrient concentrations may indirectly define megalopae habitat by serving as a proxy for preferred downwelling regimes (Hales et al., 2005;Palacios et al., 2013) outside of freshwater plumes, or may indicate the presence of food sources (such as phytoplankton and zooplankton), causing a drawdown of nutrients.
A novel approach taken by this study was to evaluate the skill of modeled variables within the specific depth range and season relevant to our study species. Since the skill of the ocean variables influences the skill of the GLMs to predict megalopae occurrence, and ultimately to model preferred habitat, differences in variable skill may help explain differences in predictive power of GLMs developed with exposure histories from unique behavior simulations. Overall, our model validation showed that ocean variable skill generally improved with depth (Table 3), consistent with prior work by Siedlecki et al. (2016), but variables with strong skill in surface waters also generated GLMs capable of good model performance (e.g. EH-S1 GLM in Table 5). Notably, the 0-30 m DVM exposure history experiment (EH-DVM30), the DVM transitioning to passive behavior (EH-D15P), and the in situ model all used ocean variables within ∼0-30 m depth range, and thus the J-SCOPE variable skill was similar for all GLMs, yet both exposure history GLMs had better fit (lower AICc) and higher predictive skill (higher AUC) than the in situ GLMs (Table 4 and Supplementary Tables 6, 7). This finding further supports our conclusion that exposure to recent environmental conditions is important to include in modeling megalopae occurrence.

Model Limitations and Future Work
In this study, the suite of potential predictor variables was limited to (1) variables that were included in the J-SCOPE historical ocean simulations, (2) variables whose skill could be assessed using observational data, and (3) variables, or threshold values for severity indices, identified in the published literature as being important for M. magister megalopae or related species. For example, only microzooplankton concentration is modeled in J-SCOPE, a class of zooplankton which is not the main food source for brachyuran crab larvae (Bigford, 1977;Harms and Seeger, 1989;Sulkin et al., 1998;Casper, 2013), so we omitted this variable from consideration. Regarding the severity indices, we used thresholds to characterize stressful conditions that may not be biologically relevant for M. magister megalopae, due to limited availability of published scientific studies (see discussions in Hettinger et al., 2012;Waldbusser et al., 2015). Given the importance of temperature and salinity on modeling megalopae occurrence, severity indices for these conditions could also be developed if critical thresholds were identified.
Here, we assembled the five best-performing GLMs into a biological ensemble with equal weighting of its members, due to a lack of information about realistic larval behaviors in wild populations. Future in situ behavioral studies of M. magister latestage zoeae and early-stage megalopae would help shape realistic larval behavior in particle tracking simulations and inform the relative weighting of member GLMs in the biological ensemble.
Application of the biological ensemble for modeling megalopae habitat is best applied during the temporal and spatial window of the megalopae observations used in this studynamely, late May (2009May ( -2012 and late June (2009-2017) over the continental shelves of Washington and Oregon. Since marine conditions typically become more stressful as the upwelling season evolves, beginning in ∼mid-April (Austin and Barth, 2002;Hales et al., 2006;Hauri et al., 2015), our model may weight exposure to more stressful conditions more heavily given the relative under-representation of the earlier (May) sampling period in recent years. Additionally, because all megalopae sampling stations were located over the continental shelf, but we applied the biological ensemble to forecast habitat over the entire J-SCOPE model domain, evaluation of the model's prediction of high-quality megalopae habitat in offshore areas would be essential to utilizing the predicted habitat fields offshore.
Finally, this study laid the groundwork for future forecasting of megalopae abundance. To model more complex dynamics, such as abundance, we will apply either generalized linear mixed models (GLMMs), delta-GLMs, or generalized additive models (GAMs), which have relaxed constraints on the types of relationships allowed between the predictor and response variables (Guisan et al., 2002;Venables and Dichmont, 2004;Brodie et al., 2019). Additionally, we will rely on J-SCOPE seasonal forecasts to predict megalopae abundance on seasonal timescales. Since megalopae abundance is correlated with recruitment into the M. magister fishery 4 years later (Shanks and Roegner, 2007;Shanks et al., 2010;Shanks, 2013), improved forecasts of megalopae abundance, generated without arduous field sampling, would extend the management time horizon from seasonal to more than 4 years in advance, potentially promoting increased long-term planning and stability in the fishery (Hobday et al., 2016;Tommasi et al., 2017).

CONCLUSION
Inclusion of environmental exposure history improved our ability to predict megalopae occurrence. Ultimately, a biological ensemble was generated from GLMs developed with multiple behaviors to encompass biologically relevant variations in megalopae dispersal. This biological ensemble showed superior predictive performance (high AUC) relative to individual GLMs. The biological ensemble identified positive correlations between megalopae occurrence and oxygen concentration, salinity, and pH, and negative correlations with temperature, nitrate concentration, and the SI for ca . When considered together, these variables indicate that megalopae habitat is characterized by downwelling conditions seaward of terrestrial inputs, such as Columbia River plume or the Strait of Juan de Fuca waters.

DATA AVAILABILITY STATEMENT
The datasets analyzed in this study are included in this article, in the Supplementary Material, or can be found online (see Supplementary Table 2).

AUTHOR CONTRIBUTIONS
EN was the principal author of the manuscript, designed and performed backtracking experiments with guidance from SS and AH and conducted habitat modeling with input from SS, IK, and CS. Model skill was evaluated by EN, SO, and SS. JF and CM collected and analyzed biological and environmental samples. SA and RF collected ocean condition observations. SA derived ocean acidification formulae. EN, SS, SO, IK, JF, CM, AH, SA, RF, CS, JN, and NB contributed to the writing and editing of the manuscript.