The Role of Climate, Oceanography, and Prey in Driving Decadal Spatio-Temporal Patterns of a Highly Mobile Top Predator

Marine mammals have been proposed as ecosystem sentinels due to their conspicuous nature, wide ranging distribution, and capacity to respond to changes in ecosystem structure and functioning. In southern European Atlantic waters, their response to climate variability has been little explored, partly because of the inherent difficulty of investigating higher trophic levels and long lifespan animals. Here, we analyzed spatio-temporal patterns from 1994 to 2018 of one of the most abundant cetaceans in the area, the common dolphin (Delphinus delphis), in order to (1) explore changes in its abundance and distribution, and (2) identify the underlying drivers. For that, we estimated the density of the species and the center of gravity of its distribution in the Bay of Biscay (BoB) and tested the effect of three sets of potential drivers (climate indices, oceanographic conditions, and prey biomasses) with a Vector Autoregressive Spatio Temporal (VAST) model that accounts for changes in sampling effort resulting from the combination of multiple datasets. Our results showed that the common dolphin significantly increased in abundance in the BoB during the study period. These changes were best explained by climate indices such as the North Atlantic Oscillation (NAO) and by prey species biomass. Oceanographic variables such as chlorophyll a concentration and temperature were less useful or not related. In addition, we found high variability in the geographic center of gravity of the species within the study region, with shifts between the inner (southeast) and the outer (northwest) part of the BoB, although the majority of this variability could not be attributed to the drivers considered in the study. Overall, these findings indicate that considering temperature alone for projecting spatio-temporal patterns of highly mobile predators is insufficient in this region and suggest important influences from prey and climate indices that integrate multiple ecological influences. Further integration of existing observational datasets to understand the causes of past shifts will be important for making accurate projections into the future.

Marine mammals have been proposed as ecosystem sentinels due to their conspicuous nature, wide ranging distribution, and capacity to respond to changes in ecosystem structure and functioning. In southern European Atlantic waters, their response to climate variability has been little explored, partly because of the inherent difficulty of investigating higher trophic levels and long lifespan animals. Here, we analyzed spatiotemporal patterns from 1994 to 2018 of one of the most abundant cetaceans in the area, the common dolphin (Delphinus delphis), in order to (1) explore changes in its abundance and distribution, and (2) identify the underlying drivers. For that, we estimated the density of the species and the center of gravity of its distribution in the Bay of Biscay (BoB) and tested the effect of three sets of potential drivers (climate indices, oceanographic conditions, and prey biomasses) with a Vector Autoregressive Spatio Temporal (VAST) model that accounts for changes in sampling effort resulting from the combination of multiple datasets. Our results showed that the common dolphin significantly increased in abundance in the BoB during the study period. These changes were best explained by climate indices such as the North Atlantic Oscillation (NAO) and by prey species biomass. Oceanographic variables such as chlorophyll a concentration and temperature were less useful or not related. In addition, we found high variability in the geographic center of gravity of the species within the study region, with shifts between the inner (southeast) and the outer (northwest) part of the BoB, although the majority of this variability could not be attributed to the drivers considered in the study.

