Using vessels of opportunity for determining important habitats of bottlenose dolphins in Port Phillip Bay, south-eastern Australia

Study area, data collection and processing

The study was conducted mostly in the southern extent of Port Phillip Bay, a large semi-enclosed marine embayment that is located on the northern coast of Bass Strait, south-eastern Australia (Sampson, Easton & Singh, 2014) (Fig. 1A). The bay is relatively shallow (mean depth 13.6 m), with a maximum depth of only 24 m in the center. The south-eastern region of the bay is characterized by steep bathymetric gradients associated with a channel that opens into Bass Strait through Port Phillip Heads (Holdgate et al., 2001). Port Phillip Bay is a busy waterway accommodating various vessel types. Each year ~3,500 cargo and container ships visit the port of Melbourne, making it Australia’s largest container and general cargo port (Sampson, Easton & Singh, 2014). Numerous recreational fishing boats, cruise ships and ecotourism vessels operate within Port Phillip Bay, alongside two passenger ferries and eight commercially licensed longline fishing vessels (Sampson, Easton & Singh, 2014; Victorian Fisheries Authority, 2023). Sightings of dolphins by staff and patrons onboard licensed ecotourism vessels within the south-eastern region of the bay were recorded opportunistically by the operators throughout the austral spring-autumn period (September to April), with these collated records forming the basis for the present study. In addition, data were obtained from the publicly available Victorian Biodiversity Atlas (VBA) database (Victorian Department of Energy Environment and Climate Action, 2024) which includes sightings from vessels of opportunity and researchers. Where multiple sightings were recorded on the same day, only those sightings with a minimum time interval of >15 min or distance >500 m were considered unique presence records. Thresholds were selected based on methods of Lacetera et al. (2023), and DEECA vessel regulations regarding marine mammal approach limits in Port Phillip Bay (Hale, 2002). To identify potential bias associated with repetitive sampling of individuals, distance between unique presence records was further examined. On the same day <20% of all records were within a 500 m distance and were therefore included in further analyses.

Since 2006, information on sightings (Global Position System (GPS) derived latitude and longitude; date and time) have been obtained from vessel logbooks, the VBA and the metadata from digital images of dolphin dorsal fins taken by patrons on ecotourism vessels. Fine-scale data collection methods (e.g., route, speed and duration of vessel tours; sea-state conditions etc.) could not be ascertained accurately for these data. However, all records used in the present study were obtained from vessels that adhered to DEECA regulations regarding marine mammal approach limits for tour operators and research vessels. These included a speed <10 knots (18.5 km·h⁻¹) being maintained when <100 m from a dolphin or pod, and remaining a maximum of 20 min when within a 100 m perimeter (Hale, 2002). In addition, tours are limited by the weather for safety, patron comfort and dolphin visibility, operating only in Beaufort Sea State Conditions <3. Nonetheless, due to the uncertainty of search methods and effort, these data were used in presence-background analyses of habitat affiliations with dolphin distribution.

Environmental variables

Processing of environmental data was conducted using ArcGIS (Version 10.8; ESRI, Redlands, USA). As dolphins are top-order predators it is assumed their distribution is highly influenced by prey availability (Eierman & Connor, 2014). Prey abundance information for Port Phillip Bay is limited and, therefore, environmental variables known to influence their distribution were selected as proxies for habitat suitability models. Dynamic environmental variables such as sea surface temperature and sea surface chlorophyll-a are available only at low spatial resolution across Port Phillip Bay (SST: 0.05° × 0.05° (Good et al., 2020); Chl-a: 0.042° × 0.042° (Sathyendranath et al., 2019)). Although they vary seasonally, the low resolution and variation suggested that these variables were not useful in the subsequent analyses. Instead, static environmental variables that have been shown to influence the distribution of bottlenose dolphins in previous studies (Beddoe et al., 2024; Hanf et al., 2022; Haughey et al., 2021) were incorporated into habitat suitability models. Temperature, light attenuation and the presence of benthic primary producers are influenced by seafloor depth, and areas of upwelling that attract fish and cephalopods are reflected through associated seafloor slope (Eierman & Connor, 2014). Therefore, bathymetry (m) and seafloor slope (% gradient) were included as predictor variables in the models. Bathymetry data for Port Phillip Bay were sourced from a Coastal Digital Elevation Model created by the Cooperative Research Centre for Spatial Information (CRCSI) (2022) with a resolution of 10 m × 10 m and depth values ranging from 0 m (mean sea level) to −80 m. The benthic terrain modeller in ArcGIS was used to calculate slope, aspect and TPI from bathymetry (Walbridge et al., 2018). To obtain perspective of fine-scale gradients, slope was calculated using a focal window of 10 m × 10 m.

Physical topographic features of the seafloor, such as its ruggedness and aspect, can influence ocean currents and tides that effect prey distribution (Bailey & Thompson, 2010). Therefore, seafloor aspect (degrees) and Topographic Position Index (TPI) were additionally included as predictor variables in the models. TPI is a measure that reflects upper, middle and lower features of the landscape, with an index of 0 or near 0 reflecting flat or continuous features, a large positive value reflecting ridges or mounts, and a large negative value reflecting the bottom of a ridge or shelf (Amatulli et al., 2018). The sine and cosine of aspect, representing the relative eastness and northness of downslope direction (Hession & Moore, 2011), were calculated to prevent problematic logistic models. As features of the study area are small-scaled due to the relatively flat landscape of Port Phillip Bay (Holdgate et al., 2001), TPI was calculated based on a cell neighborhood of 5 × 5 grid cells each 10 m × 10 m.

