Predicting cetacean abundance and distribution in a changing climate

Changes in abundance and shifts in distribution as a result of a warming climate have been documented for many marine species, but opportunities to test our ability to forecast such changes have been limited. This study evaluates the ability of habitat‐based density models to accurately forecast cetacean abundance and distribution during a novel year with unprecedented warm ocean temperatures caused by a sustained marine heatwave.

Species distribution models (SDMs) have been used to predict the abundance and distribution of seabirds, fish, sea turtles and cetaceans, and to evaluate risk from human activities (Benson et al., 2011;Forney, Becker, Foley, Barlow, & Oleson, 2015;Gilles et al., 2016;Hammond et al., 2013;Leathwick, Elith, Francis, Hastie, & Taylor, 2006;Louzao et al., 2006;Oppel et al., 2012;Redfern et al., 2013;Redfern, Hatch, et al., 2017;Torres et al., 2015). SDMs developed for cetaceans in the California Current Ecosystem (CCE) have shown forecasting ability at time-scales ranging from weeks  to months , but the range of habitat covariate values for the forecast time were generally within those used to build the predictive models. Evaluating the ability of SDMs to accurately capture responses of cetacean species during extended anomalous periods when conditions fall outside those represented in model development may help to strategically focus efforts aimed at projecting responses to climate change (Silber et al., 2017).
In 2014, waters in the CCE became anomalously warm, often 4°C above the climatological mean, as an unprecedented marine heatwave spread over the area (Bond, Cronin, Freeland, & Mantua, 2015;Di Lorenzo & Mantua, 2016;Leising et al., 2015). Indices of climate and ocean conditions in the northeast Pacific (e.g., Multivariate ENSO Index, Pacific Decadal Oscillation and North Pacific Gyre Oscillation) shifted in 2014 from conditions promoting high primary productivity to warm, less productive conditions . Large-scale warming of the northeast Pacific began in 2013, reached record high temperatures in 2014-2015 and continued into 2016 with the arrival of a strong El Niño Leising et al., 2015;McClatchie, Goericke, & Leising, 2016). In addition to the unusually warm ocean temperatures, dramatic changes in other physical and biological properties resulted in overall changes in ecosystem structure, with notable responses documented for invertebrates, fish, seabirds and marine mammal species (Bond et al., 2015;Cavole et al., 2016;Laake, Lowry, DeLong, Melin, & Carretta, 2018;Leising et al., 2015). The profoundly altered ocean conditions in the CCE from 2013-2016 provide a unique opportunity to evaluate the ability of SDMs to accurately predict changes in marine species abundance and distribution during a period with unusually warm ocean temperatures, and may increase our understanding of how warming waters affect mobile marine predators such as cetaceans.
NOAA's Southwest Fisheries Science Center (SWFSC) has conducted systematic cetacean and ecosystem surveys in the CCE since 1991, and fortuitously, the most recent survey was conducted in the summer and fall of 2014 during the anomalous conditions. To evaluate the ability of SDMs to accurately predict cetacean abundance and distribution during this novel year, we developed predictive habitat-based models of cetacean density based on seven SWFSC shipboard CCE surveys conducted during summer and fall between 1991 and 2009. Models were built for eight taxonomically diverse species with varied habitat associations: short-beaked common dolphin (Delphinus delphis delphis), northern right whale dolphin (Lissodelphis borealis), Pacific white-sided dolphin (Lagenorhynchus obliquidens), striped dolphin (Stenella coeruleoalba), Dall's porpoise (Phocoenoides dalli), blue whale (Balaenoptera musculus), fin whale (B. physalus) and humpback whale (Megaptera novaeangliae). These species also exhibit varied feeding strategies and their targeted prey span a range of trophic levels. For example, in the CCE, blue whales feed exclusively on euphausiids such as the krill species Euphausia pacifica and Thysanoessa spinifera (Croll et al., 2005;Fiedler et al., 1998). Both fin and humpback whales switch prey between krill and small schooling fish such as northern anchovy (Engraulis mordax) and Pacific sardine (Sardinops sagax) (Fleming, Clark, Calambokidis, & Barlow, 2016;Goldbogen et al., 2008;Goldbogen, Pyenson, & Shadwick, 2007;Piatt & Methven, 1992), while the odontocetes feed on a wide variety of squid and fish species (Pauly, Trites, Capuli, & Christensen, 1998). Species that feed on prey spanning a wider range of trophic levels can exploit a broad range of biophysical conditions while foraging, which may have implications for how they are impacted by climate warming (Cavole et al., 2016).
For each species, we built models with three different sets of habitat covariates based on their success in past SDMs and compared the three models' explanatory ability using established metrics. Species abundance and distribution patterns predicted by each model during the anomalously warm year 2014 were assessed using cetacean sighting data collected during 2014. The results are relevant for modelling climate change effects on marine predators and also provide novel insights into species-specific responses to warming ocean waters.