INTRODUCTION
The global mean surface temperature has increased by approximately 1 • C from pre-industrial levels (IPCC, 2019), triggering shifts in the abundance, phenology and distribution of organisms worldwide (Parmesan and Yohe, 2003;Poloczanska et al., 2013). Marine ecosystems, despite having experienced a slower warming, show comparable or even greater shift rates and vulnerability than terrestrial systems (Burrows et al., 2011;Poloczanska et al., 2013;Pinsky et al., 2019), with seagrasses, corals, cephalopods and marine mammals exhibiting the most abrupt changes (Trisos et al., 2020).
Marine mammals, as wide ranging top predators, amplify trophic information across multiple spatiotemporal scales and can therefore act as sentinels of ecosystems' responses to climate variability and change . However, assessing climate change impacts in higher trophic levels and long lifespan animals such as marine mammals is challenging, as their relationships to climate may be non-linear and affected by time lags (Simmonds and Isaac, 2007;Barlow et al., 2021). In addition, identifying spatio-temporal trends in the context of climate change requires analyzing decadal or longer time series (Thorson et al., 2016), which are rarely available for marine mammal observation data.
Combining data from multiple sampling programs can help overcome this problem Maureaud et al., 2021), but also increases the intrinsic variability related to observers' skills, sampling design and protocols, which may result in confounding species range shifts with variations in the distribution and intensity of the sampling effort (Thorson et al., 2016). For that reason, separating the observation process from the true underlying spatial distribution is essential to accurately identify range shifts over time (Chust et al., 2014b) and to identify potential drivers (Erauskin-Extramiana et al., 2019b). Recently, a species distribution function (SDF) able to distinguish between sampling variation and true geographic variability has been developed (Thorson et al., 2016). Unlike conventional estimators such as the abundance-weighted average, the SDF is applied through a Vector Autoregressive Spatio Temporal (VAST) model that allows the estimation of species distribution over predicted locations rather than sampled locations, while also estimating a standard error that allows one to distinguish between sampling variation and significant variability (Thorson et al., 2016). Although model-based approaches had been used before to estimate shifts in the distribution of species, VAST typically involves estimating a Gaussian Markov random field (GMRF) representing latent variation in density that is constant over time (a "spatial" term), as well as a GMRF representing latent variation that changes among years (a "spatio-temporal" term), which is expected to improve predictions of species density and distribution compared with using only measured habitat variables (Thorson, 2019).
Until now, this estimator has been mainly applied to commercially important fish stocks (Godefroid et al., 2019;Perretti and Thorson, 2019;Xu et al., 2019), although the fragmented and methodologically variable nature of marine mammal observations suggest the method could be highly useful for analyzing the spatio-temporal patterns of marine megafauna too. Within that context, the Bay of Biscay (BoB hereafter), located in the Northeast Atlantic, off the coasts of France and Spain (Figure 1), represents an interesting study area since numerous marine mammal species (e.g., cetaceans) cohabit there, attracted by a highly diverse and abundant community of pelagic fish species (Astarloa et al., 2019;Louzao et al., 2019).
Such productivity and diversity, however, might be altered by climate change in the near future, as rising temperatures (0.26 • C per decade; Costoya et al., 2015) are expected to increase ocean stratification and reduce primary production and zooplankton biomass in the area (Chust et al., 2014a). In recent years, losses in fisheries production have already been reported (Free et al., 2019), together with changes in the composition, distribution, and phenology of fish species (Blanchard and Vandermeirsch, 2005;Chust et al., 2019;Baudron et al., 2020). Cetacean spatiotemporal variability, in contrast, has been mainly assessed by exploring changes in their relative abundance (Hemery et al., 2007;Castège et al., 2013;Authier et al., 2018), although both abundance and distribution are considered key criteria by the European Marine Strategy Framework Directive (MSFD; Directive 2008/56/EC) aiming to assess the environmental status of species and ecosystems in European Union waters.
Advancement of both MSFD criteria in this region is therefore necessary, especially when it is known that projections of climate change impacts on cetaceans at large spatial scales (e.g., global; MacLeod, 2009) do not always match with those at regional scales (Hazen et al., 2012). In the Northeast Atlantic, for example, warm-water cetaceans were predicted to expand poleward (MacLeod, 2009;Lambert et al., 2011Lambert et al., , 2014, although the southeastward shift detected for some Northeast Atlantic fish species in the BoB could indicate the opposite pattern in this particular area (Baudron et al., 2020). Indeed, some of the fish species (e.g., horse mackerel Trachurus trachurus, anchovy Engraulis encrasicolus, and sprat Sprattus sprattus) analyzed by Baudron et al. (2020) constitute an important food resource for many cetaceans in the BoB (Meynier et al., 2008;Spitz et al., 2018), FIGURE 1 | Spatial distribution of common dolphin sightings (displayed in segments of up to 10 km) over the BoB for the 1994-2018 period. Circle sizes are proportional to group size, while solid gray lines indicate the isobaths. Sightings in yellow represent the ferry data used to check model fit.
which can heavily influence the spatial movements of their predators Giralt Paradell et al., 2019).
The hypothesis that climate change may affect top predators through climate influences on their ectothermic prey has been often suggested (Robinson et al., 2005;Simmonds and Isaac, 2007;Evans and Waggitt, 2020). Most studies, however, examine environmental conditions (e.g., temperature) as proxies of prey distribution rather than studying prey data directly (Torres et al., 2008;Giralt Paradell et al., 2019) while others focus on exploring the effects of climate indices on the grounds that they act as an integrated measure of multiple variables (Hallett et al., 2004;Hemery et al., 2007). In the Northeast Atlantic, the North Atlantic Oscillation (NAO) is the dominant mode of climate variability, although additional climate indices such as the Atlantic Multidecadal Oscillation (AMO), the East Atlantic pattern (EA), or the South Biscay Climate (SBC) have been also found to exert strong influence, direct or indirectly, on both fish and cetacean species (Guisande et al., 2004;Hemery et al., 2007;Borja et al., 2008;Evans and Waggitt, 2020) through changes in ocean temperature and salinity, vertical mixing and circulation patterns (Drinkwater et al., 2003;Hurrell and Deser, 2009).
Given the multiple drivers potentially influencing cetacean spatio-temporal patterns, understanding the role of each of them is key for a better anticipating of future responses. For that reason, in this study we used a 25-year-long temporal series  to test the effect of prey biomasses, oceanographic conditions and climate indices on the abundance and distribution of the common dolphin (Delphinus delphis), one of the most abundant cetaceans inhabiting the BoB waters (Hammond et al., 2017). We used the VAST model (Thorson and Barnett, 2017) and the spatio-temporal species data compiled by Waggitt et al. (2020) to address two main research questions: (1) Has the abundance or the distribution of the common dolphin in the BoB experienced significant changes over the last two decades? (2) If so, are changes best explained by climatic, oceanographic, or prey variables? By answering these questions, this study intends to provide insights that will help understand past and future trends in the distribution and abundance of common dolphin in the BoB while contributing to the management for this species through the development of MSFD criteria in the context of climate change.

