Defining priority areas for blue whale conservation and investigating overlap with vessel traffic in Chilean Patagonia, using a fast-fitting movement model

Defining priority areas and risk evaluation is of utmost relevance for endangered species` conservation. For the blue whale (Balaenoptera musculus), we aim to assess environmental habitat selection drivers, priority areas for conservation and overlap with vessel traffic off northern Chilean Patagonia (NCP). For this, we implemented a single-step continuous-time correlated-random-walk model which accommodates observational error and movement parameters variation in relation to oceanographic variables. Spatially explicit predictions of whales’ behavioral responses were combined with density predictions from previous species distribution models (SDM) and vessel tracking data to estimate the relative probability of vessels encountering whales and identifying areas where interaction is likely to occur. These estimations were conducted independently for the aquaculture, transport, artisanal fishery, and industrial fishery fleets operating in NCP. Blue whale movement patterns strongly agreed with SDM results, reinforcing our knowledge regarding oceanographic habitat selection drivers. By combining movement and density modeling approaches we provide a stronger support for purported priority areas for blue whale conservation and how they overlap with the main vessel traffic corridor in the NCP. The aquaculture fleet was one order of magnitude larger than any other fleet, indicating it could play a decisive role in modulating potential negative vessel-whale interactions within NCP.

www.nature.com/scientificreports/ For the endangered Eastern South Pacific (ESP) blue whale (Balaenoptera musculus) population, northern Chilean Patagonia (NCP) is regarded as its most important summer foraging and nursing ground [10][11][12] . Previous studies on blue whale occurrence and movement patterns indicated that until the onset of austral autumn/winter migration, blue whales focus most of their activities within these productive coastal waters [12][13][14][15] . However, variations in how this population utilizes this region and other areas within the ESP appear to result from changes in prevailing oceanographic conditions 16 .
Species distribution models (SDM) have shown that austral spring chlorophyll-a concentration, prior to the whales' arrival, and thermal fronts are important oceanographic proxies for describing the abundance and distribution patterns of blue whales within the NCP 16 . Krill, the primary prey of blue whales 17 , can take advantage of seasonally enhanced productivity for biomass production, with some time lag linking early life-history stages (e.g. larval recruitment) with adult densities [17][18][19][20][21] . Adult krill biomass is subsequently concentrated by thermal fronts into high-density patches which blue whales prey upon [22][23][24][25] . This prey aggregation effect driven by thermal fronts could be critical for blue whales, and other large baleen whales, given their energetically costly feeding behavior [26][27][28][29] . We hypothesize that both time-lagged distribution of primary productivity and thermal front aggregating effect generates foraging conditions for blue whales within NCP. To further test predictions from this hypothesis, here we propose that individual blue whales modify their behavior within areas of high spring chlorophyll-a concentrations and/or thermal front occurrence. As foraging behavior cannot be directly assessed solely by inspecting tracking data, we consider area-restricted search behavior (ARS, lower velocity and less directional persistence) as a proxy for this type of behavior 30,31 .
Potential local threats affecting blue whales in NCP include collisions with vessels due to intense maritime traffic 16,32 , negative interactions with aquaculture and fisheries activities [33][34][35] , direct and indirect effects from poorly regulated whale-watching operations 36 , and general disturbance from noise and acoustic pollution 37 . As such, identifying priority areas for focusing conservation actions is of utmost relevance considering a population numbering the low hundreds with a very low potential biological removal from anthropogenic origin estimated at 1 individual every 1.8 years 16 for continued growth.
Vessel collisions with cetaceans have become recognized worldwide as a significant source of anthropogenic mortality and serious injuries [38][39][40][41] . Empirical work on this issue has been conducted in a few areas and populations, mostly in the northern Hemisphere 32,39,42,43 , with little effort conducted in South America 32,44 . In most countries, unreported cases, limited monitoring and insufficiently documented incidents have precluded any accurate assessment of the true collision prevalence and trend analyses 32 .
Given the earlier results from SDMs, we considered using telemetry data as a complementary tool for improving our understanding of blue whale habitat selection process 16,17,45 and investigating overlap with vessel traffic in NCP. In fulfilling these goals, here we provide: i) a novel fast-fitting model application for data gathered from satellite-monitored Argos tags (hereafter Argos tags), ii) model-derived spatial predictions of how whales use the area based on prevailing oceanographic conditions during the tracking period, iii) spatial estimates on the relative probability of encountering blue whales, based on the integration of movement model predictions with those of a previous SDM, and iv) spatial estimates on the relative probability of whales encountering vessels as a measure of risk for four different vessel fleets operating in NCP.

Methods
Study area. The NCP (41-47°S) is characterized by an intricate array of inner passages, archipelagos, channels, and fjords enclosing roughly 12,000 km of convoluted and protected shoreline (Fig. 1). Primary biological productivity here is modulated by the mixing of sub-Antarctic waters, rich in macro-nutrients, and the abundant input of freshwater (derived from river discharges, heavy precipitation and glacier/snow melt), rich in micronutrients, particularly silica [46][47][48] . Within the NCP, several micro-basins have been described, some of them having particularly high seasonal primary and secondary production [46][47][48][49][50] , providing resources that upper-trophic level species rely on 12,17,[50][51][52][53][54] . The area also hosts one of the largest salmon aquaculture industries in the world, among other anthropogenic activities that negatively affects local biodiversity 33,34,55 . Tagging and telemetry data. Argos tags were deployed on 15 blue whales during the austral summer and early autumn at their summering grounds off the NCP (Fig. 1), following procedures described elsewhere 14  Raw Argos data included locations within NCP and outside the area after the onset of migratory movement. Because we were concerned with understanding movement patterns within the NCP, we applied a cut-off point of 24 h prior to a clear sign of migration was observed. This subset of the data was filtered using the R package "argosfilter" 57 removing relocations that comprised velocities exceeding 3 m s −1 , this upper limit was defined based on previous maximum speed assessments for this population 14 .

Oceanographic covariates.
Chlorophyll-a and sea surface temperature (SST) data were extracted using R package "rerddapXtracto" 58  www.nature.com/scientificreports/ spring (DAHCC), defined as the distance to polygons enclosing areas with an average chlorophyll-a concentration equal or higher than 5 mg/m 3 during austral spring months (September, October, November), was the best explanatory variable in a SDM applied to line-transect survey data for blue whales in NCP 16 . Here we used the same procedure to construct this covariate but used the 95 th percentile of each year´s concentrations distribution within the study area as the cut-off point for defining areas of high chlorophyll-a concentration. This was preferred because whales might select areas with the highest productivity regardless of their absolute values. Maps for DAHCC were created for each year where telemetry data were available, and their values were log transformed to reduce data overdispersion before their use in the models. For SST, data corresponded to daily averages of level-4 satellite images derived from the Multi-Scale Ultra-High Resolution (MUR) SST Analysis database (Dataset ID: jplMURSST41). MUR-SST maps merge data from different satellites, combined with in-situ measurements, using the Multi-Resolution Variational Analysis statistical interpolation 59 , in a grid size of 0.01 × 0.01 degrees (ca. 1 km 2 ). From MUR-SST maps, thermal gradients maps were generated for each day that whale locations were available using the R package "grec" v. 1.3.0 60 with the Contextual Median Filter algorithm 61 as the method for calculating gradients. MUR-SST and thermal gradients maps were used to extract the associated covariate values for each whale location.
Vessel traffic data. To characterize vessel traffic patterns in the area, daily vessel tracking information (time-stamped GPS locations for individualized vessels) was obtained from the Chilean National Fisheries and Aquaculture Service (SERNAPESCA), available at www.serna pesca .cl. This database was released by the Chilean government during 2020 and comprises data involving the industrial and artisanal fisheries, aquaculture, and transport fleets, from March 2019 to present (updated daily). According to Chilean legislation it is mandatory for these fleets to provide tracking information to SERNAPESCA, except for artisanal fishing vessels smaller than 15 m and also for those smaller than 12 m in the case of artisanal purse seiners (www.bcn.cl). Artisanal fishing fleet comprises vessels up to 18 m in length and less than 80 cubic meters of storage capacity; above these metrics fishing vessels are considered part of the industrial fishing fleet. The transport fleet comprises vessels with no size limitations, engaged solely in the transportation of fishery resources. The aquaculture fleet is the most diverse www.nature.com/scientificreports/ one, considering its different operations (e.g. staff commuting, live and processed resource transportation, and supplies and infrastructure movement) with vessel sizes ranging from 5 to 100 m. All procedures described next were conducted independently for each fleet during data analyses. We used an 8 × 8 km grid to calculate vessel density (VD i ) for each grid-cell i. Vessel data are provided daily, with data gaps occurring for some days. Therefore, VD i was calculated by summing the daily number of unique vessels crossing each grid-cell i in a month divided by the total number of days with available data (range: 25-31 days). This procedure was conducted for austral summer and austral autumn months (March-June of 2019 and January-June of 2020) and then averaged into a single layer. Potential large differences in traffic patterns between these months were visually inspected through plots, which can be found as Supplementary Figures S1-S4 online. Data from austral winter and austral spring months were not used as most of the blue whale population is absent from the study area during these months 13,14 . Modeling approach. Telemetry data analysis has motivated the development and increasing use of various state-space modeling (SSM) approaches, which deal with path reconstruction and complex latent behavioral states 30,31,62,63 . Most practical applications of SSM, however, are computationally intensive and therefore require a long time for fitting them. Recently, SSM has been implemented via Template Model Builder (TMB), a R package that relies on the Laplace approximation combined with automatic differentiation to fast-fit models with latent variables [64][65][66] . Based on "TMB" tools, we fitted a continuous-time correlated-random-walk model (CTCRW) which estimates two state variables, velocity and true locations from error-prone observed locations, and two parameters, β controlling autocorrelation in directionality and velocity and σ controlling the overall variability in velocity 62 . Variances for modelling error in locations were derived from the Argos error ellipse 67 . As the error ellipses data were not available for tags deployed in 2004, we calculated the mean error ellipse for all location classes in the newer tags (2013-2019) and assigned these values to the corresponding location classes for tags deployed in 2004.
The original version of this model (with no behavioral variation) was fitted to obtain estimates of the true locations in whale's paths and used these to extract the corresponding covariate values from DAHCC, SST and thermal gradients rasters. The mean of the covariate values within a 3 km radius from each estimated location was used to partially account for uncertainty in covariate data arising from observation error. This error radius corresponded to twice the known error for Argos location classes 3, 2 and 1 67 . Covariate data were standardized, and missing values were filled with zeros, which correspond to the mean in standardized variables. This only affected 6 whales (ID#s 1,6,7,10,11 and 12), it was restricted to SST and thermal gradient data, and except for one whale never exceeded more than 2.7% of the data (with ID#7 at 10.4% of the data). We modified the original version of the CTCRW by allowing β t and σ t to be random variables that vary in time as a function of environmental covariates.
where B0 and A0 are intercepts, A and B are vectors of slopes, X t is the corresponding design matrix holding the standardized covariates, and ε 1 and ε 2 correspond to standard deviations. In every case, the estimated standard deviation ε 2 for β t was extremely small and presented exceptionally large standard errors; therefore, instead of trying to estimate this parameter, we fixed it at 0.01. In cases where no covariate presented a significant effect on β t this variable was reduced to a single parameter β, which was estimated. Estimated values of β larger than 4 produce persistence values lower than 0.05 h, indicating that at very short time differences velocity and location are poorly correlated with previous values. Therefore, in cases where model estimates for β were higher than 4 (ID#s 5 and 10) β was fixed at 4 indicating overall poorly autocorrelated movement patterns.
Our modelling approach allowed us to quantify the influence of environmental covariates on β t and σ t , with higher values of σ t indicating higher velocities and higher values of β t indicating lower directional persistence, which might be expressed as p t = 3/β t in units of time 62 . As no discrete behavioral states were explicitly included in our model, we defined behavioral states as post hoc categories based on p t and σ t values and their medians. The expected ARS state (slower and less persistent movement) was defined for locations jointly holding values of p t and σ t below their medians and the opposite was defined as transit state. The other two logical combinations (high p t with low σ t and low p t with high σ t ) were also provided and their interpretation is further discussed below. We also calculated ν t = √ π * σ t √ β t * 2 , which corresponds to long-term velocity 68 . This variable is a function of both σ t and β t (or p t ), and hence higher ν t can be obtained by either increasing σ t or reducing β t. As ν t is a function of both σ t and β t , we considered it as a proxy for the ARS-transit continuum, with higher values of ν t representing more transit-like behavior. Expected responses of ν t to covariate variation were inspected through prediction curves.
Finally, model results were used to generate spatial predictions for ν i for each grid-cell i using a 1 × 1 km grid. These predictions indicate the expected behavioral responses for whales traversing areas not necessarily visited during the tracking period. Predictive layers were generated for individual whales and averaged across individuals for depicting an overall pattern. www.nature.com/scientificreports/ Integrating movement and species distribution models. Results from a previous SDM were used for assessing spatial overlap between blue whale distribution and marine traffic. Briefly, this model consisted of a Bayesian binomial N-mixture model used to model blue whale groups counts in line-transect data (2009, 2012 and 2014), using distance sampling techniques and oceanographic covariate data 16 . Using an 8 × 8 km grid spatial predictions of blue whale density at each grid-cell i (N i ) were generated for eight years (2009-2016) and averaged into a single layer. To integrate outputs from movement models and SDM the relative probability of encountering a whale (RPEW) was calculated as follows RPEW i assumes that the probability of encountering whales increases with predicted density 39,69 . Here we consider behavior might also be part of this function as slow and less persistent movement (ARS) will result in more time spent (1/ν i ) allocated to each grid-cell i relative to all other grid cells n. As N i , had a spatial resolution of 8 × 8 km, we resampled the ν i grid to match the coarser grid resolution prior to any calculation, using the mean of aggregated grid-cells.
Defining spatial overlap with marine traffic. A quantitative measure of risk associated to vessel traffic can be considered as a monotonic function of the number of vessels and the probability of encountering a whale 39,70 . As described above, the relative amount of time allocated to each grid-cell can be obtained from 1/ν i . Therefore, as a measure of risk we calculated the relative probability of vessel encountering whale (RPVEW) 39,69 by combining N i , ν i and VD i as follows. where corresponds to the probability of observing a whale within each grid-cell i relative to all other grid cells n, corresponds to the time allocated to each grid-cell i relative to all other grid cells n, and Pv i = VD i n i=1 (VD i ) corresponds to the observed number of vessels within grid-cell i relative to all other grid cells n. fleets. Finally, to generate quantitative estimates on the degree of overlap between blue whale distribution and vessel traffic we used the Shoener's D and Warren's I similarity statistics 71 . These statistics range from 0, indicating no overlap, to 1, indicating distributions are identical. To use these statistics, the variables N i times 1/ν i and VD i were rescaled to range between 0 and 1 and inputted to the nicheOverlap function from the R package dismo 72,73 . A schematic representation of our workflow can be found as a Supplementary Figure S5 online.

Statement of approval.
The tagging methods employed in this study were approved by the Institutional Animal Care and Use Committee of the National Marine Mammal Laboratory of the Alaska Fisheries Science Center, National Marine Fisheries Service, U.S. National Oceanic and Atmospheric Administration. All methods employed in this study were carried out in accordance with guidelines from Subsecretaría de Pesca y Acuicultura (SUBPESCA), which provided full authorization to undertake this research through resolution #2267 of the Chilean Ministry of Economy and Tourism.

Results
Tracking duration for instrumented whales while within the study area ranged from 8.1 to 105 days (mean = 52.03, sd = 29.3, median = 48.7), yielding tracks that ranged from 49 to 1,728 locations (mean = 460.27, sd = 582.36, median = 140) used for modelling (after filtering, Table 1). In general, whales tended to remain in very localized coastal areas, where high productivity occurs during each austral spring (Fig. 2). No instrumented individuals departed from NCP until the onset of austral autumn-winter months (April-July) 14 . Pearson correlation analyses showed that none of the used covariates were strongly correlated (r < 0.5, p < 0.01). Except for one instrumented whale (ID#12), all animals showed a significant positive correlation between σ t and DAHCC, six animals showed a significant negative correlation between σ t and thermal gradients (Table 1). These results imply a clear pattern of whales reducing their velocities near areas that were highly productive during spring each year and/or where higher thermal gradients occur. The relationship with SST was less clear as three individuals showed a significant negative correlation and five a significant positive one (Table 1).
Regarding correlations between β t and environmental covariates, it was expected that whenever significant, they would present the opposite sign of those that were significant regarding σ t , rendering a continuum between ARS and transit behavior. This was the case for three individuals with respect to DAHCC (ID#s1, 4 and 8), four individuals with respect to SST (ID#s 1, 4, 8 and 11) and one individual with respect to thermal gradients (ID#8 ,  Table 1). Interestingly, two individuals showed the same signal in their correlation between DAHCC and β t , as well as, between DAHCC and σ t (ID# 11 and 15). The same occurred for one individual regarding SST (ID#9) and one individual regarding thermal gradients (ID#13, Table 1).
Post hoc definition of behavioral states showed the expected occurrence of both transit and ARS behavior. However, it also showed the occurrence of intermediate behavioral states at locations associated with low speed and high persistence and vice versa (Fig. 2). These types of intermediate behaviors were more predominant in individuals tagged in 2016 and 2019.  www.nature.com/scientificreports/ Prediction curves for ν t based on covariate variation provided unrealistic predictions for individuals for which a relatively small number of locations were available (< 200 locations, Fig. 3). For this reason, we only generated spatial predictions of ν t (Fig. 4) for individuals having tracks with more than 200 locations (ID#s 5, 7, 11, 12, 13, 14 and 15). Interindividual variation was observed regarding absolute values for ν t , indicating that some whales moved, in general, faster and in a more persistent manner (Fig. 4 b,c,e) than others, and also in terms of where their lowest values (ARS behavior) were expected. Despite this individual variation, some areas were consistently depicted as having the lowest values for ν t , which are highlighted when the spatial predictions for these seven whales were averaged into an overall mean (Fig. 4h). Spatial predictions on RPEW highlighted areas of aggregation for blue whales in NCP, mainly located in the western part of Chiloe Island, Ancud Gulf, Adventure Bay and northern Moraleda Channel (Fig. 5).

Discussion
Blue whale habitat selection and priority areas for conservation. Understanding the environmental drivers of blue whale habitat selection 16,17 is paramount for defining priority areas for its conservation and developing recommendations for marine spatial planning 11,74 . In pursuing this goal, our setting combined previous SDM fit to line-transect data with a movement model fit to telemetry data in a complementary manner. Telemetry data supports the spatial pertinence of previously defined areas for assessing blue whale abundance and distribution patterns through ship-borne surveys. Although, some whales performed brief excursions to adjacent offshore waters, they tended to remain within the NCP coastal areas during most of the tracking time, which in two cases extended for up to 3 months (Table 1). Potential caveats to this approach include tagging location bias (i.e. only performed in coastal waters, Fig. 1) and sampling size, which should be overcome through the ongoing tagging program.
Previous SDM 16 showed that spring productivity and, secondarily, thermal fronts were important covariates for predicting blue whale densities. Results here show that the same covariates selected by SDM are important for understanding blue whale's movement patterns. As with the aforementioned SDM, DAHCC was the most prevalent covariate retained in our models, which combined with thermal gradients, displayed an unequivocal pattern in their correlation with σ t . This is, whales tended to reduce their velocity near areas of high primary productivity that had occurred during austral spring and where strong thermal gradients take place (Table 1, Fig. 3a). As with many other large whale species, worldwide abundance and distribution patterns of blue whales have been linked to predictable highly and seasonally productive waters associated to high chlorophyll-a, among other proxies for enhanced productivity 19,20,24,[75][76][77][78] . Nevertheless, as blue whales feed almost exclusively on krill, temporal lags are expected to occur between seasonally high primary productivity, euphausiids early life-history stage processes (e.g. larval recruitment), the peak in adult euphausiid densities and the peak in whale abundance 17,20,78 . Refining our understanding of how temporal lags relate chlorophyll-a to euphausiid spatial patterns and then to blue whale distribution remains a pending task 79,80 , especially considering that euphausiid spatial ecology in the NCP is poorly understood 49,81 .
Although spring chlorophyll-a appears to be a suitable general proxy for blue whale prey availability in the NCP, whales are expected to respond in a much more complex manner to environmental heterogeneity. Previously, blue whale density in the NCP was found to be higher near areas of thermal front recurrence 16 . By using telemetry data, we were able to refine the assessment scale and test whether blue whales responded to daily changes in thermal gradients. Despite the relatively coarse resolution of Argos data, we were able to find evidence for behavioral response in six whales while traversing thermal gradients of less than 1 °C (Fig. 3c). This may even represent an underestimation given the reported response of blue whales to gradients as low as 0.03 °C 82 . Thus, our results provide additional support on the relevance of coarse to meso-scale thermal gradients when shaping marine predator distribution 16,23,82,83 . The underlying mechanism for this pattern, however, is not clear, as thermal fronts might be responsible for increasing prey availability by boosting local productivity and/or by aggregating prey patches [22][23][24][25]83,84 . Within the NCP, both processes are likely to be tightly coupled. The influence of fresh waters rich in silicic acid, among other nutrients, from high river discharges due to glacier melt and heavy rain, fertilize the photic zone by mixing with macronutrient-loaded oceanic deep water 46,49,81,85,86 . This large fresh water input in conjunction with higher irradiance reaching the surface during spring and summer, wind stress, tide and complex bottom topography promotes alternating processes of vertical and horizontal stratification/mixing of the water column, enhancing primary production as well as plankton aggregation [87][88][89] . In this context, areas selected by blue whales in the NCP might not just be of high biological productivity, but where frontal dynamics lead to highly concentrated prey patches.
SST presented an equivocal pattern regarding blue whale movement patterns, suggesting a preference for colder waters in four individuals and the opposite in four other individuals (Table 1, Fig. 3b). This might be a temporal issue if whales in some years/seasons found their prey in colder/warmer waters. For instance, Ancud Gulf tends to present higher temperatures during spring and summer than the Corcovado Gulf as the latter represents the main entrance path for sub-superficial oceanic colder waters into the Chiloe inner sea. Alternatively, Blue whales appear to respond to dynamic water-column processes by performing continuous behavioral changes without necessarily departing from relatively discrete areas (e.g. Ancud Gulf and Moraleda Channel, Fig. 2). For instance, whales ID#11, ID#13 and ID#15 presented a higher probability of reducing their velocity nearby areas of high productivity and strong thermal gradients, a higher probability of increasing persistence nearby areas of high productivity (for whales ID#11 and ID#15, Table 1), and all three spent from one to 3 months within specific micro-basins (Ancud Gulf and Moraleda Channel). This suggests that both transit-like and ARS behaviors co-occur spatially, temporarily oscillating with the suitability of foraging conditions.
Higher blue whale densities observed in the same areas where tagged individuals presented ARS behavior in a previous study 16 could have been attributed to multiple individuals entering and leaving these areas. However, the results presented here show that instrumented blue whales concentrate in relatively discrete areas for extended periods of time (up to 3 months) searching for and exploiting available resources. The limited movement elicited by blue whales might be regarded as an indicator of low interspecific competition, considering that their population abundance is still estimated to be considerably below pre-whaling levels 16,90,91 . Other mechanisms like dominance 92 and predator avoidance 93 , have been purported to explain limited animal movement. Thus, other factors should be considered in the future for understanding other dimensions of blue whales' habitat selection process, as well as temporal variations on it. Red to blue four-color ramp indicates the percentile to which each location belongs regarding variation in σ t and 3/β t (persistence). By using the medians, the four possible combinations are presented as a posteriori behavioral state identification. Locations jointly holding values of σ t and 3/β t below their medians across all whales (low s and low p) can be considered ARS behavior, while the opposite (high s and high p) can be considered transit. Blue (far) to yellow (close) color ramp in the background indicates variation in standardized distance to areas of high chlorophyll concentration (DAHCC) in log scale, which was the most consistent covariate shaping blue whale movement patterns in this study. Data layers (including maps) were created in R ver. 4.0.2 (www.r-proje ct.org) and ensembled in QGIS ver. 3.8.0 (www.qgis. org) for final rendering. Maps were created using data on bedrock topography from the National Centers for Environmental Information (https ://maps.ngdc.noaa.gov/viewe rs/grid-extra ct/index .html). Values above 0 were considered land coverage. www.nature.com/scientificreports/  www.nature.com/scientificreports/ Independently, both SDM and movement models predictions, highlighted similar areas of aggregation for blue whales in NCP based on observed oceanographic conditions (see Supplementary Fig. S5 online). These are clearly delimited by our RPEW map and considered Ancud Gulf, the Western coast of Chiloe Island, Corcovado Gulf / Moraleda Channel (CGMC), and Adventure Bay (Figs. 1 and 5). As previous SDMs were restricted to areas within 25 km from shore, some offshore areas visited by blue whales were not considered during RPEW computation. However, as the overall tendency to remain in coastal waters by instrumented whales was clear (Fig. 2f), we consider RPEW to be adequate.
Quantifying overlap with vessel traffic. For Chile, detailed and freely available vessel traffic data as those used here are limited to recent years (2019-2020), precluding long term assessments on vessel traffic spatiotemporal variation 95 . Although limited to 10 months of data, results showed little intra-fleet variation for the transport and aquaculture vessel activities, as well as, for those occurring in the inner sea for both fishing fleets (see Supplementary Figs. S1-S4 online). This was expected as transport and logistic support operations from aquaculture operations are less variable than the shifting resource-tracking operations of fishing vessels. In addition, the inner waters concentrate obligated marine corridors for entering/leaving the area which are used similarly regardless of vessel type. Henceforth, our estimates are expected to adequately reflect general vessel traffic patterns for each fleet but inspecting possible temporal variation in these patterns should be pursued in the future.
The four different vessel fleets considered here elicited differences in VD values and their spatial use of the study area (Fig. 6). While artisanal and industrial fishing fleets utilize inner waters to the east and open waters to the west of the study area, aquaculture and transport fleets are mainly constrained to inner waters (Fig. 6). According to Chilean legislation, the artisanal fishing fleet is restricted to operate within 5 nm (9.3 km) from the coast in open and inner waters while the industrial fishing operations are to be performed beyond this area to the West. This might explain the artisanal fishing fleet´s high score on the similarity statistics, indicating the largest degree of overlap with blue whale coastal distribution. In other words, this fleet distributes the RPVEW more www.nature.com/scientificreports/ homogeneously matching blue whale distribution, while other fleets concentrate only at specific areas (lower degree of overlap). In comparison with results presented here, a study using the same overlap statistics, showed a higher degree of overlap between vessels and three species of cetaceans in the Mediterranean Sea 73 . This was expected as the Mediterranean Sea is a high intensity vessel traffic area 96 . However, most of the marine traffic recorded in that study (73.3%) corresponded to small sailing boats, suggesting low probabilities of lethal shipstrikes in general but pinpointing that shipping routes (where larger vessels navigate) might pose higher risk. This brings forward the fact that spatial overlap is just one of the factors affecting collision risk and its outcome, with vessel density, speed and size also contributing to it 39,40,97 . Although the industrial fishing fleet presents a lower degree of spatial overlap with blue whales and the lowest number of operating vessels, industrial vessels might yield a higher probability of lethal interactions if they occur, due to larger vessel size. This fleet also presented a particular pattern of high RPVEW values off Adventure Bay (Fig. 6).
With up to 729 active vessels operating per day (83% of the total) and up to 78 vessels per day crossing a single grid-cell (VD), aquaculture fleet corresponds to the largest and most densely distributed fleet in the NCP. Hence, while RPVEW predictions highlights the specific areas where interactions are more likely to occur for each vessel fleet, in absolute terms, it is possible that the aquaculture fleet represents the major driver of negative vessel-whale interactions in NCP.
When considering results from all fleets together it is clear that the inner waters largely concentrate higher VD and high RPVEW values for all fleets (Fig. 6). This area holds the largest number of human settlements in the NCP and the main port pertaining to the regional capital, Puerto Montt, raising concerns for potential collisions, behavioral disturbance and/or heavy noise exposure 38,94,[98][99][100][101] for blue whales there. Although, no systematic monitoring or registering protocol exists in this region, local authorities' statements and the local press have documented at least three large whale mortality events linked to vessel collisions in the NCP (two Maps were created using data on bedrock topography from the National Centers for Environmental Information (https ://maps.ngdc.noaa.gov/viewe rs/grid-extra ct/index .html). Values above 0 were considered land coverage. www.nature.com/scientificreports/ blue whales and one sei whale), with two occurring nearby Puerto Montt and the other one at CGMC (Fig. 5). The ability of blue whales to avoid approaching vessels appears to be limited to relatively slow descents/ascents, with no horizontal movements away from a vessel 102,103 , therefore, collision events might pose significant threats to survival and recovery 97 for this endangered population. As inner waters of NCP might be considered, at the time, the spot of higher relative and absolute probabilities of negative interactions between blue whales and vessels, management actions are urgently needed to be implemented. For now, the most effective way to reduce collision risk is to keep whales and vessels apart, either in space or time, and where/when this is not possible, other measures (such as speed regulation) can be sought and applied singularly or in combination, considering variations in vessel activity and whale´s distribution 40,102,104 , as data become available. In addition, it is important to acknowledge that all analyses performed here were restricted to vessels carrying transponders and legally mandated to submit position data. Therefore, several vessels types operating in the area that could contribute to collision risk (e.g. international cargo and tankers, cruiseliners, as well as artisanal, recreational and military vessels) are currently unaccounted for. Because widely migratory species, such as the blue whale, do not recognize political boundaries, it is of great importance to identify the location of corridors and critical areas where they perform their vital activities (i.e., feed, migrate, breed, calve) to provide baseline information for their conservation. Efforts must be implemented at the local, national and international scales if success is to be reached, as ESP blue whale population recovery might be jeopardized by the loss of even a few individuals a year 16 after being severely depleted by the whaling industry during the 20th Century.
Modelling approach. One of the main differences between our modelling approach and previously published SSMs is in that behavioral variation that arises from the dependence on time-varying parameters (σ t and βt) rather than switches in discrete pre-determined behavioral states 30,31,65,107 . While the latter approach allows formal prediction, testing on the spatio-temporal occurrence of known behavioral modes (e.g. areas where ARS is likely to occur), time-varying approaches permit investigating variation in movement patterns that cannot, or are not desired to be, categorized a priori 65,107,108 . This poses a significant advantage in cases where animal movement fails to conform to the usual transit/ARS binary view. For instance, a previous work 14 fitted a switching SSM to most of the data we analyzed here and found that transit states were very rare within the NCP. In agreement with this, our results show that ca. 75% of all whale estimated locations presented persistence values lower than 1.6 h, which is consistent with the biological expectation of whales primarily engaged in foraging related activities within NCP 12 . In this scenario, attempting to explore the effect of environmental variables on switching probability between ARS and transit states 76 would have been difficult, as very few locations and their associated covariates would have been available for the transit state. By exploring changes in movement parameters, we can assess how animals' velocity and/or persistence respond to environmental covariates without the need of further assumptions. Following the transit/ARS rationale of conventional switching SSMs, one would expect that if a covariate is correlated with σ t it also would be with βt, but with an opposite sign. That is, at certain covariate values an animal's velocity and persistence are likely to decrease indicating ARS behavior, as was the case for several individuals and variables (Table 1). However, this does not need to always be the case, as shown by whales ID#11 and ID#15, which reduced their velocity near areas of high productivity in conjunction with increased persistence (Table 1). In general, this might occur because both transit and ARS behavior co-occur in similar areas with respect to DAHCC but differ in other variables (SST and thermal gradients). Nonetheless, alternative explanations for other behaviors, apart from transit/ARS, might arise. For instance, short-lived chasing bursts (escorting-like behavior) has been described for the NCP 109 , which are expected to present high velocities but not necessarily high persistence. On the other hand, slow persistent behavior, mostly present in whales tagged in years with the highest data transmission throughput (2016-2019, Fig. 2d-e, Table 1), might be explained by the ratio of the location error relative to the scale of movement. Thus, if short time periods separate two or more locations with limited movement, high persistence might arise from negligible variation in both speed and location, as observation error increases disproportionately relative to the scale of the movement process.
Overall, our modelling approach accounted for observational error and allowed for the incorporation of environmental covariates to inform movement parameters without the need for regularization of location data into fixed time intervals 30,65 , all in one single step. By fitting the model through the R package "TMB" analysis took an average of 60.5 s to run (range: 2.6-310.6, processor: Intel Core i7-7700HQ at 2.8 GHz, RAM: 32 GB) which is a significant advantage when processing large amounts of data.

Conclusions
Blue whale movement patterns agree with previous studies on their distribution, highlighting the importance of coastal waters and reinforcing our knowledge about primary production and thermal fronts as important environmental drivers for this species´ habitat selection process in the NCP. Considering defined priority areas for blue whale conservation in the area, those located at inner waters concentrated the highest probabilities of whales interacting with vessels. Among the studied vessel fleets, the unparalleled size of the aquaculture fleet indicates this could play a decisive role in modulating potential negative vessel-whale interactions within NCP. The results of this study clearly pinpoint specific areas where management actions are urgently needed, especially considering the undetermined number of vessels strikes and levels of noise exposure in the region. This information should be considered by Governmental and International organizations to inform, design, and rapidly implement mitigation action using existing national and international conservation instruments.