Climate change disrupts core habitats of marine species

Driven by climate change, marine biodiversity is undergoing a phase of rapid change that has proven to be even faster than changes observed in terrestrial ecosystems. Understanding how these changes in species composition will affect future marine life is crucial for conservation management, especially due to increasing demands for marine natural resources. Here, we analyse predictions of a multiparameter habitat suitability model covering the global projected ranges of >33,500 marine species from climate model projections under three CO2 emission scenarios (RCP2.6, RCP4.5, RCP8.5) up to the year 2100. Our results show that the core habitat area will decline for many species, resulting in a net loss of 50% of the core habitat area for almost half of all marine species in 2100 under the high‐emission scenario RCP8.5. As an additional consequence of the continuing distributional reorganization of marine life, gaps around the equator will appear for 8% (RCP2.6), 24% (RCP4.5), and 88% (RCP8.5) of marine species with cross‐equatorial ranges. For many more species, continuous distributional ranges will be disrupted, thus reducing effective population size. In addition, high invasion rates in higher latitudes and polar regions will lead to substantial changes in the ecosystem and food web structure, particularly regarding the introduction of new predators. Overall, our study highlights that the degree of spatial and structural reorganization of marine life with ensued consequences for ecosystem functionality and conservation efforts will critically depend on the realized greenhouse gas emission pathway.