Data Collection and Standardization
Cetacean data analyzed in this study, despite focusing on the BoB, belong to a large compilation made by Waggitt et al. (2020) that included observations collected on aerial and vessel (dedicated and opportunistic) surveys conducted in the Northeast Atlantic between 1980 and 2018. Although the data analyzed here (data providers in Supplementary Table 1) is a more updated version that includes higher-resolution tracklines (meaning that fewer data were omitted due to overlap with land-masses and more accurate measurements of distance traveled were obtained), the steps taken in the data processing and standardization stage were the same as in Waggitt et al. (2020), in which they (1) assessed differences in protocols by grouping data according to the (a) survey transect design (line transects, strip transects, and an intermediate method called ESAS, European Seabirds At Sea) and (b) the platform-type (vessel vs. aircraft) and (2) fitted detection functions using platform height and Beaufort sea-state as explanatory variables to estimate the proportion of animals missed by the observers (Marques and Buckland, 2004). They also assessed response bias (when animals react to the presence of the platform) through double-platform surveys that enabled the detection of animals before responsive movements. This correction was applicable to vessel surveys and is particularly relevant to common dolphins, which typically show a positive response to vessels (Cañadas et al., 2004). Finally, they calculated the effective strip half-width (ESW) which considers the decline in the detection probability as a function of distance and covariates and serves to estimate the area effectively covered (Area covered = ESW * s * L) when including the number of observation sides (s) and transect length (L). Full details can be found in Waggitt et al. (2020).

Sampling Effort
In order to match with the spatial resolution of the environmental data that we examined in later steps (see "Identification of Main Drivers" section), we divided larger transects into 10 km segments . Then, we examined the spatiotemporal coverage of surveys by summing the effort comprised in all segments per month and per year. In addition, we checked whether compiling data had led to a non-uniform distribution of sampling in space and time by exploring the annual latitudinal and longitudinal mean distributions and the corresponding linear regression trends.

