Alejandro Ariza
email: alejandro.ariza@ird.fr

Matthieu Lengaigne

Christophe Menkes

Anne Lebourges-Dhaussy

Aurore Receveur

Thomas Gorgues

Jérémie Habasque

Mariano Gutiérrez

Olivier Maury

Arnaud Bertrand

Global decline of pelagic fauna in a warmer ocean

Pelagic fauna is expected to be impacted under climate change according to ecosystem simulations. However, the direction and magnitude of the impact is still uncertain and still not corroborated by observational-based statistical studies. Here we compile a global underwater sonar database and 20 ocean climate projections to predict the future distribution of sound-scattering fauna around the world's oceans. We show that global pelagic fauna will be seriously compromised by the end of the 21 st century, if we continue under the current greenhouse emission scenario. Low and mid latitudes are expected to lose from 3 to 22% of animal biomass due to the expansion of low-productive systems, while higher latitudes would be populated by present-day temperate fauna, supporting conclusions drawn from ecosystem simulations. We further show that strong-mitigation measures to contain global warming below 2°C would reduce these impacts to less than half.

The pelagic zone-the vast volume of open-ocean water from the surface to the seabed-is the largest living space of the planet. This three-dimensional environment holds half of the global primary production and sustains most of the animal biomass on earth [START_REF] Field | Primary Production of the Biosphere: Integrating Terrestrial and Oceanic Components[END_REF][START_REF] Bar-On | The biomass distribution on Earth[END_REF] . Dominated by small crustaceans, fish, molluscs and a variety of gelatinous life-forms, this mid-trophic level community channels primary production to the top of the food web [START_REF] Choy | Finding the way to the top: how the composition of oceanic mid-trophic micronekton groups determines apex predator biomass in the central North Pacific[END_REF] , being the fundamental forage base for apex predators such as marine mammals, seabirds, and major commercial fish stocks [START_REF] Pauly | Primary production required to sustain global fisheries[END_REF][START_REF] Bertrand | Broad impacts of fine-scale dynamics on seascape structure from zooplankton to seabirds[END_REF] . In addition, a substantial proportion of mid-trophic level fauna is involved in diel vertical migration [START_REF] Brierley | Diel vertical migration[END_REF] , the largest movement of animals on earth [START_REF] Behrenfeld | Global satellite-observed daily vertical migrations of ocean animals[END_REF] , and a fundamental pathway of biomass exchange between epipelagic (0-200 m) and mesopelagic (200-1000 m) waters [START_REF] Angel | Vertical Distribution of the Standing Crop of Plankton and Micronekton at Three Stations in the Northeast Atlantic[END_REF][START_REF] Cook | Deep-pelagic (0-3000m) fish assemblage structure over the Mid-Atlantic Ridge in the area of the Charlie-Gibbs Fracture Zone[END_REF] . By respiring, defecating, and occasionally being predated at depth, vertical migrants inject large amounts of carbon into the ocean interior, contributing to the sequestration of atmospheric CO2 [START_REF] Hidaka | Downward transport of organic carbon by diel migratory micronekton in the western equatorial Paci"c: its quantitative and qualitative importance[END_REF][START_REF] Ariza | Migrant biomass and respiratory carbon flux by zooplankton and micronekton in the subtropical northeast Atlantic Ocean (Canary Islands)[END_REF][START_REF] Saba | Toward a better understanding of fish-based contribution to ocean carbon flux[END_REF] .

Projections of ocean warming, deoxygenation and primary production declines in response to climate change [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF][START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] have led to concerns about the future evolution of mid-trophic level fauna. This fostered the emergence of global marine ecosystem models that include this critical component in their simulations [START_REF] Tittensor | A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0[END_REF] .

These models, forced by future physical and biogeochemical oceanic model outputs, are used to simulate marine ecosystems, from primary producers to top predators. Overall, ecosystem models project a global decline of animal biomass by the end of the 21 st century, primarily driven by increasing temperature and decreasing primary production [START_REF] Bryndum-Buchholz | Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins[END_REF][START_REF] Kwiatkowski | Consistent trophic amplification of marine biomass declines under climate change[END_REF][START_REF] Lotze | Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change[END_REF][START_REF] Tittensor | Next-generation ensemble projections reveal higher climate risks for marine ecosystems[END_REF] . However, because of their complex design and diverse representation of ecological processes [START_REF] Tittensor | A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0[END_REF] , these models still project a wide range of directional changes at regional scale, even within the same climate simulation experiment [START_REF] Heneghan | Disentangling diverse responses to climate change among global marine ecosystem models[END_REF] . Most importantly, ecosystems models lack ground truth validation of mid-trophic level fauna at global-scale.

