Essential krill species habitat resolved by seasonal upwelling and ocean circulation models within the large marine ecosystem of the California Current System

foundational forage species and their negative response to large-scale climate variations, such as El Niño and marine heatwave conditions. The approach can be easily transferred to other eco-systems or krill species that respond to similar ocean and climate forcing.


Introduction
Euphausiids (hereafter krill) are some of the most abundant species in the world and are directly or indirectly preyed upon by a diverse array of marine predators including commercially important fish, endangered species, and other upper trophic level species (Laws 1985, Baker et al. 1990, Siegel 2000).Given the global importance of krill to marine food webs and ecosystem function (Mangel and Nicol 2000), understanding the drivers of krill distribution and abundance is important.While krill ecosystem interactions have been widely studied in polar to temperate systems (Lawson et al. 2008, Espinasse et al. 2012, Santora et al. 2012b, Bernard et al. 2017), prediction of krill distribution remains challenging because krill form patchy yet dense aggregations that are highly variable in space and time (Okubo 1986, Levin 1994, Marinovic et al. 2002, Lawson et al. 2004, Croll et al. 2005, Benoit-Bird et al. 2013, Atkinson et al. 2014).Variability in krill populations and distributions can be driven by climate, weather, and ocean circulation (e.g.advection and retention), which can influence survival and reproduction (Brinton 1962a, b, Haury et al. 1978, Winder et al. 2001, Hofmann and Murphy 2004, Croll et al. 2005, Dorman et al. 2011, Ressler et al. 2012, Dorman et al. 2015).Quantifying the relationship between environmental drivers and krill distribution patterns can elucidate krill spatiotemporal variability, species interactions and improve prediction of ecosystem condition (Cury et al. 2008).Long-term monitoring of krill is a critical component of many ecosystem assessments globally (Smith et al. 1995, Sydeman et al. 2010, Goericke and Ohman 2015, Ralston et al. 2015, Atkinson et al. 2017), and the development of predictive models (i.e.species distribution models) from these datasets can map variability in krill habitat and inform spatial management of krill predators to achieve ecosystem-based management objectives (Hazen et al. 2018).
In the California Current Ecosystem (CCE), there are two primary krill species, Thysanoessa spinifera (TSPIN) and Euphausia pacifica (EPAC), that are important for the structure and function of the coastal upwelling food web (Brinton 1962a, b, Ruzicka et al. 2012).High connectivity in the CCE food web makes quantifying interannual variability in krill distribution and abundance important for understanding the distribution, productivity, and survival of many krill predators (Sydeman et al. 2001, Schroeder et al. 2009, Santora et al. 2014).Seabirds, whales and commercially important fish species feed on krill (Chess et al. 1988, Fiedler et al. 1998, Tanasichuk 1998, Brodeur et al. 2007, Lee and Sampson 2009, Hertz et al. 2015, Daly et al. 2017), and organisms that feed higher in the food web are rarely greater than two trophic links removed from krill (Field et al. 2006).TSPIN are concentrated mostly over the continental shelf whereas EPAC, the more numerically abundant species, are concentrated over the continental slope (Brinton 1976, Feinberg and Peterson 2003, Lu et al. 2003).Interannual variability in krill abundance was related to upwelling dynamics and primary production (Santora et al. 2011, Garcia-Reyes et al. 2014), and spatial variability was related to proximity to submarine canyons (Santora et al. 2018) and mesoscale oceanographic features (e.g.fronts, eddies; Santora et al. 2012b).Krill hotspots are often located downstream of upwelling centers and/or adjacent to submarine canyons, possibly due to food availability or lack of offshore transport (Santora et al. 2011(Santora et al. , 2018)).Further, these species favor cooler waters (Marinovic et al. 2002, Peterson et al. 2002, Brinton and Townsend 2003), and thus, respond negatively to warm water anomalies (Peterson et al. 2017) including El Niño Southern Oscillation (ENSO) conditions (Brinton andTownsend 2003, Santora et al. 2014).
Despite many investigations into the environmental determinants of krill, few studies have compared krill species-specific habitat associations to test ecological hypotheses and to inform conservation tools for krill predators.Furthermore, species distribution models that assimilate observations of krill abundance and ocean conditions provide a robust framework for evaluating processes influencing distribution patterns and examining biophysical forces controlling abundance.Here, we present a conceptual framework (Fig. 1) and hypothesize that the spatiotemporal variability in krill abundance in the spring/ summer relates to static geomorphic features, preceding winter upwelling dynamics, and mesoscale oceanographic conditions in spring (e.g.currents, temperature at scales 10-1000s km).We predict that winter upwelling conditions that minimize advection and the associated primary production influence krill population size.Further, spring oceanographic conditions likely influence spatial heterogeneity of krill distribution through transport dynamics and water column properties.To test these hypotheses, we integrated data on TSPIN and EPAC distribution with environmental data from remote-sensing and a data-assimilative Regional Ocean Modeling System (ROMS) in the central CCE from 2002 to 2018 (Moore et al. 2011, Neveu et al. 2016).To independently evaluate model results, we assess whether the presence of krill predators relates to modeled krill abundance throughout the modeling domain.The ability to model spatiotemporal patterns in krill across different ocean climates will benefit studies on food webs and climate impacts, and aid spatial management of fish and top predators.Further, the creation of a krill distribution model is a new technique that contributes to modelling the distribution of other krill species that respond to climate forcing, bathymetry or oceanography with our approach being easily transferable to different locations even if environmental variables or relationships vary.

