First Description of Migratory Behavior of Humpback Whales From an Antarctic Feeding Ground to a Tropical Breeding Ground

20 Background: Despite exhibiting one of the longest migrations in the world, half of the 21 humpback whale migratory cycle has remained unexamined; until this point, no study has 22 provided a continuous description of humpback whale migratory behavior from a feeding ground 23 to a breeding ground. We present new information on the satellite derived offshore migratory 24 movements of 16 humpback whales from Antarctic feeding grounds to South American breeding 25 grounds. Satellite locations were used to demonstrate migratory corridors, while the impact of 26 departure date on migration speed was assessed using a linear regression, and a Bayesian 27 hierarchical state-space animal movement model was utilized to investigate the presence of 28 feeding behavior en route. 29 Results : 35,642 Argos locations from 16 tagged whales from 2012-2017 were collected. The 16 30 whales were tracked for an average of 38.5 days of migration (range 10-151 days). The length of 31 individually derived tracks ranged from 645 – 6,381 km. Humpbacks were widely dispersed 32 geographically during the initial and middle stages of their migration but convened in two bottleneck 33 regions near the southernmost point of Chile as well as Peru’s Illescas Peninsula. The state space 34 model found almost no instances of ARS, a proxy for feeding behavior, along the migratory route. 35 The linear regression assessing whether departure date affected migration speed found suggestive 36 but inconclusive support for a positive trend between the two variables. No clear stratification by 37 sex or reproductive status, either in migration speed, departure date, or route choice, was found. 38 Conclusions: Southern hemisphere humpback whale populations are recovering quickly from 39 intense commercial whaling and, around the Antarctic Peninsula, are doing so in the face of a rapidly changing environment. The current lack of scientific knowledge on marine mammal 41 migration is a major barrier to cetacean conservation. This multi-year study sets a baseline against which the effects of climate change on humpback whales can be studied across years and conditions and provides an excellent starting point for the investigation into humpback whale migration.


INTRODUCTION 51
Humpback whale (Megaptera novaeangliae) migrations, with recorded one-way 52 distances of up to 8461km, are part of an annual cycle consisting of journeys between tropical 53 calving grounds in winter and high latitude feeding grounds in summer (1,2). Baleen whale 54 migrations are considered a response to the need to feed in cold waters and reproduce in warm 55 waters (1,2). Currently, NOAA recognizes 14 distinct populations of humpback whales, based on 56 breeding ground location, with seven in the Southern Hemisphere (3). These seven distinct 57 population segments (DPS) are found distributed around lower latitude coastal regions in the 58 Atlantic, Indian, and Pacific Ocean and rely on highly productive seasonal habitats in the 59 Antarctic, with several populations utilizing the Western Antarctic Peninsula, one of the most 60 rapidly warming areas in the world, as their foraging ground (4)(5)(6). 61 Humpback whales appear to generally remain loyal to their natal grounds and return for 62 breeding and calving purposes year after year. In the foraging grounds, the whales disperse 63 somewhat more broadly than in the breeding grounds, but with only limited overlap and 64 intermingling between populations that breed in different geographic areas (7). The population breeding off the western coast of South America is the Southeastern Pacific DPS. Historically, 66 animals were from the Southeastern Pacific DPS, which breeds off the Western coast of South 157 and Central America (3). Wildlife Computers (Redmond, WA, USA) SPOT5, SPOT 6,and 158 MARK 10 Platform Transmitting Terminals (PTTs) were utilized and tagging was limited to 159 adult-sized animals (>12m). Each tag was contained in a sterilized housing and was anchored in 160 the tissue beneath the blubber near the dorsal with stainless steel barbs, with the transmitting 161 antenna remaining free outside of the animal (5). Tags were deployed from a range of 3-10 m 162 from a Zodiac Mark V or a Solas ridged-hulled inflatable boat using an ARTS Whale Tagging 163 PLT compressed air system (32). 164 Satellite transmissions were activated via a salt-water switch, and locations of the whales 165 were obtained through the Argos System of polar-orbiting satellites (Argos, 1990). Tags were 166 programmed to transmit during specific hours and days. Since the tags were also being utilized 167 for other year specific projects, duty cycling varied across years. In 2012, tags were programmed 168 to transmit between 00:00-04:00 and 12:00-16:00 GMT. In 2013, tags were programmed to duty 169 cycle 3 hours on, 3 hours off, except for Sirtrack tags (identified by PTT IDs starting with 113), 170 which duty-cycled at 6 hours on/6 hours off. The 2015 tags were programmed to transmit 171 continuously, while in 2016 tags, some tags were programmed to transmit continuously, while 172 three were programmed to duty cycle at 1 day on, 4 days off. Tags deployed in 2017 were 173 programmed to duty cycle 12 hr on, 12 hr off. 174

