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

Abstract 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 centralplace 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 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 Kuhn & Costa, 2014;Lowry et al., 2017;McClatchie et al., 2016;Melin et al., 2008;Trillmich et al., 1991).
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; T A B L E 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  Melin et al., 2000;Villegas-Amtmann et al., 2011, 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 centralplace 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 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-realtime 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.

| Data sets and tagging methodology
The All data processing and analyses were carried out in the R environment, version 3.3.0 (R Core Team, 2016).

| 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 mod- 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

| Remotely sensed oceanographic data
Remotely sensed environmental data were obtained for both sea lion and CRW tracks using Xtractomatic (Simmons, 2016). The data  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).
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

| 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.

| 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) (1) Weight = 2 * (distance track − distance CRW )∕distance track

| 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/m 3 ), 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).

| 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 ( (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 exam-  (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 (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.