Tag deployment
From February to May of 2012, 2013, 2016, and 2017, we deployed 16 satellite-linked transmitting tags onto humpback whales that commenced migration in nearshore waters around the WAP. These animals were from Breeding Stock G, which breeds off the western coast of South and Central America [41]. Sirtrack tags and Wildlife Computers (Redmond, WA, USA) SPOT5, SPOT 6, and MARK 10 Platform Transmitting Terminals (PTTs) were utilized, and tagging was limited to animals > 12 m in length. Each tag was contained in a sterilized housing prior to deployment and was anchored in the muscle near the dorsal with stainless steel barbs, with the transmitting antenna remaining free outside of the animal [7]. The tags were designed to implant up to a maximum of 290 mm into the back of the whale. Tags were deployed from a range of 3–10 m from a Zodiac Mark V or a Solas ridged-hulled inflatable boat using an ARTS Whale Tagging pneumatic line thrower compressed air system [42].
Satellite transmissions were activated via a salt-water switch, and locations of the whales were obtained through the Argos System of polar-orbiting satellites (Argos, 1990). Tags were programmed to transmit during specific hours and days. Since the tags were also being utilized for other year-specific projects, duty cycling varied across years. In 2012, tags were programmed to transmit between 00:00–04:00 and 12:00–16:00 GMT. In 2013, tags were programmed to duty cycle 3 h on, 3 h off, except for Sirtrack tags (identified by PTT IDs starting with 113), which duty-cycled at 6 h on/6 h off. In 2016, some tags were programmed to transmit continuously and three were programmed to duty cycle at 1 day on, 4 days off. Tags deployed in 2017 were programmed to duty cycle 12 h on, 12 h off.
Demographic information
Skin and blubber biopsy samples were obtained from tagged whales whenever possible using standardized remote biopsy techniques [43]. Samples were obtained from the upper flank below the dorsal fin [44]. Blubber samples were used to provide life history and demographic information as covariates in models assessing migratory behavior. To determine the sex of biopsied whales, genomic DNA was extracted from these samples using a proteinase K digestion followed by a standard phenol–chloroform extraction method [45]. To assign pregnancy within sampled females, progesterone, a lipophilic steroid hormone, was quantified from a sub-sample of blubber using a progesterone enzyme immunoassay [46]. Pregnancy was then assigned by comparing the measured progesterone concentrations across a pre-validated binary logistic model developed from humpbacks of known pregnancy status sampled in the Gulf of Maine [46].
Data processing
The Argos data were Kalman processed. R (version 3.4.3, [64]) was used to filter raw observations from the satellite tags to remove points without location data, points with Argos error quality class Z (invalid location), and points with duplicate timestamps. Maps of the animals’ tracks were plotted using ggmap [47] in R [64].
Whales were determined to be migrating when they started a northward journey from the WAP without any lasting return movements. The date of departure for each whale was determined visually by graphing latitude as a function of Julian day. Trends in departure date by sex and reproductive status were assessed by creating multiple graphs of date of departure grouped by year, sex, reproductive status, year and sex combined, and year and reproductive status combined.
To determine rates of migration, speeds on the migratory route were calculated with data corrected for location error with a simple default HSSM with a 12-h timestep fitted in R using BSAM [64, 65]. Rate was the distance of the linear vector between 12 h timestep locations. Distances between locations were calculated using the function distanceTrack from the Argosfilter package [64, 66]. Individual’s mean speeds were calculated as the mean of all 12-h timestep rates for each animal. Great circle distance and speed were also calculated to allow comparison to more studies.
There were several locations where the tracks converged and allowed for a logical division of the migration corridor into three spatial sections, “WAP-Cape Horn (Drake passage)”, “Cape Horn (Chile)—Peninsula de Paracas (Peru)”, and “Peninsula de Paracas (Peru)—Zona Reserverda Illescas (Peru)”. Since not all tags transmitted for the entire migratory journey, these three discrete spatial sections allowed for a more valid estimation and comparison of speeds in some sections along the journey. Mean migratory speed was calculated for each section, as well as for the calving area. As humpback whales leave the Western Antarctic Peninsula at different times, a simple linear regression was performed using Julian day of the departure date (predictor variable) and speed (response variable) to investigate whether the timing of migration affected the speed at which the animal migrates. Because very few tags transmitted to completion of migration, we chose to look at speed in the first migratory section from the WAP to Cape Horn (latitude = − 55.9833). All data north of − 55.9833, as well as all animals that did not reach the cape, were filtered out, and the mean speed over the section was calculated for each remaining individual. To correct for issues of heteroskedasticity, speed was transformed with a log function, and the residual plot was assessed for any obvious signs of nonlinearity and heteroskedasticity. A QQ plot was used to check for the normality of residuals, and the data were tested for influential data points. To determine whether sex and reproductive status had an impact on speed, two Welch’s ANOVA tests were performed on the same speed data, using sex (male/female) as the predictor variable in the first test, and sex/reproductive status as the predictor variable in the second test (male, female-pregnant, female-not pregnant). For all tests, we report P values in the context of levels of support for the outcome rather than as a binary threshold of significance, with P-values < 0.01 indicating strong support, p-values between 0.01 and 0.1 offering suggestive, but inconclusive support, and p-values > 0.1 indicating no support [48, 49].
As coastal nations have exclusive sovereign rights for the purpose of conserving and managing marine species within the bounds of their jurisdiction [50], the amount of time the whales spent within Exclusive Economic Zone (EEZ) boundaries was calculated by summing the number of regular timestep observations from the BSAM model within each country’s national waters. While the satellite tags themselves did not collect data with great regularity, the BSAM model provides unobserved estimated locations along regular time intervals from available data, and these intervals were utilized for EEZ analysis.
Discrete behavioral modes were determined with hierarchical Bayesian state–space movement models manually constructed in JAGS. This was a departure from the simpler models used to assess locations, as it allowed for differences in transition probabilities and movement norms associated with behavioral states depending on whether the animals were in the foraging grounds, calving grounds, or migratory route. This model-associated spatial patterns of animal movement with predicted behavioral states while accounting for the significant error inherent in Argos Satellite location data. Unlike the simpler BSAM HSSM, this model did not return the coordinates of calculated unobserved locations to the user, necessitating the creation of the BSAM model to determine the earlier mentioned variable of speed.
We used a discrete-time dynamic correlated random walk model following Jonsen et al. [52] and Bestley et al. [51], where each movement stemmed from either a ‘traveling’ or ‘area-restricted search’ (ARS) state [51, 52]. When humpback whales encounter sufficient prey areas, they often engage in ARS by decreasing their travel speeds and increasing their turning angle radius and frequency; consequently, ARS behavior is defined as shorter step lengths with larger and more variable turning angles. The terminology ARS is used instead of foraging, as whales may also be engaging in other behaviors such as resting and calving in this state and our measurements are not based off a direct measure of feeding but rather use movement metrics. In humpback whales this spatial signature may persist for up to several days in one location [7, 53]. The traveling state, which is thought to occur when the animals are either actively migrating or located in habitats unsuitable for foraging, is characterized by fast travel rates and infrequent and small turning angles; in a state–space model this behavior is recognized by the presence of long step lengths with small and infrequent turning angle radius.
The first component of the state–space model was the process model, which estimates animal behavior with a first-difference correlated random walk [52]. The process model took the form:
$$d_{t} \sim N_{{2}} [\gamma_{{{\text{bt}}}} T(\theta_{{{\text{bt}}}} )d_{{t{-}1}} ,\Sigma ],$$
where dt is the difference between unobserved locations and coordinate vectors xt and xt-1 and N2 is a bivariate normal distribution with covariance matrix Σ, where \({\sigma }_{lon}^{2}\) is the process variance in longitude, \({\sigma }_{lat}^{2}\) is the process variance in latitude, and \(\rho\) is the correlation coefficient. γ is the autocorrelation of direction and speed between consecutive locations, with a value of between 0 and 1 (γ = 0 would signal a simple random walk). bt is an index used to denote behavioral mode, e.g. ARS or traveling. T(θ) is the transition matrix with mean turning angle θ which provides the rotation required to move between dt and dt-1:
$$T\left(\theta \right)=\left(\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right),$$
$$\Sigma = \left(\begin{array}{cc}{\sigma }_{lon}^{2}& \rho {\sigma }_{lon}{\sigma }_{lat}\\ \rho {\sigma }_{lat}{\sigma }_{lon}& {\sigma }_{lat}^{2}\end{array}\right).$$
This model is considered a switching model in the vein of Jonsen et al. [52], and a separate process model was run for each of the two behavioral states. As we are including two behavioral states, there were four possible transitions, two of which are calculated: α1, the probability of remaining traveling at time t if traveling at time t-1, and α2, the probability of traveling at time t given foraging at time t-1.
The second component of the state space model was the measurement equation or observation model. This equation calculated the temporally regular unobserved locations of the animals needed for the process equation from the error-prone and temporally irregular Argos location observations:
$$y_{t,i} = \, \left( {{1 }{-} \, j_{i} } \right)x_{{t \, {-}1}} + \, j_{i} x_{t} + \varepsilon_{t} ,$$
where i is an index for locations between times t and t + 1, and ji represents the proportion of the timestep at which the ith observation is made. Xt is the unobserved location of the animal at time t, yt,i is the ith observed position during the regular time interval t-1 to t, and εt is a random variable representing the error in the Argos locations. The variance in Argos observations was fixed for each Argos class error as demonstrated in Jonsen et al. [52]. Various classes of Argos errors are strongly non-gaussian, and are thus traditionally calculated with t distributions [52]. However, this can make the model so computationally complex that it cannot converge. This occurred with our model, and to counter this we ran the observation model with a multivariate normal distribution as done in Weinstein et al. [6, 7] and used the package Argosfilter in R [64, 66] to filter out implausible points that indicated speeds higher than ~ 20 km hr-1 (vmax = 6) [6, 7]. We used a timestep of 12 h, which we deemed to be a conservative balance between taking into account gaps in the data as well as ensuring behaviors did not change between locations. Although only two behavioral states were modeled, the means of the Markov chain Monte Carlo chains samples provided continuous values from 1–2. A mean behavioral mode of < 1.25 was considered traveling, whereas a value > 1.75 represented Area-Restricted Search. Estimations between 1.25 and 1.75 were treated as unknown [54].
To help address the inconsistent transmitting nature of the tags as much as possible, a joint estimation, in which estimation of behavioral states is conducted jointly across multiple animal movement datasets rather than individuals, was done. This method assumes that movement parameters may differ among individuals but are drawn from the same set of distributions, and allows the model to estimate parameters and state variables with greater precision by assuming a general range in value for all animals to borrow strength across multiple datasets, thus filling in for any animals with suboptimal data [55].
Priors for γ and θ were set to reflect the assumptions that the travelling state would have greater autocorrelation and lower mean turning angles than the ARS state. To allow for variance in transition probability between ARS and Traveling, as well as variance in behavioral state parameters as the animals switched from feeding, to migratory, and then calving areas, the variable Month was set as a random variable, allowing parameters for transition probability and autocorrelation to come from different probability distributions each month. This was different than a traditional BSAM model and allowed for potential differences in spatial characteristics of behaviors—ARS in foraging and calving grounds may present differently than ARS on the migratory route if it encapsulates breeding or calving behaviors or if the parameters of feeding bouts differ during migration (e.g. smaller turning angles for migratory ARS than foraging ground ARS). The model was fitted in R using the software JAGS [67] and the R rjags package [64, 68]. Where a gap of > 1 day existed in the raw satellite transmission data, the individual track was split and run as separate segments to avoid interpolating over long periods. Each model was run with two Markov chain Monte Carlo chains, consisting of 270,000 iterations each, the first 250,000 discarded as burn-in. The remaining 20,000 iterations were thinned, retaining every 8th sample to reduce autocorrelation and computational burden. The goodness of fit and chain convergence were assessed using the Gelman–Rubin statistic, and parameters with Gelman–Rubin (R) of less than 1.1 were considered converged as outlined by Gelman and Hill [56]. Runs were conducted on the UCSC Hummingbird computational cluster with chains running in parallel.