Data on the relative abundance and patch characteristics of prey were collected concurrently during the dolphin-habitat surveys using hydro-acoustics (Simrad ESS7 echosounder, see Brough et al., 2019b for detailed methods). The relative abundance of Hector’s dolphins showed a strong relationship with availability of potential prey represented by epipelagic schooling species (e.g., slender sprat – Sprattus antipodum, yellow-eyed mullet – Aldrichetta forsteri) and the depth of fish schools (Brough et al., 2019b). Thus, these variables were included as candidate variables in the dolphin-habitat models used to investigate the drivers of distribution.

Fine scale data on seafloor habitat types for each hotspot and reference area were available from seafloor mapping studies carried out using side scan sonar and single beam echosounders (Brough et al., 2018). Physical habitat types were classified to level 2 (abiotic habitats) using the proposed NZ thematic marine habitat classification (Dohner, 2013) and were ground-truthed using a combination of Ponar sediment grabs and drop-camera deployments. Following Sanciangco et al. (2013), a modification of the Shannon-Weiner index (Shannon and Weaver, 1963) was used to represent habitat heterogeneity calculated as an index (0-1) based on the number and aerial extent of seafloor habitat types encountered during a survey. High resolution bathymetric data (5 m x 5 m) from surveys of the Banks Peninsula nearshore environment using single-beam echosounders (Brough et al., 2018) provide fine scale information on depth, slope and aspect of the seafloor for each survey area, the latter two variables being calculated using the package raster (Hijmans, 2019) in R (R-Core-Team, 2018).

Data on the velocity of tidal currents throughout the study area were available from a regional oceanographic modelling system (ROMS) constructed for Banks Peninsula (Soutelino and Beamsley, 2015). ROMS data were used to generate mean velocity curves across the semi-diurnal tidal cycle for each survey area. A single, survey-specific tidal current velocity value was extracted using the recorded tidal state at the time and location of the mid-point of each survey. Data on local oceanographic conditions at each location at the time of a survey were derived from several sources. Sea surface temperature (SST) datasets were collected alongside hydro-acoustic data in the dolphin-habitat surveys using an onboard temperature sensor within the hydro-acoustic system. The sensor was calibrated every two weeks using an RBR Concerto (RBR Ltd). Temperature at the transducer face was down sampled to 1 Hz, and gridded into 100 m x100 m bins along the survey track. From this grid, the mean, standard deviation (SD) and horizontal gradient in SST was calculated for each survey.

Remote sensing datasets on Chlorophyll a concentration and diffuse attenuation coefficient (KD490) were accessed as level 3 data products of the MODIS Aqua sensor from NASA’s ocean color website (NASA Goddard Space Flight Center et al., 2021). Both variables are regularly used in species distribution models for marine megafauna as proxies for primary productivity (e.g., Torres et al., 2008; Hazen et al., 2012) and water turbidity (e.g., Gomez and Cassini, 2015; De Wysiecki et al., 2020) respectively. Daily datasets were downloaded at 4 km resolution and were processed in SeaDas (v 7.5.1; OBPG) and in R package raster (Hijmans, 2019) to provide pixel values for the closest cell to the centroid of each dolphin-habitat survey, at the nearest possible date. Surveys were carried out with variable, but high frequency during the summer and winter seasons; often with multiple surveys of a given survey area during a week. Due to this high frequency and the fact that conditions change rapidly in this exposed, temperate habitat, remote sensing data with high temporal resolution (i.e., daily) were necessary. The date of variable availability depended largely on local weather and sea conditions but was typically within three days of the survey.

A database that summarized the physical habitat types, bathymetric variables, current velocity, oceanography, prey field characteristics, survey effort and relative abundance of dolphins for each survey was generated. Data on the physical habitat types and bathymetry were extracted from feature and raster layers in ArcMap (v. 10.3; ESRI). A digitized track-line of each survey was imported and a buffer with a 400 m (i.e., max sighting distance; Brough et al., 2019a) radius placed around the track. The resultant ‘survey polygon’ reflected the total area surveyed for each dolphin-habitat survey. The total area and percent coverage of each physical habitat type with the survey polygon was calculated. The mean and SD of the depth and slope bathymetric surfaces that overlapped a survey polygon were calculated and exported for each survey. A velocity value for the location specific tidal state was extracted from the tidal velocity models using R package raster (Hijmans, 2019). Variables derived from remote sensing datasets were imported into the database and matched against the surveys. Every survey had in-situ measures of prey relative abundance and mean school depth. The number of dolphins (n dolphins) encountered during a survey was the measure of dolphin relative abundance used for modelling. A broad range of habitat characteristics were included as candidate variables ( Table 1 ), with multiple variables often representing single characteristics. For example, depth is represented by both mean depth and the standard deviation of depth which provides a measure of seafloor unevenness. In several instances, multiple variables from a single characteristic were correlated, which led to the removal of variables in the concurvity exercise (see below).

