On the rise: Climate change in New Zealand will cause sperm and blue whales to seek higher latitudes

Climate impacts affect marine ecosystems worldwide with island nations such as New Zealand being extremely vulnerable because of their socio-economic and cultural dependence on the marine and costal environment. Cetaceans are ideal indicator species of ecosystem change and ocean health given their extended life span and cosmopolitan distribution, but limited data availability prevents anticipating change in distribution under future climate changes. We projected the range shifts of a key odontocete and mysticete species ( Physeter macrocephalus and Balaenoptera musculus ) in 2100 relative to present day in New Zealand waters, using an ensemble modelling approach, under three climate change scenarios of different severity. The results show a latitudinal shift in suitable habitat for both whale species, increasing in magnitude with severity of sea surface temperature warming. The most severe climate change scenario tested generated 61% and 42% loss and decrease of currently suitable habitat for sperm and blue whales, respectively, mostly in New Zealand ’ s northern waters. These predicted changes will have a strong impact on the ecosystem functioning and services in New Zealand ’ s northern waters but also in coastal areas (critical for the species ’ foraging and sur-vival). Not only do these simulated range shifts help to identify future potential climate refugia to mitigate a global warming, they also generate a range of socioeconomic consequences for island nations relying on wildlife tourism, industry, and environmental protection.


Introduction
Humans have been causing severe changes in marine systems, including ocean surface warming (Bindoff et al., 2007), water acidification (Doney et al., 2009), sea level rise (Chen et al., 2017), and variation in marine primary productivity (Hammond et al., 2020). These changes affect marine organisms and communities to varying degrees, causing adjustments in dispersal and distribution patterns, species interactions, and trophic dynamics, among others (reviewed in Doney et al., 2012;Gulland et al., 2022). Based on climate change projections, these alterations and their consequences are anticipated to be even more severe in the near future (Pachauri et al., 2014). Yet, the ecological drivers controlling the responses of marine species to these upcoming climate conditions remain poorly understood.
Not only will climate impacts be more pronounced near the poles and in the tropics (Albouy et al., 2020;Solomon et al., 2007), but islands ( e. g., Pacific Islands) have become increasingly vulnerable to shifts in temperature, rainfall and sea level rises (Veron et al., 2019). As an island nation, New Zealand relies heavily on its marine and coastal environments which, besides their ecological significance, hold crucial economic, cultural, and social values (Davies et al., 2018;MacDiarmid et al., 2013). During the past century, the country has already recorded an increased frequency of marine heatwaves, along with an average air temperature increase of 1 • C and a sea level rise of ~20 cm (MfE and Stats NZ, 2017). Given the current rate of environmental change (Babcock et al., 2019;Cheng et al., 2020), monitoring ecosystems is critical, especially since the marine realm experiences faster changes than terrestrial systems (Poloczanska et al., 2013). As direct comprehensive monitoring of marine ecosystems is challenging and costly (Hazen et al., 2019), alternative strategies use indicator species (i.e., species whose status reflects the condition of their environment) and sentinel species (i.e., species reflecting changes in ecosystem and ocean health) (Hazen et al., 2019;Siddig et al., 2016). Sperm (Physeter macrocephalus) and blue (Balaenoptera musculus) whales are potential good indicator and sentinel species because of their extended life span and sensitivity to seasonal environmental shifts in their prey distribution and abundance (Hazen et al., 2019;Moore, 2008;Silber et al., 2017). Therefore, identifying how environmental shifts will shape the distribution of these indicator and sentinel species can serve as an early warning system to anticipate current or potential ecosystem changes.
However, until the cease of commercial whaling in New Zealand, both sperm and blue whales were heavily targeted by whalers and were, partly due to their low reproductive rates (Chiquet et al., 2013), subsequently decimated to within < 1% of their original population size (Branch et al., 2004;Branch et al., 2007;Sears and Perrin, 2018;Whitehead, 2018). While sperm whales are currently listed as 'Vulnerable' and blue whales as 'Endangered' by the International Union for Conservation of Nature (IUCN) Red List of Threatened Species (Cooke, 2018;Taylor et al., 2019), New Zealand's national threat classification system lists both species as 'data deficient'  since no abundance estimates exist for the local populations. This data deficiency limits current understanding of these species' habitat preferences and how they might respond to shifts in climate conditions, making predictions about future alterations of the marine system and potential socio-economic implications challenging.
Species distribution models (SDMs) are the primary tools used to forecast the impacts of climate change on biodiversity (Dawson et al., 2011) and to guide climate adaptation strategies (Guisan et al., 2013). Correlative species distribution models measure species exposure to climate change by relating observations of species' occurrence or abundance with aspects of climate using statistical or structural models (Ancillotto et al., 2020;Elith and Leathwick, 2009;Hazen et al., 2013;Ramirez-Leon et al., 2021;Randin et al., 2020;Williams et al., 2009). Due to simple and flexible data requirements, correlative SDMs provide practical advantages over alternative, more mechanistic approaches (which typically require detailed demographic and physiological data, Dormann et al., 2012), making them the principal method for predicting species distributions and community structure in space and time (Guisan et al., 2013).
Here, we investigate how environmental changes shape the presentday distribution of two indicator/sentinel species in New Zealand, and will constrain the species' respective ranges in the upcoming decades depending on the severity of the ongoing global warming. More specifically, we use an ensemble of nine of these correlative SDMs to (1) assess the present-day distribution of blue and sperm whales in New Zealand waters, and (2) model the shift in distribution of blue and sperm whales between the present day and the future (i.e., 2100) under three International Panel on Climate Change's (IPCC) Representative Concentration Pathways (RCP) of low (RCP 2.6), medium (RCP 4.5), and high (RCP 8.5) severity. We test the following hypotheses: 1) both species will show a general poleward shift in distribution which will depend on the severity of the climate change scenario; 2) given their taxonomic, trophic, and wider ecological differences, odontocetes (toothed whales) and mysticetes (baleen whales), represented by sperm and blue whales, will respond differently to climatic changes on a regional level; 3) the changes in species distribution will be strongly driven by increasing sea surface temperature. By projecting our models to 2100, we capture the largest divergence between the three climate scenarios tested (Pachauri et al., 2014) and ensure that we obtain conservative estimates of potential climate refugia and joined distribution shifts.

