Climate-driven global redistribution of an ocean giant predicts increased threat from shipping

Abstract

Climate change is shifting animal distributions. However, the extent to which future global habitats of threatened marine megafauna will overlap existing human threats remains unresolved. Here we use global climate models and habitat suitability estimated from long-term satellite-tracking data of the world's largest fish, the whale shark, to show that redistributions of present-day habitats are projected to increase the species' co-occurrence with global shipping. Our model projects core habitat area losses of >50% within some national waters by 2100, with geographic shifts of over 1,000 km (∼12 km yr^-1). Greater habitat suitability is predicted in current range-edge areas, increasing the co-occurrence of sharks with large ships. This future increase was ∼15,000 times greater under high emissions compared with a sustainable development scenario. Results demonstrate that climate-induced global species redistributions that increase exposure to direct sources of mortality are possible, emphasizing the need for quantitative climate-threat predictions in conservation assessments of endangered marine megafauna.

Keywords: Climate-change ecology; Climate-change impacts; Conservation biology.

Fig. 1. Habitat suitability for whale sharks under current and projected environmental conditions.

a,b, Regions of high (yellow) and low (blue) habitat suitability are indicated for the north Atlantic (NA), east Indian Ocean (EIO) and east Pacific (EP) based on current climatologies (2005–2019) (a) and their sum weighted latitudinal density distributions coloured by decade and scenario (b). c,d, Regions of increase (red), decrease (blue) and no change (white) are indicated for NA, EIO and EP based on 2086–2095 ssp585 climatologies (c) and their latitudinal density distributions for cells containing positive (>0.5, red) or negative (<−0.5, blue) change values separated by decade and scenario (d). Each map shows outputs from GAMs built from tracking data from the respective region projected at the ocean basin scale, and the current IUCN distribution limits are displayed in each map as greyed-out boundaries. Mapped results for other regions are given in Supplementary Figs. 1–7.

Fig. 2. Change in habitat area and quality within boundaries.

a, Shifts in mean habitat suitability within LMEs ordered from low (left) to high (right) current habitat suitability, with small grey points reflecting the present-day average (2005–2019) and predicted future averages coloured by decade and scenario. World panels show LMEs with ‘low’ (left, <0.05), ‘medium’ (centre, >0.05, <0.5) and ‘high’ (>0.5) current habitat suitability (see Supplementary Fig. 18a for LME climate zones). b, Change in total area of habitat suitability (million square kilometres) within the north Atlantic, south Atlantic (SA), northwest Indian Ocean (NIO), southwest Indian Ocean (SIO), east Indian Ocean, west Pacific (WP) and east Pacific between present and future predictions, coloured by decade and scenario, with the periods 2046–2055 and 2086–2095 shown in the left and right panels, respectively, and symbols denoting negative (blue minus) and positive (red plus) change.

Fig. 3. Temporal trends in habitat suitability.

a, Monthly habitat suitability in which high (yellow) and low (blue) means are summarized within LMEs in the north Atlantic. Upper and lower panels show predicted future for each decade and scenario and present-day annual (2005:2019), respectively, within the Southeast US Continental Shelf LME (left) and Guinea Current LME (right). Axis labels 55 and 95 refer to decadal subsets 2046–2055 and 2086–2095, respectively. b, Interannual habitat suitability metrics where high (yellow) and low (black) means and high (large) and low (small) percentage coverage of core habitat area are summarized within EEZs, which are predominately located in the Atlantic Ocean. Left, middle and right, present-day annual (2005–2019), present-day average (2005–2019) and projected future (for each decade and scenario), respectively. Red boxes denote years referenced in the text when past climatic events of note occurred.

Fig. 4. Future redistributions in the context of global shipping.

a, Projected change in habitat suitability from baseline (absolute, 2005–2019) for 14 LMEs defined as medium importance, in which the result from a Kruskal–Wallis rank-sum test is shown at top left (χ² = 32.00, P = 5.93 × 10⁻⁶). Circles denote individual LME values, the thick line denotes the median and boxes bound the interquartile range (25th to 75th percentile), with whiskers extending to the maximum and minimum values. Upper and lower boundaries of violin plots extend to the maximum and minimum values, respectively, and width represents the density of observations. b, Global distribution of areas of high (yellow) and low (purple) shipping traffic density defined as the total count of vessels from a 2019 monthly average. c–e, These areas are shown in close-up in c–e, respectively. c–e, Areas of high (yellow) and low (purple) shipping traffic density from a 2019 monthly average (left) and areas of habitat suitability gain (red) and loss (blue) predicted from GAMs (right) shown in the national waters in the United States of America, marine region identification (ID), US part of the north Pacific Ocean (c); Sierra Leone, marine region ID, Sierra Leonian part of the north Atlantic Ocean (d); Japan, marine region ID, Japanese part of the eastern China Sea (e).

