Environmental drivers of large‐scale movements of baleen whales in the mid‐North Atlantic Ocean

Understanding the environmental drivers of movement and habitat use of highly migratory marine species is crucial to implement appropriate management and conservation measures. However, this requires quantitative information on their spatial and temporal presence, which is limited in the high seas. Here, we aimed to gain insights of the essential habitats of three baleen whale species around the mid‐North Atlantic (NA) region, linking their large‐scale movements with information on oceanographic and biological processes.


| INTRODUC TI ON
Migratory marine animals are increasingly vulnerable to extinction and population declines due to the cumulative anthropogenic impacts on their environment. Understanding the distribution patterns of these wide-ranging taxa is critical for effective management (Halpern et al., 2015), but often challenging, due to the difficulty of studying their movements over vast ocean regions. Advances in satellite telemetry tags have improved the knowledge on the distribution and movements of highly mobile marine species and stimulated developments in species distribution modelling (SDM) to quantify and predict animal space use across time (Block et al., 2011). However, making population-level inferences from telemetry data poses many challenges due to the presence-only observations, such as the lack of absence data and serially correlated positions.
One solution is to incorporate pseudo-absences to the data (i.e. representatively selected locations where the animals could have been sampled but were not) and to use an error structure (i.e. correlation terms, random effects) that more accurately reflects the variability within and between different levels of the natural hierarchy of the sampling units (Aarts, MacKenzie, McConnell, Fedak, & Matthiopoulos, 2008;Vandeperre, Aires-da-Silva, Lennert-Cody, Serrão Santos, & Afonso, 2016).
Developing SDMs that can effectively forecast animals' distribution depend on the availability of ecologically meaningful environmental data at the appropriate spatial and temporal resolutions.
Marine top predators are supported by the productivity of primary and secondary consumers, and their distribution and movements are tied closely to those of their prey (Benoit-Bird et al., 2013;Boyd et al., 2015). Yet, due to the scarcity of prey data at scales relevant for marine top predators, most SDMs have relied on proxies for prey distribution such as physiographic variables and a handful of remotely sensed or in situ oceanographic measurements (Scales et al., 2017;Tobeña, Prieto, Machete, & Silva, 2016). However, the main issue with these proxies is the potential temporal and spatial lags that may occur between physical and oceanographic processes and biological responses (i.e. prey biomass), highlighting the importance of having a good understanding of the marine food web to predict the distribution of predators (Grémillet et al., 2008;Redfern et al., 2006). Recent advances in ecosystem models provide real-time and hindcast simulations of low and mid-trophic level biomass and production over large spatio-temporal scales (Lehodey et al., 2015;Lehodey, Murtugudde, & Senina, 2010), which allow the incorporation of predicted prey distributions into SDMs and a better understanding of the ecosystem components (Lambert, Mannocci, Lehodey, & Ridoux, 2014;Mannocci, Roberts, Miller, & Halpin, 2016).
Large baleen whales are believed to make seasonal migrations from high-latitude productive waters used in summer to feed, and low-latitude oligotrophic waters used in winter for breeding and calving (Stern, 2009). During feeding periods, these whales concentrate mainly in areas where their preferred prey (euphausiids, small epipelagic and mesopelagic schooling fish) occurs in predictable aggregations and at high densities (Croll et al., 2005;Laidre et al., 2010).
In the NA Ocean, studies examining the drivers of baleen whale distribution are concentrated along the US and European coasts (Kaschner, Quick, Jewell, Williams, & Harris, 2012). In the northwest Atlantic, the spatio-temporal distribution of baleen whales was strongly related to static variables and several formulations of biological productivity, such as mesoscale fronts, and eddies, which can aggregate micronekton and zooplankton biomass . In the European Atlantic shelf waters, oceanographic and physiographic features, such as sea surface temperature, depth and distance to the 2,000 m contour, were important predictors of the summer distribution of baleen whales (Macleod et al., 2009).
The drivers shaping baleen whale migration in the mid-NA are poorly understood due to the difficulty of observing their movements during migratory periods. Identifying the drivers of baleen whales' distribution is necessary to understand the dynamics of their populations and to forecast how these whales might respond to ecosystem variability associated with climatic changes . Studies characterizing the migration patterns of these whales remain largely descriptive (Prieto, Silva, Waring, & Gonçalves, 2014;Silva, Prieto, Jonsen, Baumgartner, & Santos, 2013). Therefore, our aim was to develop habitat models based on telemetry data to investigate the distribution of three baleen whale species around the mid-NA: the blue whale (Balaenoptera musculus), fin whale (B. physalus) and sei whale (B. borealis). Specifically, we implemented a comprehensive ecological modelling approach to determine the effect of oceanographic conditions and prey availability as drivers of their northward migration. Based on these models, we also predicted the spatio-temporal distribution along their northward migration in order to identify important habitats for each species in the mid-NA.  Table S1). Tagging procedures were conducted following the guidelines of the American Society of Mammalogists (Sikes, Gannon deployed from a 6-m rigid-hulled inflatable or a 12-m fibreglass boat using a compressed air gun (ARTS/RN, Restech, Norway) set at 8-10 bar pressure (Prieto et al., 2014;Silva, Prieto, Jonsen, et al., 2013). All tags were programmed to transmit on a daily basis, providing hourly data up to a maximum of 500 daily messages. The raw and unfiltered locations provided by the tags were fitted in a Bayesian switching state-space model (SSSM) to obtain positions at regular time steps, and to improve position estimates by incorporating the measurement errors provided by Argos . The SSSM was built in a hierarchical framework to estimate parameters jointly across multiple individual tracks. We tested multiple time steps (2, 5, 12 and 24 h) to find the most appropriate parameters to reduce the sample autocorrelation (see modelling approach workflow in Figure S1). For each model, we ran two Markov chain Monte Carlo (MCMC) chains for 50,000 iterations, dropping the first 45,000 samples as a burn-in and retaining every 5th sample from the remaining 5,000 assumed postconvergence samples.