Study area and species
Our study area ranges from 158 to 190 • E and 55-28 • S and includes the Tasman Sea, and parts of the South Pacific and the Southern Ocean ( Fig. 1 & Fig. S1). The oceanic floor is complex, including several basins (e.g., Tasman Basin and the New Caledonia Basin), trenches (including the Kermadec Trench) and rises (such as the Chatham Rise, the Chatham Plateau and the Challenger Plateau) (Heath, 1985), see Fig. S1 for names and locations. While the extended continental shelf reaches beyond 200 nautical miles from the coastline off Kaikōura, South Island ( Figure S1), the Kaikōura Canyon lies just 500 m off the coastline with depths of more than two kilometres at the edge of the continental shelf (Childerhouse et al., 1995).
Both sperm whale and blue whale have a cosmopolitan distribution ranging across different climatic zones including the ice edge of both hemispheres to the tropical waters near the equator (Sears and Perrin, 2018;Whitehead, 2018). As the largest odontocetes and deep-diving mesopelagic apex predators, sperm whales primary feed on squid, but also demersal fish in some regions, such as in New Zealand (Gaskin and Cawthorn, 1967;Gómez-Villota, 2007;Okutani and Nemoto, 1964;Santos et al., 1999). A known important foraging location in New Zealand water are the deep nearshore waters off Kaikōura, where sperm whales occur year-round (Childerhouse et al., 1995). Blue whales (Balaenoptera musculus) are mysticetes that feed on zooplankton, particularly on dense krill schools usually occurring in cold-water coastal or offshore upwelling systems (Croll et al., 2005;Croll et al., 1998;Gill, 2002). Once considered infrequent visitors to New Zealand (Branch et al., 2007;Olson et al., 2015), our understanding of blue whale occurrences around the island nation has evolved recently, in particular with the South Taranaki Bight (Fig. S1) being identified as an important blue whale foraging habitat, with dense aggregations of Australian krill (Nyctiphanes australis) (Torres, 2013). Analyses of photoidentification catalogues suggest a high degree of isolation in New Zealand (Barlow et al., 2018), and while the possibility of migration between New Zealand and eastern Australia and Bass Straight is still unclear, the idea that at least some animals may be resident to New Zealand waters for at least part of the year gains support (Barlow et al., 2018;Goetz et al., 2021;Torres, 2013).