Fig. 5. Global habitat reshuffling leads to increased ship co-occurrence.

a, EEZ marine regions coloured by degree of change in SCI from the 2005–2019 baseline years. Red represents an increase in SCI and blue a decrease for 2100 ssp585. b, Percentage change in SCI from the 2005–2019 baseline years within each EEZ marine region, sorted and coloured by decade and scenario combination. c, Mean SCI calculated across EEZ marine regions, coloured by decade and scenario combination where the black dotted line represents present-day baseline SCI (2005–2019) and the percentage change from baseline is shown above each bar.

Extended Data Fig. 1. Model outputs and validation.

a, Regions of high (yellow) and low (blue) habitat suitability are indicated globally where regions have been joined together at boundaries (white border) based on current climatologies (2005–2019). b, Encounter locations sourced from the Ocean Biodiversity Information System (OBIS, n = 9,379) and c, sharkbook.ai for whale sharks (n = 13,267). d, Locations of four qualitative validation regions per major ocean (Atlantic: 1–4, Indian Ocean: 5–8, Pacific: 9–12) with corresponding e, mean monthly habitat suitability trends. f, The most important core areas for whale sharks currently (2005–2019) indicated by regions within quantile bands (50^th, dark blue; 75^th light blue; 90^th green; 95^th, light green).

Extended Data Fig. 2. North Atlantic habitat shift comparison with AquaMaps.

Regions of high (yellow) and low (blue) habitat suitability in the north Atlantic (NA) generated using the AquaMaps environmental envelope algorithm based on informed parameters from our tracking dataset (n = 1,021, left panel) and independent occurrences from Global Biodiversity Information Framework (GBIF) and Ocean Biodiversity Information System (OBIS) datasets from Food and Agriculture Organization (FAO) major fishing area 31 (n = 312, centre panel) projected into 2050 (decadal average of 2046–2055) under Representative Concentration Pathway (RCP) 5.8. Right panel shows our modelled projections for 2050 (decadal average of 2046–2055) under Shared Socioeconomic Pathway (SSP) 585. All models predict a band of habitats for whale sharks across the central Atlantic basin around equatorial latitudes.

Extended Data Fig. 3. Location and area of predicted changes in habitat suitability.

a, Regions of positive (red), negative (blue) or no change (or agreement, white) identified by both Generalised Additive Model (GAM) and Bayesian Additive Regression Tree (BART) algorithms in the east Pacific (EP), north Atlantic (NA) and east Indian Ocean (EIO) with top panel showing all regions of either positive or negative model agreement, centre panel showing regions of positive or negative agreement >0.1 or <−0.1, respectively, and lower panel showing regions of positive or negative agreement >0.25 or <−0.25, respectively for 2050 ssp585. b, Area (in million km²) of predicted change in habitat importance (>0.1 or <−0.1) identified by both GAM and BART models in the NA, south Atlantic (SA), northwest Indian Ocean (NIO), southwest Indian Ocean (SIO), EIO, west Pacific (WP) and EP coloured by decade and scenario were values on the left of each panel denote total area of change and points denote area of positive (red) and negative (blue) change. Mapped results for other regions are given in the Supplementary Figs. 9–15.

Extended Data Fig. 4. Latitudinal shifts in habitat suitability.

Arrows (left panel) indicate shifts in the geographic mean of the relative 90^th habitat suitability quantile (core habitat), with the grey mapped points indicating the current geographic mean, and the arrowheads the future geographic mean coloured by decade and scenario. Lines indicate the maximum and minimum latitudes where core habitat was projected, with the grey line indicating the current geographic limits and the future geographic limits coloured by decade and scenario. The lollipop plots (right panel) indicate the corresponding northerly (right upper panel) and southerly (right bottom panel) core habitat limit shifts (distance, km) within each region. Thresholds for core habitat are the relative percentile, calculated from the habitat suitability within each projection subset.

Extended Data Fig. 5. Changing habitats in national waters.

