Seasonal variation in environmental and behavioural drivers of annual‐cycle habitat selection in a nearshore seabird

Conservation of highly mobile species often requires identifying locations or time periods of elevated vulnerability. Since both extrinsic habitat conditions and intrinsic behavioural and energetic requirements contribute to habitat use at the landscape scale, identifying spatial or temporal foci for conservation intervention requires understanding how habitat needs and distributions vary across the annual cycle. Nearshore marine birds inhabit highly dynamic systems and have widely varying habitat needs among breeding, moult and non‐breeding seasons, making them a useful case study for testing the relative contributions of individual resource requirements and environmental conditions in driving annual variation in distribution patterns.


| INTRODUC TI ON
A frequent goal in wildlife conservation is to identify locations or time periods in which management actions are likely to have the greatest positive effect (Radeloff et al., 2013). For mobile organisms, achieving this objective requires first understanding how both extrinsic habitat conditions and intrinsic habitat requirements vary across space and time (Marra, Cohen, Loss, Rutter, & Tonra, 2015).
Individuals require different habitat components for activities that vary within or among different stages of the annual cycle, including obtaining food, sheltering from predators, thermoregulating, raising young and moving or migrating among habitat patches (Börger, Dalziel, & Fryxell, 2008;Morrison, Marcot, & Mannan, 2012).
Temporal changes in environmental conditions, habitat features and individual resource requirements drive habitat selection on daily (e.g., commuting vs. resting), seasonal (e.g., breeding vs. non-breeding) and annual or even decadal scales (Garthe & Hüppop, 2004).
Thus, designing effective conservation strategies requires not only incorporating background variation in biotic and abiotic habitat features, but also understanding how temporal variation in individual resource requirements contributes to observed patterns of occupancy and distribution (Small-Lorenz, Culp, Ryder, Will, & Marra, 2013).
Opportunities to conduct an individual-based analysis of habitat requirements have vastly improved due to miniaturization and refinement of animal-borne tracking devices (Hebblewhite & Haydon, 2010). To date, however, many habitat assessments derived from individual tracking data have considered each location as an independent occupancy record, rather than considering that observed locations are samples of continuous movement processes that vary with underlying behaviour (Tremblay et al., 2009). To better understand the behavioural processes that generate observed locations, advanced computational approaches such as state-space models (Jonsen, Myers, & James, 2007) and Hidden Markov Models (Patterson, Basson, Bravington, & Gunn, 2009) have been developed to detect the underlying patterns represented by a series of locations. Rather than treating each location in isolation, these modelling frameworks interpret locations as steps in a continuous movement process, allowing for the detection of changes in behaviour based on shifts in these latent movement processes. From a conservation standpoint, incorporating dynamic behavioural processes into habitat selection analysis is an important refinement of resource selection models, as it allows the identification of locations or time periods at which either environmental conditions or individual behaviour may result in increased vulnerability (e.g., Cureton & Deaton, 2012).
Marine birds provide a particularly interesting case study for assessing effects of temporal variation in individual resource requirements on habitat selection. Seabirds have widely varying habitat requirements between the breeding season, when birds are central-place foragers breeding on land (often offshore islands) and foraging at sea close to the colony site, and the non-breeding season, when birds range widely throughout marine habitats (Weimerskirch & Wilson, 2000). Within each stage of the breeding cycle, habitat use also depends on individual characteristics and condition (Bearhop et al., 2006), phenology (Catry, Ramos, LeCorre, & Phillips, 2009), colony size and location (Lamb, Satgé, & Jodice, 2017;Lewis, Sherratt, Hamer, & Wanless, 2001), and environmental features (Tew Kai et al., 2009). These factors all contribute to variation in individual energy requirements, resulting in differences in foraging strategies and habitat preferences (Daunt, Afanasyev, Silk, & Wanless, 2006;Phillips, Bearhop, McGill, & Dawson, 2009). Moreover, since marine environments are highly dynamic, habitat features present at any given location may vary widely in space and time. Nearshore habitats contain a particularly high diversity of landscape features, including estuaries, river plumes and upwellings, which support varying prey species assemblages (Becker & Beissinger, 2003), and are structured by different oceanographic and anthropogenic processes than large marine ecosystems (Gray, 1997). Among marine birds, factors driving habitat selection differ between pelagic seabirds, which inhabit offshore marine environments, and nearshore seabirds that utilize coastal systems (Thaxter et al., 2012). Many of the same underlying individual, colony-based and environmental factors that influence habitat choice in pelagic species also affect nearshore seabirds (e.g., Erwin, 1977;Suryan, Irons, & Benson, 2000); however, since most studies of habitat selection in marine birds have focused on pelagic species, the importance of nearshore habitat features to seabirds remains poorly understood.
One of the most visible nearshore seabirds in the neotropics and subtropics is the Brown Pelican (Pelecanus occidentalis), typically found within 20 km of the nearest shoreline (Shields, 2014). Although currently abundant, the species was reduced to near extinction by DDT exposure during the mid-twentieth century (McNease, Richard, & Joanen, 1992) and continues to experience high mortality rates during oil spills (USFWS, 2011), making it an important target of regional monitoring and mitigation efforts (Jodice, Adams, Lamb, Satgé, & Gleason, 2019). Despite this conservation interest, limited information is available on habitat requirements of Brown Pelicans outside of terrestrial breeding sites. The majority of the eastern subspecies P. o. carolinensis of Brown Pelicans breeds in the northern Gulf of Mexico (GOM), where their breeding habitat occurs within three administrative planning areas (BOEM, 2012). However, since boundaries of planning areas are based on offshore geomorphology rather than species ecology, it is unclear whether they represent appropriate units for management of mobile marine species. Indeed, data from band recoveries suggest that Brown Pelicans may move widely during non-breeding (Schreiber & Mock, 1988;Stefan, 2008), meaning that current administrative boundaries may not encompass year-round habitat needs for breeding pelicans. Previous studies of fine-scale habitat use and movement patterns of Brown Pelicans in the GOM (King et al., 2013;Walter, Leberg, Dindo, & Karubian, 2014) have been limited to pelicans captured in a small portion of the central planning area. An integrated understanding of yearround movements of the species throughout the GOM would allow for more detailed analyses of spatial patterns, including identifying locations or time periods of critical habitat, quantifying vulnerability to stressors and evaluating whether current administrative boundaries effectively encompass the habitat needs of breeding populations within those regions.
Based on 3 years of tracking data from Brown Pelicans breeding in the northern GOM, we characterized year-round habitat preferences. To ensure that locations used to determine habitat preferences accurately reflected areas where individuals were interacting with marine habitat features, we filtered points based on modelled behavioural states, retaining only points that represented resident behaviour. We then evaluated habitat selection across a suite of environmental variables. Our results provide a baseline understanding of how nearshore habitat features influence patterns of use and occurrence in marine birds, and how seasonal variation in individual behaviours and resource requirements drives habitat selection on an annual scale.

