The study area extends over 4.2 million km² of the South Pacific Ocean within the New Zealand EEZ (≈ 25–57°S; 162°E–172°W; Figure 1). New Zealand’s EEZ contains highly productive zones of mixing between warm, higher salinity, nutrient poor, northern waters, and cold, lower salinity, nutrient rich, southern waters. These productive zones support high biological diversity and a variety of customary, recreational, and commercial fisheries (Bradford-Grieve et al., 2006; Leathwick et al., 2006; Stephenson et al., 2018, 2020a). The New Zealand continental shelf environment is influenced by a combination of climatic features that produce a range of nearshore conditions, upwelling, tidal mixing, and upper-ocean mixing. Surface and thermocline waters of different origin are separated by three major fronts, the Subtropical Front, the Subantarctic Front, and the Polar Front, while the flow of bottom water is predominately carried by the Deep Western Boundary Current (Chiswell et al., 2015). New Zealand bathymetry is equally as complex, reaching abyssal depths beyond 10,000 m in the Kermadec Trench and includes over 800 distinct sea features (Rowden et al., 2005; Linley et al., 2017).

Basking shark records (n = 401) were collated from various sources including records from observed and self-reported captures in commercial fisheries, public sightings, media reports, museum records, fishery surveys, aerial surveys and beach cast specimens. Records variously included information on date, number of individuals, geographic co-ordinates and source (where available) and spanned a period of 131 years (1889–2020). The only directed sampling effort for New Zealand basking sharks were aerial surveys off the east coast of the South Island over the summer months (January to March) of 2010–2011. No sharks were observed during these surveys. The data were groomed to ensure accurate identification and only records that were confirmed or probable basking shark observations within the New Zealand EEZ were retained. Where geographic co-ordinates were missing an approximate position was assigned based upon the description of the location. This was necessary for all aerial survey records, media reports and public sightings. Because of difficulties in correcting for differences in sampling methods, all catch records were converted into presence records (Elith et al., 2011; Stephenson et al., 2018). To minimize the effect of spatial bias in the occurrence data, species records were aggregated spatially to a 1 km² grid resolution (Aiello-Lammens et al., 2015; Stephenson et al., 2020b). Strandings and records without an approximate date reference (month) were removed. The final dataset included presence records of basking sharks at 369 unique sampling locations (Figure 1).

To characterize variability in the New Zealand marine environment, spatial environmental and biotic variables were collated at a 1 km² grid resolution, with each spanning the breadth of the New Zealand EEZ (Table 1 and Supplementary Table 1, further details are available in Stephenson et al., 2020a). In addition to these variables, spatial estimates of various zooplankton densities (inferred prey for basking shark) (Pinkerton et al., 2020) were used as biological predictors in the models (Supplementary Table 1). Estimates of zooplankton densities did not cover the entire New Zealand EEZ (Supplementary Figures 1.9–1.15). Areas lacking this information will simply represent the modeled relationship between basking shark records and the environmental variables. A preliminary examination of currently available zooplankton density estimates reveals these are likely to cover core areas of basking shark distribution.

Of the available environmental and biotic variables, a subset was selected to be used in the models (Table 1) based on model tuning described in section “Predictor Variable Selection.” Although most of the chosen environmental variables were static (e.g., bathymetry, Bathy), several variables were dynamic in time, representing mean monthly statistics for the past 20 years (e.g., chlorophyll-a concentration, Chl-a, “temporal resolution” column in Table 1). The environmental data spans a much shorter timeframe than the basking shark observations (131 years), and thus, long-term trends in habitat suitability could not be examined. Prior to fitting of the habitat suitability models, values for each environmental and biotic variable were extracted for locations of basking shark records by overlaying the records onto each of the environmental and biotic variable layers using the “raster” package in R (Hijmans and van Etten, 2012). For dynamic environmental variables (mean monthly statistics), recorded dates of basking shark records were used to extract respective values from the month the record was made.