To confirm differences in relative abundance between hotspots and reference areas, a generalized linear model (GLM) was fit in R with the number of dolphins encountered at each survey as a response and the survey area ( Figure 1 ) as a single categorical predictor. The GLM was fit with a negative binomial distribution and incorporated the distance of each survey as an offset. Tukey’s post-hoc tests were used to test the difference in relative abundance between each pairwise comparison of the four hotspot areas with the six reference areas (n= 24 pairwise comparisons). Tukey’s post-hoc tests were implemented in the R package multcomp (Hothorn et al., 2008).

The relationship between the relative abundance of dolphins (n dolphins) during each survey and the physical and biological characteristics of habitat were modelled under a species distribution modelling (SDM) framework using generalized additive models (GAMs; Hastie and Tibshirani, 1987). GAMs were fit using the package mgcv (Wood, 2017) in R with a negative binomial distribution. Correlation (linear and non-linear) among environmental covariates was investigated using the concurvity function in mgcv (Wood, 2006; Wood, 2017) using a threshold of 0.3 above which terms were considered correlated (He et al., 2006; Bennington et al., 2020b). When concurvity between variables was apparent, the decision of which to retain was based on the AICc values returned by single variables GAM models that fit the variable in question against the response (n dolphins, with an offset for survey distance). The variable that produced the model with the lowest AICc value was retained. A global model containing all uncorrelated covariates was then fit in mgcv. Variables included in the global model were: Prey_RA, Sand, Perc_Reef, Velocity, Slope, SST, SST_Grad, Turbidity and Env_Het (see Table 1 ). Each continuous variable was fit as a smoothed function using penalized thin plate regression splines, and the degrees of freedom for each term were limited to maximum of 5 (i.e., 4 knots; (Tepsich et al., 2014; Rayment et al., 2015) in order to minimize the chances of overfitting impacting model selection and interpretation. A model set that contained a formulation for every possible combination (and inclusion/exclusion) of candidate variables contained within the global model was constructed and each model (n=2038) fit separately. Model selection was undertaken using AICc (Hurvich, 1989) and model weight scores (Burnham and Anderson, 1998), with the ‘top-model’ being the formulation with the lowest AICc and highest model weight. The importance of habitat variables was determined in a hierarchical approach using three steps; 1) Variables retained in the top-model were considered to have some effect on the relative abundance of dolphins, 2) The strength of a variable’s effect was visually identified by viewing the values of the y-axis for areas on GAM partial dependence plots where the 95% credible interval does not overlap zero. The values on the y-axis of such plots indicates the strength of the effect on the response variable (Wood, 2006). 3) As a further test, the approximate significance of each smoothed-term was extracted from the top-model. Variables with p-values < 0.05 are considered to have a significant (α=0.95) influence on the relative abundance of dolphins.

The global and top-models were checked for the validity of assumptions using a combination of diagnostic plots and the gam.check function within mgcv. The extent of spatial auto-correlation among input observations (i.e, surveys at the same location) was investigated by generating spatial correlograms and bootstrapped 95% confidence intervals (Zuur et al., 2009). To assess goodness of fit and the predictive capacity of the top-model the deviance explained was summarized and a validation procedure established. For model validation, the full dataset was randomly split into a ‘training’ set that contained two thirds of the total dataset and an ‘evaluation’ set containing the remaining third. The top-model was fit using the training data and used to predict values on the response scale given the values of the variables in the evaluation data. The validation procedure was bootstrapped 100 times with different, randomly selected training and evaluation data. The mean Spearman’s correlation coefficient was used as a measure of correlation between the predicted and observed values (Guisan and Zimmermann, 2000; Grüss et al., 2016), with a value higher than 0.50 indicating good predictive capacity (Grüss et al., 2016; Waldock et al., 2022).