Distribution data
We compiled sighting records of sperm and blue whales in New Zealand waters from both the literature (Peters and Stockin, 2022;Stephenson et al., 2020), and the sightings database of the New Zealand Department of Conservation (DOC), ranging from 1977 to 2020. Sighting records originate from various sources, including scientific and seismic surveys, the public, fishing vessels, boat charters, ferries, and various aircraft. Out of a total of 807 (sperm whales) and 730 (blue whales) records, we removed duplicates (within and between databases), stranding records/records located on land, and sightings outside our study area, resulting in a total of 579 (sperm whales) and 477 (blue whales) data points (Fig. 1). The data showed for both species a similar clustered spatial structure (i.e., Clark and Evans index < 1, see description of the method in Table S2 (Clark and Evans, 1954)) over the entire study period. The species distribution models that we used (see 2.5 Species distribution models) require presence/absence or presenceonly data, but true absences were not available in our dataset. Therefore, we generated pseudo-absence data by randomly sampling points for each species within the study area where the focal species was not recorded ( Barbet-Massin et al., 2012). Depending on the model, performance is generally highest with a large number of pseudo-absences (e.g. 10,000) and/or with a 10:1 ratio of pseudo-absences to presences (Barbet-Massin et al., 2012). Therefore, we generated a total of 7,000 pseudo-absences for each species.
In the Southern Hemisphere, three subspecies of blue whales are currently recognized: the Antarctic blue whale B. m. intermedia and the pygmy blue whale B. m. brevicauda (Rice, 1998), and Chilean blue whales (recognized as a subspecies by the Society for Marine Mammalogy Committee on Taxonomy, but not yet named, Barlow et al., 2018). Recent research has described a genetically distinct blue whale population occurring year-round in New Zealand waters (Barlow et al., 2018), that is genetically most similar to pygmy blue whales. This newly described population has extremely low genetic diversity, which lowers their ability to respond to changing environments and thus makes them vulnerable to future climatic change (Attard et al., 2015;Barlow et al., 2018). However, since the two subspecies that occur in New Zealand are difficult to distinguish based purely on morphology, we use the generic term 'blue whale' and Balaenoptera musculus to refer to both herein.

Current environmental variables
We used eight environmental variables to build our SDMs and predict the habitat suitability of sperm and blue whales based on existing knowledge: depth (m), slope gradient (degrees), slope aspect (compass direction), distance to shore (km), distance to 200 m depth contour (km), distance to 1000 m depth contour (km), sea surface temperature ( • C) and chlorophyll a (surface concentration, in mg/m 3 , as a proxy for phytoplankton biomass in the surface layer and thus primary productivity) (Aïssi et al., 2014;Mackay et al., 2018;Pace et al., 2018;Stephenson et al., 2020). We compiled data for these environmental variables at 0.01 × 0.01 • spatial resolution (latitude × longitude) for our entire study area, as well as for each respective sighting.
We calculated depth as a negative value in metres using the bathymetry information from the marmap package (Pante and Simon-Bouhet, 2013) in R version 4.0.2 (R Core Team, 2022). Using the 'terrain' function in the raster package, we calculated the bathymetric slope gradient (defined as the degree change from one depth value to the next) and slope aspect (defined as the compass direction of the slope, ranging from 0 to 360).
For each sighting record, we calculated the distance 200 and 1000 m depth contour by calculating the shortest Euclidean distance between each data point to the 200 m and 1000 m depth contour lines (based on the bathymetry data). We downloaded monthly sea surface temperature averages for years 2003-2018 from 'https://www.ncei.noaa.go v/erddap/index.html' (Grid DAP data, dataset ID: ncdc_oisst_-v2_avhrr_by_time_zlev_lat_lon) using the rerddap package in R (version 4.0.2). We retrieved monthly chlorophyll a averaged at 4-km resolution for years 2003-2018 from the Moderate Resolution Imaging Spectroradiometer (MODIS) data from NASA's Aqua satellite system from 'htt ps://oceandata.sci.gsfc.nasa.gov/opendap/catalog.xml' using the obpgcrawler package in R (Tupper, 2018).
For both sea surface temperature and chlorophyll a, we averaged monthly values across all years and assigned the respective values to each sighting (i.e., presence data) whereas each pseudo-absence, were associated with monthly sea surface temperature and chlorophyll a values (averaged across years 2003-2018) randomly assigned. We tested for multicollinearity among climate variables using the variance inflation factor (VIF) and ensured that all variables returned VIF <10 to minimise correlation between the variables (Table S1).