Habitat Suitability Models (HSMs) were used to analyze and spatially predict the distribution of basking shark habitat suitability (measured as habitat suitability index—HSI). Acknowledging that our environmental predictors are mean (or mean monthly) averages for the past 20 years, we explored models for two time periods: 1889–2020 (all data, n = 369) and as subset of basking shark occurrence which matched the time frame of our environmental predictors 2000–2020 (n = 123). The relationship between environmental variables, biotic variables and basking shark records was explored using ensemble predictions (Ensemble HSM) from Boosted Regression Tree (BRT) and Random Forest (RF) models. This approach limits dependence on a single model type or structural assumption and may result in a more robust characterization of the predicted spatial variation and uncertainties (Robert et al., 2016). To estimate basking shark distributions, BRT and RF models require locations of both presences (occurrence records) and absences. Here, true absences (i.e., sample locations where no basking sharks were recorded) were not available for opportunistic records such as public sightings, media reports, or museum records. True absences were also unavailable for the non-opportunistic sampling methods (i.e., trawl tows, observer records, scientific surveys), particularly from commercial records which are complicated with the inclusion of multiple gear types and fishing protocols (thus affecting catchability) and issues regarding a lack of reporting of basking shark interactions (Francis, 2017). Therefore, presence only modeling approaches using pseudo-absences (i.e., locations where basking sharks were not recorded within our study area) were necessary.

We produced an ensemble model by taking weighted averages of the predictions from each model type, using methods adapted from Anderson et al. (2016, 2020) and Georgian et al. (2019). This adapted procedure derives a two-part weighting for each component of the ensemble model, taking equal contributions from the overall model performance (AUC value derived from the “evaluation”) and the uncertainty measure (SD) in each cell, as follows:

where AUC_BRT and AUC_RF are the model performance statistics; X_BRT and X_RF are the model predictions; SD_BRT and SD_RF are the bootstrap SDs; and X_ENS and SD_ENS are the weighted ensemble predictions and weighted SDs, respectively, from which maps of predicted species distribution and model uncertainty were produced.

Of the records retained for use in the models, most basking shark records (72%, n = 265) occurred in the spring and summer months (September to February). Most (69%) records came from fishing interactions (trawl = 244, set net = 7, surface longline = 3), followed by public sightings (13%, n = 47), aerial surveys (11%, n = 41), research vessels (4%, n = 14), and alternative capture methods (e.g., harpoon) or unknown sources (4%, n = 13). Since 2000, most records (84%, n = 103) have been from fishing events, with one aerial record and 19 opportunistic sightings. The most recent coastal sighting and coastal fishing interaction (set net) occurred in March 2012 and August 2014, respectively, and the last known report of a school (≥3 individuals) was from April 2013 where seven individuals were captured in one trawling event. Estimated lengths were available for 169 records (42%); 32 records were <5 m (19%), 126 records were 5–10 m (74.5%), and 11 records were >10 m (6.5%) (Supplementary Figure 4). For observations where length was recorded, sex was available for 81 of these records (20% of all basking shark observations). Most sharks (89%, n = 72) were male and nine (11%) were female (Supplementary Figure 5).

AUC and TSS scores using evaluation data were very similar between models, with the RF model performing slightly better than the BRT model (AUC: 0.92 and 0.89; TSS: 0.72 and 0.69, respectively, Table 2). Both indices indicated the models were useful in predicting basking shark occurrence (>0.7). Measures of BRT and RF model performance scores had low variability (measured by the standard deviation of the mean), suggesting the models were performing consistently across bootstrap samples. Model fits between training data and evaluation data were similar, with model fits for the evaluation data slightly lower than the training data (as would be expected). The similarity of these fits provides some indication that the training data were not overfitted in the models.

