Characterizing habitat suitability for a central-place forager in a dynamic marine environment.

Characterizing habitat suitability for a marine predator requires an understanding of the environmental heterogeneity and variability over the range in which a population moves during a particular life cycle. Female California sea lions (Zalophus californianus) are central-place foragers and are particularly constrained while provisioning their young. During this time, habitat selection is a function of prey availability and proximity to the rookery, which has important implications for reproductive and population success. We explore how lactating females may select habitat and respond to environmental variability over broad spatial and temporal scales within the California Current System. We combine near-real-time remotely sensed satellite oceanography, animal tracking data (n = 72) from November to February over multiple years (2003-2009) and Generalized Additive Mixed Models (GAMMs) to determine the probability of sea lion occurrence based on environmental covariates. Results indicate that sea lion presence is associated with cool (<14°C), productive waters, shallow depths, increased eddy activity, and positive sea-level anomalies. Predictive habitat maps generated from these biophysical associations suggest winter foraging areas are spatially consistent in the nearshore and offshore environments, except during the 2004-2005 winter, which coincided with an El Niño event. Here, we show how a species distribution model can provide broadscale information on the distribution of female California sea lions during an important life history stage and its implications for population dynamics and spatial management.

INTRODUCTION

Understanding how highly mobile marine predators select and prioritize habitats can be challenging. Despite the proliferation of animal tracking data and near‐real‐time availability of environmental data products, we still have a limited understanding of how environmental heterogeneity can influence the spatial distribution of many species and populations. Habitat models (or species distribution models, SDMs) provide correlative insight into the biophysical features that may drive habitat preference across a wide variety of taxa, scales, and environments (Guisan & Zimmermann, 2000). Habitat models have provided novel tools for assessing and predicting how animals interact with their environment and are increasingly used for ecological and conservation‐relevant research (Barbet‐Massin, Jiguet, Albert, & Thuiller, 2012; Buckley et al., 2010; Dambach & Rödder, 2011; Elith & Leathwick, 2009; Studwell et al., 2017; Zydelis et al., 2011). Most recently, marine SDMs have been used to identify critical habitat of understudied populations, improve our understanding of distributional shifts in habitat under changing ocean conditions, and support commercial and protected species management (Carvalho, Brito, Crespo, Watts, & Possingham, 2011; Eguchi, Benson, Foley, & Forney, 2017; Hazen et al., 2016; Hobday, Hartog, Spillman, Alves, & Hilborn, 2011; Hooker et al., 2011; Skov et al., 2016).

Despite advances in marine SDMs, the complex life histories of many marine species challenge our ability to understand the spatial distributions and patterns for a population during a particular stage or cycle (Ficetola, Pennati, & Manenti, 2013). Modeling habitat suitability for populations of mobile animals with stage‐specific spatial constraints, such as colonial breeders (e.g., seabirds, pinnipeds), has been particularly difficult (Pinaud & Weimerskirch, 2005). As central‐place foragers, these animals are constrained to land during breeding and provisioning stages, that is, nesting colonies, rookery haul‐out (Orians & Pearson, 1979). During this time, habitat preference is likely a function of both prey availability and proximity to the central location (Rosenberg & McKelvey, 1999). Accounting for this place‐based constraint in habitat models is important, as the need to provision offspring directly constrains foraging opportunities, which in turn can have profound implications for the behavior, energetics, and overall reproductive success of individuals.

The California Current System (CCS) is a highly dynamic eastern boundary current that is characterized by oceanographic variability on multiple spatial and temporal scales (Chavez & Messié, 2009; Checkley & Barth, 2009; Schwing, Husby, Garfield, & Tracy, 1991). Here, highly migratory species are known to associate with biophysical features that promote predictable coastal upwelling centers, fronts, and eddies year‐round (Croll et al., 2005; Scales et al., 2017; Yen, Sydeman, Bograd, & Hyrenbach, 2006). California sea lions (Zalophus californianus; hereafter, sea lions) are among the most abundant top predators in the CCS (Carretta, Forney, & Oleson, 2011; Villegas‐Amtmann, Atkinson, Paras‐Garcia, & Costa, 2012; Villegas‐Amtmann, Simmons, Kuhn, Huckstadt, & Costa, 2011), and populations have continued to grow under the Marine Mammal Protection Act (Carretta et al., 2016; Lowry, Melin, & Laake, 2017; Melin, DeLong, & Siniff, 2008; Roman et al., 2013). In the United States, sea lions breed primarily at four of the Channel Islands in southern California (Santa Barbara Island, San Clemente Island, San Miguel Island, and San Nicolas Island). San Miguel Island and San Nicolas Island are the two largest rookeries, accounting for approximately 90% of the pups produced in the United States (Lowry et al., 2017).

