Emergence of climate change signals in marine ecosystem thermal niches

Temperature is one of the most important drivers of global ocean patterns of biodiversity 1,2,3 , shaping thermal niches through thresholds of physiological thermal tolerance 4 . Because of anthropogenic global warming, lower and upper thermal niche bounds are predicted to change affecting the future distribution of marine species 5,6 . Current working hypotheses suggest an expansion of ectotherms toward their poleward boundaries 7,8 . Nonetheless, the knowledge of the timing and extent of these rearrangements across latitude and depth remains limited. Here, using daily data across the water column from both Ocean Sites network observations and novel Earth System Model, we track the emergence of thermal niches whose lower bound is warmer than their current upper bound, potentially disrupting marine habitats. We show that these developments will emerge by ~2030 in subsurface waters (~50 – 1000 m) if anthropogenic emissions continue to rise, whereas they delay several decades if emissions are substantially reduced. By 2100, thermal niches will be warmer than current counterparts. However, we further show that depending on the vertical level, concomitant changes in both boundaries will result in wider or narrower thermal niches. These results suggest that the redistribution of marine species might differ across depth, shedding light upon a much more complex picture of the impact of climate change on marine habitats. Anthropogenic climate change fingerprints in the world's oceans by warming significantly their upper layers 9,10 . This increase in heat content is projected to alter long-term well defined thermal niches driving species redistribution at a global scale 11,12 . These changes are already affecting goods and services provided by the oceans 13 , and are projected to be amplified with rising greenhouse gases emissions. A close coupling between marine organisms' physiological thermal tolerances and environmental temperature 14,15 suggests that distributional shifts can be predicted by tracking changes on the lower and upper bounds of the ecosystems’ thermal niches 6 . Few studies have yet considered the vertical structure of temperature in the ocean (e.g., ref. 16), thus current understanding on how and when climate change

just the tip of the iceberg of an ocean at change under anthropogenic pressure, and the general extent to which species distribution will be affected by changes in ecosystems' thermal niches below the surface remains an open question.
Predicting changes in the distribution of marine ecosystems due to variations in environmental temperature relies on the assumption that species' tolerance ranges reflect the magnitude of the local temperature variability 15,19 . Marine organisms tend to occupy the extension at which environmental temperatures lie within their tolerances 20 , perhaps owing to oxygen limitation of these tolerance ranges 21 . We thus consider the lower and upper limits of the thermal niches to be represented by the environmental temperature minimum (Tmin) and maximum (Tmax). At each depth, Tmin and Tmax correspond to the annual 1 st and 99 th percentile as computed using the statistical distribution of daily records, respectively. Though this approach encompasses most of the range of temperature variability, it can yield conservative projections as tolerance ranges can be wider than environmental thermal niches 22,23 , and as ectotherms display some plasticity to adapt to environmental temperatures that challenge their tolerance limits 24,25 .
Nonetheless, there is limited capacity of acclimation when long-term heating occurs 26 , especially for tropical species 27 and during reproductive stages 28 . In this regard, our overarching objective is to understand how and when global warming-induced changes in thermal niches over the water column will take place in the future affecting current marine ecosystems.
We take advantage of comprehensive data sets of daily three-dimensional ocean temperatures from both the long-term Ocean Sites (OS) network and a state-of-the-art Earth System Model (ESM). This novel data set and its exploitation represents a major novelty to improve our understanding of how future climate change will impact marine ecosystems. To provide a reliable framework upon which to derive the future evolution of lower and upper thermal niche bounds, we select six OS stations for which at least seven years of daily temperature observations from the surface to ~1000 m depth are available (Table S1). We map them into polar, temperate, and tropical ocean domains, and determine the surface area of ocean domains informed by the OS stations by computing the level of similarity in daily temperature profiles using a p-value analysis (Fig. 1a, see Methods). A pattern of alternation of cooling and warming periods is seen over the time of available observations (Fig. 1b). In general, these episodes last a few years though they last longer at the northernmost stations. Though measurements' coverage is not complete along depth and time for some stations, we consider  us to confidently compute annual Tmin and Tmax, and extract trends to compare with   the ESM simulation (see Supplementary Material).   Three-members ensemble simulations were performed with CNRM-ESM2-1 29 (see Methods) encompassing the historical period (1850-2014) followed by two future projections (2015-2100) that explore very contrasted emission pathways 30