Baseline Spatio-Temporal Model
Observations of common dolphin were analyzed by means of a spatio-temporal delta-generalized linear mixed model (delta-GLMM), referred to here as a VAST model (Thorson and Barnett, 2017) and available in R. 1 This model is a flexible variant of the classical delta models that decompose density into two components (Stefánsson, 1996): (1) the probability of encountering the species at a given location and time; and (2) the expected density of the species when encountered. This two-part approach, also known as a hurdle model, helps combat statistical problems with zero-inflation and overdispersion in the original data (Martin et al., 2005) and is therefore suitable for use with cetacean survey data that usually show patchy distributions .
Another feature of the VAST model is that it decomposes spatio-temporal patterns in available point-count data into multiple additive components: 1. A temporal main effect ("intercepts") representing changes in median abundance over time; 2. A spatial main effect ("spatial component") representing the average spatial distribution during the modeled interval; 3. An interaction of space and time ("spatio-temporal component") representing variation in distribution among years; 4. Density covariates, representing the impact of environmental conditions on expected density; 5. Catchability (a.k.a. detectability) covariates, representing the impact of environmental and/or sampling conditions on expected sampling data, but which do not reflect variation in population density and hence are "partialled out" prior to predicting densities. 1 https://github.com/james-thorson/VAST Each of these components can be included in each of two linear predictors, and these two linear predictors are then transformed via inverse-link functions to predict the value of a response variable (in this case, dolphin samples). Spatial and spatio-temporal components are estimated as a Gaussian Markov random field (GMRF) and treated as a random effect. To improve computational speed, the value of these GMRFs is predicted at a fixed set of "knots" that defines a mesh of triangles that covers the entire modeled spatial domain. The value of the GMRF at any location within this domain is then predicted from the value of three knots surrounding that location. We use the stochastic partial different equation (SPDE) approximation to calculate the probability of GMRFs (Lindgren et al., 2011), and the projection from knots to locations is accomplished using bilinear interpolation as computed using R-INLA (Lindgren, 2012). The value of fixed effects are estimated using maximum likelihood techniques while integrating across the probability of random effects (Kristensen et al., 2016), and standard errors are calculated using a generalization of the delta method (Tierney et al., 1989). For further details, please see the VAST user manual. 2 In our case, we treated year as a fixed effect (default VAST setting), such that there is no shrinkage in overall abundance across years. We modeled spatial and spatio-temporal variation as random effects to help account for multidimensional factors that are not included directly in the model but that can affect the density and distribution of the modeled species . In particular, we estimated first-order autocorrelation among years in the spatio-temporal component, such that predicted hotspots in density decay slowly over time; this treatment allows spatio-temporal patterns to be predicted (with associated uncertainty) even in locations with sporadic sampling.
Detectability covariates were not considered here, because Beaufort sea-state and platform height were included in Waggitt et al. (2020). Density covariates were also omitted for our initial investigation of trends (but see "Identification of Main Drivers" section). As a response variable, the density of common dolphin was analyzed, after truncating the highest 5% to control outliers (Buckland et al., 2001). The spatio-temporal model was fitted assuming a lognormal error distribution and a Poisson-linked delta model such that the sum of both linear predictors is predicted log-density; this structure, was selected based on the lowest Akaike Information Criterion (AIC) (Sakamoto et al., 1986). Model parameters, as well as spatio-temporal components, were estimated using 200 knots (Supplementary Figure 1) based on previous studies that applied this same resolution in bigger areas Thorson, 2019), while confirming that results are qualitatively similar when increasing the number of knots (Supplementary Table 2). Species density was predicted at each knot by multiplying the predicted probability of occurrence by the predicted density. Density estimates for each knot were then interpolated to a standard grid of 0.1 • spatial resolution (latitudinal range: 43 • -49 • N; longitudinal range: 1 • -10 • W) to match with the spatial resolution of the environmental data (see "Identification of Main Drivers" section) and multiplied by the area of the grid cell to create annual surfaces of common dolphin abundances across the BoB.
The annual abundances of common dolphin predicted for the study area were then analyzed by means of a linear regression to identify significant temporal trends and compared by means of a correlation test with an observed abundance index to check model fit. The observed abundance index was based on the encounter rate (individuals/km) of common dolphin estimated from monthly at-sea observations taken by a team of experienced observers in a constant effort-based systematic sampling scheme, i.e., the Pride of Bilbao ferry (Louzao et al., 2015;Robbins et al., 2020). This survey consistently crosses the BoB using the same route every year (Figure 1), and although it was also used as input for the baseline model, it only forms the 8% of the whole data set. Thus, we believe it can be used to compare the observed (ferry) and predicted (VAST) abundance indices and to determine whether the model predictions have been biased by differences in the effort.
An additional analysis with predicted abundances was also conducted to identify areas in which significant spatio-temporal changes occurred over the study period. For that, predicted abundances per grid cell were analyzed as a function of year by means of a linear regression. The slope and the p-value obtained in each cell, as indicators of change rate and its significance, were then plotted over the standard grid covering the study area.

Distribution Shift Metrics
Shifts in distribution were summarized by calculating the centroid of the distribution for a given year (termed center of gravity, CoG) after having predicted the density associated with every knot and year in the previous step. By means of the SDF estimator implemented in the VAST model, the CoG was calculated for the BoB population domain and standardized by the total abundance predicted for the study area, so that our analysis focused on changes in distribution after controlling for changes in total abundance (Thorson et al., 2016). Shifts in CoG were displayed in terms of "Eastings" and "Northings, " meaning km from the most western point of the study area and km from the Equator, respectively. Significant trends were identified using a linear regression against year.

