Multi-Decadal Humpback Whale Migratory Route Fidelity Despite Oceanographic and Geomagnetic Change

Understanding how organisms respond to environmental change is one of the most pressing grand challenges of organismal biology. In the vast oceans that cover 71% of Earth’s surface, remote sensing technologies have created unprecedented opportunities to create new knowledge and deliver integrated understandings of marine organism-environment interactions via long-term monitoring. Using historic whaling records and >15 years of satellite-derived data, we show that movement parameters associated with long-distance humpback whale migrations, including utilization of a south-southeast directed migratory corridor, migration path straightness, direction, timing, and velocity, have not significantly changed during a period of dynamic oceanographic and geomagnetic conditions. These findings reveal an apparent paradox: humpback whale migrations do not change in a changing ocean. Geophysical analyses of the same humpback whale movements demonstrate that these whales maintained prolonged migratory fidelity to a limited suite of spatiotemporal trajectories through gravitational coordinates, raising the possibility that migratory decisions are relatively insensitive to changing oceanographic and geomagnetic conditions. Our findings highlight the importance of filling the knowledge gaps that currently limit our ability to understand and anticipate organismal responses to rapidly changing Earth system conditions.


INTRODUCTION
Migratory animals face an uncertain future (Cotton, 2003;Wilcove and Wikelski, 2008;Hazen et al., 2013;Hof et al., 2017;Cohen et al., 2018). Global change, climate change, development, resource extraction, habitat fragmentation, and other human-induced perturbations combine to present unprecedented challenges to the sustainability of migratory populations and the ecosystem services they provide (Peñuelas et al., 2002;Harris et al., 2009Harris et al., , 2018Roman and McCarthy, 2010;Thackeray et al., 2016;van Doren et al., 2017). Ensuring that sustainable populations of migratory animals persist into the future requires knowledge of how migrants respond to Earth system dynamics (Schwenk et al., 2009;Bowlin et al., 2010).
It is widely agreed that addressing the challenges inherent to a dynamic Earth system requires long-term monitoring delivered through integrated interdisciplinary research approaches (Schwenk et al., 2009;Bowlin et al., 2010;Hays et al., 2016;Urban et al., 2016;Miloslavich et al., 2018). At present, models intending to predict the effects of Earth system dynamics on the biosphere suffer from both a lack of data and a limited understanding of the dominant biological mechanisms through which adaptations occur, including: phenology and demography; physiology; range dynamics; evolutionary potential; species interactions; and organismal responses to environmental variability (Urban et al., 2016). Robust prediction of future scenarios requires monitoring of the relevant variables at appropriately matched spatial and temporal scales (Hazen et al., 2013;Pereira et al., 2013;Urban et al., 2016;Miloslavich et al., 2018). Although significant progress in remote monitoring has occurred, relatively few studies integrate long-term datasets of both organismal behavior and the dynamic environments they inhabit (Sprogis et al., 2018;Abrahms et al., 2019). The primary goal of this research was to explore gaps in our knowledge regarding how humpback whale (Megaptera novaeangliae) long-distance migrations respond to oceanographic and geomagnetic change through long-term remote ocean monitoring.
Environmental changes should evoke biological responses, including habitat shifts in marine organisms (e.g., Hoegh-Guldberg and Bruno, 2010;Hazen et al., 2013), and large whales are no exception (Santora et al., 2020). For example, recent research suggests that ocean circulation dynamics in the North Atlantic have forced changes in right whale (Eubalaena glacialis) foraging patterns (Record et al., 2019), and in the North Pacific, humpback whale feeding behavior has been shown to respond to oceanographic conditions, including sea surface temperature, upwelling, and biomass (Fleming et al., 2016). Along the coast of California, habitat compression associated with the 2014-2016 marine heatwave likely contributed to an order of magnitude increase in annual humpback whale entanglement rate (Santora et al., 2020). In the South Pacific, Oceania's changing climate is predicted to cause distribution shifts in endangered humpback whales (Derville et al., 2019). In the Southwest Atlantic, environmental variables, including ocean currents and sea surface temperatures, were found to be strong predictors of humpback whale distribution in the Abrolhos Bank breeding area (Bortolotto et al., 2017). Ocean water temperature has also been linked to changes in humpback whale foraging behavior in the Southern Ocean (Owen et al., 2019). Although the literature clearly establishes that large whale distribution and behavior are sensitive to environmental conditions in their breeding and feeding areas, the cues and conditions that inform long-distance migratory behavior remain largely unknown.
Current evidence suggests that whale movement decisions are likely informed by some combination of magnetic, oceanographic, and gravitational cues. For example, integrating magnetic cues into whale behavior research has led to novel insights into whale stranding and navigation (Kirschvink et al., 1986;Horton et al., 2017). Oceanographic conditions, including temperature, productivity, and currents, are known to be key variables in whale distribution (e.g., Horton et al., 2011;Hazen et al., 2013;Fleming et al., 2016;Derville et al., 2019;Owen et al., 2019), and gravity is an inescapable driver of whale buoyancy and movement behavior (e.g., Clarke, 1978;Horton et al., 2017).
Yet, oceans are highly dynamic, and in the Southwest Atlantic magnetic and oceanographic cues are anything but constant. Sophisticated long-term analyses of satellite-derived oceanographic observations of the southwest Atlantic's subtropical gyre, where the Malvinas (10-88 Sv) and  currents converge, demonstrate that oceanographic conditions throughout the region, including water transport, eddy kinetic energy, sea height anomaly, sea surface temperature, and ocean current circulation patterns, are extremely variable throughout the region (Goni et al., 2011;Wu et al., 2012;Marcello et al., 2018). With respect to magnetism, the South Atlantic Anomaly (SAA) is the most rapidly changing feature in Earth's magnetic field over the past 400 years (Hartmann and Pacca, 2009). In marked contrast, latitude and bedrock dependent gravitational cues are relatively stable despite the presence of pronounced gravitational anomalies, like the Rio Grande Rise (Mohriak et al., 2010).
Because the Southwest Atlantic is known to be dynamic with respect to oceanographic and magnetic conditions, but stable with respect to gravity, it is an ideal region in which to explore possible correlations between movement behaviors and oceanographic and geophysical cues, including gravity. By integrating long-term humpback whale satellite tracking and historic whaling datasets with satellite-derived oceanographic, magnetic, and gravitational monitoring and mapping datasets, our research illuminates one of biology's most pressing questions: How do animal migrants respond, if at all, to environmental dynamics?