Ecosystem assessment survey
Krill species abundance and distribution was quantified during the NOAA-NMFS Rockfish Recruitment and Ecosystem Assessment Survey (RREAS) during May-June each year (see Sakuma et al. 2016 for survey design).The survey now spans most of the California coast; however, the most frequently sampled 'core' region ranges from Monterey Bay to Point Reyes (36-38°N) where stations were usually sampled multiple times each survey (Fig. 1).The RREAS has been conducted since 1983, but standardized species estimates began in 2002 (Sakuma et al. 2006).Therefore, we focused on the core region from 2002 to 2018.The 35 sampling stations were sampled in 15 yr on average with most stations having data for 15 of the 17 yr (n = 549).
Krill were sampled using a modified Cobb midwater trawl with a 9.5 mm cod-end liner at night (Sakuma et al. 2006, Ralston et al. 2013).At each station, 15 min net tows were conducted with a headrope depth of 30 m. Individual krill from net hauls (5% of the total haul) were sorted to species, and the total number of krill per haul was estimated by extrapolating the species subsample to the total krill volume (Sakuma et al. 2006).For each survey, we estimated the abundance of all krill species (total krill) and EPAC and TSPIN as log-transformed catch-per-unit-effort (ln(CPUE + 1)), which we refer to as CPUE throughout.To examine interannual variability, CPUE across all stations were averaged to give a single annual value for EPAC, TSPIN, and total krill.Preliminary analyses included testing for correlations between EPAC and TSPIN time series using Pearson's correlation and linear regression to test for trends over time, and temporal autocorrelation functions ('acf ' in the R stats package) determined if cycles of variability were present (R Core Team).The climatological mean of EPAC and TSPIN CPUE was determined for each station from 2002 to 2018, to visualize patterns in species distributions.
The ecosystem survey coincides with the breeding period of many resident seabirds, and the presence of some migratory marine mammals.Predator observations used standardized survey methods (Tasker et al. 1984, Yen et al. 2004).During daylight hours, an observer counted seabirds and marine mammals from the ship's bridge (vessel speed > 5 kts).The width of the survey transect was estimated with a rangefinder, and seabirds within a 300 m arc from the bow and 300 m off the side were recorded.All mammals were counted to the horizon.The density of species was calculated (number per 0.9 km 2 for seabirds and number per 3 km for mammals).We included the migrant humpback whale Megaptera novaeangliae and sooty shearwater Ardenna griseus, and locally-breeding Cassin's auklet Pychoramphus aleuticus and common murre Uria aalge in model evaluations.These species consume one or both species of krill in the region (Szoboszlai et al. 2015).While predator surveys occurred during daytime and krill trawls were conducted at night, the datasets are comparable as predator surveys often occur around and between trawls and thus, predator spatial patterns extend over a large area.

