Climate Change-Induced Emergence of Novel Biogeochemical Provinces

The global ocean is commonly partitioned into 4 biomes subdivided into 56 biogeochemical provinces (BGCPs) following the accepted division proposed by Longhurst in 1998. Each province corresponds to a unique regional environment that shapes biodiversity and constrains ecosystem structure and functions. Biogeochemical provinces are dynamic entities that change their spatial extent and position with climate and are expected to be perturbated in the near future by global climate change. Here, we characterize the changes in spatial distribution of BGCPs from 1950 to 2100 using three earth system models under two representative concentration pathways (RCP 2.6 and 8.5). We project a reorganization of the current distribution of BGCPs driven mostly by a poleward shit in their distribution (18.4 km in average per decade). Projection of the future distribution of BGCPs also revealed the emergence of new climate that has no analog with past and current environmental conditions. These novel environmental conditions, here named No-Analog BGCPs State (NABS), will expand from 2040 to 2100 at a rate of 4.3 Mkm2 per decade (1.2% of the global ocean). We subsequently quantified the potential number of marine species and annual volume of fisheries catches that would experience such novel environmental conditions to roughly evaluate the impact of NABS on ecosystem services.


INTRODUCTION
The biosphere is partitioned according to well-defined biomes corresponding to unique sets of physical, chemical and biological conditions. On land, biomes such as "desert, " "mountain, " or "forest" can be delineated by environmental variables such as temperature, precipitation and altitude (Peel et al., 2007;Sayre, 2014). Similarly, the ocean has been subdivided in various ways, using either expert opinion or data driven methods based on both biotic and abiotic components of the marine realm (Kavanagh et al., 2004;Oliver and Irwin, 2008;Reygondeau and Dunn, 2018). The widely accepted system of ocean provinces elaborated originally by Longhurst (2007) covers both open ocean and coastal zones. Two levels structure Longhurst's province system. The first one distinguishes coastal and open ocean areas and subdivides them into temperate, tropical and polar biomes. The second sub-divides each of the coastal and oceanic biomes into oceanographically, ecologically and topographically homogeneous regions. This results in 56 distinct biogeochemical provinces (BGCPs) that have been shown to be closely related to global patterns in marine biodiversity (Pauly, 1999;Beaugrand et al., 2000), ecosystem functions and services such as fisheries productivity (Chassot et al., 2011;Demarcq et al., 2012;Reygondeau et al., 2012).
While the relevance of Longhurst's BGCPs to partition the ocean into homogeneous environmental and ecological regions has been established , their delineation has been modified to improve their representativeness of the biogeography of diverse biological communities (Reygondeau et al., 2013). Specifically, Reygondeau et al. (2013) have shown how seasonal and inter-annual climate variability -such as such as such as El Niño Southern Oscillation (ENSO) event-modifies BGCPs' boundaries and drives their reorganization.
Climate change is expected to lead to drastic changes in ocean conditions. Oceans have already absorbed more than 93% of the heat resulting from the accumulation of greenhouse gases. They are rapidly getting warmer, less oxygenated (Gattuso et al., 2015) and changes in the seasonality of oceanographic conditions have already been reported. These oceanographic changes will alter the delineation and location of the BGCPs. Notably, some variables are expected to reach levels that are beyond those observed in recorded history (Froelicher et al., 2016). As a result, existing BGCPs might not be able to characterize regions where changes are projected to exceed the range of past observations. Such "novel" climatic zones have already been identified in the terrestrial realm. They have also been reported during marine heatwaves (Frölicher et al., 2018;Oliver et al., 2018).
In this study, we test the hypothesis that climate change might result in a large-scale bio-geographic re-organization of the oceans accompanied by the emergence of novel BGCPs. In particular, we examine the extent to which these novel BGCPs might overlap with marine biodiversity and fisheries (Cheung et al., 2016b;Cheung, 2018). For this purpose, we use the numerical approach described in Reygondeau et al. (2013) to analyze projected ocean conditions from three Earth System Models (ESM; Institut Pierre Simon Laplace, Max plank Institute and Geophysical Fluid Dynamics Laboratory models) under two climate change scenarios (the "strong mitigation" Representative Concentration Pathway, or RCP 2.6, and "business-as-usual" RCP 8.5). For each ESM and RCP scenario we calculate the projected annual average distribution of each BGCP from 1950 to 2100. We then identify regions that are characterized by novel combinations and levels of oceanographic variables and that have therefore never been observed historically. We name these novel BGCPs No Analog Biogeographical State regions (NABS). We quantify the timing of NABS emergence and evaluate the potential consequences on marine biodiversity and key marine ecosystem services.