To link the habitat covariates deemed important by the GAMs to the existence of hotspots, we used two complementary approaches. Firstly, a ‘GAMvelope’ (Torres et al., 2008) was used to extract the values of the environmental covariates that had a positive relationship with the abundance of dolphins. A GAMvelope describes the area on a GAM curve where values of a habitat variable (x-axis) have a positive effect on the response (y-axis) (Torres et al., 2008; Tepsich et al., 2014; Correia et al., 2015). Extracting these data and examining how they are represented spatially among survey locations allows an appraisal of the location of important habitat (Torres et al., 2008). We extracted all surveys with positive GAMvelope values and used a Chi-squared test to investigate the null hypothesis of no difference in the occurrence of positive GAMvelopes between hotspots and references areas. The test was undertaken for each variable retained in the top-model.

Secondly, principal component analysis (PCA; Jolliffe, 2002) was undertaken with each dolphin-habitat survey as an input observation and the combination of all environmental covariates retained in the top-model used to calculate principal components. The PCA was undertaken using the prcomp function in R and plotted to provide a visualization of any apparent clustering of observations among hotspots and reference areas. The alignment of each variable vector (arrows) with each PCA axis shows the correlation between a variable and the principal components, while the length of each vector is indicative of the contribution of each variable (Jolliffe, 2002). Further, correlations of clusters of observations with the magnitude and direction of vectors of environmental variables in a PCA bi-plot can be used to distinguish the unique habitat characteristics of observations in PCA space (Pitcher et al., 2011; Stephenson et al., 2018).

We integrated intensive surveys of Hector’s dolphins (Cephalorhynchus hectori) at Banks Peninsula with data on environmental characteristics and prey to identify the statistical correlates of dolphin distribution. Among all the variables tested, the relative abundance of potential prey had one of the strongest correlations with relative abundance of dolphins. When prey were more abundant, so were the dolphins. This finding confirms our hypothesis of prey being an important driver of distribution in this species. Many studies have shown the spatial and temporal overlap between marine top-predators and their prey (Baumgartner et al., 2003; Fauchald, 2009; Benoit-Bird et al., 2013; Barlow et al., 2020). Especially for small-bodied marine mammals inhabiting temperate areas, the need for frequent foraging to satisfy high energetic requirements (Harding et al., 2005; Wisniewska et al., 2016) means that these taxa are likely to exhibit particularly high overlap with prey. In studies that combine both prey and environmental variables as drivers of marine mammal habitat use, Embling et al. (2012) and Hazen et al. (2011) found that both prey and habitat covariates were important drivers of abundance, while Torres et al. (2008) found habitat models performed best with environmental covariates only. Such differences among studies are likely due to variation in the ways that habitats and prey are characterized, the foraging patterns of individual species, and the heterogeneity of the underlying community (Torres et al., 2008; Bennington et al., 2020b).

The relative abundance of dolphins had a strong positive relationship with the proportion of sandy substrate in their habitat ( Figure 2 ). Mud and sand were the dominant habitat types in all survey areas and thus the positive association with sand is also a negative association with mud. At fine-scales, harbor (Phoca vitulina) and grey seal (Halichoerus grypus) were also less common where bottom types had a high proportion of muddy sediment (Bailey and Thompson, 2010). Preference for certain bottom types is probably related to the preference of certain prey for these habitats (Torres et al., 2008; Embling et al., 2012), though preferences of Hector’s dolphin prey for specific habitat characteristics are poorly known.

The habitat models suggest that reef habitat was not favored by Hector’s dolphins. In temperate areas, reef habitat and the macroalgae communities that it sustains, are important local sources of productivity, and provide heterogeneity that promotes high biodiversity (Smale et al., 2013; Rees et al., 2014). The strong negative relationship between Hector’s dolphins and reef shows that this species does not respond to these ecological values, at least on the scale of this study; the species rarely targets reef-associated prey (Miller et al., 2013). The abundance of dolphins showed a negative relationship with areas of high habitat heterogeneity ( Figure 2 ), in contrast to the only other study that included heterogeneity as a factor influencing marine mammal distribution (Eierman and Connor, 2014). This finding suggests that Hector’s dolphins may not be a good indicator of habitats with high diversity of seafloor communities. Rather, they appear to favor specific habitat-types (e.g., high-current, sandy habitat) that may support high productivity and prey abundance.

