Fine‐scale spatial and temporal distribution patterns of large marine predators in a biodiversity hotspot

Large marine predators, such as cetaceans and sharks, play a crucial role in maintaining biodiversity patterns and ecosystem function, yet few estimates of their spatial distribution exist. We aimed to determine the species richness of large marine predators and investigate their fine‐scale spatiotemporal distribution patterns to inform conservation management.


| INTRODUC TI ON
Large marine predators, such as cetaceans and sharks, often occupy the top trophic level in marine food webs. They play a crucial role in maintaining biodiversity patterns and ecosystem function, transferring energy as they travel both vertically and horizontally through the ocean, often over large spatial and temporal scales (e.g. Catano et al., 2016;Hindell et al., 2020;Sequeira et al., 2018;Wang et al., 2019). It is therefore important that the distribution and abundance of marine predator populations are understood; yet, given their K-selected life history and high mobility, this requires largescale and long-term monitoring. However, few examples of comprehensive monitoring programmes exist in marine ecosystems and are often retrospective analyses (e.g. Hindell et al., 2020;Sequeira et al., 2018) despite the anthropogenic pressures that threaten the persistence of biodiversity, including large marine predators Schipper et al., 2008;Worm et al., 2006). Anthropogenic impacts are well-illustrated by the poor conservation status of cetacean and shark species, with a decline of more than 70% of sharks and rays since the 1970s  and 37% of marine mammals at risk of extinction (Davidson et al., 2012). Given the importance of top predators in shaping community structure and composition (Bowen, 1997), conservation interventions are required to ensure populations and ecosystems are maintained (Albouy et al., 2017).
Large marine predators that depend on productive coastal habitats are at risk because their habitat use often overlaps with areas of economic, social or cultural value Block et al., 2011;Pichegru et al., 2009). Area-based management tools, such as marine protected areas (MPAs), are a popular conservation strategy for managing human activities in ecologically significant marine environments (Davidson & Dulvy, 2017;Grorud-Colvert et al., 2021). MPAs are generally spatially explicit, static units, with varying levels of protection ranging from no-take marine reserves to multi-use MPAs (Kelleher, 1999;Sciberras et al., 2013). Regardless of their conservation aims, the successful use of MPAs for the management of ocean ecosystems depends on knowing what species inhabit the area and their spatial and temporal distributions (McClellan et al., 2014;Morzaria-Luna et al., 2018).
Information on large predator distribution and density is required for ecosystem-based management (Bossart, 2011;Zacharias & Roff, 2001), yet despite their iconic status and potential to act as ecosystem indicators, distribution data for many marine predators is poorly known. Taxa such as cetaceans and sharks are challenging to study because of their vagile and cryptic nature and the financial and logistical difficulties of large-scale and long-term studies.
Even where data exist, large marine predators are difficult to manage spatially as they can cover vast distances between patches of critical habitat that vary at multiple temporal and spatial scales Hindell et al., 2020;Sequeira et al., 2019). In addition, critical habitat patches may be relatively stable or highly dynamic depending on physical, climatic and oceanographic features determining the availability and distribution of resources (Cox et al., 2018;Lambert et al., 2017;Murphy, 1995). As a result, static ocean management tools may not be an effective approach to ensuring the protection of highly mobile marine predators and their dynamic habitat (e.g. Birkmanis et al., 2020;Hartel et al., 2014;Lodi et al., 2020).