Emergence of changes in current thermal niches
With the shape of current thermal niches established, we track when and where concomitant changes in the lower and upper boundaries of the thermal range will take place. Concomitant changes in Tmin and Tmax can be seen as a compound event 34 since they can result in several developments of the thermal niches, possibly leading to a redistribution of marine ecosystems or driving to their collapses (Fig. S2). Thermal niches will shift toward warming or cooling temperatures if both the temperature of the lower and upper bounds of the thermal range increases or decreases, respectively. In addition, thermal niches can also expand if Tmax warms more rapidly than Tmin, or narrow if Tmin warms more rapidly than Tmax. As depicted by the observations and as simulated by CNRM-ESM2-1, trends in Tmin and Tmax over the recent years ( Fig. S3 and S4) lead to various developments in thermal niches. In general, warming trends were stronger than cooling trends thus resulting in warmer thermal niches. The pace at which these changes occurred varies between stations, but it has been shown to push marine organisms to acclimate to a warmer environment 11,35,36 (see Supplementary Material).
Overlapping to this warming, an imbalanced warming of Tmax results in wider thermal niches.
These expanded thermal niches may host more marine species generating interspecific competition for resources and space. By contrast, excess warming of Tmin shrinks the thermal range, possibly leading to the loss of ecosystems' habitability. Though the responses of marine organisms to these changes can be difficult to predict 37 , as global warming trends are likely to increase 16,38,39 , it is key to understand the time at which these changes may occur.
To assess how and when future thermal niches may emerge from warming-induced changes in the thermal range bounds, we modify the canonical approach of the Time of Tmin surpasses a first threshold (Tmidpoint), we consider that the shift of the thermal range may represent a warning to current ecosystems since Tmidpoint has been observed to align well with the temperature of maximum ecological success of a specific species (see ref. 5). Furthermore, when an ecosystem will be exposed to a Tmin that is warmer than the current Tmax, we consider that organisms should deal with a completely new thermal niche (see Fig. 3a). Considering the approach discussed above, we find a rather good agreement across all OS stations in the first 200 m of the water column (Fig. 3b), with only noticeable changes in the thermal range within lower epipelagic waters (50 -200 m). This feature results from both the shape of the current thermal ranges (Fig. 2), and the rather homogeneous warming of the ocean in these layers (see

End-of-the-century thermal niches
We now focus on the changes in thermal ranges by the end of the century (2080 -2100). This analysis is relevant for tracking the impact of climate change on marine ecosystems as the lower and the upper boundaries of the thermal range will continue to change even after crossing ecosystems thresholds. As changes in thermal range boundaries can take place in both directions and at different pace because of physical and dynamical ocean processes, they can result in complex rearrangements of the thermal niches by the end of the century.
Under the high emission scenario, end-of-the-century thermal niches differ from those estimated over the historical period (1990 -2014) ( Fig. 4 and S6). In general, both the lower and upper  Consequently to these changes, complex thermal niche developments take place (Fig. 4).
Depending on the station and the layer considered, end-of-the-century thermal niches will be warmer and narrower as a result of a quicker warming of Tmin (yellow rectangles), or warmer and wider as Tmax warms more rapidly than Tmin (violet rectangles). In the former case, new conditions will challenge local adaptation of inhabitant organisms. In the latter, possible spread of species from neighbour habitats may generate additional stresses by changing species interaction 12 . Wider thermal niches will result above 200 m at all stations under a high-emission scenario with the only exception at station CIS-1 due to an excess warming of Tmin (Fig. S6a).
Below 200 m, the pace of warming of both bounds are comparable, generating both wide or narrow thermal niches depending on the station. Considering a low emission scenario, more heterogeneous developments appear. All stations show warming anomalies (Fig. S6b), except at station CIS-1 where both bounds will generally be cooler than the historical mean. Only developments at stations FRAM and K276 are consistent across emission scenarios.

Concluding remarks
Marine ecosystems are already being threatened by numbers of anthropogenic stressors like fishing 41 , acoustic pollution 42 or plastics 43  immediate emission reduction consistent with a high mitigation scenario. In response to ocean warming, the shape of thermal niches will also change. We project that by 2100, depending on the vertical layer, thermal niches will be warmer and wider, potentially increasing the interspecific pressure; or warmer and narrower, forcing organisms to adapt to new local conditions or migrate. Assuming organisms are adapted to current environmental conditions, such changes may lead to rearrangements of marine habitats in the decades to come across latitude and, if possible, depth 47 . These results shed light on a much more complex picture (see

Ocean Sites observations
The Ocean Sites (OS) network constitutes a worldwide effort to monitor ocean parameters through high-quality data extracted from long-term, high-frequency observations at several locations of the World ocean. Six OS stations, listed in Table S1, were selected because of the availability of continuous daily measurements of ocean temperature and salinity across the water column for more than seven years, allowing a robust computation of thermal range boundaries (see below). All of the six stations provide data from the surface to about 1000 m depth, that have been resampled daily at each depth, and then interpolated into the vertical grid of CNRM-ESM2-1.

Simulations
This work exploits simulations from a state-of-the-art Earth system model, CNRM-ESM2-1 29 , that has been developed by the CNRM-CERFACS climate group for the sixth phase of the In this study, we performed four simulation with CNRM-ESM2-1: a 250 year-long pre-industrial control simulation (without anthropogenic forcing) to estimate the model's internal variability; and a historical simulation from 1850 to 2014 followed by two future scenarios from 2015 to 2100, which are used to derive present and future variations in temperature minimum and maximum.
These simulations were produced using the external forcing as recommended by CMIP6 for the pre-industrial state and the historical period. For the future scenarios we used contrasting pathways: a low (SSP1-2.6) and high (SSP5-8.5) emission pathways as described in ref. 30.
All simulations provide daily outputs from the ocean surface to 4000 m for ocean temperature and salinity as well as oxygen, pH and net primary productivity. Here, we exploit only the first 47 vertical layers for ocean temperature and salinity in order to describe the first 1000 m depth of the water column. Finally, to ensure both observations and model data (historical + SSP5-8.5 simulation) cover exactly the same period, model data was selected to begin and end at the same date as observations.

Ocean domains informed by Ocean Site stations
Though OS networks are located throughout the World ocean, our selection of OS stations is disproportionately located in the North Atlantic and North Pacific oceans (see Fig. 1 and Table   S1). To assess how large is the surface area of ocean domains informed by our six selected OS stations, we compute the level of similarity between daily profiles as provided by observations and the model hindcast over the current period (1990-2020) using the statistical approach presented in ref. 50. This approach compares simultaneously the mean and the daily variations of OS daily profiles with a neighbour grid-point model profile using a Chi-squared-based test.
The test consists in comparing the cumulative sum of the Welch's tz 2 51 across depth levels to an empirical Chi-squared distribution with 47 degrees of freedom (i.e., the number of depth levels).
We use 10,000 random samples of this Chi-squared distribution to estimate the empirical distribution of the Chi-squared law. The distribution is then used to compute an empirical 'integrated' p-value that represents an objective metric to determine how far the two profiles are consistent between each other within the depth interval.
The empirical 'integrated' p-value allows us to quantify the match between profiles. We establish a threshold of 0.90 to consider a profile over a grid-cell consistent with the OS profile. For further analysis, stations were grouped into three ocean domains: polar, temperate and tropical waters.

Estimation of thermal niche boundaries
The working definition of the thermal niche employed in this work assumes that the species' tolerance ranges reflect the magnitude of the local temperature variability 15,19 . As a consequence, we infer the vertical structure of the thermal niche from the lower and upper limits of the thermal range that is captured by the minimal and maximal environmental temperature across the water column.
Thanks to high-frequency data, we provide a robust yearly estimate of these bounds using the  For illustration purposes, we illustrate how our approach works for four depth levels of HOT-01 station (Fig. S5).

End-of-the-century thermal niches
End-of-the-century thermal niches provide a snapshot of the concomitant changes in thermal range boundaries resulting from climate change. We compute the end-of-the-century thermal niche from daily data over the 2080-2100 period. Fig. 4 displays the end-of-the-century thermal range anomalies of both Tmin and Tmax with respect to the mean over 1990-2014, corresponding to the last years of the historical simulation. To compare with the historical profiles, we also include their anomalies. At each depth level, we assess the magnitude of the changes between the end-of-the-century and the historical profiles.
As changes in thermal range boundaries can evolve in both directions, and with a different pace, they may result in a re-arrangement of the vertical shape of the thermal range. To track if end-of-the-century thermal niches are also wider or narrower, we compute the difference between Tmax and Tmin anomalies at each depth level (Fig. S6). If the difference is positive, thermal niches will be wider, i.e., Tmax warms more rapidly. If the difference is negative, thermal niches will be narrower, i.e., Tmin warms more rapidly. If differences are < 0.05 ºC (i.e., level of uncertainty informed from the analysis of the internal variability of thermal range profiles in Fig. S4), we consider no changes in the shape of thermal niches will take place, i.e., only shifting toward warming or cooling is projected. Both emission pathways are displayed in Fig.   S6.
In Fig. 4 we have assigned a colour code for each of these developments. A depiction of these developments is provided in Fig. S2.

Model internal variability
As a consequence of the chaotic nature of processes in the Earth systems being simulated (ocean-atmosphere-land-biosphere-cryosphere), one of the main sources of uncertainties in climatic future projections is their internal variability 52 . One way to isolate these uncertainties is