Identification of Main Drivers
To understand spatio-temporal patterns, three main groups of drivers were analyzed (Table 1), classified into local and regional covariates according to their spatio-temporal structure (a local covariate varies across space while a regional covariate is a univariate time series representing the covariate over the entire study area; Thorson, 2019): (1) Local oceanographic conditions integrated at 100 m depth, specifically temperature and chlorophyll a concentration (Chl-a), based on their direct relationship with climate change and their importance for predicting top predators distribution (Hazen et al., 2012;García-Barón et al., 2020).
(2) Regional climate indices, specifically North Atlantic Oscillation (NAO), East Atlantic Pattern (EA), and Atlantic Multidecadal Oscillation (AMO) climate indices (details in Table 1), due to their ability to extract the leading pattern in weather and climate variability over the North Atlantic and their relationship to cetacean and prey populations (Simmonds and Isaac, 2007;Borja et al., 2008;Evans et al., 2010;Evans and Waggitt, 2020). (3) Regional biomasses of potential prey species, based on the assumption that climate change will affect cetaceans distribution through changes in their prey (Robinson et al., 2005;Simmonds and Isaac, 2007;Evans and Waggitt, 2020).
Temperature and Chl-a values were sourced from the Iberian Biscay Irish Ocean Reanalysis Model available at the Marine Environmental Monitoring Systems, 3 providing values at a 0.08 • spatial resolution, a 1-month temporal resolution and at 22 discrete depth intervals ranging from surface to 100 m depth. To test their effect on the annual estimates predicted by the baseline spatio-temporal model, the annual mean of both temperature and Chl-a was estimated integrating the data available in the first 100 m of the water column and then resampled with the raster package (Hijmans et al., 2017) at 0.1 • (∼10 km) resolution . The three climate indices were downloaded from the National Oceanic and Atmospheric Administration (NOAA) at a monthly scale and averaged to obtain annual values, 4 while the biomass of prey species was acquired from the International Council for The Exploration of Seas (ICES) website at annual scale. 5 We selected prey species based on their relative importance in the common dolphin's diet in the BoB (Meynier et al., 2008;Santos et al., 2013) as well as data availability and suitability because not every potential prey species (e.g., sprat, myctophids) was available for the spatio-temporal scale defined in this study. European anchovy (Engraulis encrasicolus) was the only prey species whose biomass had been estimated exclusively for the BoB. Horse mackerel (Trachurus trachurus) estimates were for the Northeast Atlantic, Atlantic mackerel (Scomber scombrus) and blue whiting (Micromesistius poutassou) for the Northeast Atlantic and adjacent waters and sardine (Sardina pilchardus) estimates for the Cantabrian-Atlantic Iberian waters (for information on the extent of stocks see Table 1). Although there is an assessment for the sardine stock of the BoB, data were only available from 2000 onward (ICES, 2019c), so we decided to use the biomass estimations from the Cantabrian sea and Atlantic Iberian waters instead after having checked that both indices were highly correlated (r = 0.87) and followed similar trends (Supplementary Figure 2). Finally, the biomasses of all species were summed and used as a proxy for total prey biomass available in the BoB.
For modeling purposes, local temperature and Chl-a variables were included as quadratic forms in the model to allow for non-linear responses (Perretti and Thorson, 2019). Regional climate indices were included as "spatially varying coefficients" as in Thorson (2019), which means that instead of estimating a single slope parameter presenting the effect of an oceanographic index on density, the model estimates a separate slope parameter for every modeled location (every knot). The biomass of each 1 | Summary of the local oceanographic, regional climatic and regional prey variables used in this study accompanied by a little description and the source from which they were obtained. prey species, as well as the total biomass index, were first log transformed and then included as spatially varying coefficients since they were also available as a single regional time-series. As a preliminary analysis, potential drivers were correlated with the abundance and CoG of common dolphin obtained in the previous baseline spatio-temporal model. Then, covariatesbased modeling was performed in two different ways to identify the most parsimonious drivers and to uncover the relative contribution of covariates: (1) Univariate spatio-temporal models were fitted for each variable using the same configuration as in the baseline spatio-temporal model. Univariate models were then compared with the baseline model by means of the AIC (Sakamoto et al., 1986). Only a decrease in the AIC > 2 in relation to the baseline spatio-temporal model was considered an improvement. When models differed by less than 2 units of AIC ( AIC ≤ 2), they were considered statistically equivalent (Arnold, 2010). The way in which covariates were related to the spatio-temporal patterns of common dolphin was also explored by plotting the functional relationships from the model parameters.
(2) Univariate models were fitted for each variable after setting the spatio-temporal variation (i.e., spatio-temporal random effects) to 0. This was done to remove the contribution of random effects and isolate the effect of the covariates since in VAST, random fields can also account for changes in distribution over time by capturing the residual spatial patterns that cannot be attributed to the fixed effect . The abundances and CoG obtained from these models were then compared with those from the baseline spatio-temporal model to determine the amount of variation attributable to covariates.

Sampling Effort
A total of 1728 sightings of common dolphin collected in 21 different surveys were analyzed (Figure 1 and Supplementary Table 1). Those surveys mainly covered spring-summer months and showed a peak of maximum effort between the 2007 and 2012 period (Supplementary Figure 3). The mean latitude of sampling also varied and shifted significantly south over time (p = 0.001), while no significant change was observed in the mean longitude of sampling (Figure 2).

Common Dolphin
The common dolphin abundance estimated by the baseline spatio-temporal model showed a significant increase (p < 0.001) throughout the study period, accompanied by high variability (Figure 3 and Supplementary Table 3). This increase was most pronounced over the more recent years (2011-2017) and mainly occurred in the southeast corner of the BoB (Figure 4). These results agreed with the ferry data, which also showed an increasing trend and a significant correlation (r = 0.7, p = 0.003) with the predicted abundances ( Supplementary  Figures 4, 5).
The CoG also showed a high interannual variability, but no significant trend was found over time in either of the two axes (Figures 5A,B). In contrast, the correlation between eastings and northings showed as significant pattern (p = 0.005) in the direction of the shift, indicating that the distribution of common dolphins generally varied between the inner (southeast) and the outer part (northwest) of the BoB (Figure 5C).