Models capable of forecasting suitable habitat show promise in mov-
ing towards time-varying management tools, such as dynamic MPAs, to protect large marine predators, especially highly mobile species (Barlow & Torres, 2021;Hausner et al., 2021;Hobday et al., 2010).
Spatial Distribution Models (SDM) are one tool that can be used to explore species habitat preferences and predict distributions.
SDM can identify relationships between species occurrences and environmental covariates to describe habitat preferences, then translate the output into geographical space to predict distribution within and beyond the study area . Some predator distribution studies include prey (biotic) data as covariates in SDMs Torres et al., 2008). However, as these data are often challenging to collect at the same temporal and spatial scale as the study species, most SDM studies include environmental covariates only (e.g. depth, salinity), with the assumption that environmental conditions drive resource distribution (Cañadas et al., 2002;Finucci et al., 2021;Skov et al., 2008;Virgili et al., 2021). Understanding fine-scale details about the habitat in which the study species live and the drivers that may influence their distribution is near impossible for most ocean environments, but broad-scale environmental variables can inform models, as seen in our study site.
The Hauraki Gulf/ Tīkapa Moana/ Te Moananui-ā-Toi, Aotearoa New Zealand (henceforth the Gulf), is a highly productive marine environment that supports a diverse community of resident and migratory large marine predators (Gostischa et al., 2021;Whitehead et al., 2019). It is a seasonally stratified wind-driven coastal upwelling marine ecosystem that sits on a wide shelf (~60 km ;Stevens et al., 2019). The topography of the Gulf is relatively homogenous in depth (mostly 20-70 m) and is characterized by a gentle slope to the edge of the continental shelf. The presence of multiple islands, embayments and harbours supplies the diversity necessary to host 46 habitat types (Jackson & Lundquist, 2016). The Gulf has some of the highest primary production levels in New Zealand due to the joint effect of seasonal upwelling, a western boundary sub-tropical current (East Auckland Current) and the interactions between local surface winds, currents and rather complex geomorphology of the embayment (Zeldis et al., 2004). The Gulf is an area of high economic, social and cultural value with recreational and commercial fisheries and the country's largest commercial port. Indicators of ecological quality have shown ongoing environmental degradation, despite the declaration of the Hauraki Gulf Marine Park in 2000 (Hauraki Gulf Marine Park Act, 2000;Hauraki Gulf Forum, 2020), and actions to protect this habitat are urgently needed. This requires an understanding of the biodiversity and drivers behind organisms' habitat use and distribution. Despite the importance of the Gulf for marine biodiversity and large marine predators, few estimates of the spatial distribution of these animals exist (Dwyer et al., 2016), with no distribution studies of large shark species.
To fill this knowledge gap, we conducted a replicate systematic aerial survey of the Gulf over a full year to determine the species richness of large marine predators, the Bryde's whale (Balaenoptera edeni brydei), common dolphin (Delphinus delphis), common bottlenose dolphin (Tursiops truncatus), bronze whaler shark (Carcharhinus brachyurus), juvenile smooth hammerhead shark (Sphyrna zygaena) and the species grouping of pelagic sharks, and investigate their fine-scale, spatiotemporal distribution patterns. Specifically, we used boosted regression tree (BRT) models, a flexible machine learning method, to identify dominant environmental drivers and biotic drivers (prey) of species distributions and generate monthly predictive maps, illustrating the variation in the likelihood of occurrence and species richness across the study area and months. We then highlight the utility of this information for conservation planning, particularly the importance of understanding spatiotemporal distribution patterns for dynamic ocean management.

| Survey design and protocol
Twenty-two double-observer line-transect aerial surveys were conducted typically twice a month from November 2013 to October 2014, following MacKenzie and Clement (2014). Surveys followed a systematic grid of eight parallel transects evenly spaced apart every 10 km, orientated offshore at 62° from the north (Figure 1). Under this design, it was possible to sample across a wide range of environmental conditions. The survey followed a simple, unstratified design due to a wide range of species with different or unknown distributions (Dawson et al., 2008). To ensure equal coverage probability (i.e. a random sample of the total habitat area in the survey area; Buckland et al., 2004) the start point of each survey was randomly chosen from two locations, either the transect line starting at the top northwest of the grid or at the bottom southwest of the grid, using the striplet method (Fewster, 2011). For cetaceans, we collected data following distance sampling protocols (Buckland et al., 2004) as this study was part of a larger project that also estimated cetacean abundances (Hamilton et al., 2018). For all other species, sightings data were collected using the strip transect methodology (Eberhardt, 1978). Surveys commenced when the Beaufort sea state (BSS) was ≤3 across the study area, there was no rain or fog and there was sufficient light, i.e. 1 hour after sunrise and before sunset. records of sightings within 80-90°; however, none of the sightings made by observers on the opposite sides of the plane were a potential duplicate as they never recorded the same sighting. While the bubble window allowed the rear observers to search for animals on the track line, the reality is that their viewing area did not overlap, and it is therefore highly unlikely a single animal would be detected by both observers. The observers operated independently with no communication during the survey. They logged all sightings of large marine predator species, i.e. cetaceans, pinnipeds, sharks, rays and oceanic fish observed during the flights.
For each sighting, the observer recorded (1) time of the sighting to the second; (2) species or the highest taxonomic level possible; (3) number of individuals or, where they could not do an absolute count, an estimate of the minimum, best and the maximum number of individuals; (4) for cetaceans, the perpendicular angle of the sighting to the plane was measured using an inclinometer; and (5) environmental conditions: the BSS, glare direction, glare coverage and water colour.