Conceptual overview and environmental variables
Our conceptual framework builds upon previous studies in the CCE (Santora et al. 2011, 2012a, 2014, 2018, Schroeder et al.We hypothesize that krill distribution and abundance is related to winter upwelling dynamics, static geomorphic features as well as transport dynamics and water column properties during the survey period. 2013, Friedman et al. 2018), and thus provides an integrative examination of biological and physical factors that influence krill (Fig. 1).We modeled the distribution of krill abundance using three initial models focusing on key hypothesized drivers: 1) geomorphology (Geomorphic model), 2) preceding winter upwelling dynamics (Winter model), and 3) ocean conditions during the survey (May model) (Fig. 1).Following these models and variable selection (described below), we modeled krill abundance using the combined important variables from these three models, referred to as the Full model (models are listed in Supplementary material Appendix 1 Table A2).Environmental variables were separated according to these three types of drivers (geomorphic, winter and May conditions) to determine the type of variables that best explain krill abundance, to put our results in the context of past studies that found one of these types drivers was important, to aid in creating the most parsimonious model, and finally, to determine if a model containing a combination of these variables is more explanatory.Our models used static, remotely-sensed and ROMS variables (Supplementary material Appendix 1 Table A1).All variables were interpolated to the ROMS grid (0.1°, 30-48°N from the coast to 134°W) and monthly averages were calculated (described further in Supplementary material Appendix 1).
In summary, the Geomorphic model included bottom depth (z), bottom rugosity (z_sd), distance from shore and distance to canyon.For the Winter model, the indicators of winter upwelling were sea surface temperature (SST), wind stress curl, isothermal layer depth (ILD), and upwelling indices determined from vertical transport and nitrate flux near the coast.Four winter models were tested to determine the time period that best related to krill CPUE: 1) January physical variables and February chlorophyll-a concentration (CHL) (referred to as model A), 2) February physical variables and March CHL (model B), 3) March physical variables and April CHL (model C), 4) January-February-March (JFM) physical variables and February-March-April (FMA) CHL (model D) (Supplementary material Appendix 1 Table A1, A2).The four separate models were tested using the three different months and the 3-month mean because upwelling and spring bloom phenology varies interannually and it is unclear how this timing relates to krill CPUE during our surveys.For the May model, ocean conditions included water column structure, SST and its spatial variability, advective or transport conditions, and mesoscale ocean variability during May as the majority of the surveys were conducted during this time (Supplementary material Appendix 1 Table A1, A2).