Seafloor habitat type (i.e. sediment, seagrass beds, macroalgae on sediment, biogenic reef, rocky reef and sessile invertebrate beds) directly influences the presence of a diverse range of invertebrate, fish and cephalopod species (Eierman & Connor, 2014). Therefore, benthic habitat was included as a categorical predictor variable in the model. Benthic habitat data were sourced from a Statewide Marine Habitat Map (2023) (10 m × 10 m spatial resolution), within the DEECA CoastKit database (Mazor et al., 2023).

Human activities occurring in Port Phillip Bay have been suggested to influence bottlenose dolphin behavior and habitat relationships (Filby et al., 2017; Scarpaci, Corkeron & Nugegoda, 2003). Hence, two variables reflecting the environment in relation to urbanization and vessel (commercial and recreational) activity were included as predictor variables in the models as Euclidean distance to coastline (m) and Euclidean distance to shipping channels (m), respectively. As previous studies have found these factors to reduce the likelihood of bottlenose dolphin presence in habitat areas (Filby et al., 2017; Zanardo et al., 2017), it is hypothesized that probability of presence will increase with increasing distance (m). These two predictor variables were calculated using geospatial vectors and spatial analyst in ArcGIS.

Raster layers were generated for the extent of Port Phillip Bay with grid cell resolution of 10 m × 10 m. This resolution was chosen to ensure fine-scale sensitivity when analyzing bathymetric features and Euclidean distances.

Ensemble habitat suitability modelling

Statistical modelling of habitat associations with dolphin distribution was conducted in R statistical environment (version 4.4.0) (R Core Team, 2024). An ensemble presence-background habitat suitability model was used to predict potential habitat areas in Port Phillip Bay most suitable for bottlenose dolphins and to assess the importance of environmental variables.

As true/inferred absences and search effort information were lacking in the present study random background points were generated. To reduce the effect of sampling bias, background points were extracted within the model calibration zone (area sampled by ecotourism and other vessels of opportunity). Environmental variables were clipped at a 2 km radius of presence points and 3,000 random background points were sampled within this radius. The model calibration zone was defined considering visibility range relative of pseudo-absence and height of the observation platform. Three background point sample replicates were generated to account for variation within each modelling algorithm. If multiple presence or background points were recorded in a single 10 m grid cell, these were combined to represent a single presence or background point to avoid pseudo-replication in the model (Correia et al., 2021; Elith et al., 2011). Once having removed NA (not available) values from the presence-background data, statistical assumptions could be checked and met.

Before running habitat suitability models, collinearity (high correlation) between predictor variables was assessed using ‘cor’ from the car package (Fox & Weisberg, 2019), and ‘vif’ from the usdm package (Naimi et al., 2014). Variance inflation factors (VIF) and correlations between continuous variables were computed to allow the exclusion of highly correlated variables from the same models, using a linear correlation threshold of 0.7 and VIF threshold of 3 (Haughey et al., 2021). As VIFs and coefficients were below the thresholds for the dataset, all seven predictor variables were included in the same models.

An ensemble habitat suitability model was generated using the BioMod2 package (Thuiller et al., 2009). Predictions from five different presence-background modelling algorithms; flexible discriminant analysis (FDA), generalized additive model (GAM), generalized boosted model (GBM), multivariate adaptive regression splines (MARS) and maximum entropy (MAXENT), were combined. These algorithms were selected as they perform fairly well in their predictive accuracy (Phillips, Dudík & Schapire, 2004; Thuiller et al., 2009) and allow an overall comparison between regression (GAM, GBM, MARS), classification (FDA) and machine learning (MAXENT).

Default algorithm parameters were used with each FDA, GAM, GBM and MARS, as recommended by Thuiller et al. (2009). However, MAXENT parameters were refined using ‘ENM evaluate’ from the ENMeval package (Muscarella et al., 2014), with preliminary analyses having indicated linear and quadratic features to provide the greatest area under the receiver operating characteristic curve (AUC) score and smallest ΔAICc. Model algorithm settings are listed in Material S2. Habitat suitability models were built with a binomial error distribution and ‘logit’ link function. Data was randomly split into model calibration (75%) and evaluation (25%) data, following a 10-fold cross-validation method (Thuiller et al., 2009). To prevent negative prevalence effects, presence and background points were equally weighted with a set prevalence of 0.5.

The AUC metric was used to assess model predictive performance. An AUC score < 0.5 suggests that a model performs worse than random (i.e. counter-predictions). Model performance was determined following (Araújo et al. 2005): AUC < 0.5 = counter-predictions; 0.5 < AUC < 0.6 = fail; 0.6 < AUC < 0.7 = poor; 0.7 < AUC < 0.8 = fair; 0.8 < AUC < 0.9 = good; 0.9 < AUC = excellent. Only algorithms with a mean AUC value > 0.7 were included in further analyses to ensure sufficient prediction accuracy. The importance of environmental variables was determined using a randomized 10-permutation run procedure within BioMod2 (Thuiller et al., 2009). This approach calculates the Pearson’s correlation between standard predictions using all predictor variables and predictions where the variable being evaluated is randomly permutated (rearranged). Low correlation (i.e. significant difference between the two predictions of variable importance) indicates the predictor variable is of high importance, and high correlation indicates the variable is of low importance. According to the mean correlation coefficient, variables are then ranked from 0 to 1, with the highest-ranking variable the most influential and the lowest the least influential (Thuiller et al., 2009).

To generate the final ensemble model, the mean AUC score of all algorithms if >0.7 was calculated. The importance of each environmental variable across the different model algorithms was averaged using a mean of means value. An ensemble habitat suitability map illustrating mean predicted probability and associated variance from model algorithms was generated within the model projection zone. Response curves describing the variation in probability of bottlenose dolphin presence as a response of each environmental variable were plotted.