MATERIALS AND METHODS
We present satellite-monitored essential biodiversity and ocean variables during a multi-decadal period of humpback whale migratory route fidelity in the Southwest Atlantic Ocean. Our study integrates published whaling records with several satellite-derived remote sensing datasets, including: (1) the latitude, longitude, and date on which 243 humpback whales were killed by the Soviet Yuri Dologorukiy fleet in the Southwest Atlantic Ocean between 1965 and 1973 (Zemsky et al., 1996; Figure 1A); (2) platform transmitting terminal (PTT) satellite telemetry data that document the spatial and temporal locations of 20 humpback whales as they migrated through a <350 km wide and >3000 km long south-southeast corridor (SSEC) in the South Atlantic Ocean between 2003 and 2018 ( Figures 1B,C); (3) magnetic field parameters determined using the Swarm satellite derived Enhanced Magnetic Model (Chulliat et al., 2015); (4) Terra satellite derived Moderate Resolution Imaging Spectroradiometer (MODIS) estimates of near-surface temperature (NASA, 2018) and chlorophyll-a concentrations (NASA, 2018); (5) Ocean Surface Current Analysis Real-time (OSCAR) estimates of near-surface ocean  1965-1966; 189 whales killed between the Abrolhos Bank calving ground and Rio Grande Rise during austral Spring 1967; and 2 whales killed in South Georgia Basin on February 15-16, 1973(Zemsky et al., 1996, (B) 2003-2009 tracks of 9 whales (individual platform transmitting terminal numbers as indicated in the legend), (C) 2010-2017 tracks of 11 whales (individual platform transmitting terminal numbers as indicated in the legend). White lines depict the 200 m bathymetric contour. Gray-scale basemap depicts 1 arc-minute resolution (Amante and Eakins, 2009) bathymetry. White shaded polygon corresponds with the 350 km wide south-southeast corridor (SSEC) used by the 20 tracked whales. Highly productive surface waters (i.e., mesotrophic waters with >2.6 g m −3 chlorophyll-a) for the month of December, when the tracked whales first arrived in the South Georgia Basin feeding grounds, are shown as colored polygons as indicated in the legends.
currents based on data collected by Jason-1/Jason-2 and Gravity Recovery and Climate Experiment (GRACE) satellites (ESR, 2009); (6) latitude and bedrock dependent gravitational acceleration data derived, in part, from TOPEX/Poseidon satellite observations (Balmino et al., 2012;Götze, 2014). We integrated these datasets through geospatial and time-series analyses that collectively reveal how humpback whales respond, and fail to respond, to changing Earth system conditions. The specific hypothesis we tested was: Southwest Atlantic humpback whale migratory movements describe systematic and highly reproducible sinusoidal trajectories when plotted in spatiotemporal gravitational coordinates .

Humpback Whale Locations and Movement Variables
As part of a long-term monitoring program, 138 PTT tags were deployed on humpback whales seasonally residing off the coast of Brazil between 2003 and 2018 (Zerbini et al., 2006this study). Humpback whales were tracked using published methods (Zerbini et al., 2006Horton et al., 2011Horton et al., , 2017. In brief, SPOT radio-frequency PTT satellite tags (Wildlife Computers, Redmond, WA, United States) were transdermally implanted into the upper flank of each whale near the base of the dorsal fin using a carbon-fiber pole or a modified pneumatic line-thrower. In all references to PTT tag numbers in the current study, the two digits to the right of the decimal point correspond with the abbreviated Julian calendar year in which the tag was deployed.
Raw humpback whale location estimates were assigned a location class (i.e., 3, 2, 1, 0, A, B, Z) by the Argos-CLS system based on the estimated location error and the number of messages received (Argos-CLS, 2016). All locations that passed a 20 km h −1 velocity filter were used in this study. Velocityfiltered locations were combined to determine single average daily locations for each whale using PAST (v. 3.26;Hammer et al., 2001) to ensure that individual tracks were uniformly distributed with respect to time. The humpback whale tracking research we report was performed in accordance with research approvals granted by the Brazilian Environmental Agency (IBAMA), permit #009/02/CMA/IBAMA and process #02001. 000085/02-27.
The date and location of 243 humpback whales killed by the Soviet Yuri Dologorukiy fleet between 1965 and 1973 were derived from the records reported by Zemsky et al. (1996;see also: Zerbini et al., 2006). Our understanding is that the Yuri Dologorukiy fleet humpback whale kills we report were hidden, in hard copy form, in the potato cellar of Dimitry Tomorsov until the end of the Cold War.
Humpback whale movement variables, including movement direction, distance traveled, straightness index, and movement velocity, were determined using the equations (1-4) as reported below. Movement direction was determined by: where: α is the whale's movement direction in degrees; λ 1 and ϕ 1 are the whale's longitude and latitude, respectively, in decimal degrees at location 1; λ 2 and ϕ 2 are the whale's longitude and latitude, respectively, in decimal degrees at location 2; e is the eccentricity of the spheroid (i.e., 0.081819791). Distance traveled was determined by: where: S is the distance traveled between locations; a is the length of the major semiaxis in km (i.e., 6378.137 km); a is the whale's movement direction in degrees; e is the eccentricity of the spheroid (i.e., 0.081819791); ϕ 1 is the whale's latitude in decimal degrees at location 1; ϕ 2 is the whale's latitude in decimal degrees at location 2. Since the Argos-CLS system provides point locations rather than movement paths, the distance traveled data are minimum estimates of the true distance traveled between sequential locations. Straightness index was determined by: where: SI is the straightness index as defined by Batschelet (1981); D is the finite rhumb line distance between a starting location and an end location; L is the length of path followed between the starting location and the end location. Movement velocity was determined by: where: v is the velocity; S is the distance traveled (see: equation 2); t 2 -t 1 is the time it took for the whale to travel from location 1 to location 2. Since S is necessarily a minimum estimate of the true distance traveled, the movement velocity values we report will also be minimum estimates of true movement velocity. Kernel density estimation (Silverman, 1986) was used to determine the number of modes present in the humpback whale movement velocity data distribution.

Spatiotemporal, Geophysical, and Astronomical Variables
Spatiotemporal, geophysical, and astronomical variables were determined for the humpback whale locations we report. The time at which each whale initiated its continuous southward migration away from Abrolhos Bank was determined via piecewise linear regression breakpoint analysis (PLR-BPA) using the 'segmented' package in [R] (Muggeo and Muggeo, 2017).
We used PLR-BPA, rather than state-space switching models, as PLR-BPA enables quantification of the time at which movement behavior likely changed between sequential PTT locations. These piecewise linear regression breakpoint analyses allowed us to determine both when and where major changes in humpback whale movement behavior occurred, including: (1) migratory departures, represented in the tracking data as the time and place where each whale initiated continuous directional movement; (2) migratory stop-overs; (3) migratory re-starts following stop-overs. Geomagnetic field parameters, including inclination (I) and intensity (F) at each PTT location, were determined using the Enhanced Magnetic Model (Chulliat et al., 2015). Latitude and bedrock dependent gravitational accelerations at each PTT location were determined using the International Gravity Formula (Götze, 2014) and spherical harmonic analysis of the Earth's topography-bathymetry ETOPO1 (Amante and Eakins, 2009) dataset up to degree and order 10,800 as reported in the International Gravimetric Bureau's 2 × 2 arc-min (i.e., ∼3.7 × ∼3.7 km) World Gravity Map (Balmino et al., 2012). Astronomical variables, including moon illumination, Earthmoon distance, and lunar declination, were calculated using astronomical algorithms (Meeus, 1991). The magnitude of the tidal gravity vector at the Abrolhos Bank migratory departure site of PTTs 24641.05, 121203.17, and 120942.17 were determined using ETIDE (Fisahn et al., 2012(Fisahn et al., , 2015.
Using these geophysical data we determined the latitude and bedrock dependent gravitational accelerations experienced by each whale. We then normalized the sum of these two spatially dependent gravitational cues to the gravitational acceleration present at the individual humpback whale migratory departure sites identified using PLR-BPA. This departure site normalization was done by: where: g L is the latitude dependent gravitational acceleration at the whale's location; g B is the bedrock dependent acceleration (i.e., Bouguer gravity anomaly); g Ld is the latitude dependent gravitational acceleration at the whale's migratory departure site identified by PLR-BPA; g Bd is the bedrock dependent gravitational acceleration at the same migratory departure site for which g Ld was determined. We plotted these departure sitenormalized gravity data against moon illumination, a temporally dependent visual and gravitational cue. We analyzed the resulting migratory trajectories using PAST's (v. 3.26;Hammer et al., 2001) sinusoidal regression. In this study, we restricted our least-squares regression analyses to a single sinusoid of the general form, where: a is the amplitude of the sinusoid, which we did not restrict to the natural range in moon illumination values (i.e., 0 to 1); g L , g Ld , g B , and g Bd are the same as in equation (5); Frontiers in Marine Science | www.frontiersin.org is the sinusoidal period; is the sinusoidal phase. PAST minimizes the residual sum of squares via an ordinary least squares regression process by adjusting the fitted sinusoid's amplitude, period and phase.

Oceanographic Variables
Oceanographic variables, including sea surface temperatures, ocean currents, and Chlorophyll-a concentrations, were extracted from satellite monitored raster datasets at humpback whale locations using ArcGIS Desktop (ESRI, 2011). Sea surface temperature and chlorophyll-a concentration rasters (4 km pixel size) were obtained from the Goddard Space Flight Center (NASA, 2018); zonal and meridional ocean surface currents (1/3 • pixel size) were estimated using the ocean surface current analysis real-time (OSCAR) model provided by JPL Physical Oceanography DAAC and developed by Earth and Space Research (ESR, 2009). In order to determine if any significant trends were present in the oceanographic conditions the tracked whales experienced, we performed Mann-Kendall trend tests (Gilbert, 1987), a non-parametric test for the presence of significant temporal trends in time-series data, on equally spaced sea-surface temperatures in the SSEC, total areal extent of mesotrophic ocean surface waters (i.e., >2.6 g m −3 chlorophyll-a concentration) in the southwest Atlantic/Southern Ocean basin, and ocean surface current direction within the SSEC every December during the 2002-2017 satellite tracking period. Mann-Kendall trend tests were also performed on the 2003-2018 humpback whale movement parameters: straightness, swimming direction, date of migration onset, and swimming velocity.

RESULTS
The data we report include six primary findings. First, the historic Soviet whaling and satellite telemetry data demonstrate that Southwest Atlantic humpback whales have utilized a spatially restricted ∼1.0 million km 2 migratory corridor (i.e., SSEC) for >50 years (Figure 1). Second, fidelity to the SSEC manifests as relatively constant movement parameters, including straightness (Batschelet, 1981), movement direction, timing, and velocity when analyzed with respect to time. Third, humpback whale movement velocities in the SSEC was bimodal. Fourth, orientation cues derived from Earth's magnetic field changed significantly (p << 0.05) during the 2003-2018 monitoring period. Fifth, essential oceanographic variables (Miloslavich et al., 2018), including sea surface temperature, ocean current direction, and productivity, were highly variable during the 2003-2018 monitoring period. Sixth, the humpback whale migrations we report describe highly significant and reproducible sinusoidal gravitational coordinate (i.e., g-space) trajectories.

Soviet Whaling and Humpback Whale Satellite Telemetry
Between 1965 and 1973 the Soviet Yuri Dologorukiy fleet killed 243 humpback whales in the Southwest Atlantic Ocean (Zemsky et al., 1996). Fifty-two of these whales were killed on their feeding grounds in the South Georgia Basin during the austral summer of 1965-1966. Not long after, during a 2-week period between 31 October and 13 November, 1967, the same fleet killed 189 humpback whales on the Abrolhos Bank calving grounds or within the SSEC (Figure 1A). Prolonged utilization of this same south-southeast directed migratory corridor (i.e., SSEC) is confirmed by the 2003-2018 humpback whale tracking dataset reported below.
The humpback whale tracking dataset we analyzed includes 45 PTT tags (i.e., 33% of all tags) that successfully recorded the onset of individual southward migrations away from the continental shelf and migratory movement south of −24 • S latitude. Of these 45 migratory movements, 20 whales followed open-ocean paths that fell within the 350 km wide and 3000 km long SSEC located between −20 • S and −45 • S latitude (Figures 1B,C). Of the 20 whales that utilized the SSEC, only 6 were successfully tracked south of −45 • S latitude and into the population's feeding grounds within and adjacent to the South Georgia Basin (Figures 1B,C). Satellite telemetry tracking data for 5 of these 20 whales (PTTs: 27259.03; 10946.05; 26712.05; 24641.05; 87771.09) have been published previously (Zerbini et al., 2006Horton et al., 2011).
The PTT satellite tracking data demonstrate that 44% of migrating whales (i.e., 20 of 45 tracked whales) utilized less than 10% of the available ocean area during a >15 year-long period, consistent with the finding that multiple species of marine megafauna are capable of prolonged spatial and temporal fidelity to well-defined migratory domains . Movement parameters, including straightness (Figure 2A), swimming direction (Figure 2B), the timing of migratory onset (Figure 2C), and swimming velocity (Figure 2D), for the 20 whales migrating through the SSEC, showed no significant trends through time (p > 0.05; non-parametric Mann-Kendall trend test).
However, kernel density estimation revealed that the swimming velocity data distribution was bimodal (Figure 3), with a slow mode peaking at 3.23 km h −1 and a fast mode peaking at 4.54 km h −1 . The swimming velocities we report are consistent with PTT derived humpback whale swimming velocities reported by others, including: (1) 3.4 and 3.6 km h −1 in the Southwest Pacific Ocean (Hauser et al., 2010;Riekkola et al., 2018); (2) 3.6 km h −1 in the Indian Ocean (Trudelle et al., 2016); (3) 4.1 km h −1 in the North Atlantic (Kennedy et al., 2014); (4) 4.5 and 6.25 km h −1 in the North Pacific (Mate et al., 1998). Due to the bimodal distribution we observed, we classified individual whales as either slow (<4.0 km h −1 ) or fast (≥4.0 km h −1 ), based on their average migratory swimming velocities ( Table 1).
The bimodal swimming velocity distribution is likely real, and not an artifact of the potential PTT underestimation of true swimming velocity, because of the straightness of the migratory paths followed by these whales (Table 1 and Figures 1, 4). To demonstrate this point, we calculated the cumulative migratory distance traveled, using both raw Argos PTT location estimates and our interpolated whale locations, for both a slow (PTT 87775.08) and a fast (PTT 121203.17) whale during the first ∼1000 km of their southward migrations ( Figure 4C). As expected, the more frequent Argos PTT location estimates yielded higher  (Gilbert, 1987) probabilities (i.e., p-values) that a significant (i.e., p < 0.05) temporal trend is present.
cumulative distances traveled at any point in time during the satellite monitored migratory movements. However, there is an extremely significant difference (p < 0.05; t-test on the linear regression slopes) between the fast and the slow swimming velocities for both the raw Argos PTT data and the interpolated daily locations. Although interpolation of higher temporal resolution PTT location estimates minimizes movement distances and velocities, these effects do not obscure significant differences in movement parameters in directional long-distance migration data. Noad and Cato (2007) demonstrate that singing humpbacks swim significantly slower than non-singing humpbacks during long-distance migration, suggesting a potential sex-related or behavioral driver for the bimodal swimming velocity distribution we report, and further investigation of the role acoustics play in long-distance movement behaviors is warranted. Piecewise linear regression breakpoint analysis identified significant breakpoints associated with migratory departures for 15 of the whales that migrated through the SSEC (Figure 5A). The satellite tracking data for the remaining 5 whales that used the SSEC did not include significant breakpoints associated with the onset of southward migration because either: (1) the whale was already migrating southward when its PTT tag was deployed (e.g., Figure 5B); or (2) the PTT dataset included a >3 daylong gap spanning the onset of continuous southward movement (e.g., Figure 5C).
Humpback whale PTT 120942.17 was the only whale to perform a multi-day stop-over with associated decisions to stop and (re)start directed migratory movements ( Figure 5A). The complexity of 120942.17's track is noteworthy as it demonstrates spatiotemporal fidelity to a well-defined migratory path at the ecological expense of both time and energy. Following ∼15 days and >1700 km of continuous south-southeast (159 • ) directed swimming away from Abrolhos Bank, 120942.17 stopped over above Rio Grande Rise. Unexpectedly, 120942.17 resumed continuous directed swimming several days later along an overall north-northwest bearing (339 • ), antithetical to its initial southsoutheast directed movement ( Figure 2B and Table 1). Following ∼10 days and ∼1000 km of swimming along this "reverse" trajectory, 120942.17 again stopped its directed swimming for a single overnight period, only to return to its original southsoutheast directed migration 2 h (±1 h) before sunrise on 11 November, 2017.
In total, PLR-BPA identified 20 significant changes in latitudinal time-series for the 20 whales migrating through the SSEC. Of these significant changes in movement behavior, 15 occurred within ±4.5 h of sunrise and 5 occurred within ±5 h of sunset ( Figure 5D). A significant local peak in the data distribution (p = 0.0007, two-tailed exact binomial probability), centered on dawn twilight, includes nine of the twenty breakpoints identified. Our breakpoint analyses demonstrate that crepuscular navigation is a common component of humpback whale movement behavior.

Geomagnetic and Oceanographic Conditions
In contrast to the prolonged stability of the humpback whale movement parameters, the magnetic field these whales swam through changed significantly during the study period. For example, geomagnetic inclination (I) significantly changed (p < 0.05; non-parametric Mann-Kendall trend test) by >12.5% at the humpback whale wintering grounds on Abrolhos Bank between October 2003 and October 2017 ( Figure 6A).
The relationship between this significant change in geomagnetic coordinates and geographic coordinates is notable. Magnetic inclination angles, equivalent to those present at Abrolhos Bank in 2003, were located >340 km to the northwest and >130 km inland of the Brazil coast by 2017.
Secular variation in magnetic orientation cues were more severe at the southern end of the SSEC where equivalent magnetic inclinations (I) shifted >400 km to the west between 2003 and 2017. Equivalent magnetic field intensity-inclination polar bi-coordinate locations (Lohmann et al., 1999;Brothers and Lohmann, 2018) present in 2003 no longer existed in 2017. The severity of these temporal changes in the magnetic field are due to the South Atlantic Anomaly (SAA), the most rapidly changing feature in Earth's magnetic field over the past 400 years (Hartmann and Pacca, 2009).
In addition to a changing magnetic field, the whales we tracked also experienced dynamic oceanographic conditions. Although no significant trend was detected (p > 0.05; nonparametric Mann-Kendall trend test), between November 2003 and November 2017, average sea surface temperatures in the SSEC increased by approximately 0.3 • C (Figure 6B), roughly three times larger than the contemporaneous global trend (Hausfather et al., 2017).
It is not surprising that the humpback whales we tracked through the SSEC experienced dynamic ocean surface currents, annually varying across a >90 • range of predominantly headcurrents during late austral spring humpback whale migration period ( Figure 6C). Such considerations are important as animals moving through a flowing medium, such as air or ocean, are affected by the direction and magnitude of currents, and these currents have the potential to impact animal movement trajectories (Gaspar et al., 2006;Chapman et al., 2011). To consider this possibility, we analyzed both uncorrected and current-corrected (Gaspar et al., 2006) humpback whale tracks. Our analyses demonstrate that: 1) uncorrected and currentcorrected humpback whale locations were not significantly different (average drift = 9.2 km day −1 ; paired two-tailed t-test, p > 0.05; see: Supplementary Figures S1, S2), there is no significant difference between uncorrected and currentcorrected humpback whale movement trajectories (two-tailed paired t-test, p > 0.05; see: Supplementary Figure S2). Despite the dynamic nature of near-surface ocean currents in the SSEC, our analyses demonstrate that these currents have    Figure 1), (B) humpback whale PTT 87778.10's latitude time-series profile, wherein no significant breakpoint was detected due to continuous southward movement, (C) humpback whale 87775.12's latitude time-series profile, wherein no significant breakpoint was detected due to a multi-day gap in the satellite telemetry dataset, and (D) histogram showing the temporal distribution of the 20 significant breakpoints identified using PLR-BPA relative to sunrise (light shading) and sunset (dark shading).
relatively minor effects on migrating humpback whale trajectories due to the factor 10 difference between average current velocity (0.43 km h −1 ) and average whale movement velocity (4.6 km h −1 ).
With respect to feeding, the primary motivation of poleward humpback whale migrations, the areal extent of highly productive (i.e., >2.6 g m −3 chlorophyll-a) surface water decreased significantly (p < 0.05; non-parametric Mann-Kendall trend test)

Movements in Gravitational Coordinates
A central goal of this research was to test the hypothesis that humpback whale migratory movements describe highly significant sinusoidal gravitational coordinate (i.e., g-space) trajectories. The results of these analyses demonstrate that humpback whales utilizing the SSEC followed a limited range of temporally modulated gravitational coordinate trajectories (Figure 7). These highly significant sinusoidal correlations (sinusoidal regression; F-test; α = 0.05; p < 0.05) support the hypothesis as posed and suggest that the observed g-space trajectories are a non-random consequence of navigational FIGURE 7 | Humpback whale movements in gravitational coordinates (i.e., g-space). (A) g-space trajectories for slow swimming (i.e., <4 km h −1 ) humpback whales (n = 9, symbols as in Figure 1), (B) g-space trajectories for fast swimming (i.e., >4 km h −1 ) humpback whales (n = 11) and the 243 Southwest Atlantic humpback whales killed by the Soviet Yuri Dologorukyi fleet (Zemsky et al., 1996(Zemsky et al., ) between 1965(Zemsky et al., and 1973 (symbols as in Figure 1). Normalized gravity values in (A,B) were determined by dividing the sum of the latitude dependent gravity (Götze, 2014) and the bedrock dependent gravity (Balmino et al., 2012) at each whale location by the value present at each individual's Abrolhos Bank migratory departure site. Moon illumination depends solely on time, and it cyclically evolves from no illumination (i.e., New Moon, 0) to full illumination (i.e., Full Moon, 1) to no illumination (i.e., New Moon, 0) over a 29.5 day average period (i.e., synodic month).
behavior (Figure 7). Although the twenty whales that migrated through the SSEC followed distinct geographic-Julian calendar coordinate paths (e.g., Figure 5A), these same movements describe overlapping temporally phased trajectories when plotted in gravitational coordinates (Figure 7). Yet, two different motifs, consistent with the two modes in the humpback whale movement velocity distribution, are apparent in the g-space trajectories: a slow motif ( Figure 7A) and a fast motif ( Figure 7B).
The slow motif includes 9 whales that swam through the SSEC at an average swimming velocity of 3.6 ± 0.9 km h −1 , whereas the 11 whales in the fast motif swam significantly faster (p = 0.008; two-tailed t-test), at an average velocity of 4.8 ± 0.9 km h −1 ( Table 1). The data demonstrate that eight of nine whales in the slow motif initiated their southward migrations in the days immediately prior to first-quarter or last-quarter moons (Figure 7A), while the ninth whale in the slow motif, PTT 27259.03, initiated its migration 2 days after first-quarter moon. In contrast, the eleven whales in the fast motif initiated their southward migrations near full or new moons ( Figure 7B). Both the slow and fast motifs include significant in-phase and out-of-phase sinusoidal correlations between departure site-normalized gravitational acceleration and moon illumination (Figure 7), similar to significant correlations present in the g-space trajectories followed by diverse megafaunal species, including white sharks (Carcharodon carcharias) and northern elephant seals (Mirounga angustirostris) . Based on these results, we accept the hypothesis that humpback whale migratory movements describe highly significant sinusoidal g-space trajectories.
The results we report demonstrate that, despite changing oceanographic and geomagnetic conditions, humpback whale long-distance migrations in the Southwest Atlantic Ocean are characterized by prolonged (i.e., >50 years) fidelity to a heavily utilized south-southeast directed migratory corridor connecting calving grounds on the Abrolhos Bank to feeding grounds in the South Georgia Basin. When analyzed in spatiotemporal gravitational coordinates, these same humpback whale migratory movements describe richly patterned and highly reproducible trajectories in either a fast or a slow motif. Thus, the historic Soviet whaling and humpback whale satellite telemetry data we report reveal an apparent paradox: humpback whale longdistance migrations do not change in a changing ocean.

DISCUSSION
Understanding organismal responses to environmental change remains one of the most pressing challenges in ecology (Schwenk et al., 2009), and as relatively high trophic level predators, it is particularly important to determine how large whales respond to dynamic oceanographic conditions given the role cetaceans play in circulating marine nutrients, trophic architecture, trophic cascades, and carbon turn-over rates (e.g., Croll et al., 2006;Roman and McCarthy, 2010;Witteveen and Wynne, 2016). In this context, it is somewhat surprising that Southwest Atlantic humpback whale movements do not appear to have responded to a significant decrease in the areal extent of mesotrophic surface waters and dynamic oceanographic and geomagnetic conditions over a period of several years. This paradox of migratory fidelity despite environmental change raises important questions, including: (1) How do humpback whales maintain migratory stability despite oceanographic and geomagnetic change? (2) What cues and rules inform humpback whale migratory decisions? One possible answer to both questions invokes the spatially and temporally dependent, yet stable, gravitational cues ubiquitous to the Earth system .
The recurrent pattern of phased gravitational coordinate trajectories in our long-term satellite tracking dataset supports FIGURE 8 | Geographic and gravitational coordinate migratory trajectories of PTT 24641.05's slow motif humpback whale migration (average migratory velocity = 3.5 ± 1.0 km h −1 ). (A) Geographic coordinate (plate-carée projection) track map of 24641.05's single stage south-southeast directed migration, (B) latitude time-series plot for 24641.05's geographic coordinate Julian calendar migratory trajectory, (C) g-space trajectory of 24641.05's migration track, including the actual and mirror image g-space trajectories. The mirrored trajectory (dashed black curve) is plotted on the reverse direction moon illumination and normalized gravity axes. Actual migratory trajectories in (A-C) are shown as red-hued triangles that are color-coded and sized according to 24641.05's average daily velocity as shown in the legend. Normalized gravity values in (C) were determined by dividing the sum of the latitude dependent gravity (Götze, 2014) and the bedrock dependent gravity (Balmino et al., 2012) at each average daily location by the value present at 24641.05's Abrolhos Bank migratory departure site. Colored basemap in (A) depicts bedrock gravity anomalies across a 150 to 660 mGal range. the interpretation that humpback whale migrations are both stable and non-randomly distributed in space and time (e.g., Figure 7). In the following sections, we discuss three separate examples from our tracking dataset that describe the conditions under which different humpback whales were able to maintain migratory fidelity to a limited range of spatiotemporal trajectories through gravitational coordinates.
Several rare conditions were met during 24641.05's migration that collectively determined the gravitational accelerations it experienced. First, 24641.05 initiated its migration away from Abrolhos Bank on a quarter moon neap tide, when: (1) the tidal gravity vector achieved an anomalously low daily range in acceleration values (∼100 µGal; Supplementary Figure S3); FIGURE 9 | Geographic and gravitational coordinate migratory trajectories of PTT 121203.17's fast motif humpback whale migration (average migratory velocity = 5.0 ± 1.3 km h −1 ). (A) Geographic coordinate (plate-carée projection) track map of 121203.17's single stage south-southeast directed migration, (B) latitude time-series plot for 121203.17's geographic coordinate Julian calendar migratory trajectory, (C) g-space trajectory of 121203.17's migration track, including the actual and mirror image g-space trajectories. The mirrored trajectory (dashed black curve) is plotted on the reverse direction moon illumination and normalized gravity axes. Actual migratory trajectories in (A-C) are shown as blue-hued triangles that are color-coded and sized according to 121203.17's average daily velocity as shown in the legend. Normalized gravity values in (C) were determined by dividing the sum of the latitude dependent gravity (Götze, 2014) and the bedrock dependent gravity (Balmino et al., 2012) at each average daily location by the value present at 121203.17's Abrolhos Bank migratory departure site. Colored basemap in (A) depicts bedrock gravity anomalies across a 150 to 660 mGal range.
(2) the moon was within 72 h of its quasi-monthly (i.e., the 27.32 day duration lunar sidereal month; Supplementary Figure S3) lunar declination maximum of ∼28 • ; (3) the moon was within 24 h of its quasi-monthly (i.e., the 27.55 day duration anomalistic month; Supplementary Figure S3) apogee of ∼404,000 km. Second, 24641.05's migration describes a highly symmetrical gravitational coordinate trajectory with mirror planes through both the gravitational midpoint of its migration (i.e., 100.134% of its departure site's gravitational acceleration; Figure 8C) and the 29.5 day duration synodic month's 50% moon illumination position. Maintenance of this highly symmetrical trajectory required 24641.05 to swim both faster and slower when it experienced relatively higher or lower bedrock derived gravity anomalies, respectively ( Figure 8A). Third, 24641.05 ended its south-southeast migration, on a neap tide first-quarter moon, as it approached 46 • S latitude in the end of the first week of December, 2005. Similar patterns are present in the gravitational coordinate trajectories of all 9 slow motif humpback whales (Figure 7A), irrespective of the year in which the migration occurred or the prevailing geomagnetic and oceanographic conditions (Figure 6).

Fast Migration: PTT 121203.17
Humpback whale PTT 121203.17's south-southeast directed migration was relatively fast, averaging 5.0 ± 1.3 km h −1 ( Table 1 and Figure 9). Like the 10 other fast motif whales, 121203.17 started its migration in the days surrounding new or full moon, and it continued its non-stop and highly directional (straightness index = 0.95; average heading 161 ± 11.7 • ; Table 1) swimming for ∼3900 km during a ∼1.25 month long period. Humpback whale 121203.17 ended its migration on 26 November 2017 within 24 h of first-quarter moon ( Figure 9C) as it approached 54 • S latitude ( Figure 9B).
As was the case for PTT 24641.05, several rare conditions were met during 121203.17's relatively fast and highly directional migration that collectively determined the gravitational accelerations 121203.17 experienced. First, PLR-BPA suggests 121203.17 initiated its migration on 24 October 2017: (1) when the moon was <25% illuminated and the tidal gravity vector achieved a relatively moderate ∼180 µGal daily range (Supplementary Figure S4); (2) within 24 h of the sidereal month's lunar declination minimum of −19.4 • (Supplementary Figure S4); (3) within 24 h of the moon's anomalistic month's apogee (i.e., ∼405,000 km; Supplementary Figure S4). Second, 121203.17's migration describes a highly symmetrical gravitational coordinate trajectory with mirror planes through both the gravitational midpoint of its migration (i.e., 100.158% of its departure site's gravitational acceleration; Figure 9C) and the synodic month's 50% moon illumination (i.e., neap tide) position. As was the case with 24641.05, maintenance of this highly symmetrical trajectory required 121203.17 to swim both faster and slower when it experienced relatively higher or lower bedrock derived gravity anomalies, respectively ( Figure 9A). Similar patterns are present in the gravitational coordinate trajectories followed by all 11 whales in the fast motif ( Figure 7B), irrespective of the year in which the migration  Figure 6A (i.e., 9 slow motif whales) and Figure 6B (i.e., 11 fast motif whales), respectively. All four sinusoidal regression fits are highly significant (i.e., p << 0.05). occurred or the prevailing geomagnetic and oceanographic conditions (Figure 6).

Dynamically Paced Migration: PTT 120942.17
The migratory movements of one of the whales we tracked, PTT 120942.17, were particularly anomalous. Our results indicate that 120942.17 was the only whale in our PTT dataset to both: (1) perform a multiple day-long migratory stop-over; and (2) retrace its migration route with north-northwest swimming over a ∼10 day-long period (Figure 10). PTT 120942.17 is also the only whale to clearly shift from one velocity motif (slow) to the other (fast) during its staged southward migration ( Table 1).
Humpback whale PTT 120942.17's movements can be divided into four separate stages (Figure 10). Stage 1 includes its initial south-southeast directed migration from Abrolhos Bank to the northern flank of Rio Grande Rise (Figures 10A,C). Stage 2 includes its stop-over above Rio Grande Rise, where ocean floor altitudes ∼550 m below sea level are ∼4000 m higher than adjacent abyssal plains (Amante and Eakins, 2009). Stage 3 includes 120942.17's reverse north-northwest directed migration back toward Abrolhos Bank (Figures 10B,D), and stage 4 includes 120942.17's recovery of its initial south-southeast directed migration path (Figures 10B,D).
Despite the anomalous and staged structure of 120942.17's movements, the g-space trajectories described by these movements coincide with the phased sinusoidal trajectories followed by several other whales that also migrated through the SSEC at distinctly different times between 2003 and 2017 (Figures 7, 10). Specifically, stage 1 of 120942.17's migration ( Figure 10C) initially coincides with the phased movements of other whales in the slow motif (Figure 7A), including the g-space trajectory followed by PTT 24641.05, approximately 12 years earlier ( Figure 8C). Despite initially following a g-space trajectory similar to the one followed by PTTs 24641.05 and 96380.10 ( Figure 7A), 120942.17 subsequently departs from this slow motif trajectory as it approaches Rio Grande Rise and transitions into stage 2 of its migration (Figures 10A-C).
Importantly, PTT 120942.17 continued in an overall southward direction during stage 2 of its migration ( Figure 10A). Yet, the gravitational acceleration 120942.17 experienced during this ∼10 day period did not change ( Figure 10D). This apparent discrepancy between 120942.17's geographic and gravitational coordinate movements is easily explained. Although 120942.17 continued swimming southward (Figure 10B), the gravitational effects of this continued southward movement were counter balanced by local changes in gravity associated with the bedrock underlying Rio Grande Rise ( Figure 10A). These facts resulted in 120942.17 experiencing a constant gravitational acceleration throughout stage 2 of its migration ( Figure 10D). As a direct consequence of its movement behavior during stage 2 of its migration, 120942.17 is moving in geographic coordinates, but not moving in gravitational coordinates for precisely one-half of a tidal gravity cycle (i.e., from new moon to full moon; Figure 10D).
Finally, stage 4 of 120942.17's migration begins in the dawn twilight on a last-quarter moon. The associated g-space trajectory for stage 4 coincides with the phased sinusoidal trajectories followed by the remaining 5 whales in the fast motif (PTTs: 10946.05, 26712.05, 50687.07, 87775.12, 121192.12; Figures 7B, 10D).
Despite making multiple movement decisions that both delayed its southward migration by more than 24 days and added more than 2000 km of swimming distance, 120942.17's migration describes highly symmetrical and phased gravitational coordinate trajectories that coincide with segments of the g-space trajectories followed by 13 of the other 19 whales that migrated through the SSEC. The symmetrically patterned and reproducible structure of these gravitational coordinate trajectories implies that the benefits of swimming through a well-defined and heavily utilized spatiotemporal corridor outweigh the energetic and temporal costs of swimming farther for longer periods. The risks associated with such a strategy are not minor and include exposure to changing oceanographic conditions, shifts in prey availability (e.g., Figures 1, 6), heightened competition, and predictability.
Indeed, the predictability of the humpback whale migrations we report may have already hurt the population (Rosenbaum et al., 2009). The fact that Soviet whale ships killed humpback whales at the southeast corner of Abrolhos Bank, within the SSEC, and above Rio Grande Rise but crucially not elsewhere, over a 2-week period suggests there has been >50 years of spatiotemporal fidelity to the migratory corridor revealed by our satellite tracking research despite dynamic oceanographic and geomagnetic conditions. Plotting the time and location of the Western South Atlantic Soviet humpback whale kills in gravitational coordinates supports this interpretation: the 1967 Soviet whaling kills coincide with the sinusoidal gravitational trajectories followed by the humpback whales we tracked between 2003 and 2018 ( Figure 7B).

CONCLUSION
Using long-term satellite remote sensing data, our research provides empirical support for what several models have predicted: changes in Earth system conditions will challenge the sustainability of some populations of marine megafauna (MacLeod, 2009;Barbraud et al., 2011;Hazen et al., 2013;Silber et al., 2017). The prolonged spatiotemporal fidelity of humpback whale movements, despite contemporaneous oceanographic and geomagnetic change, suggests humpback whale movement decisions include mechanistic responses to stable and predictable exogenous cues, including gravity.
The relative stability of gravitationally derived cues helps explain the apparent paradox of humpback whale migratory fidelity despite pronounced oceanographic and geomagnetic change. Yet, the extent to which dynamic oceanographic conditions evoke changes in marine megafaunal long-distance migratory movement decisions, and the environmental and biogeophysical thresholds that trigger specific movement behaviors, remain unknown. Navigation during long-distance migration is likely to be informed by a diverse suite of cues, and if we are to ever truly know how whales navigate, all reasonable mechanisms must be considered as part of inclusive and integrated research on diverse species across a variety of environmental and biogeophysical contexts over prolonged periods. Documenting the extent to which animal migration routes, destinations, and movement decisions change, or perhaps fail to change, in a changing environment is essential to preserving biodiversity and sustainably managing diverse ecosystems and the services they provide.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
TH conceived the study, compiled all satellite remotesensing data, prepared the figures, performed the geophysical, astronomical, and statistical analyses, and wrote the initial manuscript. AZ, AA, DD, and FS performed all of the field research, including satellite tag preparation, deployment, and data management. All authors contributed to the writing and revision of the initial manuscript. The scientific results and conclusions, as well as any views or opinions expressed herein, are those of the author(s) and do not necessarily reflect those of NOAA or the US Department of Commerce.

FUNDING
Funding from Shell Brasil and CGG Brasil supported collection of the humpback whale satellite telemetry data we report.