| INTRODUC TI ON
Climate change and biodiversity change are intricately linked on land and in the ocean (Díaz et al., 2019). Understanding this link is key to evaluate risks for ecosystem functioning and nature's contribution to people as well as to conceive mitigation strategies to sustain the increasing demand for natural resources (Cheung et al., 2016;Jouffray et al., 2020). Recent observational studies on climate change impacts on biodiversity show that (i) marine biodiversity changes faster than terrestrial biodiversity (Blowes et al., 2019) and (ii) marine species track temperature changes more closely than terrestrial ones (Antão et al., 2020). The tighter link between climate and biodiversity change in the ocean can be explained by lower constraints on dispersal and colonization in the oceans (Lenoir et al., 2020;Poloczanska et al., 2013), lower thermal safety margins in marine species distributions (Pinsky et al., 2019), and the lack of microclimatic thermal refuges (Sunday et al., 2014).
The challenge remains to cast these observations into biodiversity scenarios for the future under different CO 2 emission pathways. One way is to use the observed poleward distribution shifts (Poloczanska et al., 2013 and the predictions of contractions or expansions of thermally suitable habitat (Cheung et al., 2009;Khan et al., 2018;Martinez et al., 2018) to predict global biodiversity change based on climate envelopes. Yet, changing temperatures are only one out of many factors driving species occurrence patterns (McHenry et al., 2019). More complex modeling approaches combine population models with climate models (Hare et al., 2010;Hollowed et al., 2009;Jones & Cheung, 2015) to predict future species distribution patterns. However, similar to trait-based vulnerability assessment frameworks (Hare et al., 2016;Jones & Cheung, 2018) these require a level of information on species-specific traits, which are not available for the vast majority of marine species and can consequently only deliver limited insights with regard to taxonomic coverage and spatial extent, that is, global marine biodiversity change.
Here, we project changes in global marine biodiversity toward the year 2100 through the integration of geo-referenced species occurrence data for 33,518 marine species based on global climate model simulations and an environmental niche model (Kaschner et al., 2006). Our model assigns probabilities of occurrence for every species and for 0.5° grid cells of the global oceans, resulting in habitat suitability maps that describe the potentially suitable habitat for a specific species under current and future environmental conditions. The occurrence probabilities are estimated by constructing environmental envelopes based on seven environmental parameters. These are water depth, water temperature (sea surface, sea bottom), salinity (sea surface, sea bottom), oxygen concentration (sea bottom), integrated primary productivity over the euphotic layer, sea ice cover, and distance to the coast. For our analyses, we focus only on potential core habitat area, that is, cells with probabilities of occurrence > 0.5. With this approach, we not only outpace previous efforts for marine biodiversity change scenarios  by enhancing taxonomic coverage by factor 2.6 but do so at a very fine scale in a multidimensional environmental envelope. We report on the effects of future climatic conditions on the areal extent of potentially suitable habitat space at a global scale under three contrasting Representative Concentration Pathways (RCP): strong mitigation (RCP 2.6), moderate mitigation (RCP 4.5) and no mitigation, that is, high-emission scenario (RCP 8.5). In addition, we evaluate the potential for ecological consequences with respect to trophic interactions and disruptions in distributional patterns.

| MATERIAL S AND ME THODS
We used the AquaMaps species distribution model (v.10/2019) to determine current and future habitat suitability maps. AquaMaps combines occurrence data with species ecological information from the literature to derive thresholds of species' environmental preferences and tolerances, referred to as environmental envelopes.
Environmental envelopes are matched to local environmental conditions to predict a species' natural distribution range in terms of probabilities of occurrence based on habitat suitability. The environmental parameters used to determine habitat suitability are sea surface temperature, salinity, integrated primary productivity over the euphotic layer, sea bottom oxygen concentration, sea ice concentration (SIC), and distance to coast.

| Species occurrence data
Geo-referenced species occurrence data were obtained for 33,518 marine species, primarily from global open-access biodiversity databases (Global Biodiversity Information Facility-GBIF, the Ocean Biodiversity Information System-OBIS) by data dump based on a list of currently accepted species names, validated using the Catalog of Life: 2018 Annual Checklist (https://catal ogueo flife.org/annua lcheck list/2018/). Point data were also gathered from museum collection records and scientific literature made accessible via FishBase and SeaLifeBase. We did not consider other ecological information (e.g., life stage, sex) from these sources because it is not consistently available for all species. However, it is expected that most species records in the data set comprise later-stage form observations. We also used proxy occurrence data, called country points, based on the reported presence (as native or endemic) of a species in a country. These consist of geographic coordinates of a fixed pair of points (one each over shallow and deep water), assigned to a specific country and marine FAO area combination, as a possible approximation of a species' locality of occurrence in that country. Either or both proxy points are used depending on the country of occurrence and depth distribution of the species. The proxy points were used to augment occurrence data required to model the distribution of a species, particularly for capturing environmental conditions in countries with reported species occurrence but no actual georeferenced data (see the section on Construction of environmental envelopes) and to adjust for biases from non-representative species distribution coverage. We also included additional species occurrence data (Teschke et al., 2019;Teschke & Brey, 2019a, 2019b, to close data gaps in high latitudes.
We corrected for biases resulting from species misidentifications and sampling efforts. Independent validation of the occurrence data was performed by filtering only for points that fall within a species' known range using bounding box coordinates, or in the absence thereof, the geographic limits of FAO areas where the species naturally occur from FishBase and SeaLifeBase. This filter eliminates occurrence data from beyond a species' native range from the model.
The filtered points were then assigned to a global grid of half-degree cells. The resulting data set consists of a unique set of cell center latitude and longitude coordinates, that is, a half-degree cell intersected by the point data is only registered once, regardless of how many times the same occurrence points fall within. This data set, referred to as "good cells," is used in the succeeding process of building species' environmental envelopes.

| Environmental data
The AquaMaps algorithm included georeferenced 0.5° cell resolution values of seven environmental parameters (see above) representing key physical and biological factors structuring the habitat of many marine species at large scales (Kaschner et al., 2019). The present-day state of habitat suitability was for RCP2.6 (CMIP5 range of 0.2 to 1.7°C), 1.8°C for RCP4.5 (CMIP5 range 1.1 to 2.6°C), and 4°C for RCP8.5 (CMIP5 range of 2.6 to 4.8°C, Collins et al., 2013;Giorgetta et al., 2013).
In addition, for the other employed environmental variables and on a regional scale, MPI-ESM-MR projected changes for 2045-2054 and 2090-2099 periods are close to the CMIP5 ensemble mean for most of the 66 large marine ecosystem (LME) defined by the Global Environment Facility. We also note that while we are using environmental data from just one model, we do use different scenarios, which at least for temperature cover a wider range of environmental conditions than the respective scenario model ensembles.

| Construction of environmental envelopes
Species-specific environmental envelopes were constructed using four thresholds: absolute minimum (Min A ), preferred minimum (Min P ), preferred maximum (Max P ), and absolute maximum (Max A ) to describe the environmental tolerances of species in terms of depth, sea temperature, salinity, primary productivity, dissolved oxygen, SIC, and distance to coast. Environmental envelopes for all parameters were calculated algorithmically, except for the depth envelope, which was based on absolute and common depth ranges recorded in FishBase and SeaLifeBase. Long-term average local conditions were extracted from the present-day environmental data layers at locations corresponding to the geographic coordinates of "good cells." Species environmental envelope thresholds were calculated for sea temperature (surface and bottom), salinity (surface and bottom), primary production, bottom dissolved oxygen, and distance to coast based on the following rules: 1. Min A = 25th percentile − 1.5 × interquartile or absolute minimum in extracted data (whichever is lesser) 2. Max A = 75th percentile + 1.5 × interquartile or absolute maximum in extracted data (whichever is greater) 3. Min P = 10th percentile of observed variation in an environmental predictor 4. Max P = 90th percentile of observed variation in an environmental predictor Species envelope thresholds for SIC were computed as follows: 1. Min A = absolute minimum in extracted data. For all species where Min A = 0, set Min A to SIC mean value + (−1). This extends Min A to avoid the exclusion of species from all non-ice-covered areas.
2. Max A = absolute maximum in extracted data 3. Min P = 10th percentile of observed variation in an environmental predictor 4. Max P = 90th percentile of observed variation in an environmental predictor Additional rules were used to avoid the use of illogical values and ensure conformity to basic biological concepts. In particular, environmental envelopes for temperature and salinity were based on annual means for surface layers for species with depth Min A ≤ 200 m, whereas bottom layers were used to compute the envelope thresholds for species having deeper Min A . For species with temperature Max P ≥ 25°C, temperature Max A was set to Max P + 4.2°C which sets the upper limit to around 34°C, approximating the lethal limit for tropical marine aquarium fishes. A minimum temperature distance of 0.25°C was also set between temperature Min P and temperature Max P for species with temperature Max A ≤ 5°C (polar and deepwater species), and 1.0°C where temperature Max A > 5°C. In addition, Min A was set to Min P -1 where primary production Min P -Min A <1.

| Habitat suitability maps
Species-specific environmental envelopes were used to predict relative habitat suitability (0 to 1), by scoring how well environmental attributes (i.e., local conditions) in a given 0.5° cell match a species' tolerances. First, the probability of species occurrence for each environmental parameter per 0.5° cell is calculated. The AquaMaps model assumes a probability of occurrence equal to 1 for average local conditions within a species' preferred range (Min P to Max P ) and decreases linearly as conditions approach the absolute minimum and maximum parameter thresholds of the species Furthermore, for marine mammals, the probability of occurrence in relation to depth is calculated using the mean depth in the cell (rather than the cell depth range).
Using a multiplicative approach to get P c allows each environmental parameter to act as a "knock-out" criterion, such that cells, where average conditions are outside a species' tolerance in any of the predictor parameters used, will yield a P c = 0.00 and be excluded as predicted suitable habitats.
To represent present-day core habitat distribution, all 0.5° cells with P c > 0.5 within a species' natural range, defined by its bounding box extents or the limits of its FAO area allocations, were used.
Mapping extents for species ranges based on FAO areas were extended into adjacent poleward FAO areas to allow for the emergence of natural range boundaries that could be artificially constrained by a species' FAO area limits.

| Statistical analyses
The code for all analyses described in the following sections is accessible under https://github.com/Dorot heeHo dapp/Code-Marin e-Biodi versi ty-Change.

| Analysis of biodiversity change
Estimated probabilities of occurrence were converted to presence/ absence (1/0) information using a cut-off value of 0.5. This value is arbitrarily chosen but has only a minor impact on global biodiversity patterns, which affects ranges of cosmopolitan or generalist species rather than distribution patterns of specialist or endemic species (Jones & Cheung, 2015). Hereinafter, we use the term core-habitat distribution to refer to the species distribution comprising cells with for each grid cell as well as relative species loss.
Current and future species range sizes were estimated by transforming the corresponding subset of coordinates and species occurrence probabilities in a raster object and adding up the individual cell sizes (km 2 ) obtained with the function area() from the "raster" R package (Hijmans, 2020). The computed current and future ranges were compared for the estimation of species range size changes until 2100.
We also computed several landscape metrics, for example, the number of patches, the mean size of patches, the average distance between patches, and the patch density using the according command from the R package "landscapemetrics" (Hesselbarth et al., 2019).
Current and future north and south poleward (leading) and equator-bound (trailing) species distribution edges were defined as the 95th percentile (poleward) or 5th percentile (equatorward) of the latitudinal range of a species' core habitat extent for each hemisphere. The distance between the species current and future range edge (2100 under the three RCPs) was then calculated using Vincenty's formulae, which assumes an ellipsoidal shape of the earth (distm(), R package "geosphere" (Hijmans, 2019)).

| Equatorial gap analysis
For every species that has current and projected distributional ranges on both hemispheres under all three RCP scenarios (n = 21,648), we determined the 0.01-quantile of its latitudinal distribution in the northern hemisphere to be the southernmost edge of the core-habitat area (i.e., probability of occurrence >0.5) with a value of 0 indicating distributions crossing the equator. The northernmost edge in the southern hemisphere of the distributional range of each of these species was identified accordingly.
We then calculated the distance between these equator-bound distribution limits under current and future projected conditions again assuming an ellipsoidal shape of the earth (Vincenty, 1975) as implemented in the distm() function of the R package "geosphere" (Hijmans, 2019).

| Top predator loss and invasion
We defined fish species with a trophic level >3.9 (according to FishBase, https://www.fishb ase.de/) as fish top predators, resulting in 2343 species. We then determined how many of these species disappeared from or emerged in marine ecoregions (Spalding et al., 2007) in the future projections under the three RCP scenarios.

| Impact of environmental variables on compositional change
To assess impact strength across the environmental variables on future changes in species composition (Jaccard turnover), we ran random forest regression models (Breiman, 2001) (R packages "randomForest" (Liaw & Wiener, 2002), "ranger" (Wright & Ziegler, 2017), and "caret" (Kuhn, 2020)). As predictors, we included the current mean values of the environmental variables depth, mean annual sea surface temperature (SST), surface salinity, bottom oxygen, primary productivity, annual SIC, and distance from land as well as the difference between their current and projected future values. Only changes in SIC were excluded from the analysis as they were highly correlated with current SIC.
We applied hyperparameter tuning to find the best values for the following parameters: (i) the number of explanatory variables to randomly sample from at each node (6) and (ii) the minimum number of samples within the terminal nodes (3). We also tested for optimal sizes of training (0.8) and test data subsets (0.2) to find the set of parameter values for a best possible trade-off between bias introduction and overfitting the model.
The impact of each predictor variable across all trees and nodes was evaluated by applying permutation feature importance, which assesses the effect of re-shuffling each predictor variable's values on the overall model performance as implemented in the 'ranger' R package (Wright & Ziegler, 2017).

| Poleward range shift and loss of core habitat area
Our projections reveal considerably different outcomes among the three CO 2 emission scenarios. Under RCP 2.6 and 4.5, poleward expansions of the leading edges of species distributions prevail over equatorward retractions of the trailing edges (Figure 1a). In other words, for most species the latitudinal range remains the same or even increases as species move poleward. This coincides with empirical findings on recent marine species range shifts (Fredston-Hermann et al., 2020; Poloczanska et al., 2013Poloczanska et al., , 2016. However, under RCP 8.5, the velocities of leading and trailing edges change, resulting in approximately equal latitudinal shifts of leading and trailing edges or even greater trailing edge shifts compared with the leading edges. In consequence, species would lose low-latitude range faster than gain highlatitude range. The predicted larger shifts in the northern hemisphere compared with the southern hemisphere have been attributed to greater warming in northern hemisphere oceans (Lenoir et al., 2020).
Importantly, our analyses on the total extent of all suitable habitat patches reveal that despite the commonly reported range shifts toward the poles and latitudinal range expansions, the majority of species will suffer net range loss of core habitat area (Figure 1b,c) in all three CO 2 emission scenarios. The reduction in the number and size of suitable habitat patches as well as growing isolation ( Figures S7-S9) between them increase with projected greenhouse gas emissions, leading to a 50% range contraction in suitable core habitat area for 1%, 4%, and 48% of species under RCP 2.6, 4.5, and 8.5, respectively.
The extent of core habitat loss shows some variation across hemispheres and CO 2 emission scenarios. For exclusively northern hemisphere species we find much higher variability in changes of suitable habitat area and slight median gains under RCP 2.6 and RCP 4.5, in contrast to median range size losses of southern hemisphere species.
In addition to the different rates of warming in the southern and northern hemispheres, one possible explanation for this pattern could be the Antarctic circumpolar current (Böning et al., 2008) whose sharp gradients in temperature and other environmental parameters act as a dispersal barrier for many species. However, under RCP 8.5, median range size loss is projected to increase dramatically for northern hemisphere species outpacing range size losses of species in the southern hemisphere ( Figure 1b). Part of the reason for this pattern could lie in the higher fraction of area with substantial temperature increases in the northern compared with the southern hemisphere (see Figure S10).

| Changes in biodiversity
Overall, projection medians show net local increases in species richness under RCP 2.6 and 4.5 scenarios, but net local decreases under RCP 8.5 with higher variability in subpolar regions and semi-enclosed tropical and temperate seas (Figure 2). Species losses are expected F I G U R E 1 Median and confidence intervals of latitudinal range shifts and changes in range sizes of marine species suitable habitat area taking place between the beginning and the end of the 21st century under three CO 2 emission pathways (RCP2.6, RCP4.5, RCP8.5). (a) Poleward shifts of marine species with cross-equatorial distributions (left), and species with only northern hemisphere (center) or southern hemisphere (right) occurrences. (b) Percentage range size change. (c) Median relative range size change for different marine organism groups. Error bars represent lower (25%) and upper (75%) quartiles. n represents the sample size of the respective species group.
F I G U R E 2 Species turnover maps showing changes in marine species composition between the beginning and the end of the 21st century under three Representative Concentration Pathways. Side panels represent latitudinal gradients of the median number of local species loss (red) and gain (blue) according to 25 and 75% quantiles (shaded area). The light bands in the eastern Pacific are an artifact of hard longitudinal distribution bounds applied to species with limited occurrence data, or to species whose range limits stop at islands in the central and eastern Pacific (see Methods for further details).
to be highest at the equator and in the (sub-)tropics and will be particularly pronounced in the Caribbean Sea, eastern Mediterranean, Red Sea and Indo-Pacific region ( Figure S5).
Net species gains are achieved through invasion without loss in higher latitudes, including Arctic and Antarctic regions. However, in the tropics, where the probability of newly invading species is low, projected changes in habitat suitability could lead to extreme reductions in the number of species (Figures S11-S14). Overall, 4%, 7% and 19% of the global ocean area will become unsuitable as core habitat for at least a third of their current species inventory under RCP2.6, RCP4.5 and RCP8.5, respectively ( Figure S13). Cells with very little (< 10%) or no change in community composition by 2100 represent 36% under RCP2.6 but only 15% and 1% under RCP4.5 and 8.5. According to our model, some areas will experience high rates of compositional reorganization much sooner than 2100, that is, assuming RCP8.5 emission pathways our projections show drastic changes in the polar regions and parts of the Mediterranean and Red Sea already by 2050 ( Figure S1).
Results of a random forest analysis assessing the impacts of environmental parameters and changes thereof on the amount of compositional turnover show that the highest rates of biodiversity change are driven by SST, primary production, and SIC ( Figure S2).
The results concur with the high turnover rates in polar and tropical regions shown in Figure 2. The relative importance of these parameters differs across RCP scenarios ( Figure S2). In the low-emission scenario (RCP2.6), our model predicts the biggest changes in species composition for areas with low and further declining primary production (see Figure S3 and maps in S20 and S22) and areas with low temperature and high SIC, that is, polar regions ( Figures S3 and maps in S20 and S23). With increasing emissions, the effects of SIC and SST increasingly prevail over primary production and other environmental variables like salinity or sea bottom oxygen ( Figures S2, S24 and S25). Under RCP8.5, we see the highest predicted compositional turnover in areas of very low and very high current SST as well as in regions with currently high SIC ( Figure S3).

| Distributional gaps around the equator
Poleward shifts and range contractions of suitable habitat can affect the spatial integrity of species distributions. An illustrative example is species with cross-equatorial distributions. Our model predicts the emergence of areas of unsuitable environmental conditions around the equator for 8% (RCP2.6), 24% (RCP4.5), and 88% (RCP8.5) of the 11,853 species in our data set that currently show continuous cross-equatorial distribution ranges (Figure 3).

| Loss and invasion of top predators
Top predators in general and large ocean predators, in particular, constitute important top-down food web controls in the ocean that affect marine ecosystems across a broad range of spatiotemporal scales (Hazen et al., 2019). Therefore, the future extent of the loss, displacement and/or replacement of these species in a particular region will have broad-scale implications for community integrity and functional and compositional stability (Estes et al., 2011;Myers et al., 2007). We estimate that even under the strong mitigation scenario (RCP2.6) > 65% of all fish top predator species in our data set will extend their ranges to new ecoregions. Under RCP8.5 this fraction of species increases to 92%. At the same time, 1% and 14% of marine ecoregions could lose more than half of their current fish top predator species under RCP4.5 and 8.5 due to their shifting ranges.
Both processes can lead to substantial changes in food web structure and trophic dynamics (Fossheim et al., 2015;Smith et al., 2017).
Especially in the polar regions, top predator species experience anon average-higher turnover than the remaining species ( Figures S4   and S6). In the Arctic, projected changes in habitat suitability will allow for increasingly frequent fish top predator invasions, while in the Antarctic, the picture is more differentiated across regions. Here, areas are projected to mainly experience either loss, gain or replacement of current fish top predator species assemblages ( Figure S5).

| Limitations of the modelling framework
The vertical lines on the maps showing projected species composition turnover (Figure 2) are an artefact of the AquaMaps modelling approach. The root cause of these is potentially data-deficient pelagic-oceanic or bathypelagic species whose range maps include predicted highly suitable habitat at their currently known or published range limits of the native range, defined by either a species' bounding box coordinates or, in the absence of such, the limits of FAO areas where a species is reported as native. These limits are used in the AquaMaps model when plotting predicted probabilities of species occurrence.
The vertical bands east of the Hawaiian islands around 154°W may also be indicative of a real dispersal barrier, called the East Pacific Barrier (Briggs, 1961), an uninterrupted expanse of ocean between 5000 and 8000 m deep and nearly 6500 km wide, that separate distinct fauna in the central and Indo-West Pacific from the eastern Pacific. The vertical band seen in this area is formed by the 10° outward extension of current range limits that are applied in the AquaMaps model when plotting the future range extent of a species. Thus, for species with an eastern limit in Hawaii, the hard edge is moved eastward to around 144°W in the future range, resulting in the banding we see. Similar cases appear to result from species whose eastern limit produces hard edges in the area of Ducie Island (~124°W, extending to 114°W) and Easter Island (109°W, extending to 99°W). Similar banding occurs from this 10° outward future range extension for species whose eastern limit is set to either FAO area 61, 71, or 81. Thus, faint bands are also seen around 175°W-165°W (FAO 61 and 71), and 120°W-110°W (FAO 81).

| DISCUSS ION
Our projections of future suitable habitat ranges come with all the limitations and uncertainties that have frequently been discussed for environmental niche modelling approaches (Araújo et al., 2019).
Limitations include for instance the equilibrium assumption, that is, that species are found in all suitable areas with regard to environmental and climate conditions, whilst being absent from unsuitable habitats, which does not always hold (Araújo & Pearson, 2005). Additionally, marine species occurrences are not evenly reported across the globe, resulting in much more complete species inventories in highly sampled areas such as coastal areas in the northern hemisphere, whereas polar regions, ocean gyres and the tropics are clearly under-sampled. We would, therefore, like to stress that AquaMaps provides estimates of habitat suitability, that is, it assigns probabilities of occurrence for each species in each cell given the local environmental conditions and a species' ability to occur under these conditions based on its speciesspecific environmental envelope. Despite these limitations, we argue that our estimates offer an informed picture of global biodiversity dynamics and are, in fact, conservative due to (i) ignoring species interactions which affect a species' chance to establish at a new location in addition to habitat suitability (Wisz et al., 2013), (ii) the consideration of annual average temperatures as species extinctions are often driven by extreme values rather than annual means of environmental variables (Cheung & Frölicher, 2020;Román-Palacios & Wiens, 2020;Vasseur et al., 2014), and (iii) our assumption that species will be able to track climate velocities to adjust distributional ranges. The latter assumption is based on findings from metaanalyses (Gregory et al., 2009;Poloczanska et al., 2013), but there are studies on single species or species groups indicating a mismatch between range shift and changing environmental conditions (Alabia et al., 2018;Gaudin et al., 2018). AquaMaps relies on single model estimates as opposed to model ensemble means, which are commonly used for climate projections (see Material and Methods section for details). Nevertheless, we argue that a multimodel average approach would not necessarily provide a more robust assessment of future climate conditions, since most models have common biases and similar large uncertainties in their projections (as discussed in Power et al., 2021). Regarding the more recent CMIP6 model ensemble, which is driven by a different but similar in terms of radiative forcing set of scenarios, the RCP8.5-derived warming in our analysis for the LMEs is about 0.5°C weaker than the unconstrained SSP85 ensemble mean at the end of the 21st century, indicating that our analysis is on the conservative side of potential future climate projections. However, if constraints from observations of past climate change are applied to the CMIP6 ensemble, the overall warming is reduced for all SSP scenarios compared with the unconstrained scenarios, bringing warming rates back close to the CMIP5 ensemble (Lee et al., 2021;Tokarska et al., 2020). Furthermore, simulation's improvement from CMIP5 to CMIP6 are overall relatively small and not always evident (e.g., Dieppois et al., 2021;Séférian et al., 2020) and while ocean model resolution has increased for a range of CMIP6 models, it is still coarser than the AquaMaps grid (1.0° vs. 0.5°).
Our results show that ongoing climate change will cause many species to suffer range contractions and disruptions accompanied by a general loss of suitable habitat space, despite expected expansions of latitudinal ranges toward the poles. According to our predictions on habitat range, latitudinal distribution shifts will be accompanied by substantial loss of potential habitat, which has so far been overlooked, as previous global scale studies concentrated on assessing latitudinal range shifts (Cheung et al., 2009;García Molinos et al., 2016) and differences between leading and trailing edges (Poloczanska et al., 2013). However, changes in habitat suitability may illustrate ecological change better than trends in latitudinal range shifts. Adverse local environmental conditions together with physical barriers can lead to the loss of suitable habitat patches within latitudinal ranges. In addition, on a sphere, the same latitudinal range translates to smaller areas at higher latitudes, resulting in loss of potential habitat if distributions are shifted poleward. Our results concur with multifactor habitat suitability model projections (McHenry et al., 2019;Roberts et al., 2020) and individual species examples of observed and predicted marine habitat loss as a consequence of global warming (Assis et al., 2016;Dahlke et al., 2018;Filbee-Dexter et al., 2016;Hoegh-Guldberg et al., 2007;Martinez et al., 2018).
With respect to the high numbers of species with complete loss of suitable core habitat in tropical areas, two aspects should be noted: (1) Our calculations are based on core habitat area, that is, cells with a probability of occurrence > 0.5 for a given species. Hence, even though the equatorial regions in large parts do not satisfy the core habitat requirements of marine species under RCP8.5 toward the end of this century, this does not imply that marine life will completely disappear from these areas. (2) The definition of upper physiological limits of the thermal tolerance (i.e., environmental envelope) of marine species inhabiting the warmest parts of the oceans is difficult and can only be estimated based on experimental results, which do not exist for the vast majority of species. Our model does permit potential adaptation of tropical species to further temperature increases (see method section for further details), but that permission may be too optimistic and any interpretation of model projection in these areas should consider these uncertainties.
The ecological consequences of the depicted habitat loss are likely far-reaching as many habitat-forming species (e.g., kelp, sea grass, corals; Figure 1c) are among the organism groups with very high median contractions in potential range size. Contractions, disruptions or shifts in their distributions may seriously compromise habitat conditions for many more associated species that depend on food or shelter provided by the habitat-forming species (Steneck et al., 2014). We identify further ecological consequences of the ongoing alterations in local marine species inventories with regard to trophic structure, loss of keystone species and global connectance of populations. Indeed, recent studies have observed disruptions of for example, circumglobal fish trans-equatorial migration (Carlisle et al., 2017), central equatorial Pacific coral reef cover and seabirds' presence and reproductive activity (Brainard et al., 2018), as a consequence of increasing hypoxia and thermal stress during extreme climatic events. This together with the general climate change induced decline in species number around the equator (Chaudhary et al., 2021) supports our projections. Ecological implications of emerging distributional gaps are, therefore, changes in spatial dynamics of habitat use, novel or irregular trophic interactions and reduced effective population sizes. The latter may increase the susceptibility of species to localized depletion and reduce genetic pool size and exchange, which could eventually lead to the emergence of new species and bipolar distribution patterns that are a common phenomenon now and in the past (Lindberg, 1991).
Such fundamental spatial and structural reorganization of marine life calls for renewed efforts to increase the area, connectivity, and integrity of natural systems and to reduce the number of threatened species. It also emphasizes the relevance of implementing dynamic management plans and international collaborations for effective global climate change mitigation and biodiversity conservation. International agreements like the UN Convention on Biological Diversity (CBD/COP/15/1/Add.3; www.cbd.int) and the "30×30" target (i.e., expand protected areas to cover 30% of the planet's land and sea by 2030) or the EU Biodiversity strategy for 2030 (European Commission, 2020) can be very valuable political frameworks to reach these goals. Indeed, natural and fully protected marine landscapes, when carefully planned and managed, can be vital aspects in mitigating the ongoing changes in climate and biodiversity (Jacquemont et al., 2022). Our study highlights that the degree of range contraction and loss of suitable habitat will critically depend on the realized greenhouse gas emission pathway. In consequence, the realized scenario will also determine the costs, effort and success of future conservation measures implemented to meet global conservation goals to sustain marine ecosystem functionality.