| Survey data
Cetacean sighting data used to build the SDMs were collected in the CCE during the summer and fall (July through early December) of 1991, 1993, 1996, 2001, 2005, 2008 and 2009 using systematic line-transect methods (Buckland et al., 2001). The 1991-2008 surveys covered waters from the west coast of the United States to approximately 556 km offshore (Figure 1a; Barlow & Forney, 2007), while the 2009 survey focused on waters off Southern California (Figure 1b;Carretta, Chivers, & Perryman, 2011). Transect lines were arranged in a systematic grid to provide even coverage of the survey region over the course of each survey. The survey protocol was the same for all years (see Barlow & Forney, 2007;Kinzey, Olson, & Gerrodette, 2000). Research vessels travelled at approximately 18.5 km/hr along predetermined transect lines. Two experienced observers searched with pedestal-mounted 25X binoculars on the flying bridge of the ship; a third observer searched by eye or with hand-held 7X binoculars and recorded cetacean sightings and survey conditions.
The 2014 survey was conducted from 5 August to 9 December and used the same methods as the 1991-2009 surveys, and transect coverage was roughly uniform throughout the study area (Barlow, 2016; Figure 1c).

| Habitat covariates
Continuous portions of survey effort were divided into approximate 5-km segments to create samples for modelling using the approach described by Becker et al. (2010). Species-specific sighting data were assigned to each segment (total number of sightings and average group size), and habitat covariates were derived based on the segment's geographical mid-point. To maintain consistency with the species-specific effective-strip-width estimates derived by Barlow, Ballance, and Forney (2011) and used in this study to estimate cetacean densities, sighting data were truncated at a distance of 5.5 km perpendicular to the track line for the delphinids and large whales, and at 3.0 km for Dall's porpoise (Buckland et al., 2001). Models were built with three different sets of predictor variables to compare their ability to predict changes in abundance and distribution during the novel year. Each set included a different combination of dynamic, bathymetric and geographic covariates as described below.

| Dynamic variables
Environmental variables from a data-assimilative CCE implementation of the Regional Ocean Modeling System (ROMS), produced by the University of California Santa Cruz Ocean Modeling and Data Assimilation group , were used as dynamic predictors as they have proven effective in similar habitat models in this study area . We used daily averages for each variable served at the 0.1 degree (~10 km) horizontal resolution of the ROMS output. Given the broad temporal span of our survey data (1991-2014), we used output from both a historical reanalysis (1980-2010Neveu et al., 2016) and a near-real-time data assimilation system (2011-present; Moore et al., 2013). These two systems differ in the specific data used and in assimilation details, but both provide data-constrained state estimates for our study area. We limited the predictors to those consistent between the two sources : sea surface temperature (SST) and its standard deviation sd (SST), calculated for a 3 × 3-pixel box around the modelling segment mid-point), mixed layer depth (MLD, defined by a 0.5°C F I G U R E 1 Completed transects for the Southwest Fisheries Science Center systematic ship surveys conducted between 1991 and 2014 in the California Current Ecosystem study area. The green lines show on-effort transect coverage in Beaufort sea states of 0-5 for (a) surveys conducted late July through early December in 1991December in , 1993December in , 1996December in , 2001December in , 2005 and 2008 off the US west coast, (b) a smaller scale survey conducted September through December of 2009 and (c) the 2014 survey conducted from 5 August to 9 December deviation from the SST), sea surface height (SSH) and sd(SSH). A simple offset (+0.035 m) was applied to the near-real-time SSH data to match the historical reanalysis dataset, which had a different reference level (Scales et al., 2017).