The relative importance of each predictor and their influence on basking shark habitat suitability were consistent across BRT and RF models (Supplementary Figures 6, 7). The three most important variables in predicting basking shark habitat suitability were vertical flux (POCFlux, with 26.0% influence on the response), slope (Slope, 14.1%), and turbidity (Turbidity, 10.6%) (Figure 2). Bathymetry (Bathy, 9.7%) and BPI broad (BPI broad, 9.6%) were also moderately important variables. There was a strong positive relationship of predicted basking shark HSI with vertical flux, highest in areas where vertical flux was 20 mgC m^–2 d^–1 or greater than what would be expected for the given depth. High HSI was predicted in gently sloping and less complex seafloor topographies with moderate turbidity. Two depth strata had high HSI—nearshore depths and depths between 200 and 550 m. A less clear relationship was observed between HSI and SST and mixed layer depth (MLD), with low HSI occurring between temperatures of 12.5 and 15°C and in areas where the mixed layer depth was approximately 75 m. There was a weak relationship between HSI and copepod (Copepoda) densities, with low HSI occurring with low levels of copepod densities, a peak in HSI at moderate copepod densities (10–20 counts per 5 nautical miles), and a plateau in HSI values at the highest levels of copepod densities (>25 counts per 5 nautical miles). HSI was lowest at moderate levels of chl-a concentration (Chl-a) (0.5–1.0 mg m^–3) and highest at high chl-a concentration (>1.2 mg m^–3).

Areas of predicted high habitat suitability for basking sharks in New Zealand waters occurred along the continental slope, particularly along the 250 m contour along the North and South Islands, Mernoo Bank, Pukaki Rise, Puysegur Ridge, and around New Zealand’s offshore islands (Chatham Islands, Stewart Island, Bounty Islands, and Auckland Islands) (Figures 3, 4). Within the spatial coverage of species occurrence, areas of moderate uncertainty (SD > 0.2) included most offshore waters north of 40°S, the deeper depths (>500 m) of the Hokitika Canyon, northern Chatham Rise, coastal waters off east coast of the South Island (Canterbury Bight), Foveaux Strait and Puysegur Ridge (Figure 5). The northern North Island and features further from the continental shelf, including the eastern half of the Chatham Rise were outside of the estimated spatial coverage (95% KDE) of species occurrence. In addition, moderate—high uncertainty (SD > 0.2) was reported along deep sea features north of New Zealand, including the Kermadec Ridge and Trench, the Lau-Colville Ridge, and the Norfolk Ridge (Figure 5).

Basking shark habitat suitability was largely influenced by variables representing ocean processes. The environmental predictors used in this work were comprehensive and many were dynamic (i.e., monthly means were available). Overall, areas with high levels of vertical flux of particulate organic matter at the seabed had high habitat suitability. This is likely indicative of higher levels of primary production in the surface ocean and higher prey density in the mesopelagic zone and at the seafloor and may be a suitable proxy when prey data is unavailable. The biotic predictive layers included here were found to have lower influence on habitat suitability compared to some of the environmental predictors. However, prey availability is highly patchy and temporally variable; thus, it is likely that a static variable reflecting prey abundance was unable to accurately represent the spatial distribution of prey. The inclusion of biotic predictors in the model is important in understanding species’ relationship with the marine environment in unobserved space and has been identified as a potential link in understanding effects of climate change. In the Northeast Atlantic, basking sharks are often observed in shallow, highly productive coastal waters during spring and summer months where they feed on zooplankton blooms (Sims, 2008). Historically, basking sharks were commonly observed in similar environments (e.g., east coast of the South Island) around New Zealand (Francis and Duffy, 2002). Although prey preference for New Zealand sharks is poorly understood, there is some relationship between New Zealand basking shark distribution and copepod abundance, as seen in the North Atlantic (Sims and Merrett, 1997).

The inclusion of dynamic (mean monthly) environmental variables here may allow the models to capture temporal change in patterns of basking shark distribution, including seasonal changes and interannual variability. Despite the availability of dynamic environmental predictors, true temporal changes in distribution are difficult to confirm due to the limited amount of biological data (basking shark observations) available. However, in our results, both inshore and offshore regions were highlighted as areas of high habitat suitability. This is particularly evident in the bimodal effect of the bathymetry predictor, where basking shark habitat suitability was observed to be highest in very shallow depths (<100 m), and again at depths between 200 and 500 m. This result is consistent with previous work where basking sharks have been shown to exhibit seasonal vertical space use in the Northeast Atlantic, with tagged individuals occupying shallow depths (<100 m) in the summer months and depths greater than 1,000 m in late winter/early spring (Doherty et al., 2019).