Primary producers and other plankton components can be remotely sensed from space or continuously sampled from moorings and ships of opportunity [START_REF] Reid | Mesopelagic-boundary community in Hawaii: Micronekton at the interface between neritic and oceanic ecosystems[END_REF][START_REF] Ben Mustapha | Automatic classification of waterleaving radiance anomalies from global SeaWiFS imagery: Application to the detection of phytoplankton groups in open ocean waters[END_REF] , however, research on mid-trophic level fauna requires large sampling nets and complex deep water operations [START_REF] Pakhomov | Report of the Advisory Panel on Micronekton Sampling Intercalibration Experiment[END_REF][START_REF] Kaartvedt | Efficient trawl avoidance by mesopelagic fishes causes large underestimation of their biomass[END_REF] . As a result, knowledge on its distribution, diversity, and response to environmental factors has long been limited to regional research or to global reviews based on variable methodological approaches [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Catul | A review on mesopelagic fishes belonging to family Myctophidae[END_REF] . Observations with ship-borne acoustic echosounders are an alternative approach that, while not taxonomically precise, provide consistent and large-scale data on predominantly mid-trophic level fauna from pelagic ecosystems [START_REF] Benoit-Bird | Ecological Insights from Pelagic Habitats Acquired Using Active Acoustic Techniques[END_REF][START_REF] Annasawmy | Micronekton diel migration, community composition and trophic position within two biogeochemical provinces of the South West Indian Ocean: Insight from acoustics and stable isotopes[END_REF] . These observations are increasingly available over large ocean regions [START_REF] Haris | Sounding out life in the deep using acoustic data from ships of opportunity[END_REF][START_REF] Irigoien | The Simrad EK60 echosounder dataset from the Malaspina circumnavigation[END_REF] and have recently permitted to revise estimates of global mesopelagic fish biomass [START_REF] Irigoien | Large mesopelagic fishes biomass and trophic efficiency in the open ocean[END_REF] as well as to describe sound-scattering and migration in major ocean basins [START_REF] Klevjer | Large scale patterns in vertical distribution and behaviour of mesopelagic scattering layers[END_REF][START_REF] Proud | Fine-scale depth structure of pelagic communities throughout the global ocean based on acoustic sound scattering layers[END_REF] . Linear statistical modelling between environmental drivers and the intensity of sound-scattering layers also allowed for an assessment of future mesopelagic fish biomass in response to climate change [START_REF] Proud | Biogeography of the Global Ocean's Mesopelagic Zone[END_REF] .

This study indicates that ocean warming and decreased primary production is expected to cause an increase of fish biomass by the end of the 21 st century [START_REF] Proud | Biogeography of the Global Ocean's Mesopelagic Zone[END_REF] , in contradiction with projections from all ecosystem models Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 forced under comparable climate scenarios [START_REF] Bryndum-Buchholz | Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins[END_REF][START_REF] Kwiatkowski | Consistent trophic amplification of marine biomass declines under climate change[END_REF][START_REF] Lotze | Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change[END_REF][START_REF] Tittensor | Next-generation ensemble projections reveal higher climate risks for marine ecosystems[END_REF] . Such discrepancies between statistical and mechanistic modelling illustrate that there is still no consensus on the fate of pelagic fauna under multiple global climate stressors.

To shed light on this issue, we constructed a global atlas of sound-scattering mid-trophic level fauna, based on the largest acoustic dataset available to date. This database includes fundamental oceanic biomes that were previously unrepresented, such as coastal upwelling, oxygen-minimum, and equatorial zones. The structure and distribution of mid-trophic level fauna is projected worldwide for the present period, and under different future climate change scenarios, based on the latest generations of Earth system models from the Coupled Model Intercomparison Project Phase 5 and 6 (CMIP5 [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF] and CMIP6 [START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] ). These projections relied on twenty different climate simulations, an optimized suite of environmental predictors, and nonlinear statistical modelling based on machine learning. With this design, we project strong declines of midtrophic level fauna in a warmer and less productive ocean by the end of the 21 st century, corroborating ecosystem model projections with an observationally-based approach.

Acoustic seascape classification. We compiled a database with publicly available acoustic data at 38 kHz, the most widely used frequency for observing mid-trophic level fauna at epipelagic and mesopelagic depths. The database covered more than 350,000 km around the globe, extended from 20 to 750 m depth, and spanned from year 2001 to 2020 (Supplementary Table 1). The ocean was gridded into 4x4 degree cells and average day and night profiles were computed for each cell (Extended Data Fig. 1). In total, we obtained 459 cells, which represented 16% of the ocean surface, from open to coastal-boundary systems and from 66ºS to 50ºN. We transformed the vertically discretized data of profiles into continuous curves and explored the variation of these curves-day and night together-using functional principal component analysis [START_REF] Ramsay | Functional data analysis[END_REF] . The principal components that accounted for more than 90% of the total backscatter variance of profiles were used to classify the data through hierarchical clustering (see Methods and Extended Data Figs. 2 and3).

We identified six acoustic seascape classes associated with subpolar, gyre, subtropical, temperate, upwelling, and low-oxygen upwelling waters (Fig. 1a). Vertically-integrated backscatter of acoustic profiles-a proxy of pelagic biomass-was minimum in subpolar and oligotrophic-gyre waters, moderate in subtropical waters, and higher in rich temperate and upwelling waters, including those in oxygen minimum zones (Fig. 1b-m). With the exception of subpolar waters, this distribution was consistent with patterns of pelagic biomass worldwide, according to net-sampling studies [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Moriarty | Distribution of mesozooplankton biomass in the global ocean[END_REF] . The most distinctive feature in the classification was the shape of acoustic profiles, which illustrated the vertical structure of marine life across oceanic systems. While subpolar profiles were characterized by a smooth depth-increasing backscatter signal (Fig. 1b,c), those in rich upwelling and low-oxygen areas presented prominent peaks in Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 epipelagic and mesopelagic waters (Fig. 1j-m). These peaks denoted the presence of strong soundscattering layers, which were also present-albeit with lower intensity and deeper depths-in gyre, subtropical and temperate systems (Fig. 1d-i). The classification also showed the main factors modulating the distribution of sound-scattering layers at global-scale. First, increased light penetration in oligotrophic ocean-gyre areas push mesopelagic fauna to deep waters [START_REF] Aksnes | Light penetration structures the deep acoustic scattering layers in the global ocean[END_REF] , while in productive systems, the presence of subsurface oxygen minimum zones constrains organisms to thinner layers near the surface [START_REF] Bertrand | Acoustic Observation of Living Organisms Reveals the Upper Limit of the Oxygen Minimum Zone[END_REF][START_REF] Bianchi | Intensification of openocean oxygen depletion by vertically migrating animals[END_REF] . Second, the intensity of these animal layers is expected to be higher in productive oceanic areas such as temperate, upwelling, and low-oxygen-upwelling systems [START_REF] Irigoien | Large mesopelagic fishes biomass and trophic efficiency in the open ocean[END_REF] , as shown by our classification. The differences between day and night profiles-in particular the exchange of backscatter maxima between epipelagic and mesopelagic waters-illustrates that diel vertical migration is more intense in productive areas. In principle, the integrated backscatter of profiles was expected to be higher at night due to the incorporation of migrants from waters below our sampling range, however, that was not the case for upwelling and low-oxygen regions (Fig. 1j-m). While vertical migration itself is the main source of day and night differences in watercolumn backscatter, changes in the swimbladder condition of fish during these excursions are also known to have important effects on how sound is reflected back to surface [START_REF] Godø | Diel migration and swimbladder resonance of small fish: some implications for analyses of multifrequency echo data[END_REF] Present-day global seascape. We tested ten supervised learning algorithms, as well as ten potential environmental predictors, to extend our seascape classification beyond the observation locations. Crossvalidation tests and recursive elimination of weak predictors revealed that the best performance was achieved with a random forest classifier trained with satellite-derived chlorophyll concentration, subsurface dissolved oxygen, and sea surface temperature (78% predictive accuracy, see model selection and performance assessment in Methods and Extended Data Figs. 456). We extended the classification with this model at global scale, using the mean environmental conditions for the period 2000-2020. Seascape classes were coherently distributed along large-scale oceanographic systems associated with subpolar, gyre, subtropical, temperate, upwelling, and low-oxygen upwelling waters (Fig. 2a). As this partitioning relied on the combined day and night structure of acoustic profiles, their boundaries reflected not only differences in water-column biomass but in qualitative community properties such as the vertical distribution or the migration patterns of pelagic fauna. Indeed, net-based studies on the distribution of mesopelagic fishlikely the major contributor to backscatter [START_REF] Agersted | Mass estimates of individual gas-bearing mesopelagic fish from in situ wideband acoustic measurements ground-truthed by biological net sampling[END_REF] -show comparable community transitions between oceangyre, equatorial, coastal-boundary, and mid-latitude temperate systems [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Catul | A review on mesopelagic fishes belonging to family Myctophidae[END_REF][START_REF] Backus | Pelagic faunal provinces and sound-scattering levels in the Atlantic Ocean[END_REF] broader pelagic fauna inventories also identified fairly similar oceanographic biomes [START_REF] Longhurst | Ecological Geography of the Sea[END_REF][START_REF] Spalding | Pelagic provinces of the world: A biogeographic classification of the world's surface pelagic waters[END_REF][START_REF] Sutton | A global biogeographic classification of the mesopelagic zone[END_REF] . We therefore consider these acoustic seascape classes as "echobiomes", as they represent sound-scattering communities with comparable structural and functional properties but not necessarily geographically connected or sharing the same species composition [START_REF] Spalding | Pelagic provinces of the world: A biogeographic classification of the world's surface pelagic waters[END_REF] . The vertical integral of the acoustic profiles predicted by our model can also be regarded as an indicator of pelagic biomass around the globe (Fig. 2b). Day and night biomass proxies presented similar spatial patterns, with minimum values in ocean-gyre and subpolar systems, and maximum in equatorial, temperate, and coastal-boundary systems (see night projections in Supplementary Fig. 1). Differences between low and high backscatter areas were however larger when using daytime data. These extended global projections were also consistent with large-scale biomass charts derived from net sampling [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Moriarty | Distribution of mesozooplankton biomass in the global ocean[END_REF] . The only exception occurred in subpolar echobiomes, where net based estimates were not as low in comparison to other regions [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Moriarty | Distribution of mesozooplankton biomass in the global ocean[END_REF] .

Future global seascape. Future changes in the distribution of echobiomes were then calculated by rerunning our model with average environmental conditions projected for 2080-2100. These were initially obtained from simulations of 13 different CMIP6 climate models forced by the "business-as-usual" highemission scenario SSP5-8.5 (see Methods). Our results point to a massive redistribution of the acoustic seascape by the end of the 21 st century, largely driven by the expansion of low-productive ocean-gyre and subtropical echobiomes towards current equatorial, temperate, and coastal-boundary rich systems (Fig. 2ce). Overall, this rearrangement of the seascape would modify 20% of the total ocean surface analysed in comparison to present-time values (Fig. 3e). The expansion of ocean-gyres against surrounding systems is expected to promote a generalized biomass decline between 40°S and 40°N, while the poleward migration of temperate systems would result in biomass gains at higher latitudes. These regional biomass changes could locally reach up to 40% (Fig. 2f). According to our projections, biomass changes will be significant and consistent across latitudes and ocean basins, for most climate forcings considered, and regardless of whether day or night backscatter is used as proxy of biomass (Fig. 3a). Considering the day backscatter biomass proxy, expected biomass losses from 40°S and 40°N would range between 3 and 22%, with major losses in the tropics. Biomass gains at subpolar latitudes would range between 9 and 16% in the northern hemisphere, and between 8 and 31% in the southern hemisphere (Fig. 3a). Despite the opposing changes between low and high latitudes, projected losses would far outweigh projected gains as the ocean's surface between 40°S and 40°N represent 70% of the total. The picture was qualitatively similar when the night backscatter proxy is considered, but the magnitude of biomass change would be lower (see Fig. 3a and Supplementary Tables 2 and3). Both day and night biomass changes were largely driven by mesopelagic fauna and, to a lesser extent, by migratory epipelagic fauna at night (see Supplementary Fig. 1). Further, the displacement of echobiomes is expected to affect the vertical distribution of biomass around the globe. The Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 most significant and consistent change would take place in tropical upwelling regions, where epipelagic and mesopelagic fauna is expected to deepen from 10 to 100 meters. This process would result from an expansion of the deep sound-scattering fauna, currently in warm low-productive regions, into tropical upwelling systems (see Supplementary Fig. 2).

Likely future. Mid-trophic level fauna projections presented above are based on changes in environmental conditions projected by the average of the latest generation of Earth system models, CMIP6, and considering the high carbon emission scenario SSP5-8.5. While CMIP6 is the current baseline for climatebased socioeconomic and ecosystem assessments [START_REF]Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] , comparative studies with the previous generation of climate models indicate that CMIP6 projects stronger warming but weaker and less consistent declines in primary production than CMIP5 [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF][START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] . This raises the question of the sensitivity of our projections under such differences in climate forcing. To address this, we extended our analysis and forced our projections with a wide variety of CMIP5 and CMIP6 individual model outputs, under equivalent high carbon emission scenarios, RCP8.5 and SSP5-8.5 (see Supplementary Figs. 345678910). The most significant finding of this intercomparison was that 19 out of 20 different climate forcings provided results consistent with the average CMIP6 projection shown here. Virtually, all future oceans projected by Earth system models over the last decade produce the same outcome: biomass losses in low and mid-latitudes and gains in subpolar ecosystems. In fact, our model projects very similar average biomass changes across latitudes when forced by CMIP5 or CMIP6 average environmental conditions (maximum median differences of 3 percentage points, see Fig. 3b). However, larger differences across CMIP6 simulations resulted in more uncertain biomass change projections. We further conducted an assessment on how environmental drivers would affect mid-trophic level fauna if they acted separately. Such analysis showed that the major source of uncertainty in CMIP6-forced biomass projections comes from chlorophyll. It also shows that biomass losses in tropical and subtropical zones would be primarily driven by a combination of warming and chlorophyll decline, while biomass gains at higher latitudes would be fundamentally promoted by warming (Fig. 3c). This contrasting response of mid-trophic level fauna is explained by the different primary production regimes that comes along with warming in low and high latitude systems. Warming is assumed to increase metabolic rates such as feeding, respiration, or growth [START_REF] Heneghan | Disentangling diverse responses to climate change among global marine ecosystem models[END_REF] . In a food-limited scenario, as indicated by our low-chlorophyll projections in tropical and subtropical waters, this would promote a decline in the consumer biomass stock [START_REF] Kooijman | Dynamic Energy Budget Theory for Metabolic Organisation[END_REF] . On the other hand, warming where primary production is not limiting would promote species expansions, growth and reproduction, as it has been already observed in temperate-polar transitions zones affected by ocean warming [START_REF] Cheung | Signature of ocean warming in global fisheries catch[END_REF][START_REF] Fossheim | Recent warming leads to a rapid borealization of fish communities in the Arctic[END_REF] . Temperature and chlorophyll persist as the major drivers of biomass change in both CMIP5 and CMIP6 projections (Extended Data Figs. 7 and8), but most Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 importantly, these changes would still take place even if only warming occurred, the most likely and unequivocal ocean change projected by Earth system models [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF][START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] .

Alternative future. We additionally projected biomass changes of mid-trophic level fauna, forced by CMIP6 environmental conditions in 2080-2100 under the low-emission scenario SSP1-2.6 (Fig. 3d). That is, if nations achieved strong mitigation measures to keep global warming below 2°C [START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] . Our results show that such measures would largely halt the redistribution of echobiomes in comparison to the high-emission scenario. The contraction of echobiomes in upwelling and subpolar regions would be strongly mitigated, while oligotrophic gyres would remain virtually unaffected (Extended Data Figs. 9-10). In terms of biomass, these strong mitigation measures would contain changes below 10% in virtually all climate zones, and near zero in temperate systems (Fig. 3d). According to our calculations, all biomass changes projected under the CMIP6 low-emission scenario, either positive or negative, would be mitigated by 4 to 11 percentage points with respect to the high-emission scenario. This implies that the impact will drop by at least two thirds at all latitudes. Similar mitigation rates were estimated under CMIP5 forcing or when using night backscatter as proxy of biomass (impact reductions from 3 to 13 percentage points, see Supplementary Tables 2 and3).

Interpreting climate-forced acoustic-based projections. Our projected biomass changes are consistent

with virtually all marine ecosystem models under comparable climate change scenarios, both in magnitude and direction [START_REF] Bryndum-Buchholz | Twenty-first-century climate change impacts on marine animal biomass and ecosystem structure across ocean basins[END_REF][START_REF] Kwiatkowski | Consistent trophic amplification of marine biomass declines under climate change[END_REF][START_REF] Lotze | Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change[END_REF][START_REF] Tittensor | Next-generation ensemble projections reveal higher climate risks for marine ecosystems[END_REF] , reconciling acoustic-based and ecosystem model projections for the first time [START_REF] Proud | Biogeography of the Global Ocean's Mesopelagic Zone[END_REF] . The way our model works stands out for its simplicity. Understanding species physiology or ecological interactions is not required to determine their response to ocean warming and shifts in primary production. This is often a difficult task in the implementation of ecosystem models and a major source of uncertainty in marine fauna projections [START_REF] Tittensor | A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0[END_REF][START_REF] Heneghan | Disentangling diverse responses to climate change among global marine ecosystem models[END_REF] . Our approach simply links sonar observations of sound-scattering fauna with current environmental conditions and projects the future distribution of this fauna under new conditions. As such, it assumes that current echobiomes will follow their environmental niches. This assumption, which is ecologically coherent, has further been corroborated by studies that already detected large-scale shifts in the distribution of multiple species following new thermal niches due to ocean warming [START_REF] Cheung | Signature of ocean warming in global fisheries catch[END_REF][START_REF] Fossheim | Recent warming leads to a rapid borealization of fish communities in the Arctic[END_REF] . Yet, as any observational approach, our method presents limitations that must be taken into consideration.

Echosounders detect organisms that are efficient sound reflectors at certain frequencies according to their physical properties. For instance, the frequency used here, 38 kHz, is known to maximize the signal of gasbearing organisms such as small mesopelagic fish or siphonophores [START_REF] Proud | From siphonophores to deep scattering layers: uncertainty ranges for the estimation of global mesopelagic fish biomass[END_REF] . Weaker fluid-like organism might hence be under-represented and this should be taken into consideration when comparing our results with more specific fauna compartments from ecosystem models. In relation to this issue, the acoustic signal in subpolar echobiomes was the only one that misrepresented pelagic biomass charts elaborated with fishing Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 gears [START_REF] Gjøsaeter | A review of the world resources of mesopelagic fish[END_REF][START_REF] Moriarty | Distribution of mesozooplankton biomass in the global ocean[END_REF] . It is unlikely that rich subpolar systems host less animal biomass than low-productive ocean gyres. We therefore believe that the subpolar decay of signal observed is related to documented latitudinal changes in the backscatter-to-biomass ratio [START_REF] Chapman | Deep scattering layer spectra in the Atlantic and Pacific Oceans and adjacent seas[END_REF] , and due to the reduction of gas-bearing organisms in subpolar sytems [START_REF] Dornan | Swimbladder morphology masks Southern Ocean mesopelagic fish biomass[END_REF][START_REF] Escobar-Flores | Estimates of density of mesopelagic fish in the Southern Ocean derived from bulk acoustic data collected by ships of opportunity[END_REF][START_REF] Dornan | Large mesopelagic fish biomass in the Southern Ocean resolved by acoustic properties[END_REF] . Because of this, the expansion of rich temperate echobiomes over subpolar echobiomes might amplify our projected biomass gains at high latitudes. We attribute therefore a greater uncertainty to biomass change projections in subpolar waters, not only for the aforementioned reasons, but also because this echobiome lacked observations in the North hemisphere. Future assessments will therefore require to investigate on basin-scale acoustic properties of pelagic fauna, and increase observations at high latitudes, including the Arctic. Another limitation relates to the model's inability to predict new echobiomes in regions where projected environmental changes have no analog with historical observations. New climate will emerge in tropical waters, where existing species are expected to reduce in size and biomass due to extreme warming and low productivity conditions [START_REF] Reygondeau | Climate Change-Induced Emergence of Novel Biogeochemical Provinces[END_REF] . As these processes are not captured, this might constitute another source of underestimation in our projected declines of ocean fauna. Finally, our projections, as any other, inherit the uncertainties of climate simulations [START_REF] Bopp | Multiple stressors of ocean ecosystems in the 21st century: projections with CMIP5 models[END_REF][START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF] . Our results should then be revisited under future simulation exercises, especially those projecting primary production and low trophic level compartments, as they currently represent the largest source of uncertainty in mid-trophic level fauna projections [START_REF] Kwiatkowski | Twenty-first century ocean warming, acidification, deoxygenation, and upperocean nutrient and primary production decline from CMIP6 model projections[END_REF][START_REF] Tittensor | Next-generation ensemble projections reveal higher climate risks for marine ecosystems[END_REF] .

Concluding remarks. We use the last-generation of Earth system models, CMIP5 and CMIP6, along with sonar observations, to project the future distribution of ocean fauna under multiple scenarios of climate change. As such, this study constitutes a ground-truth reference for ecosystem model simulations and provides an alternative empirical approximation to envisage and anticipate the effects of climate change with greater confidence. [START_REF] Heneghan | Disentangling diverse responses to climate change among global marine ecosystem models[END_REF] Our results indicate that ocean warming and shifts in primary production will result in 1) substantial reorganization of biogeographical provinces, 2) significant reductions of mid-trophic level fauna in low and mid latitude systems, and 3) massive expansion of temperate species towards highlatitudes systems. Nearly 70% of the global ocean is projected to lose animal biomass. Such loss is expected to bring food insecurity worldwide [START_REF] Blanchard | Linked sustainability challenges and trade-offs among fisheries, aquaculture and agriculture[END_REF] and declines of carbon export mediated by vertical migration or the sinking of matter produced by animals [START_REF] Bianchi | Estimating global biomass and biogeochemical cycling of marine fish with and without fishing[END_REF] . Mid-trophic level fauna may hence constitute a new feedback loop component in the global climate system which may accelerate the risks and time horizons projected by the IPCC [START_REF]Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF] . We further show that strong mitigation measures to contain global warming below 2°C would reduce our projected biomass changes by at least two thirds in comparison to the high-emission scenario. Paradoxically, now that ambitious actions need to be taken to halt the worst effects of climate change, mesopelagic fish, likely the largest unexploited stock of mid-trophic level fauna worldwide [START_REF] Irigoien | Large mesopelagic fishes biomass and trophic efficiency in the open ocean[END_REF] , is in the spotlight of the fishing industry [START_REF] Grimaldo | Investigating the potential for a commercial fishery in the Northeast Atlantic utilizing mesopelagic species[END_REF][START_REF] Olsen | Can mesopelagic mixed layers be used as feed sources for salmon aquaculture? Deep Sea Research Part II: Topical Studies in Oceanography[END_REF] . Based on these results, we call for caution in managing this ). The vertical integration of these curves-a proxy for pelagic biomass-is indicated in grey at the bottom of each profile (nautical area scattering coefficient, m 2 nmi -2 ). Colours indicate the acoustic seascapes identified, which were associated to subpolar (brown), ocean-gyre (olive), subtropical (cyan), temperate (orange), upwelling (red), and low-oxygen upwelling (purple) waters.

  . Additionally, part of the night backscatter signal of migrants might be hidden above our sampling depth range. Due to this diel variation of backscatter , we ultimately used day and night acoustic estimates as complementary indicators of pelagic biomass in the remainder of this study. Both indicators together provided a more complete view of the global variation of biomass, by including the two possible backscattering scenarios in the water column.