| Bathymetric variables
Bathymetric data were derived from the ETOPO1 1-arc-min global relief model (Amante & Eakins, 2009). Water depth (m) was obtained for the mid-point of each transect segment. The standard deviation in water depth, a proxy for seafloor slope, was calculated based on values within a 3 × 3-pixel box around the modelling segment mid-point.

| Alternative predictor sets
Three separate models were built for each species based on different combinations of predictor variables. The first included the full set of dynamic and bathymetric variables ("Model 1"; Table 1). These habitat covariates have been used successfully in previous SDMs in the study area .
The second set of predictors included the same variables as the first set but SST was replaced with a bivariate interaction term between SST and latitude ("Model 2"; Table 1). Since the study area spans a range of latitudes, including temperate to subtropical and shallow to deep-water habitats with water masses of different origins (McClatchie, 2014;Reid, Roden, & Wyllie, 1958), a bivariate SST:latitude predictor would account for these effects, and other modelling studies have used similar approaches to account for latitudinal temperature gradients within large study areas (Forney, 2000;Yuan et al., 2017).
The third model used an identical set of predictors as Model 1 but included longitude and latitude as a bivariate term to explicitly account for geographical effects ("Model 3"; Table 1). The inclusion of a spatial interaction term prohibits predictions outside of the study area, but given the prevalence of longitude and latitude in many cetacean modelling studies (e.g., Cañadas & Hammond, 2008;Forney et al., 2015;Hedley & Buckland, 2004;Hedley, Buckland, & Borchers, 1999;Tynan et al., 2005;Williams, Hedley, & Hammond, 2006), we included this set for comparison.
Year was also included as a potential predictor in all three models to capture population trends for species whose abundance has changed substantially during the time period considered in our analyses, including fin whale (Moore & Barlow, 2011), humpback whale (Barlow, Calambokidis, et al., 2011), blue whale (Monnahan, Branch, & Punt, 2015) and short-beaked common dolphin (Barlow, 2016).

| Modelling framework
Generalized additive models (Hastie & Tibshirani, 1990) were developed in R (v. 3.3.2; R Core Team, 2014) using the package "mgcv" (v. 1.8-17;Wood, 2011). Methods largely followed those described in Becker et al. (2016), but we used different (daily rather than weekly) ROMS predictors that were consistently available for the 2014 forecast. We used restricted maximum likelihood (REML) to optimize the parameter estimates and a variable selection process that uses a shrinkage approach to modify the smoothing penalty, allowing the smooth to be identically zero and removed from the model (Marra & Wood, 2011). To ensure that models were not overfit, we also removed variables that had pvalues > 0.05, and then refit the models to ensure that all remaining variables had p-values < 0.05 (Redfern, Moore, et al., 2017;Roberts et al., 2016).
We used two different species-specific modelling frameworks, depending on group size characteristics. For large whales and Dall's porpoise, species that generally occur in small groups, we fit a single response model using the number of individuals per transect segment as the response variable and a Tweedie distribution to account for overdispersion (Miller, Burt, Rexstad, Thomas, & Gimenez, 2013).
For the delphinids, species that occur in groups with large and variable sizes, we fit separate encounter rate and group size models. All transect segments were used to build the encounter rate models, using the number of sightings per segment as the response variable with a Tweedie distribution. Group size models were built using only those transect segments that included sightings, using the natural log of group size as the response variable and a Gaussian link function.
Geographical differences in group size have been observed for many delphinids (Barlow, 2015;Cañadas & Hammond, 2008;Ferguson, Barlow, Fiedler, Reilly, & Gerrodette, 2006), so we modelled group size using a tensor product smooth of longitude and latitude for use in density estimation (Wood, 2003). For both the single response and encounter rate models, the natural log of the effective area searched (described below) was included as an offset to account for TA B L E 1 Potential predictor variables offered in each of the three models developed for each species considered in this analysis