Drivers and Covariate Contributions
Neither the annual temperature nor the Chl-a concentration integrated at 100 m depth revealed a significant (p > 0.05) temporal trend across the full BoB (Supplementary Figure 6). The climate index AMO has remained in a positive phase since 1997, whereas NAO and EA indices have shown a higher variability with alternation between positive and negative phases (Supplementary Figure 7). Both anchovy and mackerel biomasses showed a significant (p ≤ 0.05) recovery after a period of low abundance, while sardine and horse mackerel underwent a severe decline (p ≤ 0.001). In contrast, blue whiting did not show any significant temporal trend (p = 0.2). The prey biomass index, on the other hand, exhibited a significant increase (p = 0.003), despite the large variability (Supplementary Figure 8).
The correlation between the potential drivers and the CoG (easting and northings) of common dolphin only showed weak relationships. In contrast, predicted abundance revealed several strong relationships (r > 0.5) with prey species, specifically mackerel and anchovy (positive correlation), and sardine and horse mackerel (negative correlation) (Figure 6). After prey species, only EA and NAO climate indices showed a moderate correlation with abundance (r∼0.40). Blue whiting was not significant (p > 0.05), while temperature, Chl-a, AMO and the prey biomass index showed weak relationships (r∼0.20) (Figure 6).
For covariates-based models, the AIC score showed that the most substantial decrease was for the NAO index and regional prey species biomasses (especially anchovy and sardine). Local Chl-a concentration, as well as horse mackerel and mackerel, only contributed slightly, while remaining drivers (temperature, AMO, EA, blue whiting and prey species biomass index) were not relevant in terms of AIC ( Table 2). Functional relationships of those important drivers revealed positive responses for NAO, anchovy, mackerel and negative for Chl-a, horse mackerel and sardine (Supplementary Figure 9).  Similarly, covariate-only models (with no random effects) showed that the NAO index and prey species biomasses were able to explain the increase in regionwide abundance of common dolphin (Figure 7). Chl-a concentration, despite having shown a decrease in AIC score ( Table 2), did not contribute to explain the variability in the relative abundance (Figure 7), and neither did temperature, AMO index, or blue whiting (Supplementary Figure 10). EA and biomass indices did show a higher contribution in terms of variability, but they were not identified as important drivers according to AIC score (Supplementary Figure 10).
In the case of CoG, only Chl-a and temperature contributed to explain the observed variability but, even then, only in a very small proportion (Figure 8). In fact, the variation in the CoG explained by these variables only accounted for about 10-20 km, while the spatio-temporal model suggested variation of 100-300 km.

DISCUSSION
The evaluation of the spatio-temporal patterns of common dolphin in the BoB agrees with the MSFD aiming to assess the abundance and distribution of species in European waters. Surveys providing information on species distribution and abundance in this region, however, have shown significant shifts in the spatial distribution of observations, which make necessary the application of methods such as VAST to account for uneven sampling effort.

Spatio-Temporal Trends in Common Dolphin Abundance
The modeling of common dolphin sightings revealed a significant increase in abundance, which is in agreement with previous studies conducted in the BoB (Hemery et al., 2007;Authier et al., 2018;Saavedra et al., 2018) and in the wider Northeast Atlantic (Hammond et al., 2017;Evans and Waggitt, 2020) that also reported an increasing trend. In addition, data from ferry surveys, known to perform the same route every year, showed the same pattern and confirmed that the results were not biased by the detected latitudinal shift in effort.
In addition, the predicted abundance estimates were found to be quite coherent with those obtained in previous surveys conducted in summer 2012 in the BoB (Laran et al., 2017) and in summer 2016 in the Northeast Atlantic (ICES, 2020), in which 490,000 (95% CI: 340,000-720,000) small delphinids (common and striped dolphins) and 634,000 (95% CI: 353,000-1,140,000) common dolphins were estimated, respectively. Although it is not possible to make a direct comparison with our predictions, the ratios for common/striped dolphins and Northeast Atlantic/BoB estimated from Hammond et al. (2017) would lead to an approximate abundance of 360,000 (95% CI: 250,000-526,000) and 425,000 (95% CI: 237,000-764,000) individuals of common dolphin in the BoB for 2012-2016, respectively. These numbers were similar to our predictions in those years (359,000 ± 49,000 and 376,000 ± 71,500 individuals, respectively; Supplementary Table 2), and would indicate that, overall, abundance estimates from VAST were consistent with previous studies. This good agreement is remarkable, given the heterogeneity of the data used in this study that comprised 21 datasets, and emphasizes the importance of applying methods that are robust to shifts in sampling effort. In addition, the concordance between our results and those estimates made on summer also suggest that the spatiotemporal patterns obtained in this study should be interpreted as spring-summer trends, as this was the period of the year when most data were collected (Supplementary Figure 2B).
The increasing trend in abundance found in this study for the BoB, however, does not necessarily imply an overall population increase at the Northeast Atlantic level (i.e., species whole distribution range), and instead, could be due to the arrival of individuals from unsampled areas. That is why the results found in this study should be treated with caution and never be used to downplay the effects of incidental capture on common dolphin, especially when recent estimates suggest that the bycatch in the BoB is unsustainable for the population as a whole (ICES, 2020).