The relative abundance of dolphins was positively correlated with higher tidal current velocity ( Figure 2 ). Several studies have related tidal state to distribution or foraging in marine predators (Johnston et al., 2005; Bailey and Thompson, 2010; Lin et al., 2013; Nuuttila et al., 2017). Brough et al. (2020) found tidal state to influence the onset of foraging bouts in Hector’s dolphin. Harbor porpoise (Phocoena phocoena), an ecologically similar species to Hector’s dolphin, were also found to be more abundant at high tidal stream velocities (Marubini et al., 2009). Tidal currents are hypothesized to promote prey aggregation via the creation of hydrological features such as wakes and eddies (Johnston et al., 2005) and/or enhance catchability of prey (Bailey & Thompson, 2010).

Oceanographic covariates including surface temperature, temperature gradient and turbidity all had some influence on the relative abundance of dolphins. Of these, the influence of temperature was the most pronounced; dolphins were substantially less abundant at temperatures cooler than 12°C (i.e., winter temperatures). Thus, the influence of temperature likely describes the seasonal preference for this near-shore component of the dolphins wider habitat (Rayment et al., 2010). However, the indication of an ‘upper-limit’ for the correlation between high temperatures and higher abundance ( Figure 2 ) should be considered in the context of warming sea conditions in this area (Shears and Bowen, 2017). The GAMs suggest the dolphins are more abundant in less turbid conditions ( Figure 2 ). While it is generally accepted that on a broad scale, Hector’s and Maui (C. hectori maui) dolphins show a preference for nearshore habitat (that is incidentally more turbid) (Bräger et al., 2003; Derville et al., 2016), our study suggests that at fine scales the species selects locations that are less turbid – a result that raises questions about the transferability of modelled relationships between scales of analysis. Turbidity may play a role in the detection of prey, threats or conspecifics (Karczmarski et al., 2000; Heithaus and Dill, 2006). However, turbidity may also act as a proxy for correlated variables including primary productivity and swell exposure (Cloern, 1987; Dittmann et al., 2016).

Several variables deemed important by the SDMs were also more common at hotspots, including habitat with a high proportion of sand, high prey abundance and current velocity ( Figure 3 ). While not statistically significant in the GAMs, the preferred depth range was more frequent at hotspots compared to reference areas ( Figure 3 ). In contrast, preferred values of reef habitat, temperature and temperature gradient, turbidity and heterogeneity were not more common at hotspots ( Figure 3 ). SST, SST gradient and turbidity may be correlated with the dolphins’ broader seasonal distribution in this nearshore habitat, rather than hotspots being locations with particularly high or strong gradients in temperature, or low turbidity. Hector’s dolphins are more common inshore during summer, when waters are warmer and have lower turbidity (Bräger et al., 2003; Rayment et al., 2009), However, turbidity does distinguish a large cluster of observations from reference areas in PCA space ( Figure 4 ) that are likely related to reference areas on the northern coast of Banks Peninsula (LY, ME) and inner harbor locations (WA) that are more turbid. While independently unrelated to hotspots ( Figure 3 ), the combination of reef cover and environmental heterogeneity along with depth, does seem to contribute to areas that are not hotspots; as seen by a discrete cluster in the PCA analysis ( Figure 4 ). It is likely these characteristics combine to form a habitat ‘type’ that is not routinely used by the dolphins and confirms the importance of multi-variate analyses when appraising habitat selection.

The discrepancy between the frequency of occurrence of positive GAMvelopes among hotspots/reference areas and the outputs of the models are likely attributable to the reference area ‘LB’. While this area was not designated a ‘hotspot’ by Brough et al. (2019a), our surveys indicated the area had high dolphin abundance, that was not significantly different from the abundance encountered at hotspots ( Table S1 ). The environmental signature of this location was also more similar to hotspots than reference areas, particularly with regards to prey abundance ( Figure 4 ). It may be that habitat at LB has improved in quality over time. A dataset spanning 30 years was used to initially designate hotspot areas, and thus any very recent preference for this area would be obscured in the analysis by Brough et al., 2019a.