Model 1 Model 2 Model 3
Dynamic both varying segment lengths and the different detection conditions recorded during the surveys.
Density (D; number of animals per km 2 ) was estimated by incorporating the model results into the standard line-transect equation (Buckland et al., 2001): where i is the segment, n is the number of sightings, s is the average group size, and A is the effective area searched: where L is the length of the effort segment, ESW is the effective strip half-width, and g (0)  Given that the greatest source of uncertainty in these models is the interannual variability in population density due to movement of animals within or outside of the study area (Barlow et al., 2009;Boyd et al., 2018;Ferguson et al., 2006), this source of uncertainty was used to produce lognormal 90% confidence intervals for the spatial density estimates . We clipped the prediction grid to the boundaries of the approximate 1,141,800km 2 study area to ensure that predictions were not extrapolated outside the region used for model development, particularly since longitude and latitude were included as potential predictors in one of the three alternative models.

| Model comparison
We used a variety of established metrics to compare the performance of the models built with the three different sets of predictors, including Akaike's information criterion (AIC; Akaike, 1973), REML score, the percentage of explained deviance, root-meansquared error (RMSE), ratios of observed to predicted density, and visual inspection of predicted and observed distributions during the 1991-2009 SWFSC cetacean surveys (Barlow et al., 2009;Becker et al., 2010Becker et al., , 2016Forney et al., 2012). Each of these models was then used to predict on daily 10-km resolution grids of environmental conditions during the 2014 SWFSC survey and averaged as described above for the 1991-2009 models to produce a final 2014 prediction. A previous study estimated how species abundance changed from 1996 to 2014 based on the same survey data used here (Barlow, 2016; Table 2). To assess the models' 2014 forecast ability, we compared overall study area ratios of observed to predicted density in light of the major changes noted by Barlow (2016) and inspected predicted 2014 distribution patterns as compared to the 2014 survey observations. While we can make predictions using all three models, in a true forecast situation, we do not know beforehand which will have the best predictive ability on a future year.
This mock forecast (i.e., simulated forecast of a past event) provided a unique opportunity to help us understand issues associated with a true forecast, for example, when a model with relatively high explanatory power (i.e., model goodness-of-fit) might have relatively low predictive power (i.e., performance on a novel dataset), and vice versa.  (Table 3).
Functional plots for the habitat variables common among the three model types were also generally consistent (see Figure S1.1 in Appendix S1 in Supporting Information). The percentage of explained deviance varied the most among the three species-specific models and tended to be higher for Models 2 and 3 (Table 4). Model 3 exhibited the best overall explanatory performance for five of the eight species and Model 2 for the remaining three species (Table 4).
For all species, the multiyear average density predictions from the best models captured observed distribution patterns as indicated by actual sightings during the 1991-2009 surveys ( Figure 2).
The areas of highest predicted densities varied among species, reflecting ecological differences within the broad study area ( Figure 2).
The upper 90% confidence interval plots for some species revealed major differences from the average density predictions (e.g., Dall's porpoise, fin whale), while others were more similar (e.g., striped dolphin, humpback whale).
Despite the similarity of the three different models' explanatory performance, their ability to predict density during the novel 2014 year varied substantially, particularly for blue and humpback whales (Table 5). For five of the eight species, the models with the best explanatory performance were also those that most accurately predicted study area abundance during the novel year (Table 5).
For short-beaked common dolphin, Pacific white-sided dolphin and striped dolphin, models that included longitude and latitude showed the best explanatory power for 1991-2009; however, models without these static terms exhibited superior performance when forecasting on the 2014 anomalous environmental conditions (Table 5).
With the exception of striped dolphin, the best models' 2014 abundance predictions were within approximately 20% of estimates based on the observed data, with observed-to-predicted density ratios ranging from 1.02 to 1.22 (Table 5). The most accurate density prediction (1.02 ratio) was for short-beaked common dolphin, whose 2014 abundance in the study area more than doubled in comparison with past line-transect estimates (Barlow, 2016; Table 2). This increase in abundance likely reflects an influx of short-beaked common dolphins into the study area, as interannual shifts in distribution with changing ocean conditions have been documented previously (Anganuzzi, Buckland, & Cattanach, 1993;Forney, 2000;Forney & Barlow, 1998;Forney et al., 2012;Heyning & Perrin, 1994). The dynamic short-beaked common dolphin model captured this influx quite well as is evident from the 2014 density plot that shows high predictions throughout major portions of the study area compared to the 1991-2009 average ( Figure 3a). In contrast, the best model-predicted abundance estimate for striped dolphin was approximately half the estimate based on observed data ( Table 5). The line-transect abundance estimate for striped dolphin in 2014 was more than three times higher than the 1996-2008 average (Table 2), and 10 times higher than the 2008 estimate (Barlow, 2016). All three models for this species underestimated the substantial influx of animals into the study area in 2014 by a factor of 2-3 ( Table 5).
The best predictive models also were able to capture the sub-   Black dots in the average plots show actual sighting locations from the SWFSC summer/fall ship surveys for the respective species. The abrupt discontinuity that runs east-west at 40°N that appears in many of the multiyear average density and uncertainty plots reflects the location of the Mendocino Escarpment, a bathymetric feature evident in the majority of the models that included depth or slope as significant predictors  Table 2) and model-predicted density estimates (Table 5).