| Tracking data
Model convergence and sample autocorrelation were visually checked in the diagnostic plots, and a time step of 12 h was selected as the most suitable due to its lower sample autocorrelation.

| Simulated correlated random walks
As the track locations obtained through the SSSM only provide information on the presence of whales, we created simulated tracks (i.e. "pseudo-absences") for each original track (case points) to characterize environmental conditions potentially available to individuals (Aarts, Fieberg, & Matthiopoulos, 2012). The simulated tracks (control points) were built using a correlated random walk (CRW) model (Kareiva & Shigesada, 1983), randomly taking the first or last transmitted position of each satellite track as the starting point to avoid bias towards areas more intensively sampled (Phillips et al., 2009; see Figure S2). For each real track, we

| Environmental and prey-related variables
Candidate variables were selected based on species' ecological preferences and regional oceanographic and physiographic characteristics (Redfern et al., 2006). Potential prey biomass distributions were obtained from the mid-trophic level spatial ecosystem and population dynamics model (SEAPODYM-MTL, Lehodey et al., 2015;Lehodey et al., 2010; Table 1). The SEAPODYM-MTL model simulates spatial and temporal dynamics of production and biomass of micronekton between the surface and 1,000 m depth.
The model is driven by physical and biological variables, using a system of advection-diffusion-reaction equations (Lehodey et al., 2010). Physical variables, current and temperature (T) data (Table 1) and euphotic depth (ZEU) data, were derived from ocean colour satellite data (http://www.scien ce.orego nstate.edu/ocean ), using the vertically generalized production model (VGPM) of Behrenfeld and Falkowski (1997). Both GLORYS-2v1 outputs and primary production data were interpolated onto a regular 0.25° × 0.25° grid with a weekly time step to be used as forcing for the SEAPODYM.
In the model, three vertical layers are defined in relation to the euphotic depth, as micronekton diel vertical migration (DVM) is mainly induced by daylight variations, validated with biomass estimates from micronekton sampling cruises and acoustic backscatter data (Lehodey et al., 2015(Lehodey et al., , 2010. The SEAPODYM-MTL model simulates six functional groups of micronekton according to the dynamics of their vertical distribution: epipelagic, migrant upper mesopelagic, upper mesopelagic, highly migrant lower mesopelagic, migrant lower mesopelagic and lower mesopelagic ( Table 1). As baleen whales mainly forage on prey available within the epipelagic layer (Goldbogen et al., 2017), functional groups remaining during day and night within the lower mesopelagic layers (>400 m deep) were not included in the analysis.
We also incorporated the outputs of the lower trophic level SEAPODYM (SEAPODYM-LTL), which follow the same modelling framework of the mid-trophic level, and considered a single functional group including all mesozooplankton organisms (i.e. both holo-and mero-zooplankton). This model also had two state variables, production (LTL-prod) and biomass (LTL-bio), integrated over the epipelagic layer defined by the euphotic depth (Conchon, 2016) ( Table 1). In addition to these prey-related variables, we included the surface chlorophyll-a concentration (Chl-a) used to compute the primary production and eddy kinetic energy (EKE) calculated as EKE = ½ × (u2 + v2), with u and v, respectively, the meridian and zonal geostrophic current components that are used as forcing variables for the epipelagic layer of SEAPODYM.
Regarding static environmental variables, depth was obtained from ETOPO1 Global Relief Model at a spatial resolution of 1/60°, provided by NOAA's National Geophysical Data Center (NGDC, https ://www.ngdc.noaa.gov/mgg/globa l/). Derived static variables were slope, calculated using the elevation values of the four immediately adjacent depths (Ritter, 1987), and distance to coast (See ST1).
Following the identification of time-lagged Chl-a concentration as an important predictor of the spatio-temporal distribution of baleen whales around the study area , we tested a time-lag of 1 month prior to the point location for the following variables (indicated by variable name_1m): LTL-bio, LTLprod, Chl-a and NPP. All environmental and prey-related variables were extracted for each original and simulated track point using the R software and the packages "raster" (Hijmans & van Ettern, 2012), "ncdf4" (Pierce, 2019) and "rgdal" (Bivand et al., 2019). Data points with incomplete information were removed from the final dataset, corresponding to 0.02%, 0.06% and 1.60% of the total dataset for fin, blue and sei whale, respectively. Prior to running the models, we examined collinearity between pairs of covariates (Zuur, Ieno, & Elphick, 2010) using the Pearson's correlation coefficient and selected those variables less correlated (<0.7) and ecologically more relevant (Dormann et al., 2013).

| Habitat modelling
We applied generalized additive mixed models (GAMMs) to the original and simulated track points to predict the habitat preference of each whale species. GAMMs were fit with a binomial distribution and a logit link function, using the restricted maximum likelihood (REML, mgcv 1.8.7; Wood, 2006) as an optimization method and individual IDs as a random effect. To model our data as a binary response variable, original points were assigned a value of 1 and simulated points a value of 0 (Aarts et al., 2012). In addition to single covariates, we defined two biologically interpretable interactions between covariates using two-dimensional tensor product smooths: T × LTL-bio and EKE × LTL-bio (Wood, Scheipl, & Faraway, 2013).
These interactions were chosen based on current knowledge of the oceanographic conditions, as well as the relationship between these variables to understand the habitat preference of marine top predators within our study area (Vandeperre et al., 2016).
The GAMMs were built using a backward selection of relevant individual smooth terms and tensor product smooths. Nested models were compared using chi-square tests. To select the tensor product smooths, the original structure was tested to the equivalent additive structure while controlling for the degrees of freedom using the k-parameter (i.e. the basis dimension) (Vandeperre et al., 2016;Wood, 2006). Once the best model was chosen based on the lowest corrected AIC, the smoothing splines were limited to a maximum of 5 degrees of freedom to avoid overfitting and preserve the ecological interpretability of functional relationships. The effect of this limitation on the number of degrees of freedom was evaluated by comparing the final outputs and predictions from models with constrained and unconstrained smooth terms. No differences were detected between constrained and unconstrained models on variable selection and predictive performance (results not presented).
Following this habitat modelling approach, we evaluated the performance of models built with both environmental and prey-related variables, with that of models including only environmental 0.25°SEAPODYM ecosystem model* Production of mesozooplankton within the epipelagic layer Note: The SEAPODYM defined three vertical layers in relation to the euphotic depth: Epipelagic layer as <1.5 of the euphotic depth (~150 m), upper mesopelagic layer as >1.5 and <4.5 of the euphotic depth (between 150 and 400 m) and lower mesopelagic layer as >4.5 of the euphotic depth (>400 m). Examples of species composition of functional groups are provided in Lehodey et al. (2010).
*is to describe the seapodym ecosystem model variables, in order to assess the contribution of modelled mid-and lower trophic level biomass and production.
A sensitivity analysis was also carried out to determine whether the simulated CRW tracks produce outputs that were robust to parameter uncertainty. We randomly selected two of the 200 CRWs tracks per whale and assessed whether the predictors from the best-fitting model were also significant for the different simulated points, repeating this procedure 40 times (Hazen et al., 2016). For each of the 40 runs, we calculated the area under the receiver operating characteristic (AUC) curve and the variance explained (R 2 ). The AUC index ranges from 0 to 1, with values below 0.6 indicating a discrimination ability no better than random, values between 0.7 and 0.9 considered as reasonable and values >0.9 as excellent (DeLong, DeLong, & Clarke-Pearson, 1988 (Žydelis et al., 2011). GAMMs fit to the subsampled dataset were compared with the standard models through variable importance and AUCs.

| Model predictions
Monthly predictions of habitat preferences (scale from 0 to 1) were produced for each species based on the best GAMMs, at the spatial resolution of the fitted models (0.1°). Environmental and preyrelated variables were compiled at a monthly resolution during the tracking period to obtain mean predictions for the study area and calculate the standard deviation (SD), in order to measure the stability of the predicted distribution, with stable and unstable habitat represented by low and high SD.

| Tracking data
Tracking data were collected for the three target species between March and July, from 2008 to 2016. Two fin whale tracks shorter than 3 days were excluded from the analysis (final dataset of n = 29;

| Fin whale
The best-fitting GAMM included two additive covariates (Depth and log_EKE) and one tensor product smooth (LTL-bio × T) ( Table 3; Table   S2). Fin whales preferred water depths < 2,500 m with high EKE (> 0.005 cm 2 /s 2 ) ( Figure 2). The surface of the tensor product smooth of T and LTL-bio was plotted on the scale of the response. The tensor product smooth indicated an association with cold (<10°C) productive habitats (LTL-bio > 100 gWW). Habitat preferences decreased sharply at T values between 10 and 15°C, but increased again for warmer waters with low LTL-bio. The environmental and prey-related variables of best model were significant in the large majority of the 40 models ran with different simulated points, showing the robustness of these variables to control point selection (Table 1)  of uncertainty across the study region from March to May, with this uncertainty remaining in the southern region of the study area in June and July ( Figure S3).

| Blue whale
The best model predicting blue whales' habitats included one tensor smooth product (LTL-bio × T) and two additive covariates (Depth and LTL-prod_1m) ( Table 3; Table S2). The tensor smooth product of T and LTL-bio defined a complex interaction surface (Figure 4) Note: Models were run 40 times, and the significance of their variables was evaluated with p values of less than 0.01. Variances explained (R 2 ), Akaike information criterion (AIC) and area under the curve (AUC) were obtained for each run, calculating a mean and SD. Variables included in the best models below: lower trophic level biomass (LTL-bio), lower trophic level production on the previous month (LTL-prod_1m), eddy kinetic energy (log_EKE) and depth (Depth). in June, and in the south-west coast of Iceland by July ( Figure 5).

TA B L E
The probability of occurrence of blue whales in the Azores decreased from April to July, reaching almost zero in July. Predictions showed a fair amount of uncertainty in the northern part of the study area in April, but decreased drastically in the following months ( Figure S4).

F I G U R E 3
Predicted distribution of fin whales from March to July. Habitat preferences were scaled from 0 (low preference) to 1 (high preference). Spatial resolution was set to 0.1°F I G U R E 4 Generalized additive mixed model (GAMM) smoothers for blue whales. (a) Surfaces of the tensor product smooths (temperature, T, and lower trophic level biomass, LTL_bio) plotted at the scale of the response. (b) and (c) GAMMs smother for log_LTL_prod_1m and Depth, respectively. Dashed lines represent the 95% confidence intervals for the fitted relationships

| Sei whale
The final model for sei whales contained one tensor smooth product (LTL-bio × T) and two additive covariates (Depth and LTL-prod_1m) ( Table 3; Table S2). Sei whales were associated with deep waters (>1,500 m) and intermediate values of LTL-prod from the previous month (25-65 gWW/day) ( Figure 6). The tensor smooth product of T and LTL-bio indicates that sei whales used relatively cold waters (~10°C) with average productivity values (10-100 gWW), as well as colder areas (~5°C) with high productivity (LTL-bio > 100 gWW).
Preference decreased sharply for warmer areas (>15°C), even when LTL-bio was high. The tensor smooth product was significant in 100% of the 40 models ran with different control points, in 80% of the models with the covariate Depth and in only 20% of models with the LTL production from the previous month. The average AUC was 0.83 (SD = 0.04), and the average deviance explained was 31% (SD = 0.09). The best model with environmental variables included four covariates (temperature, EKE, Depth and Slope) and had a lower predictive performance and variance explained than models with prey-related variables (Table S3).
Latitudes above 45°N, mainly around the south-west of the Irminger Sea and in the Labrador Sea (except the shallow waters along the Greenland coast), were identified as important habitats for sei whales (Figure 7). In May, the main hotspot was near the Charlie-Gibbs Fracture Zone (CGFZ) in the mid-NA. In June and August, the predicted distribution shifted northwards and was concentrated in the entire Labrador Sea and in a few areas along the Greenland coast. Low occurrences were predicted in the area around Azores.
Predictions indicated a low uncertainty across the whole study region in all sampled months ( Figure S5).

| Serial correlation evaluation
Models based on subsampled and full datasets yielded similar results for all three species in terms of importance of predictor variables and robustness to control point selection (Table S4). In the case of fin and blue whales, subsample models had slightly lower AUC values than full dataset models but AUC values were similar for sei whale models. Nonetheless, habitat preference predictions from these models were comparable to those from full datasets for all the species (Figures S6-S8). Taken together, these results indicate that our analyses and predictions were robust to serial correlation.

| D ISCUSS I ON
This study offers novel insights into the distribution of highly migratory species across large geographical scales, and fills in a gap in knowledge of baleen whale distribution in the mid-NA. We provide spatially and temporally explicit information on the distribution of the three species of baleen whales through the combination of satellite tracking data, remotely sensed oceanography and ecosystem models, which could contribute to implement appropriate conservation and management measures. However, it should be noted that our analysis is based on a relatively small sample size, especially for sei whales, and that some of the tracks were short and did not cover the full extent of movements to the whales' presumed northern feeding grounds. Thus, our habitat preferences may not adequately reflect the relationship of whales with conditions at their feeding grounds.
Moreover, inter-monthly variations may not be fully representative F I G U R E 5 Predicted distribution of blue whales from April to July. Habitat preferences were scaled from 0 (low preference) to 1 (high preference). Spatial resolution was set to 0.1° of the seasonal movements of the species due to the limited number of tracks per month within each sampling year.
Regarding prey-related variables, SEAPODYM data also had some issues related to data estimation and validation. First, various sources of uncertainty are associated with the observations used to validate the model, due to the different techniques employed to collect zooplankton biomass data and its spatio-temporal distribution (Conchon, 2016). Second, the LTL data are based on mesozooplankton species, which constitutes only a portion of the diet of our targeted baleen whale species. While further developments are required to optimize final output parameters, this study shows promising results on the application of modelled zooplankton biomass to study marine top predators at large scale.

| Habitat preferences
Our models show that the distribution of blue, fin and sei whales is strongly associated with prey availability and with static and F I G U R E 6 Generalized additive mixed model (GAMM) smoothers for sei whales. (a) Surfaces of the tensor product smooths (temperature, T, and lower trophic level biomass, LTL_bio) plotted at the scale of the response. (b) and (c) GAMMs smother for log_LTL_prod_1m and Depth, respectively. Dashed lines represent the 95% confidence intervals for the fitted relationships F I G U R E 7 Predicted distribution of sei whales from May to August. Habitat preferences were scaled from 0 (low preference) to 1 (high preference). Spatial resolution was set to 0.1° dynamic environmental variables closely linked to productivity.
However, habitat preferences varied among species and across the distribution range of each species. Fin whales exhibited a bimodal response to environmental variables along the northward migration route, with a marked preference for cold waters associated with high LTL-bio levels and for warmer waters, and avoidance of intermediate temperatures. Such bimodal response reflects differential exploitation of distinct water masses and productivity regimes at the ends of their studied range, one characteristic of high productivity in cold waters and the other of low productivity in temperate waters. Sei whales also showed a clear preference for productive subpolar waters in favour of temperate ones, displaying a two-phase distribution pattern during the study period, consistent with limited available observations (Sigurjónsson, Gunnlaugsson, & Payne, 1989;Waring, Nottestad, Olsen, Skov, & Vikingsson, 2008). In comparison, blue whales had a more homogenous response across the study area, targeting areas with intermediate LTL-bio values that were productive on previous months. The retention of T and LTL-bio, and the different response obtained by the three species suggest that interactions between environmental and prey-related covariates may help to predict species distribution. Past studies showed the utility of proxies such as SST and productivity to predict the occurrence, abundance and species richness of marine top predators. Although these predictors might indicate important foraging areas for these species, they were probably decoupled in time and space, and from prey availability at higher trophic levels (Grémillet et al., 2008). Incorporating the novel LTL data allowed us to shorten the temporal lags between the predictors and the predators, helping to reduce the possible spatial mismatch between predictor and response variables and representing a more direct biological relationship.
The habitat preference for the three species was significantly constrained by depth, which has been shown to be a good predictor of baleen whale distribution in multiple regions (Scales et al., 2017), including in the NA (Víkingsson et al., 2015). Blue and sei whales preferred oceanic areas and avoided coastal regions, whereas fin whales used both coastal and oceanic areas.
In addition, fin whale distribution was strongly influenced by EKE, a measure of mesoscale activity. Upwelling/downwelling processes are common in areas with intense mesoscale activity, such as eddies and fronts, enhancing marine productivity and leading to the formation of predictable prey patches (Worm, Sandow, Oschlies, Lotze, & Myers, 2005). Numerous studies have linked these features to the presence of marine top predators, including fin whales (Reisinger et al., 2018;Scales et al., 2017). While fronts and eddies are very dynamic in the mid-NA and can have significant interannual variations, their persistence around the south-west Irminger Sea and CGFZ has been well documented (Volkov & Oceanic, 2005). Due to persistent thermal fronts, these latest regions are well known for being extensively exploited by a wide range of marine top predators .
Differences in habitat preferences among whale species should partly reflect differences in the distribution of their main prey. In the NA, sei whales primarily feed on calanoid species, especially on Calanus finmarchicus, and to a lesser extent on euphausiids (Sigurjónsson & Vikingsson, 1997). Blue whales feed almost exclusively on krill (e.g. Meganyctiphanes norvegica and Thysanoessa spp) (Gavrilchuk et al., 2014;Víkingsson, 1990), and fin whales feed on zooplankton and a variety of small schooling fish (Sigurjónsson & Vikingsson, 1997). Therefore, the wider dietary breadth of fin whales likely contributes to higher plasticity in foraging behaviour and enables this species to exploit a broader set of environmental conditions across their migration route.

| Predicting spatio-temporal distribution across the NA
Unsurprisingly, our models showed that high latitudes were the preferred habitats for all three baleen whale species during summer months, in accordance with known seasonality and abundance patterns. However, our results also highlight that their predicted summer habitats are not totally overlapping due to distinct habitat preferences. The prediction maps for fin whales dis-  Waring et al., 2008). These findings agree with our predictions for the species during this month. In July, the preferred habitat for blue whales expanded north, showing a higher preference for the eastern NA. Passive acoustic monitoring off north-western Britain and Ireland (55°N 20°W) showed few blue whale detections in April and June, with a gradual increase in mid-July until reaching a peak in November-December (Charif & Clark, 2009). In recent years, visual surveys revealed the presence of large baleen whales during these months in the area, including blue whales, mainly around the highly productive waters of the Porcupine Seabight and Rockall Trough (Baines & Reichelt, 2017;Wall et al., 2013;Wall, O'Kelly, Whooley, & Tyndall, 2009). Our predictions in northern latitudes, close to west Iceland, largely coincide with the known distribution of this species in summer, where a population of around 1,000 blue whales was estimated. In the Azores, the predicted shift towards northern latitudes closely matches the seasonality in the sighting data, which shows a decrease from April to July . During these months, blue and fin whales have a similar environmental niche with a high degree of spatial coincidence in this region, strongly influenced by the higher productivity in the region .
The highest habitat predictions for sei whales occurred in areas above 45°N, from May to August. These findings agree with previous reports of the known occurrence of feeding sei whales in these waters (Anonymous, 2009;Mitchell, 1975;Sigurjónsson et al., 1989). In the eastern NA, the deep waters of the Irminger Sea, between east Greenland and west of Iceland, appear to be an important habitat for the species, while the Faroe-Shetland Channel and Norwegian Sea are used less intensively (Cattanach, Sigurjónsson, Buckland, & Gunnlaugsson, 1993). In the western NA, sei whales seem to concentrate in the northern banks at 66°N and off the southern tip of Greenland, as well as along the continental shelf waters of the Labrador and Newfoundland coasts.
Recently, Roberts et al., (2016)  Our monthly predictions for sei whales displayed a northward expansion of their habitat preference from May to August. In May, the preferred habitat for the species occurred around the CGFZ.
In this region, the interaction of the NA current with the mid-Atlantic Ridge generates intense persistent fronts leading to predictable foraging grounds for marine top predators (Doksaeter, Olsen, Nøttestad, & Ferno, 2008). A high concentration of sei whales reported just north of the CGFZ was associated with the fine-scale processes interacting in the upper 100 m of the water column Waring et al., 2008). Higher zooplankton biomass, particularly of Calanus finmarchicus, was found in the northern region of the CGFZ compared to the southern region, which is believed to be related with retention/advection processes concentrating large zooplankton over steep topography (Gaard et al., 2008). From June to August, sei whales were predicted over the same habitats in the deep waters of the Irminger Sea and Labrador Sea, including a few hotspots along the west Greenland. Laidre et al., (2010) found that krill biomass, mainly Meganyctiphanes norvegica and Thysanoessa sp, was the most important variable explaining the presence of baleen whales in west Greenland. Although the number of sei whale sightings (n = 17) was limited, sighting aggregations matched well the areas with the highest densities of krill. In Azores, sei whales are seen every year, with most sightings from spring to early summer (Silva, Prieto, Cascão, et al., 2013). However, sei whale presence has not been linked to local productivity. This matches well with the lack of association with LTL biomass and the low preference for habitats in these latitudes. Although the species feeds occasionally around Azores, our results indicate that the archipelago is not an important feeding ground for sei whales, at least during spring and summer.

| Management and conservation implications
Highly migratory species are exposed to increasing impacts from anthropogenic activities, which can lead to the alteration of their habitats, disruption of natural behaviours or increased mortality (e.g. from ship strikes) (Maxwell et al., 2013). In the mid-NA, a major international shipping route located mainly between 35° and 45°N overlaps with the predicted preferred habitat of the three baleen whale species. Specifically, blue whales have the highest proportion of the predicted habitat overlapping with this main route, potentially increasing the risk of ship collisions, as well as exposing them to increased noise levels. Alongside these chronic forms of noise pollution, the NA has also experienced an intense warming over the last decades, which is projected to continue over this century (IPCC, 2007). These physical changes are anticipated to lead to a poleward shift of several species of copepods (Villarino et al., 2015). As a consequence, baleen whales might move northwards to track the movements of their prey, similar to predicted habitat shifts in the North Pacific . These shifts may decrease usage of middle latitude habitats by baleen whales. Moreover, species with a specialized diet, such as the blue whale, may have less capacity to adapt to these climate-induced changes compared to generalist predators, such as fin whales, which may be able to switch to other prey species (Pauly, Trites, Capuli, & Christensen, 1998

| CON CLUS IONS
Our study provides the most detailed investigation into the spatiotemporal distribution of large baleen whales in the mid-NA. Our results demonstrate that the movements of all three baleen whale species towards and at high latitudes were closely linked to prey availability, although species exhibited distinct ecological niches.
These results support previous work for these species, suggesting that baleen whales adjust their distribution during their migrations to ocean resources (Hazen et al., 2016). In our study, the habitat preferences of baleen whales were strongly related to the LTL biomass and production within the epipelagic layer, consistent with the diet of the species, and highlighting the strength of combining spatially explicit and vertically resolved prey information with oceanographic and telemetry data. Previous studies suggested the importance of including LTL biomass to develop climatological habitat models in the north-west NA using visual observations . Our study confirms the relevance of this explanatory variable for modelling the distribution of cetaceans at the scale of the NA. Future research should focus on overcoming some of the limitations of this study to develop cetacean distribution maps in near real-time to provide a more accurate and efficient conservation tool.

ACK N OWLED G EM ENTS
We are very grateful to Cláudia Oliveira, Irma Cascão, Maria João Cruz, Miriam Romagosa and many volunteers, skilled skippers, crew and spotters that participated in the tagging fieldwork. This work was supported by Fundação para a Ciência e Tecnologia