Environmental Data
The ocean properties that we used to delineate the BGCPs from 1950 to 2100 were based on outputs from three ESMs. They include annual average surface (average between 0 and 10 m)and bottom sea water temperature ( • C), oxygen concentration (ml.L −1 ), salinity, pH, surface net primary production (mgC.km 2 .year −1 ),particulate organic carbon concentration (mgC.km 2 .year −1 ) and sea ice coverage (%). These data have been simulated by the Geophysical Fluid Dynamics Laboratory Earth System Model (GFDL ESM2M), the Institut Pierre Simon Laplace Climate Model (IPSL-CM5-MR) and the Max Planck Institute for Meteorology Earth System Model (MPI-ESM). They were made available through the Coupled Model Intercomparison Project phase 5 (CMIP5; Taylor et al., 2012). The ESM outputs were re-gridded onto a regular grid of 0.5 • of latitude by 0.5 • of longitude using the nearest neighbor method and values in coastal cells were extrapolated using bilinear extrapolation. In addition, annual climatologies of euphotic depth (Morel et al., 2007), mixed layer depth (de Boyer Montégut et al., 2004), and bathymetry (Smith and Sandwell, 1997) were also used to delineate BGCP boundaries. These variables were gathered from observations and are here used to optimize the quantification of the environmental envelopes of specific BGCPs such as frontal oceanic systems, coastal regions or boundary of oceanic subtropical gyres.
The ESM outputs for the future period (2006-2100) were simulated under RCP 2.6 and RCP 8.5. The RCP 2.6 is a lower emission "strong mitigation" scenario under which the radiative forcing trajectory peaks at 3 W.m −2 before 2100 and then is followed by a decline to 2.6 W.m −2 by 2100. The RCP8.5 is a high emission "business-as-usual" scenario with a rising radiative forcing pathway leading to 8.5 W.m −2 by 2100.

Modeling the Geography of BGCPs
To delineate the BGCPs' boundaries for each year from 1950 to 2100, we applied the numerical approach described in detail in Reygondeau et al. (2013). The approach quantifies the environmental envelopes of each BGCP based on the variables described above and projecting their future spatial distribution (Figure 1). First, the spatial coordinates of the boundaries for each BGCP as defined by Longhurst (2007) were retrieved ( Figure 1A). Second, to include the monthly variation of each variable and characterize the discrete set of environmental conditions that characterize each BGCP, we obtained monthly values of the ESM outputs for each variable and spatial cell in each BGCP for the time period January 1970 to December 2000. This time period represents the timeframe over which the BGCP delineation was defined originally. We here used all environmental variables (surface and bottom) for BGCP located in the coastal biome and only surface environmental variables for BGCP located in the tropical, temperate and polar biomes (see Supplementary Table 1). We obtained 56 matrices of environmental properties (one for each province) named X n1,p,z (n1 = 360 months, p = x environmental variables, z = number of geographical cells of the selected province). Third, we quantified the environmental envelope of each BGCP using three environmental niche models (ENM): Non-Parametric Probabilistic environmental niche model , Maxent (Phillips et al., 2004), and boosted Regression Trees (Elith et al., 2008). We treated spatial cells included in each distinct BGCPs as a set of "occurrence records" analogous to species distributions and applied each ENM to model the environmental envelopes of each BGCP. We applied a multimodel ensemble approach to capture the uncertainty associated with the different numerical methods. Finally, we computed the time-dependent spatial distributions of the probability of occurrence (values ranging between 0 and 1) of each BGCP from the annual average ENM outputs for each ESM and each RCP scenario. Global BGCP division is then identified by attributing the geographical cell to the BGCP with the highest probability of occurrence for each year in each geographical cell and year.