F I G U R E 2 Continued
Study area 2014 abundance estimates for humpback whales were higher than estimates based on past survey efforts (Barlow, 2016; Table 2) and the best model successfully predicted the documented population increase (observed-to-predicted density ratio = 1.15). The 2014 line-transect abundance estimate for humpback whale was substantially higher than past estimates (Table 2)

| D ISCUSS I ON
The unusual ocean conditions during the 2013-2016 marine heatwave offered a unique opportunity to examine the performance of SDMs for cetaceans in the CCE during novel habitat conditions (Table 3). The three different types of models examined in our study varied in their explanatory and predictive capabilities, providing insights and guidance for the development of robust marine-species SDMs in a changing climate. Below we provide details on the models' performance relative to the ecology of the species included in our study.

| Model comparison: explanatory performance
The 1991-2009 average density predictions captured observed distribution patterns and had generally good performance metrics (Table 4), similar to models developed by Becker et al. (2016). Model 1 (dynamic/bathymetric variables only), which included similar covariates to previous models that have been validated extensively (Barlow et al., 2009;Becker et al., 2012Becker et al., , 2014Calambokidis et al., 2015;Forney et al., 2012), was not selected as the best explanatory model for any of the species considered here. This suggests that the spatial interaction terms included in both Models 2 and 3 picked up additional unexplained variance that resulted in improved performance metrics relative to Model 1.
The Model 1 functional plots for short-beaked common dolphin, blue whale and fin whale showed distinct bimodal distributions for water depth (see Figure S1.1 in Appendix S1 in Supporting Information). Although bimodality can suggest separate ecotypes for a given species, it can also suggest prey switching or the ability to exploit multiple niches, for example, blue whales feed on Thysanoessa spinifera in shallow environments (<100 m water depth) and Euphausia pacifica in shelf-edge and open-ocean environments. In the absence of genetic data, we cannot presently determine whether indeed we are modelling separate ecotypes or providing evidence that these species have plastic distributions.
Model 2, which included the SST:latitude interaction term, provided the best explanatory performance for northern right whale dolphin, Dall's porpoise and fin whale. The interaction term serves as a proxy to account for potential latitudinal changes in the relationship between SST and species densities over the broad study area, which has previously been recognized as important for large study areas Forney, 2000;Yuan et al., 2017).
Model 3, which included a spatial interaction term, had the highest explained deviance for the other five species (Table 4), identifying fixed geographical patterns in the distributions of these species.
This result was expected because models that include longitude and latitude can improve explanatory performance in SDMs (Hedley & Buckland, 2004;Hedley et al., 1999;Tynan et al., 2005); however, they do not necessarily increase our understanding of dynamic species-environment relationships Wood, 2006) or allow for accurate predictions outside of the original covariate space or study area. Therefore, we expected that Model 3 would not be able to capture large shifts in species distributions resulting from the unusually warm conditions.