The distribution, foraging ecology, and reproductive strategies of sea lions have evolved under the influence of short‐term (i.e., upwelling) and long‐term (i.e., El Niño Southern Oscillation) changes in ocean conditions (Melin et al., 2008; Weise & Harvey, 2008). Throughout the summer breeding months, animals remain close to the rookeries. During the nonbreeding season (August–May), demographic groups spatially segregate. Males disperse from the rookery and are free to exploit productive habitats throughout the CCS (Melin, Delong, Thomason, & Vanblaricom, 2000), while adult females are central‐place foragers for the entirety of the 10‐ to 11‐month lactation period. During this time, females must attend to their young, alternating periodic trips at sea (1–7+ days) with time on land to nurse their pups (McHuron et al., 2016; Melin et al., 2000). Because they are nonmigratory, females are vulnerable to prolonged changes in their foraging environment (Melin et al., 2008). The ability to locate suitable habitat close to the rookery is critical to pup survival (Costa, 2007; Melin et al., 2000; Ono, Boness, & Oftedal, 1987). Suboptimal conditions reduce prey availability, requiring females to alter foraging and attendance patterns. In extreme events, protracted changes to prey distribution and abundance have led to reproductive failures and long‐term population affects (DeLong et al., 1991; Kuhn & Costa, 2014; Lowry et al., 2017; McClatchie et al., 2016; Melin et al., 2008; Trillmich et al., 1991).

Moreover, a majority of sea lion prey items are commercially important species (e.g., northern anchovy (Engraulis mordax), sardine (Sardinops sagax), Pacific hake (Merluccius productus), jack mackerel (Trachurus symmetricus), Pacific mackerel (Scomber japonicus), and market squid (Doryteuthis opalescens)) (Lowry & Carretta, 1999; Lowry, Stewart, Heath, Yochem, & Francis, 1991; Orr, VanBlaricom, DeLong, Cruz‐Escalona, & Newsome, 2011; Weise & Harvey, 2008). A depletion of foraging resources near the colony can lead to increased spatial overlap with fisheries, leading to both direct competition and indirect biological interactions (NMFS, 1997; Weise & Harvey, 2005). Such conflicts include human‐related injuries (Goldstein, Johnson, Phillips, Hanni, & Fauquier, 1999), depredation (loss of commercial and recreational fish), incidental capture in fisheries or bycatch, and entanglement in fishing gear (Beeson & Hanan, 1996; Carretta & Chivers, 2004; Stewart & Yochem, 1987). Many of these interactions are complex and have been difficult to manage (Arthur et al., 2017; Lewison et al., 2015; Maxwell et al., 2015). In particular, bycatch of marine mammals in fisheries that overlap with their foraging grounds has been identified as a critical management issue for the United States fisheries (Beeson & Hanan, 1996). Bycatch mitigation, however, requires a robust understanding of sea lion habitat use and distribution.

While a robust body of literature has documented the biology, ecology, and physiology of female California sea lions (Antonelis, Stewart, & Perryman, 1990; Costa, 1991; Feldkamp, DeLong, & Antonelis, 1989; McDonald & Ponganis, 2013; McHuron et al., 2016; Melin et al., 2000; Villegas‐Amtmann et al., 2011, 2012), these studies focused on foraging characteristics such as dive behavior, trip length, and duration of individuals from a single rookery or in response to El Niño events (Antonelis et al., 1990; Costa, 2007; Melin, Orr, Harris, Laake, & DeLong, 2012; Melin et al., 2008; Ono et al., 1987; Sydeman & Allen, 1999; Trillmich et al., 1991). Although aspects of at‐sea habitat use have been explored (Kuhn & Costa, 2014), our ability to broadly identify suitable foraging habitat of this central‐place forager has been limited. This is in part because the importance of prey species in sea lion diet fluctuates seasonally and annually, making direct observations that coincide with prey distribution difficult (Lowry et al., 1991; Melin et al., 2012; Orr et al., 2011). Because we currently lack the fine‐scale resolution required to correlate foraging habitat with prey distribution, we must rely on the oceanographic processes that serve as proxies to prey distribution (Arthur et al., 2017; Bost et al., 2009).

Here, we couple a multiyear tracking data set with near‐real‐time environmental data to quantitatively characterize and predict the spatial extent of habitat suitability of lactating female sea lions from the two main rookeries in the CCS. We examine broadscale habitat use using satellite tracking data from 72 lactating sea lions to elucidate the biophysical relationships associated with habitat preference. We then develop predictive models of habitat to explore how accessibility and use changes among years. Our findings demonstrate the utility of a marine SDM to identify changes in habitat use of a central‐place forager. Given the importance of habitat use of breeding and provisioning females on population‐level processes, our model can be used to connect at‐sea distribution shifts to changes in population trends, as well as inform species protection and fisheries bycatch management efforts.

MATERIALS AND METHODS

Data sets and tagging methodology