| Study area
From 2013 to 2015, we collected data on breeding adult and nestling pelicans as part of a region-wide study from six colonies in the northern GOM between 83° and 98° W and 27° and 31°N ( Figure 1).
All colony sites were within the same marine ecoregion (Spalding et al., 2007), although local environmental conditions and breeding population size varied among sites (Lamb et al., 2017). Individuals of both sexes were captured at all colonies, although sex ratios varied between sites since sex could not be determined prior to capture and therefore could not be included in the experimental design (Table 1).

| Pelican locations
To track movement patterns of adult pelicans, we used 65 g solar GPS Platform Terminal and Cellular Terminal transmitters (GeoTrak Inc; NorthStar Science and Technology) with a backpack-style Teflon ribbon harness attachment (Dunstan, 1972). To elevate the transmitters and prevent feathers from covering the solar panels and antenna, we mounted each device on a 6 mm thick neoprene pad that also extended 6 mm beyond the perimeter of the transmitter in all directions. Transmitters were programmed to collect 12 fixes/day during breeding (April-August; every 90 min from 1030 to 0130 GMT), 10 fixes/day during pre-and post-breeding (September-October and February-March; every 90 min from 0700 to 0100 GMT) and 8 fixes/day during winter (November-January; every 120 min from 0700 to 0100 GMT). We obtained an average error estimate for GPS points from transmitters at known locations (N = 220) of 4.03 ± 2.79 m.
Over 3 years, we fitted 85 individual pelicans with GPS transmitters, 67 of which recorded sufficient data for subsequent analysis (i.e., transmitted over at least 30 days; Table 1). Adults were captured at nests using leg nooses in either the late incubation or early chick-rearing stage of breeding. All captured adults were weighed, measured, banded and sampled for blood and feathers. We also calculated adult body condition index (BCI) as the residual of the linear relationship between culmen length and mass (Eggert, Jodice, & O'Reilly, 2010). Since morphology is not always sufficient to determine sex in Brown Pelicans, adults were later sexed via PCR using DNA obtained from blood samples (Itoh et al., 2001). Total handling time from capture to release averaged 19 min (±6.5 min standard deviation). To avoid subsequent influence of disturbance on observed movement patterns, we avoided entering or trapping in areas of the colony where we had already captured nesting pelicans. Since individual characteristics (e.g., sex, size, physical condition) may influence pelican foraging movements during breeding (Walter et al., 2014), we used nested ANOVAs to determine whether individual characteristics of tracked adults differed systematically among planning areas or between colonies within each planning area. Unless otherwise specified, all statistical manipulation of spatial data was conducted using the "adehabitat" family of packages (Calenge, 2006) in R 3.2.3 (R Core Team, 2014).
For each individual, we removed all locations between when the individual was first released and when it returned to its nest site (typically 2-4 hr after release). We manually identified and removed outliers using a speed cut-off of 65 km/hr between successive points, which is the maximum travel speed recorded for Brown Pelicans (Schnell & Hellack, 1978). Cleaned locations for each individual were then interpolated to regular 90-min intervals. Since location data were not collected overnight, we chose not to interpolate tracks between successive days, and we differentiated each day as a separate trajectory by cutting tracks between each set of two successive points separated by a gap of greater than 6 hr.

| Movement states
To distinguish resident from commuting behaviour, we fit a twostate Hidden Markov Model (HMM; Patterson et al., 2009) to the regularized movement trajectories using the moveHMM R package (Michelot, Langrock, Patterson, & Rexstad, 2015). Hidden Markov Models are a particularly flexible and efficient way of characterizing behavioural states from precise and regularized tracking data (Langrock et al., 2012), and thus, are a good fit for GPS tracking locations. Briefly, the model assumes a priori that observed movement data are driven by underlying movement "states," characterized by a distribution of step lengths (distance between successive points) and turning angles. Each point therefore represents a sample drawn from these underlying distributions and can be assigned to one or the other based on probability. A Markov chain is used to describe the state parameters and classify data according to its most probable state membership.
To reduce spatial autocorrelation and characterize patterns of movement between rather than within days, we fit the model to a reduced data set of one location per day, calculated as the centroid of all locations for that day. We began with the assumption that local movement would be characterized by short step lengths and sharp turning angles (State 1) and commuting movement by long step lengths and wide turning angles (State 2). Therefore, initial step lengths were set at 5 (±5) km for State 1 and 10 (±10) km for State 2.
Initial turn angles were set at π radians for state 1 and 0 radians for State 2. Angle concentration for each state was initially set at 1. In subsequent analyses, all points along the trajectory for a given day were assigned to the movement state associated with that day.
We further assigned resident locations to one of three stagesbreeding, migration or wintering-based on their positions within the annual cycle. We defined the breeding stage as the period during which individuals returned to the colony site at least once every three days (i.e., successive locations at the colony site were separated by no more than 72 hr). We considered the winter stage to begin when an individual arrived at the end of its final period of commuting movements. All movements between the end of the breeding stage and the beginning of the wintering stage were classified as migratory. Individuals that did not migrate (i.e., whose breeding and wintering ranges overlapped) progressed directly to the wintering stage at the conclusion of breeding.

| Environmental variables
We measured environmental characteristics of resident locations using seven habitat variables, four of which were fixed for any given point (distance to coastline, distance to river outflow, bathymetry and bottom substrate) and three of which we measured monthly (net primary production, sea surface salinity and sea surface temperature; Table 2). We chose these variables to represent a suite of likely drivers of nearshore habitat variation, particularly the distribution of pelican prey which tends to include Gulf Menhaden (Brevoortia tyrannus) and other schooling fish (e.g., Deegan, 1990). Since limited data are available on fine-scale variation in oceanographic features such as currents and eddies, and since these features have a high degree of short-term variability in coastal areas (Kaltenberg, Emmett, & Benoit-Bird, 2010), we used the distance to physical features that influence the movement of water in coastal areas (coastline, major river outflow) as proxies for these processes. Following the North American Rivers and Lakes map, we defined major rivers as orders 4 and higher in the Horton classification system. Depth and bottom substrate can influence both prey distributions and oceanographic characteristics. Net primary production uses remotely sensed ocean colour and its known relationships to primary production throughout the water column to calculate daily phytoplankton fixation of chlorophyll per unit of ocean surface (see Behrenfeld & Falkowski, 1997 for full details of these calculations). Compared to sea surface chlorophyll concentrations, net primary production provides a more robust index of total oceanographic productivity, which in turn influences the distribution of consumers at higher trophic levels. Salinity and temperature also influence the distribution of aquatic prey species depending on their osmotic and thermal tolerances. Since some data were reported at finer spatial resolutions than others (Table 2), we standardized all variables to a resolution of 0.1-degree grid squares (approximately 100 km 2 ). Distance values were calculated as the distance from the grid square centroid to the feature of interest. For all other variables, we resampled the data using the mean value for each 0.1-degree grid square.

| Habitat suitability and distribution
We mapped preferred habitat characteristics in ecological space using a multivariate ordination of all habitat variables using a Hill-Smith principal components analysis (Hill & Smith, 1976), which allows the inclusion of both categorical and continuous variables.
For each grid square, we calculated habitat suitability as the squared Mahalanobis distance of that point from optimal location TA B L E 2 Environmental data layers used for habitat analysis To characterize individual responses to the measured habitat variables, we used an Outlying Mean Index (OMI) analysis (Dolédec, Chessel, & Gimaret-Carpentier, 2000). Briefly, OMI is an ordination technique that characterizes available sites based on a set of environmental variables, sets the mean of all conditions at zero in n-dimensional space then determines the axis that describes the maximum amount of marginality (difference from the mean: i.e., similar to residuals) of individual animals or species in ecological space. Thus, the first axis of the OMI is the combination of environmental characteristics that best explains the position of animals across available resources. Similarly, the position of each habitat characteristic on the first axis of the OMI represents that variable's contribution to animal distributions; that is, the strength of selection on that characteristic. Outlying Mean Index does not assume specific resource selection functions and allows differences in individual niche selection to be taken into account when describing the distribution of a group of animals. We conducted OMIs for each month on all individuals and habitat variables for each behavioural state, then averaged the scores of individuals on the first OMI axis to calculate niche location and breadth for groups within the population. We considered scores between −1 and 1 to represent relatively weak associations, while scores <−1 or >1 represented relatively strong associations. We also examined spatial distribution of breeders from different administrative planning areas (western, central and eastern planning areas for the Bureau of Ocean Energy Management; BOEM) that contained differing levels of oil and gas activity in marine waters. We calculated 95% kernel density estimates for all individuals from each breeding area using the "ks" package (Duong, 2015) with a plugin bandwidth estimator (Gitzen, Millspaugh, & Kernohan, 2006;Wand & Jones, 1994).
We then used an Albers Equal-area Conic Projection to calculate the areas included within each planning area's 95% kernel contour and to estimate the intersection areas between kernels from different planning areas.

| Pelican locations
After cleaning and interpolating all collected locations (N = ca. Overall, 61.5% of bird-days were classified as resident and 38.5% as transient (Appendix S1 in Supporting Information, Figure S1.2).

| Habitat suitability and distribution
Resident locations were concentrated in shallow nearshore waters (Table 3). The habitat variables most strongly associated with resident locations year-round were net primary production (positive) and sea surface salinity (negative) (Figure 2). Sea surface temperature was negatively associated with residency during non-breeding, but the association diminished to near zero during the breeding season. Compared to dynamic variables, fixed factors were less strongly associated and less variable in their relationship to pelican habitat use and did not vary during the year. Bathymetry had a positive relationship with residency (i.e., pelicans were more likely to occupy shallower waters), while distance to coastline and distance to river outflow were both negatively associated with use by pelicans (i.e., pelicans preferentially occupied habitats closer to coastlines and river mouths). Overall, areas of highest habitat suitability yearround were located in the northern Gulf, particularly the central and western planning areas (Figure 3). During summer, the total area F I G U R E 3 Suitability scores of available habitat for Brown Pelicans in the Gulf of Mexico based on Mahalanobis distances, 2013-2015 (Albers Equal-Area Conic projection). Darker colours indicate higher suitability of preferred habitat was narrowly restricted to coastal areas of the northern Gulf; however, during the fall and winter, suitable habitat characteristics also occurred from the nearshore region out to ca. 600-km offshore.

| Seasonal habitat associations
The average values of dynamic habitat characteristics in occupied habitats were generally similar between the breeding and staging periods, but differed significantly during winter (Figure 4). In general, pelicans utilized lower-productivity, higher-salinity and lowertemperature waters during winter compared to either breeding or migration. Selection for high productivity and low salinity relative to mean values peaked during breeding and migration periods, while selection for low temperatures was highest during winter ( Figure 2).

| Regional variation and overlap
Patterns of association with dynamic habitat variables varied be-

| D ISCUSS I ON
Our study provides insights into the habitat associations of an abundant and highly visible nearshore seabird, the Brown Pelican, with fine-scale environmental characteristics in the northern GOM. We documented strong associations between pelicans and several dynamic habitat features, particularly primary production and salinity. These associations peaked during the early chick-rearing and moult stages of the annual cycle, suggesting that highly productive Although extensive work has described the environmental factors driving at-sea habitat use by seabirds in pelagic waters (e.g., Haney, 1985;Pinaud & Weimerskirch, 2005;Tew Kai et al., 2009), relatively little is known about the factors driving marine habitat use in nearshore seabirds. For the most part, prior studies of habitat preferences in nearshore-foraging species have been conducted in northern temperate waters (e.g., Becker & Beissinger, 2003;Day, Nigro, & Prichard, 2000;McLeay et al., 2010;Yen, Sydeman, Bograd, & Hyrenbach, 2006). Similar to results from these systems, we found that marine productivity was the most significant driver of habitat selection of Brown Pelicans in the GOM. Also, in concordance with previous results (Becker & Beissinger, 2003;Day et al., 2000), we found that the influence of sea surface temperature on at-sea distribution was significant but highly variable over time. In a departure from previous assessments of habitat use of nearshore seabirds, which generally found little effect of salinity on habitat use, we found that salinity strongly influenced habitat use for Brown Pelicans. Although the effects of salinity have not been extensively documented on this species or other coastal seabirds, recent studies (e.g., Zamon, Phillips, & Guy, 2014) have suggested that river plumes can be important nearshore-foraging habitat for seabirds, concentrating prey in a manner analogous to oceanic fronts in pelagic systems. While distance to river outflow was only weakly related to pelican habitat suitability in this study, pelicans were often located in relatively large estuarine complexes and therefore may ultimately be responding to salinity gradients that exist even at comparatively greater distances from river mouths.
Since the scale of movement that we observed was relatively small (on the order of tens of kilometres per day, rather than hundreds of kilometres as is commonly observed in pelagic seabirds), we chose environmental variables likely to relate to the distribution of prey rather than those that might facilitate long-distance movement (e.g., prevailing winds) or visual identification of foraging areas (e.g., ocean colour). The influence of salinity in particular is correlated with the abundance and distribution of prey items. Brown Pelicans in the GOM forage primarily on Gulf Menhaden (Lamb et al., 2017), which concentrate during the spring and summer in low-salinity estuarine environments (Deegan, 1990). Both summer and winter distribution of preferred pelican habitat corresponded closely with Gulf Menhaden distributions, indicating that pelicans select habitat principally as a function of prey concentrations. We did not find that the spatially fixed metrics we tested (distance to coastline, distance to river outflow, bathymetry or bottom substrate) had a strong influence on habitat suitability. Previous studies (e.g., Suryan, Santora, & Sydeman, 2012) have suggested that such metrics tend to provide a more consistent predictor of the distribution of seabirds than sea- 3) of monthly averages for each variable, and grey areas are 95% confidence intervals of regression lines. Coloured bars along the x-axis indicate the approximate timing and duration of breeding (red), wintering (blue) and migratory staging (green) periods environmental data available (10 km 2 ) and the temporal resolution of the GPS data collected (90-min intervals) did not allow us to distinguish fine-scale foraging areas from commuting or resting habitat. Thus, we focused on mesoscale movement patterns and habitat selection on a monthly timescale. The fact that seasonally varying parameters were more strongly related to habitat selection than physical oceanographic features is consistent with previous observations that mesoscale habitat use is likely to be driven by primary productivity, while physical features become more important at the micro (<10 km) scale (Becker & Beissinger, 2003). Habitat selection likely also occurs at finer temporal scales than those described by this study (Kristan, 2006) and may vary with daily or weekly changes in estuarine dynamics that alter distribution and concentrations of prey.

| Seasonal variation
We found that pelicans showed a particular preference for high-productivity, low-salinity habitats during the energetically demanding breeding and post-breeding moult periods. Correspondingly, the availability of suitable habitat was much lower during the summer (May-July) and early fall (August-September), which are peak times for chick rearing and post-breeding dispersal respectively, than during the remainder of the year. Moreover, the habitats that contained these preferred features tended to be located in areas of high anthropogenic activity. Preferred foraging areas included highly productive coastal estuaries, which were often located near shipping hubs and ports in the comparatively cooler waters of the northern Gulf. During moult, breeding pelicans from throughout the northern Gulf used habitats in and around the Mississippi Delta ( Figure 6, inset), which is also a hot spot for shipping, oil extraction and concentrations of agricultural and industrial contaminants from throughout central North America.
The early chick-rearing and moult periods therefore may represent hot moments (McClain et al., 2003), in which a confluence of habitat requirements and availability may increase exposure to risk factors for breeding Brown Pelicans; however, further analysis is required to assess the potential extent and consequences of this elevated exposure. Although our study focuses on Brown F I G U R E 6 Annual 95% kernel density estimates for locations of Brown Pelicans originally captured at breeding colonies in the eastern (blue), central (orange), and western (green) planning areas, 2013-2015 (Albers Equal-Area Conic projection). Areas used by breeders from two or more planning areas are shaded in purple, and areas used by breeders from all planning areas are shaded in red. The inset map shows only the area of overlap between breeders from all three planning areas Pelicans, the variation in both habitat features and annual-cycle energetic needs that drive pelican occupancy is likely to be experienced by a suite of nearshore waterbird species, which also rely on prey-rich nearshore habitats and experience increased energy demands during chick-rearing and moult.

| Annual-cycle movements and spatial population structure
In addition to biophysical habitat associations, we examined specific To date, studies of Brown Pelican non-breeding movements have been limited to information on band recoveries, typically from birds banded as juveniles (Schreiber & Mock, 1988;Stefan, 2008) and tracking data from individuals captured during non-breeding (King et al., 2013;Poli, 2015). This has limited the possibility of linking non-breeding birds to breeding colonies outside the breeding season. Our study incorporates individual data on year-round movements of pelicans from known breeding locations. Measuring overlap between breeding populations in different planning areas of the GOM used for large-scale marine spatial planning helps to refine current understanding of the distribution of environmental risk among management units and to better identify which segments of the overall breeding population are affected by spatially explicit threats in the marine environment. Moreover, our results highlight the wide-ranging movements of marine vertebrates, which regularly travel between administrative planning areas across the northern GOM. For these species, an integrated approach to risk assessment and restoration activities may be needed in order to effectively manage the regional population.

ACK N OWLED G EM ENTS
Funding for this study was provided by the Bureau of Energy

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this study are available on Movebank (movebank.org, Movebank Study ID 296027617) and are published in the Movebank Data Repository with https ://doi.org/10.5441/001/1.7856r086.