| Model comparison: predictive performance
For three species (short-beaked common, Pacific white-sided and striped dolphins), Model 3 exhibited the best explanatory performance but had the poorest performance when predicting on the anomalous environmental conditions, and substantially underestimated the absolute number of animals in the study area in 2014 (Table 5). For short-beaked common and striped dolphins, two warm-temperate/tropical species whose distributions shifted northward during 2014, Model 3 was not able to capture the additional influx of animals into the study area, most likely due to the static spatial term. For both of these species, Model 1 (dynamic/bathymetric variables only) predicted the absolute 2014 abundance most TA B L E 5 Model comparison based on the overall metrics presented in Table 4 for the 1991-2009 models and their ability to accurately predict study area abundance for the novel year (2014). The Observed:Predicted (Obs:Pred) 2014 Density Ratios reflect total study area density values. Model-based abundance values were calculated as the sum of the individual model-predicted grid cell abundances throughout the study area based on the model with the best 2014 density ratio. Differences in the model-based abundance estimates shown here and the Barlow (2016) line-transect estimates shown in Table 2 [1991][1992][1993][1994][1995][1996][1997][1998][1999][2000][2001][2002][2003][2004][2005][2006][2007][2008][2009]) is shown in the third panel. Blues represent predicted 2014 density values that were lower than the 1991-2009 average (i.e., <0), light yellow represents density values that were similar to the 1991-2009 average, and reds represent density values that were substantially higher than the 1991-2009 average accurately, particularly for short-beaked common dolphin (Table 5).
Model 1 not only captured the northern shift in distribution of shortbeaked common dolphins, but was also able to forecast the higher densities within the study area (Figure 3a).
Model 1 also had the best predictive performance for striped dolphin, but all three models for this species severely underestimated the substantial influx of animals into the study area in 2014 (Table 5) (Barlow, 2016;Peterson et al., 2006). Although Model 1 included a linear SST function (see Figure S1.1.D in Appendix S1 in Supporting Information), which should have captured the increase in the number of striped dolphins in 2014, the overall predictions were too low.
This could indicate that the magnitude of the species-environment response for striped dolphins changed during 2014, similar to the pattern identified by Boyd et al. (2018) for Dall's porpoise when the amount of suitable habitat within the study area changed substantially. Our study area represents the northern portion of the range of striped dolphin (Mangels & Gerrodette, 1994;Perrin, Scott, Walker, & Cass, 1985),  (Table 5). Interestingly, the 2014 estimated change in abundance from the 1996-2008 average varied among these three species; for northern right whale dolphin it was more than twice the average, for Pacific white-sided dolphin it was about 20% lower than the average, and for Dall's porpoise it was more than half the average (Table 2; Barlow, 2016). For these three species, the best novel predictions were made by Model 2 (SST:latitude interaction term), indicating that accounting for the spatial structure of a dynamic variable could be important when making forecasts during changing climate conditions (Rind, 1998).
Model 2 also made the best novel prediction for fin whale and predicted absolute density to within 5% of what was observed (Table 5). Although a large migratory species, fin whales occur year-round off Southern California (Barlow, 1994;Campbell et al., 2015;Carretta, Forney, & Barlow, 1995;Dohl, Norris, Guess, Bryant, & Honig, 1980;Forney & Barlow, 1998). They also have large home ranges, fast travel speeds, and feed in both coastal and pelagic waters (Calambokidis et al., 2015;Falcone, Diehl, Douglas, & Calambokidis, ). Given their widespread distribution and more dynamic foraging behaviour, it is not surprising that Model 3 (longitude:latitude interaction term) exhibited the worst predictions for the novel year.
Model 3 exhibited both the best explanatory and predictive performance for blue and humpback whales (Tables 4 and 5). For both species, predicted density estimates were within about 15% of what was observed, in contrast to the other two model types that underestimated observed density by a factor of 2-6 ( Table 5) (Table 2; Barlow, 2016). Blue and humpback whales were the only species for which Model 3 made the best novel predictions, consistent with observations that these large migrators return to the same largely coastal foraging grounds every season, perhaps due to learned behaviour and/or the persistence of trophic hotspots (Baker et al., 2013;Barlow, Calambokidis, et al., 2011;Calambokidis et al., 2008Calambokidis et al., , 2015Irvine et al., 2014;Santora et al., 2017).