| Presence records
Due to the double-observer design, records partly consisted of duplicate sightings, i.e. two records of the same animal or group made by observers seated on the same side of the plane.
Duplicate sightings were reconciled post-hoc by comparing all records of the same species made by observers seated on the same side of the plane during the same survey. For cetaceans, sightings made within 7 s and within 10° of each other were considered a duplicate sighting. Sharks and rays were generally encountered at low densities allowing for easy identification of duplicate records.
We therefore applied a narrower time frame of 5 s (and within 10° of each other) to classify records of sharks and prey and avoid duplication.
Sightings of sharks exceeding an estimated total length of 2.5 m (Last & Stevens, 2009) that share a similar ecological role were combined to form the species group 'pelagic sharks' as they typically inhabit similar habitats (pers. comm. Clinton Duffy, Department of Conservation). These included rarely-sighted species, i.e.

| Generation of pseudo-absence data
To estimate the distributions of the large marine predators, binomial BRT models require locations of both presences and absences (Manly et al., 2007). Presence records consisted of sightings recorded during aerial surveys (Table 1 for details). Given that our sightings only reflect a portion of the animals' habitat use (i.e. when the animal is at, or close to, the surface of the water), presence records will likely be an underestimate of true distribution. In line with this, absence information collected from our aerial surveys should therefore also not be considered true absences. Instead, pseudo-absences (artificial absence points) are generated (Lobo & Tognelli, 2011;Phillips et al., 2009).
When data are collected during systematic surveys (as was the case here), pseudo-absences are distributed within the absence zones, i.e. where no sightings were made (Derville et al., 2016).
To capture habitat use at a fine temporal (1-day) scale (Derville et al., 2016), 22 unique sets of pseudo-absences were generated for each species, corresponding to each survey. Pseudo-absences were distributed in the on-effort portions of transects, i.e. where observers were in search mode and within the viewing strip as per distance sampling methods (Hamilton et al., 2018; 420 m transect width for sharks and dolphins; 570 m for Bryde's whale).
To ensure that pseudo-absence points reflected available but unused habitat in a nonbiased manner, sightings collected over the entire period were pooled to create a kernel density map for each taxon using the Spatial Analyst kernel density tool in ArcMap Bryde's whales) to prevent environmental overlap between presences and pseudo-absences-the length approximates to the 1-km spatial resolution of the study (Torres et al., 2008). We generated stratified random pseudo-absence points with numbers inversely related to the density contours. Pseudo-absences were distributed with a minimum of two times the transect width to avoid serial autocorrelation (Derville et al., 2016). To weight pseudo-absences according to the stratified design, we created a standardized, inverted kernel density map and extracted the kernel density estimate (KDE) of the cell within which a pseudo-absence point fell to assign a weight. The process produced approximately 30 times more pseudo-absences than presence points for sub-sampling during the modelling process.

| Environmental and biotic variables
Spatial estimates of 13 high-resolution environmental and biotic variables (1-km grid resolution) that influence the distribution of sharks and cetaceans (Redfern et al., 2006;Schlaff et al., 2014;Stephenson et al., 2020;Tardin et al., 2017) were collated (Bostock et al., 2019a(Bostock et al., , 2019bKozmian-Ledward, 2014;Pinkerton et al., 2018;Stephenson et al., 2020; Table S1). Some variables were static (no temporal variation, e.g. bathymetry), whereas others were dynamic (monthly temporal variation, e.g. sea surface temperature-SST). Two oceanographic variables-seafloor sediment disturbance from wave action and tidal currents-were only available as annual averages and therefore were considered static, despite knowing that these vary through time. Biotic variables included the distance of species' sightings to zooplankton and fish sightings, representing prey, acknowledging that a complete lack of spatial coverage across the study area can result in less certain spatial predictions where observations were not made (Bowden et al., 2021). To account for changes over time in the occurrence of the large predators not captured by the available environmental and biotic variables, 'month' was included as a factor in the BRT models (Compton et al., 2012).

| Boosted regression tree model fitting and evaluation
Relationships between large marine predator presence / pseudoabsence and environmental and biotic variables were investigated using BRT models. BRT modelling combines many individual regression trees (models that relate a response to their predictors by recursive binary splits) and boosting (an adaptive method for combining many simple models to give improved predictive performance) to form a single ensemble model (Elith et al., 2008). Detailed descriptions are available in Ridgeway (2007) and Elith et al. (2008). All statistical analyses were undertaken in R (R Core Team, 2020) using the 'dismo' package (Hijmans et al., 2017).
Habitat models were built using a portion of the pseudo-absence dataset. Points were randomly selected until the combined weight of the pseudo-absences was equal to the combined weight of the presences, where presences were given a weight of 1 and pseudoabsences were given the weight of the KDE. BRT models were fitted with Bernoulli error distribution with a bag fraction of 0.6, a learning rate (the shrinkage parameter) of 0.001-0.0001 (to fit between 1000 and 2500 trees for each species' model) and a tree complexity of 3 (the number of interactions fitted) using a random 5-fold cross-validation.
In general, BRTs are reasonably robust to correlation between predictor variables (Charlène et al., 2020;); though highly correlated variables can complicate the interpretation of the results, and their inclusion may not greatly increase the predictive power (Leathwick et al., 2006). Here, multiple steps were undertaken to produce parsimonious models for each species (i.e. select which predictor variables would be used in the final models). First, variables with Pearson's correlation coefficient >0.85 were removed (i.e. as per Elith et al., 2010). Distance to 40-m depth contour was consistently correlated with bathymetry and slope (i.e. >0.85) and was therefore removed (Figures S1-S6). For each species and species group, a BRT model was then initially fitted using all environmental variables (bar Distance to 40-m depth contour), which was then subjected to a simplification process whereby predictor variables were TA B L E 1 Total number of sightings of species or species group by month recorded during aerial surveys from November 2013 to October 2014. removed from the models, one at a time, using the 'gbm.simplify' function (Elith et al., 2008). Model simplification can be beneficial for small datasets as including more predictor variables than necessary may degrade the quality of model performance (Elith et al., 2008).
This process first assesses the relative contributions of each variable in terms of deviance explained, with the lowest contributing variables removed from the model. The model is then refitted with the remaining environmental variables. Model simplification was stopped when model performance (deviance explained) was decreased by >1% (Leathwick et al., 2006;Stephenson et al., 2020). Predictor variables for the final 'simplified' models varied by species (see Section 3).
Finally, the effect of spatial autocorrelation (SAC) was evaluated BRT models were assessed using measures of model performance, including the deviance explained, area under the receiveroperating characteristic curve (AUC) and Boyce index. The explained deviance provides a measure of the goodness-of-fit between the predicted and raw values (total deviance) (Compton et al., 2012).
Acknowledging that there are uncertainties associated with using pseudo-absences for assessing model fits, we use the Boyce Index that only requires presences and measures how much model predictions differ from the random distribution of the observed presences across the prediction gradients (Boyce et al., 2002).
The relative influence of each environmental variable was assessed in the models was the number of times it was selected for splitting, weighted by the squared improvement to the model as a result of each split (using in-bag data) (Friedman & Meulman, 2003).
The association between species and species group presence/ pseudo-absence and environmental predictor variables was illustrated using partial dependence plots (i.e. predicted response curve of species probability of occurrence across the gradient of the variable of interest when all other variables are held at their means, and the median factor level is used for factors).

| Spatial predictions and measures of uncertainty
BRT models were bootstrapped 100 times for all species using the same parameter settings as those used in the simplified models. For each bootstrap iteration, a random sample of sightings was drawn with replacement and random samples of pseudo-absences were selected until the weighting of the pseudo-absence (from the KDE) was the same as the overall weighting of sightings available for each species. Sampling of pseudo-absences was drawn without replacement.
At each bootstrap iteration, the model fits were estimated using both 'training data' (the randomly selected data used to tune the model) and 'evaluation' data (the remaining data not randomly selected from the present data and a random pseudo-absence selection as detailed above). Model fits calculated using evaluation data are presented here because these are considered a more robust and conservative method of evaluating goodness-of-fit (Friedman et al., 2001).
For each species, bootstrapped BRT models were predicted geographically for each month using mean monthly estimates (derived from daily estimates) for those variables that were considered dynamic and the static variables. The mean probability of occurrence and the coefficient of variation (Anderson et al., 2016) were calculated for each 1 km 2 cell to produce 12 monthly habitat prediction maps and 12 spatially explicit maps of uncertainty, respectively.
Monthly distribution of species richness was calculated by summing the monthly occurrence probability predictions (ranging from 0 to 1) from individual models (six BRT models; Calabrese et al., 2014;Ferrier & Guisan, 2006). Finally, a mean annual occurrence probability was produced for each species and species group by averaging each species' and species groups' monthly predictions.

| RE SULTS
Following data quality control, 449 sightings of cetaceans and sharks were collected during 22 aerial surveys conducted from November 2013 to October 2014 (see Table S2 and Figures S7-S14 for details).
These data were used to generate BRT models for Bryde's whale, common dolphin and bottlenose dolphin, bronze whaler shark, the species grouping of pelagic sharks and the juvenile hammerhead shark.

| Model performance
Based on mean model fit measures (averaged across bootstrapped models), the deviance explained by BRT models ranged from a mean of 18% (standard deviation of the mean: ±0.09%) for Bryde's whales to a mean of 55% (±0.06%) for immature hammerhead sharks (

| Predictor variable selection and contribution
The number of predictor variables retained for species' models ranging from six (Bryde's whale) to nine (common dolphin) (

| Bronze whaler shark
The highest probabilities of occurrence for bronze whaler sharks were in January, with the most suitable habitat predicted throughout the inner Gulf, the Colville Channel, around the outer Gulf islands and in pockets along the Coromandel Peninsula coastline (Figure 2).
The probability of occurrence was lower during other months, other than a peak in the cooler months of April, May and June. Uncertainty in predictions was lowest in January and highest around prey patch occurrences (Figures S17 and S18).

| Pelagic sharks
Probability of pelagic shark occurrence was high throughout the majority of the outer Gulf in January, particularly in the northwest and on the east coast of the Coromandel Peninsula ( Figure 2). Probability of occurrence was moderately high in May and June, but these areas were not as extensive as those in January. Prediction maps mostly showed a low to moderate likelihood of occurrence during the remaining months. Uncertainty around spatial predictions was TA B L E 2 Mean estimates (± standard deviation, SD) of model performance (using evaluation data) for the bootstrapped boosted regression tree models fitted with presences/pseudo-absences per species or species group.

TA B L E 3
Mean environmental and biotic predictor contribution (as a percent ± standard deviation) for the 100 bootstrapped boosted regression tree models fitted with presences/pseudo-absences per species or species group. Environmental and biotic predictors are grouped into static (orange), substrate (grey) and dynamic (blue) categories.

Common dolphin
Bottlenose dolphin

| Immature hammerhead sharks
Immature hammerhead shark probability of occurrence was highest in January and March, and highly suitable habitat was predicted in coastal habitats within the inner Gulf and around Great Barrier-Aotea and Little Barrier Islands-Hauturu (Figure 2). Like the other shark models, the likelihood of occurrence was low in February.
There were low predicted occurrences from April to December across the study area, broadly spanning the cooler water months.
Uncertainty in spatial predictions of immature hammerhead distribution was low in coastal regions with high predicted probability of occurrences, i.e. in January and March. However, uncertainty in spatial predictions was highest in coastal regions (compared with the deeper, mid-inner Gulf) for the remaining months with few or no sightings of juvenile hammerhead sharks (Figures S17 and S18).

| Bryde's whale
There was variability in the seasonal patterns of Bryde's whale monthly distributions with the highest predicted occurrence observed from August to January and lowest from March to July October, a large portion of the study area was predicted to be suitable habitat. Probability of occurrence was lowest from May to July (<0.4). The degree of uncertainty around the spatial predictions was relatively similar across months; however, uncertainty around spatial predictions was dynamic, likely due to prey (Figures S17 and S18).

| Common dolphin
Common dolphins were predicted to occur year-round with moderate to high probabilities of occurrence (0.5-1.0); the highest probability of occurrence was predicted in August to November (Figure 2) when common dolphins were predicted to occur in the eastern region of the inner Gulf and around Great Barrier Island. The geographical extent was variable across months, ranging from a highly localized to broad distribution. Higher mean occurrence probabilities were associated with a lower standard deviation ( Figures S17 and S18).

| Bottlenose dolphin
Spatial predictions of bottlenose dolphins were concentrated in coastal habitats along the mainland and islands throughout the Gulf ( Figure 2). The likelihood of occurrence peaked from November to March and was lowest from April to July. Prediction maps showed an expansion and contraction in habitat preferences and a higher degree of connectivity between coastal habitats when the overall likelihood of occurrence was high. The degree of uncertainty in spatial predictions was relatively similar across months (Figures S17 and S18). However, there was a higher degree of uncertainty around spatial predictions associated with prey patches and at greater distances from the coast.

| Mean annual distributions
Across study taxa, the mean annual probability of occurrence was moderate to low, reflecting that most taxa are only predicted to be present during a portion of the year (Figure 2). For taxa with relatively consistent distributions of moderate to high probability occurrence over time (e.g. Bryde's whale, common dolphin, bottlenose dolphin), the distributions of mean annual probability were the highest in similar 'core' areas as those identified in the monthly predictions, albeit with much lower probability occurrence values. By contrast, for taxa with either very few months with high predicted probability occurrence or which were more variable in the locations of these high probability areas over time, the distributions of mean annual probability were very low and spread out across large parts of the study area (e.g. bronze whaler shark, pelagic shark and immature hammerhead shark).

| Monthly richness
Monthly predicted richness of large marine predators was the highest in the late spring and summer months: November-March and peaking in January (mean richness >2, Figure 2); areas with the highest richness were predicted in coastal habitats along the mainland and islands in the inner and outer Gulf during these months. The remaining months were predicted mainly to have low richness (<1), although large parts of the study area had predicted richness values 1-2 primarily in central parts of the study region.

| DISCUSS ION
This investigation of a known global biodiversity hotspot for large marine predators, the Hauraki Gulf (Hauraki Gulf Forum, 2020), gives unique insights into sharks and cetaceans' fine-scale spatiotemporal distribution patterns. The comprehensive coverage of the Gulf via aerial surveys enabled us to compile a year-round understanding of the importance of these waters to a range of large marine predators. Temporal changes in distribution are not routinely accounted for in species distribution studies, and even fewer studies also account for temporally variable biotic variables (prey) (see Barlow et al., 2020;Bennington et al., 2020;Torres et al., 2008). We explicitly examined temporal changes of species distribution patterns; spatial and temporal overlaps of our study taxa across large areas in summer illustrate the importance of the Gulf for large marine predators. However, there were also apparent differences in distribution over time suggesting that a generic management plan would be inadequate for protecting these species. Our approach highlights the value of multi-taxa surveys and the importance of considering temporally dynamic abiotic and biotic drivers for understanding biodiversity patterns and informing conservation management plans.

| Static drivers of predator occurrences
Static variables were important in all habitat models and defined the general ecological niches of various species, and predictions were generally consistent with current knowledge (Dwyer et al., 2016(Dwyer et al., , 2020Francis, 2016;Gostischa et al., 2021). For example, bronze whaler sharks and bottlenose dolphins are classified as coastal species, reflected by their predicted preference for shallow waters close to shore. By contrast, common dolphins were predicted to occur at depths between 30 and 50 m, reflecting their more offshore distribution. While food is hypothesised to be the primary driver of large marine predators (both predation and fear of predation, Gaynor et al., 2019), there are other factors of importance, such as reproduction, which may drive distribution patterns and may be better explained by physical variables (e.g. Rayment et al., 2015).
Sharks show different habitat preferences with ontogenetic shifts, as highlighted by the contrasting predictions of adult and immature hammerhead sharks where there was a clear preference for pelagic and coastal waters, respectively. Pelagic sharks, including adult hammerheads, were predicted to occur in areas of low seabed disturbance and a low percentage of mud substrate, characteristic of the outer Gulf and their pelagic nature. By contrast, immature hammerhead sharks were predicted to prefer shallow (<20 m) waters in areas of high seabed disturbance and mud, which is characteristic coastal habitat and in between inner Gulf islands (Francis, 2016).
These differences in hammerhead sharks highlight the need to model life stages separately for conservation and management (Kinney & Simpfendorfer, 2009).

F I G U R E 2
Predicted monthly distributions of Bryde's whales, common dolphins, bottlenose dolphins, bronze whaler sharks, pelagic sharks, immature hammerhead sharks and richness of large marine predators in the Hauraki Gulf, New Zealand, for the months (a) January to June and (b) July to December. Probability of occurrences for species or species groups ranges from 0 (blue) to 1 (red) in increments of 0.1 (see legend in a). Richness (last row) ranges from 0 (blue) to 3.5 (red) in increments of 0.25 (see legend in a).

| Dynamic drivers of predator occurrences
The SDMs included temporally dynamic environmental and prey variables, providing evidence of the key biotic drivers of the distribution of top predators. The importance of prey in all predator models provides a causal link between predator occurrence and prey distribution (Wirsing et al., 2020). Month was the most important predictor of distribution patterns for all six species or species groups and allowed us to somewhat account for unobserved (latent) processes. This is a common approach using BRT models and other SDM modelling approaches (e.g. Compton et al., 2012).
While not an environmental or biological covariate, it proxies for a wide range of ecosystem variations and was important in explaining the occurrences of both sharks and cetaceans. These patterns reflect increased productivity and availability of prey associated with seasonal changes in climatic and oceanographic processes. Winter and spring see an upwelling of nutrient-rich waters (Sharples & Greig, 1998;Zeldis et al., 2004). The combination of wind-mixing and upwelling produces some of the highest spring chlorophyll-a stocks on the New Zealand continental shelf (Murphy et al., 2001), with a significant bottom-up effect on higher trophic levels (Jillett, 1971;Sharples & Greig, 1998).
Predictions of all sharks were highest from November-December through to March, and as the water cools in May-June for bronze whaler and pelagic sharks. This sharp increase aligns with the response of sharks to warmer SST and low chlorophyll-a concentrations, as found in Australia (Birkmanis et al., 2020). In summer warm waters (Shears & Bowen, 2017) converge with cold nutrientrich waters from the continental shelf and bring a new assemblage of zooplankton and fishes, increasing additional prey sources to large predators .
Immature hammerhead sharks were predicted to occur in extreme coastal habitats, especially in the southern reaches of the inner Gulf in just a few months of summer, which may experience greater fluctuations in salinity, temperature and purity due to freshwater runoff and pollutants (Stevens et al., 2019;Zeldis et al., 2004). Mangroves are nursery habitat in the Gulf, but these have been degraded due to removal for coastal development, which highlights the need to consider pressures on land and in the ocean in the management of coastal marine predators (Hauraki Gulf Forum, 2020;Sievers et al., 2019). In addition, seagrass meadows were once extensive throughout the shallow harbour regions (Morrison et al., 2014), which may have been critical nursery habitat for juvenile bronze whaler sharks, as found in southern Australia (Drew et al., 2019). Therefore, the conservation and management of coastal shark species will need to consider the restoration and protection of critical nursery and breeding habitats in the Gulf.
By contrast, pelagic sharks were predicted to occur in the outer Gulf mainly, in the same region as major fisheries including trawling and set longlines. The smooth hammerhead sharks are prone to the effects of overfishing (Rigby et al., 2019) and have a high postrelease mortality rate due to stress (e.g. Gallagher et al., 2014). Smooth hammerhead sharks are not a target species under the New Zealand Fisheries Act 1996, yet it is often reported as bycatch (Francis, 2016 Gladics et al., 2014). The ecosystem structure and function in the Gulf have changed as a result of human arrival (MacDiarmid et al., 2016). The Gulf ecosystem has been highly modified with a reduction in the historical biomass of seabirds mainly due to introduced predators, and a reduction in marine mammals and fish species due to fishing, alongside the increased phytoplankton productivity with land-borne nutrients (Pinkerton et al., 2015). This reduction in biomass and biodiversity almost certainly affects its food-web dynamics although the nature of the coupling between the lower/middle and upper levels is not well understood (Pinkerton et al., 2015). These centurieslong anthropogenic pressures highlight the need to understand contemporary trophic interactions of large marine predators to identify key prey resources, how sensitive they are to changing resource availability and their resilience to ongoing change (Colbert, 2019).

| Management considerations
Area-based management tools, such as MPAs, aim to protect biodiversity by regulating ocean activities. However, MPAs are predominately static, despite marine habitats and species' distributions often being highly dynamic Kavanaugh et al., 2016). While monthly distribution estimates between our study taxa showed some commonalities, there were also apparent differences suggesting that a generic management plan for coastal predators would be inadequate for protecting these species. Dynamic management is a potential alternative to staticbased protection measures by allowing management decisions to be updated in response to changing environmental or socioeconomic conditions (Lewison et al., 2015;Maxwell et al., 2015).
Accounting for dynamic changes is particularly important for protecting top predators because static management strategies typically offer partial protection to these highly mobile species or require vast areas to be protected to cover their distributions (Dunn et al., 2016). The (relatively) fine-scale spatial-temporal species distributions presented here lend themselves well to informing dynamic management measures (Foley et al., 2010) but come with the caveat that there is increasing inter-annual variation that must be considered in the future. Dynamic habitat management may offer more flexibility than static approaches for managing large marine predators, although there are likely to be conflicts between economic, socio-cultural and conservation objectives, which will need attention to make it work (Gaines et al., 2010). Despite these difficulties, dynamic MPAs may be an effective way of protecting these widely distributed marine predators without the need to designate large areas as protected (i.e. in this case, most of our study area) (e.g. Hausner et al., 2021). In addition, the monthly spatially explicit maps of uncertainty (which also varied over time) provide important additional information required for spatial management (Azzellino et al., 2017;Moilanen & Wintle, 2006). That is, uncertainty layers can allow trade-offs between biological quality and the certainty of that information in spatial planning tools (e.g, see Moilanen & Wintle, 2006;Stephenson et al., 2021).
Sharks and cetaceans use the entirety of the Gulf at all times of the year and their distributions are dynamic both in time and space.
Nevertheless, there appears to be a period of around six months where the occurrences of cetacean and shark species peak, and it is strongly recommended that conservation planning focuses on implementing temporal restrictions of harmful activities within spatially defined boundaries to maximize protection. Marine predators are equally as difficult to manage as they are to study; however, there are excellent examples where conservation goals are balanced with economic objectives and a diverse community of large predators can thrive (Chilvers et al., 2005;Handley et al., 2020).
Finally, any management measures should also consider exploring the spatial overlap between human pressures (e.g. fishing) and habitat preferences of predators to ensure appropriate measures are implemented given the substantial decline in marine biodiversity in the Gulf since human arrival (MacDiarmid et al., 2016). The shifting baseline is a common phenomenon in studies of marine ecosystems, and the Gulf has considerable anthropogenic pressures causing a decline in system function that requires urgent attention (Hauraki Gulf Forum, 2020). Our research has provided a broader focus on multiple large marine predators and their environment, including prey, and future studies should focus on multispecies approaches to better understand marine communities and ecosystem function.

| CON CLUS ION
Few distribution studies assess a wide range of large predators simultaneously, instead of focussing on single families, such as sharks, cetaceans or seabirds, that are easier to access or of high conservation concern. Yet, large predators can be sentinels of the environment because they integrate many different trophic and environmental layers of data that are reflected in their distribution patterns . Assessing and combining the distribution of several top predators may therefore provide information on the wider ecosystem, the impact of anthropogenic stressors and possibly on the distribution of resources, and critical habitat for species not monitored directly . This multispecies data therefore favour an ecosystem-based management approach, where the goal is to protect biodiversity and important habitats rather than individual species (Hyrenbach et al., 2000). Here, we provide fine-scale monthly predictions for a range of top predators simultaneously. The dynamic relationships highlighted between environmental and biotic conditions in this study are essential for understanding large marine predator distribution patterns and possibly the wider ecosystem. Our multi-species approach may result in more effective management if environmental managers and decision-makers dynamically prioritize the protection of these critical ecosystem components.

ACK N O WLE D G E M ENTS
We thank the survey crew: Our extraordinary pilot, William

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors have no conflict of interest to declare.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/api/gatew ay/wos/peer-revie w/10.1111/ ddi.13705.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in Dryad https://doi.org/10.5061/dryad.wm37p vmrn.