Future environmental variables
For the projection to 2100, we selected the climatic variables from three Representative Concentration Pathways (RCPs, Pachauri et al., 2014): RCP 2.6 (high effort to curb greenhouse concentration that sees emissions peak early, then fall due to active removal of atmospheric carbon dioxide), RCP 4.5 (medium effort to curb greenhouse concentration, where radiative forcing level stabilizes by 2100 using a range of technologies and strategies for reducing greenhouse emissions), and RCP 8.5 (worst-case scenario due minimal effort to curb greenhouse concentration with a failure to reduce warming by 2100). As the confidence intervals of the different IPCC scenarios only considerably diverge later in this century, we chose to project habitat suitability to the year 2100, as this is the point in time where we will see the largest divergence between the three scenarios.
Given the nature of our environmental variables, we assumed that only chlorophyll a and sea surface temperature will change within the next 80 years, while depth, slope gradient, slope aspect, distance to shore, distance to 200 m depth contour, and distance to 1000 m depth contour will remain constant. Thus, we used these two climatic variables  (Rickard et al., 2016;Taylor et al., 2012), as these models cover a range of CMIP5 model performances and estimate the changes in both sea surface temperature (Fig. S4) and chlorophyll a (Fig. S5). For each scenario and model, sea surface temperature and chlorophyll a variables are based on long-term average data (2080-2100) and were applied as anomalies relative to present day (averaged for 2003-2018).

Species distribution models
We used correlative species distribution models to predict the habitat suitability of sperm and blue whales by estimating the statistical relationship between the in situ species occurrences (and pseudo-absences) and the environmental variables of those locations. Amongst the broad range of available statistical algorithms to predict species distributions, we used an ensemble modelling approach based on nine widely used algorithms: Artificial Neural Networks ( Table S2). This approach integrates statistical models of different complexities and statistical properties when projecting a species' distribution through time (Araújo and New, 2007;Elith et al., 2011) and ensures that several possible projections are considered as it maps both the main trend (i.e., mean, median, or some other percentile) and the overall variation (and thus uncertainty) across all models.
For each species, we first randomly split our dataset (including pseudo-absences) into a 70% subset that we used to train each of the nine models using present-day environmental variables. To account for the stochasticity in the pseudo-absence generation, we repeated this process 20 times, thus generating 20 different training and evaluation datasets. Using 20 iterations proved to be sufficient to provide a robust estimate of species distributions and subsequent analyses because the variation in prediction as a function of the number of repetitions is very low (Figs. S2 and S3). Then, we computed each of the nine models independently and used the remaining 30% of our original dataset to evaluate each model's performance (i.e., k-fold cross validation, Fielding and Bell, 1997). We evaluated model performance for each repetition using the area under the receiver operating characteristic curve (AUC) and the true skill statistic (TSS), two intuitive metrics to assess the predictive performance of species distribution models transposed into presence-absence mapping (Allouche et al., 2006;Swets, 1988). From the relative suitability map generated by each model for each repetition, we determined a threshold maximizing TSS (which includes both sensitivity and specificity, Guisan et al., 1998) below which we considered the species as 'absent'. This threshold method is commonly used to transform continuous probabilities of suitability into probabilities of presence/absence in species distribution models (Nenzén and Araújo, 2011).
We projected on the complete study site and averaged predictions for each model across the 20 repetitions. We then generated the final 'ensemble' projection averaging the predicted occurrences across all models, while weighting each model contribution to the average based on its respective TSS score (Thuiller et al., 2009), assuming that TSS is more reliable than AUC as a measure of accuracy when using dichotomous presence/absence data (Allouche et al., 2006). Models with higher TSS thus had a greater contribution to the final averaged estimate.
To predict future habitat suitability for each species (i.e., sperm and blue whales) under three climate change scenarios (i.e., RCP 2.6, RCP 4.5, and RCP 8.5), we followed a similar approach and projected each of the nine models independently (with 20 repetitions per model) accounting for the future change in sea surface temperature and chlorophyll a (see '2.4 Future environmental variables') generated by the four General Circulation Models (i.e., HadGEM-ES, IPSL-CM5A-LR, MIROC-ESM-CHEM and MPI-ESM-LR) per climate scenario (i.e., 20 repetitions per species distribution model × 3 climate change scenarios × 4 climate models = 240 simulations). For each climate change scenario simulated by each of the four General Circulation Models, we first averaged each of the nine species distribution models' outputs across the 20 repetitions, and calculated the ensemble predictions using the TSS-weighted average (based on their respective TSS score previously calculated) of the nine species distribution models. Ultimately, we returned and mapped the mean probabilities of species presence along with its standard deviation for each grid cell, for each RCP climate scenario, across the four different generalised circulation models.