While bathymetry (and slope) were also found to be important predictors, their effect may be partially influenced by basking shark availability to fisheries (see below). Basking sharks have been shown to dive as deep as 1,264 m and have been regularly documented at depths of 600–1100 m (Francis and Duffy, 2002; Gore et al., 2008; Doherty et al., 2017). The species has also been shown to follow distinct water masses at depth, remaining at depths of 250 m or more for months without coming to the surface (Braun et al., 2018; Dewar et al., 2018). Basking sharks are known for complex diel vertical movements, which are thought to be influenced by shifts in prey availability and oceanography (Sims et al., 2005; Dewar et al., 2018). In well-stratified deep waters, basking sharks exhibit normal diel vertical movements (shallow depths at night, deeper depths during daylight), while sharks occupying inshore, inner-shelf areas near thermal fronts conduct reverse diel vertical movements (shallow depths during the day, deeper depths at night) (Sims et al., 2005). This may explain, at least in part, why more contemporary reports of basking sharks in deepwater fisheries have been made during daylight hours (when sharks are occurring at their preferred deeper depth range) but does not offer insight into the disappearance of inshore observations or why individuals are no longer seen at the surface in the region.

Water temperature had relatively minimal influence on basking shark distribution. Basking sharks appear to have a broad thermal range and are therefore relatively unrestricted by temperature (Sims et al., 2003). They are also known to cross tropical regions by submerging into deeper, colder water (Skomal et al., 2009), and one individual was encountered in tropical waters off Indonesia (Fahmi and White, 2015). While gradual changes in sea temperatures may have minimal effect on basking sharks, ocean heat waves and extreme processes associated with SSTs might be more relevant and are expected to become more prevalent with climate change. By 2100, climate change projections predict SST will increase by 2.5°C, which in turn is predicted to lead to a deepening in the surface mixed layer depth (by 15%), declines in primary production (4.5%) and particle flux (12%); the largest changes in macronutrients are predicted in eastern Chatham Rise and southern Sub-Antarctic waters (Law et al., 2018). Such changes may alter food availability (prey distribution and abundance) for basking sharks resulting in shifts in their distribution.

Predictors found to positively influence basking shark HSI could be further explored to better understand historic and future basking shark distribution. Predictors including chl-a concentration and vertical flux are often used as an index of phytoplankton abundance (primary production) and are strongly linked to primary consumers such as copepods. In recent decades, dramatic shifts in chl-a concentration have been reported in the South Pacific and the Southern Oceans (Del Castillo et al., 2019). Significant declines in chl-a concentration were observed in spring and summer months in the South Pacific from 1979 to 2000, while significant increases were linked to extreme summer marine heatwaves in the Southern Ocean between 2002 and 2018 (Gregg and Conkright, 2002; Montie et al., 2020). Highly variable oceanographic conditions around New Zealand are strongly correlated with climate phenomenon such as El Niño–Southern Oscillation (ENSO) (Bradford-Grieve et al., 2006). In the Northern Hemisphere, basking shark movement patterns have been linked to shifts in prey availability and oceanography (Sims et al., 2005; Gore et al., 2008; Dewar et al., 2018). One tagged individual was shown to remain in an area with putative upwelling and high abundance of phytoplankton in the Western Atlantic for up to a month (Gore et al., 2008). In the Northeast Atlantic, a northward shift in basking shark distribution in response to long-term zooplankton declines was found to correspond with declines in basking shark catch in Irish fisheries to the south from 1948 to 1975 (Sims and Reid, 2002). Similar models used in this project could be explored to predict basking shark distribution response to future climate change forecasting and events such as ENSOs [e.g., Earth System Models, Anderson et al. (2020)].