| Did the models accurately predict changes in species abundance and distribution during an anomalously warm year?
For seven of the eight species considered in this study, the best predictive models developed using 1991-2009 survey data successfully captured shifts in distribution and change in absolute abundance during an anomalously warm year. The average 2014 SST value was substantially higher than in 1991-2009 (Table 3), consistent with multiple studies that revealed dramatic changes in SST and other oceanic properties in 2014 (Bond et al., 2015;Cavole et al., 2016;Leising et al., 2015). The anomalously warm temperatures observed in 2014 are projected to be commonplace by the end of this century, as average SST is predicted to rise from 1 to 6°C by 2,100 Our results suggest that the more dynamic models are most successful at novel predictions for odontocetes and fin whale, while models that include longitude and latitude may be most successful for blue and humpback whales. That said, the inclusion of longitude or latitude in any species habitat model may ultimately limit its forecasting ability under potential climate change scenarios. For example, although it is unclear how trophic hotspots may respond to extreme climate changes, it is likely that upwelling patterns will change (Rykaczewski et al., 2015), and there could be shifts in the locations of such hotspots (Santora et al., 2017) and consequently the foraging grounds for blue and humpback whales. that the concept is more complicated when considering diverse cetaceans at the species level, and deserves further study.
New periods of climatic and oceanic extremes could result in unexpected interactions, state changes in species-environment relationships, or changes in the intensity of the response (e.g., Boyd et al., 2018), and these could foil predictions in novel situations (Myers, 1998). For example, during the 1982-83 El Niño, the coastal stock of bottlenose dolphins (Tursiops truncatus) expanded its range northward into central California and has remained there ever since Wells et al., 1990). This is a process we cannot predict with our current habitat models, as the underlying relationships changed post-El Niño, and we are explicitly modelling stable species-environment relationships. If indeed the intensity of the response of striped dolphins changed in 2014 as we suspect, our ability to accurately predict the large influx of this tropical species into the study area was compromised.
While our findings demonstrate the models' ability to predict changes in distribution and abundance as temperatures increase, it would be useful to know if they can also predict what happens if temperatures decrease. While we are not able to address this in the current study, it warrants investigation, particularly since in addition to warming waters, climate change projections also include greater variation in ocean conditions (Rosenzweig et al., 2008), and thus periods when temperatures may decrease.
This study shows promise for predicting changes in abundance and distribution of mobile marine predators during anomalously warm ocean conditions that are likely to become more common as the climate changes (Intergovernmental Panel on Climate Change, 2007;Oliver et al., 2018;Rosenzweig et al., 2008). However, our predictive accuracy will likely be limited as ocean conditions deviate more and more from past conditions, particularly if underlying species-environment relationships or the intensity of species responses change. These challenges underscore the utility of methods such as those developed by Boyd et al. (2018) that can explicitly take such changes into account. It is also important to recognize the limitations of predictive models for use in marine spatial planning, because risk assessments, mitigation measures and designated marine protected areas may no longer perform as intended (Hobday & Hartog, 2014;Hobday et al., 2018). For effective management, ongoing studies of species distribution and abundance are imperative to allow habitat models to be updated, re-evaluated and improved, and to allow a better understanding of ecosystem dynamics in a changing climate.

ACK N OWLED G EM ENTS
The survey data used in this analysis were collected by a dedicated team from NOAA's Southwest Fisheries Science Center (SWFSC) and we thank everyone who contributed to collecting these data.
Chief Scientists for the surveys included Tim Gerrodette, Susan Chivers and two of the coauthors (JB and KAF). We thank the UCSC Ocean Modeling group for providing ROMS output. The ROMS near-real-time system is supported by NOAA through a grant to the CeNCOOS Regional Association. This project was funded by the Navy, Commander, U.S. Pacific Fleet and by NOAA/ NMFS/SWFSC. We thank Paul Fiedler, Chris Edwards and Chip Johnson for their thoughtful comments on an earlier version of the manuscript.

DATA ACCE SS I B I LIT Y
NOAA's Southwest Fisheries Science Center cetacean sightings and survey effort can be accessed freely on OBIS SEAMAP (https:// seamap.env.duke.edu/). Bathymetric data were derived from the ETOPO1 arc-minute global relief model (Amante & Eakins, 2009) and are available for download at https://www.ngdc.noaa.gov/ mgg/global/relief/ETOPO1/data/. Regional Ocean Modeling System (ROMS) data were derived from the data-assimilative CCE implementation of ROMS produced by the University of California Santa Cruz Ocean Modeling and Data Assimilation group  and are available for download at https://oceanmodeling.ucsc.edu/ ccsnrt.