Range shift analyses
Based on the outcomes of the species distribution models (i.e., maps of species distributions), we ran the following analyses: (i) quantifying both sperm whale and blue whale range shifts between present day and 2100, (ii) identifying the environmental variables that mainly drive these species' range shifts, and (iii) analysing the main response curves of each species to highlight the relationship between the probability of presence for the species and the main environmental variables driving their shift in distribution.
We identified the changes in habitat suitability for sperm and blue whales in the forecasted scenarios relative to the present day by calculating gained areas (newly suitable habitat in 2100 that were unsuitable during present day), lost areas (suitable habitat in present day that are projected unsuitable in 2100), areas with an increase in probability of presence (higher suitability in 2100 than present day) and areas with a decrease in probability of presence (lower suitability in 2100 than present day). We also identified the common trends in habitat suitability changes and areas of potential importance for both species (i.e., location showing the same shifts such as: gain, loss, increase, or decrease) for both sperm and blue whales under each future climate scenario. Identifying and highlighting areas in which both species exhibit the same changes can be useful for the placement of marine protected areas and the legislation of oil and gas exploration.
We estimated the individual contribution of all variables in the species distribution models (Thuiller et al., 2009). For each of the nine statistical algorithms used, we first used their present-day projection as a benchmark. Subsequently, we ran these algorithms with one environmental variable changed (by randomly reshuffling the values of the variable) while the others remained constant. We then calculated a (Spearman) correlation score between the new prediction and the benchmark prediction to use as a metric of relative variable importance (a high correlation score shows that the randomised variable has little importance for final predictions). We repeated this process until each of the eight environmental variables had been randomised (one at a time) using all 20 training datasets, resulting in 20 iterations per variable. We subsequently calculated the mean and standard deviation variable importance for each variable across the 20 iterations per algorithm, and we finally calculated the ensemble predictions using the TSS-weighted average (based on their respective TSS score previously calculated) of the nine species distribution model algorithms.
We evaluated the responses of the species distributions to the gradients of explanatory variables based on the response curves derived from each model. The response curves were generated following this calculation: N-1 variables were held constant at their mean value whilst the variable of interest contains 100 points varying across the maximum and the minimum of its range. Variation in predictions, made to these 100 cells, only reflects the effects of the one selected variable. Thus, a plot of these predictions allows visualisation of the modelled response to the variable of interest, contingent on the other variables being held constant. This was done subsequently for all the selected variables.
Detailed methods are reported in Table S2 following the protocol outlined in Zurell et al. (2020).

Present-day predictions and variable importance
The ensemble present-day (2003-2018 average) projections for sperm whales showed the highest probability of occurrence in waters east of the country (e.g., Kaikōura Canyon) and in the southwest (Fig. 1a), whereas ensemble blue whales prediction showed its highest (> 0.7) probability of habitat suitability in the South Taranaki Bight and in the Hauraki Gulf (Fig. 1b). These results are supported by a high predictive power of the ensemble models (sperm whale: AUC = 0.95, TSS = 0.77, with only < 2.7% of total sperm whales presences misrepresented by the models in the northwest area; blue whale: AUC = 0.99, TSS = 0.90, Fig. 2 & Fig. S2).
The estimates of variable importance showed that distance to shore, distance to the 1000 m depth contour, depth, and sea surface temperature contributed most to predicting sperm whale presence, while depth, slope, distance to the 200 m depth contour, and sea surface temperature contributed most to predicting blue whale presence (Table 1).
Sperm whales showed highest probability of habitat suitability in waters < 200 m depth, between 50 and 200 km distance to shore, and between 20 and 200 km distance to the 1000 m-depth contour ( Fig. 3a-e). Probability of habitat suitability was high (i.e., > 0.9, Fig. 4h) for sea surface temperature between 7 and 24 • C, but decreased to almost zero outside this range. The probability of habitat suitability for blue whales increased (from 0.25 to 1, Fig. 3) with increasing sea surface temperature and decreasing depth but decreased steeply (i.e., < 0.25, Fig. 3 h) when sea surface temperature was above 21 • C and when depth exceeded ~200 m (Fig. 3a). Blue whales' probability of habitat suitability was also the highest (~1, Fig. 3c) when slope was close to zero and distance to the 200 m-depth contour was less than ~200 km (Fig. 3d).