The movement and distribution of California sea lions were examined using ARGOS tracking data from 72 adult lactating females, tagged in November of 2003–2009. Seventeen females were tagged from San Miguel Island (34.0°N, −120.4°W) and 55 from San Nicolas Island (33.3°N, 119.5°W; Table 1). The tracking period for each year lasted between November and February, although the tracking duration of individual animals did not necessarily span this entire time period. All analyses presented were restricted to this time period. Description of animal capture and instrumentation are provided in Kuhn and Costa (2014) and McHuron et al. (2016).

Table 1

Biometric data for the 72 adult female California sea lion trips used in the GAMM, from 2003 to 2009. A foraging trip was defined as the time at sea between haul‐out periods (Villegas‐Amtmann et al., 2008). Maximum trip distance from colony refers to the distance between the colony and the farthest away an animal traveled

Tag ID	Tagging location	Start date of trip	Date end of trip	Trip duration (days)	Max distance from colony (km)	Total trip distance (km)
2103020	San Nicolas	11/16/03	11/24/03	8	119.43	223.3
2103029	San Nicolas	11/14/03	11/20/03	6	80.24	80.24
2103030	San Nicolas	11/18/03	11/22/03	4	37.92	79.15
2103033	San Nicolas	12/16/03	12/26/03	10	102.39	144.4
2103035	San Nicolas	12/3/03	12/9/03	6	78.78	95.11
2103036	San Nicolas	11/14/03	11/16/03	2	8.08	8.08
2103037	San Nicolas	11/18/03	11/23/03	5	64.58	92.85
2104001	San Nicolas	11/20/04	11/28/04	8	57.24	77.36
2104003	San Nicolas	11/2/04	11/16/04	14	456.95	895.12
2104004	San Nicolas	11/3/04	11/12/04	9	135.78	178.93
2104005	San Nicolas	10/31/04	11/10/04	10	100.06	181.96
2104006	San Nicolas	2/23/05	3/8/05	13	84.37	112.27
2104007	San Nicolas	11/24/04	12/10/04	16	168.62	230.69
2104008	San Nicolas	11/13/04	11/22/04	9	208.75	432.93
2104010	San Nicolas	11/1/04	11/9/04	8	95.77	160.51
2104011	San Nicolas	11/9/04	11/18/04	9	49.08	65.22
2104012	San Nicolas	11/2/04	11/12/04	10	126.25	238.12
2105019	San Miguel	12/17/05	12/23/05	6	163.08	227.6
2105020	San Miguel	12/14/05	12/19/05	5	50.13	101.04
2105021	San Miguel	1/5/06	1/15/06	10	369.27	746.64
2105022	San Miguel	11/18/05	11/25/05	7	136.5	273.04
2105023	San Miguel	1/30/06	2/7/06	8	210.3	264.3
2105024	San Miguel	12/29/05	1/5/06	7	105	222.76
2105025	San Miguel	1/8/06	1/13/06	5	165.75	273.12
2105026	San Miguel	1/11/06	1/19/06	8	355.4	705.26
2105027	San Miguel	1/5/06	1/11/06	6	47.8	77.29
2105028	San Miguel	12/21/05	12/29/05	8	321.62	537.21
2105029	San Nicolas	12/6/05	12/11/05	5	40.75	125.31
2105030	San Nicolas	12/13/05	12/25/05	12	255.98	447.27
2105031	San Nicolas	12/26/05	1/3/06	8	194.02	394.32
2105032	San Nicolas	12/2/05	12/7/05	5	49.22	129.57
2105033	San Nicolas	12/5/05	12/10/05	5	13.05	75.43
2105034	San Nicolas	12/6/05	12/12/05	6	8.02	62.97
2105035	San Nicolas	11/17/05	11/23/05	6	57.54	115.6
2105036	San Nicolas	12/25/05	12/31/05	6	143.47	263.84
2105037	San Nicolas	12/3/05	12/9/05	6	78.43	164
2105038	San Nicolas	12/28/05	1/6/06	9	100.79	177.77
2105039	San Nicolas	12/5/05	12/9/05	4	57.93	67.09
2105040	San Nicolas	12/7/05	12/23/05	16	226.93	275.41
2106001	San Nicolas	11/11/06	11/20/06	9	255.41	493.04
2106002	San Nicolas	11/3/06	11/10/06	7	124.96	223.36
2106003	San Nicolas	11/22/06	12/10/06	18	340.42	582.47
2106004	San Nicolas	11/4/06	11/14/06	10	96.23	228.58
2106005	San Nicolas	11/7/06	11/19/06	12	158.29	352.43
2106006	San Nicolas	11/20/06	11/25/06	5	34.05	69.54
2106007	San Nicolas	11/21/06	11/28/06	7	127.12	253.85
2106008	San Nicolas	11/5/06	11/10/06	5	5.02	8.31
2106009	San Nicolas	12/12/06	12/21/06	9	117.92	231.82
2106010	San Nicolas	11/15/06	11/29/06	14	218.31	484.83
2106011	San Nicolas	11/8/06	11/14/06	6	65.82	65.82
2106012	San Miguel	12/17/06	12/24/06	7	133.73	238.52
2106014	San Miguel	12/15/06	12/22/06	7	195.4	349.16
2106015	San Miguel	11/27/06	12/11/06	14	387.18	796.15
2106016	San Miguel	11/11/06	11/19/06	8	92.87	179.56
2106018	San Miguel	12/3/06	12/8/06	5	67.09	178.68
2106020	San Miguel	11/21/06	11/27/06	6	80.47	176.58
2106021	San Miguel	11/14/06	11/19/06	5	137.39	251.58
2107009	San Nicolas	11/21/07	11/27/07	6	121.95	240.23
2107010	San Nicolas	12/22/07	1/4/08	13	108.12	213.83
2107011	San Nicolas	12/31/07	1/6/08	6	41.3	84.56
2107012	San Nicolas	1/4/08	1/15/08	11	502.15	876.47
2107013	San Nicolas	11/28/07	12/8/07	10	243.75	498.97
2107014	San Nicolas	1/12/08	1/18/08	6	158.89	299.49
2107015	San Nicolas	11/11/07	11/19/07	8	74.07	134.17
2107016	San Nicolas	12/28/07	1/9/08	12	456.12	838.08
2107017	San Nicolas	11/25/07	12/5/07	10	248.24	509.06
2108001	San Nicolas	12/12/08	12/18/08	6	94.72	151.28
2108002	San Nicolas	12/14/08	12/21/08	7	102.9	218.35
2108003	San Nicolas	11/18/08	11/30/08	12	415.27	765.57
2108005	San Nicolas	11/19/08	11/29/08	10	185.83	397.52
2108006	San Nicolas	11/18/08	11/26/08	8	123.1	249.3
2108010	San Nicolas	12/18/08	1/10/09	23	120.26	157.86