Comparing Predictions With Observations
We compared the spatial distribution of and temporal fluctuation in the BGCPs predicted from ESM outputs with the ones gathered using observed environmental variables in Reygondeau et al. (2013). To harmonize the two sets of BGCP predictions, the BGCP distributions were averaged annually over the same period (1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007) and regridded using the nearest neighbor method on a 1 • ×1 • spatial grid (original spatial grid used in Reygondeau et al., 2013). We first performed a spatial correlation of the average probability of presence from 1998 and 2007 to evaluate the level of congruence between observed and modeled BGCPs. We then performed a temporal correlation of the standard deviation of the probability of occurrence for the period 1998-2007 for each BGCP. Results from the analyses were used to evaluate the ability of ESM outputs to represent the temporal variability of the probability of a given BGCP (Supplementary Figure 1).

Geographical Trends
For each BGCP, we calculated changes in two geographical indices over time and compared the results between climate change scenarios. First, we calculated the annual centroid of each BGCP as the average coordinates of the center of the spatial cell belonging to a BGCP weighted by the value of the probability of occurrence of the same BGCP in the analyzed cell (Figure 1). Second, we calculated the total area covered by each BGCP based on the sum of area of the spatial cells belonging to the a given BGCP for each year. We calculated these indices using the ensemble average BGCP distributions across ESM and ENM under each RCP from 1950 to 2100. We then evaluated the latitudinal shifts and changes in total area of the BGCP from 1950 to 2100 (relative to the average of 1970-2000) under RCP 2.6 and RCP 8.5.

Identification of No Analog Biogeographical State
No-Analog BGCPs State (NABS) represents a multienvironmental range of condition that has not been encountered during the training set period  across the BGCPs. More precisely, a NABS set of environmental conditions is significantly different from all possible environmental combinations encountered from 1950 to 2000 and used to inform the multi variable environmental envelopes of the 56 BGCPs. Numerically, NABS are characterized by spatial cells exhibiting a null probability of occurrence for all of the 56 Longhurst BGCPs. We then quantified NABS distribution and coverage by delineating their annual boundaries and total area for each ESM and RCP. We identified the year at which a spatial cell first became a NABS and attributed a confirmed NABS location if the environmental conditions are similar (Null probability of original and spatially surrounding BGCPs) for more than five consecutive years for a given ESM and RCP pathway. We subsequently mapped the distribution of NABS as the 2/3 agreement between the 3 ESMs over time.

Assessing the Potential Impacts of NABS on Biodiversity and Ecosystem Services
The exposure to novel ocean conditions in NABS might seriously challenge the viability of marine species currently inhabiting those areas, thus leading to substantial species turnover, changes in catch composition, and declines in potential fisheries production (Cheung et al., 2016b). To evaluate the risk that the emergence of NABS may pose to marine biodiversity and ecosystems services, we calculated the number of marine species and the annual volume of fisheries catches that will be exposed to NABS conditions. For marine exploited biodiversity, we collated the global gridded marine species richness dataset (Gagné et al., 2020) and added several other exploited species distributions from the Sea Around Us 1 and species used in Asch et al. (2017). 1 www.seaaroundus.org This dataset included 1,105 species ranging from invertebrate to top predator. We extracted average annual total fisheries catch (average between 2001 and 2015) from the Sea Around Us database. Both datasets were in the 0.5 • latitude × 0.5 • longitude grid consistent with the ocean conditions and BGCPs data. We then computed the number of species and total annual catch in spatial cells characterized as NABS for each year between 2000 and 2100 across ESMs and RCPs.