Future projections for different IPCC scenarios
Variation in predicted habitat suitability among the four different climate models was low for both species (Table S3, Fig. S6). Sperm whales' habitat around the northern offshore islands (Norfolk Islands and Raoul Island) and at the outer edges of their current northern distribution becomes increasingly unsuitable with each scenario (Fig. 4a, c, & e), with a five-fold increase in habitat loss in RCP 8.5 compared to RCP 2.6. Newly suitable habitat remains relatively constant across all three scenarios (9.45 -12.95%), with most new habitat gained around the southern offshore islands (Auckland Islands and Campbell Island, Fig. 4a, c, & e). While in RCP 2.6, increase outweighs decrease in habitat suitability five-fold, this is almost even in RCP 8.5 (43.86% increase vs 43.48% decrease). Most of North Island's waters decrease in habitat suitability, as do large areas surrounding the South Island, including Kaikōura. Increase in sperm whales' habitat suitability is largely restricted to southern and eastern offshore islands (Auckland, Campbell and Chatham Islands).
Blue whales' suitable habitat loss commences at the northern end of its present distribution in RCP 2.6 ( Fig. 4b), where it doubles in area moving south in RCP 4.5 (Fig. 4d). New habitat is gained mainly in RCP scenarios 4.5 and 8.5 (16.31% and 31.88%, respectively), mainly at the western and south-eastern edges of the current distribution, as well as around the Chatham Islands, whereas RCP 8.5 generated a total habitat loss of a quarter of the blue whales' present habitat, mostly surrounding the North Island (Fig. 4f). While increase and decrease are more variable in scenarios RCP 2.6 and RCP 4.5, a clear geographical shift in habitat suitability is visible in RCP 8.5, with decrease all around the North Island, and increase around the South Island and offshore islands. Interestingly, each RCP scenario showed some decrease in blue whales' habitat suitability in the north-eastern part of the South Taranaki Bight (Fig. 4b, d, & f).
The results showing that RCP 8.5 scenario returns the worst prediction with more habitat loss (12.95% for sperm whale, 24.95% for blue whale, compared to the current suitable habitat) than RCP 2.6 are largely consistent with the results for joined responses (i.e., changes in predicted habitat suitability that were the same for both species, Fig. 5). Loss of both species' suitable habitat as well as decrease in projected Fig. 2. Mean true skill statistic (TSS) and mean area under the curve (AUC) for training and test datasets for a) sperm whales (Physeter macrocephalus) and b) blue whales (Balaenoptera musculus). Averaged model runs using training and test datasets obtained AUC scores of > 0.7 and > 0.8 and TSS scores of > 0.4 and > 0.6 for all nine algorithms used for sperm whale and blue whale data, respectively. These results indicate very good predictive power of the fitted models. For both species, RF, BRT and MaxEnt had the best performance. suitability occurred mainly around the northern tip of the North Island, extending southwards with increasing scenario severity. Habitat suitability is projected to increase in coastal and some offshore areas all around New Zealand in RCP 2.6 ( Fig. 5a) and 4.5 (Fig. 5b), and our ensemble models predict habitat gains and increase in suitability (of 89% relative to their total joined changes in RCP 2.6, and 93 % in RCP 4.5, Fig. 5a and b respectively) in coastal and offshore areas surrounding New Zealand as well as its eastern and southern offshore islands ( Fig. 5a and b). However, under RCP 8.5, 40% of joined changes are loss of suitable habitat and decrease habitat suitability, predominantly in coastal and offshore waters surrounding the North Island. The increase in habitat suitability (56% of joined changes) and gain of newly suitable Vertical histograms represent, for each panel, the standardized latitudinal distribution of probability of occurrence for the present day (purple histogram) and each future climate scenario (yellow histogram). We also provide for each panel, the percentage of areas gained, lost, showing an increase and a decrease in probability of presence in 2100 relative to present day (Fig. 1). White areas represent the land mass.
habitat (4% of joined changes) for both species under this most severe scenario is largely limited to coastal waters of the South Island, as well as around eastern and southern offshore islands (Auckland and Chatham Islands, Fig. 5c).