Figure 1 :

 1 Figure 1: Acoustic seascape classification. a, The shape of day and night acoustic profiles was used to classify the ocean seascape through hierarchical clustering. b-m, Day and night profiles for each seascape class are shown in top and bottom panels, respectively. Central and side curves represent the median and interquartile range (Sv, dB re m -1). The vertical integration of these curves-a proxy for pelagic biomass-is

Figure 2 :

 2 Figure 2: Biogeography and acoustic-based biomass. Distribution of subpolar (SP), gyre (G), subtropical (ST), temperate (T), upwelling (UW), and low-oxygen (LO) echobiomes and daytime water column backscatter (m 2 nmi -2 ) as a proxy of animal biomass. Results are shown for 2000-2020 in a and b, for 2080-2100 projections in c and d, and as the present-to-future change between these two periods in e and f. Biogeographic change in e is represented with light and bold colour shades, indicating the present distribution and future expansion of echobiomes, respectively. Backscatter change in f is indicated as percentage, relative to present-time values. Future projections are forced by the average climate outputs from 13 CMIP6 Earth system models. Contour dashed and continuos lines indicate areas where the standard deviation of projected changes exceeded values of 15 and 20, respectively.

Figure 3 :

 3 Figure 3: Biomass changes by 2080-2100 under variable climate forcing. a, Day (red) and night (black) biomass projected latitudinal changes forced by CMIP6 climate projections under high-emission scenario. b, Day biomass forced by CMIP6 (red) and CMIP5 (black) climate projections under high-emission scenario. c, Day biomass projected changes forced by CMIP6 climate projections under high-emission scenario if sea surface temperature (SST, pink), subsurface dissolved oxygen (SDO, blue), or chlorophyll (CHL, green) projected environmental drivers occurred in isolation. d, Day biomass projected changes forced by CMIP6 under high (red) and low (black) emission scenarios. e, Percentage of the global ocean's area for each latitudinal range. Fixed and changing features of biomass projections are indicated at the top and at the bottom of each panel, respectively. Boxes show the median and the interquartile range while whiskers indicate the 5-95% percentile range. They represent the spread of biomass changes due to different climate forcing (7 CMIP5 or 13 CMIP6 climate forcings). Low and high emissions respectively designate SSP5-8.5 and SSP1-2.6 IPCC greenhouse emission scenarios.

  Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2

  . Regionalisations based on Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2

  Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 fundamental component of pelagic ecosystems and urge to place mid-trophic level fauna at the very centre of global climate policies and research.

  Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2

	a	Hierarchical clustering	100	b	Subpolar	d	Gyre	Class-specific mean profiles Subtropical f Temperate h Upwelling j	Low-oxygen l
			200
			300 400										Day
			500
	Latitude 20°S 0°2 0°N 40°N 60°N	Depth (m)	500 400 300 200 600 700 100	c	217 [109, 374]	e	248 [173, 398]	835 [558, 1288] g	1672 [1146, 2490] i	1994 [1416, 2834] k	1395 [725, 2536] m	Night
	40°S 60°S		600 700 308 [154, 515]		591 [385, 842]	1179 [752, 1801]	1709 [1116, 2563]	1762 [1297, 2416]	1038 [395, 2392]
		Longitude 180°90°W 0°90°E		95 85 75 65	95 85 75 65	95 85 75 65 Sv (dB re 1 m 1 ) 95 85 75 65	95 85 75 65	95 85 75 65