RESULTS
We projected that 26-39% of BGCPs will expand their area, while between 61 and 74% will shrink in size by 2100 under RCP 2.6 and 8.5, respectively (Figures 1, 2 and Supplementary Tables 1, 2), due to changes in environmental conditions (Supplementary Figure 2). Specifically, the tropical BGCPs were projected to expand while the temperate BGCPs were projected to contract by losing area on their tropical side, without compensation via poleward expansion. The two polar BGCPs were projected to shrink in size. This redistribution of BGCPs was moderate by 2050, but projected to be substantial by 2100 under RCP 8.5. The level of such redistribution was found to be more important over the Indo-Pacific basin with a high spatial perturbation located around the warm pool and monsoon provinces (WARM and MONS; Figure 1 and Supplementary Table 1), altering the position of all neighboring BGCPs that tend to move poleward. The change in the distribution can be attributed to changes in SST, pH, salinity and oxygen concentration in the open ocean regions and to mainly NPP and Sea Bottom Temperature in the coastal biome (Supplementary Figure 2). Notably, we projected the emergence of a large area of NABS in tropical regions by the end of the 21st century under RCP 8.5 (Figure 1).
We projected significantly higher rates of poleward centroid shifts and changes in total area of the BGCPs under RCP 8.5 compared to RCP 2.6 ( Figure 2). The average projected rate of poleward shift across the 56 BGCPs was 0.36 ± 0.66 km.year −1 (average and standard error) and 1.83 ± 1.99 km.year −1 under RCP 2.6 and RCP 8.5, respectively. ARAB (Arabian sea, 10.99 ± 0.81 km.year −1 ), NADR (North Atlantic Drift Region, 6.28 ± 4.72 km.year −1 ), and SATL (South Atlantic gyral province, 5.63 ± 3.45 km.year −1 ) were the BGCPs exhibiting the strongest shift in their centroid (Figure 2). The ARAB shift would translate into a move north up to a 1,180 km by 2100. In contrast, the ISSG (Indian ocean South Subtropical Gyre province)and BRAZ (Brazilian coastal province) BGCPs did not shift, keeping their centroid locations constant (rate of shift of 0.02-0.05 km.year −1 ). Biogeochemical provinces that were projected to expand the most were SARC (Sub atlantic ARCtic province; 1.98%.year −1 ), CNRY (1.52%.year −1 ) and NECS (NorthEast Atlantic Continental Shelve province, 1.08%.year −1 ) -effectively doubling or more in size by 2100 -while PEQD (Pacific EQuatorial Divergence province, -1.10%.year −1 ), ARCT (ARCTic province, -11.09%.year −1 ) and NADR (North Atlantic Drift Region, -11.09%.year −1 ) were expected to shrink substantially.  3). The probability of occurrence of the NADR province was projected to decrease sharply at the southern portion of its distribution towards the end of the 21st century under RCP 8.5, with the area between the Greenland Sea and Norwegian Sea displaying more NADR-like conditions over the course of this period. Hence, the overall location of this BGCP is projected to shift northward. For NATR (North Atlantic Tropical Region), the BGCP boundary is projected to expand northward and southward into the Atlantic Ocean, while the probability of occurrence of the NATR province shows only a slight decline by the end of the century under RCP 8.5. For the CHIL BGCP, the projected boundary and probability of occurrence remained relatively stable under all scenarios.
We identified the emergence of a NABS in large areas of the western tropical Pacific Ocean, eastern tropical Indian Ocean and eastern tropical Atlantic Ocean (Figure 4). The earliest emergence of a NABS region (across all 3 ESMs) was projected to occur in the tropical south Pacific Ocean in the 2040s, under RCP 8.5. This NABS region was projected to grow very rapidly during the 2050s and 2060s, extending into the eastern Indian Ocean and emerging off the coast of tropical west Africa in the 2070s and in the Caribbean in the 2080s. In contrast, under RCP 2.6, the emergence of NABS occurs on average one decade later than under RCP 8.5 in the same general areas. NABS projected size by the end of the 21st century is on average 11.9% of projections under RCP 8.5. NABS regions are characterized by high SST, low oxygen concentration and low NPP or low pH and low NPP conditions compare to the reference period (Supplementary Figure 3).
We projected that a substantial proportion of global marine diversity and fisheries catches are presently located in regions where NABS are projected to develop. Under RCP 8.5, 18.97 ± 6.05% in 2050 and 59.49 ± 4,08% in 2100 of the total number of exploited marine species considered here (N = 1,105 species) will be exposed to the expansion of NABS regions. These NABS regions are indeed projected to develop in areas with high seafood dependence [in terms of calories and nutrition (Golden et al., 2016)] with relatively low adaptive capacity (Blasiak et al., 2017) and high exposure to simultaneous losses in terrestrial food production capacity. The proportion of global marine fisheries catches (average between 2001-2011) caught in areas of predicted NABS in 2050 and 2100 are 7.87 ± 3.94% and 30.39 ± 5.32%, respectively, under RCP 8.5. The projected exposure to NABS was much lower under RCP 2.6, with only 15% of exploited species and 5% of total catch in NABS areas by the end of the 21st century. The difference in trajectories in regards to NABS between RCP 2.6 and RCP 8.5 became significantly larger than the inter-ESM variability after the end of the 2040s.