All data processing and analyses were carried out in the R environment, version 3.3.0 (R Core Team, 2016).

Track filtering and trip identification

Bayesian state–space modeling techniques were used as a filtering method to account for location error (Bailey et al., 2008, 2012; Breed, Costa, Goebel, & Robinson, 2011; Breed, Costa, Jonsen, Robinson, & Mills‐Flemming, 2012; Jonsen, Flemming, & Myers, 2005). In order to reduce autocorrelation of at‐sea locations, final position estimates along each track were generated at 24‐hr intervals resulting in one position per day (Austin, Bowen, & McMillan, 2004; Kuhn & Costa, 2014). All points over land were removed from final track locations. A foraging trip was defined as the time at sea between haul‐outs from the rookery (Villegas‐Amtmann, Costa, Tremblay, Salazar, & Aurioles‐Gamboa, 2008). Kernel density analyses of at‐sea locations were used to determine the home range (95%) and core (50%) for each rookery during the entire tracking period (adehabitatHR package, Calenge (2011).

Quantifying space use

The efficacy of SDMs often depends on the quality and quantity of presence data points as well as the method of selection of absence points or background data (“pseudoabsences”; e.g., random, environmentally, or spatially stratified; Barbet‐Massin et al., 2012). Here, we created suitable habitat models using presence‐only tracking data and generated pseudoabsences from correlated random walk models (CRWs; Aarts, MacKenzie, McConnell, Fedak, & Matthiopoulos, 2008; Hazen et al., 2016; Willis‐Norton et al., 2015). Such absences represent a theoretical null model where sea lions would travel independent of environmental parameters. Comparison of environmental conditions along sea lion tracks and CRWs can test whether animals are selecting habitat based on specific oceanographic variables (Willis‐Norton et al., 2015).

Owing to the difficulties of quantifying habitat suitability for central‐place foragers (Aarts et al., 2008; Matthiopoulos, Harwood, & Thomas, 2005) and of parameterizing correlated random walks that can accurately approximate their movements, this approach has not been widely implemented. Here, we explore the utility and parameterization of CRWs for a central‐place forager in modeling habitat suitability over broad spatial and temporal scales. Random walk trajectories were simulated using the adehabitatLT package in R (Calenge, 2015). Due to the nature of female sea lion movements, CRWs were simulated by trip and by trip phase (i.e., incoming and outgoing portions of each foraging trip). A minimum of 10 CRW simulations were generated per trip and were allowed to move unconstrained except for on land, in which a new location along the trip length was sampled with replacement. Each simulation started at the first observed trip latitude/longitude location and was built iteratively so that the simulated movement was sampled from a normal distribution (e.g., Figure S1). Each simulation maintained the same relative distance, turning angle, and duration in time between successive locations (Calenge, Dray, & Royer‐Carenzi, 2009).

Each simulated trip was weighted based on how closely it resembled the actual sea lion trip. The weight value was calculated as the normalized difference between the actual trip and simulated trip length distance, summed with the normalized distance in net angular displacement of the sea lion and CRW track (Hazen et al., 2016).

The higher the weight value, the more dissimilar the CRW to the actual trip. Such weighted values provided a means of ensuring that the CRWs were at an appropriate distance and direction as compared to the actual movements of sea lions. Simulated trips with weighted values in the upper quartile, and those that crossed land, were removed to ensure that the CRWs were representative of possible movements, distributions, and habitats that sea lions could have encountered (Hazen et al., 2016).

Remotely sensed oceanographic data

Remotely sensed environmental data were obtained for both sea lion and CRW tracks using Xtractomatic (Simmons, 2016). The data sets included time series of sea surface temperature (SST), surface chlorophyll‐a concentrations (Chl‐a), surface winds (v‐component, for upwelling‐favorable conditions), mean sea‐level anomaly (SLA), SST standard deviation (SST SD), SLA standard deviation (SLA SD), bathymetry, and rugosity (bathymetric standard deviation, bathymetry SD; Table 2, see Table S1 for data references). For each oceanographic parameter, a mean value was calculated based on the mean latitude and longitude error 1° longitude × 1° latitude × 1–8 day intervals and centered at the position of each daily SSM‐interpolated sea lion position (Willis‐Norton et al., 2015). The distance of each satellite location from the colony was calculated using great circle distances to account for the Earth's curvature (Kappes et al., 2010).

Table 2

List of environmental variables and hypothesized influence on adult female sea lion habitat selection

Variable	Hypothesized mechanistic link
Sea Surface Temperature (SST)	Description of thermal regime
Sea Surface Temperature Standard Deviation (SST SD)	Mesoscale thermal structure
Chlorophyll‐a (Chl‐a)	Proxy for primary productivity
Eddy kinetic energy (EKE)	Index of mesoscale convergence and divergence, prey retention
Mean sea‐level anomaly (SLA)	Index of mesoscale features
Mean sea‐level anomaly Standard Deviation (SLA SD)	Index of mesoscale variability
Wind (v‐component)	Upwelling‐favorable winds
Bathymetry	Depth to seafloor
Bathymetry Standard Deviation (Bathymetry SD)	Roughness of seafloor
Distance from colony	Index of movement from rookery

We also explored mesoscale structure in surface currents using eddy kinetic energy (EKE), which was calculated from geostrophic current components as follow (Cayula & Cornillon, 1992):

Transformations of variables were explored to ensure data were normally distributed. A logarithmic transformation was required for Chl‐a and EKE. A square root transformation was applied to bathymetry SD.

Generalized additive mixed models

Generalized Additive Mixed Models (GAMMs) were used to quantify the statistical correlation between oceanographic parameters and sea lion spatial distribution (Redfern et al., 2006). GAMMs allow for multiple nonlinear relationships between a response variable and its covariates in a semiparametric manner (Hastie & Tibshirani, 1990; Su, Sun, Punt, Yeh, & DiNardo, 2011; Wood, 2006). The GAMMs link the environmental covariates to animal presence/absence with individual as a nested variable. Specifically, GAMMs were fit with a binomial distribution, logit link function, and random effect of individual sea lion. To avoid pseudoreplication, only one trip per individual (the trip with the most number of locations) and one CRW of that trip (randomly selected) were used in the models. GAMMs were run using the gamm4 package (Wood & Scheipl, 2013).

Because the main focus of this study is on the broadscale habitat use of lactating female sea lions within the CCS, we chose not to run separate models by rookery, but rather to run one model for all individuals. Results provide information on the population‐level habitat associations for the two largest California sea lion rookeries in the CCS.

Model performance metrics

Candidate models were generated based on hypothesized combinations of environmental covariates (Table 2). All variables were tested for multicollinearity using Generalized Variance Inflation Factors (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). The model with the lowest Akaike's Information Criterion (AIC) and highest receiver operating curve (ROC) area under the curve (AUC) statistic was run 40 times with a 1:1 ratio of randomly chosen simulated tracks for each foraging trip (outgoing and return) to examine variability. Model validation using ROC curves and AUC statistics was calculated using the ROCR package in R (version 1.0‐7; Sing, Sander, Beerenwinkel, & Lengauer, 2015).

Habitat models and predictive surfaces

Predictive surfaces were generated daily and were fit over a set of time‐matched environmental data that corresponded to each November–February satellite tracking period between 2003 and 2009. The spatial resolution of each predictive surface was set to 0.25°, the lowest common resolution of environmental data (Table S1). Daily surfaces were averaged for each November–February period, generating a total of six winter habitat maps. Relative habitat suitability was scaled from 0 (unsuitable) to 1 (highly suitable). Cumulative mean and standard error (SE) suitability maps show the variability associated with model predictions.

RESULTS

General habitat use

From 2003 to 2009, adult female sea lions were tracked for 14–131 days (mean 56.9 days ± 24.6 SD). The average number of foraging trips per female was 13.3 ± 6.9 SD trips, with a mean trip duration of 8.4 days (±3.6 days SD) (Table 1). The maximum straight‐line trip distance from the rookery ranged from 5.0 to 502.2 km (mean 149.4 km ± 115.0 km SD). Dispersal primarily extended north/northwest of each colony (San Miguel and San Nicolas islands, Figure 1a); however, core areas of use remained closest to the colonies (50% UD, Figure 1b,c). While most individuals favored nearshore habitat from the Southern California Bight and along the mainland coast, several individuals from each colony were tracked offshore into waters greater than 500 m depth (Figure 1b,c).

Figure 1

(a) Daily locations from all foraging trips of lactating female California sea lion tracks (Zalophus californianus, n = 72), from November to February 2003–2009 (color‐coded by year of deployment) displaying one location per day. Distribution of tracks by colony: (b) San Miguel (n = 17) and (c) San Nicolas Islands (n = 55). Deployment locations shown as black circles. The 500‐m isobaths are shown in light gray, and the 50% and 95% kernel density utilizations for each colony are shown in dark gray

Model predictions

The best‐fitting model included SST, Chl‐a, bathymetry, EKE, SLA, and SLA SD and explained 46% (R 2 = .46) of the total deviance (Table 3) with an AUC = 0.91 (Figure S2). Partial response curves (Figure 2) show that the probability of sea lion occurrence was greatest with cool SST (<14°C), productive waters (i.e., Chl‐a ranging from −0.5 to 1.0 mg/m3), and shallow depths (<500 m below sea level). Sea lions were also associated with increased SLAs (0.05–0.1 cm) and EKE, while SLA SD (i.e., index of mesoscale variability) was negatively associated with sea lion occurrence. Distance from the colony was considered in a competing candidate model, but surprisingly was a less important predictor of sea lion habitat than bathymetry. Overall, SST, EKE, and bathymetry were the most consistently significant predictors of sea lion habitat (Table 3).

Table 3

Selection diagnostics from the final Generalized Additive Mixed Model (GAMM). Model was run 40 times to examine the number of times a variable was significant (n‐significant). All variables represent a p value <.001. Values listed as follow: mean (min – max)

Variable	Effective degrees of freedom (edf)	Chi‐squared	n‐significant (n/40)
SST	2.7 (1.0–3.8)	160.4 (120.6–198.8)	40
Chl‐a	2.9 (1.0–3.9)	20.0 (1.5–42.3)	24
EKE	3.3 (1.0–3.8)	57.3 (25.1–80.7)	40
SLA	3.1 (1.0–3.9)	39.3 (19.0–70.1)	38
SLA SD	1.2 (1.0–3.5)	27.0 (4.9–71.7)	33
Bathymetry	3.6 (3.3–3.9)	146.2 (106.5–183.2)	40

R 2 = .46 (.36–.54).

AIC = 864.92 (712.61–1014.59).

AUC = 0.91 (0.88–0.93).

Figure 2

GAMM partial plots showing the relative habitat suitability for adult female California sea lions in response to: sea surface temperature (SST, °C), Chlorophyll‐a (Chl, log mg/m3), eddy kinetic energy (EKE, cm2/s2), sea‐level anomaly (SLA, cm), sea‐level anomaly standard deviation (SLA SD, cm), and bathymetry (m). Gray shading represents the 95% confidence intervals for the fitted relationships

The combined influence of these biophysical parameters is evident in the broadscale spatial habitat predictions. Throughout each winter, highly suitable habitats were evident in the near to offshore environments. Available habitat was identified along the California coast, from the northern Channel Islands, up through Monterey Bay. Two exceptions were the winters of 2004–2005 and 2008–2009. During both time periods, habitat suitability diminished within the Southern California Bight (Figure 3b,f). In the winter of 2004–2005, habitat suitability shifted offshore and north of Monterey Bay (Figure 3b). In the winter of 2008–2009, suitability was most concentrated along the nearshore central California coast. Prediction errors across all sampled years were highest in offshore and within the Southern California Bight (Figure 4b).

Figure 3

Spatial habitat predictions of adult female California sea lions by year. Maps show relative habitat suitability for female California sea lions during foraging trips, based on environmental data from November to February, from 2003 through 2009. Suitability is scaled from 0 (unsuitable) to 1 (highly suitable)

Figure 4

Spatial habitat predictions of adult female California sea lions, averaged over all tracking periods (November–February), from 2003 to 2009. (a) Mean spatial prediction of relative habitat suitability for female California sea lions, from November to February 2003–2009, based on a suite of environmental data. Suitability is scaled from 0 (unsuitable) to 1 (highly suitable); (b) Standard Errors in spatial prediction from November to February 2003–2009. The spatial resolution of predictive surfaces was set to 0.25°, the lowest common resolution of environmental data

DISCUSSION

Characterizing habitat associations for a mobile marine predator can be challenging for animals whose movement patterns change among life stages. Modeling habitat suitability for spatially constrained foragers is a particularly complex exercise because the options to respond to their environment are limited (Kappes et al., 2010; Pinaud & Weimerskirch, 2005). Using multiyear satellite tracking data and near‐real‐time environmental data, we developed a species distribution model to identify the oceanographic conditions that characterize and predict foraging habitat for lactating sea lions from the two largest rookeries in the CCS.

The foraging behavior of many central‐place foragers has evolved to repeatedly exploit areas within proximity to breeding grounds where resources may be, to some degree, spatiotemporally predictable, as they offer higher, more efficient levels of energy acquisition and thus maternal provisioning (Baylis, Page, McKenzie, & Goldsworthy, 2012; Bonadonna, Lea, Dehorter, & Guinet, 2001; Chilvers, 2008; Irons, 1998; Lowther, Harcourt, Hamer, & Goldsworthy, 2011). For female California sea lions, we found relationships with static and dynamic environmental covariates that suggest they repeatedly target areas of enhanced productivity. Sea lions preferentially selected cold SST (<14°C), shallow depths, and elevated chlorophyll‐a values, common proxies for upwelling along the continental shelf, where tightly coupled biophysical processes drive the development of a robust food web along central and southern California (Ainley, Sydeman, Parrish, & Lenarz, 1993; Sydeman & Allen, 1999). Coastal upwelling offers seasonally predictable resources on broad spatial scales and has been shown to influence the distribution and abundance of top predators, including several pinniped species off the California coast (Sydeman & Allen, 1999). As California sea lions are known to be shallow water, epipelagic foragers (Feldkamp et al., 1989; Kuhn & Costa, 2014), the combination of shallow depths and upwelling of cold, nutrient waters along the shelf would provide seasonally reliable prey resources. Results from our spatial predictions identified a high degree of suitability close to the breeding colonies and along the California coast, suggesting a strong preference for the nearest predictable and most profitable areas for nursing females.

Interestingly, we found positive associations with metrics of mesoscale activity (e.g., eddy kinetic energy and sea‐level anomalies), suggestive of shelf‐break and offshore foraging, where physical convergence processes can serve as useful foraging patches for top predators due to the aggregation of resources (Bailleul, Cotté, & Guinet, 2010; Hyrenbach, Forney, & Dayton, 2000). Previous top predator studies have documented associations with eddies [e.g., seabirds, Yen et al. (2006); sea turtles, Polovina et al. (2006); and other pinniped species, Fadely, Robson, Sterling, Greig, and Call (2005); Ream, Sterling, and Loughlin (2005)]. While these features are more ephemeral in nature, their persistence in the CCS has been well documented (Batteen, 1997; Lynn & Simpson, 1987; Strub & James, 2000). To our knowledge, this is the first study to detect an association between California sea lion foraging habitat and mesoscale activity. Sea lions are known to display extensive intraspecies variability in their at‐sea movements, behaviors, and distributions (Kuhn & Costa, 2014; McHuron et al., 2016; Melin et al., 2008). This is especially true during lactation, as females must continually expand and adjust their foraging behavior in response to prey movements (Melin et al., 2008) and direct competition for resources under a restricted range and large population size. While finer‐scale studies are necessary to explore an association with eddies in the CCS, it is possible that sea lions may utilize these mesoscale features for foraging opportunities while at sea.

Our modeling approach demonstrates preferred habitat for lactating sea lions is broadly associated with the continental shelf, but exhibits a degree of spatial modulation across years, most notably during basin‐wide environmental perturbations. Considerable studies have documented interannual changes in lactating female foraging behavior in response to El Niño Southern Oscillation (ENSO) events (Costa, Antonelis, & DeLong, 1991; DeLong et al., 1991; Kuhn & Costa, 2014; Melin et al., 2008; Trillmich et al., 1991); however, our ability to broadly identify spatially explicit changes has been limited. For example, higher than normal upwelling and subsequent high productivity during the 2007–2008 La Niña winter provided the most spatially robust habitat availability for sea lions throughout our study period. By the winter of 2008–2009, La Niña conditions had weakened; however, central and northern California maintained stronger than normal upwelling (McClatchie et al., 2009). This most likely drove prey resources north, thereby weakening habitat suitability closest to the breeding colonies. As a result, females would have had to travel farther north to reach access profitable habitat.

Our findings suggest that, during El Niño events, females may face serious limitations in habitat accessibility due to the spatial constraints of the breeding colony. Warmer sea surface temperatures, reduced productivity, and elevated SLA have been associated with reduced food availability in the Southern California Bight (Schwing et al., 2006; Thomas & Brickley, 2006; Trillmich et al., 1991). Suboptimal conditions require females to alter foraging and attendance patterns, with females moving further offshore and north of Monterey Bay to find food (Kuhn & Costa, 2014; Melin et al., 2008), which may result in pup abandonment (Melin et al., 2008). Our model predictions during the 2004–2005 winter reflected this pattern, as a weak El Niño appeared to reduce productivity and prey resources within the nearshore environment (Kuhn & Costa, 2014), shifting habitat suitability away from the southern California coast and north of Monterey Bay. This displacement is consistent with previous studies that identified an increase in adult female abundance in central and northern California during El Niños in response to better foraging conditions (Lowry & Forney, 2005; Sydeman & Allen, 1999).

Although our model predictions extend from 2003 to 2009, the distribution patterns and feature associations we describe here can be used to track shifts in distribution under future ocean conditions. In a warming ocean, alterations to suitable habitat—smaller in extent and further from the breeding colonies—may have severe consequences on the population dynamics of this species. Prolonged environmental change will result in increased energetic costs and decreased reproductive success with potentially long‐term population declines (Hazen et al., 2012; Lowry et al., 2017; McHuron, Mangel, Schwarz, & Costa, 2017; Melin et al., 2008, 2012). Recent studies have documented an increasing number of pups born on Año Nuevo Island, an important central California haul‐out site that supports a small breeding population (Lowry et al., 2017; McHuron, Block, & Costa, 2017). However, there is limited space for the population to expand northwards, as the Channel Islands represent the only islands with enough available space to support such large breeding rookeries under a growing population.

Habitat loss and large sea lion concentrations may lead to increased pressure on coastal fisheries, due to overlap with commercially important prey species, with a potential to impact to top‐down food web dynamics (Lowry & Forney, 2005; Lowry et al., 1991; Weise, Costa, & Kudela, 2006). Furthermore, the relationships we identified have been found in other mobile marine predators that utilize dynamic biophysical features (e.g., upwelling centers, fronts, and eddies) and have important implications for overlap with human use (Maxwell et al., 2013; Scales et al., 2014). Recent research that considers what dynamic habitat use means for ocean resource management has found that while place‐based static management approaches can be used to define general areas of overlap between protected species and human use, our ability to monitor and manage animals in relation to human activities will likely require management approaches that are also dynamic in space and time (Gregr, Lessard, & Harper, 2013; Lewison et al., 2015; Maxwell et al., 2015; Wedding et al., 2016). The predictive models used in this study offer an increased understanding of when and where potential conflicts may arise and reflect the importance of dynamic, spatially explicit conservation and management initiatives for many other marine top predators (Louzao et al., 2011).

Our broadscale modeling approach presented provides information on the population‐level habitat associations for females from the two largest California sea lion breeding colonies in the California Current System. However, some caveats must be considered, as foraging habits and at‐sea distribution can vary by season and life history stage. For lactating females, distribution may be closer to the rookery at different times of year, specifically during the breeding season (June–July). During this time, females give birth and their foraging trips are even more constrained by pup age. Our satellite observations did not include this period; therefore, we caution extrapolation of model predictions to other data‐limited times of the year. Second, separate species distribution models should be constructed for nonlactating females, as are free to disperse away from the rookery to exploit productive areas (Melin et al., 2000). Finally, habitat selection may be colony‐specific. While core residency for both colonies indicated proximity of lactating females to breeding grounds and some overlap in spatial distribution, San Miguel Island females were more likely to move north near Monterey Bay, whereas only San Nicolas Island individuals used the Southern California Bight (Figure 1b,c). This spatial foraging segregation may be a mechanism to reduce intraspecific competition between two breeding colonies (Kuhn & Costa, 2014; McHuron et al., 2016), but also may reflect spatial constraints associated with lactation, as San Nicolas is approximately 120 km southeast of San Miguel (Kuhn & Costa, 2014). Future work may include finer‐scale, colony‐specific, or behavior‐specific models that may capture habitat preference more relevant to specific colonies or individual behaviors. Additional tagging studies may help capture suitable habitat for all colonies within the Southern California Bight, including Santa Barbara Island and San Clemente Island. However, for the purposes of understanding the environmental drivers that influence distribution on a population level, this study represents an important first step in modeling spatially explicit habitat use for this species. Our findings demonstrate the utility of a marine species distribution model as a novel approach for identifying changes in central‐place forager habitat, with important implications at the population level. An increased understanding of habitat use can not only improve our ability to monitor and predict future shifts in distribution as a function environmental variability, but also serve in the context of species protection and fisheries management.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTION

DB, SF, KS, EH, SB, EM, and EB contributed to research design and data analysis. SM, EM, PR, and CK conducted fieldwork and collected the data. DC, SM, EH, SB, LC, and RL were responsible for funding application. DB wrote the first draft. All authors contributed to draft revisions and approved the final version of the manuscript.

Click here for additional data file.