Acknowledgments. We acknowledge to the Australian Integrated Marine Observing System (IMOS), the French National Research Institute for Sustainable Development (IRD), the British Antarctic Survey (BAS), the Peruvian Marine Institute (IMARPE), the Pierre and Marie Curie University (UPMC), and the Spanish National Research Council (CSIC) for their generous and invaluable contribution to the public acoustic databases used in the present study. A.A. was funded by a post-doctoral IRD fellowship. This work is a contribution to and was supported by the International Joint Laboratory TAPIOCA and the Horizon 2020 UE projects PADDLE (gran agreement No. 73427) and TRIATLAS (grant agreement No. 817578).

We would like to acknowledge the time and effort devoted by reviewers to improve the quality of this work.

Data availability

All data used in the present study are publicly available. Acoustic data can be accessed through diverse internet repositories indicated in the Supplementary Table 1. Observations of sea surface temperature and dissolved oxygen are available in the 2018 World Ocean Atlas (https://www.ncei.noaa.gov/products/worldocean-atlas), and satellite chlorophyll observations can be accessed through the Ocean Climate Change Initiative data portal (http://www.esa-oceancolour-cci.org). CMIP5 and CMIP6 simulations are publicly available in the Earth System Grid Federation data portal (https://esgf-node.llnl.gov). We also provide in supplementary data the global acoustic atlases elaborated for the present study.

Code availability

Raw acoustic data from the Malaspina circumnavigation expedition were processed using the open-source software Matecho v6.7, following the standard procedures detailed in Methods. The rest of acoustic repositories were already available as processed data. Data analysis was conducted with custom-made analysis routines in Python v3. [START_REF] Angel | Vertical Distribution of the Standing Crop of Plankton and Micronekton at Three Stations in the Northeast Atlantic[END_REF] and diverse open-source Python packages indicated in Methods. All analysis routines used in the present study are available upon request.

Methods

Data collection and processing. Georeferenced single-beam acoustic data at 38 kHz, spanning from 2001 to 2020, was collated from public databases, sourced by Australian, British, French, Peruvian and Spanish research programs in the Pacific, Atlantic and Indian Oceans (Supplementary Table 1). Most of data were already available as quality-controlled, calibrated processed data, with the exception of the Malaspina circumnavigation expedition data, which were processed for this study following standard procedures. This processing included background noise correction [START_REF] De Robertis | A post-processing technique to estimate the signal-to-noise ratio and remove echosounder background noise[END_REF] , automatic removal of corrupted signal [START_REF] Ryan | Reducing bias due to noise and attenuation in open-ocean echo integration data[END_REF] , manual cleaning of persistent corrupted signal and aliased seabed, and resampling to 10 m height and 1000 m wide data cells, using Matecho v6.7 [START_REF] Perrot | Matecho: An Open-Source Tool for Processing Fisheries Acoustics Data[END_REF] . Acoustic profiles were required to present consistent depth range and resolution, and continuous water column signal, in order to be transformed into functional curves (see section below). We therefore limited all profiles to the common largest depth range available in the dataset, from 20 to 750 m depth, and they were interpolated to a consistent vertical resolution of 10 m. To ensure the continuity of the water column signal, we remove from the dataset any profile with either missing or extremely weak backscatter, the later likely introduced when cleaning low signal-to-noise regions [START_REF] De Robertis | A post-processing technique to estimate the signal-to-noise ratio and remove echosounder background noise[END_REF] . The threshold below which a profile was considered to be altered by this pre-processing step was set to -130 dB, since no biological backscatter is expected below this level [START_REF] Stanton | Review and recommendations for the modelling of acoustic scattering by fluid-like elongated zooplankton: euphausiids and copepods[END_REF] . This threshold permitted to remove profiles where the continuity of the signal was compromised by noise conditions without the risk of removing backscatter from animals. Profiles where the seabed was above 1000 m depth were removed to exclude continental shelves from the analysis, and profiles where the sun was between 0° and 18° below the horizon were removed to avoid vertical migration events during dawn and dusk, following the astronomical definition of twilight. GPS time and position of profiles along with global bathymetry data [START_REF] Gebco | A continuous terrain model of the global oceans and land[END_REF] were used to estimate the solar azimuth angle and seabed at the time of data acquisition. Profiles that passed the filtering process (equivalent to more than 350,000 km around the globe, see Supplementary Table 1) were grouped and averaged into geographical cells of 4 degrees longitude and 4 degrees latitude. Only cells with at least 10 km of acoustic data recorded were accepted to compute average profiles. This was done separately for day and night acoustic data, i.e, we generated two global datasets with day and night profiles. As the resulting averaged profiles might result from different sampling missions across the annual cycle, the monthly sampling effort for each averaged profile was recorded. This was used afterwards to correct the environmental inputs used for training the acoustic seascape predictive model (see below). In total, mean profiles were computed for 465 cells, which represented the 16% of the ocean surface (Extended data Fig. 1). Functional data analysis. Functional data analysis (FDA [START_REF] Ramsay | Functional data analysis[END_REF] ) was used to describe and analyse the variation of the acoustic signal along the depth dimension. The analysis was implemented using the FDA module available in EchoPY v1.1 [START_REF]Fisheries acoustic data processing in Python[END_REF] . This included turning the discretized data of profiles into functions and use the descriptors of these functions to explore the main modes of signal variation in the water column through functional principal components analysis (fPCA [START_REF] Ramsay | Functional data analysis[END_REF] ). In the present study, each observation consisted of a 4x4 cell with two variables: day and night Sv profiles. Each profile was decomposed into a sequence of 18 basis functions (see Extended Data Fig. 2). We used cubic splines as the most suitable basis functions to represent non-periodic functional data [START_REF] De Boor | A Practical Guide to Splines[END_REF] . With 18 basis functions from 20 to 750 m depth, each cubic spline was adjusted to the data every ~40 m. At ocean-basin scale, 38 kHz layers are usually 100 m wide or larger [START_REF] Klevjer | Large scale patterns in vertical distribution and behaviour of mesopelagic scattering layers[END_REF] , we therefore considered this resolution appropriate to outline the vertical structure of 38 kHz profiles without over-fitting. In a functional scenario, and for the application of fPCA, the counterparts of the Sv values are the function coefficients [START_REF] Ramsay | Functional data analysis[END_REF] , which were concatenated from day and night profiles. In total, 36 coefficients (18+18) for each of the 465 geographical cells were used to build a cross-covariance matrix and find a suitable number of principal components (PCs) explaining most of the dataset variance (Fig.

Version accepted by

Extended Data Fig. 2). One of the most significant advantages of fPCA with respect to the conventional one is the possibility to access the variance contained within each principal component as a function of depth [START_REF] Ramsay | Functional data analysis[END_REF] . This allowed diagnosing the main modes of profile variance in the system (Extended Data Fig. 3), simultaneously during day and night. This was essential for later interpretation of the acoustic seascape classification. Further details on the implementation of fPCA can be found in Ramsay and Silverman (2005) [START_REF] Ramsay | Functional data analysis[END_REF] .

Seascape classification. Day-night profiles averaged for each geographical cell were classified using agglomerative hierarchical clustering from the scikit-learn Python module [START_REF] Scikit-Learn | Clustering module[END_REF] . We used the function coefficients decomposed in 8 principal components (PCs) as descriptors for the classification, which accounted for more than 90% of the dataset variance. Agglomerative hierarchical clustering was chosen because it is suitable to classify globularly distributed data [START_REF] Scikit-Learn | Clustering module[END_REF] , provides valuable information on the hierarchical similarity of profiles between clusters, and showed fairly consistent results in comparison with other classification tools (see Extended Data Fig. 5) . The same first 8 PCs used in the classification were used to reconstruct the acoustic profiles, that is, keeping more than 90% of profile variance in the dataset (Extended Data Fig. 3). These profiles were then averaged for each acoustic seascape class (median and interquartile range, see Fig. 1).

Future environmental conditions. We initially analysed simulation outputs from 13 climate models extracted from the CMIP6 archive [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF] . For each model, we use the first member of the historical experiments and the corresponding member of a low and high emission future scenario from the Shared Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 Socioeconomical Pathways (SSP), namely the SSP126 and SSP585. While the SSP585 scenario assumes accelerated globalization and rapid economic and social development of developing countries coupled with the exploitation of abundant fossil fuel resources, SSP126 is an optimistic scenario designed with the aim of simulating a development that is compatible with the 2°C target, assuming climate protection measures being taken. The model selection procedure was based on the availability of temperature, chlorophyll and oxygen data for the three experiments selected for a given model. Future mean anomalies for each of these models were calculated as the averaged annual difference between the 2080-2100 and the 1995-2014 historical periods. All data were then compiled over the same 4° × 4° grid. For the sake of comparison with the previous CMIP exercise and based on the same selection procedure, we also analysed simulations outputs from 7 CMIP5 69 climate models: for each model, we used one member of the historical experiment and the corresponding member of the high emission scenario RCP85. The RCP 8.5 scenario has an identical radiative forcing level to SSP5-8.5 (i.e., 8.5 W m -2 at 2100). Future mean anomalies for these models were calculated as the averaged annual difference between the 2080-2100 and the 1986-2005 periods.

Seascape predictive model. We used a machine learning approach trained with clustering results, in order to project the acoustic seascape beyond our observations 70,71 . A Random Forest learning algorithm was used to predict the acoustic profiles at global scale. This method was used among others after testing the prediction accuracy of up to 10 models from the Python Scikit-learn classifier module [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF] by means of crossvalidation analyses (Extended Data Figs. 456). The model was trained with sea surface temperature, subsurface dissolved oxygen (150-500 m depth), and chlorophyll as inputs, while acoustic seascape classes were used as the output. We used temperature and dissolved oxygen data from the 2018 World Ocean Atlas [START_REF] Locarnini | World Ocean Atlas[END_REF][START_REF] García | Dissolved Oxygen, Apparent Oxygen Utilization, and Oxygen Saturation[END_REF] and chlorophyll from the Ocean Colour Climate Change Initiative [START_REF]Ocean Colour Climate Change Initiative dataset, Version 5.0[END_REF] . The data period chosen was between years 2000 and 2020 in order to encompass the time coverage of the acoustic dataset. A wider set of predictors were tested initially, such as temperature and dissolved oxygen at several depth intervals, the depth of mixed layer, or primary production. The optimal selection of predictors was decided using a backward elimination approach, starting with all candidate variables and recursively removing the one with less involvement in the classification result. In order of importance, subsurface oxygen, chlorophyll, and sea surface temperature were the variables that-combined-provided the best prediction accuracy during cross-validation tests. These tests were performed using a 75% of the acoustic classification as the training dataset and the remaining 25% for validating the results. The operation was repeated 100 times with random subsets for training and validation to obtain average prediction accuracy metrics (F1 score). Each time, a random shuffle operation was selected to improve randomization (see Extended Data Fig. 4). To overcome the problem of the uneven season coverage in the acoustic dataset, the model was trained with monthly-weighted environmental predictors for the period 2000-2020, based on the monthly sampling Version accepted by Nature Climate Change on 19 August 2022 Published on 29 September 2022 in https://doi.org/10.1038/s41558-022-01479-2 effort of each acoustic profile in the dataset (see Extended Data Fig. 1). For instance, if the mean acoustic profile of a given geographical cell was computed with a sample of profiles biased to the summer season, temperature, dissolved oxygen, and chlorophyll predictors for training the model would be biased to the same direction. For that particular case, only mean summer environmental conditions for the period 2000-2020 would be used as predictors.

Seascape and biomass projections. Once the Random Forest classifier was trained, acoustic seascape classes were projected at global scale using mean environmental conditions of temperature, subsurface dissolved oxygen, and chlorophyll as inputs, for the periods 2000-2020 and 2080-2100. Contemporary environmental variables were obtained from the same sources as for training the model. For future variables, we computed historical-to-future mean anomalies from CMIP6 or CMIP5 Earth System model ensemble projections and we added these anomalies to contemporary environmental observations (see "Future environmental conditions" section). Our projections include day and night acoustic seascapes in the epipelagic and mesopelagic domains (Supplementary Fig. 1), global maps with the depth distribution of epipelagic and mesopelagic sound scattering (Supplementary Fig. 2), and vertically integrated backscatter maps, as an indicator of water-column biomass across multiple CMIP5 and CMIP6 Earth system models, and under both high-emission and strong-mitigation scenarios (Fig. 2 and Supplementary Fig. 345678910). The Random Forest projection provided the probability of each acoustic seascape class to be allocated in a given location. The most voted ones were used to build the global acoustic seascape and the probabilities were used to compute weighed mean profiles for these locations, using the profiles of each acoustic seascape class (see Fig. 1). The backscatter of these projected profiles was vertically integrated and used as a metric of animal biomass in the water column. To provide measures of uncertainty due to multiple climate forcing, this biomass proxy was computed with backscatter projections forced by all the environmental outputs within the CMIP6 (n=13) and CMIP5 (n=7) Earth system model ensembles (Supplementary Figs. 3-10). Median and percentile ranges (25 th -75 th and 5 th -95 th ) were provided from these computations (see Fig. 3 and Supplementary Tables 2 and3).