Humpback whale migrations to Antarctic summer foraging grounds through the southwest Pacific Ocean

Humpback whale (Megaptera novaeangliae) populations typically undertake seasonal migrations, spending winters in low latitude breeding grounds and summers foraging in high latitude feeding grounds. Until recently, a broad scale understanding of whale movement has been derived from whaling records, Discovery marks, photo identification and genetic analyses. However, with advances in satellite tagging technology and concurrent development of analytical methodologies we can now detail finer scale humpback whale movement, infer behavioural context and examine how these animals interact with their physical environment. Here we describe the temporal and spatial characteristics of migration along the east Australian seaboard and into the Southern Ocean by 30 humpback whales satellite tagged over three consecutive austral summers. We characterise the putative Antarctic feeding grounds and identify supplemental foraging within temperate, migratory corridors. We demonstrate that Antarctic foraging habitat is associated with the marginal ice zone, with key predictors of inferred foraging behaviour including distance from the ice edge, ice melt rate and variability in ice concentration two months prior to arrival. We discuss the highly variable ice season within the putative foraging habitat and the implications that this and other environmental factors may have on the continued strong recovery of this humpback whale population.


PTT
Spatial gradient in the rate of ice retreat.
For each season a spatial raster was created of the day of ice melt (i.e., the day of season any pixel was last covered by ice of at least 15% concentration). Data were derived from the daily NSIDC SMMR-SSM/I sea ice concentration data and reprojected to a 0.23° longitude-latitude grid. To capture negative gradients (e.g. early melt close to the Antarctic continent in the south) a first difference calculation 25km was applied north-to-south. Values are square root transformed (with sign retained).
Ice cv lag 2 (ice cv lag 2) Coefficient of variation (CV) in the sea ice concentration two months prior We used the coefficient of variation (CV = σ/µ), i.e. the ratio of the standard deviation to the mean, as an indicator of the ice edge variability. Where there is generally low mean ice concentration and high variance there will be a high CV giving a spatial representation of the marginal ice zone. The CV for each month was calculated from AMSR-E daily sea ice concentration at 6.25km resolution. So for example a tracking location in March would be associated with a CV (lag 2) value for that spatial pixel from January.

Mean sea ice concentration (%) one month prior
Mean ice concentration (%) for the month prior. Calculated from AMSR-E daily sea ice concentration at 6.25km resolution 6.25km Table S2.2. Exploratory model series examining environmental predictors of whale movement behaviour (n = 1442). The discrete behavioural state estimate (b, where 1 = 'transit', 2 = 'search') is modelled as a binomial response (i.e. the GAMM was fit using a binomial family with a logit link), with individual whale (id) included as a random effect. The GAMM allows for the relationships between behaviour and environmental predictors to be flexible and non-linear. The final GAMM assessed during the resampling procedure is highlighted in bold. All statistical models were fit using the gamm4 library in R (Wood and Scheipl, 2017). Initially the variables describing mean ice concentration at both one and two months prior (ice mn lag 1, ice mn lag 2) as well as the CV of ice concentration one and two months prior (ice cv lag 1, ice cv lag 2) were included in exploratory analyses. The one and two month terms were highly correlated (r > 0.8) therefore only one term describing ice concentration lag and one term describing ice variability lag could be included in the final model. The AIC values guided model selection, although the log CHLa (∆AIC +0.5) was retained in order to include at least one predictor representing biological productivity during the resampling procedure.

S2.4. Examining for effects relating to spatial and temporal autocorrelation.
Spatial. Spatially autocorrelated data can affect the smoothness estimation of model terms; i.e., the estimation procedure in GAM is likely to under-smooth. To further examine for effects relating to spatial autocorrelation, we therefore also implemented the same model including a spatial term interaction. Spatial coordinates (i.e. longitude and latitude) were transformed prior using a Lambert equal area transformation, and included in the model as explanatory variables via a full tensor product smooth on x and y (i.e., t2(x,y)). In a second implementation these were included as separate smooth terms (i.e., s(x) and s(y)).
In both cases the GAMMs retained the significance of the ice variables representing distance to the ice edge and ice melt rate; despite, as to be expected, spatial smooths absorbing much of the variation. There was some exchange in term significance between the lagged ice variables such that the mean ice concentration prior was prioritised rather than the ice variability in both spatial models. This indicates that the influence of the ice variability term likely mimics/is mimicked by the spatial term; i.e., it is likely to be associated with the focussed spatial area 145 -175°E around the Balleny Islands (66°55′S 163°45′E) where the majority of foraging behaviour was identified. Overall, these results reinforce the prominent role that features of the marginal ice zone play in influencing the search behaviour of humpback whales.
i) Full tensor product smooth on x and y (t2(x,y)
As to be expected, with more destructive resampling comes a loss of information and a degradation of the strength of the signal across all environmental variables. Even so, the ice related variables continue to be the most commonly retained across models. As such, the environmental interpretation that inclusion of the ice-related variables was the most highly resilient to uncertainty, appears robust.
Subsampling results. P-values are given for the original GAMM fit to the most probable discrete behavioural state estimate from the hSSSM, and indicate the approximate significance of GAMM smooth terms based on Chi-sq. statistics. Second column (N sign. /100) indicates the number of times each environmental variable was determined to be significant (based on p-value < 0.05) under the resampling procedure (n = 100 iterations), which resamples randomly from the posterior the estimates for both behavioural state and location.