Demographic Information 175
Skin and blubber biopsy samples were obtained from tagged whales whenever possible 176 using standardized remote biopsy techniques (33). Samples were obtained from the upper flank 177 below the dorsal fin (34). Blubber samples were used to provide life history and demographic 178 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 180 followed by a standard phenol-chloroform extraction method (35). To assign pregnancy within 181 sampled females, progesterone, a lipophilic steroid hormone, was quantified from a sub-sample 182 of blubber using a progesterone enzyme immunoassay (36). Pregnancy was then assigned by 183 comparing the measured progesterone concentrations across a pre-validated binary logistic 184 model developed from humpbacks of known pregnancy status sampled in the Gulf of Maine 185 (36). 186 187 Data Processing 188 R (version 3.4.3, R Core Team, 2017) was used to filter raw observations from the 189 satellite tags to remove points without location data, points with Argos error quality class Z 190 (invalid location), and points with duplicate timestamps. In addition, clearly implausible points 191 (e.g. on land or hundreds to thousands of kilometers from expected location) were visually 192 inspected and removed. Maps of the animals' tracks were plotted using ggmap (37) in R (R Core 193 Team, 2017). 194 Whales were determined to be migrating when they started a northward journey from the 195 WAP without any significant or lasting return movements. The date of departure for each whale 196 was determined visually by graphing latitude as a function of Julian day and assessing at which 197 point the animal moved northward without any return movements. The static nature of the 198 environmental data combined with the mobile nature of the humpback data's mobile nature 199 precluded us from statistically evaluating the potential environmental trigger of light. Instead, we 200 matched daylight hours in the WAP to tagging data and graphed this against the animals' 201 latitudes in the same fashion that we determined departure dates. 202 To determine rates of migration, speeds on the migratory route were calculated with data 203 corrected for location error with the simple default Hierarchical State Space Movement Model 204 with a 12-hour timestep fitted in R using BSAM (Jonsen 2016, R Core Team 2017. Rate was the 205 distance of the linear vector between 12-hour timestep locations. Distances between locations 206 were calculated using the function distanceTrack from the Argosfilter package (Freitas 2012, R 207 Core Team 2017). Average rates were calculated as the average of all 12-hour timestep rates for 208 each animal. 209 As coastal nations have exclusive sovereign rights for the purpose of conserving and 210 managing marine species within the bounds of their jurisdiction (38), the amount of time the 211 migrators spent within EEZ boundaries was calculated by summing the number of regular 212 timestep observations from the BSAM model within each country's national waters. While the 213 satellite tags themselves did not collect data with great regularity, the BSAM model calculates 214 true unobserved locations along regular time intervals from available data, and these intervals 215 were utilized for EEZ analysis. 216 There were a number of locations where the tracks converged and allowed for a logical 217 division of the migration corridor into three spatial sections, "WAP-Cape Horn (Drake passage)," 218 "Cape Horn (Chile) -Peninsula de Paracas (Peru)," and " Peninsula de Paracas (Peru)-Zona 219 Reserverda Illescas (Peru)." Since not all tags transmitted for the entire migratory journey, these 220 3 discrete spatial sections allowed for a more valid estimation and comparison of speeds in some 221 sections along the journey. Average migratory speed was calculated for each section, as well as 222 for the breeding area. As humpback whales leave the Antarctic peninsula at different times, a 223 simple linear regression was performed using Julian day (predictor variable) and speed (response 224 variable) to investigate whether the timing of migration affected the speed at which the animal 225 migrates. Because very few tags transmitted to completion of migration, we chose to look at 226 speed in the first migratory section from the WAP to Cape Horn (latitude = -55.9833). All data 227 above -55.9833, as well as all animals that did not reach -55.9833, were filtered out, and the 228 average speed over the section was calculated for each remaining individual. To correct for 229 issues of heteroskedasticity, speed was transformed with a log function, and the residual plot was 230 assessed for any obvious signs of nonlinearity and heteroskedasticity. A QQ plot was used to 231 check for the normality of residuals, and the data were tested for influential data points. To 232 determine whether sex and reproductive status had an impact on speed, two Welch's ANOVA 233 tests were performed on the same speed data, using sex (male/female) as the predictor variable in 234 the first test, and sex/reproductive status as the predictor variable in the second test (male, 235 female-pregnant, female-not pregnant). For all tests, P-values <.01 indicated strong support, p-236 values between .01 and .1 offered suggestive, but inconclusive support, and p-values>.1 237 indicated no support (39,40). 238 Discrete behavioral modes were determined by manually constructed hierarchical 239 Bayesian state-space movement models. This was a departure from the simpler models use to 240 assess true locations, as it allowed for differences in movement norms associated with behavioral 241 states depending on whether the animals were in the foraging grounds, breeding grounds, or 242 migratory route. This model associated spatial patterns of animal movement with predicted 243 behavioral states while simultaneously accounting for and correcting the significant error 244 inherent in Argos Satellite location data. 245 We used a discrete-time dynamic correlated random walk model following Jonsen et al. 246 (2005) and Bestley et al. (2013), where each movement stemmed from either a 'traveling or 247 'area-restricted search' (ARS) state (41,42). When humpback whales encounter sufficient prey areas, they often engage in ARS by decreasing their travel speeds and increasing their turning 249 angle radius and frequency; consequently ARS behavior is defined as shorter step lengths with 250 larger and more variable turning angles. The terminology ARS is used instead of foraging, as 251 whales may also be engaging in other behaviors such as resting and breeding in this state and our 252 measurements are not based off of a direct measure of feeding but rather use movement metrics. 253 In humpback whales this spatial signature may persist for up to several days in one location (43). 254 The traveling state, which is thought to occur when the animals are either actively migrating or 255 located in habitats unsuitable for foraging, is characterized by fast travel rates and infrequent and 256 small turning angles; in a state-space model this behavior is recognized by the presence of long 257 step lengths with small and infrequent turning angle radius. 258 The first component of the state space model was the process model, which estimates 259 animal behavior with a first-difference correlated random walk (42). The process model took the 260 where dt is the difference between true unobserved locations and coordinate vectors xt and xt-1 263 and N2 is a bivariate normal distribution with covariance matrix Σ, where 2 is the process 264 variance in longitude, 2 is the process variance in latitude, and is the correlation coefficient. 265 γ is the autocorrelation of direction and speed between consecutive locations, with a value of 266 between 0 and 1 (γ=0 would signal a simple random walk). bt is an index used to denote 267 behavioral mode, e.g. ARS or traveling. T(θ) is the transition matrix with mean turning angle θ 268 which provides the rotation required to move between dt and dt-1. in Weinstein et al. (2017aWeinstein et al. ( , 2017b (4,5). We used a timestep of 12 hours, which we deemed to be 292 a conservative balance between taking into account gaps in the data as well as ensuring 293 behaviors did not change between locations. Although only two behavioral states were modelled, 294 the means of the MCMC samples provided continuous values from 1-2. A mean behavioral 295 mode of <1.25 was considered traveling, whereas a value > 1.75 represented Area-Restricted 296 Search. Estimations between 1.25 and 1.75 were treated as uncertain (44). 297 To help address the inconsistent transmitting nature and duty cycling of the tags as much 298 as possible, a joint estimation, in which estimation of behavioral states is conducted jointly 299 across multiple animal movement datasets rather than individuals, was done. This method 300 assumes that movement parameters may differ among individuals but are drawn from the same 301 set of distributions, and allows the model to estimate parameters and state variables with greater 302 precision by assuming a general range in value for all animals to borrow strength across multiple 303 datasets, thus filling in for any animals with suboptimal data (45). 304 Priors for  and  were set to reflect the assumptions that the travelling state would have 305 greater autocorrelation and lower mean turning angles than the ARS state. To allow for variance 306 in transition probability and behavioral state characteristics as the animals switched from 307 feeding, to migratory, and then breeding areas, the variable Month was set as a random variable, 308 allowing parameters for transition probability and autocorrelation to come from different 309 probability distributions each month. This is different than a traditional BSAM model and 310 important as it allowed for potential differences in spatial characteristics of behaviors -ARS in 311 foraging and breeding grounds may present differently than ARS on the migratory route. This 312 model was fitted in R using the software JAGS (Plummer, 2013) and the R rjags package 313 (Plummer, 2016; R Core Team 2017). Where a gap of >1 day existed in the raw satellite 314 transmission data the individual track was split and run as separate segments to avoid 315 interpolating over long periods. Each model was run with two MCMC chains, consisting of 316 270,000 iterations each, the first 250,000 discarded as burn-in. The remaining 20,000 iterations 317 were thinned, retaining every 8 th sample to reduce autocorrelation and computational burden. 318 The goodness of fit and chain convergence were assessed using the Gelman-Rubin statistic, and 319 parameters with Gelman-Rubin (R ) of less than 1.1 were considered converged as outlined by 320 Gelman and Hill (2006)

Demographic information
Of the 16 animals that initiated migration, four were pregnant females, four were resting 336 females (one juvenile), four were males, and four did not have biopsy samples and were thus of 337 unknown sex. None of the animals were accompanied by calves at the time of tagging. 338 339

Individual data analyses 340
The start of migration, end of migration tracked, duration of migration tracked, number of 341 transmissions during migration, and length of migration tracked were found for each animal 342 (Table 1). The animals showed differences in regards to their migratory speeds, the start of 343 migration, and geographic routes. A summary of each of the 16 animals' individual movements 344 is provided in Table 1, and their routes can be seen in Figures 1 & 2.  (Table 2). 371 The average speed for all the animals was 5.88 km hr -1 (SD=1.31). In general, average 372 speeds followed a slow-fast-slow trajectory by track segment, with the average speed calculated 373 for the animals highest during the middle section of migration from Cape Horn to Peninsula de 374 Paracas, and lowest in the breeding area (Table 1, Figure 3). 15 migrators had tracks reaching to 375 Cape Horn, and their average speeds over the distance can be seen in Table 1. The regression 376 results showed suggestive but inconclusive support for the hypothesis that whales have faster 377 migratory speeds the later they leave the peninsula (F(1,13)= 4.117, p=.06346). There was no 378 relationship between speed and sex (F(2, 3.11=0.003, p=.96)) or speed and sex / reproductive 379 status (F(2, 4.8=.37, p=.71)). 380 The animals appeared to be almost exclusively traveling during their northward were documented in these areas, nor did animals deviate from their northbound migration to 416 enter the Straits of Magellan (47,48). It is worth noting that one individual recorded "Unknown" 417 behavior near Northern Chilean Patagonia. 418 Our migratory tracks tentatively identify the area around Zona Reservada Illescas, Peru, 419 as the start of the breeding area based on abrupt route change and the transition from transiting to 420 ARS in animal PTT=123224. This delineation of the breeding ground is more in agreement with 421 Guzman (21) than Rasmussen (2), which placed the border close to the equator in Salinas, Ecuador, more than 550 km away. Tagged whales in our study reached as far north as Panama,423 which was in agreement with Rasmussen's findings regarding the geographical extent of the 424 breeding grounds. 425 One tagged whale, PTT 123232, provided information on the complete migratory cycle 426 from the Antarctic to the tropical breeding ground and back to the Antarctic. While the tag 427 stopped transmitting for a significant portion of the northward migration, this deployment 428 represents the first tagged humpback to provide data for a continuous annual cycle. The 429 southward route lined up closely with the northward route, indicating that humpbacks may use 430 the same routes, regardless of migratory direction. 431 Interestingly, none of our migratory parameters lined up with one of the most touted 432 characteristics of migrationsegregation along sex, reproductive, and age classes. While our 433 sample size was not large (n=16), it was much larger than most similar cetacean telemetry 434 studies of migration, and the lack of stratification is notable. It may be possible that segregation 435 by sex and reproductive status has been overemphasized in past literature, that this pattern varies 436 by DPS and is not adhered to in the Southeastern Pacific DPS, or that there are additional 437 parameters that have not been accounted for. Our sample's nature can also explain some of the 438 discrepancies -Felix and Guzman (2014) looked only at southbound migration and hypothesized 439 that the coastal route vs oceanic route differed by whether the animal was a single adult or 440 mother with calf. By the time of the northbound migration, calves had already been weaned and 441 none of our tagged females were accompanied by offspring; therefore, the lack of observed 442 coastal route does not contradict their findings. 443 It is also of note that our findings seemingly oppose those of Avila et al (2020), which states that whale arrival in the breeding grounds is becoming consistently earlier, with an average 445 arrival date of the last week of May (10). Of our 16 animals, 8 had not even commenced 446 migration by the last week of May, let alone made it to the breeding grounds. 447 Our study supported Dawbin's (1966) conclusions on migratory foraging, which stated 448 that the animals did not forage on their northward migration. While telemetry data cannot 449 conclusively rule out foraging behavior, only 1% of our recorded locations on the migratory 450 route indicated ARS and all of these points belonged to one animal and occurred in the Drake 451 Passage. Without more detailed data (e.g. dive parameters) it would not be possible to determine 452 if this ARS included actual feeding behavior versus the myriad other reasons that an animal may 453 cease transiting for a short period of time. A few cases of behavior were classified as unknown 454 on the route, but the majority of points in this category were found in the breeding or foraging 455 areas. As previously stated, certain instances of feeding bouts have been recognized on the 456 migratory route in recent years (13,21,23,25,27,49). However, all of these recorded instances 457 have occurred while individuals were migrating from breeding to feeding grounds. It is possible 458 that supplementary feeding is a phenomenon relegated only to the route from breeding to feeding 459 groundsperhaps because there is less of a definitive date that whales need to reach their 460 destination by, or because energy stores are running low while on the journey from foraging to 461 breeding grounds whales have just replenished their food stores. 462 The average migratory rate for our animals was 5.88 km·h −1 ± 1.31 and 5.88 km·h −1 ± .59 463 for the complete tracks of the 5 animals that completed migration. Our animals completed the 464 migration in 41-54 days and traveled between 33°-43° per month. These speeds were 465 significantly faster than Dawbin (1966), who recorded south to north speeds of 15°per month, 466 with approximate migration durations of 60-120 days. They were slightly higher than previously 467 recorded telemetry speeds of 4.04± 1.08 km·h −1 (21), 4.3 ± 1.2 km· h −1 (13), 4.5 km·h -1 (14), and 468 3.83 and 3.48 km·h −1 (15,16). It is possible that the whales in our study utilized coastal currents, 469 such as the Humboldt Current, along the west coast of South America, to increase their traveling 470 speeds without incurring additional energetic costs. It is also possible that the Southeastern 471 Pacific DPS experiences slightly higher migratory speeds than other populations or that, 472 alternatively, migratory rates in the direction of the breeding ground are higher than that of the 473 return route given that the whales are at their maximum energy storage and are motivated to 474 establish themselves on breeding grounds. 475 The telemetry data also revealed that our humpback whale speeds, on average, were not 476 constant and tended to be highest in the middle of migration. If this is a typical pattern, it could 477 mean that many of the telemetry estimations in different studies of average migratory rates could 478 be biased if calculations are based on only a short portion of the route. 479 We found no evidence that migration was triggered by daylight hours. There was no 480 number of daylight hours at which all whales initiated migration. Instead, the whales departed 481 from the Antarctic in conditions ranging from two to eight hours of sunlight. Suggestive support 482 was offered for a positive relationship between migratory speed and departure date. This increase 483 of speed with a later departure date could indicate that animals feel compelled to make up for 484 lost time, presumably to arrive at the breeding ground in a coordinated manner. should be made to pick a duty cycling regime implemented consistently across years and 497 specifically with state-space model timesteps in mind. In addition, the JAGS model should 498 include the ability to fill in smaller gaps, as seen in BSAM and Jonsen (2007) (44). 499 500

Management Implications 501
The conservation of migratory species requires a knowledge of migratory routes' 502 geographical locations, which can highlight areas of particular importance to a species (29,31). 503 The humpback whales in this study spent the vast majority of their migratory time in territorial or 504 exclusive economic zone waters of several nations, and knowledge of the jurisdictions in which 505 the animals migrate can be taken into account when determining management policies as coastal 506 nations have exclusive sovereign rights for conserving and managing marine species within the 507 bounds of their jurisdiction (38). 508 To maximize conservation resources, the concept of site conservation, specifically 509 focusing resources on sites particularly important to a species' life history, has been developed 510 (50). Bottleneck sites, as well as breeding areas, are considered key areas (50). This study 511 identifies two bottleneck regions off Chile's coast and from Peru's Peninsula de Paracas up into Panama ( Figure 1B). These two areas represent regions to concentrate conservation resources 513 and pass legislation, and this information can be shared with the appropriate national 514 organizations to advance efficient and effective conservation measures such as Marine Mammal 515 Protected Areas (MMPA) (51). In addition, our data has been contributed to the Migratory 516 Connectivity in the Ocean project (MiCO), which is currently developing a system to aggregates 517 and generated actionable knowledge to support worldwide conservation efforts for numerous 518 migratory species (52). 519 520 521

CONCLUSION 522
Understanding humpback whale migratory behavior and routes gives us a greater context 523 to make effective and efficient conservation decisions in the face of the animals' changing 524 environment. This study is a starting point for the long-term monitoring of the animals in an era 525 of climate change. In the coming years, a significant challenge in the conservation of migratory 526 species will be migrants' potential to shift routes in response to their changing environment. 527 Long-term monitoring programs will allow conservationists and management specialists to 528 monitor and anticipate these changing behaviors (29), identify conservation priorities, and 529 provide baseline data against which the impacts of climate change on ecosystems and migratory 530 species can be highlighted (19,29). Future studies should continue to grow the sample size and 531 investigate routes, behaviors, sex, and reproductive segregation of migration. In particular, 532 emphasis should be given to the bottleneck region between Magellan and Northern Patagonia's 533 strait, to research whether or not our animals are feeding in this location on Antarctica's return 534 route. The information presented here currently defines the behavior of humpback whale 535 migratory behavior from feeding to breeding grounds and can serve as a baseline for future work 536 on the species to compare and contrast how different environmental conditions and populations 537 impact this behavior.

Availability of data and materials 574
The humpback whale datasets generated and or analysed during the current study are available in the 575 WhalePhys repository, https://github.com/bw4sz/WhalePhys/tree/master/Data/Humpback 576 The Palmer Station datasets analyzed during the current study are available in the Palmer Station 577 Weather -Daily Averages repository,     Average speeds of humpback whales by segment of migratory route Figure 4 ARS, traveling, and unknown behavior exhibited by satellite tagged humpback whales on their northward migration from Antarctica