DISCUSSION
This study supports the hypotheses that climate change will deeply affect the biogeography of the global ocean, and lead to the emergence of biomes with no historical analog. At the global scale, we predict that trade wind provinces will expand (Staten et al., 2018), with the distribution of several westerly wind provinces migrating poleward. This evolution will be driven mostly by the poleward extension of warm surface waters, pH, oxygen concentration, and NPP (Supplementary Figure 2). Warming of the tropics will drive a poleward migration of the BGCPs (Figure 1 and Supplementary Figure 2) with the polar provinces progressively shrinking over time and being concentrated at the highest latitudes. These biogeographical changes appear to be corroborated by the rapid modifications in biogeochemical processes, species composition and food web dynamics already documented for these regions (Beaugrand et al., 2002;Dulvy et al., 2008;Polovina et al., 2008;Stramma et al., 2012;Fossheim et al., 2015). Such patterns are also in accordance with other projections performed using modeled distributions of marine species (Cheung et al., 2009;Lenoir et al., 2011), biomes (Sarmiento et al., 2004) and climate velocity computed using only SST (Burrows et al., 2011). However, in areas where province boundaries are constrained by the presence of land, we expect a significant decline in the original area of several provinces (ARAB, GUIN, GUIA, CARB, CHIN), with the potential to significantly impact the ecological assemblages characterizing those provinces.
Our results show that NABS regions will cover areas where a substantial proportion of global marine biodiversity presently occurs and with a crucial dependence on seafood production. The emergence of wide NABS areas in the tropical ocean will exacerbate the high vulnerability of populations living in developing coastal nations and small island states (Lam et al., 2016;Blasiak et al., 2017). Biodiversity, structure and function of marine ecosystems, and fisheries catches are closely related to the characteristics of the BGCPs (Reygondeau et al., 2013). The climate-induced changes to BGCPs that we project provide additional support to the idea that climate change will deeply alter the distribution and function of marine ecosystems as well as the benefits derived from them by people. Predicting the evolution of marine biodiversity in the NABS areas raises significant challenges, as the knowledge we currently hold regarding the organization of ecosystems in existing biomes is likely to be invalidated in no-analog environmental conditions (Fitzpatrick and Hargrove, 2009). While we cannot definitively state that no species stand to benefit from these new conditions, as the NABS predominantly occur in tropical regions where many species are already close to physiological maxima (Walther et al., 2002;Sunday et al., 2011Sunday et al., , 2012, it is likely the emergence of NABS will substantially elevate the risk of impacts on marine biodiversity and fisheries. The NABS were characterized by very warm mean annual temperatures, high salinity, low oxygen concentration and low NPP compared to the reference period (Supplementary Figure 3). Most marine species will be physiologically stressed under such conditions, which could impact their survival rate (Sunday et al., 2012). This explains the high rate of local extinction predicted by previous studies in the regions where NABS environmental conditions will not fit the environmental tolerance range of the majority of endemic species (Beaugrand et al., 2015;Jones and Cheung, 2015). Paleontological records have shown that the apparition of novel types of climate increases the rate of migration impacting the local biodiversity pool (Jablonski, 2008). Also, such changes in the local environmental condition will impact the biological processes of marine species. In NABS regions, we therefore expect a decrease in the size of endemic species as well as a decline in local biomass, as suggested by the literature. However, it is possible that some species will manage to survive in NABS-like environments, particularly small organisms with high turnover rates, such as microbes, that have good potential for rapid adaptation, as well as species that are physiologically more flexible, such as widely distributed marine species that occur in the Pacific Warm Pool. However, species that prove capable of adapting in NABS may only represent a small fraction of the currently high biological diversity occurring in these regions (Donelson et al., 2019;Palumbi et al., 2019).
Our results can be linked to recent studies on marine extreme events (Frölicher et al., 2018;Oliver et al., 2018) that have shown, using satellite observations, that recurring and permanent marine heatwaves already occur in the regions where we expect NABS to appear by the middle of the 21st century. We can reasonably expect that the impacts that these marine heatwaves already have had on marine life (Le Nohaïc et al., 2017;Oliver et al., 2017) will be similar or even amplified under NABS conditions. The emergence of NABS would be limited in most regions of the global ocean if CO 2 emissions can be held to the level required to achieve RCP 2.6 (Figure 1). However, if this is not the case (which currently appears likely), we expect the earliest NABS to emerge in about twenty years from 2010 in the tropical Pacific region and a decade later in the Indo-Pacific basin. There is therefore an urgent need to develop adaptation policies anticipating the risk of rapid ecosystem collapses in these regions, in addition to supporting efforts for strong mitigation.
The original methodology developed in Reygondeau et al. (2013) was applied on observed environmental data and used the ocean partition proposed by Longhurst (2007) as framework to quantify the environmental envelopes of BGCPs. Consequently, the application of the methodology to ESM have large uncertainty in correctly representing the ocean properties in several regions of the Ocean and hence distribution of BGCPs. The natural provinces (or regions of similar environmental conditions) derived from the models may not spatially match the prescribed boundaries found in Reygondeau et al. (2013). We subsequently tested the spatial and temporal agreement between our predictions based on different ESM projections and results from Reygondeau et al. (2013) (Supplementary Figure 1). Overall, the spatial division and internal variation of BGCPs are well captured in ESM outputs (Supplementary Figure 1), but several regions need to be considered with caution. First, the dynamic of coastal regions are not well captured by this generation of ESMs in the CMIP5 repository (Stock et al., 2011;Cheung et al., 2016a) and consequently, the distribution and dynamics of coastal BGCPs are not well represented for current or future projection. Second, semi enclosed sea (Red sea, Persian Gulf, Baltic sea, Mediterranean sea) have similar biases as coastal regions and results for these areas need to be taken with caution. Third, while the division of the southern ocean is well captured by the methodology using ESM outputs, several biases in the representation of the region biogeochemical dynamics are known (Bindoff et al., 2019). Consequently, future projection of the BGCPs in the southern Ocean are sensitive to thse biases. We hope that the finer resolution of the new generation of ESM models that are under development may resolve these issues in a future re-analysis.
Our results indicate that the environmental changes that would occur in the global ocean along a "no mitigation" RCP 8.5 scenario would lead to a drastic reorganization of global marine biogeography, associated biodiversity and trophic networks. These changes would include the emergence of wide regions with no environmental analog compared to current observations. If the global climate is not kept below 2 • C warming, NABS areas can be expected to emerge, as early as 20 years from the 2010s. It would affect 19% of the total number of exploited species in 2050 and 59% in 2100 and would cover regions that are currently responsible for 8% of global marine fisheries catch in 2050 and 30% in 2100, under RCP 8.5. These numbers would change to only 15% of exploited species and 5% of total fisheries catches in NABS areas by the end of the 21st century under the RCP 2.6 scenario. Mitigating anthropogenic pressures at a level sufficient to reach the Paris agreement targets would therefore substantially reduce the risk of emergence of large NABS regions in the global ocean, and the dramatic consequences that such large-scale ecological changes would entail for tropical marine biodiversity, associated fisheries and the human communities that they support.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
GR performed all the analyses and wrote the original of the manuscript. Idea originated from OM and GR. All other coauthors contributed to data sharing, idea, analysis, and writing of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
WC, GR, CW, and VL acknowledge funding support from the Nippon Foundation-the University of British Columbia Nereus Program. WC acknowledges support from the Natural Sciences and Engineering Research Council of Canada. OM acknowledges the support of the French ANR project CIGOEF (grant ANR-17-CE32-0008-01). TF thanks the Swiss National Science Foundation (PP00P2_170687) and the European Union's Horizon 2020 Research and Innovation Programme under grant agreement no. 820989 (project COMFORT, Our common future ocean in the Earth system -quantifying coupled cycles of carbon, oxygen, and nutrients for determining and achieving safe operating spaces with respect to tipping points) for financial support.

ACKNOWLEDGMENTS
This study is dedicated by GR to the memory of Madeleine Remy.