1. Study domain
We used the most-recent global marine regions dataset to define the spatial extent of the high seas (v11, https://www.marineregions.org/). We defined the high seas as marine areas outside of EEZs, which represent 95% of habitat on Earth by volume, 65% of the surface of the ocean, and 46% of the surface of the Earth22. We classified the high seas into three depth layers using the ETOPO1 bathymetry dataset41: epipelagic (0–200 m), mesopelagic (200–1,000 m) and bathyabyssopelagic (>1,000 m). We treated each depth layer as a separate planning domain. For every planning domain, we created an equal-area hexagonal grid of planning units of 2,620 km2 (~0.5° at the equator). This yielded a total of 90,065 planning units for the epipelagic planning domain, 88,528 planning units for the mesopelagic planning domain, and 87,170 planning units for the bathyabyssopelagic planning domain; there are fewer planning units by depth because of seamounts and underwater mountain ranges. Information about climate metrics, species, and fishing data were assigned to each planning domain (described below), each of which could be prioritised for inclusion in a climate-smart MPA network.
2. Climate change metrics: sources and processing
Climate-change metrics (i.e., climate velocity and the relative climate exposure index) were estimated using future ocean temperatures from a multi-model ensemble mean derived from 11 general circulation models (GCMs) from the Coupled Model Intercomparison Project Phase 6 (Earth System Grid Federation, https://esgf.llnl.gov; see Supplementary Table 2). We used models under three IPCC Shared Socioeconomic Pathways (SSPs)42: SSP1-2.6, SSP2-4.5 and SSP5-8.5. Pathway SSP1-2.6 represents an optimistic scenario, characterised by a shift to a more sustainable economy and reduction in inequality resulting in a peak in radiative forcing of ~3 W m-2 before 2100. SSP2-4.5 represents an intermediate scenario, with a stabilisation of radiative forcing levels at ~4.5 W m-2 by 2100. SSP5-8.5 is characterised by a continued increase of greenhouse gas emissions resulting from a fossil-fuel-based economy and increased energy demand, with a radiative forcing >8.5 W m-2 by 2100, rising thereafter.
For construction of the multi-model ensemble, we followed published methods18. Briefly, for each of the 11 climate models, we re-gridded from the original grid to a uniform 0.5° spatial grid using an area-weighted bilinear interpolation43. Then, we extracted depths according to three different ocean layers: epipelagic (0–200 m), mesopelagic (200–1,000 m) and bathyabyssopelagic (>1,000 m)44. For each depth layer, we averaged temperatures using a volume-weighting approach, with volumes of 0.5° grid squares in each depth layer. Finally, to avoid artefacts caused by inconsistent numbers of grid cells available by depth for different models, we included only grid cells common to all models within each depth layer18. All analysis was undertaken using software tools Climate Data Operators45 and R46 (Supplementary Table 2).
3. Climate velocity and the relative climate exposure metric calculations
To prioritise a climate-smart MPA network, we focused on areas across ocean depth layers that accomplish two objectives: 1) high retention of biodiversity, and 2) low levels of exposure to future climate warming. We used two metrics of climate change to represent these objectives: climate velocity and a Relative Climate-Exposure (RCE) index. Climate velocity is a metric that gives expectations for species’ range shifts under projected future ocean warming15,25. Thus, it is expected that in areas of slow climate velocity, species’ distributions are likely to shift less, promoting their retention within a given area. We estimate local climate velocity at 0.5° resolution for the second half of the century (2050–2100) at each ocean depth layer of the multi-model CMIP6 ensemble. The temporal trend (i.e., numerator of climate velocity) was calculated as the slope of a linear regression of mean annual temperatures (°C yr-1) for the corresponding climate scenario time period. The spatial gradient (i.e., denominator of climate velocity) was calculated from the vector sum of the latitudinal and longitudinal pairwise differences of the mean temperature across the corresponding climate scenario and time period at each focal cell using a 3×3 neighbourhood window (°C km-1)15. All calculations were performed using the VoCC R package47.
RCE is a metric we developed to obtain information about the amount of exposure to climate warming that local populations of a species would face relative to its experience of variation in seasonal temperatures. We calculated RCE as the ratio of the slope of a linear regression of projected mean annual temperatures (°C yr-1 2050–2100) to the current mean seasonal temperature range (°C from 2015–2020):
Therefore, it is expected that by prioritising areas with a low climate-exposure index we will select MPAs that will be more likely to minimise the potential exposure of species to future warming. Although some metrics of exposure have been already incorporated in marine spatial planning48–50, our metric considers not only information on warming, but also the relative local change relative to seasonal variation, and therefore, presumably, the relative vulnerability of resident biodiversity to projected future temperatures. For each planning unit in each depth domain, we calculated climate velocity and RCE separately for all three SSPs.
4. Conservation features
To solve the minimum-set conservation prioritisation problem, it is necessary to set a protection target a priori for each conservation feature, which indicates the minimum amount of each feature (i.e., species, habitat) to be included within the final prioritised network51. Here, our conservation features were marine species distribution maps from AquaMaps52 (v2019). The AquaMaps dataset predicts marine species distributions using a probability of occurrence (0–1) derived from an environmental niche model based on depth, temperature, salinity and oxygen at 0.5° spatial resolution. It includes 33,518 marine species, 23,700 of which we considered, as their environmental envelopes were generated using at least 10 observations 18,53. We used a minimum threshold of 0.5 probability of occurrence to select range maps for our high seas planning region as a conservative threshold defining core distribution ranges of the species10,18,53. Since most biodiversity in the ocean is located in coastal regions, our selection criteria yielded 11,701 species distribution maps in the high seas planning region (Supplementary Table 3).
We assigned each species distribution to each depth layer planning domain based on the species’ depth range derived from AquaMaps depth envelopes. This yielded a total of 8,538 conservation features in the epipelagic planning domain (0-200 m), 4,329 conservation features in the mesopelagic planning domain (200–1,000 m), and 1,472 conservation features in the bathyabyssopelagic planning domain (>1,000 m) (Extended Data Fig. 1d-f, Supplementary Table 3).
To obtain a better representation of conservation features in the global climate-smart MPA network across latitudes and ocean depth layers, we ensure that each is represented in every marine biogeographical province that it overlaps with. We used different biogeographical provinces for different depth layers. For the epipelagic planning domain, we used Longhurst Provinces54 and for the mesopelagic and bathyabyssopelagic planning domains, we used Glasgow Provinces55. Since most conservation features appear in multiple provinces, this categorisation process yielded a total of 56,996 conservation features: 19,761 for the epipelagic planning domain, 23,360 for the mesopelagic planning domain, and 13,875 for the bathyabyssopelagic planning domain.
5. Opportunity cost of fishing
When solving the minimum-set problem in spatial prioritisation, the objective is to identify areas where MPAs can meet conservation targets at a minimum cost. In our global prioritisation analysis, we used the value of fisheries as a cost layer because it is most common cost layer used in marine planning as resulting priority areas avoid valuable fishing grounds, where possible56. To estimate the cost ($US) of each planning unit in each depth layer planning domain, we multiplied an estimate of the total catch (kg) by the price per kg for each species caught ($US kg-1). The data were obtained from the Sea Around Us project dataset, which compiled fishing records for 1,242 species of fish and invertebrates from different publicly available databases from 1950-201430. This dataset includes landings at 0.5° spatial resolution, it is interpolated to account for missing values, and includes an estimate of illegal, unreported and unregulated fishing recorded as discards30. We calculated the total catch from 2005–2014 in each 0.5° cell (Extended Data Fig. 1a-c) and obtained the mean price of each species ($US)31. We then used the FishBase database (https://www.fishbase.se/) to categorise each species by depth (minimum and maximum preferences). This yielded a total of 1,099 prices for species ($US), 669 prices for the epipelagic layer, 322 prices for the mesopelagic layer and 118 prices for the bathyabyssopelagic layer. To obtain a final cost layer for each ocean depth layer, we calculated a total price by adding prices for every species within each 0.5° cell. We overlapped each cost layer generated with the corresponding high seas depth layer planning domain to obtain an area-weighted mean total cost ($US) for each planning unit (Extended Data Fig. 1a-c).
6. Spatial conservation prioritisation
We used integer linear programming (ILP) to find climate-smart MPA networks across ocean depth layers that minimise the overall cost (i.e., fishing value, US dollars) of the MPA network, and achieve the representation targets for each conservation feature. We solved the minimum-set objective51 using the Prioritizr R package32 and Gurobi optimisation software57. We set Gurobi to achieve a solution within 10% of the optimal solution. This is a relative arbitrary number that establishes the difference between the upper and lower bounds of the objective function58. For example, a value of 0.10 will result in the optimizer stopping and returning solutions when the difference between the bounds reaches 10% of the upper bound. The optimal solution is that which achieves the coverage target with the lowest possible cost.
Targets for conservation features were generated as follows (Extended Data Fig. 6). First, each conservation feature was intersected with both climate velocity and RCE maps to obtain its corresponding climate-change condition in each planning unit. Then, for each conservation feature, we selected only the planning units in which the conservation feature experienced slow climate change. We defined a slow climate change as values in the lowest quartile for each climate velocity and RCE metric (i.e., slow climate velocity and low RCE) within the range of the species under consideration (Extended Data Fig. 6). This process reduced the distribution size of each conservation feature (i.e., number of cells represented in the high seas) to one quarter of its initial distribution (Extended Data Fig. 6), yielding “restricted” conservation features. Next, for each of these restricted conservation features, we assigned targets based on the conservation status reported by the IUCN Red List59. For this, we used the rredlist R package60 to obtain the IUCN classification for each conservation feature. For taxa reported as threatened (i.e., VU = vulnerable, EN = endangered, or CR = critically endangered), every associated restricted conservation feature was assigned a fixed target of a 100%10 (i.e., 100% of its distribution within the first quartile of our climate-smart metrics for its overall distributional range) (Extended Data Fig. 6). For conservation features not reported as threatened in the IUCN Red List, we assigned relative targets based on the size of its global distribution (i.e., the restricted conservation feature) represented in each ocean depth planning domain, with a minimum of 10% for features with broad ranges, and 100% for features with limited ranges (10–100% of its restricted distribution for each conservation feature) (Extended Data Fig. 6, see Sensitivity Analysis section). We calculated these relative targets following the equation:
where Targetmax and Targetmin refer to the maximum and minimum target of protection (%) for restricted conservation feature i, PUs(i) represents the number of planning units represented for slow conservation feature i, and PUstotal refers to the total number of planning units for each high-seas depth planning region. By using relative targets instead of fixed targets, we ensured both the representation of multiple areas where widely-distributed restricted conservation features were conserved, as well as the protection of restricted conservation features within limited distribution ranges10.
We created different prioritisation planning scenarios to determine how the incorporation of climate-change metrics (i.e., slow climate velocity and low RCE) drives the selection of a climate-smart MPA network under alternative climatic futures. We ran three scenarios: one where we included restricted conservation features under SSP1-2.6; one with restricted conservation features under SSP2-4.5; and one with restricted conservation features under SSP5-8.5. Each prioritisation scenario locked in protection in existing MPAs (data extracted from www.protectedplanet.net) and Vulnerable Marine Ecosystems (VMEs, data extracted from www.fao.org).
To determine the spatial similarity of selected planning units among prioritisation scenarios and ocean depth planning domains, we calculated the Cohen’s Kappa coefficient. Cohen’s Kappa is a pairwise statistic that indicates the degree of agreement among scenarios, ranging from -1 to 1, where -1 represents complete disagreement, 0 represents agreement due to chance, and 1 represents complete agreement61.
7. Sensitivity analyses
The total area protected in a prioritisation analysis is highly correlated with the initial target assigned to each conservation feature. To test the sensitivity of our analysis to each conservation feature, we tested different targets of protection for the climate-smart prioritisation planning approach, with a minimum of 10% and adding 10% to a maximum of 90% (10–90% of the restricted distribution for each conservation feature). This gave a total of seven new and different planning scenarios (Supplementary Table 1). For each set of targets (i.e., planning scenarios), we performed a prioritisation analysis using the Prioritizr R package and the Gurobi software. As in the main design, we set Gurobi to achieve a solution within 10% of the optimum solution.
8. Caveats
There are several key sources of uncertainty in the current study. First, the climate-smart prioritisation planning approach taken to the conservation of individual species in this study assumes that every cell considered habitat in the AquaMaps dataset is of equal value. Scaled up, this implies conservation of any area-based target of a species distribution confers equal protection for the species. Recent analyses made the same assumption62 because these distributions are the best available data we have that cover the geographic and taxonomic scope of these studies. However, this approach does not take into account metapopulation structure, life-history, or how a population uses different habitat patches (e.g., for breeding, foraging, or migrating). It is clear that all parts of a species’ range are not equally important and that conserving any target for a species range will not necessarily confer the protection necessary to conserve metapopulation structure or even the species as a whole. We strongly encourage further research to improve our understanding of how species use and connect across the ocean, and the incorporation of this knowledge to improve global species distribution models like AquaMaps.
Second, our climate-smart prioritisation planning approach did not consider the three main aspects of a well-connected marine MPA network in its design63 (i.e., structural connectivity, functional connectivity, climate connectivity) and only represents the most cost-efficient climate-smart solution in a global scale analysis – meeting all the proposed conservation feature targets by minimising the fishing cost. But connectivity in the marine realm might be even more challenging to address, especially considering the three-dimensional nature of the ocean. This would require consideration of horizontal and vertical connectivity within and across ocean layers. These aspects were beyond the scope of this study but could be considered in further research. Third, there are limitations associated with the static cost layer of current fish prices in our analysis. The main objective of this study is to identify current climate-smart areas that will retain biodiversity and be least exposed to climate change to proactively plan for long-term ocean protection64. However, the cost layer used could also change as a result of changes in climate velocity and RCE which can derive in changes in the global distributions of fishes as the climate warms. Incorporating dynamic conservation features and cost layers using projections for biodiversity and fisheries could be a promising way forward for future research.
In this study, there is uncertainty associated with the climate models used to estimate the climate velocity and RCE. Although CMIP6 GCMs models have finer spatial resolution than CMIP5 models for both horizontal and vertical grid cells, uncertainty can still vary with depth18. These uncertainties might also generate differences in estimates of climate velocity and RCE metrics . This is an important consideration, especially for climate-smart planning approaches that rely on climate-change projections.
Finally, both climate velocity and RCE are generic metrics that do not directly link climate warming with species’ thermal preference. The RCE metric developed here indirectly relates vulnerability of marine biodiversity to climate warming (i.e., how sensitive they are to that level of exposure), but the analyses might be strengthened by better linking the hazard (i.e., climate warming), the exposure to that hazard (the amount of warming at the spatial unit of analysis), and the vulnerability of local marine communities to that hazard (how sensitive they are to that level of exposure). Again, these aspects were beyond the scope of the current study that included thousands of marine species, but could be considered in further research, should such data become available.