The PCA analysis suggests that not all hotspots and reference areas have the same characteristics. Observations at some reference areas were defined by reef cover, heterogeneity and depth, and others by turbidity ( Figure 4 ). The large spread of the observations from hotspots suggests that these important areas do not have a well-defined ‘recipe’ that describes their constituent characteristics, rejecting our initial hypothesis for a unique environmental signature of these important areas. However, there are some clear contributors to good quality habitat for this species, with high sand cover, prey abundance and tidal velocity clearly separating observations among hotspots and reference areas ( Figure 4 ). While foraging behavior is more common at these hotspots (Brough et al., 2020), it is possible that different prey with different habitat requirements are targeted in different locations. Additionally, other behaviors known to influence habitat selection (e.g., nursing young) may contribute unique habitat preference signals to this appraisal (Weir et al., 2008; Webster et al., 2009). In any case, the complexity in distinguishing a unique environmental signal for hotspots confirms the importance of targeted investigations into hotspot characteristics when habitat-based management is required.

There is good reason to believe that the contributors to hotspot occurrence for Hector’s dolphins may also reflect habitat of high ecological value more generally. The characteristics of seafloor sediment are known to have a significant influence on benthic communities (Thrush et al., 2004). Mud content has been shown to impact benthic invertebrate abundance (Thrush et al., 2003), diversity (Dauer et al., 2000; Ellingsen, 2002) and ecosystem structure and function (Thrush et al., 2003; Jones et al., 2011; Pratt et al., 2013). At Banks Peninsula, mud content has a strong influence on benthic invertebrate fauna, with sand dominated habitat having higher richness and biomass (Fenwick, 2004). The avoidance of areas with high mud content by this apex predator may be due to the impacts of sedimentation of seafloor communities (Thrush et al., 2004) altering the availability of Hector’s dolphin prey (e.g., through benthic foraging fish – red cod, Pseudophycis bachus). If this hypothesis holds true, Hector’s dolphins may be useful indicators of seafloor community health, however further sampling at lower trophic levels is required to confirm such a claim. Sedimentation is a known threat to this nearshore habitat, with some areas experiencing up to 1 m of deposition of fine, muddy sediment over the last 70 years (Hart et al., 2009). If sedimentation is adversely impacting habitat of high ecological value that support Hector’s dolphins, management should target a reduction in this stressor.

High and dynamic tidal flow can promote localized productivity increases due to enhanced mixing and the creation of features such eddies, wakes and fronts, that in turn influence the abundance of mid-trophic species and their catchability to predators (Johnston et al., 2005; Bailey and Thompson, 2010; Embling et al., 2012; Pirotta et al., 2013). Thus, high tidal flow may be a key component of ecologically important habitat in this area. The flow of tidal currents can be impacted by coastal development such as land reclamation and wharf construction (Jefferson et al., 2009), aquaculture (Grant and Bacher, 2001) and dredging (van Maren et al., 2015). These stressors occur at Banks Peninsula and are prevalent throughout coastal New Zealand. Resource management should consider the importance of tidal current velocities for Hector’s dolphins and for contributing to habitat of high ecological value, especially due to the possible increases in coastal development over coming decades.

Mid-trophic level, epipelagic fishes perform a critical role in ecosystems by connecting primary consumers with top predators in so-called wasp-waist ecosystem structuring (Cury et al., 2000). Variability in the abundance of these taxa has been related to distinct changes in ecosystem function and health (Cury et al., 2000; Griffiths et al., 2013). Thus, locations that are routinely used by epipelagic fish, such as the hotspots for Hector’s dolphins investigated in this study, may be considered ecologically important habitat due to their perceived role in supporting high biomass of these taxa. Threats to mid-trophic level fishes are well known and include overfishing, habitat modification, disease, and climate change (Ward, 2001; Bertrand et al., 2004; Roux et al., 2013). There is limited commercial harvest of the small epipelagic species in New Zealand, however changes in oceanographic conditions, via terrestrial impacts (e.g., suspended sediment) or climate change are a distinct possibility.

Increased turbidity from suspended sediments derived from terrestrial runoff is a common feature of coastal habitat internationally (Galy et al., 2015; Halpern et al., 2015). High turbidity can alter rates of primary productivity in phytoplankton (Cloern, 1987; May et al., 2003) and macroalgae (Desmond et al., 2015; Blain et al., 2021), and can influence the health, abundance and survival of juvenile fish (Wilber and Clarke, 2001; Lowe et al., 2015). Thus, similar to the previous environmental variables, the effects of turbidity on marine mammals are likely indirect via prey, though further analysis is required to confirm relationships between prey abundance and variable turbidity.