Regional vs. Locally Estimated Environmental Variables
Local environmental variables, such as temperature and Chla used in this study, are often unable to capture complex associations between environment and ecological process due to time lags in species responses coupled with the non-linear intrinsic nature of population dynamics (Hallett et al., 2004).
This can be particularly true for Chl-a and cetaceans species that feed on zooplanktivorous fishes, since the abundance of the latter has been related to a period of zooplankton grazing and FIGURE 6 | Pearson correlation among the common dolphin's predicted abundance, CoG and potential drivers. Circle sizes are proportional to the correlation coefficient, which is indicated inside the circles. Non-significant correlations (p > 0.05) are shown without a circle. a phytoplankton decay . Under such circumstances, many researchers working with cetaceans often apply time-lagged Chl-a concentration for one and/or 2 months prior to the sighting month (Tobeña et al., 2016;Prieto et al., 2017;Pérez-Jorge et al., 2020;Barlow et al., 2021).
In this study, however, predictors were introduced at an annual scale to match the available temporal scales of both prey and climatic indices, which prevented its incorporation in a lagged phase and likely led to the low contribution of Chl-a in explaining the spatio-temporal patterns of common dolphin. Similarly, the lack of importance shown by temperature could be also a consequence of this annual resolution or could instead suggest that, within the core of the species range, temperature is not such an important variable to explain its abundance and distribution.
On the contrary, regional indices of climate, spanning several months and considering wider areas of influence, are less disturbed by local variability and very often outperform locally estimated environmental variables (Hallett et al., 2004). In addition, they usually hold information about several environmental factors (e.g., temperature, storms and precipitation, mixed layer depths or circulation patterns), which make them act as an integrated measure of meteo-oceanographic conditions that tend to explain more of the variability of the system than just, for example, ocean temperature (Hurrell and Deser, 2009;Thorson, 2019).
The results found in this study are a good example of this, as the NAO climate index was found to be the best predictor explaining the abundance of common dolphin according to AIC scores. Specifically, results showed a positive relationship between both, meaning that common dolphin abundance is enhanced during positive phases of NAO, which are characterized by colder and drier conditions over Mediterranean regions, central and southern Europe (e.g., BoB), and warmer and wetter conditions in northern Europe (Visbeck et al., 2001;Aravena et al., 2009;Hurrell and Deser, 2009).
Although the NAO index and similar climate indices have been previously related to the abundance of wide ranging predators in the BoB (Hemery et al., 2007;Louzao et al., 2015), responses are likely mediated through the influence of the climate indices on food resources rather than directly on higher trophic predators such as cetaceans (Drinkwater et al., 2003;Lusseau et al., 2004). Indeed, the NAO climatic index has been related to some biologically important phenomena, such as upwelling (Pérez et al., 2010), river run-off (Dupuis et al., 2006) and Ekman transport (Guisande et al., 2004), which are known to influence the recruitment of some of the main prey species (i.e., anchovy, sardine) of common dolphin (Guisande et al., 2004;Borja et al., 2008;Planque and Buffaz, 2008). We could therefore hypothesize a potential bottom-up process, in which NAO affects common dolphins through its influence on prey. In fact, bottom-up control has been suggested for the continental shelf food web of the BoB, where a highly diverse and abundant community of forage fishes regulates higher trophic levels (Lassalle et al., 2011).