Areas of high basking shark habitat suitability included the east and west coasts of the South Island, Puysegur Ridge, and the Auckland Island slope. Some areas of Chatham Rise, specifically around Mernoo Bank (including Mernoo Saddle) and off the southern slope of Pitt Island (Chatham Islands), were also identified as areas of high habitat suitability. Much of Chatham Rise, however, was outside the spatial coverage of species occurrence and thus habitat suitability predictions hold a higher degree of uncertainty. Chatham Rise is a known hotspot for chondrichthyan diversity in New Zealand waters (Wetherbee, 2000), but interestingly, basking sharks have very rarely been reported from the area. Chatham Rise, as well as Puysegur Ridge, have relatively low densities of copepods (see Supplementary Figure 1.9) and may not be optimal feeding grounds for basking sharks. However, in international waters east and north-east of Chatham Rise, 15 juvenile basking sharks (180–310 cm total length) were reported by Japanese drift net vessels operating at shallow depths (10 m) (Yatsu, 1995). This report suggests juvenile sharks may inhabit epipelagic waters in the open ocean (Francis, 2017).

Given the long temporal span of the data, model predictions may be more representative of past, rather than current habitat suitability, particularly some inshore parts of the predicted distribution. Basking sharks are occasionally recorded from northern New Zealand and were reported to be regular visitors to the Hauraki Gulf during spring in the late nineteenth century (Cheeseman, 1891). The predictions presented here are smoothed over time, as there is a mismatch between the availability of basking shark records (131 years) and environmental data (approximately 20 years). It is possible that the models have highlighted seasonal patterns of distribution by indicating both inshore and offshore regions as areas of high habitat suitability with the presence of bimodality in the bathymetry HSI. However, these patterns could also be indicative of a shift to deeper, offshore habitats in recent decades.

There were several areas where the spatially explicit uncertainty (measured as the SD) was relatively high, indicating the relationship between basking sharks and the environment was more uncertain. Our understanding of basking shark use of the pelagic habitat remains relatively unknown, largely due to the spatial bias in observations (e.g., lack of open ocean pelagic research surveys). In areas with high uncertainty, such as Cook Strait, the northern Chatham Rise, and Foveaux Strait, few basking shark sightings where available and uncertainty might be linked to low sample size. Uncertainties regarding the most northern predictions of habitat suitability (north of 40°S) may, in part, be explained instead by a lack of information on copepod density north of 40°S (Pinkerton et al., 2020), but is more likely associated with a lack of basking shark data. Despite a lack of surface sightings around northern New Zealand, there is some evidence of basking shark captures in deepwater trawl fisheries off central Hawke’s Bay. Basking sharks may be present along the upper slope around northern New Zealand but there is no fisheries-independent data (research or public sightings) to indicate their presence. In the Kermadec Islands where there is minimal human presence, there is one unconfirmed stranding of a basking shark from Raoul Island (Morton, 1957).

The estimate of spatial coverage of species occurrence (top 95% of the KDE of basking shark occurrences) provides a representation of the likely geographic (and in turn environmental) space occupied by basking sharks within New Zealand waters. Predicted distribution outside of this area, should be treated with caution as this represents prediction into largely unsampled space. In this study, the environmental threshold reflects the distribution of presences only—and thus retains any spatial biases associated with these datasets. In particular, the spatial distribution of presences is related to the distribution of fishing effort and human population centers (for opportunistic sightings) and may not be an accurate representation of hotspots. However, using the top 95% of the KDE of basking shark occurrences provides a more conservative estimate of this species’ spatial distribution, which can be useful in determining when modeled predictions are occurring outside of sampled environmental space. This measure provides a meaningful threshold with which to classify broad areas as “uncertain.”