a, Percentage of Exclusive economic zones (EEZs, n = 200) summarised by continent where positive or negative changes in habitat suitability are projected for the low (blue, ssp585) and high (orange, ssp126) climate mitigation scenarios. b, EEZs ordered by the mean habitat suitability showing reshuffling from the baseline (2005–2019) under high mitigation/ sustainable development (ssp126, left to centre panel) and low mitigation/ high emissions (ssp585, centre to right panel) scenarios in 2086–2095. Individual EEZ trajectories are shown with Asia (orange) and Europe (blue) highlighted, demonstrating reshuffling under climate change with examples on the right (for example, where the Portuguese waters increase in relative habitat suitability scores, Papua New Guinean waters decrease).

Extended Data Fig. 6. Global known whale shark distribution.

The current global extant distribution of whale sharks as defined by the International Union for Conservation of Nature (IUCN) (grey) with expanded range (orange) mapped to identify current predicted and future projected regions where the species is likely to occur based on our models.

Extended Data Fig. 7. Change in ship co-occurrence index within national waters.

Exclusive economic zone (EEZ) marine regions coloured by degree of change in shipping co-occurrence index (SCI) from the 2005–2019 baseline years. Red reflects an increase in SCI and blue a decrease. Scenarios ssp126 and ssp585 are shown for the 2046–2055 (rows 1 and 2) and 2086–2095 (rows 3 and 4) decadal averages.

Extended Data Fig. 8. Ship co-occurrence index within national waters.

Exclusive economic zone (EEZ) marine regions coloured by shipping co-occurrence index (SCI, relative units where yellow is high and black is low) in the 2005–2019 baseline years (top panel) and 2086–2095 decadal average for ssp585 (bottom panel).

Extended Data Fig. 9. Model controls and performance.

a, Map of simplified whale shark (presences, yellow) and background (blue) location dataset in the north Atlantic used in the algorithm control procedure detailed in Supplementary Fig. 27 0a. b, Regions of high (yellow) and low (blue) habitat favourability in the north Atlantic identified by each algorithm. c, Internal validation performance metrics coloured by algorithm with the chosen method (Generalised Additive Models, GAM) in yellow. Note that after high Random Forest (RF) scores were disregarded due to concerns with data overfitting in mapped outputs (b), GAMs had high correct classification rate (CCR), precision, kappa and specificity scores second only to Bayesian Additive Regression Trees (BART), whereas Generalised Linear Models (GLM), Generalised Boosted Models (GBM) and Maximum Entropy (MXT) were better at correctly predicting presences (sensitivity and recall). d, Cross-fold validation performance metrics based on five spatial folds coloured by test where the black line denotes the median, boxes bound the interquartile range (IQR) (25^th to 75^th percentile) and whiskers extend to the smaller quantity of data extremes or medians ± 1.5× IQR with outliers shown as open circles and a red asterisk positioned above the chosen method (GAM). GAMs performed well across measures and outperformed BART (which showed better internal validation scores (c)), when measuring the true skill statistic (TSS) and Millers calibration statistic (MCS), with the closest MCS value to 1 across all models tested. e, Map of simplified whale shark (presences, yellow) and background (blue) location dataset in the north Atlantic used in the background sampling selection procedure detailed in Supplementary Fig. 27 0b. f, Internal validation performance metrics coloured by sampling method with the chosen method (MCP) in yellow. MCP had the highest score across tests. g, Cross-fold validation performance metrics based on five spatial folds coloured by test where the black line denotes the median, boxes bound the interquartile range (IQR) (25^th to 75^th percentile) and whiskers extend to the smaller quantity of data extremes or medians ± 1.5× IQR with outliers shown as open circles and a red asterisk is positioned above the chosen method (MCP). The MCP method performed well across measures. h, Regions of high (yellow) and low (blue) habitat favourability in the north Atlantic identified by each sampling method. The abbreviations in a–d refer to the algorithms tested and in e–h to the background sampling methods tested.

Extended Data Fig. 10. Projection method framework.

Steps taken to prepare global climate model (GCM) data for use in the study. Steps 1 to 6 were undertaken first to prepare the data for section 5 in Supplementary Fig. 27. Steps 7 and 8 are a summarised version of section 0a to 2b in Supplementary Fig. 27. Steps 1 to 5 were undertaken for each essential ocean variable (EOV) which were then stacked in steps 5 and 6, all of which was repeated for each decade and scenario combination. See Supplementary Information 3.5 for detailed description and equations used.