Discussion
The present-day predictions of habitat suitability for both sperm and blue whales generally concur with previous findings (Stephenson et al., 2020), highlighting the high predictive power of our ensemble modelling approach (Fig. 2, Guisan et al., 2017). Our results place blue whales predominantly in areas associated with upwelling such as the South Taranaki Bight (Branch et al., 2007;Chiswell et al., 2017;Rennie et al., 2009;Torres, 2013;Warren et al., 2021), and sperm whales in deeper waters over canyons and trenches (Gannier and Praca, 2007), such as the Kaikōura Canyon. Areas with high habitat suitability for sperm whales (> 0.7, Fig. 1a), besides waters off Kaikōura, were coastal waters surrounding the southern tip of the North Island following the edge of the Hikurangi Trough, and some areas along the south-western edge of the South Island at the end of the Puysegur Trench, indicating the species' preference for productive waters along the shelf edge (López and Methion, 2019;Whitehead, 2018). For blue whales, besides the South Taranaki Bight, areas of high habitat suitability were the Cook Strait, the Hauraki Gulf, and waters off the Karamea Bight (Fig. 1, see Fig. S1 for place names), which are known areas of blue whale abundance (Barlow et al., 2018;Goetz et al., 2021;Torres, 2013).
We showed that sea surface temperature is one of the main predictor variables of habitat suitability for both sperm whales and blue whales in New Zealand waters (Table 1) which echoes the difference in thermal constrains for both species (Fig. 3). These differences in the effect of temperature on habitat suitability between the two species is likely due to their differing ecology and habitat requirements. Both species are cosmopolitan and occur widely throughout both tropical sub-polar waters, but sperm whales are known to have a wider temperature tolerance (Evans, 1997;Pirotta et al., 2011;Whitehead, 2018). Similar to other baleen whales such as humpback whales (Derville et al., 2019), sea surface temperature has been shown to be an important predictor for blue whales in other regions (e.g., Chile, Buchan and Quiñones, 2016;California, Fiedler et al., 1998;southern Australia, Gill et al., 2011;Shabangu et al., 2019), which is likely due to the relationship between sea surface temperature and primary production, as well as upwelled nutrients, which provide optimal conditions for krill (Whitehead et al., 2010). However, chlorophyll a showed a surprisingly minor impact on both species' habitat suitability (Table 1) in particular for blue whales (Stephenson et al., 2020). Being highly mobile predators, cetacean presence and movements are strongly influenced by their prey Palacios et al., 2013). We used chlorophyll a as an indicator for primary productivity, which is a relatively suitable measure of food abundance for blue whales because their main prey (i.e., krill) is generally found in areas of high primary productivity (Barlow et al., 2020;Brinton, 1962;Mauchline, 1969;Torres, 2013). For sperm whales, we used chlorophyll a as a crude measurement for their prey abundance (generally large deep-water organisms such as demersal and mesopelagic fish and cephalopods, Clarke and Roper, 1998;Gaskin and Cawthorn, 1967;Gómez-Villota, 2007;Whitehead, 2018), due to a limited knowledge of their prey's ecology and potential response to climate change (de la Chesnais et al., 2019). However, the temporal lag between primary production and cephalopods, caused by the time it takes for primary productivity to work through the food web (Jaquet, 1996), could explain the low chlorophyll a contribution in sperm whales habitat suitability. Time-lag analyses could potentially capture such ecological processes but are mostly applied at local scales (e.g., New Zealand South Taranaki Bight region) over small temporal resolution (e. g., daily projections for a maximum of three weeks) (Barlow and Torres, 2021), which is challenging to apply across New Zealand waters for the next 80 years (using 20-yrs mean climate data).
The expected changes in sea surface temperature and chlorophyll a (Figs. S4 & S5) lead to southward shifts in habitat suitability for both species by 2100 that will increase with the severity of the climate change and, more importantly, vary regionally. Sperm whales typically showed higher gain and increase in habitat suitability than loss and decrease, because the expected changes in future sea surface temperature remain in the thermal tolerance for the species (Fig. 3 & Fig. S4) except in the far north of the study site where habitat loss becomes predominant. However, in RCP 8.5 (Fig. 4), the decrease in suitability in nearshore waters is of concern, in particular off Kaikōura, since this is an important sperm whale foraging ground year-round (Childerhouse et al., 1995;Giorli and Goetz, 2019;Guerra et al., 2017;Guerra et al., 2020;Jaquet et al., 2000). Foraging is one of the most critical behaviours for survival for any animal, but particularly for large predators with high energy expenditure (Kershaw et al., 2021;Spitz et al., 2012). Similarly, the South Taranaki Bight is an important foraging ground for New Zealand's blue whales (Torres, 2013). While we estimated an increase in suitability in this area by 2100 for RCP 4.5, we also project a decrease in suitability in RCPs 2.6 and 8.5, Fig. 4 d & f) in its northern end range, mostly driven by an increase in sea surface temperature in waters surrounding the North Island (> 21 • C, Fig. S4), which is where the foraging ground is located (Torres, 2013). Models for Pacific populations of five mysticete species, including blue whales, predicted declines and local extinctions due to warming-induced prey reduction and increased interspecific competition (Tulloch et al., 2019), highlighting the vulnerability of these large predators to changes in prey availability. Since there are only few known blue whale foraging grounds outside of the Antarctic (Gill, 2002;Gill et al., 2011;Rennie et al., 2009), a change in habitat suitability in this areas could potentially have population-level consequences for this endangered species.
We acknowledge that by using a static correlative modelling approach, we limit the ecological interpretation of our results since the processes are implicit in the model and we assumed that the species distribution is in equilibrium with its environment and that the speciesenvironment relationships remain constant throughout space and time (Dormann et al., 2012). First, the ensemble modelling approach could be improved by using pseudo-absence differently for each SDM, but the way the prediction of each model is weighed to obtain the final ensemble predictions makes this procedure challenging (Barbet-Massin et al., 2012) and would only result in minor improvements of already strong model performances (AUC > 0.9 for both species, Fig. S2 ). Second, the marine environment being highly dynamic, particularly regarding mobile species (Mannocci et al., 2017), which could potentially be captured better by process-based dynamic SDMs (Fernandez et al., 2017;Reside et al., 2010). Using a process-based approach (including demographic processes, species adaptation, etc.) could also better account for novel or no-analogue species responses that may occur with future conditions (Williams and Jackson, 2007). However, dynamic process-based models require environmental predictor variables with a high spatio-temporal resolution, as well as temporally and spatially unbiased sightings, which are often unavailable (see discussion in El-Gabbas et al., 2021). Finally, we suggest that including species migration in our model would be at this stage irrelevant because of (i) the lack of conclusive understanding regarding the migratory behaviour Mӧller et al., 2020;Szesciorka et al., 2020;Warren et al., 2021), and movement (Whitehead, 2003) of our species for the New Zealand region, and (ii) the mobility of our species, which would be high enough to allow them to track the changes in their habitat suitability . Despite some regional differences between their projected future distributions, both species are expected to shift their range southward which opens a whole range of socioeconomic implications for New Zealand. For example, if sperm whale numbers off Kaikōura should decrease considerably, or change from year-round to more seasonal abundance, the tourism industry would likely be negatively affected due to fewer and less reliable sightings affecting the known direct and indirect employment demographics of the town in respect to whale watching activities (Fumagalli et al., 2021). Similarly, while no whale watch operations currently exist for blue whales in New Zealand owing to their relative inaccessibility to tour vessels, the current and projected habitat of blue whales remains of significant interest and accessibility to industry. With the Taranaki Bight at the centre of significant pressure for seabed mining, the international importance of the region as a foraging region for New Zealand blue whales (Torres, 2013) has not prevented the mining application to be reconsidered by the Environmental Protection Agency. Should future consent allow seabed mining in the southern regions, predicted here to have increasing importance for foraging blue whales because of climate change, this may have severe consequences on the future size of the blue whale populations.
Given the stark decline in whale populations due to commercial whaling, many ecosystem services have been severely degraded already. Great whales (baleen whales and sperm whales) are also known to be marine ecosystem engineers, as they facilitate nutrients transfer from deep waters to the surface, as well as across latitudes via migration from feeding to calving areas (Roman et al., 2014). These whales further sequester large amount of carbon to the deep sea, thus contributing to natural climate-change mitigation (Lavery et al., 2010;Martin et al., 2021;Roman and McCarthy, 2010) and maintain iron levels in the surface layer, which promotes krill abundance (Nicol et al., 2010;Pershing et al., 2010;Roman et al., 2014;Roman and McCarthy, 2010;Willis, 2014). Our results indicate that, in particular at local and regional scales, a further reduction in these services due to whales moving to other, more suitable areas, could affect wider ecosystem functioning and destabilize ecological processes (Roman et al., 2014).
However, the areas showing increases in habitat suitability and newly suitable habitat for both species in the South Island and offshore islands could also potentially be seen to serve as climate refugia (Ban et al., 2016;Verdura et al., 2021) for both species, which might warrant their increased protection in the future. While this study focusses on sperm and blue whales only, it is possible that distribution models with future projections for other species will yield similar results, particularly for cetaceans with similar habitat requirements as sperm whales (e.g., beaked whales, Visser et al., 2021) and blue whales (e.g., sei and Bryde's whales, Kato and Perrin, 2018;Romagosa et al., 2020;Silva et al., 2019). Sperm and blue whales are cosmopolitan species with a global distribution. It is likely that their response to climatic changes presented here, will follow a similar pattern in other places. Given that increases in sea surface temperature are likely to be amplified in the northern hemisphere (Albouy et al., 2020), whale populations may be even more affected by climate change.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.