The Role of Prey
Common dolphins are assumed to be opportunistic predators that feed on a wide variety of species, although a preference for energy-rich species, such as the anchovy, sardine, mackerel and horse mackerel investigated in this study, has been suggested (Meynier et al., 2008). Atlantic mackerel, however, is only present in large quantities during the first half of the year in the BoB, coinciding with its spawning period (Uriarte and Lucio, 2001), while Atlantic horse mackerel and the Iberian sardine are currently in serious decline (ICES, 2018(ICES, , 2019b. European anchovy, in contrast, has been at a sustainable level since 2010, with an overall increasing trend that reached its maximum in 2019 (ICES, 2019a). The importance of prey species in common dolphin diet has been found to be related to their availability in terms of abundance (Santos et al., 2004;Meynier et al., 2008), which could explain the negative responses shown by species with low abundances (e.g., Iberian sardine and Atlantic horse mackerel) and the positive and larger contribution in terms of AIC made by those species with higher abundance (i.e., European anchovy). Blue whiting, on the other hand, did not seem to be relevant in explaining the variability of common dolphin over the study period, despite being more abundant than, for example, anchovy or mackerel. Evidence of blue whiting in the diet of the common dolphin was found in the BoB in the 1980s (Desportes, 1985), which could mean that it was important in the past but less so now, or that it is only important, given its poorer energetic condition (4.4 kJ g −1 ), in the absence of other remarkable prey species (Santos et al., 2013). However, not all potential prey species were included and differences in the distribution of stocks may have also affected the results. In fact, only anchovy's biomass had been estimated exclusively for the BoB. Remaining species biomasses were either estimated using adjacent areas (i.e., Iberian sardine) or distribution areas that extended considerably the observations range of common dolphin (i.e., blue whiting, mackerel and in a lesser extent horse mackerel), which could have contributed, for example, to the higher prominence of anchovy detected in this study.

Distributional Shifts
The common dolphin is considered a warm-temperate species, and accordingly, its range is expected to expand in response to increasing water temperature (MacLeod, 2009). This northward expansion seems to be already happening, at least at the northern limit of the species range, as evidenced by a higher frequency of strandings and sightings in northern Britain and southern Scandinavia (MacLeod et al., 2005;Evans and Waggitt, 2020). The BoB, however, does not constitute a range edge within common dolphin's distribution, which can explain why we did not find a northward shift in its CoG, but instead, switches between the inner (i.e., southeast) and the outer (i.e., northwest) part of the BoB. This pattern has also been detected when forecasting the future distribution of anchovy's egg density in the BoB for spring (Erauskin-Extramiana et al., 2019a) and was associated to the contraction (southeast) and expansion (northwest) of anchovy population (Motos et al., 1996). A prey driven distribution was already suggested for albacore tuna in the area (Lezama-Ochoa et al., 2010), so we could hypothesize that the distributional shifts of common dolphins in the BoB are also driven by the distribution of their main prey. Similarly, the increase in common dolphin abundance detected in the southeast corner of the BoB could be also related to a higher prey availability. Indeed, other important prey species of the diet of common dolphin (e.g., horse mackerel, sprat) also shifted to the southeast of the BoB in the past 30 years (Baudron et al., 2020).
The prey variables considered in this study, however, could not explain much of the observed spatio-temporal variability of the CoG as a result of being introduced as a biomass index that changed across time but not across space, and hence, could not confirm or reject the hypothesized prey-driven distribution. Whether top predator abundance and distribution is driven by the environment or prey is a much debated question in ecology (Grinnell, 1917;Elton, 1927;Torres et al., 2008). However, acquiring co-occurring top predator and prey data in space and time to test these hypotheses is challenging. In this study, we have taken advantage of a large spatio-temporal compilation of top predator sightings, but in contrast, we have only been able to incorporate annual, non-spatial biomass indices of prey. Future work, therefore, should focus on improving prey data inputs to better understand their role in driving top predator distributional shifts in the BoB, a question that remains open. Climate indices, as for prey biomasses, were regional time-series rather than spatio-temporal datasets (i.e., changed across time but not across space), so their effect on the CoG is also difficult to understand. Local oceanographic variables did account for spatio-temporal changes, but even so, only explained a very small proportion of spatial shifts, which means that most of the distributional shifts occurred due to unidentified sources. This inability to attribute a source to distributional shifts was also found in previous studies with fishes Perretti and Thorson, 2019), and suggests that more effort must be made to understand when distributional shifts can be attributed to covariates in spatial random effects models (Hodges and Reich, 2010).

CONCLUSION
Climate change is believed to affect marine mammals through changes in their physical environment but also in their prey. However, many studies aimed at understanding climate impacts often employ environmental characteristics as proxies for prey distribution. In this study, we incorporated both environmental and prey variables estimated at local and regional scale and explored the relative importance of each of them in explaining the spatio-temporal variability in common dolphin data. Although we could not attribute much of the detected distributional shifts to the variables considered in this study, we could conclude that, in the BoB, climate indices and prey species biomasses can play an important role in driving the abundance patterns of the common dolphin. Further research on climate change effects on common dolphin, however, should focus on comprising the whole distribution range of the species, given the increasingly feasible possibility for combining surveys across areas and regions provided by methods such as those used here. This way, we could address important knowledge gaps that have not been solved here, for example, if the increasing trend found in abundance is due to the arrival of new individuals or it is the result of an overall population growth. Answering to this question will undoubtedly help understand population dynamics and bycatch implications, but meanwhile, we reiterate our call for caution when interpreting the abundance patterns predicted in this study.

DATA AVAILABILITY STATEMENT
Any requests for survey data should be addressed to their owners. In future, some survey data may become open access. Please contact PE (peter.evans@bangor.ac.uk) for further details.

AUTHOR CONTRIBUTIONS
AA, MLo, and GC conceived and designed the study. AA applied the models and wrote the main text of the manuscript. JW collated and standardized survey data while MP and JT helped with the modeling approach. The remaining co-authors contributed survey data and revised the manuscript. All authors contributed to the article and approved the submitted version.   Supplementary Table 1), to the Bangor University, James Waggitt and Peter Evans for allowing a short visit, and to Rutgers University and the Pinsky Lab for hosting an academic exchange. We also thank Megan Ferguson and Leire Citores for their feedback and help as well as the editor and the two reviewers for their valuable suggestions and comments. This is contribution 1076 from AZTI Marine Research.