The lack of basking shark records in New Zealand waters during recent years highlights the need to better understand their overall abundance, distribution, movements, and habitat use in the South Pacific Ocean as well as changes to the sighting effort. Currently data collection on the species in the Southern Hemisphere is heavily reliant on interactions with fisheries, especially trawl fisheries (Francis and Smith, 2010; Hernández et al., 2010; Lucifora et al., 2015; Francis, 2017). Because of this, our knowledge of New Zealand basking shark distribution is essentially limited to areas of relatively high historic and current trawl fishing effort (Clarke et al., 2017; Baird and Wood, 2018). Most basking shark interactions occur during the spring-summer months, corresponding to when fishing vessels target commercially important species, such as spawning aggregations of arrow squid (Nototodarus sloanii) (Hurst et al., 2012). As a protected species, it is mandatory to report basking shark interactions with fisheries. However, there is uncertainty in the levels of reporting, and observer coverage is relatively low in some fisheries (e.g., inshore fisheries), particularly around northern New Zealand, so that presence records are likely underestimated (Francis, 2017). Understanding habitat use will assist in assessing risk to fishing activities and could be incorporated into management frameworks such as the spatially explicit risk assessment that New Zealand has in place for other protected species (Large et al., 2019).

Identifying areas of high habitat suitability could also assist in decision making processes for future research efforts. Previous research has identified the need to tag free-swimming basking sharks to better understand species movement, habitat use, and interactions with fisheries (Francis, 2017). This will require the ability to find individuals at the surface, and at an accessible location. Meeting these requirements has prevented any efforts of tagging New Zealand basking sharks to date. In the Atlantic Ocean, basking sharks have been successfully tagged off Cape Cod, Massachusetts and the west coast of Scotland and Isle of Man (Skomal et al., 2009; Doherty et al., 2017; Hawkes et al., 2020) and in the Pacific, off San Diego and Monterey Bay, California (Dewar et al., 2018). By identifying areas of high habitat suitability, research efforts can be directed to specific areas of interest to increase the tagging success. For example, the Auckland Islands has been identified as an area of high habitat suitability for basking sharks where basking sharks were historically sighted at the surface (Parrott, 1958). This area is also a known hotspot for other large filter-feeding vertebrates that feed along the Subtropical Front (STF), a continuous feature within the Southern Tropical Convergence at latitudes 39°–42°S, characterized by elevated primary productivity (Murphy et al., 2001; Rayment et al., 2015; Mackay et al., 2020). Southern Ocean oceanographic fronts have been identified as important foraging areas for a range of marine predators (Bost et al., 2009) and may also be important for basking sharks.

Differences in habitat suitability among sexes or size classes, a common observation among sharks, were not examined at this time due to the relatively small sample size of basking sharks across the region and low availability of size and sex data for most records. This information is becoming more readily available through fisheries observer data collection and should be explored further in the future. New Zealand has a considerable number of records of small (<5 m) sharks, accounting for nearly 20% of records where length has been recorded. Some of the most recent reports of basking sharks in New Zealand waters have been juvenile individuals, including a 3 m female captured east of the Auckland Islands in January 2018 and a 3.3 m male captured off the west coast South Island in August 2020. Both individuals were released alive. Continued collection of biological data on basking sharks is essential for understanding differences in habitat use across life history stages, particularly for juvenile basking shark as they are globally rare and their habitat preference is unknown.

More data on at-sea distribution of basking sharks is required to understand habitat use, threat overlap, and population status throughout the New Zealand and South Pacific region. The total South Pacific basking shark population size is unlikely to be high; in the Northeast Atlantic, basking shark numbers likely do not exceed 10 000 individuals (Lieber et al., 2020). With little population differentiation across global regions, it is plausible that basking sharks observed in New Zealand waters engage in large scale intra- and inter-oceanic migrations (Lieber et al., 2020). This may make species’ detection more difficult in the vast marine space of New Zealand’s EEZ. In the past, aerial surveys have been successful in detecting basking sharks distributed over relatively large spatial scales in New Zealand, and elsewhere for estimating regional population sizes (Francis and Duffy, 2002; Westgate et al., 2014). However, surveys conducted in 2010–2011 off Banks Peninsula failed to locate any basking sharks (Chapman and Duffy, 2011). These surveys have not been repeated since. Basking sharks may travel or feed in subsurface habitat, and therefore go undetected in aerial surveys. Alternative means of tracking these animals, such as autonomous underwater vehicles (AUV) (Hawkes et al., 2020) offer great potential for the future.