Study area
We conducted our study in the South Atlantic Bight, USA, which extends from the Cape Fear River Basin to approximately Cape Canaveral (Fig. 1). The coast here is characterized by a complex geomorphology of barrier islands, estuaries, and salt marshes. The area supports ca. 15 brown pelican colonies annually (active breeding from April – September) and many of the beaches and islands are used as migratory stopover, staging, or wintering grounds for this species and others [33].
Satellite transmitter deployments
Nesting pelicans were outfitted with GPS satellite transmitters (GeoTrak Inc., North Carolina, USA) at four colonies in coastal South Carolina (Bird Key Stono, 32° 38′ N, 79° 58′ W, n = 21; Castle Pinckney, 32° 46′ N, 79° 54′ W, n = 12; Marsh Island, 32° 59′ N, 79° 33′ W, n = 7; Deveaux Bank, 32° 32′ N, 80° 10′ W, n = 5). Colony size ranged from ca. 50–2000 pairs. Deployments commenced during the chick-rearing stage (May–July) of the 2017 and 2018 breeding seasons. Transmitters weighed ~ 65 g (10 × 3.5 × 3 cm) and constituted ≤3% body mass of instrumented individuals (range = 2475–4350 g), the recommended threshold for large seabirds [34]. Briefly, nest-attending adults were captured via either neck or leg noose and equipped with a solar GPS Platform Terminal Transmitter dorsally using a backpack-style harness system. For a description of specific attachment procedures, see [35]. During the post-breeding stage of deployment (September – November), units were programmed to record 10 locations per day at 90 min intervals between the hours of 01:00–23:30 GMT and were duty-cycled on an 8 h on to 36 h off activity schedule. Unit error was assumed to be similar to that of [32], i.e. 4.03 ± 2.79 m.
Hurricane events
Our opportunistic analysis of pelican movement in relation to hurricane activity includes three storm events. On 10 September 2017, Hurricane Irma made landfall in southwestern Florida, USA, as a Category 4 tropical cyclone. Over the subsequent 1.5 days, Irma proceeded north along the coast of western Florida before weakening and degenerating near the central Georgia-Alabama border. Although the storm was centered mainly along the Gulf coast of Florida, much of the southeastern Atlantic seaboard was affected by the outer cyclonic bands (Fig. 1).
Hurricane Florence made landfall on 14 September 2018 in southern North Carolina, USA, as a reduced Category 1 tropical cyclone, having been a Category 4 cyclone 4 days prior. Florence tracked inland in a southeasterly direction as it weakened, degenerating over West Virginia, USA, three days after landfall, affecting predominantly the coastal Carolinas (Fig. 1).
Less than 1 month later, Hurricane Michael made landfall in the panhandle of Florida on 10 October 2018 as a Category 4 tropical cyclone. Michael followed a northeasterly trajectory after landfall, weakening incrementally over the southeastern United States before restructuring as an extratropical cyclone 2 days later off the Mid-Atlantic coast (Fig. 1). Similar to Irma, Michael impacted much of the Atlantic seaboard due to the trajectory, strength, and spatial extent of the storm.
Meteorological data
A kernel density analysis was used to identify the core spatial area utilized by instrumented pelicans during each hurricane event. Subsequent utilization distributions (UDs) were used to determine a representative location for assessing pelican response to meteorological indices. This approach allowed for the acquisition of meteorological data that would represent shared conditions for the greatest number of individuals throughout the tracking period. We used only locations recorded during the calendar month of the respective hurricane event, which corresponded with peak cyclonic activity but limited seasonal changes in weather. Distributions therefore reflected core use areas during the entire passage of the cyclone as well as the remainder of the month in which the cyclone occurred. Erroneous locations were identified and removed through a combination of visual inspection (e.g. consecutive locations separated by unrealistic distances) and a speed filter of ≥65 km per hour [36]. Kernel bandwidth was determined using R statistical software (v 3.4.2.) through a plug-in bandwidth selector in package ks [37]. Locations within the 25% UD (i.e. core range) identified in the kernel density output during the month of each respective hurricane (grid = 400, extent = 0.4°) were then used to assess movement patterns in relation to storm events. Roughly, the area of highest use by pelicans during these time periods paralleled the coastline from central South Carolina to north-central Georgia (Fig. 1). Individual pelicans located outside of the prior 25% UD at the time of hurricane passage (e.g. in Chesapeake Bay) were manually excluded from further analysis, as well as individuals for whom movement data was not complete for the entire time period.
Meteorological data were obtained via the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information from the Hunter U.S. Army Airfield, Savannah, Georgia (station 74780413824), to represent conditions experienced during Hurricane Irma, and from the Marine Corps Air Station Beaufort, Beaufort, South Carolina (station 72208593831), to represent conditions during Hurricanes Florence and Michael (https://www.ncdc.noaa.gov/). These sites were within the 25% UD in the kernel density analysis. Although spatially similar, multiple weather locations were required as neither station had complete data for all three hurricane events in totality. Meteorological data were collected hourly and spanned the entire month of each cyclonic event. Data were requested 04 November 2017, 28 November 2018, and 12 December 2018, respectively.
Behavioral clustering
We used an Expectation Maximization binary Clustering (EMbC) algorithm to derive biologically-relevant behavioral states for individual brown pelicans [38]. EMbC uses unsupervised relationships between successive locations incorporating path distance and tortuosity (i.e. velocity and turning angle) to infer underlying behavioral processes. EMbC is particularly appropriate for remotely-sensed location data as it accounts for spatial and temporal correlations and uncertainties in the input features and is robust to spatial data collected at relatively long intervals [39]. Critically, EMbC is capable of producing biologically-relevant classifications for locational data recorded at timescales relevant to the current study (e.g. [40]). Each point within individual tracks was clustered into one of four categories: low velocity/ low turning angle (LL), low velocity/ high turning angle (LH), high velocity/ low turning angle (HL), and high velocity/ high turning angle (HH) (Fig. 2). These four behavioral nodes were biologically interpreted as corresponding to inactive, localized search, commuting, and dispersive search behaviors, respectively. Following [38], a post-processing smoothing procedure was applied based on consecutive behavioral correlations to manage temporally-irregular data. This smoothing procedure searches for clusters of the same behavioral assignment that contain a single point of a different classification, and adds additional likelihood weight to that single point belonging to the larger cluster, a feature explicitly implemented in state-space models. In this way, the smoothing procedure favors homogenized bouts of behavior instead of single-point behavioral switches during clusters of equal assignment. We also calculated mean step length (distance between successive points) and net displacement (maximum distance from the first location in the series) for descriptive purposes. Each point was finally matched temporally to the closest hourly meteorological variable for statistical analysis.
Statistical analyses
We assessed the effects of meteorological drivers on pelican behavioral state with multinomial logistic regression following [41]. To simplify model interpretation and to examine activity patterns more accurately matched to the temporal resolution of the data, models were conducted on a reduced set of two behavioral nodes classified as either active (including localized search, commuting, and dispersive search; LH, HL, and HH, respectively) or inactive (LL). Environmental variables of interest (barometric pressure and wind velocity) were chosen a priori based on data completeness, relevance to cyclonic activity, and probability of being sensed by individual pelicans [14].
Both tracking and meteorological data were further subset to exclude other potentially confounding anomalous conditions. We defined an anomalous event as a barometric pressure reading ≥1 SD from the monthly mean. Only data collected from the end of the last pressure anomaly pre-cyclone to the first pressure anomaly post-cyclone were therefore included in our regression analysis, thus creating a temporal segment of activity that was exclusively characterized by ‘baseline’ conditions with the exception of the cyclonic event. Significant differences of barometric pressure and wind velocity between study periods were assessed via Kruskal-Wallis chi-squared tests, with Wilcoxon rank sum tests used when significant differences were found.
Four multinomial logistic regression models were fit to the data using R package mlogit [42], including a null model, single-effect wind velocity model, single-effect barometric pressure model, and global model including both wind velocity and barometric pressure. Model selection was performed within each set using Akaike’s Information Criterion (AIC), with the best-performing model indicated by the lowest AIC value. Given low AIC similarity between models, we did not model average. Environmental variables were interpreted as having a significant effect on individual behavioral states at p < 0.05. We further assessed transition probabilities using the top-performing model, with the null state (i.e., reference level) defined as inactive (i.e., the probabilities are reflective of transitioning from inactivity to activity).