Krill species distribution model
We used a species distribution model to relate environmental parameters to EPAC, TSPIN, and total krill CPUE (Supplementary material Appendix 1 Table A1, A2).We matched krill CPUE from the 35 core stations each year to the nearest value for each environmental variable.Separate models were run for EPAC, TSPIN, and the total krill CPUE due to observed differences in species distributions (Fig. 1).
The total krill model was used to determine what information would be lost if all krill species were not distinguishable, such as in acoustic studies.
We used a popular machine learning and predictive modeling technique, boosted regression trees (BRTs), to model the krill species distribution in relation to environmental parameters.BRTs outperform other regression-based approaches (e.g.generalized linear/additive models), can elucidate complex non-linear relationships, and can handle missing data, outliers, multi-collinearity, irrelevant predictors, and violations of traditional statistical assumptions (Leathwick et al. 2006, De'ath 2007, Elith et al. 2008, Oppel et al. 2012).Despite the utility of BRTs, a disadvantage to using a correlative approach as opposed to a mechanistic approach is the challenge of interpreting results and distinguishing between correlation and causation, and BRTs are prone to overfitting.
BRTs combine many models (each classification tree) and include stochasticity to decrease model variance and improve predictive performance.We followed established protocols for fitting BRTs and used the 'brt.functions'package in R (R Development Core Team, Leathwick et al. 2006, Elith et al. 2008).BRT models were built using a gaussian family appropriate to krill CPUE (the response variable), a tree complexity of 3, bag fraction of 0.6, and a learning rate such that at least 1000 trees were included in the model (Elith et al. 2008).

Model evaluation and variable selection
To assess model robustness and sensitivity, we ran 50 iterations of four-fold cross-validation where 75% of the data was used for training and 25% for testing and report the mean metrics of the 50 model outputs (Brodie et al. 2018).Model results were evaluated using a set of diagnostic metrics.We quantified the relative influence of predictor variables and constructed partial dependency plots to show the influence of a variable after accounting for the average effects of all other variables in the model (Elith et al. 2008).The root mean squared error (RMSE) expressed as a percentage of the maximum abundance value was calculated to compare predictions to observations along with the percentage of deviance explained.
Model results and performance from the three initial models were evaluated before running the Full models.To choose the appropriate time period for the Winter model, we compared model performance and the partial plots from models A-D (Supplementary material Appendix 1 Table A2).We chose the model with the highest explained deviance/lowest RMSE (Supplementary material Appendix 1 Fig.A1a, b) to include in the Full models as the variable partial plots from models A-D were generally similar (Supplementary material Appendix 1 Fig.A3-A5).For TSPIN, model D (the 3-month mean) was used, and for EPAC and total krill, model C was used, but the four models performed similarly for each species (Supplementary material Appendix 1 Fig.A1a, b), and there was often a high correlation between each variable measured over the three different months.Model D likely also accounts for variability in phenology (e.g.upwelling or the spring bloom).Only variables that had a variable percent contribution > 10% were included in the Full models (Supplementary material Appendix 1 Table A2).Two exceptions were made, distance_from_shore and distance_to_canyon were substituted for z in the EPAC and total krill models because the surveys did not extend into deep waters (> 2400 m), and thus, the full offshore range of EPAC was not adequately sampled.With z in the model, the relationship between krill CPUE and depth was not completely accurate (Supplementary material Appendix 1 Fig.A2b), leading to erroneous spatial predictions of high krill CPUE at all deep depths.

Model prediction and extrapolation
For each year from 2002 to 2018, we predicted the Full models onto environmental data from the core sampling region and the U.S. West Coast to compare the distribution and abundance of krill species.Additionally, for each species, we averaged krill predictions across all trawl stations to compare the interannual variability in predictions and observations and computed Pearson's correlations.For each species, we also computed Pearson's correlations between all observations and predictions using station-level data for all years as well as mean predictions and observations at each station to assess spatial capabilities of the model.We also examined krill CPUE under El Niño and La Niña conditions using the Oceanic Niño Index (ONI, Supplementary material Appendix 1).El Niño (La Niña) events are defined when the ONI has 5 consecutive months with an SST anomaly ≥ 0.5 °C (≤ −0.5°C).We focused on moderate (1.0-1.4°C),strong (1.5-1.9°C) and very strong (≥ 2.0°C) SST anomalies.During our study, 2003, 2010and 2016were considered El Niño years, and 2008, 2011and 2012 were considered La Niña years.In addition, 2015 was the year of peak CCE warming associated with the North Pacific marine heatwave (Bond et al. 2015, Di Lorenzo and Mantua 2016, Jacox et al. 2018).To understand the impacts of these climate events, krill CPUE anomalies for 2015 as well as for composites of the El Niño and La Niña years were computed by subtracting the 2002-2018 climatology from the Full model predictions.

Model evaluation and realism assessment with observed predator distribution patterns
Predator abundance by species was averaged from 2002 to 2018 onto our environmental grid.We only included grid cells that were sampled in at least three years to remove anomalous sightings and irregular survey locations.We matched predator abundance to predicted TSPIN and EPAC abundance to test for predator-prey associations using the climatological means from 2002 to 2018 for both predator observations and krill predictions.We compared the kernel density distribution of TSPIN and EPAC abundance for cells where each predator was absent (abundance = 0) and present (abundance > 0).Using a nonparametric twosample Kolmogorov-Smirnov (K-S) test, we compared two one-dimensional probability distributions and quantified the distance between the two distribution functions (Venables and Ripley 2002).A two-sided K-S test was performed to identify significant differences (p < 0.05) between distributions, and an alternative K-S test was used to identify if a distribution was greater or less than the other.Differences in sooty shearwater distributions were not tested due to the low sample size (n = 8) of absence locations.

Temporal and spatial variability in krill observations in the central CCE
Spatially, the mean proportional abundance of krill at each sampling station showed EPAC was proportionally more abundant offshore over the continental slope and TSPIN was more abundant nearshore (Fig. 1, 2).Temporally, there was high variability in mean krill CPUE, where mean CPUE was often higher for EPAC than for TSPIN, and the means were not correlated over time (r = 0.003, p = 0.99).Positive linear trends in krill CPUE over time were significant at the 90% confidence level for TSPIN (p = 0.056, adjusted R 2 = 17%) and total krill (p = 0.069, adjusted R 2 = 15%) but not for EPAC (p = 0.35).There was no significant temporal autocorrelation present in any time series (p > 0.05).

Performance of krill species distribution models
Model performance metrics showed Full models had an explained deviance and correlation between observations and predictions that was higher than the other three models (Supplementary material Appendix 1 Fig.A1).All models tended to overpredict lower abundance values and under predict higher values (Supplementary material Appendix 1 Fig.A1c-h).Though the Geomorphic model for EPAC had a similar deviance explained (35%) to the Full model, the correlation between predictions and observations was lower (Supplementary material Appendix 1 Fig.A1a, c).Additionally, the EPAC Geomorphic model predicted high abundance in all regions with deep bottom depths (Supplementary material Appendix 1 Fig.A2b), and thus, the predictions were inaccurate.The Full model for total krill had an explained deviance that was roughly half (19%) of the individual TSPIN and EPAC models (35%).For TSPIN, the Geomorphic and May models had more explanatory power (21%) than the Winter models on average (mean 14%) but model D was similar (21%).For EPAC, the Geomorphic model (34%) had more explanatory power compared to the Winter (mean 12%) and May models (21%).Winter models A-D had an explained deviance within ~10% of each other within the TSPIN, EPAC and total krill models, with the best performing model being model D for TSPIN and model C for EPAC and total krill.The root mean squared error (RMSE) generally decreased in order from the Geomorphic, Winter, May to the Full model (Supplementary material Appendix 1 Fig.A1b), with the TSPIN and EPAC models having similar RMSE and the total krill RMSE being ~4% lower.
For the three initial models, variable contributions were similar for TSPIN, EPAC, and total krill (Fig. 3), though partial dependence plots revealed that TSPIN and EPAC responded differently to environmental variables (Supplementary material Appendix 1 Fig.A2-A6); this is described by model type below.

Geomorphic model
Bottom depth (z) contributed the most, followed by rugosity (z_sd), distance from shore and distance to canyon (Fig. 3a, Supplementary material Appendix 1 Table A2).TSPIN was more abundant in nearshore waters (shallower, low rugosity and closer to shore), and EPAC and total krill were more abundant offshore (deeper waters, higher rugosity (i.e.slope regions) and further from shore) (Supplementary material Appendix 1 Fig.A2).

Winter model
The variables tested had similar contributions to each species model (Fig. 3b) with the top three variables including a combination of CHL, wind stress curl, SST and ILD (Fig. 3b, Supplementary material Appendix 1 Table A2).Across all months investigated TSPIN CPUE had negative correlations with ILD and SST, and a positive correlation with CHL (Supplementary material Appendix 1 Fig.A3a,  e, f ).EPAC CPUE had positive correlations with ILD, SST, and wind stress curl, and a negative correlation with CHL (Supplementary material Appendix 1 Fig.A4).For total krill, the functional responses were more variable, but there were positive correlations with ILD and wind stress curl, and negative correlations with SST and CHL (Supplementary material Appendix 1 Fig.A5).For the total krill model, multi-modal response curves were more common, suggesting models may be combining species-specific responses into a single multimodal curve.

May model
The most important variables in the May models were the strength of stratification (BV), surface northward wind stress (SVSTR), surface eastward velocity (SU), and sea surface height (SSH) (Fig. 3c, Supplementary material Appendix 1 Table A2).For all species models, there were negative correlations with SSH, surface eastward wind stress (SUSTR), ILD and BV, and a positive correlation with spatial SST variability indicative of fronts (SST_sd) (Supplementary material Appendix 1 Fig.A6).For TSPIN, there were positive correlations with SU and SVSTR, and a negative correlation with surface northward velocity (SV) while the opposite was seen for EPAC.

Full models
For the most important variables, the functional relationships showed TSPIN CPUE was higher in areas with shallower depths, lower SSH, and 10-30 km from shore; EPAC CPUE was higher in areas with higher bottom rugosity, more southward surface wind stress, higher March wind stress curl, westward surface flow and lower April CHL; total krill CPUE was higher in areas with higher March wind stress curl, more southward surface wind stress, intermediate southward surface flow and westward surface flow (Fig. 3d-i, Supplementary material Appendix 1 Table A2, Fig. A7).Again, multi-modal responses in the total krill model were apparent.

Prediction of regional distribution patterns
The Full model predictions at the trawl stations agreed with observations (Supplementary material Appendix 1 Fig.A1fh, R 2 > 70), including both interannual (p < 0.05, EPAC R = 0.88, TSPIN R = 0.96, total krill R = 0.91, Fig. 2) and spatial variability in krill CPUE (p < 0.05, EPAC R = 0.89, TSPIN R = 0.97, total krill R = 0.86, Fig. 4a-c).Model predictions within the sampling region (Fig. 4a-c) and along the U.S. West Coast (Fig. 4d-f ) revealed differences in TSPIN and EPAC distribution and abundance.TSPIN was predicted to be more abundant nearshore with highest abundances within the inner Gulf of the Farallones and lower abundances north of California and south of the Channel Islands (~33°N, 119°W).The predicted distribution of EPAC followed the continental slope with notably high abundances in Monterey Bay, and lower abundances within the inner Gulf of the Farallones.The total krill model generally showed a high abundance of krill from nearshore to over the continental slope.

Influence of ocean climate
During warm years (2015 marine heatwave and El Niño), models predicted lower than average CPUE, and during cold years (La Niña), there was higher than average CPUE, though there were subtle differences between species (Fig. 5).For TSPIN, during 2015, there was lower CPUE within the Gulf of the Farallones, while there was higher CPUE in northern Figure 3. Radar charts of the variable percent contribution (a-f ) and partial dependence plots (g-i) for boosted regression tree models for TSPIN, EPAC, and total krill, including the (a) Geomorphic, (b) Winter, (c) May, and the Full models for (d) TSPIN, (e) EPAC, and (f ) total krill.The 10% threshold is shown by the thick gray line.The partial dependence plots include the top two variables from the Full models with variable percent contribution shown.From left to right, (g) log10 bathymetry (z) and sea surface height (SSH) for the TSPIN model; (h) rugosity (z sd) and surface northward wind stress (SVSTR) for the EPAC model; (i) wind stress curl in March (curl Mar) and SVSTR for the total krill model (also see Supplementary material Appendix 1 Fig.A7).Each line is the mean response curve for the 50 model iterations.Variable acronyms and models are described in Supplementary material Appendix 1 Table A1, A2.
California and southern Oregon (Fig. 5a).During El Niño (La Niña), there was overall lower (higher) TSPIN abundances but the opposite pattern was observed in northern California (Fig. 5b, c).For EPAC, while there was generally lower (higher) than average CPUE during 2015 and El Niño (La Niña) years, there were patches where the opposite patterns were observed (Fig. 5d-f ).The total krill model anomalies were similar to the anomalies in the individual species models (Fig. 5g-i), but interestingly, the positive anomaly in total krill in northern California/southern Oregon in warm years was driven entirely by TSPIN (Fig. 5g, h).Notably, the magnitude of the anomalies was greatest in 2015.

Coherence of modeled krill and predator distributions
The distribution of modeled krill CPUE relative to predator occurrence patterns provided further confirmation of the structural realism of the krill models.Humpback whales were present along the California coast, with the highest numbers from San Francisco to Monterey Bay (Supplementary material Appendix 1 Fig.A8a).Humpback whales were present in areas of higher EPAC (K-S test, D = 0.25, p < 0.0001) and TSPIN CPUE (D = 0.19, p = 0.0011) compared to where whales were absent (Fig. 6a).The resident breeder Cassin's auklet was abundant on the continental shelf from Monterey to north of Point Reyes and offshore of their Farallon Islands breeding colony (Supplementary material Appendix 1 Fig.A8b), and were present in regions of significantly higher TSPIN CPUE (D = 0.19, p = 0.0031; Fig. 6b).Common murres were abundant in the Gulf of Farallones near their breeding colonies with their full distribution spanning from Monterey Bay to north of Point Reyes (Supplementary material Appendix 1 Fig.A8c).Common murres were present in areas of significantly higher EPAC (D = 0.20, p = 0.00071) and TSPIN CPUE (D = 0.47, p < 2.2e-16; Fig. 6c).Sooty   El Niño (2003, 2010, 2016;center) and La Niña (2008, 2011, 2012;right) anomalies in predicted krill ln(CPUE + 1) from the mean CPUE climatology for (a-c) TSPIN, (d-f ) EPAC, and (g-i) total krill.Red (blue) indicates higher (lower) than average CPUE.
shearwaters were the most prevalent species and were observed in almost all of the surveyed grid cells (Supplementary material Appendix 1 Fig.A8d) and thus, no statistical tests were done (Fig. 6d).

Discussion
By applying a conceptual model and habitat modeling framework, we found krill abundance was related to geomorphic features and ocean conditions during winter and spring.Krill species had different spatial distributions, habitat associations, responses to ocean physics and associations with predators.However, both species responded negatively to large-scale climate events, suggesting they share some common population drivers.Our models were evaluated using krill (~36-39°N) and predator observations (~32-41°N) but we acknowledge that predictions beyond this range have yet to be confirmed.Nonetheless, our krill distribution models can aid in fine-tuning essential krill habitat, identifying ecologically important areas and informing ecosystem-based management.Further, integrating remote-sensing and dataassimilative oceanographic model output improved predicted distribution patterns and subsequent evaluation with observations of krill predators.In creating a habitat model for krill, a species often thought to have unpredictable distributions, we provide a framework for modeling krill species in other ecosystems where populations respond similarly to physical forcing (Brinton and Townsend 2003).

Drivers of krill species-specific habitat associations
Static and dynamic physical factors often influence the distribution and abundance of marine species, where similar species partition available habitat to avoid competition (MacArthur 1958).Our modeling results showed differences in EPAC and TSPIN distributions, abundances, and habitat associations, which reinforced our conceptual framework (Fig. 1) and a historic geographic atlas of krill (Brinton 1962a).The spatial habitat partitioning of the two species and the poorer performance of the total krill model suggested it is important to separate krill species before analysis.Further, the multimodal response curves in the total krill model may result from combining species that have different preferences, lending further support for the individual species models.Our modeling results indicated the combination of variables from the three initial models contributed to the most accurate predictions of krill distribution.We discuss our interpretations of species-habitat relationships in combination with past findings below, but acknowledge that our modeling approach cannot distinguish between correlation and causation, and that bimodal biological responses were apparent for some variables as would be expected in upwelling-driven systems (Jacox et al. 2016).
Both species were influenced by geomorphic features, such as bottom depth, rugosity and distance to shore, confirming coastal associations (Santora et al. 2011(Santora et al. , 2012b(Santora et al. , 2018)).Interestingly, distance from canyon was not an important predictor, but information contained within this variable was likely represented by another feature.However, the previous association of krill hotspots and canyons was based on acoustic data, which track finer-scale habitat and aggregation patterns, which are not detected by trawl surveys (Santora et al. 2011(Santora et al. , 2018)).As shown in acoustics studies, our models indicate krill abundance hotspots were oriented in the alongshore direction, associated with the shelf break/ slope and submarine canyons, and were largely attributed to EPAC (Santora et al. 2011(Santora et al. , 2018)).
Krill abundance and distribution were also related to winter upwelling, suggesting an important role of these conditions.Interestingly, EPAC and TSPIN had different responses to ecosystem variability.TSPIN was associated with shallower isothermal layer depth, cooler surface waters and high CHL, while EPAC was associated with the opposite.These opposing relationships could be driven by different spatial distributions (neritic versus oceanic) and associated water properties as a function of distance from shore rather than active habitat selection (i.e.correlation and not causation).For example, isothermal layer depths are deeper offshore where EPAC was abundant and shallower nearshore where TSPIN was abundant.However, the relationship between TSPIN and cold water and high CHL suggests an affinity for coastal upwelling environments, particularly within upwelling shadows.Satellite observations of surface CHL are typically higher on continental shelves and therefore the functional relationship between EPAC and CHL offshore may need closer examination with CHL measurements throughout the water column.Concurrent hydrographic conditions (e.g.surface currents/ wind stress, stratification) likely impacted krill spatial patterns through advection or other transport dynamics.EPAC and TSPIN were associated with upwelling favorable conditions, yet TSPIN was related to onshore surface flow while EPAC was associated with offshore flow, revealing likely impacts of advection or retention processes.EPAC was more strongly associated with southward wind stress, which leads to stronger offshore advection, but also relates to productivity (Chelton et al. 1982, Dorman et al. 2015) and regional mesoscale distribution patterns.

Impacts of climate variability on krill
Shifts in CCE oceanographic conditions and epipelagic forage community structure are impacted by ocean variability, including ENSO (Chavez et al. 2003, Jacox et al. 2016).For both EPAC and TSPIN, the models predicted warm El Niño (cool La Niña) conditions resulted in lower (higher) than average krill abundances along much of the U.S. West Coast, as previously suggested (Ainley et al. 1996, Marinovic et al. 2002, Peterson et al. 2002, Brinton and Townsend 2003).Krill predictions showed fine-scale deviations from this overall pattern, which need to be confirmed with observations (e.g.opposite anomalies for TSPIN between San Francisco and Eureka).The marine heatwave negatively impacted both krill species and also caused unprecedented northward shifts in other species (Friedman et al. 2018, Sanford et al. 2019).Notably, the negative TSPIN anomalies along the coast were almost twice the magnitude of El Niño anomalies (Fig. 5).Decreased krill abundance may be attributed to below average coastal upwelling and primary production, along with unprecedented warm SST (Bond et al. 2015).El Niño events and marine heatwaves have downstream effects on seabirds and marine mammals, where reproductive failures coincided with warm and weak upwelling years (Ainley et al. 2005, Sydeman et al. 2006).This study highlights the potential to use krill spatial predictions as a regional metric of climate effects on ecosystem function with the ability to assess climate impacts before May-June surveys.

Spatial organization of predators in relation to predicted krill abundance
We found significant overlap between observed predator occurrence and areas of higher modeled krill abundance, suggesting predator distributions can be indicative of krill patches or seascapes.Cassin's auklets were observed in areas of high TSPIN abundance, which has been related to higher breeding success (Ainley et al. 1996, Abraham and Sydeman 2006, Sydeman et al. 2013).Humpback whales and common murres, both known to prey switch and shift distributions (Fleming et al. 2016, Wells et al. 2017), were observed in high EPAC and TSPIN CPUE regions.Sooty shearwaters were nearly ubiquitous, being seasonal visitors that were not constrained to any location (Briggs and Chu 1987, Veit et al. 1997, Santora et al. 2012a).Krill-predator habitat requirements should account for life history strategies, with a single krill model likely being insufficient for all predators.For example, migratory species may select seasonal foraging hotspots, while central place foragers target the most profitable prey habitat within their foraging range, but may be less selective on species due to offspring feeding commitments.
A method to define ecologically important areas and essential krill habitat around the world Krill are the foundation of many food webs, essential to ecosystem function and beacons of climate change/variability (Richardson 2008).While the importance of krill is established and mesoscale hotspots are known (Santora et al. 2011(Santora et al. , 2018)), critical aspects of krill ecology remain unclear.In the CCE, this is portrayed by both the ban on krill fishing off the U.S West Coast and the vague definition of krill habitat, 'the area between the coastline and the 2000 m isobath' (PFMC 2008).Elsewhere, the importance of krill is partly conveyed by small fisheries in the northern hemisphere (i.e.British Columbia, Japan) and the largest fishery in the Southern Ocean (Nicol andEndo 1999, Nicol andFoster 2016).Antarctic krill E. superba are a keystone species and important to biogeochemical cycles (Cavan et al. 2019) but are substantially commercially harvested.Their vast spatial distribution makes synoptic surveys to estimate biomass for annual catch limits challenging (Atkinson et al. 2009, Cavan et al. 2019), and even though krill removals are at precautionary levels, there is potential for fishery-related effects on predators, food webs and other ecosystem ramifications (Pauly et al. 1998, Nicol et al. 2012, Plagányi and Butterworth 2012, Hinke et al. 2017, Watters et al. 2020).
The location of krill hotspots around the world are generally known (Dalpadado and Skjoldal 1991, Santora et al. 2018, Atkinson et al. 2019), but our results suggest identifying variability in hotspots may be equally important, as changes in intensity likely result from large-scale processes and may drive trophic cascades.Therefore, for any ecosystem or krill species, the creation of time-varying annual spatial predictions (Hazen et al. 2018) of krill abundance could be useful not only for identifying variations in krill habitat and focusing attention to spatial anomalies but also to guide catch limits or dynamic fishing zones as all krill populations have high interannual variability.Additionally, this would be relevant for understanding the demography of krill predators, avoiding ship strikes or entanglement of whales in fishing gear (Santora et al. 2020), informing offshore development (e.g.wind farms, oil exploration), and guiding adaptive sampling strategies for further research.We provide a framework that is transferable to any location or species, where similar habitat characteristics and other variables (e.g.sea ice) could be incorporated to provide a greater understanding of the long-term spatiotemporal dynamics of krill and the drivers of their distribution and abundance.

Figure 1 .
Figure 1.The mean proportional abundance of the two dominant krill species, T. spinifera (TSPIN) and E. pacifica (EPAC), along the central California coast during the RREAS in May-June from 2002 to 2018.Bathymetry is shown in the background in 600 m intervals to a maximum depth of 4000 m.A conceptual framework for the factors driving krill distribution and abundance is shown on the right.We hypothesize that krill distribution and abundance is related to winter upwelling dynamics, static geomorphic features as well as transport dynamics and water column properties during the survey period.

Figure 2 .
Figure 2. Interannual variation in mean ln(CPUE + 1) of T. spinifera (TSPIN), E. pacifica (EPAC), and total krill from 2002 to 2018 within the core central California survey region shown in Fig. 1.Linear trends are shown by the straight dashed lines.Krill species distribution model predictions from the Full models are shown by the dashed lines with open circles.The Oceanic Niño Index (ONI, blue and red bars) is shown over the same period, where 5 consecutive month periods with values > 0.5(< −0.5) are indicative of El Niño (La Niña) events.

Figure 4 .
Figure 4. Climatology (temporal mean) of predicted ln(CPUE + 1) for krill from the Full model along the central California coast (top panel) and along the U.S. West Coast (bottom panel) from 2002 to 2018 for (a, d) TSPIN, (b, e) EPAC, and (c, f ) total krill.Points in (a-c) are mean observations from the mid-water trawl sampling stations.Bathymetry line contours (contour interval of 500 m) are shown in black.The red box in (d-f ) represents the region shown in (a-c).

Figure 6 .
Figure 6.Kernel density distribution of TSPIN (orange) and EPAC (green) CPUE where predators were present (abundance > 0; solid) and absent (abundance = 0; dashed) in the climatological mean from 2002 to 2018 (Supplementary material Appendix 1 Fig.A8), including the (a) humpback whale, (b) Cassin's auklet, (c) common murre, and (d) sooty shearwater.The 95% confidence interval is shown around each kernel density estimate except sooty shearwater absence locations (d) due to the small sample size.