**The Interactive Role of Hydrocarbon Seeps, Hydrothermal Vents and Intermediate Antarctic/Mediterranean Water Masses on the Distribution of Some Vulnerable Deep-Sea Habitats in Mid Latitude NE Atlantic Ocean**

**Luis Somoza 1,\*, José L. Rueda 2, Olga Sánchez-Guillamón 2, Teresa Medialdea 1, Blanca Rincón-Tomás 3, Francisco J. González 1, Desirée Palomino 2, Pedro Madureira 4, Enrique López-Pamo 1, Luis M. Fernández-Salas 5, Esther Santofimia 1, Ricardo León 1, Egidio Marino 1, María del Carmen Fernández-Puga <sup>6</sup> and Juan T. Vázquez <sup>2</sup>**


**Abstract:** In this work, we integrate five case studies harboring vulnerable deep-sea benthic habitats in different geological settings from mid latitude NE Atlantic Ocean (24–42◦ N). Data and images of specific deep-sea habitats were acquired with Remoted Operated Vehicle (ROV) sensors (temperature, salinity, potential density, O2, CO2, and CH4). Besides documenting some key vulnerable deep-sea habitats, this study shows that the distribution of some deep-sea coral aggregations (including scleractinians, gorgonians, and antipatharians), deep-sea sponge aggregations and other deep-sea habitats are influenced by water masses' properties. Our data support that the distribution of scleractinian reefs and aggregations of other deep-sea corals, from subtropical to north Atlantic could be dependent of the latitudinal extents of the Antarctic Intermediate Waters (AAIW) and the Mediterranean Outflow Waters (MOW). Otherwise, the distribution of some vulnerable deep-sea habitats is influenced, at the local scale, by active hydrocarbon seeps (Gulf of Cádiz) and hydrothermal vents (El Hierro, Canary Island). The co-occurrence of deep-sea corals and chemosynthesis-based communities has been identified in methane seeps of the Gulf of Cádiz. Extensive beds of living deep-sea mussels (*Bathymodiolus mauritanicus*) and other chemosymbiotic bivalves occur closely to deep-sea coral aggregations (e.g., gorgonians, black corals) that colonize methane-derived authigenic carbonates.

**Keywords:** seafloor mapping; vulnerable deep-sea habitats; deep-sea corals; chemosynthesis-based communities; vulnerable marine ecosystem; Atlantic Ocean

## **1. Introduction**

Maintaining the sustainable functioning of the global biosphere requires protection of deep-sea ecosystems, particularly because they face major changes related to human and climate-induced impacts [1]. For an effective protection, an improvement in knowledge and research is essential, because the deep sea is truly the Earth s last frontier. Although

Sánchez-Guillamón, O.; Medialdea, T.; Rincón-Tomás, B.; González, F.J.; Palomino, D.; Madureira, P.; López-Pamo, E.; Fernández-Salas, L.M.; et al. The Interactive Role of Hydrocarbon Seeps, Hydrothermal Vents and Intermediate Antarctic/Mediterranean Water Masses on the Distribution of Some Vulnerable Deep-Sea Habitats in Mid Latitude NE Atlantic Ocean. *Oceans* **2021**, *2*, 351–385. https://doi.org/ 10.3390/oceans2020021

**Citation:** Somoza, L.; Rueda, J.L.;

Academic Editor: Antonio Bode

Received: 13 December 2020 Accepted: 16 April 2021 Published: 26 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

oceans cover more than 70% of the Earth's surface, less than 0.0001% of the deep sea has been explored to date at a high resolution [2]. Seafloor bathymetry and habitat mapping are basic tools used to improve management in some areas of conservation importance and scientific interest [1]. High-resolution seabed mapping is essential for improving the current knowledge and understanding of deep-seafloor morphology, key geological processes (including sedimentary depositional and erosional processes), habitat distribution, and typology as well as availability of mineral deposits and energy resources, among other aspects. Seabed mapping is also valuable for more applied purposes such as the deployment of submarine cables and pipelines, detection of marine geohazards, development of early warning systems for volcanic eruptions, tsunamis, and earthquakes, prevention of pollution, management of fisheries, analysis of potential environmental impacts associated with deep-sea mining and for safer shipping. In summary, high-resolution seabed mapping studies provide data and information for researchers, contractors, industry, mining companies, regulators and Non-Governmental Organizations to manage and mitigate impacts on marine environments.

The main aim of the present study is to integrate bathymetry and seafloor geomorphology, water masses properties, and the distribution of vulnerable deep-sea habitats (and their associated biota) along the southernmost NE Atlantic Ocean. These vulnerable deep-sea habitats include those conformed by deep-sea corals (in this study including deep-sea scleractinians, stony octocorals, bamboo corals, soft corals, gorgonians. and antipatharians), sea-pens, sponges and chemosynthesis-based communities. This integration has been applied to five case studies located between 24◦ and 42◦ N at 90–2500 m below sea level (mbsl) including: (1) Deep-sea scleractinian reefs (dominated by *Desmophyllum pertusum*) still thriving on the Galicia Bank; (2) Chemosynthesis-based communities of hydrocarbon seeps in the Gulf of Cádiz colonized by sponges, gorgonians, scleractinians, and black corals; (3) Deep-sea stony octocoral aggregations occurring on volcanic submarine ridges located along the Passage of Lanzarote; (4) New habitats formed after the Tagoro submarine eruption in 2011–2012; and (5) Aggregations of scleractinians, antipatharians, gorgonians, bamboo corals, and sponges on ferromanganese crusts of seamounts of the subtropical Canary Island Seamount Province (CISP). These five case study areas have been mapped in detail with MBES and ground truthed by means of ROV surveys and/or coring/dredging samples. The water mass properties for each specific habitat were obtained with in-situ ROV-mounted CTD (Conductivity–Temperature–Depth sound) and Niskin bottles. The comparison between vulnerable deep-sea habitats of these case studies allows us to analyze: (i) the potential relationship between water mass properties (salinity, oxygen, potential density, oxygen and methane concentration) for each specific habitat as a function of water depth and latitude, (ii) the influence of methane seeps on the distribution of some of these deep-sea habitats; and (iii) the influence of abrupt submarine eruptions and low-temperature hydrothermal vents in the Macaronesia volcanic archipelagos.

## **2. Geological and Oceanographic Settings**

The five case-study areas include important types of world-wide recognized geomorphological seabed features such as banks, mud volcanoes (MVs), coral mounds and reefs, oceanic passages, contourite drifts, volcanic ridges, salt domes, seamounts, and submarine volcanoes [3,4].

## *2.1. The Galicia Bank*

The first case study concerns banks and reefs of colonial scleractinians (mainly *D. pertusum*) located at 700–1800 mbsl on the Galicia Bank, NE Atlantic (Figures 1 and 2A). The Galicia Bank, with an area of 2117 km<sup>2</sup> (Figure 3A), constitutes the main structural high located on the NW Iberian continental margin, 200 km west off the Galician coast. The summit of the bank yields a contourite depositional system associated with the Mediterranean Outflow Water (MOW) and consists of erosive furrows and depositional bottom-current features such as drift and sediment wave fields [5,6]. Large pavement-like plates composed of phosphorite slabs and ferromanganese crusts dominate the top of the Galicia Bank [7]. Several water masses have been recognized to influence key processes of the Galicia Bank. Water masses above 500 mbsl are dominated by the Eastern North Atlantic Central Water current (ENACW). Below the minimum salinity located near 500 mbsl, the ENACW gradually mixes with the MOW reaching a peak in salinity values (36.1–36.5 psu) at approximately 1000 mbsl. The MOW, in this area, is characterized by an increase of the salinity values at 800 and 1200 mbsl and acts as a contour current above and around the Galicia Bank reaching velocities of 5–10 cm s−<sup>1</sup> [8].

**Figure 1.** Location of case studies areas along the North East and Central Atlantic Ocean and schematic oceanographic circulation affecting the cases studies. CISP: Canary Islands Seamount Province; PoL: Passage of Lanzarote. Main currents from [8–12] are labelled as: North Atlantic Current (NAC), Azores Current (AC). Regional bathymetric map extracted from the EMODnet Project (http://www.emodnet.eu/bathymetry, accessed on 1 December 2020).

#### *2.2. Gulf of Cádiz and Moroccan Atlantic Margin*

The second case study refers to chemosynthesis-based habitats and other vulnerable deep-sea habitats (aggregations of sponges, gorgonians, antipatharians, among others) in mud volcanoes (MVs) located in the Moroccan Atlantic margin of the Gulf of Cádiz (GoC) between 34◦ and 35◦ N (Figures 1 and 2B). In the GoC, extensive geophysical research and seabed sampling have resulted in the discovery of 84 MVs and MVs-mud diapir complexes since 2000 [13–21]. Mud volcanism in the GoC is primarily associated with high concentrations of thermogenically-derived methane [22]. Different sources have been invoked for mud volcanism in the GoC: Pliocene and Late Miocene shales, the so-called Allochthonous

Unit of the GoC and even deeper sources of Mesozoic age [23,24]. Carbonates and ferromanganese rocks such as methane-derived authigenic carbonates (MDACs) [25,26] formed by microbial-mediated anaerobic oxidation of methane (AOM) [27] and Fe-rich nodules formed by later exhumation, oxidation and diagenesis of buried AOM products [28]. Three types of intermediate water masses are currently considered in the case study of the GoC (Figure 1): the Eastern North Atlantic Central Water current (ENACW), the Mediterranean Outflow Water (MOW), the North Atlantic Deep Water current (NADW), and the Subarctic Intermediate Water current (SAIW).

**Figure 2.** E–W regional bathymetric profiles along the North East and Central Atlantic Ocean showing the morpho-structural context and water depth distribution of the studied areas. (**A**) Galician margin; (**B**) Gulf of Cádiz; (**C**) Passage of Lanzarote (PoL), (**D**) Tagoro volcano, El Hierro Island, and (**E**) Canary Island Seamount Province (CISP). Profiles extracted from bathymetry of the Global Multi-Resolution Topography (GMRT) synthesis data set [29]. See Figure 1 for location.

**Figure 3.** Seabed distribution of scleractinian reefs and mounds (M) in the Galicia Bank: (**A**) 3D shaded relief model using multibeam bathymetry data; Sunlight from the NW; (**B**) Morphological map including morpho-sedimentary features and scleractinian reefs identified with multibeam and backscatter data; (**C**) Ultra high-resolution sub-bottom profile crossing scleractinian reefs were (**D**,**E**) living *Desmophyllum pertusum* and *Madrepora oculata* appear.

## *2.3. The Passage of Lanzarote-Canary Islands-NW African Margin*

The third case study focuses on the Passage of Lanzarote (PoL), between the Canary Islands (Fuerteventura-Lanzarote Islands) and the West African Continental Margin (Figure 1). The bathymetric profile transverse to the PoL is asymmetrical with the passage narrowing towards the south, showing a steeper western flank (Figure 2C). The PoL has maximum water depths at its center, ranging between 1240 and 1460 mbsl (Figures 1 and 2C). Two groups of seafloor morphologies stand out in the submarine relief of this passage [30,31]: (a) seamounts and hills related to volcanic ridges and salt domes and (b) bottom current-related features. Three main types of water masses have been identified circulating through the PoL (Figure 1) [10,32,33]: The upper thermocline North Atlantic Central Water current (NACW), which spans from the surface to the neutral density (Υn) value of 27.3 kg/m<sup>3</sup> (approximately 600 mbsl); and the Antarctic Intermediate Water current (AAIW) reflected by a clear S minimum and found below the NACW (roughly 700–1600 mbsl).The AAIW is well identified along the African continental margin with a predominantly northward flow which often exceed 0.02 m s−1, appearing particularly intense (up to 0.06 m s−<sup>1</sup> near the bottom) during the summer seasons [33]; the third, the Mediterranean Modified Water current (MW), a modified branch of the MOW identified in individual profiles by high-salinity peaks flowing southwards between 900 mbsl to the bottom with a mean southward flow of −0.05 Sv [33]. In winter and spring, the formation of active meddies (Mediterranean Outflow Water eddies) have been reported at the north of the PoL between 12◦ W and 15◦ W [33]. The occurrence of submarine mounds along this deep-water passage promotes the intensification of the turbulence at the interface between the AAIW and MW bottom currents [31,32].

## *2.4. El Hierro Island, Canary Islands*

The fourth case study highlights the natural successional changes in newly formed volcanic substrates and habitats related to recent submarine eruptions that occurred in 2011–2012, south of El Hierro Island, Canary Islands (Figures 1 and 2D). Dead fish and patches of pale colored water at the sea surface indicated the onset of a submarine eruption on 10 October 2011 at 5 km distance from the town of La Restinga. The El Hierro eruption continued during late February, when seismicity decreased drastically until 6 March 2012. As a result of this shallow submarine eruption, a new submarine volcano grew from 375 to 89 mbsl [34]. Since 2016, the new volcano appears on the official hydrographic charts as "Tagoro" (meaning "stone circle for meeting place" in the aboriginal language of Canary Islands). The eruption consisted of two main phases of edifice construction intercalated with collapse events. After the cessation of the volcanic eruption, emissions of large quantities of iron and carbon dioxide from submarine hydrothermal vents were reported [35,36]. Several oceanographic expeditions were carried out during and after the eruption of the Tagoro volcano (Supplementary Table S1). These expeditions have allowed the identification of the development of new habitats after the eruption, associated with the emission of large quantities of iron and carbon dioxide from hydrothermal submarine vents.

#### *2.5. The Canary Island Seamounts (23*◦*48 N–26*◦*30 N)*

The fifth case study comprises the southern seamounts of the Canary Island Seamount Province (CISP), located southwest of the Canary Islands and 100 km west off the African coast (Figures 1 and 2E). The CISP is composed of more than 100 seamounts and submarine hills rising from ~5000 to 200 mbsl [37]. The CISP extends from the Lars-Essauira seamount, approximately 400 km from the north of Lanzarote (32.7◦ N, 13.2◦ W) to the Tropic seamount, 500 km at the southwest of the archipelago (23.8◦ N, 20.7◦ W). This case study focuses on the southern CISP area which is characterized by the presence of several large seamounts, including The Paps, Echo, Drago and Tropic, rising from ~5000 m to 200 mbsl (Figures 1 and 2E). These are the oldest seamounts of the CISP, dated between 91 and 119 Ma [38]. The seabed of these seamounts is mainly composed of ferromanganese crusts and nodules growing on basalt-sedimentary rock substrates containing high average

amounts of Fe (23.5 wt%), Mn (16.1 wt%), and trace elements like Co (4700 μg/g), Ni (2800 μg/g), V (2400 μg/g), and Pb (1600 μg/g) [39–41]. The oldest age of initiation of the ferromanganese crusts growth was estimated at 90 Ma [40]. The CISP seamounts are located beneath the northern end of the Oxygen Minimum Zone (OMZ), which extends from 100 to 700 mbsl, from near the equator to 25◦ N [42]. The core of the OMZ is located at 400–500 mbsl, reaching very low oxygen values (<50 μmol/kg). The dissolved oxygen is also depleted by the high biological productivity in proximity to the Canary Islands due to upwelling and nutrient-rich currents [42]. At intermediate depths, the low oxygen waters are ventilated by the ENACW and by the South Atlantic Central Water current (SACW) above ~700 mbsl and by the AAIW at ~1000 mbsl. The North Atlantic Deep-Water current (NADW) flows at depths below 1500 mbsl and the Antarctic Bottom Water current (AABW) flows below ~4000 mbsl [33] (Figure 1). These deep bottom currents accelerate around the seamounts, increasing the supply of nutrients related to high inputs of Fe and P sourced from the Sahara dust [39].

## **3. Materials and Methods**

Mapping the seafloor and its specific habitats (and Vulnerable Marine Ecosystems, VMEs) at regional scale requires interdisciplinary research teams using a wide suite of surveying techniques. Data of this work were collected mainly during the SUBVENT-2 cruise aboard RV *Sarmiento de Gamboa* [43]. Other data included in this work were collected during a set of extensive scientific expeditions listed in the Supplementary Table S1. A summary of the cruise details for each of the five case studies including data sets collected, MBES systems and ground truthing data sets are listed in Supplementary Table S1.

The multibeam bathymetry echosounder (MBES) Atlas DS 1 × 1 on board the RV *Sarmiento de Gamboa* has been used for seabed mapping of the distinct habitats and to plan ROV tracks (see Supplementary Table S1 for further details on MBES systems). MBES data were processed using CARIS HIPS&SIPS™ version 9 software (Teledyne CARIS, Frederiction, NB, Canada). The MBES data were used to generate Digital Terrain Models (DTM) at spatial resolutions of 30–50 m (depending on water depth). The multi-resolution DTM was used to produce regional sun-shaded image renders, perspective views and to extract margin-wide bathymetric profiles using Fledermaus™ version 7 software (QPS, Porstmouth, NH, USA) to interpret the submarine landscapes. It was also used to generate derivative products such as slope angle maps by means of ArcGIS™ desktop v. 12.4.

The ATLAS P-35 used aboard the RV *Sarmiento de Gamboa* is a parasound echosounder with two frequencies: a Primary High Frequency (PHF) at 20 kHz and a Secondary Low Frequency (SLF) at 4.5 kHz. PHF was used to record anomalies within the water column, whereas SLF was used as a sub-bottom profiler to record sediment/rock features (Supplementary Table S1).

Submersible observations and sampling of the seafloor were carried out with the ROV *Luso* (from Estructura de Missão para a Extensão da Plataforma Continental, Portugal), equipped with a high definition video camera Argus HD-SDI Camera (ARGUS Remote Systems AS, Laksevag, Norway), two robotic manipulators to recover biological and geological samples, a CTD (conductivity, temperature, and depth measurements) with fluorescence, turbidity and CO2-, CH4- and O2- sensors, and four Niskin bottles for water sampling (Supplementary Table S1). Water samples (two replicates of 20 mL) were taken for methane determination using a gas chromatograph in cold seep areas. Rock-mineralization samples and associated biota were taken by means of conventional rectangular benthic dredges (0.8 m in width by 0.6 m in height) towed on the seafloor during 10 to 20 min at 1 knot speed. ROV tracks were planned after analysis and interpretation of previous MBES mapping (DTM of the bathymetry) and backscatter models. Gravity cores with a maximum length of 300 cm were taken on the summits and/or slopes of the mud volcanoes. The cores were photographed and a visual lithostratigraphic description was made in each case.

Regional oceanographic data were extracted from Conductivity-Temperature-Depth (CTD) of the World Ocean Database 2018 [44] of the National Centers for Environmental Information (NOAA) (http://www.nodc.noaa.gov/OC5/indprod.html, accessed on 1 December 2020). We use Ocean Data View 4.0 software (Alfred-Wegener-Institut Helmholtz-Zentrum für Polar- und Meeresforschung, Bremerhaven, Germany) [45] for displaying regional oceanographic variables.

#### **4. Results**

#### *4.1. Description of Vulnerable Deep-Sea Habitats of Each Case Study*

In this section, the morphology and the type of habitats and communities in the five case studies from the north Atlantic (42◦ N) to subtropical Atlantic (23◦ N) have been described. Full details of substrate types, water depth ranges, geomorphology, habitat types, and the associated biota for each specific site are listed in the Supplementary Table S2.

## 4.1.1. Case Study 1: Scleractinian Reefs in the Galicia Bank (42◦15 N–42◦55 N)

A large number of live scleractinian colonies of *Desmophyllum pertusum* (previously named as *Lophelia pertusa*) and *Madrepora oculata* have been identified conforming reefs on the summit and flanks of the Galicia Bank at water depths between 620 and 1175 mbsl (Figure 3). The distribution of these scleractinian reefs has been mapped by means of MBES backscatter mosaic images, ultra-high resolution and high-resolution multichannel seismic reflection data and seabed sampling (Figure 3). A complete sequence of stepped elongated reefs (M1 to M5 in Figure 3), known as Breogham mounds and characterized by high backscatter strengths with a seabed expression, were identified on multibeam bathymetry along the western flank of the Galicia Bank (Figure 3A,B).

These semi-buried reefs display heights up to 70 m and widths of 450 m and are lined up in along-slope trending ridges from 826 to 1126 mbsl. Furthermore, on the top and at the eastern flank of the bank, a series of mini-reefs, 2–4 high and 80–100 m wide appear on a flat erosion surface at 780–750 mbsl. Mini reefs with living scleractinian colonies appear as acoustically transparent domes in very high-resolution sub-bottom profiles (Figure 3C). We interpret their low acoustic response (~1520 ms<sup>−</sup>1) close to values for the water column, as due to the high contents of seawater within their skeletons, as well as the open space created by the skeletal growth patterns. Samples of these still thriving scleractinian reefs yielded living *D. pertusum* and *M. oculata* (Figure 3D,E). In contrast, older and semi-buried scleractinian reefs were partially cemented forming elongated mounds along the flanks of the bank, being characterized by higher backscatter responses [46]. In some cases, scleractinians, mollusc monoplacophora (*Laevipilina rolani*), and bryozoans [47] occurred on encrusted phosphorites and ferromanganese nodules substrates at 750–1400 mbsl in the southeast Galicia Bank.

4.1.2. Case Study 2: Chemosynthesis-Based Communities and Other Vulnerable Deep-Sea Habitats on Hydrocarbon Seeps of the Gulf of Cádiz (34◦ N–35◦ N)

Six MVs located between 350 and 3000 mbsl (Mercator, Algacel, Yuma, Madrid, Las Negras, and Bonjardim), were explored with the ROV "Luso" during the oceanographic expedition SUBVENT-2 [43] (Location shown in Figure 4A). Chemosynthesis-based communities and high methane concentrations were detected in these explored MVs. Reduced grey mud breccias were recovered by gravity cores in all cases, when testing the nature of extrusion of these cones.

**Figure 4.** (**A**) Location of mud volcanoes explored in the present study in the Gulf of Cádiz along 35.5◦ N and the Pompeia Coral Mound Province. (**B**) Dead scleractinian skeletons at the top of the mounds. (**C**,**D**,**F**) Stony octocorals (*Corallium tricolor*) overlying live and dead colonies of *Desmophyllum pertusum* on mounds at the Pompeia Coral Province; (**E**) Ultrahigh-resolution sub-bottom profile showing giant coral mounds (transparent bodies) up to 42 m in height beneath the dead scleractinian reefs at the Pompeia Coral Mound Province.

The type of habitats identified in these MVs can be resumed into seven types (See details in Supplementary Table S2):


(7) Graveyards of scleractinians colonized by stony octocorals.

Chemosynthesis-based communities were mainly detected on the summits and sometimes on the flanks (e.g., Algacel MV), of the surveyed MVs (Figure 5) and contained different types of taxa and rates of methane flux. Regarding this, in MVs with high methane fluxes such as Algacel MV (methane concentrations up to 97.60 nM) (Figure 5A–C), harbored living deep-sea mussel *Bathymodiolus mauritanicus* beds, containing both small- and large-size specimens were found. These beds formed large linear and circular mussel clumps (up to 10 m in diameter) surrounding emissions of intermittent gas bubbling with methane concentrations of 97.60 nM (Figure 6A–E).

**Figure 5.** Vulnerable deep-sea habitats in mud volcanoes (MVs) of the Moroccan margin (Gulf of Cádiz). (**A**,**B**) 3D models of the Algacel and Mercator MVs as well as the Northern Pompeia Coral Ridge at the Pompeia Coral Mound Province, where ROV dives were carried out. Distribution of deep-sea habitats and organisms across (**C**) Algacel and (**D**) Mercator MVs. Methane concentrations are displayed in red letters. Location of MVs in Figure 4A. See further explanation in the text.

**Figure 6.** Underwater images obtained with the ROV "Luso"during the SUBVENT-2 expedition in mud volcanoes (MVs) of the Moroccan margin. (**A**): *Beggiatoa*-like bacterial mats and two cerianthids at the summit of the Mercator MV; (**B**): Accumulation of shells of the chemosymbiotic bivalve *Isorropodon megadesmus* on muddy substrates at the summit of the Bonjardim MV; (**C**): Shells of the chemosymbiotic bivalve *Acharax gadirae* at the summit of the Yuma MV; (**D**): Beds of the living chemosymbiotic *Bathymodiolus mauritanicus* in the Algacel MV; (**E**): Extensive graveyards of the chemosymbiotic bivalve *Lucinoma asapheus* at the Northern Pompeia Coral Ridge (Pompeia Coral Mound Province); (**F**): Sea-pen (*Pennatula* sp.) on muddy sediments at the Las Negras MV; (**G**): An iron bar of anthropogenic origin colonized by small gorgonians, solitary corals (*Desmophyllum dianthus*) and crinoids at the Yuma MV; (**H**): Living scleractinians(*D. pertusum*) and stony octocorals (*Corallium*) colonizing dead scleractinian skeletons in the Pompeia Coral Mound Province.

> In MVs with moderate methane fluxes (30–50 nM), including the Mercator, Las Negras, and Madrid MVs, living chemosymbiotic bivalves such as *A. gadirae* and *Solemy aelarraichen-*

*sis* occurred with high densities of shells of *Thyasira vulcolutre*, *S. elarraichensis, Isorropodon megadesmus* and *B. mauritanicus*. Dense populations of frenulate polychaetes (*Siboglinum* sp.) also occured in some of these MVs with moderate methane fluxes (Figure 5B–D). Sulfur-oxidizing bacterial (SOB) mats forming macroscopically visible cohesive white patches of 10–30 cm in diameter were also detected (Figure 6A).

In MVs with low methane fluxes such as the Bonjardim MV (20.18–28.12 nM), the past occurrence of chemosynthesis-based communities is revealed by abundant shells of the chemosymbiotic bivalve *I. megadesmus* (Figure 6B).

Crater-like pockmarks identified on the MVs seafloor, 1–3 m in width, provide evidence of former blow-outs that are nowadays colonized with non-chemosynthetic organisms such as sea-pens (*Funiculina quadrangularis*), cerianthids, and annelids (*Hyalinoecia tubicola*).

Aggregations of gorgonians (*Swiftia*), antipatharians (*Bathypathes, Leiopathes, Stichopathes*), bamboo corals (*Chelidonisis*), and scleractinians (*D. pertusum*) colonize MDACs blocks (up to 1 m) formed by anaerobic oxidation of methane (AOM) on the summit and flanks of the MVs with high rates of diffuse methane release such as the Algacel MV (Figure 5C). Otherwise, large desmosponges aggregations (*Geodia* sp., *Phakellia* sp.) colonized blocks of MDACs formed on the summit and flanks of MVs with lower methane fluxes (Figure 5D). Along the flanks of some MVs, muddy sediments covered with scattered patches of dead scleractinians (coral rubble) and small MDAC gravels are colonized by hexactinellid sponges (*Pheronema, Hyalonema*) and gorgonians (*Radicipes* on soft substrates and *Swiftia* on MDAC gravels).

Other non-chemosymbiotic deep-sea species such as sea-pens (*Kophobelemnon* sp., *Pennatula*, *Anthoptilum*), cerianthids, holothurians, and decapod species of commercial interest (*Aristeomorpha foliacea, Plesiopenaeus edwardsianus*) were also detected on the soft substrates of MVs (Figure 6F).

Ridges and mounds with abundant scleractinian skeletons (mainly *D. pertusum* and *M. oculata*), but with scarce living *D. pertusum* colonies and stony octocoral colonies (e.g., *C. tricolor*) were detected westwards of the Algacel MV in the Pompeia Coral Mound Province (Figure 4). Patches of chemosymbiotic bivalve shells (mainly *L. asapheus* with some *Thyasiravulcolutre)*, and scattered bacterial mats (*Beggiatoa*-like sulfur oxidizers) were also detected (Figure 6E,H). Ultra-high-resolution profiles below the scleractinian skeletons show that they represent the top of giant buried mounds that reach up to 42 m in height (equivalent to 50 ms two-way travel time, TWT, assuming an average sonic velocity of 1750 m s−<sup>1</sup> for recent sediments) (Figure 4E). The potential causes of the regression of these giant coral mounds in the GoC will be considered in the Section 5. Methane concentrations on the seafloor of the ridge at the Pompeia Coral Mound Province range from 41.93 to 43.24 nM, indicating the occurrence of moderate seepage. The E-W orientation of these ridges at the westward side of the Algacel MV (Figure 5A) points to the influence of strong bottom currents causing strong turbulences on their lee face.

4.1.3. Case Study 3: Stony Octocorals, Gorgonians, Soft-Corals and Sponge Aggregations along the Passage of Lanzarote (28◦30 N–29◦ N)

The third case study focuses on the hard-substrate and sedimentary habitats of the Passage of Lanzarote (PoL) located between the Canary Islands (Fuerteventura-Lanzarote ridge) and the West African Continental margin (Figure 1). Three main types of morphologies can be described based on the bathymetric map along the PoL: Elongated mounds (M), volcanic ridges, volcanic cones and sub-circular sea-floor depressions (Figure 7)

**Figure 7.** Bathymetric map obtained along the SUBVENT-2 expedition on board RV *Sarmiento de Gamboa*. The main characteristics of the seafloor such as sub-circular domes, volcanic ridges, marginal channels and circular to elongated depressions are highlighted. Locations of dives (from SV2\_D12 to SV2\_D15) are also shown as circles. M1: Dome; WTP: Western Twin Pool; ETP: Eastern Twin Pool.

Mounds are cylindrical or sub-rounded highs, 2–10 km in diameter with flat summits (e.g., M1 in Figure 7). The base of the mounds along the PoL ranges between 1225 and 1530 mbsl and their summits between 828 and 1336 mbsl. The volcanic ridges are 6 km long, 150 m high, E-W trending linear ridges detached from the mounds (e.g., the *Volcanic Ridge* in Figure 7). Volcanic cones are single structures 15–70 m height and 350–700 m in diameter crossing the flat summits of the mounds (e.g., the *Volcancito* in Figure 7). Six circular or slightly elliptical depressions have also been identified in the continental slope of the West African continental margin at the PoL [48]. They are located between 750 and 1415 mbsl, and their reliefs vary between 40 and 240 m. The length and width of their axes range from 2 and 4.5 km to 2 and 3.3 km, respectively. The surveyed deep depressions have named as the Twin Pools [48] (Western and Eastern Twin Pool, WTP and ETP respectively in Figure 7). Otherwise, the basin morphology of the PoL is dominated by channels and moats related to the interaction of the strong bottom currents with the seabed (Figure 7). Therefore, a NE-SW central channel 100 km in length and 1–5 km in width is located at 1290 mbsl in the central sector and deepens toward the NE (1460 mbsl) and towards the SW (1320 mbsl). Rimmed depressions around the mounds classified as moats are 0.5–2 km widths incising to 180 m depth. Moats are better developed along the east and south sides of the mounds reaching up to 1270–1530 mbsl at their bases.

The vulnerable deep-sea habitats identified within the PoL (Figure 8) can be resumed into the next types (see details in the Supplementary Table S2):


(4) Deep-Sea Desmosponges aggregations on deep giant circular depressions (*Western and Eastern Twin Pools*) covered by thin layers of fine sediments at 1200–1300 mbsl.

**Figure 8.** Underwater images obtained with the ROV "Luso"during the SUBVENT- 2 expedition along the PoL. (**A**): The hexactinellid sponge *Pheronema carpenteri* on soft substrates of the northern flank of M1 mound; (**B**): A Venus fly-trap anemone on soft substrates of the northern flank of M1 mound; (**C**): Brissingid starfishes colonizing small volcanic rocks at the summit of M1 mound; (**D**): The starfish *Nymphaster arenatus* on volcanic rocks at the summit of M1 mound; (**E**): Large colonies of the stony octocorals *Corallium* spp. (*C. niobe* and *C. tricolor*) on basaltic rocks of the volcanic ridge located at the northeastern part of M1 mound; (**F**): Basaltic rock colonized by different glass sponges, octocorals (*Corallium*, *Swiftia*) and hydrozoans at the Volcancito Hill. See Figure 7 for location.

> Mounds along the PoL harbor a mixture of fauna from typical sedimentary and rocky habitats, including the hexactinellid sponge (*Pheronema carpenteri*), venus flytrap anemones, isolated stony octocorals (*Corallium*) and bamboo corals (*Acanella arbuscula*), decapod crustaceans (*A. foliacea, P. edwardsianus*), starfishes (*Hymenodiscus, Nymphaster*), holothurians (*Mesothuria*), and leather sea urchins (*Araeosoma*).

> The volcanic ridges and volcanic cones are colonized by a wide variety of suspensionfeeding species, such as stony octocorals (mainly *C. tricolor* and *C. niobe*), black corals (*Leiopathes glaberrima*), bamboo corals (*A. arbuscula*), and small and large gorgonians (*Swiftia*, Plexaurid gorgonians), together with crinoids and hexactinellid and lithistid sponges, among a wide variety of benthic species. Some areas covered by abundant coral rubble were colonized by isolated bamboo corals (*A. arbuscula*). The occurrence of strong bottom

currents supports the fact that tree-like branches of some octocorals are oriented along an E-W direction (Figure 8E)

Giant circular depressions, such as WTP and ETP (Figure 7), contain a mixed seabed with soft and carbonated flagstones. The sandy and muddy substrates are colonized by typical sedimentary bathyal species such as sponges (*Thenea*), cerianthids, sea-pens (*Pennatula*), echinoderms (*Ceramaster, Mesothuria, Araeosoma*), decapods (*A. foliacea, P. edwardsianus*, *Plesionika* spp.) and large scalpellid barnacles. The rocky substrates are covered by a layer of sediment that is poorly colonized by fauna, including some massive demosponges (Pachastrellid-like sponges), small encrusting sponges, scattered hydroids and solitary scleractinians (*Caryophyllia*).

4.1.4. Case Study 4: Low-Temperature Hydrothermal Habitats at Tagoro Volcano, El Hierro Island (27◦35 N)

Recent submarine eruptions in El Hierro Island from 2011 to 2012 have provided the opportunity to map the development of a new submarine volcano, named as Tagoro, and the identification of habitats and communities associated with emissions of large quantities of iron and carbon dioxide from hydrothermal submarine vents.

The El Hierro eruption started on 10 October 2011, within a pre-existing submarine valley at 250–350 mbsl (Figure 9A). Initial bathymetric data collected on 25 October 2011, [49] identified a single cone around 205 mbsl (Figure 9B). MBES data from 29 November 2011, indicated that this volcano developed four vents with summits at 220, 195, 180 and 165 mbsl [34]. On 22 December, new bathymetric data showed a drastic collapse of the main edifice that was not present in previously acquired bathymetry. On 6 March 2012, the extrusive magmatic activity ceased 137 days after the onset of the eruption, leaving the summit of the new submarine volcano at 89 mbsl. Hydrothermal activity with intensive bubbling released from active vents continued until at least June 2012.

ROV observations in 2014 showed that the base of the Tagoro volcano is dominated by an area with extensive accumulations of scoriaceous bombs at 280 to 380 mbsl formed from basanitic lava emissions (Figure 9C,E,F). When observed closely, these bombs appeared to be either whole intact lava balloons or broken fragments of larger balloons. These bombs correspond to lava balloons floating on the sea surface during the eruption (Figure 9E). Echosounder images of the water column taken during the eruption overlie on the 3D multibeam bathymetry model showed high-reflective bright spots interpreted as these floating bombs sank along the flanks of the volcano (Figure 9A).

Two main habitats were recognized after two years of the cessation of the eruption (see details in Supplementary Table S2):


ROV images of the summit of the volcano show 5 m-high chimneys or hornitos with numerous degassing conduits punctuating their walls and marked by yellow mats, linked to sulfur-related bacteria (Figure 9D). Temperatures measured during the 2014 survey (Supplementary Table S1) showed abrupt increases in temperature up to ~2.69 ◦C, related to active hydrothermal chimneys or hornitos located at the volcano summit [50]. Peaks of CO2 emissions were detected with the ROV sensors around these hornitos [50]. On the contrary, the flanks of the volcano showed recolonization of fauna such as small oysters (*Neopycnodonte cochlear*) and serpulids, shrimps (mainly *Plesionika*), and eels (*Conger conger, Gymnothorax*) living on caves and crevices generated by the cooling of the plastic bulbous lavas flowing downslope at 250–300 mbsl.

**Figure 9.** Sequence of MBES during and after El Hierro submarine eruption superimposed on previous bathymetry (in white colour): (**A**) Bathymetry from November 2011 during the eruption when the summit was located at 220 mbsl. Overlay on the bathymetry, a 2D image of parasound echosounder showing the volcanic eruption rising throughout the water column; (**B**) Bathymetry after the eruption in June 2012 (summit located at 89 mbsl). ROV images taken in April 2014 showing distinct habitats developed after the eruption: (**C**) Pillow-lavas flowing downslope with abundant shrimps (mainly *Plesionika*), serpulid worms and small oysters (*Neopycnodonte cochlear*, some of the small white dots on the rocks); (**D**) Hornitos (hydrothermal vents) developed at the summit of the volcano draped by yellow and orange iron-sulfide bacterial mats; (**E**) Lava balloons along the slope of the recent volcano showing demersal fauna such as moray eels (*Gymnothorax*); and (**F**) Accumulations of lava bombs colonized by shrimps (*Plesionika*).

4.1.5. Case Study 5: Ferromanganese Crust-Bearing Seamounts of Southern Canary Islands (23◦48 N–26◦30 N)

In this fifth case study, we present the data from four ferromanganese crust-bearing seamounts (Echo, The Paps, Drago, and Tropic) located in the west-southwest region of the CISP at water depths ranging from 300 to 4300 (Figure 10). The Echo Seamount is a

sub-circular seamount of 10 km in diameter with the flat shallow summit located at only 300 mbsl and the base at −3700 mbsl (Figure 10A). The Paps seamount has the shape of a NW-SE ridge, 40 km length in the SE direction with the summit at 1600 mbsl and the base at 4300 mbsl (Figure 10B). The Tropic seamount is a star-shaped guyot, with its base located at 4200 mbsl and its summit at about 1000 mbsl (Figure 10C). The Drago Seamount is elliptical in shape with major axis NW-SE oriented. The summit is situated at 2200, whereas the base at 3000 mbsl (Figure 10D).

**Figure 10.** Location of the studied seamounts in the southwestern part of the Canary Islands Seamount Province (CISP) including the positions of the benthic dredge samples (DR) in (**A**): Echo Seamount; (**B**): The Paps Seamount; (**C**): Tropic Seamount; (**D**): Drago Seamount. See Supplementary Table S2.

The full list with the ROV data including substrate, morphology, water depth, oceanographic values, and associated taxa for each of the seamount surveyed can be found in the Supplementary Table S2. The habitats identified within the surveyed seamounts (Figure 11) can be resumed into the next types (Supplementary Table S3):


The summits of shallow guyot-like seamounts (e.g., Echo Seamount) are colonized by habitat-forming species from hard substrates such as large sponges (*Pachastrella, Poecillastra*), scleractinians (mainly *Dendrophyllia cornigera*) and antipatharians (*Stichopathes*) (Figure 11A,B). Other dominant fauna includes molluscs (*Asperarca nodulosa, Clelandella*, different muricids), small hydroids and echinoderms (mainly cidarids and the crinoid *Thalassometra* cf. *lusitanica*). The southern flank of this shallow seamount (DR03, 1757–1949 mbsl, in Figure 10A) showed a higher abundance of sessile species that included large gorgonians (*Metallogorgia melanotrichos*), small encrusting sponges as well as antipatharians (*Bathypathes*) and different serpulids (Figure 11C,D). Solitary and colonial scleractinians (*Desmophyllum*-like corals, *Solenosmilia variabilis*), bamboo corals, and large brachiopods were also collected in that area (Figure 11E).

**Figure 11.** Images of some samples obtained during the DRAGO0511 expedition [48] in the Echo (**A**–**E**), The Paps (**F**–**J**), Tropic (**K**) and Drago (**L**) seamounts at the Canary Islands Seamount Province. (**A**): Rocks, sponges and corals; (**B**): A living colony of the scleractinian *Dendrophyllia cornigera* collected at the summit; (**C**): A gorgonian (*Metallogorgia melanotrichos*) harboring an euryalid ophiuroid collected at the flank; (**D**): A small antipatharian (*Bathypathes*) collected at the flank; (**E**): Coral rubble retrieved from the SE flank; (**F**): Holaxonian gorgonian with abundant epibionts captured at the summit; (**G**): *Metallogorgia* gorgonian collected at the summit; (**H**): Fe-Mn crusts and hexactinellid sponges retrieved from the summit; (**I**): Fe-Mn crust samples obtained at the flank; (**J**): Goniasterid starfish collected at the flank; (**K**): Coral rubble of scleractinians (mainly *Solenosmilia variabilis* and *Desmophyllum*-like cup corals) collected at the eastern; (**L**): Thick Fe-Mn crust on volcanoclastic rock substrate retrieved from the summit.

Ferromanganese crusts of The Paps Seamount (Figure 10B) are colonized by large and fragile deep-sea gorgonians (*Metallogorgia melanotrichos, Chrysogorgia, Paramuricea*) (Figure 11F–H), bamboo corals, and small soft corals (*Anthomasthus*), hexactinellid sponges (*Aphrocallistes*), together with crustaceans (*Munidopsis*) and echinoderms that are generally attached to the gorgonians (mainly the ophiuroids *Ophiocreas* and *Asteroschema,* as well as Antedonidae crinoids). The ferromanganese crust substrate of the northern flank of this seamount (DR10, 2839–3010 mbsl in Figure 10B) is colonized by hexactinellid sponges, gorgonians (*M. melanotrichos*) and echinoderms (brisingid and goniasterid starfishes, Figure 11I,J, ophiuroids on the gorgonians, e.g., *Ophiocreas* and *Asteroschema*).

The ferromanganese crusts and carbonate-phosphorite pavements of the guyot-like Tropic Seamount (Figure 10C) are colonized by scleractinians (mainly *Solenosmilia variabilis* and *Desmophyllum*-like solitary corals) (Figure 11K), large colonies of bamboo corals (*Acanella*, *Keratoisis* up to 60 cm long), hexactinellid sponges (mainly *Poliopogon amadou*), polychaetes (mainly eunicids), and stalked crinoids (*Endoxocrinus wyvillethomsoni*). The

associated fauna included hexactinellid sponges, antipatharians (*Bathypathes*), ophiuroids, polychaetes and crustaceans (*Munidopsis*, balanids, pandalid shrimps).

The deepest seamount, the Drago Seamount (Figure 10D), harbored colonies of stony octocorals (*Corallium*) and the gorgonian *M. melanotrichos*, together with hexactinellid sponges (*Aphrocallistes, Regadrella*) and small unidentified gorgonians at its summit (DR13, 2290–2426 mbsl, Figure 10D) colonizing thick ferromanganese crusts (Figure 11L).

## **5. Discussion**

## *5.1. Potential Drivers of Deep-Sea Habitat Distribution from Subtropical North Atlantic*

This section discusses the potential drivers affecting the distribution and biodiversity of vulnerable deep-sea habitats, including both chemosynthesis and non-chemosynthesisbased habitats. Some of these drivers are: (a) the seafloor water mass properties as temperature, salinity, dissolved oxygen and potential density [12]; and (b) the active geological processes affecting the seafloor such as (i) methane seeps [51–55] and (ii) submarine volcanic eruptions followed by low-T hydrothermal degasification [50]. Furthermore, an overview of the effects of methane seeps and low-T hydrothermal vents following recent submarine eruptions on the distribution of deep-sea habitats at regional scale along the northeast Atlantic Ocean from 24◦ N to 42◦ N is provided.

#### *5.2. Influence of Water Mass Properties on Vulnerable Deep-Sea Habitats*

The relationships between water depths and water mass properties (temperature, salinity, dissolved oxygen, and potential density) for specific vulnerable deep-sea habitats highlight important considerations on suitable environmental conditions for them. The distribution of undercurrents along the NE Atlantic from the NW African Margin and Canary Islands to the north Iberian Margin is mainly driven by the outflow of the MOW at the Gibraltar Strait and by the intrusions of the Antarctic-derived AAIW currents at intermediate waters (Figure 12). The ROV-mounted CTDs data are in-situ measurements of the benthic layer, where the deep-sea habitats allow us to compare the regional distribution of water masses (Figure 12) with the local variables affecting the specific habitats (Figures 13–15). A full list of oceanographic measures in-situ for each site surveyed with the ROV-mounted CTD is shown in Supplementary Table S3.

**Figure 12.** Distribution of water masses from the Subtropical Atlantic to the North East Atlantic based on regional salinity values. Mediterranean Outflow Waters are mainly flowing northwards from the Gibraltar Strait affecting the Gulf of Cádiz and west Iberian Margin. Otherwise, Antarctic Intermediate Waters are flowing northwards along the NW African Margin throughout the Passage of Lanzarote. Data from World Ocean Database 2018 [44].

**Figure 13.** T-S diagram for the in-situ water mass conditions measured with the ROV for each of the type of habitats (full data listed in Supplementary Table S3). Background grey pattern corresponds to CTD data and dashed lines (from top to bottom) to constant potential density (isopycnal) contours = 26.51, 27.30 and 27.82 dividing four regions corresponding to the main different water masses defined by [33]: NACW, AAIW, MOW and NADW defined by [33]. Legend: FAO Taxa Classification and Codes for Vulnerable Marine Ecosystems (VMEs) (http://www.fao.org/in-action/vulnerable-marineecosystems/vme-database, accessed 1 December 2020).

**Figure 14.** Seafloor water mass parameters vs. water depths for each of the identified vulnerable deep-sea habitats and their exposure to the main bottom currents: (**A**) Temperature; (**B**) Salinity; (**C**) Potential Density, and (**D**) Dissolved Oxygen. Color legend for deep-sea habitat types is the same as Figure 13. AAIW: Antarctic Intermediate Water current; MOW: Mediterranean Outflow Water; NADW: North Atlantic Deep-Water current. Data are full listed in Supplementary Table S3. Further explanation in the text.

Figure 13 shows the T-S diagram for the in-situ water mass conditions measured for each type of habitat (data can be found in Supplementary Table S3). In this T-S diagram, most of vulnerable deep-sea habitats are located within the field of the intermediate waters: AAIW and MOW bounded between 27.30 and 27.82 isopycnals.

The three types of deep-sea coral habitats are located between these two intermediate water masses. Therefore, the aggregations of stony octocorals and bamboo corals along the PoL are bathed by AAIW, which is reflected by a clear S minimum (S = 35.4). Otherwise, aggregations of gorgonians, scleractinians, and antipatharians of the GoC as well as scleractinian reefs of the GB are bathed by the MOW, characterized by warmer and saltier waters (S = 35.8–36.2 psu). On the contrary, the sponge aggregations seem to follow a distinct pattern, with those of Desmosponges related to the NACW at uppermost levels in the PoL, whereas those of Hexactellinids in the GoC related to intermediate waters MOW-AAIW. The main difference in the water masses of these two types of sponge aggregations is the T values (Figure 13).

**Figure 15.** Seafloor water mass parameters vs. latitudes for each vulnerable deep-sea habitat detected in this study: (**A**) Temperature; (**B**) Salinity; and (**C**) Potential Density. Colour legend is the same as in Figure 13.

Chemosynthesis-based habitats of the GoC are also located under the influence of the MOW, except for the deepest mud volcano (e.g., Bomjardim MV), which is located under the influence of the NADW. The values for low-T hydrothermal habitats of the summit of the Tagoro volcano with high S and T are disturbed by the influence of hydrothermal fluids not reflecting its location at the thermocline. In contrast, benthic habitats of volcanic caves at the flanks of Tagoro are bathed by the NACW reflecting the major influence of water masses (Figure 13).

#### 5.2.1. Water Mass Temperatures

Water temperature has been traditionally considered as one of the most important factors influencing the distribution of deep-sea organisms due to their different physiological tolerances [12,56]. Thermal thresholds for deep-sea coral habitats (including those conformed by scleractinians, stony octocorals, gorgonians, and black corals) are well constrained to temperature ranges of 7.3 and 11.5 ◦C within the depth range of 800–1400 m (Figure 13A). Within this threshold, all these deep-sea corals appear distributed into three levels of seafloor temperatures: (i) scleractinians conforming reefs (*Desmophyllum pertusum*, and *Madrepora oculata*) located at the maximum temperature tolerance of 11–11.5 ◦C; (ii) Gorgonians, scleractinian and antipatharians conforming aggregations at intermediate temperatures of 8.8–10 ◦C; and finally iii) Stony octocoral aggregations at the lowest temperatures of 7.3–7.7 ◦C (Figure 14A).

The temperature range of the distribution of the reef-forming scleractinians *D. pertusum* and *M. oculata* (4–12 ◦C) across the five case studies matches the temperature range of cold North Atlantic Intermediate waters and the temperature range where those reefforming species have been detected in the past [56,57].

On the other hand, some aggregations of fragile octocorals (*Acanella, Radicipes*) appear at the lowest temperatures of 7.3–7.7 ◦C (Figure 14A) in the case studies but at a higher temperature range (1.5–6.1 ◦C) than that detected previously for the North Atlantic, so the present study provides new data for case studies located southwards of those considered [58]. Moreover, the low temperature values could point out that stony octocoral gardens living along the PoL are somehow influenced by the cold Antarctic Intermediate Waters (AAIW) that flows between the easternmost Canary Islands and the African margin (Figures 12 and 13), which was totally unexpected previously in this study.

Otherwise, desmosponges aggregations dominated by *Geodia* sp. and *Phakellia* sp., occurring within the crater of Mercator MV (Figure 5D) at shallow depths of 350–370 mbsl, tolerate high sea-floor temperatures up to 12.5 ◦C (Figure 14A) as also detected in the northern sector of the GoC, where the MOW influence is higher [54]. On the contrary, hexactinellid sponge aggregations with *Pheronema* and *Hyalonema* detected along the flanks of the Algacel MV at 820–840 mbsl (Figure 5E) live within the thresholds of the aggregations of gorgonians, scleractinians, and antipatharians, colonizing MDACs around methane seeps (Figure 14A). Indeed, these sponges were found at similar depths and environmental conditions in the northern part of the GoC [59,60]. On the other hand, deep-sea chemosynthesis-based habitats do not show seafloor temperature constraints, varying from 3.2 to 12.5 ◦C, as function of the seeps location. In the case of the large mussels *Bathymodiolus mauritanicus*, the seafloor temperature of 10.1 ◦C may fluctuate as function of intensification of the hydrocarbon seeps. No information on the effect of temperature on the distribution of methane seeps fauna has been found, but the fact that the same chemosymbiotic-species may occur in cold seeps with different oceanographic conditions in methane seeps of the northern and southern GoC may indicate a lower dependence on water mass properties for these chemosymbiotic organisms than for the above mentioned cnidarians and sponges [60,61]. In the case of the hydrothermal-related habitats (V1 = hornitos-like hydrothermal chimneys and V2 = lava cavities in Figure 14A), temperatures reflect thermal anomalies associated with venting hydrothermal fluids. Therefore, maximum measured temperature of expelled fluids in Tagoro volcano was 39 ◦C, while the

ambient seawater range between 15 and 19 ◦C with local temperature anomalies of +2.7 ◦C associated with CO2 emissions [50].

#### 5.2.2. Salinity and Potential Density

Although salinity variations are very low in the deep sea, the distribution of deep-sea scleractinians and octocorals in some areas of the NE Atlantic has been related to potential density, which represents a parameter defined by salinity and temperature [62]. Indeed, reef-forming scleractinians of the selected sites are controlled by a narrow salinity range of 35.35–36.1 psu (Figure 14B). As seafloor temperatures, habitats conformed by deep-sea corals appear also distributed into three levels of salinity: (i) Scleractinian reefs (*D. pertusum* and *M. oculata*) occurred at salinity ranges between 35.8 and 36.1 psu; (ii) Gorgonian aggregations at salinity ranges between 35.6 and 36.1 psu; and (iii) Stony octocoral aggregations at the lowest salinity values between 35 and 35.2 psu. These data point out that there is an influence of the MOW waters bathing the scleractinian reefs, and of the mixing between the saltier MOW waters and the lower salinity AAIW waters for some gorgonian and stony octocoral aggregations (Figure 13). Indeed, the MOW flow seems to play an important role in the distribution and connectivity of deep-sea corals across the NE Atlantic, with a high diversity of deep-sea corals along the MOW pathway [63].

Deep-sea desmosponges aggregations at 1400 mbsl appear associated with salty water masses of 36 psu (Figure 13B) related to the export of bathyal benthos to the Atlantic through the MOW outflow as detected in other areas of the GoC [60]. On the contrary, shallow water desmosponges grounds associated with methane seeps of the Mercator MV (Figure 5C), located at 350 mbsl above the MOW influence, show the lowest levels of salinity of 35.75 psu (Figure 14B). As occurred for temperature values, large hexactinellid sponge aggregations occurred within the same salinity range as for deep-sea corals (Figure 14B). The role of the MOW exporting deep-sea sponge larvae from the Mediterranean to the Atlantic Ocean and some NE locations close to the GoC has been demonstrated recently, and this may partly explain the occurrence of some sponges in salty water masses related to the MOW [64].

The occurrence of scleractinian reefs along the NE Atlantic margins has been related to values of potential density, a parameter defined by salinity and temperature. Thus, *D. pertusum* and *M. oculata* reefs, were firstly constrained to a narrow potential density range (potential density; δ<sup>θ</sup> = 27.35–27.65 kg m<sup>−</sup>3) [62] which further, in other North Atlantic regions, proved to be larger than initially suggested as ranging δ<sup>θ</sup> = 27.10–27.84 kg m−<sup>3</sup> [63]. This range is similar for the *D. pertusum* and *M. oculata* reefs identified in the Galicia Bank, which falls within this suggested range with values of δ<sup>θ</sup> = 27.353–27.431 kg m−<sup>3</sup> (Figure 14C). In the case of the Galicia Bank, these conditions seem to be caused by the turbulent mixing between the MOW and the ENACW water masses promoted by the MOW current impinging the topography of the Galicia Bank [46]. The formation of turbulent eddies of MOW flow facilitates the downward transport of nutrients, influencing the present distribution of the scleractinian reefs and mounds mapped along its flanks and on the summit (Figure 3B).

On the contrary, the gorgonian and stony octocorals aggregations do not fit with this potential density threshold (Figure 14C). Thus, stony octocoral aggregations in the subtropical latitudes occur at potential density values δ<sup>θ</sup> = 32.725 kg m−<sup>3</sup> , whereas gorgoniandominated aggregations occur at values of δ<sup>θ</sup> = 32.050–34.900 kg m−3, much higher than those for scleractinian reefs. Thus, according our data, the potential density range proposed for scleractinian reefs [62,63,65] would be only valid for the northernmost NE Atlantic margins, and more data on the presence of these species southwards is then needed for estimating the appropriate density range. One explanation is the potential influence of intrusion of AAIW flow along the NW Africa margin towards the GoC (Figure 12).

Desmosponges and hexactillenid aggregations also show different potential densities. Therefore, desmosponge aggregations show higher potential density ranges (δ<sup>θ</sup> = 32.050–34.900 kg m<sup>−</sup>3) than the studied hexactillenid sponge aggregations (δ<sup>θ</sup> = 31.058 kg m−3) (Figure 14C) Hydrography plays an important role in shaping the distribution of sponge aggregations in the Atlantic, and several authors have noted the association of sponges with particular water masses due to their temperature and salinity characteristics or hydrodynamic conditions, such as tides or internal waves which enhance the food supply. Thus, the influence of the NACW (Figure 13) in the study area allows the presence of desmosponges assemblages with a wide Atlanto-Mediterranean distribution [60,64].

#### 5.2.3. Dissolved Oxygen Concentrations

Dissolved oxygen (DO) values are close to saturation in most deep-sea areas, ranging from 2.5 mL·L−<sup>1</sup> off SW Africa to 6 mL·L−<sup>1</sup> in the NW Atlantic [65]. The values of DO for deep-sea corals in the study area are limited to values ranging from 4.3 to 5.5 mL·L−<sup>1</sup> (Figure 14D). As occurred with temperature and salinities values, these habitats show three levels of DO for deep-sea coral habitats living at 800–1400 mbsl: (i) High oxygenation values of 5.44–5.55 mL·L−<sup>1</sup> for scleractinian reefs; (ii) Intermediate oxygen values of 4.73–5.25 mL·L−<sup>1</sup> for gorgonian aggregations; and (iii) low oxygen values of 4.32 mL·L−<sup>1</sup> for stony octocorals aggregations. Otherwise, desmoponges sponge grounds show similar values of DO to aggregations of gorgonians and stony octocorals (Figure 14D). Deep-sea corals and sponges, like other organisms, are not very sensitive to small variations in DO unless it drops a certain threshold [12].

Otherwise, hypoxic or anoxia conditions are mainly related to submarine areas with methane seeps or hydrothermal vents. Seafloor anoxic conditions may be generated by intense activity of methane seeps, hydrothermal vents or submarine eruption plumes. Fluid emissions are not continuos but may be reactivated periodically on focused submarine vents or remain diffuse during long time [66]. Very low DO values of 0.64 mL·L−<sup>1</sup> observed on deep-water mud volcanoes as Bonjardim in the GoC are closely to anoxia suggesting the intense activity of methane seeps (Figure 14D). If methane seeps are intense and/or undercurrents are weak in such deep waters, then low oxygen conditions prevail above seafloor and chemosynthesis-based habitats are only composed by colonies of Siboglinidae tubes [67]. However, in areas of methane seeps influenced by strong intermediate undercurrents (Figure 14D), co-occurrence of chemosynthesis-based and scleractinian-gorgonian aggregations has been detected as in this study and other MVs of the Gulf of Mexico [68]. This is the case for the extensive beds of deep-sea mussels Bathymodiulus mauritanicus beds living around bubbling methane seeps (Figure 5D) and MDACs patches colonyzed by gorgonians, antipatharians, and scleractinians in Algacel MV (Figure 5D). In such methane seeps, the siboglinids tubeworms play an important role in the connectivity between methane seeps and deep-sea corals by filtering sulphur reducing bacteria (SRB) allowing the predominance of anaerobic oxidation of methane (AOM) by archaeas building up carbonates patches and allows deep-sea coral larvae to grow isolated from the surrounding highly acidic methane muds [69].

#### *5.3. Distribution of Vulnerable Deep-Sea Habitats from Subtropical to North Latitudes*

The biodiversity of deep-sea habitats and their biogeographical affinities can vary over space and time due to regional and basin-scale changes in the oceanographic conditions [65]. However, active geological process affecting the ocean s seafloor such as methane seeps, hydrothermal vents or submarine eruptions can modulate the spatial distribution of deepsea habitats as detected in this study. The water properties for deep-sea habitats show trends with latitude which are related to oceanographic conditions, but also anomalous values associated with methane seeps and hydrothermal vents (Figure 15). Therefore, deep-sea coral habitats show a notable decreasing trend in the suitable temperatures from north to south (Figure 15A): (i) 11–11.5 ◦C for scleractinian reefs habitats in the Galicia Bank; (ii) 8.8–10 ◦C for gorgonian and scleractinian aggregations associated with mud volcanoes of the GoC; and (ii) 7.4–7.7 ◦C for stony octocoral gardens of the PoL. Otherwise, the sponge grounds show a similar decreasing trend but at higher temperature ranges. Thus, desmosponges habitats of the GoC are related to temperatures thresholds as high as 12.5 ◦C, whereas hexactillined sponge grounds of the Pol show lower temperatures of 7.3–7.33 ◦C (Figure 15A). This negative trend identified in the sensitive temperatures for these habitats may be caused by the cooling of the AAIW intermediate waters and it may affect the composition of these communities as identified on T-S diagram (Figure 13). The AAIW flows along the NW African coast and extends to the GoC, where it interacts with the MOW outflow, perhaps restricting the spreading of the later [70].

The values of salinity from 24◦ to 44◦ N show a different trend than the temperature ones (Figure 15B). In this case, the deep-sea habitats supporting maximum salinity values are located on the mud volcanoes of the GoC, that are strongly influenced by the MOW undercurrent (Figures 12 and 13). Therefore, gorgonian and scleractinian aggregations growing on MDACs around methane seeps show maximum salinity values up 36 psu. Both to the north and to the south, the salinity threshold for deep-sea corals diminishes. Thus, the scleractinian reefs of the Galicia Bank (42◦ N) show slightly lower salinity values of 35.9 psu and stony octocoral aggregations located on the PoL (28.30◦–29◦ N) show the lowest salinity values of 35.35 psu (Figure 15B). Considering this, the low salinity and temperature values bathing the stony octocoral aggregations of the PoL are due to the influence of the AAIW flowing along the NW African (Figures 12 and 13).

The salinity threshold for the sponge aggregations follows the same descending trend from 34 to 28◦ N at regional scale (Figure 15B). Desmosponge grounds living close to methane seeps in the GoC are living under similar high salinity values to gorgonian and scleractinian aggregations (Figure 15B). However, hexactillenid sponge grounds living at shallower depths on the summits of the PoL mounds support higher salinity values than stony octocoral habitats (Figure 15B). This is interpreted as a result of the influence of periodic southwards intrusions of salty MOW waters into the PoL at water depths shallower than the AAIW undercurrents which may influence the development of those habitats.

## *5.4. Temporal Variability of AAIW Latitudinal Extension along the Northern Atlantic Ocean Might Have Caused the Massive Mortality of CWC Reefs?*

At present live *D. pertusum* and *M. oculata* frameworks are very scarce in the GoC and generally substituted by other deep-sea corals (e.g., different species of gorgonians, antipatharians, etc.). Besides, a geographical boundary between scleractinian reefs and aggregations of other deep-sea corals is observed (Figure 15B,C). A distribution possibly related to the presence of the AAIW waters (Figure 13). Evidence for such a relation is provided by the occurrence of giant dead mounds up to 30 m in height conformed by dead scleractinians (*D. pertusum* and *M. oculata*), such as in the Pompeia Coral Province in the GoC (Figure 4), pointing to a massive mortality of scleractinian reefs off northwest Moroccan [69,71–73] and along the Mauritanian margins [74].

Based on radiocarbon data of scleractinians along the Western Mediterranean Sea during the last deglaciation times [75], we hypothesize that the beginning of intrusion of AAIW waters into the GoC during the last deglacial times may have been one of the causes of a massive mortality of scleractinians by shifting northwards the biogeographical boundary between some scleractinians and stony octocoral aggregations (Figure 15B). This point out that stony octocorals (*Corallium tricolor* and *C. niobe)*, antipatharians, and gorgonians are more suitable to AAIW conditions than the scleractinians *D. pertusum* and *M. oculata.*

#### *5.5. The role of Methane Seeps Driving Distribution of Chemosynthesis and Non-Chemosynthesis-Based Habitats*

Seabed features formed by seafloor venting of hydrocarbon-enriched fluids are generically termed cold seeps, which are associated with high geological, geochemical and, biological activity [51,52]. Submarine MVs are one of the main seabed morphological expressions of cold seeps formed by violent eruptions, followed by progressive degassing of hydrocarbon-enriched fluids and sediments onto the seafloor [66].

## 5.5.1. Drivers Controlling Distribution of Habitats in Methane Seeps: Acidic Muds vs. Carbonates

Based on the habitat types identified, ROV-mounted CTD parameters and CH4 water analyses along MVs of the GoC, two main drivers controlling their formation and distribution should be highlighted: (i) the rate of release of deep-seated methane-enriched fluids, and (ii) the formation of hard substrates such as MDACs (chimneys, slabs or pavements) by AOM processes [27].

The upward migration of methane is transformed by AOM into large amounts of sulphide on the surface [76]. This extremely acid seafloor is buffered by large populations of sulphur-oxidizing siboglinid tubeworms and sulphur-oxidizing bacterial mats, allowing formation of hotspots on the surface by consuming the AOM-sourced sulphur [69]. The occurrence of siboglinid tubeworms in the anoxic-oxic zone seems to be essential for the formation and the non-dissolution of carbonates at the seabed and, moreover for the progressive colonization by deep-sea corals and other suspension feeders on these hard carbonated substrates. Thus, some scleractinians, gorgonians, antipatharians, bamboo corals, and demosponges can colonize the upper parts of MDACs blocks, slabs and pavements, up to 1 m in diameter. In areas with extensive hydrocarbon seeps, massive MDACs are the dominant substrate for coral colonization and reef formation as detected in other MVs of the GoC [54,59].

#### 5.5.2. CWC Mounds and Methane Seeps

The formation of giant carbonate mounds up to 30 m high built up by colonial scleractinians has been identified related to methane seeps in the GoC [69,71,72,77]. In the case of the Northern Pompeia Coral Ridge (westwards Algacel MV), we identified values of methane ranging from 41.93 to 43.24 nM and patches of shells of chemosymbiotic bivalves, mainly *L. asapheus* with some *Thyasira vulcolutre*, and scattered bacterial mats (sulfur-oxidizing-like *Beggiatoa*) indicating active methane seeps throughout the carbonate mounds. It has been proposed a direct relationship between carbonate reefs and fluid seepages by fertilizing waters sourced from methane seeps [78]. Recently, a new hypothesis on the formation of carbonate mounds conformed by scleractinians has been proposed in relation to methane seeps [69]. Thus, scleractinian larvae may colonize MDACs hotspots only if siboglinid tubeworms previously reduce the concentration of sulfide in the anoxicoxic zone allowing skeletal calcification of scleractinians in MDACs surrounded by the highly acidic muds of methane seeps. Transition from active methane seeps, carbonate mounds and deep-sea coral colonization stages has also been proposed as evolutionary stages of MVs in the GoC [79].

Deep-sea coral colonizing MDACs hard substrates in areas of methane seeps have been reported throughout all margins of the world: The Gulf of México [80,81], Hikurangi Margin in New Zealand [82], Brazilian margin [83], the Darwin Mounds in the northern Rockall Trough [63], the Norwegian shelf [78], and the Porcupine Basin, west of Ireland [84]. Furthermore, the co-existence of chemosymbiotic vestimentiferan worms and bacterial mats with deep-sea corals in cold seeps has been reported in the Gulf of Mexico [81]. The co-occurrence of MDACs with annelids has also been reported from hydrocarbon seeps along the US northern and mid-Atlantic margin [85]. Both MDACs and chemosymbiotic deep-sea mussels in the U.S. Atlantic margin seeps shows average <sup>δ</sup>13C signature of −47‰, a value consistent with microbially driven anaerobic oxidation of methane-rich fluids occurring at or near the sediment-water interface [85]. In these seep areas, macrofaunal densities dominated by annelid families Dorvilleidae, Capitellidae, and Tubificidae were four times greater than those present in deep-sea mussel beds habitats. These differences between chemosynthesis-based habitats have also been observed in the GoC [59,61]. In this way, it has been suggested that the difference between habitats in such extremophile environments is driven by quality and source of organic matter [86].

## 5.5.3. Type of Habitats and Methane Concentration

Based on this study, a relationship between the type of habitat and methane concentration in the water column may occur: (i) High concentrations of methane 97.60 nM with bubbling: extensive beds of living chemosymbiotic bivalves *B. mauritanicus, Lucinoma asapheus, Acharax gadirae,* or *Thyasira vulcolutre* fueled by bubbling gas through fissures and cracks; (ii) Intermediate methane concentrations ranging between 65 and 92 nM: AOM-formed MDACs accumulations colonized by aggregations of gorgonian and antipatharians (*Bathypathes, Leiopathes, Stichopathes*), bamboo corals (*Chelidonisis, Acanella*), stony octocorals (*Corallium*) and, sometimes even scleractinians (*D. pertusum*) as well as of demosponges; (iii) Low concentration ranging 20–30 nM allowing colonisation of typical non-chemosymbiotic bathyal soft bottom species such as sea-pens (*Kophobelemnon* sp., *Pennatula*, *Anthoptilum*), cerianthids, holothurians and decapods.

The influence of strong undercurrents bathing the flanks of the MVs favours the occurrence of non-chemosymbiotic fauna such as large sponge grounds (*Geodia* sp., *Phakellia* sp.) that colonize the scattered MDAC slabs, as well as sea-pens (*Funiculina quadrangularis*), cerianthids and annelids (*Hyalinoeciatubicola*) on the soft, muddy sediments.

#### *5.6. Potential Ecological Restoration of Deep-Sea Habitats after Submarine Eruptions in the Macaronesia Region*

Volcanism and associated hydrothermal systems are relevant processes for the evolution of the ocean basins, due to their impact on the geochemistry of the oceans and their potential to form significant deposits [87]. Low-T hydrothermal vents after violent submarine volcanic eruptions generate long-term CO2 inputs to oceans due to the continuous degasification of the magmatic systems mainly placed on hot-spot volcanic islands like Hawaii or Canary Islands. This is due to the high contents in C bearing in the thick oceanic sediments below the submarine volcanoes that are expulsed by low-T hydrothermal vent systems [88].

In the Tagoro volcano, recolonization of pioneering fauna such as small oysters, serpulids and mobile species (e.g., shrimps, eels) was detected, representing first colonizers of the newly formed substrates in this volcanic environment. Newly formed habitats were also detected [89] together with the burial of extensive areas with aggregations of antipatharians and some deep-sea corals. Some authors indicated that the first colonizers at Tagoro volcano included a large proportion of common suspension feeders and predators of circalittoral and bathyal hard bottoms of the Macaronesian fauna, which could have exploited the uncolonized hard bottoms and the post eruptive fertilization of water masses [35,36]. Along the flanks of this volcano, caves show high values up 6.44 mL·L−<sup>1</sup> of DO related to active CO2 degasification two years after the volcanic eruption [50]. The trend representing the DO range for deep-sea complex habitats from 40◦ N to 26◦ N (Figure 16) shows the potential pathway from the present parameters to reach suitable DO conditions for the settlement and development of slow-growing organisms such as octocorals (gorgonians and soft corals) as well as antipatharians as detected in other areas of El Hierro Island [89].

We suggest that the comparison between habitats growing after submarine eruptions at different ages might be used to infer the rate of biological colonization and the natural ecological succession of habitats after a violent submarine eruption. Previously, a similar approach has been successfully done for understanding geological and biological evolutionary stages in methane seeps displaying different scenarios of rates of fluid venting [79] and hydrodynamics [21]. Therefore, within the Macaronesia volcanic archipelagos, several recent submarine eruptions have taken place in Capelinhos (west Faial Island, Azores) in 1958–59 or south of El Hierro Island (Canary Islands) in 2011–2012 [34]. The recent discovery of soft coral gardens dominated by *Alcyonacea* in the Azores Archipelago [90] colonizing the volcanic substrate created from the eruption of Capelinhos in 19581–959, opens new studies of succession and survivorship of habitat-forming species in the Macaronesia volcanic areas as they have already been developed in the Hawaiian Islands [91].

**Figure 16.** Seafloor dissolved oxygen values vs. latitude for each deep-sea habitat. Grey pattern shows the oxygen threshold for CWC reef and gardens, and deep-sea sponges along the northeast Atlantic Ocean. Dashed black lines show the potential pathways from conditions of cold seeps and recent submarine eruptions to be suitable for non-chemosynthesis-based habitats according the level of dissolved oxygen on the benthic layer. Data are fully listed in Supplementary Table S3.

This fact would support the hypothesis that the type of seabed substrate and the occurrence of fluid flow (CO2, methane, sulfide, iron) together with some water mass properties (e.g., current speed, salinity, and temperature) might control the evolution of the type of habitats in such volcanic environments, i.e., (i) soft octocoral gardens developed after recent volcanic eruptions (~50 years ago) with latent hydrothermal CO2 degassing and Fe fertilization but causing seabed acidification unsuitable for carbonates; (ii) Coral gardens and sponge aggregations with a wide variety of sessile species (e.g., octocorals, hexactinellid and lithistid sponges, black corals, bamboo corals and large gorgonians) supported by non-degassing volcanic rocks (up thousands of years old) as the submarine ridges of PoL allowing occurrence of carbonate-bearing benthic fauna.

#### *5.7. Tools for Future Management of Vulnerable Marine Ecosystems*

The integration of high-resolution bathymetry, measurements of in-situ water mass properties, and identification of deep-sea habitats ground truthing by ROV at regional scale represents an important tool for increasing the knowledge and protection of deep-water marine ecosystems, especially those areas with VMEs, affected by the ongoing climate changes and/or contemporary anthropogenic impacts. The future scenario for the impact in deep-sea species and habitats, such as scleractinian reefs, coral gardens or deep-sea sponge aggregations, with potential increases in CO2 and CH4 emissions or rising in the seawater temperature of the oceans at global scale, might be compared with the changes produced, at local scale, by natural emissions of methane in cold seeps or carbon dioxide in recent submarine eruptions.

Otherwise, these integrated studies are also of paramount importance for the efficient management of natural marine resources, such as the ferromanganese crusts bearing seamounts with high contents in Co, Ni, and other Rare Earth Elements (REE), which are presently targets for mineral exploration in the Area under contracts with the Inter-

national Seabed Authority (ISA). These seamounts host VMEs, as detected in this and previous studies [92] as well as mineral resources [93] and, thus, spatial management plans are needed to address potential conflicts between deep seabed mining interests and the conservation of the deep-sea biodiversity presented in this study. The same applies for the hydrocarbon-bearing fluid habitats of the studied methane seeps, which are potential target areas for hydrocarbon exploitation in the subsurface and contain singular habitats with endemic species such as the chemosymbiotic bivalves of the GoC.

The five case studies presented here might be considered Ecologically or Biologically Significant Marine Areas (EBSAs) as considered by the United Nations Convention of Biological Diversity (according to different criteria of Annex I, decision IX/20). Unfortunately, significant ecosystem components of those areas (mainly benthic ones) may be threatened by human activities (e.g., fishing activities, hydrocarbon exploitation, or seabed mining), at present or in a near future. Therefore, these selected areas comprising a wide variety of deep-sea habitats in the North Atlantic should be protected and used as pilot areas for further mapping and monitoring of habitats and VMEs, contributing to the environmental management of the North Atlantic Ocean seafloor.

## **6. Conclusions**

Five case studies located between 90 and 2500 mbsl depth, along the NE and Central Atlantic Ocean from 23◦ N to 42◦ N latitudes, are presented in this work. These were targeted with MBES mapping, ROV-mounted CTD data of the benthic layer (temperatures, salinity, potential density, dissolved oxygen, potential density, CO2, CH4 concentration) and identification and sampling of deep-sea habitats. Based on these datasets, we conclude that:


**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/oceans2020021/s1. Table S1: Cruise details, MBEs system and configurations, oceanographic data and sampling methods; Table S2: Seabed substrate, morphology, oceanographic variables, VME Indicator Taxa and seafloor micro-morphologies and habitat types; Table S3: ROV-mounted CTD mass water data.

**Author Contributions:** L.S. was the chief scientist of the SUBVENT-2 cruise and led the overall organization of the text. J.L.R. directed the collection of biological observations and the biological text of the manuscript. O.S.-G. contribute to the overall text of the manuscript and biological descriptions. P.M. lead the ROV "Luso" operations. B.R.-T. directed the microbiological analyses. D.P., R.L. and L.M.F.-S. processed bathymetry, provide multibeam maps and contributed to geomorphological description of the manuscript. T.M. and M.d.C.F.-P. contributed to the geology of the manuscript. F.J.G. and E.M. contributed to the petrological parts of the manuscript. E.L.-P. and E.S. contribute to physic-chemical data. J.T.V. was the chief scientist of the DRAGO511 cruise. All authors contributed to the writing of a draft of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** Support for this work comes from the EXPLOSEA (CTM201675947-R), SUBVENT (CGL2012- 39524-CO2), and EXARCAN (CTM2010-09496-E) projects funded by the Ministerio de Ciencia e Innovación of Spain. This study also benefits from the Atlantic Seabed Mapping International Working Group (ASMIWG) as part of the Atlantic Ocean Research Alliance Coordination and Support Action (AORA-CSA). This study is also a contribution to the EMODNET-Geology project (EASME/EMFF/2018/1.3.1.8-Lot 1/SI2.811048) and the European project H2020 GeoERA-MINDeSEA (Grant Agreement No. 731166, project GeoE.171.001).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** The authors thank captains, crew and ship parties of all scientific cruises aboard R/V "*Hespérides*", "*Sarmiento de Gamboa*" and "*Miguel Oliver*". Special thanks to the ROV "Luso"pilots team Antonio Calado, Andreia Afonso, Miguel Souto and Renato Bettencourt that have been the basis for this paper. Special thanks to the personnel of the Instituto Hidrográfico de la Marina (IHM) for their cooperation in the processing of multibeam bathymetry data.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Review* **Offshore Geological Hazards: Charting the Course of Progress and Future Directions**

**Gemma Ercilla 1,\*, David Casas 1, Belén Alonso 1, Daniele Casalbore 2, Jesús Galindo-Zaldívar 3, Soledad García-Gil 4, Eleonora Martorelli 5, Juan-Tomás Vázquez 6, María Azpiroz-Zabala 7, Damien DoCouto 8, Ferran Estrada 1, Ma Carmen Fernández-Puga 9, Lourdes González-Castillo 3, José Manuel González-Vida 10, Javier Idárraga-García 11, Carmen Juan 12, Jorge Macías 13, Asier Madarieta-Txurruka 3, José Nespereira 14, Desiree Palomino 6, Olga Sánchez-Guillamón 6, Víctor Tendero-Salmerón 3, Manuel Teixeira 15,16,17, Javier Valencia <sup>18</sup> and Mariano Yenes <sup>14</sup>**


**Abstract:** Offshore geological hazards can occur in any marine domain or environment and represent a serious threat to society, the economy, and the environment. Seismicity, slope sedimentary instabilities, submarine volcanism, fluid flow processes, and bottom currents are considered here because they are the most common hazardous processes; tsunamis are also examined because they are a secondary hazard generated mostly by earthquakes, slope instabilities, or volcanic eruptions. The hazards can co-occur and interact, inducing a cascading sequence of events, especially in certain

**Citation:** Ercilla, G.; Casas, D.; Alonso, B.; Casalbore, D.; Galindo-Zaldívar, J.; García-Gil, S.; Martorelli, E.; Vázquez, J.-T.; Azpiroz-Zabala, M.; DoCouto, D.; et al. Offshore Geological Hazards: Charting the Course of Progress and Future Directions. *Oceans* **2021**, *2*, 393–428. https://doi.org/10.3390/ oceans2020023

Academic Editor: Pere Masqué

Received: 16 December 2020 Accepted: 24 May 2021 Published: 31 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

contexts, such as tectonic indentations, volcanic islands, and canyon heads close to the coast. We analyze the key characteristics and main shortcomings of offshore geological hazards to identify their present and future directions for marine geoscience investigations of their identification and characterization. This review establishes that future research will rely on studies including a high level of multidisciplinarity. This approach, which also involves scientific and technological challenges, will require effective integration and interplay between multiscale analysis, mapping, direct deep-sea observations and testing, modelling, and linking offshore observations with onshore observations.

**Keywords:** seismic faults; slope instabilities; submarine volcanism; fluid-flow processes; bottom currents; tsunamis; canyon heads; tectonic indentation; multidisciplinary approach

## **1. Introduction**

Geological processes occurring within or at the surface of the Earth may lead to natural disasters causing loss of life, environmental damage, and major impacts on the economy and food security. Human behavior may also trigger natural disaster processes where no hazards existed before or increase the risks where they do exist. Knowledge of the geological elements likely to produce a disaster and their distribution, as well as the understanding of the mechanisms (conditioning factors and triggers), is critical for the prevention and mitigation of catastrophic events [1–4]. The time scales of those mechanisms are highly variable, from geologic to human scale, which contrasts greatly with the sudden impact of such events when they are triggered.

Like in any other Earth domain, the marine environment is associated with potentially hazardous geological processes that may represent serious threats to society, the economy, and the environment (http://www.esf.org/fileadmin/Public\_documents/Publications/ Natural\_Hazards.pdf, accessed on 28 May 2020) (Figure 1). Such processes are present in all physiographic domains; even features and events that occur on the continental shelf and in deep-sea areas may have catastrophic effects on large areas in coastal environments [5]. Coastal areas are highly populated around the world and are the sites of most megacities [6]. A growing population, currently approximately 2.4 billion people, lives within 100 km of the coast (oceanconference.un.org, accessed on 28 May 2020). The expansion of urban coastal areas and coastal and offshore industries (e.g., communications, energy, and mineral extraction) has greatly increased the exposure and risk of large subaerial and submarine infrastructure. Contrary to the social perception, submarine areas host different and active features that present frequent geohazards.

**Figure 1.** Sketch showing the main offshore geological hazards. Inspired by [7].

Avoidance of hazardous areas is not always possible. Therefore, to reduce the risk of an area, conceptualized as the product of hazards and vulnerability, it is imperative to reduce vulnerability and obtain better knowledge of hazardous processes through new approaches and the development of new techniques and tools [8,9].

The field of offshore geohazards is wide because it covers different topic areas, such as identification and mapping of the hazards, risks, vulnerability, predictions, and warnings; covering all of these topics is beyond the scope of this paper. Specifically, this paper aims to provide the summary of the course and progress in the research of geoscientists regarding the main offshore hazardous geological processes, namely: seismicity, slope sedimentary instabilities, submarine volcanism, fluid flow processes, bottom currents, and tsunamis. Some geological settings, where interactions arise from cascading effects and thus create multihazard scenarios, are also shown. We analyze the key characteristics and main shortcomings of those geohazards, enabling the identification of present and future directions for marine geosciences focusing on offshore geohazards.

#### **2. Definition and Classification of Offshore Geological Hazards**

There are several definitions of geological hazards (e.g., [10–12]). As a synthesis, geological hazards can be considered all phenomena or conditions (on land and offshore), natural or induced by human activity, that can produce damage and that geological information can be used to predict, prevent, or correct.

The classification of offshore geological hazards may vary depending on their implications. From an engineering point of view, they are generally classified based on the problems they can generate during the exploration, installation, and operation of structures [13]. In contrast, marine geoscientists are more interested in understanding the features and processes that define a hazard; therefore, they classify them according to their causes.

Offshore geohazards arise from different geomorphological and geological features that produce scenarios in which diverse processes may act alone or in combination with others, triggering a chain of events. Morphological characteristics, such as relief (negative or positive) or overstepped slopes, may be indicators of processes that can generate hazards; however, depending on the activity planned, they may be considered a threat themselves.

The most widespread offshore geohazards are seismicity, slope instabilities, submarine volcanism, and processes related to fluid flow and bottom currents (e.g., [5,7,9,14]; https:// niwa.co.nz/natural-hazards/marine-geological-hazards;https://marineboard.eu/marinegeohazards-blue-economy, accessed on 3 October 2020) (Figure 1). Tsunamis deserve special mention because they are a secondary hazard derived from or generated by another event, especially earthquakes, slope instabilities, or volcanic eruptions ([5,15,16], among others). Additionally, many of these processes are associated with intense submarine erosion that may be responsible for the general topography and microtopography of the seafloor (e.g., [17]). Seismicity has a direct impact on the ground due to vibrations that affect infrastructure and buildings. Moreover, seismicity can also produce indirect effects, such as liquefaction and slope instabilities [18]. Earthquakes are related to the presence of seismogenic faults.

Submarine slope instabilities, resulting in products collectively referred to as landslides, mass-transport deposits, or mass-transport complexes, are capable of damaging infrastructure resting on or fixed to the seafloor, such as vertical foundations, communication cables, and pipelines, due to the associated impact, dragging, excessive burial, or undermining effects [19,20]. Active submarine volcanoes are significant geological hazards because of their violent and explosive eruptions and related earthquakes, collapses of their summits (i.e., caldera) and fluid emissions; moreover, volcanic activity can trigger secondary hazards such as tsunamis and landslides [5].

Fluid flow processes, such as seepage of light hydrocarbons, migration of overpressurized muddy fluids forming volcanoes and diapirs (Figure 1), and gas hydrate formation and breakdown, constitute a main type of potential geohazard [7,21]. They

are commonly a threat to navigation and offshore infrastructure, during both installation and operation, because they can cause damage or uncontrolled release of gas that in turn induces explosions or landslides.

The persistent action of bottom currents over the seafloor may create large areas affected by erosion and scouring and highly active seafloor conditions (Figure 1). Similarly, reworked or transported marine sediment is subsequently deposited, forming mounds in areas with high sedimentation rates. Bottom currents may also affect the sedimentological and geotechnical properties of seafloor and subsurface sediments, affecting their stability [22]. Thus, bottom currents may represent a hazardous process to subsurface and seafloor installations and infrastructure crossing the water column, because surface and intermediate currents can induce stress on them.

#### *Human Activities in Submarine Environments*

Human activities on the seafloor have increased sharply since the last half of the last century, accompanied by significant technological advances (Figure 2). Therefore, the understanding of offshore geological hazardous processes is fundamental for seafloor management. The main important offshore activities potentially exposed to offshore geohazards are as follows:


**Figure 2.** Some of the main human activities in deep-sea areas (>200 m water depth) in the NE Atlantic and western Mediterranean. The figure was used as idea and base to create a new one including more information from [25].

#### **3. Some Prehistorical and Historical Cases of Offshore Geohazard Events**

Offshore geological hazard events are recognized in all cultures and in all seas and oceans but are most common on highly active continental margins associated with tectonic plate boundaries (https://www.emdat.be/, accessed on 28 October 2020). Table 1 presents a small percentage of them, with only some of the widely known case events. All these disasters generated by prehistoric and historic marine geological events affected coastal populations, infrastructure, and the environment, highlighting the great vulnerability of these areas as well as the scarce knowledge regarding the triggering mechanisms of the events in particular or of the hazards in general. They have also revealed that active geological features do not recognize political frontiers and that hazard assessment must cross national boundaries. Furthermore, their occurrences have revealed that even small events may have a large impact on intensively exploited coasts (e.g., industry or tourism).


**Table 1.** Some prehistorical and historical cases of marine geohazard events.


**Table 1.** *Cont.*

## **4. Offshore Geohazards and Their Main Key Questions**

*4.1. Tectonic Earthquakes: Seismogenic Faults*

Internal geodynamic processes constitute the main engine that determines the presentday Earth's configuration. Most of the plate boundaries are located in offshore areas where seismically active faults constitute a principal marine threat [50,51]. In addition, some of them continue onshore, providing a good chance for direct observations (e.g., San Andreas Fault, [52]; the Alpine Fault, [53]). Progressive and continuous plate motion is accommodated by deformation along the plate boundaries where most of the active faults are located [54,55]. Some faults may undergo creep and are aseismic [56]. However, in sectors where two large fault blocks become coupled by asperities on the fault surface, elastic deformation, and stresses may increase, reaching the strength of the rocks [57]. The above factors drive sudden slip, producing an earthquake and, as a result, seafloor shaking, liquefaction, and permanent deformation [58].

Thrust faults related to subduction zones, including those associated with the Ring of Fire surrounding the Pacific Ocean and in the Indian Ocean, produce the most intense earthquakes [59–61] (Figure 3a). In these regions, although deep seismicity occurs, the most devastating earthquakes are located at shallow depths close to the coastlines. Examples include the Alaska (1964, Mw 9.2), Sumatra (2004, Mw 9.0–9.3.1), Chile (20101960, Mw 8.89.5), and Japan (2011, Mw 9) earthquakes and the related tsunamis [62] (https://earthquake.usgs.gov/, accessed on 11 September 2020). In addition, transcurrent faults can accumulate high stresses, eventually resulting in major strike-slip earthquakes, such as in Cape Mendocino in California (1992, Mw, 7.2) [63]. However, such pure strikeslip events generally do not produce vertical displacements of the seafloor in flat areas, although they may be significant when steep slopes are affected [64]. Regardless, seafloor deformation occurs in transpressional (e.g., Shackleton fracture zone, Figure 3b, [65,66]) and transtensional faults (e.g., Incrisis-Al Idrissi faults, Figure 3c, [29]). A particular setting occurs at the tips of such faults, where vertical displacement may trigger tsunamis [67]. Seismogenic normal faults are relatively scarce in marine environments and are most likely related to the isostatic response of the Earth's crust to ice loads close to coastal areas [68].

**Figure 3.** Examples of active seismogenic faults: (**a**) Subduction earthquakes in convergent plate margins have the highest magnitudes and produce seafloor shaking, submarine landslides, and tsunamis; (**b**) oblique-view perspective of the Shackleton fracture zone (Antarctica) that constitutes the transpressive sinistral active margin between the Scotia and Antarctic plates. Bathymetric features are determined by the permanent deformation of the seafloor; more details on this fracture zone can be found in [65,66]; (**c**) high-resolution seismic profile of the Incrisis and Al Idrisi sinistral fault zones in the central Alboran Sea, western Mediterranean. The tectonic activity of the Al Idrisi fault zone, now mostly covered by recent sediments, has been transferred to the new Incrisis fault zone, where the active faults affect the most recent sediments; more details on these fault zones in [29].

Detailed knowledge of fault features and their seismic characteristics is essential to preventing the effects of these geological hazards. The main target is to constrain the fault geometry, kinematics, dynamics, and seismic behavior to determine the related maximum magnitude of the seismic events and their recurrence interval. The numerical modelling of seafloor deformation and the propagation of seismic waves constitute one of the main tools to establish the direct impact of seismic waves on coastal populations and to determine the potential to produce submarine landslides and tsunamis. The determination of these key fault features and seismicity is developed either by coring, which allows the study of deposits related to the main events [69], or by geophysical methods. Seismological observations in areas surrounding active seismic zones are generally far from land, and the coverage of active seismic zones is limited, decreasing the accuracy of the locations of active seismogenic faults. Ocean bottom seismometers (OBSs) can record seismicity over long periods, but there are large delays in data recovery. Seismic reflection and acoustic techniques highlight the geometry of faults and improve the analysis of their activity. However, present-day standard techniques fail to produce a complete and detailed image of faults.

To improve the accurate analysis of the location and geometry of active faults, seafloor mapping of large parts of the oceans with new multibeam and sonar equipment is mandatory. The study of the continuity of faults at depth will require new seismic acquisition and data processing techniques, including 3D seismic methods for complex areas reaching depths of up to 12 km for crustal faults. The accuracy of seismological observations will increase with the installation of seafloor seismological observatories in wired networks to

provide real-time data. In addition, a denser network of OBSs will be necessary to achieve good coverage in active regions [50]. The characterization of fault behavior over long periods is also highlighted as new insight. In this sense, the improvement of submarine geomorphological indexes will likewise contribute to better identification and analysis of the active faults. Shallow coring techniques and detailed high-resolution seismic parametric profiles will improve the analysis of the geometry and age of the related deposits. These observations will help to determine the fault slip and the recurrence periods of the main submarine palaeoearthquake records. Deep coring techniques may enhance the significance of fluids in fault activity. In addition, the study of marine faults with offshore to onshore continuity will be essential for collecting direct observations of fault zones (fault surfaces, fault striations, fault gouges, and related deposits) and for determining their kinematics. Geodetic networks (global positioning system (GPS) and high-precision levelling (HPL)) and satellite-based Earth observations will also be mandatory to measure their present-day activity and constrain isostatic rebound in glacial margins close to the coastlines.

The integration of these new data will improve the probabilistic seismic hazard models in offshore tectonically active areas and the evaluation of their importance as secondary triggering factors of other hazards, such as slope sedimentary instabilities and tsunamis.

#### *4.2. Submarine Slope Instabilities*

Sedimentary instabilities are common processes in all submarine environments, where the largest slope instabilities in the Earth occur [70] (Figures 1 and 4). They may be classified according to different approaches, such as mechanical behavior, particle support mechanisms, sediment concentrations, and longitudinal changes in their deposits, or according to the relationship between source areas, dimensions, and geometries of deposits [71–76]. Based on the mechanical properties and rheology of the processes, two main groups can be defined: (i) slides/slumps/spreads and (ii) gravity flows. These two groups, with important differences in their pre- and post-failure behavior, occur in all physiographic environments and are efficient transporters of sediment, organic carbon, nutrients, and pollutants [77–81]. They are scale-invariant processes that range greatly in size from the meter scale to many km across (Figure 4).

Slides/slumps are movements of sediment or rock along a surface of rupture that develops in a layer with low shear strength or a weak layer [82]. They are elastoplastic movements that include translational and rotational movements (Figure 4a–c). Spreads, the submarine characterization of which has increased during the last decade, are sediment or blocks of consolidated sediment moving over liquefied underlying material and not a basal shear plane [83]. Depending on the mechanical behavior or the energy available, all these mentioned movements of sediments may evolve into a sediment flow, but a flow may also develop directly if the sediment is completely remolded.

A wide range of flow types can occur because of the interplay of rheology, grain size composition, and concentration (Figure 4). Flows in general have viscoplastic behavior and can be divided into cohesive (e.g., mud flows and debris flows) (Figure 4d) and non-cohesive flows (grain flows), depending on the amount of fine-grained matrix [84]. One type of cohesionless flow involving large volumes of failing masses is debris/rock avalanches. Usually, such failures originate from deep rotational failures on high-gradient slopes and in volcanic environments [85]. Turbidity currents are a type of Newtonian flow in which fluid turbulence is key to supporting the sediment and keeping it in suspension. Turbidity currents can transport up to hundreds of cubic kilometers of sediment at high velocities (up to 19 m/s) over thousands of kilometers [33,86] (Figure 4c).

**Figure 4.** Examples of submarine slope instabilities: (**a**) multibeam bathymetry displayed in the Gebra Valley area, Bransfield Basin (Antarctic); more details regarding this valley can be found in [87]. The dashed lines show the cross-sections of seismic records in b, c, and d; (**b**) Airgun seismic record displaying the repeated large-scale slope failure events that were responsible for the cut-and-fill features forming the Gebra Valley; (**c**,**d**) parametric seismic records showing the different slope instabilities affecting the external margins of the Gebra Valley; modified from [88].

Field monitoring of turbidity currents has increased in recent years. However, these processes and failures are still primarily recorded in nature, and most of the knowledge acquired is through the interpretation of the resulting geomorphology and the study of the final deposits. This situation leaves unanswered key questions and uncertainties about all the mechanisms involved in sedimentary instabilities, which need further study. These questions will have to be addressed at different scales through repeated very highresolution bathymetric surveys, high-resolution (and 3D) seismic surveys, new direct deep-sea monitoring and mobility sensors, in situ geotechnical tests, and experimental and numerical models. Better field observations and models will help to achieve an improved understanding of numerous aspects, such as the rates of seafloor changes, the

role of preconditioning factors, the impact of triggers, how rapidly a slope failure can develop, and the volumes of sediment involved and reworked. Improved knowledge of the governing mechanisms, evolution, and transformation of these submarine sedimentary instabilities will be crucial to understanding the hazard they represent [89–94]. Their modelling will also help to assess their consequences.

When landslide hazard assessment is considered, concepts such as distribution, time, and magnitude must be considered [95–97]. Regional inventories and magnitude-frequency relationships, including events triggered over a long period of time or almost instantaneously, could provide critical information [98]. Nevertheless, accurate chronological constraints (ideally combining biostratigraphic and radiometric techniques) will be essential for hazard evaluation.

#### *4.3. Submarine Volcanism*

Volcanoes are vents in the Earth's crust through which molten rock, hot rock fragments, and gases can erupt. Magma can rise along conduits to the surface, forming lava that either continuously flows out or shoots upward. Furthermore, the lava can break into pieces that are thrown into the air or into the sea due to decompression of the gases it contains (https://www.britannica.com/science/volcano, accessed on 3 September 2020).

From a marine point of view, volcanic eruptions affecting the sea water column can be grouped into four basic types (Figure 5):


**Figure 5.** Main types of eruptions examples: (**a**) Historical image of surtseyan eruption plumes from Surtsey (Iceland) in 1963 (image from National Oceanic and Atmospheric Administration, from https://commons.wikimedia.org/wiki/File: Surtsey\_eruption\_2.jpg, accessed on 03 June 2021); (**b**) emplacement of the 2007 lava delta at Stromboli (image from Rolf Cosar, downloaded from https://commons.wikimedia.org/wiki/File:Stromboli-Lavadelta-2007-03-10.JPG, accessed on 03 June 2021); more details on the 2007 submarine eruption at Stromboli can be found in [107]; (**c**) pyroclastic flow occurred on 16 December 1902, during a volcanic eruption of Mt. Pelée in Martinique (image from A. Lacroix, https://commons.wikimedia.org/wiki/File:Photographie\_du\_nuage\_noir\_du\_16\_d%C3%A9cembre\_1902\_lors\_ de\_l%27%C3%A9ruption\_de\_la\_Montagne\_Pel%C3%A9e\_en\_Martinique.jpg, accessed on 3 June 2021); (**d**) floating volcaniclastic materials and gas emissions during the Tagoro volcanic eruption offshore El Hierro in 2011 (image courtesy of Eugenio Fraile, IEO); (**e**) 3D bathymetric map of the Tagoro volcano showing the main morphological characteristics; more details on this submarine eruption can be found in [108,109]; (**f**) 3D simplified sketch of the main types of eruptions affecting insular volcanoes.

Another important hazardous phenomenon associated with volcanic eruptions that can become a hazard is gas emissions. The most common volcanic gases are water vapor, carbon dioxide, carbon monoxide, sulfur dioxide, and hydrogen sulfide; to a lesser extent, methane, hydrogen, helium, nitrogen, hydrogen chloride, or hydrogen fluoride can also be emitted (https://www.britannica.com/science/volcano, accessed on 1 October 2020). These gas emissions can cause loss of life due to suffocation if they reach the surface and can kill fauna present in the environment surrounding the underwater eruption. Such eruptions can also induce strong acidification of seawater, resulting in the subsequent loss of habitats around the eruption. For example, in the eruption of the Tagoro volcano, the pH dropped to 5, and the partial pressure of dissolved carbon dioxide increased almost 1000 times [108,110].

Despite the large number and volume of submarine volcanic eruptions, our understanding of such processes and associated landforms is still limited, especially compared with their subaerial counterparts ([111–113] and reference therein). Many concepts are still based on the interpretation of ancient deposits and on theory (e.g., [114] and references

therein). However, the growing availability of detailed digital elevation models (also used to depict seafloor changes associated with eruptions through repeated surveys) integrated with hydroacoustic monitoring and in situ observations of volcanic settings will exponentially increase their detection (considering that volcanic eruptions only occasionally reach the sea surface) and our knowledge of submarine eruptive processes (e.g., [115–120]).

A challenge-based study will provide knowledge to understand the processes that take place in the evolution of a submarine volcano at different depths. The main effort should focus on monitoring these processes using a variety of instrumentation (including on-land seismometers and marine stations with OBSs, hydrophones, pressure sensors, CTD (conductivity, temperature, and depth) instruments, and geochemical parameter sensors to control emissions) to allow study of the eruptive pulses and the content of emissions. Some of this instrumentation will be able to be connected by optical cables to laboratories onshore for online monitoring (e.g., [121]), and profiles can be made with a towed oceanographic rosette (tow-yo). Another challenge will be understanding the changes in seafloor morphology through time-lapse high-resolution bathymetry surveys, taking advantage of the use of autonomous underwater vehicles (AUVs) for deep volcanoes.

## *4.4. Fluid Flow Processes*

Seepage is a global process that occurs in different geodynamic contexts in both active and passive continental margins. Generally, this process includes the leakage of hydrocarbons (particularly methane as both dissolved and free gas), water, and/or sediment [21,122] (Figure 6). The gas in shallow marine sediments [123,124] is mainly composed of methane, and its origin is attributed to either biogenic or thermogenic processes. The escape of fluid from the sediment may occur as micro-seeps or as sudden violent escapes (cold seeps), producing diverse types of morphologies on the sea floor (Figure 6a) or in the subsurface [21]. Some features have positive relief (e.g., mounds, methane-derived authigenic carbonate, gas hydrates, mud volcanoes), and others have negative topographies on the seafloor (e.g., pockmarks, collapses) (Figure 6). Gas can migrate through unlithified sediments along bedding planes, faults, and fractures (Figure 6b) driven by buoyancy forces and pressure gradients [125,126]. Glacial-isostatic and tectonic events may reactivate fractures and faults, producing temporal variability in spatially heterogeneous fluid flow [127].

**Figure 6.** *Cont*.

**Figure 6.** Examples of hazardous features related to fluid flow processes: (**a**) At the top: very high-resolution seismic profile (3.5 kHz) in the Ria de Vigo showing an acoustically blank zone related to the presence of shallow gas accumulations near the seabed (less than 1 meter below the seabed). It also shows the presence of pockmarks and the main paths of gas escape. In the middle: multibeam echosounder images from Ría de Vigo displaying depressions related to gas escapes (pockmarks) and mounds formed by the accumulation of debris (mud and bivalve shells) from mussel rafts. At the bottom: methane-gas bubbles escaping from the sandy and muddy seabed of the Ría de Vigo; (**b**) mud diapirs related to the compressive regime (left figure) and to listric faults (right figure); (**c**) high-resolution bathymetric map of the northeastern Gulf of Mexico showing salt diapirs (rounded positive structures) piercing the seafloor and their associated seep anomalies associated with oil and gas seepages, slides and slumps, and gas hydrates to a lesser extent (images from https://www.boem.gov/oil-gas-energy/mapping-and-data/ map-gallery/boem-northern-gulf-mexico-deepwater-bathymetry-grid-3d; accessed on 15 September 2020; https://www.boem.gov/oil-gas-energy/mapping-and-data/mapgallery/seismic-water-bottom-anomalies-map-gallery, accessed on 16 September 2020; https://metadata.boem.gov/geospatial/WBA\_Metadata.xml, accessed on 16 September 2020); more details on the northeastern Gulf of Mexico salt diapirs in [128,129]). Enclosed seismic records show: a detailed view of salt-related normal faults acting as drainage pathways for deep fluids and lateral slides associated with overburden salt diapirs in a 2D seismic line located in the Gulf of Lion, western Mediterranean Sea.

The mud fluidization and degassing processes associated with overpressure contribute to the formation of mud volcanoes. This fluidization is mainly due to the overpressure generated by tectonic stresses or by lithostatic pressure in regions with high sedimentation rates [130]. These types of structures are one of the most important methane sources in the hydrosphere and atmosphere [131–133]. Diapirs are gravitational/tectonic structures (Figure 6b) produced mainly by salt, clays, or a mix of these lithologies that, currently, have little consideration in submarine hazard models. These intrusive bodies of relatively low density tend to migrate upward, deforming and piercing the overlying sedimentary sequences. They can appear with other fluid migration structures, such as mud volcanoes or pockmarks. In an extensive context, the development of diapirs commonly occurs close to deep listric faults (Figure 6b) that act as escape migration routes. The formation of mud diapirs is more frequent in compressive settings, as highly over-pressurized diapiric material (e.g., [134]) can move upward from subsurface depths up to 3–4 km to the seafloor. In various tectonic regimes, diapir activity may also trigger slope instabilities (Figure 6c) because the deformation and elevation of these structures favor seafloor oversteepening and seismicity [135].

Gas hydrate is an ice-like crystalline solid form of water and low molecular weight gas (e.g., methane, ethane, and carbon dioxide) [136]. Methane hydrates can form at any depth where the geothermal conditions are colder than the hydrate stability curve. This occurs in the upper hundred meters of marine sediment at water depths greater than 500 m. Nevertheless, certain conditions, such as the presence of saline pore waters or clays, can inhibit gas hydrate formation, while a high fluid flux can promote gas hydrate formation [137,138]. Seismic reflection techniques are used to determine the areal extent of gas hydrates in marine locations mainly by the identification of bottom-simulating reflectors (BSRs) (e.g., [139]). However, gas hydrates have been recovered from sites without BSRs. This is because the saturation of methane hydrate in the pore space must exceed approximately 40% for the seismic velocity to be altered significantly enough to generate a BSR [140]. Gas hydrates may be a significant hazard because they alter sea floor sediment stability and can lead to collapses and landslides that may trigger tsunamis [141–145], and their breakdown and release to the water column and atmosphere may have a strong influence on the environment and climate [146].

The presence of all these fluid flow features usually denotes subsurface hydraulic activity, over-pressurization, fluidization, and degassing processes, as well as sudden fluid (gas and/or liquid) release that may produce gas explosions, slope sedimentary instabilities, and an uplifting/subsiding seafloor [122,147–149] (Figure 6c). These processes can have major impacts on seabed infrastructures and on those requiring piles that are driven into the seafloor. Therefore, it will be necessary to extend systematic investigations to identify the locations of fluid dynamic processes in areas where their activity remains unknown currently. Thus, heat flow studies will need to be increased in order to detect and map new subseafloor marine fluid flows and understand their regimes. High-resolution 3D seismic surveys will also allow an accurate acoustic characterization and distribution assessment of the different fluid dynamic features and definition of their origins. They also have the potential to document and characterize in more detail the different types and timing of deformation patterns in areas close to diapirs and related mud volcanoes, with the goal of accurately determining the timing of fluid flow processes. Additionally, studies on microseismicity would allow the detection of fluid injection. Monitoring of fluid flows should also increase in active cold seeps; systematic sediment and gas sampling for biogeochemical analysis will aid the understanding of the general physical and geochemical characteristics of the escaping gas. Likewise, improved numerical models of gas hydrate formation, stability, quantification, and role in the shear strength of the host sediment will lead to progress in understanding the impact of gas hydrates on safety and seafloor stability.

## *4.5. Bottom Currents*

The term "bottom currents" refers to all persistent currents flowing near the seabed that resuspend, transport, and/or control sediment deposition [150]. Although bottom currents are semipermanent features, they are characterized by high variability over a range of time scales (from daily to geological timescales; [151]). They may occur in shelf, slope (Figure 7) and deep basin settings. In shelf settings, wind and tidal forcings are common and produce different hydraulic regimes (e.g., tide- and current-dominated regimes). In deep-water settings, bottom currents are mostly related to thermohaline circulation. These currents are driven by density gradients (e.g., the North Atlantic Deep Water) and typically flow subparallel to the bathymetric contours with velocities of 1–20 cm/s [152]. In particular settings (e.g., narrow gateways), bottom currents strongly intensify, reaching velocities of 50–300 cm/s (e.g., the Mediterranean outflow water that spills over the Gibraltar sill, e.g., [153,154]). In deep-water settings, submarine canyons can be swept by focused tidal currents with up and down flows and velocities of 25–50 cm/s [142]. Moreover, an increasing number of studies (e.g., [142,155–158]) have highlighted that several intermittent oceanographic processes can affect the seabed: eddies, gyres, helical flows, benthic storms, cascading dense shelf water, internal waves, and currents related to extreme, cyclonic, and tsunami waves.

As most of these processes can produce complex flow conditions, an increase in the velocity of bottom currents and additional shear stress may produce considerable sediment resuspension and seabed erosion [155] (Figure 7f–g). All these observations highlight that a variety of oceanographic processes are able to exert a significant impact on shaping the seafloor when bottom currents are active for a prolonged period of time (Figure 7a). At small spatial scales, they generate various erosional and depositional bedforms, ranging in size from centimeters to kilometers (see, e.g., [151,159]) (Figure 7c,f–h) and whose identification is particularly relevant for geohazard assessments of seabed infrastructure. At a larger scale, bottom currents with persistent activity on a geological time scale (e.g., thermohaline-induced currents) may form regionally extensive contourite depositional systems [158,160–163], including a variety of depositional elements (contourite drifts) (Figure 7a) and erosional elements (contourite channels, furrows, moats, and erosive terraces).

Evaluating the action of bottom currents is crucial for hazard assessment because intense seabed erosion may locally favor slope instability [164] (Figures 7a–e,i). Moreover, migrating bedforms (e.g., sand waves) and erosion (Figure 7f–h) can have major impacts on seabed infrastructure [7]. This can be extremely relevant in narrow straits swept by powerful tidal currents and by internal waves (e.g., the Messina Strait) that can create a dangerous setting for submarine cables and pipelines [165]. Other crucial areas are canyons swept by strong tidal bottom currents, topographic highs (e.g., seamounts and ridges) where bottom currents interact with topography [163,166–168], areas affected by tidal forcing and associated internal waves [169], areas of local upwelling [170,171], seasonal fluctuations in the main circulation pattern [172], or areas of sinking dense water [173], which may trigger slope sedimentary instabilities (e.g., [174,175]). Finally, it is noteworthy that contourite deposits can be prone to becoming unstable (e.g., [176]), as several predisposing factors (e.g., mounded morphology on steep slopes and the low shear strength related to high sedimentation rates) may favor slope failures [22] (Figure 7i).

**Figure 7.** Examples of hazardous features related to bottom current action of the Mediterranean outflow water around Iberia: (**a**) Sines contourite drift (Portugal margin, NE Atlantic) is affected by slide scars on the seafloor and failure surfaces in the subsurface. Additionally, upslope-migrating sediment waves occur in the near-surface sediments; (**b**) details of slide scars and potential failure surfaces; (**c**) details of sediment waves; (**d**) details of the bathymetric map with the track of the seismic profile shown in (**a**). The red dot marks the headscarp of a landslide scar; (**e**) detail of the slope gradient map; (**f**–**h**) ROV images showing bedforms in the continental slope of the Gulf of Cadiz (NE Atlantic): (**f**) Indurated muddy outcrop (grey in color) indicates intense current flow erosion; the mudstone is covered partially by sand with ripples. The two laser lines are separated 50 cm; (**g**) erosive furrows excavated on the muddy seafloor, covered partially by sandy sediment with starved rectilinear to sinuous asymmetric ripples. This starvation allows for the exposure of the mudstone surface over which the ripples are moving; (**h**) sinuous sand wave with superimposed linguoid to sinuous asymmetrical ripples on the stoss side and rectilinear to sinuous ripples in the trough area. The two laser lines are separated 50 cm; (**i**) diagram showing the intrinsic and hazardous properties of contourite deposits that may contribute to slope failures. Figure 7a–e,i was adapted with permission from [158]. Copyright 2019, Elsevier.

Recognition of the potential hazard of deep-water bottom currents is increasing because new and large seafloor areas of contourites have intensively eroded during the last 15 years due to mobile seafloors and slope sedimentary instabilities, which have been mapped (e.g., [150,155,158,160,162,177]). The assessment of the role of bottom current activity as a hazardous process is challenging for geoscientists due in part to the need to establish a dialogue with physicist oceanographers (e.g., [162]); as it is necessary to define flow conditions, induced bed shear stress and effects on morpho-sedimentary processes affecting the seabed, as well as their evolution over time. Thus, a multidisciplinary approach, including oceanographic, morphologic, sedimentologic, stratigraphic, and geotechnical studies, should be used. This will require multidisciplinary surveys and the integration of complex datasets, including oceanographic data (CTDs, acoustic Doppler current profilers (ADCPs), and transmissometers to measure the water properties and velocity not only at the near-bottom but also throughout the water column), multibeam bathymetry data, sub-bottom profile data, seafloor samples, sediment cores, and remotely operated vehicle (ROV) videos. Data integration enabled characterization of the different oceanographic processes, their interactions and timing, and their influences on the near-bottom flows acting on the seafloor. AUV surveys will also be required for high-resolution geophysical surveys in deep-water environments. Repeated bathymetric surveys and seafloor observatory systems will be required to define seabed evolution over time. Finally, hydrosedimentary modelling will be very helpful in assessing seabed changes and bed shear stress over a defined time period (e.g., the lifetime of the infrastructure).

## *4.6. Tsunamis*

The main mechanisms for tsunami generation are earthquakes caused by seismogenic fault movement (96% of events), slope sedimentary instabilities, and volcanic eruptions (Table S1). In addition, with the release of large volumes of gas from seafloor sediment, atmospheric disturbances (meteotsunamis), or even cosmic impacts can also produce tsunamis [178] (Figure 8). Finally, anthropogenically induced submarine slope sedimentary instabilities have triggered local tsunamis [179]. To perform tsunami hazard assessment, three main components of the phenomena must be addressed: the generating mechanism, wave propagation in the open sea, and coastal inundation.

Seismotectonic tsunamis are triggered by the coseismic vertical displacement of the seafloor impacted by an earthquake and the transmission of this movement to the water column [180,181] (see Video S1 in the Supplementary Materials). Normal fault movement causes the water masses to sink toward the formed depression, generating an initial large sine wave at the surface of water mass (Figure 8a). In contrast, reverse faults move the seafloor upward, and the water column is pushed upward (Figure 8b), forming an initial large crest wave on the sea surface. Tsunami wave generation by seismotectonics is controlled by the rupture velocity (mostly slow velocities), fault type, slip and average vertical displacement, width, length, and segmentation of the rupture zone of the fault in the seafloor [181].

Submarine slope instabilities (such as slumps, slides, debris/rock avalanches, debris flows) generally involve large volumes of sediments and rocks, and the associated movement within the water body generates a dipole-like water wave that can eventually generate major tsunami waves [15,182–184]. The initial shape of the tsunami wave is defined by a depression–uplift pair in the water surface (Figure 8c). The water depression is due to the sudden sediment vacuum that occurs at the slide scar (Figure 8c, number 1), which becomes occupied by sea water, and the uplift (Figure 8c, number 2) is due to the pressure force exerted upwards by the fast-moving slide material. Several aspects of submarine slope instabilities as tsunami sources are currently being discussed, namely, the rheology (because it influences the deformation of the sliding material during their runout) [84,85], acceleration, and volume. These aspects are considered to be key factors controlling the geometry (depression and uplift) of the generated tsunami wave [185].

**Figure 8.** Sketch of the main tsunamigenic sources of geological origin: (**a**) Normal fault activity; (**b**) reverse fault activity; (**c**) submarine landslide. Numbers 1 and 2 refer to the time-sequenced development of the initial tsunami shape in the water surface; (**d**) collapse flanks of volcanoes; (**e**) volcanic explosion related to eruptions and caldera collapses.

Volcanic eruptions can also induce subaerial and submarine slides, slumps, debris/rock avalanches, or debris flows on the flanks of volcanoes that, in turn, can produce tsunami waves (Figure 8d) and they can also generate caldera explosions that can cause the complete collapse of the edifice [186–188]. These explosions produce waves (Figure 8e) that are generated first by the explosion itself and then by the sinking of the volcanically mobilized material.

Tsunami waves move in all directions from the source area, affecting the entire water column [189]. Their initial propagation, especially their direction, is conditioned by the geometric and deformation characteristics of the source structure. Tsunami waves are initially characterized by their large wavelengths (tens or hundreds of kilometers), small heights in the open sea (on the scale of centimeters), and large velocities. The wave velocity is greater in deeper waters (700 km/h at depths > 4000 m), and when depth decreases, the velocity also decreases to 30–50 km/h at the coast. Simultaneously, the wavelength decreases and the wave increases in height to balance the kinetic energy with potential energy (e.g., [190]).

Tsunamis can have significant impacts on coastal communities, depending on regional and local bathymetry and coastal geomorphology variability [191,192]. The occurrence of reefs, human infrastructure, the geometry of the coastline and beaches, and the presence of bays, estuaries or deltas at river mouths can influence the size, appearance, and impact of tsunamis when they arrive at the coast. Typically, a tsunami reaches the coast as a series of successive crests and valleys, sometimes separated by several or tens of minutes, and can reach the coast as a rapid flood, more rarely as a wall of water, or, sometimes, as an initial withdrawal of the sea (e.g., [178]). Thus, the destruction caused by a tsunami on the coast can be very different at relatively short distances. The long wavelength of tsunamis gives them more momentum such that they can flood areas hundreds of meters and even kilometers from the coast. The maximum height above sea level that a tsunami reaches on the coast is known as the runup and mostly ranges from 1 m to 30 m, with extreme heights > 500 m, as in the case of the tsunami that occurred in Lituya Bay [184], when the coast is very close to the source area or where the coastal geomorphology amplifies the tsunami effects (resonant effects) (e.g., [181]).

The challenges in tsunami hazard research should focus on several aspects. First, an accurate definition of tsunamigenic sources is needed because it will help reduce the uncertainty in the triggering mechanism and will be important for studying the events themselves (e.g., the frequency of reoccurrence, potential areas to be impacted). Other aspects should include establishing their recurrence intervals by identifying past events in sediment cores (palaeoearthquakes, palaeoslides, and palaeotsunamis), and applying in situ measurements plus long-term monitoring. For seismotectonic sources, faults have to be described in terms of their tectonic style, dynamics, present-day activity, fault zone geometry, fault offsets of sedimentary units, and fault surface. For slope sedimentary instability-generated tsunamis, knowledge of the seafloor geometry, slope failure processes, and their early post-failure evolution is fundamental to determining their triggering potential. To analyze these concerns, high- and ultrahigh-resolution bathymetric data and 3D seismic reflection profiles, in situ seismicity measurements and observations, long-term monitoring, and longer sediment cores will be fundamental.

Additional challenges that will also be important for studying tsunami events and their impacts include the development of increasingly realistic mathematical models of the tsunami generation process, propagation through the water masses, and the impact on the coastal zones and the establishment of tsunami early warning systems (TEWSs). This work will include increasing the model resolution, developing more efficient and faster than realtime (FTRT) codes, and using future exa-scale computational architectures. Probabilistic tsunami hazard analyses (PTHAs) will have to be conducted in different areas of the world at global, regional, and even local scales with the aim of understanding tsunami hazards and developing tsunami risk reduction activities. PTHA increases the knowledge of the potential tsunamigenic threats at different scales by estimating the probability of exceeding specific levels of tsunami metrics, such as the maximum inundation height or runup within a certain period of time, around determined locations. Furthermore, probabilistic tsunami forecasting (PTF) attempts to address the uncertainty in tsunami forecasts by formulating a probability density function (PDF). The use of PTF in the context of rapid hazard assessment and in TEWSs is also a major challenge.

#### **5. Scenarios with Multiple Geological Hazards**

Following the above arguments and characteristics, the understanding of the different geohazard factors also needs to recognize the distribution of the main hazardous features, how they can interact, and their potential to generate cascading events. The most common of this type of event comprises an earthquake that triggers a landslide, both of which can produce a tsunami. Additionally, bottom currents can scour an overstepped seafloor, thereby reducing the shear strength of the unaffected sediments upslope and leading to their failure, forming a landslide that may produce a tsunami. Furthermore, the breakdown of sub-bottom gas hydrates can increase the pore pressure of the sediment bearing the

released gas, which may lead to tsunamigenic slope sedimentary instability. Despite the highly scattered distribution of these factors along continental margins, they commonly coincide in certain specific environments or geological contexts, which should be monitored by the scientific community. Diverse settings, such as fjords, active river prodeltas, canyonfan systems, subduction areas, or even high-latitude open slopes, may be critical. Among them, three settings are highlighted for their multiple hazardous features.

#### *5.1. Tectonic Indentation Areas*

Tectonic indentation areas in marine settings are significant because the tectonic structures developed in a framework of continental collision. The consequences of this tectonic activity are the presence of source areas that can produce multiple geohazards, such as earthquakes, sedimentary instabilities and, to a lesser extent, tsunamis. The continental indentation structures in marine areas occur in convergent continental margins related to plate corners (Taiwan, [193]) and accretionary wedges (Manila Trench, [194]) and in areas of early continental collision, such as the westernmost Mediterranean (Alboran Sea, [195]; Aguilas Arc in the Gulf of Vera, [196]). Particularly, the central Alboran Sea and the Aguilas Arc/Gulf of Vera (Figure 9) are key areas for understanding the link between indentation and geological hazards in a land–marine transition context; this is because although both exhibit similar hazardous features (seismic faults and slope instability deposits), their degree of development is different. Continental indentation influences the tectonics of the adjacent oceanic areas, as occurs during the northward indentation of the Arabian plate in Eurasia, which determines the westward motion of the Anatolian Block related to the development of the Aegean Sea [197]. The indenter blocks are generally bounded by lateral seismogenic strike-slip faults that permit displacement during the process of collision [198].

**Figure 9.** Example of multiple geological hazards in a tectonic indentation area in the western Mediterranean: (**a**) Geologic map showing the main tectonic features of the westernmost Mediterranean and the tectonic indentation zones in the central Alboran and Aguilas Arc in the Gulf of Vera; (**b**) detailed tectonic features of the Central Alboran Basin (for more details see [195]) and areas of generalized slope instability deposits (red areas). Legend: the red arrow indicates the direction of tectonic indentation and the red shading indicates the location of tectonic indentation.

Indentation structures simultaneously develop fault sets, folds, and block tilting, which can generate submarine slope sedimentary instabilities. These structures require integrated analysis by 3D analogue modelling, which has improved over time from early models [198] to recent models [199]. Future research will need to determine the stage of development of tectonic indentations. Additionally, the future development of new generations of numerical modelling is required. Another key question to address in the study of marine indentation zones is the tsunamigenic potential of strike-slip-related faults. In general, this type of fault is not considered tsunamigenic because it does not

significantly displace the seafloor. However, new data in the central Alboran Sea contradict this theory and indicate the need to investigate other strike-slip faults in similar geological frameworks [195].

## *5.2. Canyon Heads Close to Coast*

Submarine canyons, especially their shallower parts, are commonly very active geomorphological features that should be highlighted because of their association with multiple hazards (e.g., [200–202]). In general, canyons are located on the edge between a continental shelf and the continental slope, but some excavate the shelf to the point that their heads are only a few hundred meters from coastal towns, for example, the Garrucha (Figure 10a,e) and Gioia canyons in the Mediterranean [97–198,203,204] and the Congo and Capbreton canyons in the eastern Atlantic [19,205]. In these scenarios, changes can occur due to interactions among coastal processes (deposition and erosion) (e.g., [205]), river discharge (e.g., [206]), oceanographic processes (e.g., [207]), and seismicity related to tectonic processes (e.g., [208]). These activities can also produce favorable conditions for sedimentary instabilities (Figure 10a–d), which may produce tsunamis. Likewise, canyon heads close to the coast strongly influence tsunami propagation and runup (e.g., [209]).

**Figure 10.** Example of multiple geological hazards in the Garrucha Canyon, SW Mediterranean: (**a**) Bathymetric map displaying a canyon head affected by intense gullying of its two main tributaries; (**b**,**c**) ROV (remotely operated vehicle) images displaying slope failures affecting the canyon walls; (**d**) seismic profile illustrating the occurrence of mass-transport processes that contribute to the erosion of the canyon walls; (**e**) photograph of the Garrucha port, which is located at the canyon head, where the sedimentary instability processes that contribute to canyon-head retrogradation affect its pier.

Therefore, key scientific issues to be addressed include obtaining a better understanding of each submarine and coastal geological process and the oceanographic and climatic processes that govern the retrogradation, incision, and enlargement of canyon heads (e.g., [210]). To achieve this goal, detailed 2D, 3D, and 4D geomorphological visualization will be needed. This work should be carried out not only at the head of the canyons but also in the adjacent areas of the continental shelf, open continental slope, and infralittoral zone (Figure 10a). This visualization will allow us to typify and map with high precision the different morphological elements (both erosive and depositional) of the integrated canyon-head-margin system [101]. Within the canyon heads, the morphometric, chronostratigraphic, sedimentological, and geotechnical characterizations of submarine slope instabilities will be crucial; the integrated results will allow us to estimate the recurrence of events and to assess and model the potential canyon-head stability [204]. New insights will be fundamental to establishing the spatiotemporal relationship between slope sedimentary instabilities, tectonics, and oceanography and to defining areas that may be prone to failure (e.g., [211]).

#### *5.3. Volcanic Islands*

Volcanic islands and their submarine portions merge multiple geohazards, mainly associated with volcanic eruptions, flank collapses, slope instabilities (Figures 5a–e and 11), and associated tsunamis, although strong volcano-tectonic subsidence [212], retrogressive erosion at canyon heads [213,214], and earthquake swarms (e.g., [215]) deserve special attention for hazard assessment. Flank collapses and slope instabilities ([216] and reference therein) (Figures 5g and 11) represent a common hazardous process during the evolution of many insular volcanoes, which are often able to mobilize volumes up to thousands of cubic kilometers (e.g., [186,217–219]). For instance, the 100–400 × 106 m3 flank collapse affecting Anak Krakatau in 2018 generated a tsunami with a runup of up to 13 m along the Sunda Strait. Despite the high tsunamigenic potential associated with these large-scale events, their hazard is relatively low because they have recurrence times on the order of thousands of years. In contrast, small- and medium-sized slope instabilities affecting active volcanic flanks are more hazardous because they have markedly shorter recurrence times and are able to generate local but devastating tsunamis [220,221]. One of the best examples is recognizable at Stromboli Island, where five tsunamigenic landslides over just the last century have been reconstructed [222]. In the case of highly explosive eruptions, the entrance of pyroclastic currents into the sea can also generate tsunami waves or travel (their upper and dilute parts) over the sea for distances of tens of kilometers before impacting surrounding coastal communities, as described during the 1883 Krakatoa eruption [223].

Considering the multiple often closely related hazards affecting volcanic islands, the key scientific recommendation for an effective hazard assessment in such areas should include the use of an integrated and multidisciplinary approach encompassing both the submarine and subaerial flanks of the island. High-resolution mapping will be fundamental to understanding the variability in volcanic edifices and associated landforms (Figure 11) and to performing systematic parametrization to provide insights into the complex interplay between the volcanic, tectonic, erosive-depositional, and eustatic processes (for shallowwater areas) that control the genesis of volcanic islands ([113] and references therein). This mapping (Figure 11) will also be the basis for planning more detailed surveys with seismic methods, ROV dives, and seafloor sampling and for successive bathymetric comparisons aimed at understanding what occurs during eruptive crises. In this regard, the availability of multiple time-lapse bathymetric surveys has been proven to be a very effective tool for monitoring seafloor changes associated with volcanic and/or failure events occurring in both shallow water (e.g., [107,115,224]) and deep water (e.g., [225,226]). In particular, the integration of repeated bathymetric surveys with acoustic monitoring and/or ROV dives will increase our ability to detect and understand eruption dynamics in submarine environments [120,227,228].

**Figure 11.** Main characteristics of Terceira Island (Azores; North Atlantic Ocean) and the bathymetry surrounding the island. The map also shows the area where lava balloons clusters and volcanic ash were observed at the sea surface during the 1998–2001 Serreta eruptions; more details on these eruptions can found in [118,229].

## **6. Conclusions: A Distinctive Multidisciplinary Approach to Study Offshore Geological Hazards**

Offshore geological hazards include convulsive and persistent geological processes and are mainly represented by seismicity, slope sedimentary instabilities, submarine volcanism, fluid flows, and bottom currents; tsunamis are also mentioned because they are commonly a secondary hazard generated mostly by earthquakes, slope instabilities, or volcanic eruptions. They can occur in any domain or environment in the oceans and seas and represent a real and serious threat to society, the economy, and the environment.

Despite the progress in data acquisition and the establishment of evolutionary models for the different hazardous features, each dedicated section has identified knowledge gaps and how these gaps can be addressed. We also note that hazardous processes can interact and potentially generate cascading events.

This review establishes that the challenges for improving outcomes in offshore geohazard research can be addressed with multidisciplinary approach studies. This approach requires cross-disciplinary research to bring together multiscale analysis, mapping, direct deep-sea observations and testing, and modelling in scenarios with individual, but mainly multiple geohazards. This approach will lead to multicriteria decisions for understanding hazardous processes and their causative factors.

A qualitative step in the multiscale analysis involves the acquisition of long-term geological records, such as seismic profiles with different degrees of resolution and penetration and geophysical data (e.g., magnetometer and gravimeter data) (Figure 12), all acquired simultaneously in surveys using emerging technology and applying advanced tools for processing and geophysical modelling (Figure 12). The long-term records also provide the opportunity to study seismic profiles and sediment cores (the longer the better) (Figure 12) to improve our understanding of the magnitude and frequency of hazardous processes. Advancing techniques in sediment core analysis and age dating will contribute to reducing the uncertainty between stratigraphic correlations and increase their temporal resolution. This will facilitate the attainment of more accurate information about the sediment age and the recurrence interval of hazardous events. The success of these observations from these conventional but continuously advancing techniques will be closely tied to seafloor mapping, direct deep-sea observations and testing, and modelling.

**Figure 12.** A distinctive methodological approach to study offshore geological hazards. The main future directions for studying offshore geohazards will be the implementation of multiscale and multidisciplinary approaches joining conventional and emerging tools for monitoring, mapping, direct observations, in situ testing, and modelling. Scenarios can be affected by multiple hazardous features, some in a land-marine transition context, and the integration of offshore and onshore observations is essential. The figure was used as idea and base to create a new one including more information from [79].

Mapping is not a new approach but is still needed to fill the existing gaps along many continental margins and basinal plains and to provide higher-resolution seafloor maps from shallow to deep-sea areas, thus providing more details on the occurrence of hazardous features. This task is mandatory for taking the next step in geological hazard assessment. The capacity to perform repeated high-resolution multibeam bathymetric surveys is also a very effective tool for monitoring seafloor changes, and it has to be implemented for a better understanding of hazardous processes. In this sense, the use of AUVs with multibeam sonars and sub-bottom profilers in repeated surveys and of ROVs for direct observations are essential to map the active hazardous seafloor features with maximum resolution, even in deep sea environments (Figure 12); the temporal resolution of their mapping will be important to report their dynamic evolution. Seismic record acquisition also plays an important role in mapping prior hazardous processes. In this sense, a greater use of 3D seismic data is expected to offer new and unprecedented sub-bottom geomorphologic information, enabling an accurate delimitation and characterization of active structures, especially in complex geological settings. However, the success of future efforts to map hazardous seafloor features will require confident geomorphological models and mapping standards for the correct understanding and recognition of features. There is still a long way to go before the scientific community reaches an agreement on standards for marine geohazard mapping, but this is a requirement for future multiple hazard maps and catalogues.

Direct observations and testing are also technical challenges because both are linked to the development of new sensors, techniques, protocols, and infrastructures, such as seafloor observatories (Figure 12). The development of new seafloor observatory systems with capabilities and facilities for remote and real-time recording over long periods of time will lead to qualitative advances. Direct observations of active structures and measurements of smaller-scale, but highly recurrent, events will enhance our understanding of larger processes and also provide important data for small-scale models. Moreover, direct observations and measurements from the water column via the optimization of

mooring systems, CTDs, ADCPs, and transmissometers will be essential to understand the physics of the environment (e.g., bottom currents, turbidity, etc.). Additionally, direct seafloor monitoring by seafloor network systems, including OBSs with longer standing periods than are available in the present and DART (Deep-ocean Assessment and Reporting of Tsunamis) buoy systems, or the use of submarine cables to detect earthquakes, will provide better seismic data with which to define active faults, surface ruptures, volcanic activity, and tsunamis (Figure 12). Direct seafloor observations and measurements, together with inland seismic stations, will allow us to define seismogenic and aseismic faults and estimate realistic peak ground acceleration (PGA) values, which depend on the epicentral distance and earthquake magnitude (Figure 12). Seismic loading is critical to defining the factor of safety (F) and the susceptibility to slope failure of different seafloor areas. Integrated seafloor and inland monitoring will be a key element to increase the reliability and timeliness of the information used by early warning systems (Figure 12). However, future generations of AUVs and ROVs, which have become more widely accessible to the scientific community and easier to manage, will also provide new opportunities for in situ sampling of sediments, monitoring the rate of seafloor mobility, fluid seepage characterization, and measuring in situ geotechnical parameters. The in situ geotechnical properties (Figure 12) involve the goal of obtaining contact measurements (seafloor and sub-bottom) using advanced static (e.g., cone penetrometer, pressuremeter, heat flow), dynamic (e.g., XBP) and combined systems. The measurement of the dynamic effect of seismic events or other cyclic sources, such as storm waves and internal waves on the sediment, especially on the pore pressure, is another rarely performed technical approach that will need to be enhanced (e.g., dynamic simple shear tests). Special in situ tests can reveal the real effect of high sedimentation rates or fluid flow dynamics (gas emissions) on the sediment and its geomechanical characteristics. In addition, in situ geotechnical property measurements will be closely tied to the assessment of potential seafloor stability by the application of probabilistic methods, for which GIS is an adequate and very powerful tool.

The success of future efforts for an effective seismic hazard assessment will be reliant on new and better geological models taking advantage of developments in artificial intelligence. Consequently, the future of the field of geohazards is coupled to enhanced computational capabilities. This is also because this field will face a massive volume of datasets (i.e., the so-called big data problem). Datasets will be generated from new hullmounted and towed methodological instruments as well as autonomous and permanent observational systems recording multiple hazard datasets with higher temporal and spatial resolution. This means that in the future, a common but key challenge will be the ability to efficiently manage and analyze (also in real time) massive data; therefore, the need to train geoscientists and build the capacity to operate advanced computational systems are needed. Massive data linked to the advances in artificial intelligence will open new lines of research for the use of advanced deep learning and machine learning algorithms, for example, for automatic detection and classification of different variables (e.g., slope gradients, roughness, backscatter signal amplitudes, grain size, density) involved in the identification of geohazard features. The development of complex neural networks trained to detect variables of interest (e.g., in seismograms, bathymetries, cores) will offer an important advancement both qualitatively (elements could be detected automatically that could go unnoticed by the most trained analyst) and quantitatively (the analysis of multiple datasets as well as their spatiotemporal relationships will increase exponentially). Such advances will lead to a much greater understanding of hazardous processes and will have significant effect on the probabilistic methods for assessing geological hazards with more robust models from which early warning systems will benefit.

Furthermore, there is also the need to enhance multidisciplinary studies in the geological scenarios with multiple hazardous processes that can interact and generate cascading events. These scenarios can be affected by hazardous features that are connected from sea to land. Therefore, the integration of on-land information, e.g., Global Navigation Satellite System (GNSS), high-precision levelling, seismicity monitoring, multiple geophysical

datasets (including magnetometer, gravimeter, magnetotelluric (MT), electrical resistivity tomography (ERT), and HPL data, and field, drone, and satellite Earth observations), with submarine results is critical to realizing the correct assessment of hazards (Figure 12). In this sense, scientific and technical coordination of the research community working on subaerial and submarine hazardous structures with different datasets is a major task that is still in its infancy, and reinforcement is required in the future.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/oceans2020023/s1, Table S1: Okada's parameters [3] describing the two fault segments involving of the Horseshoe and Marques de Pombal faults as source for the tsunami modelling. Video S1: Video of tsunmai simulation in the Gulf of Cadiz.

**Author Contributions:** Conceptualization, G.E. and D.C. (David Casas); Writing—original draft: all authors based on their geological hazard expertise; Introduction by G.E. and D.C.; Human activities and historical cases and by G.E. and B.A.; Definition and classification of geohazards by G.E. and D.C. (David Casas); Tectonic earthquakes by J.G.-Z., F.E., L.G.-C., A.M.-T., and V.T.-S.; Landslides by D.C. (David Casas), M.A.-Z., and J.I.-G.; Volcanism by D.C. (Daniele Casalbore), O.S.-G., and J.-T.V.; Tsunamis by J.-T.V., J.M., J.M.G.-V., and D.P.; Fluid-flow by S.G.-G., D.D., M.C.F.-P., and G.E.; Bottom currents by E.M., M.T., and C.J.; Scenarios involving multiple geological hazards by G.E., D.C. (David Casas)., D.C. (Daniele Casalbore), F.E., and J.G.-Z.; Conclusions by G.E., D.C. (David Casas), J.N., J.V., and M.Y.; Integration—review-reediting: G.E. and D.C. (David Casas). All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was supported by the following scientific projects: FAUCES (MINECO/FEDER, CMT2015-65461-C2-R), DAMAGE (MINECO/FEDER, CGL2016-80687-R); B-RNM-301-UGR18, P18- RT-3275 and RNM148 from Junta de Andalucía; GRACE (Eurofleets: EFP\_SEA02-024); VULCANO (CTM2012-36317); RIGEL (IEO); MEGAFLOW (MEC/FEDER, RTI2018-096064-B-C21); and ChEESE (UE (H2020-INFRAEDI-2018-1, PROJECT ID: 823844). This study was also carried out in the framework of the IGCP 640–S4LIDE. In addition, the research was also conducted in collaboration with the Andalusian PAIDI Research Group RNM 328 (Coastal and Marine Geology and Geophysics).

**Acknowledgments:** The ICM-CSIC authors acknowledge Severo Ochoa funding from the Spanish government through the "Severo Ochoa Centre of Excellence" accreditation (CEX2019-000928-S). We also thank IHS-Markit for the Kingdom Suite educational license. This article is dedicated to the memory of Albert Figueras, who devoted the last years of his life as director of the UTM-CSIC.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Tamiji Yamamoto 1,\*,†, Kaori Orimoto 1, Satoshi Asaoka 1, Hironori Yamamoto <sup>2</sup> and Shin-ichi Onodera <sup>3</sup>**


**Abstract:** Although the water quality in Hiroshima Bay has improved due to government measures, nutrient reduction has sharply decreased fisheries production. The law was revised in 2015, where the nutrient effluents from the sewage treatment plants were relaxed, yet no increase in fishery production was observed. Herein, we investigate the distribution of C, N, S, and P within Hiroshima Bay. Material loads from land and oyster farming activity influenced the C and S distributions in the bay sediments, respectively. Natural denitrification caused N reduction in areas by the river mouths and the landlocked areas whose sediments are reductive. The P content was high in the areas under aerobic conditions, suggesting metal oxide-bound P contributes to P accumulation. However, it was low in the areas with reducing conditions, indicating P is released from the sediments when reacting with H2S. In such reductive sediments, liberated H2S also consumes dissolved oxygen causing hypoxia in the bottom layer. It was estimated that 0.28 km<sup>3</sup> of muddy sediment and 1.8 <sup>×</sup> <sup>10</sup><sup>5</sup> ton of P accumulated in Hiroshima Bay. There remains conflict between the 'Legacy of Eutrophication' in the sediment and 'Cultural Oligotrophication' in the surface water due to 40 years of reduction policies.

**Keywords:** carbon; eutrophic; Hiroshima Bay; phosphorus; nitrogen; sulfur; sediment

## **1. Introduction**

Estuaries function as filters/traps or removal devices for nitrogen (N) and phosphorus (P) to the offshore [1]. At the river mouth, flocculation of colloidal iron oxides occurs upon the mixing of riverine freshwater and seawater, and orthophosphate is trapped as Fe(III) (oxyhydr)oxides-P and sedimented on the bottom [2]. In addition to this chemical process, particulate matter tends to be sedimented at the innermost area of estuary due to estuarine circulation, which is induced by river discharge [3,4]. Phytoplankton grows by taking up nutrients supplied from land through rivers and other point sources and the sea. Phytoplankton will sink to the bottom as dead cells or feces of filter feeders along with other particulate matter. In the succession of mineralization on and in the sediments, organic-bound P can be liberated, and orthophosphate is taken up by phytoplankton and bacteria again in the water column. This is one main framework of the filtering or trapping function of P in estuarine systems.

However, the P cycle in the sediments and overlying water is much more complex. Iron phosphate is formed in rich Fe and phosphate conditions. Fe(III) (oxyhydr)oxides form at the oxic/anoxic boundary between sediment and overlying water [5]. During oxygen consumption via decomposition of organic matter supplied from the upper water column, Fe and P are released to the bottom of the water column. In anoxic sediments under non-sulfide conditions, such as lakes, ferrous phosphate may significantly contribute

**Citation:** Yamamoto, T.; Orimoto, K.; Asaoka, S.; Yamamoto, H.; Onodera, S.-i. A Conflict between the Legacy of Eutrophication and Cultural Oligotrophication in Hiroshima Bay. *Oceans* **2021**, *2*, 546–565. https:// doi.org/10.3390/oceans2030031

Academic Editor: Michael W. Lomas

Received: 30 December 2020 Accepted: 10 August 2021 Published: 16 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to P retention in the sediments [6]. However, the production of H2S in reducing sediments in estuaries counteracts the functioning of Fe(III) (oxyhydr)oxides in trapping P by forming iron sulphides (FeSx) instead [7]. Simultaneously, Fe(III) (oxyhydr)oxides-P may undergo reductive dissolution in the presence of S2<sup>−</sup> and release orthophosphate [8]. Subsequently, the increased liberation of orthophosphate, Fe2+ and Mn2+ over long periods forms authigenic phosphorites known as apatite [9–11] and vivianite [12–15]. Thus, in general, P is permanently removed from the water column and buried in the sediment, and not all sedimented P returns to the water column [16–18]. The importance of each diagenetic process is different in different areas and sites, depending on the depositional environment, sedimentation rate, and prevailing redox conditions [5].

N is removed in estuaries. The main pathway of N removal is denitrification, accounting for ~50% of N loaded through rivers [19–22], while anaerobic ammonium oxidation (anammox) and burial in sediments are negligible [23]. Nitrification can be the major nitrate source for denitrification, supplying >80% of the nitrate demand in the St. Lawrence Estuary [24,25]. When the bottom is always reductive with no nitrate supply, denitrification depends on riverine nitrate supply or oxic seawater intrusion, as observed in Lakes Shinji and Nakaumi, Japan [26]. Under anoxic conditions and during bottom water hypoxia, dissimilatory nitrate reduction to ammonium (DNRA) competes with denitrification, and it may recycle the N in the system and contribute to an increase in primary production [27–32]. This varies across different areas and sites depending on the level of eutrophication; the ratio of nitrate to organic matter is a primary factor which overcomes denitrification or DNRA in the system [33]. Furthermore, the presence of reduced sulfur in the sediments also influences denitrification, DNRA [34], and coupled nitrification-denitrification [35]. Furthermore, a larger fraction of nitrate can sometimes be retained in the system in the presence of sulfides [27,35]. Overall, estuaries can work naturally as N removal devices and phosphorus traps.

The Japanese government established the 'Law Concerning Special Measures for Conservation of the Environment of the Seto Inland Sea' in 1978, following its predecessor, the 'Law Concerning Tentative Measures for Conservation of the Environment of the Seto Inland Sea' in 1973 [36]. As part of the measures taken for the implementation of this law, reductions in P and N contents in discharge water from all workplaces/factories, which discharge >50 ton day−<sup>1</sup> of water, began in 1980 and 2005, respectively. Thus, the long-term reduction of nutrients (40 years for P and 15 years for N) has successfully contributed to improving the water quality, as evidenced by the apparent increase in water transparency [37]. However, such simple reductions of P and N do not render the desired effect to the Seto Inland Sea ecosystem because the response of the ecosystem often shows non-linearities characterized by hysteresis [38–41]. Fishery production in the Seto Inland Sea has sharply decreased with decreasing primary production since the late 1980s without tracing the same trajectory as it had traveled during the eutrophication period [41].

In Hiroshima Bay, which is in the western part of the Seto Inland Sea, oyster culture has been conducted intensively for more than 60 years, accounting for 60% of the total production in Japan [42]. However, the oyster production in the bay is in a crisis due to the lack of feed phytoplankton, and recently, it has been observed that the oyster larvae cannot survive due to the lack of small-sized phytoplankton (<5 μm in diameter) suitable for their feed [43]. In contrast, oysters, through their feeding activity, stimulate biogeochemical processes in the sediments by increasing the sedimentation of organic matter from the water column [44,45]. Therefore, the sediment quality is muddy and contains vast amounts of refractory organic matter which have been historically deposited from the oyster culture via both oyster feces and dead oyster meat, and organisms attached on the shells sink to the bottom. In the sediment of such reduced conditions, sulfide is produced along with the anaerobic decomposition of refractory organic matter. Therefore, acid volatile sulfide (AVS), water content (WC), and ignition loss (IL) of the sediment are higher in the Hiroshima Bay compared to the neighboring area with less aquaculture activity [46]. Oxygen depletion in the bottom layer (hypoxia) is still observed even after 40 years of

the P reduction measure. This can be attributed to the oxygen consumption by reductive substances such as sulfide produced in the sediment, not by freshly produced organic matter decomposition. This can be referred to as a 'Legacy of Eutrophication' or 'Legacy of Hypoxia', which has been referenced in the previous studies [47–50]. That is, H2S and other reductive matter produced in the anaerobic sediment, in which an abundance of refractory organic matter accumulated in the bottom sediments during the former eutrophication period, consume the bottom dissolved oxygen (DO) and liberate phosphate into the bottom water even several decades after the peak eutrophication period.

The fisheries production has been decreasing since the 1980s when the measure was implemented [41]. All kinds of bioproduction, not only seaweed culture, bivalve culture, and capture fisheries, have apparently decreased [41,51,52]; the fish catch and shellfish production in the Seto Inland Sea in 2014 was one-third and one-sixth of their peaks in 1980s and 1970s, respectively. Therefore, the central government revised the abovementioned law in 2015 [53]. The major revision was the relaxation of nutrient removal from sewage treatment water. However, no signs of improvement in any aquaculture production and capture fisheries were observed until 2021. The challenge in the Seto Inland Sea including Hiroshima Bay is a 'conflict' between the 'Legacy of Eutrophication' of the sediments and 'Cultural Oligotrophication' of the surface water by the excess nutrient reduction, taken as the measure to improve eutrophication. The sediment quality of the Seto Inland Sea, including those in Hiroshima Bay, have been surveyed thrice thus far by the government [54–56] with the standard methods [57]; however, there are few analyses and discussions in relation to fishery production.

In this study, we investigate the sediment quality in Hiroshima Bay, describe the differences in the distribution patterns of sediment carbon (C), N, P, and sulfur (S) contents, and discuss the possible causes of the distribution patterns of these elements in the sediments from physical, chemical, and biological process perspectives. A discussion on the conflict between the eutrophic sediment state and the oligotrophic surface water condition, which has been induced by the long-term nutrient reduction measures is also presented.

## **2. Materials and Methods**

## *2.1. Characteristics of Hiroshima Bay*

Hiroshima Bay is in the western part of the Seto Inland Sea, Japan. A strict nutrient reduction measure of P and N for the Seto Inland Sea has been mandated by the government since 1980 and 1995, respectively. Both the total phosphorus (TP) and total nitrogen (TN) discharge loads in 2014 were approximately 1/2 of their peak loads of in 1976 and 1994, respectively [52]. Oyster culture is intensively conducted in the bay, producing ~60% of the total production in Japan [42]. Although the nutrients transported from the land well cherished the oyster production through feed phytoplankton growth in the bay, the nutrient reduction measure strongly damaged the oyster production, resulting a decrease in oyster production of fresh meat to 1.6–2.0 × <sup>10</sup><sup>4</sup> ton yr−<sup>1</sup> in recent years from the peak production of 3 × 104 ton yr−<sup>1</sup> in 1986 [42]. While cultured oysters may suppress excess phytoplankton blooms by their feeding activity, the biodeposits deteriorate the sediment quality [44]. In addition to feces, pseudofeces, dead oysters, and other creatures grown on the raft culture system sink onto the bottom when farmers are handling the oysters, deteriorating the sediment quality. Thus, the sediment quality of the bay is poor compared to that of the neighboring area in terms of oxidation-reduction potential (ORP), IL, and AVS [46].

Hiroshima Bay has an area of 1043 km2, average depth of 25.8 m, volume of 26.9 × 109 m3, and catchment area of 3743 km2 [36]. It is geographically divided into two areas: the northern Hiroshima Bay (nHB), north of the narrow strait (Nasabi Strait) between Itsukushima Island and Etajima Island, and the southern Hiroshima Bay (sHB), south of the channel (Figure 1). The nHB area is 141.2 km2, excluding Etajima Bay and Kure Bay, and the remaining is that of sHB [58]. In nHB, pronounced stratification is established due to large freshwater input from the Ohta River system; the average total discharge including the other small rivers

is 7.5 × 106 m3 day−<sup>1</sup> and the contribution of the Ohta River is ~90% [58]. The freshwater residence time in nHB is estimated to be 27 d in the study [58]. Thus, the stratification is a 'halocline,' ultimately reducing vertical mixing [58], thereby contributing to the expansion of hypoxia in the bottom layer [59]. The freshwater from the Ohta River system entrains approximately seven times volume saline bottom water on average from the sHB [58], which may transport both suspended and resolved matter from the sHB. Meanwhile, stratification of the water column is not as strong in sHB, except the western part where there is freshwater input from the Nishiki River and the Oze River (Figure 1). Kure Bay and Etajima Bay are stagnant; therefore, hypoxia with DO concentration <2 mg L−<sup>1</sup> in the bottom water is observed in September in these bays in addition to the innermost area of nHB [59].

**Figure 1.** Map showing the monitoring lines of sediment thickness (dotted line) and core sampling stations (filled circle). 5–8 October 2010.

## *2.2. Field Observations*

Field observations were conducted during a research cruise at 14 stations in nHB and 18 stations in sHB during 5–8 October 2010 (Figure 1). The thickness of muddy sediments in Hiroshima Bay was monitored throughout the cruise at a ship speed of 4–5 knots, using a sounding device (SH-20, Senbon Electric Co., Ltd., Numazu, Japan) emitting two acoustic frequency beams (15 and 200 kHz) [60]. As the device is graduated in 10 cm intervals, we measured 1 cm between the tick by the eyes. The sediment thickness data was used to estimate the accumulated P amounts in the sediments in Section 2.4.

At every station, five sediment core samples were collected using a KK-type core sampler with an acrylic tube inner diameter of 4.8 cm. Three of the sediment cores were selected for analyses based on their appearance and lack of artificial disturbance. First, the overlying water temperature was measured with a bar thermometer. The temperature range was small, between 22.7 (Stations 24 and 32) and 25.8 ◦C (Stations 10 and 19) because the water column had already mixed vertically due to surface cooling. After measuring the overlying water temperature, the core samples were sliced at 5-cm intervals from the surface to 20 cm depth, and at 10-cm intervals below this until the depth that the sediment was collected with the device was reached. The slicing of the sediment core samples was rapidly performed on the ship deck after the sediments were collected at each station. Although the surface of the sediment may have partially oxidized during the slicing process, we aimed to not oxidize the sediments by storing each sediment sample in a gas-tight plastic container 7.8 cm in diameter and 4.6 cm tall and placed them in a dark, cool (0–5 ◦C) container filled with nitrogen gas until various analyses are conducted. Sediment core sample collection failed at Station 5 because the bottom was composed of sandy sediments.

#### *2.3. Chemical Analyses*

Immediately after sediment core sampling, the ORP was measured for the samples collected from 0–5 cm using an ORP meter (PS-112C, DKK-TOA Corp., Tokyo, Japan). Temperature was recorded using the ORP meter. The Eh value was calculated from the ORP and temperature values, where Eh (mV) = ORP (mV) + 206 − 0.7 (t − 25). The AVS concentration was determined with a detection tube (Hedrotech-S, Gastech Corp.) for the sediment samples collected from 0–5 and 5–10 cm [61]. By adding acid to the sediment sample, H2S gas was detected as the color of the reagent in the test tube changed (brown); H2S + Pb(CH3COO)2 (white) → PbS (brown) + 2CH3COOH. The detection limit of the detection tube was 0.0002 mg, and the Coefficient of variation (CV) was 5%. Here, AVS is an operationally defined property that is the H2S evolved from the acid leaching of sediment. A detail review on what compounds are measured by the acid leaching method of AVS is described in the literature [62]. Dividing sediment into aqueous phase and solid phase, the major dissolved species in AVS are HS−, H2S, and the aqueous FeS under the condition of reduced sediments with pH < 7. In the solid phase of sediment, FeS2 (pyrite) may not contribute to AVS. Notably, part of the FeS in the aqueous phase consists of nanoparticles of mackinawite due to sample handling usually through a 0.2 μm pore size filter. The hydrogen sulfide (H2S) concentration was determined with a detection tube (Kitagawa-type, Komyo Rikagaku, 200SA, 200SB, Kawasaki, Japan; detection limit 0.1 mg L−1, accuracy ± 15%) [63] for the interstitial water of sediment core samples collected from 0–5 and 5–10 cm after centrifugation at 3000 rpm for 15 min. The head space of the centrifugation tube was filled with nitrogen gas, and both AVS and H2S analyses were conducted within 3 min to minimize oxidation of the sample. This quick measurement procedure does not give significant results [64]. The detection tube was calibrated beforehand by H2S standard solution by dissolving an aliquot of Na2S9H2O (Nacalai Tesque, Kyoto, Japan) in 3% NaCl solution to correct a salinity error. The H2S oxidation rate is dependent on the initial concentration of H2S and salinity but not on water temperature and pH [64].

All sediment samples were stored in a cool (0–5 ◦C) dark box and brought back to the land laboratory. After thoroughly mixing each sample with a spatula, the mud wet density (g cm−3) was determined using a weighing balance. The sample dry weight was determined after drying the wet mud samples for 12 h at 110 ◦C in an oven. The WC (%) of each mud sample was determined as the difference between the wet and dry weights. The dried mud samples were combusted at 600 ◦C in a muffle furnace for 4 h and weighed again. IL was determined as the difference in weights before and after combustion [65].

After adding 2N HCl solution to a portion of the dried mud sample to remove inorganic carbon, samples were dried, and total organic carbon (TOC), TN, and total sulfur (TS) were determined using a CHN analyzer (CHNS/0 2400II, Perkin Elmer, Waltham, MA, USA). The TP was determined using an inductively coupled plasma-atomic emission spectroscopy (ICP-AES, Optima 7300DV, Perkin Elmer; The Ministry of the Environment of Japan [57]) with a detection limit of 0.1%. The measurements were conducted in duplicate in the first run, and then the second run was taken place in case the two values in the first run are apart. We adopted the average value of the two closest values.

## *2.4. Estimation of the Sediment P Content*

First, we created a depth contour map showing the thickness of muddy sediments. All the depths of the surface muddy sediment measured with the acoustic device were characterized by WC (average: 61.4% +/− standard deviation: 3.6%) and IL (average: 9.90% +/− standard deviation: 0.91%). Here, the data were selected for 18 stations where the demarcation between the soft, muddy surface sediment and the hard sediment below it was clear. The demarcation level was used to estimate the thickness of the surface muddy sediment in the other stations.

Then, the mud volume (m3) was calculated by multiplying the area (m2) of every 5-cmthick sample. A horizontal distribution map of each element concentration was produced at 5-cm depth intervals, and areas (m2) with the same concentration were summed. Then, total amount of P in the Hiroshima Bay sediments was estimated using the area (m2), average concentration (mg g<sup>−</sup>1), mud wet density (g cm−3), and thickness of each layer.

#### **3. Results**

#### *3.1. Sediment C, N, S, and P Contents*

Vertical profiles and horizontal distributions of TOC, TN, TS, and TP are illustrated in Figures 2–6. The vertical profiles of each element showed high values in the surface layer and decreased gradually with depth (Figure 2). Only the TS content was high in the depth of 10–20 cm at the Ohta River mouth (Station 27) and Kure Bay (Stations 31 and 32). The average values, standard variations, and coefficient of variations are summarized in Table 1. In the average value in the table, only TS has a peak value at 10–15 cm depth. Small variation in the vertical profile of TP values was clear (Figure 2), which is evidenced by the low CV (6–11%) as summarized in Table 1.


**Table 1.** Average (Avg) and standard deviation (Stdev) of TOC, TN, TS, and TP in the Hiroshima Bay sediments observed in 5–8 October 2010. Coefficient of variation (CV, %) is also shown.

TOC was high in nHB, particularly at the Ohta River mouth (Stations 27) with the highest value (>20 mg g dry<sup>−</sup>1) followed by Etajima Bay (Station 20) and Kure Bay (Station 32) with values ~15–20.0 mg g dry−<sup>1</sup> (Figures 2 and 3). In sHB, TOC was slightly high in the western part, where the Nishiki River water enters, with values of 10–15 mg g dry−<sup>1</sup> (Figure 3) which were high at the surface (Figure 2).

The highest TN values were observed in the central area, both in nHB (2.9 mg g dry−<sup>1</sup> at Station 22) and sHB (3.0 mg g dry−<sup>1</sup> at Stations 10) as shown in Figure 4. However, it showed lower values (1–2 mg g dry−1) below the surface (Figure 2). In contrast, the values at the Ohta River mouth (Station 27) and the Nishiki River mouth (Station 9) were

<1.5 mg g dry<sup>−</sup>1, as shown in Figures 2 and 4, which were different from the distribution of TOC, as shown in Figure 3.

The highest TS value in the surface 0–5 cm layer, 7.3 mg g dry−1, was observed at the mouth of Etajima Bay (Station 21) as shown in Figure 5; however, the concentration was low in the deeper layer at this station (not shown in Figure 2). The highest values (9.3–9.4 mg g dry−1) were detected at 10–20 cm depth at Station 20 in Etajima Bay (Figure 2). Except those high values in and around Etajima Bay and the western part of sHB (5.9 mg g dry−<sup>1</sup> at Station 9) in the surface layer, the TS concentration, including in the deeper layer, did not differ much between nHB (2.0–7.1 mg g dry−1) and sHB (1.1–5.9 mg g dry<sup>−</sup>1).

**Figure 2.** Vertical distributions of total organic carbon (TOC), total nitrogen (TN), total sulfur (TS), and total phosphorus (TP) at selected stations. 5–8 October 2010. The number in the figure indicates the station number of the selected stations; 27 (the Ohta River mouth), 9 (the Nishiki River mouth), 22 (the center of nHB), 10 (the center of sHB), 31, 32 (Kure Bay), and 19, 20 (Etajima Bay), respectively.

**Figure 3.** Spatial distributions of the total organic carbon content of muddy surface sediments (0–5 cm) in Hiroshima Bay in 5–8 October 2010.

**Figure 4.** Spatial distributions of the nitrogen content of muddy surface sediments (0–5 cm) in Hiroshima Bay in 5–8 October 2010.

**Figure 5.** Spatial distributions of the sulfur content of muddy surface sediments (0–5 cm) in Hiroshima Bay in 5–8 October 2010.

**Figure 6.** Spatial distributions of the phosphorus content of muddy surface sediments (0–5 cm) in Hiroshima Bay in 5–8 October 2010.

The spatial distribution of TP was somewhat homogeneous and was different from the other elements, with values of 0.37–0.67 mg g dry−<sup>1</sup> in both nHB and sHB (Figures 2 and 6). The values at the northeastern part of nHB along with the Ohta River mouth (Station 27) and the Nishiki River mouth (Station 9) were slightly low, as shown in Figure 6.

As shown in Figure 7, while the thickness of muddy sediments in nHB was 0.2–0.4 m at the mouths of the Ohta River and the Yahata River, it was 0.3–0.5 m thick in Kure Bay, where the water is stagnant due to the landlocked configuration (Figure 1). In sHB, the muddy sediment thickness was also high (0.4–0.5 m) in the southernmost area north of Yashiro Island and locally in western area. The volume of muddy sediments was estimated to be approximately 2.82 × <sup>10</sup><sup>8</sup> <sup>m</sup><sup>3</sup> for the entire HB. The amount of P accumulated in the Hiroshima Bay sediment, estimated using thickness data of the muddy sediment and the P content, was 1.8 × <sup>10</sup><sup>5</sup> ton. Accumulation processes will be discussed later.

**Figure 7.** Spatial distributions of the thickness of muddy sediments in Hiroshima Bay estimated for the observed data in 5–8 October 2010.

## *3.2. Subsidiary Parameters Relating to C, N, S and P Cycles*

Observed parameters (Eh, WC, IL, AVS, and H2S) for the soft, muddy surface sediment whose depths are different at every station, which relate to C, N, S and P cycles, are summarized in Table 2. The Eh in nHB was from −150 to −27 mV. The lowest area was Etajima Bay, with values ranging from −150 to −120 mV, indicating that sediment in Etajima Bay was under severe reducing conditions. In sHB, the sediment conditions were also generally reducing (−153 to −25 mV), except for some sandy areas (Stations 1 and 3: +23 and +45 mV, respectively). The WC ranged 69.8–82.5% in nHB and 46.6–81.8% in sHB. The IL values were high in nHB (10.9–12.7%), and particularly high in Etajima Bay, Kaita Bay, and Kure Bay, compared to those in sHB (4.8–11.3%). The AVS was higher in nHB, with values of 0.17–1.14 mg g−1, than in sHB (0.04–0.35 mg g−1). The highest AVS value was observed at the mouth of the Ohta River (Station 27). H2S in the interstitial water of the sediments was detected at most stations. The highest value (23 mg S L<sup>−</sup>1) was observed both in Kure Bay (Station 32) and the center of sHB (Station 13), and the second highest value (19 mg S L<sup>−</sup>1) was in Etajima Bay (Station 20).


**Table 2.** Summary of observed parameters (Eh, water content, ignition loss, acid volatile sulfide, and hydrogen sulfide) in the surface (0–5 cm) sediments of Hiroshima Bay. 5–8 October 2010. H2S is the concentration in the interstitial water of sediment. N and S mean northern and southern Hiroshima Bay, respectively.

#### **4. Discussion**

The C content was generally high at the mouths of Ohta River, Yahata River, and Seno River (Kaita Bay). This may be due to large inputs of organic matter from the cities of Hiroshima and Kure, which have large populations (Hiroshima: 1.2 million, Kure: 0.23 million). Even in sHB, the C content in the western part was high, which can be explained by the material load of the Nishiki River flowing through Iwakuni City (population: 0.14 million). One of the major sources of particulate organic matter is the direct load from the rivers. According to Kittiwanich et al. [66], the 10-years (1991–2000) average riverine particulate organic nitrogen (PON) and particulate phosphorus (PP) loadings are 14.5 mg N m−<sup>2</sup> day−<sup>1</sup> and 2.73 mg P m−<sup>2</sup> day−1, respectively for nHB, and 0.6 mg N m−<sup>2</sup> day−<sup>1</sup> and 0.27 mg P m−<sup>2</sup> day−<sup>1</sup> for sHB, respectively. Notably, other large PON and PP loads are those conveyed by the estuarine circulation, which is driven by the riverine discharge; those for nHB are 21.0 mg N m−<sup>2</sup> day−<sup>1</sup> and 2.59 mg P m−<sup>2</sup> day<sup>−</sup>1, respectively, and for sHB they are 1.0 mg N m−<sup>2</sup> day−<sup>1</sup> and -0.26 mg P m−<sup>2</sup> day−1, respectively. The minus PP value in sHB indicates PP was transported from sHB to outside the bay. These values indicate the material load by the estuarine circulation is comparative to the riverine direct N and P loads or more. Although the carbon loads from the rivers and the estuarine circulation were not estimated in their study, it is inferred that the same tendency in POC load being high in nHB and low in sHB.

The second largest source of particulate organic matter in Hiroshima Bay is biodeposits from oyster culture, which are intensively conducted particularly in nHB. According to the material budget calculation by Kittiwanich et al. [66], the amount of PON and PP produced as feces and pseudofeces by cultured oysters are 18.1 mg N m−<sup>2</sup> day−<sup>1</sup> and 2.09 mg P m−<sup>2</sup> day−<sup>1</sup> for nHB, respectively, and 1.3 mg N m−<sup>2</sup> day−<sup>1</sup> and 0.15 mg P m−<sup>2</sup> day−<sup>1</sup> for sHB, respectively. The particulate matter load by the biological process, particularly PON, is larger than those of direct riverine load both in nHB and sHB. The biodeposits from the cultured oysters contribute approximately 2/3 of the total sedimentation both in PON and PP in nHB. The remaining total biodeposits produced in the water column mainly comes from the feces of zooplankton. In nHB, including Etajima Bay, oyster farming is intensive; approximately 12,000 oyster rafts (of them, 1700 rafts in Etajima Bay) exist in nHB, equivalent to 5760 million individual oysters. According to Yamamoto et al. [44], 27.6 kg DW deposits raft−<sup>1</sup> day−<sup>1</sup> are estimated to be produced. Biodeposits of ~330 ton DW raft−<sup>1</sup> day−<sup>1</sup> sink down onto the seafloor of HB. This estimate does not include oyster meats, regardless if they were dead or alive, and other attached creatures, which may have fallen during handling by farmers, which may reflect almost the same amount of the feces and pseudofeces deposits.

Conversely, the distribution of N differed substantially from that of C, showing low values at the mouth of the Ohta River (Figure 4). Consequently, this area had high sediment C/N molar ratios (>11). In sHB, the C/N molar ratio was also high (12–23) in the western part near the Nishiki River than in the central part (5–11). The low N values in these river mouth areas can be explained by denitrification processes. The occurrence of denitrification is also evidenced by the low N/P molar ratio. The N/P value 5.6 recorded at the Ohta River mouth was much smaller than the Redfield ratio (16). It was reported that severe hypoxia (DO concentration <2 mg L−1) of the bottom layer is observed in the north-east part of nHB usually in September and recovers in October every year [59]. The hypoxia is limited in the innermost north end of nHB (Stations 27 and 28) where riverine input is received, and the area landlocked by islands with narrow connections to the neighboring areas (Stations 19, 20, 23, 26, 29, 30, 31, and 32). It is reported that the annual scale of hypoxia in the main nHB, except in the landlocked area, is dependent on the DO supply accompanied by the advective lateral flow driven by the estuarine circulation, although the seasonal recovery of hypoxia occurs due to ventilation by vertical mixing of the water column as the season progresses [59].

Long-term monitoring data on the bottom DO concentration as well as the other water quality data are available from the homepage of Hiroshima Fisheries and Marine Technology Center [67]. As the number of stations monitored by the prefectural institution is limited, here we use the data collected at two stations (Stations 17 and 21 in the original data set) whose locations are close to our sampling stations (Stations 22 and 31) for discussion. Therefore, in Figure 8, we used our station numbers for these sites. The data indicate that the bottom DO concentration at Station 31 in Kure Bay is much lower than that at Station 22 (Figure 8). This is due to the water stagnancy in Kure Bay, while DO is supplied by the estuarine circulation at Station 22 as mentioned above. There is no fixed long-term trend in the bottom DO concentrations even though the nutrient reduction measure has been implemented for 40 years. The most serious eutrophication period was during the 1970s–1980s and before in terms of the number of harmful algal blooms; the record high was 1974 (14 times yr<sup>−</sup>1) and the second record high was 1986 (10 times yr−1) in Hiroshima Bay [67]. However, no harmful algal blooms are observed in recent couple of years. This can be the 'Legacy of Eutrophication', which has been mentioned in previous studies [47–50]. That is, H2S and other reductive matter produced in the anaerobic sediment, in which an abundance of organic matter accumulated in the bottom sediments during eutrophication period, continue to consume the bottom DO and release phosphate into the bottom water, even several decades after the peak eutrophication period.

In Figure 8, we selected the lowest 10 values for DO and the highest 10 values for each nutrient species at Station 31 (Kure Bay) during 1988 to 2010 since the sampling interval is not regular before 1987. The bottom NH4-N concentration was high during June–July, except in August 1987; meanwhile, the bottom NO3-N was high during August– December with several exceptions (January 1989 and June 2004). This may imply that organic matter decomposition progresses in early summer (June–July) and is followed by nitrification (August–December). During the nitrification process followed by H2S generation, the DO consumption proceeds and forms hypoxia in September. Taking the difference between the NO3-N peak value to the following lowest value which occurred 1–2 months later, the denitrification rate is ~6.6 mg N m−<sup>2</sup> day−<sup>1</sup> for Station 31 on average over the three decades. Correspondingly, it is estimated to be 4.6 mg N m−<sup>2</sup> day−<sup>1</sup> for Station 22. This is only a first-order estimation, and the values may be overestimated because the decrease in NO3-N concentration is also caused by the other processes such as physical dilution and ammonification, etc. However, it is in the range of the reported values (0–23.2 mg N m−<sup>2</sup> day<sup>−</sup>1) determined by the acetylene inhibition method for the sediments collected at Station 22 in 1994–1995 [68], and close to the value estimated for the entire Hiroshima Bay using the data from April 1991 to March 1992 (7.9 mg N m−<sup>2</sup> day<sup>−</sup>1) [69] and slightly lower than the value estimated for the nHB using the data set of 1991–2000 (14 mg N m−<sup>2</sup> day−1) [70]. The latter two were those estimated by a material budget calculation using models. It is not so much different from the values estimated for the sediments of the eutrophic Stockholm archipelago, Baltic Sea (1.26–24.1 mg N m−<sup>2</sup> day<sup>−</sup>1, in the original study, 90–1723 μmol N m−<sup>2</sup> day<sup>−</sup>1) [50].

**Figure 8.** Temporal variations of bottom water (**a**) dissolved oxygen (DO), (**b**) PO4-P, (**c**) NH4-N and (**d**) NO3-N concentrations. The samples were collected at 1 m above the bottom. Data are cited from Hiroshima Fisheries and Marine Technology Center [67]. Month names were given to the top 10 annual peak values for each nutrient species (lowest 10 for DO) at Station 31 (Kure Bay) during 1988 to 2010 since the data set was not complete before 1987.

The S and P distribution patterns differed substantially from those of C and N. For the S and P cycles, we should consider whether the sediment conditions were aerobic or anaerobic, similar to N. Sulfate ions are plentiful in seawater and are the source of sulfate reduction under anaerobic conditions in addition to organic matter. In other words, H2S is the byproduct of sulfate reduction. As described above, the H2S concentration and AVS, in which H2S is contained, were high in Etajima Bay, where S originates not only from sulfate ions, but also from oyster feces. This may be the cause of the high S contents in the sediments of Etajima Bay (Figure 5). In Japan, the standard permissible values of AVS and H2S for aquatic life are 2 mg g dw−<sup>1</sup> and 0 mg L−1, respectively [65]. Compared

to the standards, all stations except Station 27 cleared the AVS standard, and Stations 1, 3, 8, 11, and 14 in sHB cleared the H2S standard. In comparison to the values between the samples from 0–5 cm depth and from 5–10 cm depth, the AVS was higher at 0–5 cm depth (0.34 ± 0.21 mg g dw−1) than 5–10 cm depth (0.21 ± 0.15 mg g dw−1). This may be attributed to the formation of FeS and other compounds, which can be detected as AVS [62]. Meanwhile, the H2S was higher in 5–10 cm (6.1 ± 5.4 mg L−1) than 0–5 cm (4.6 ± 5.6 mg L−1). This may indicate H2S dissolved in the pore water is diffused to the overlying water and causes the consumption of DO. The volume-based percentage content of AVS in TS and that of H2S in AVS were 7.3 ± 5.0% (0.8–28.6%) and 1.1 ± 1.3% (0.0–6.4%), respectively. Although the fraction of H2S in the sediment TS is somewhat small, it should be continuously produced as long as the vast organic matter remains in the sediments and sulfate is supplied from seawater.

The observed P content was comparatively low at the stations located at the mouth of rivers and stagnant areas (Stations 19, 20, 26, 27, 29, and 31 in Figure 6), contrary to the high C content at these stations (Figure 3). In these areas, material load is high either from rivers or oyster culture. The concentration of the elements in the sediment is determined by the balance between the sedimentation flux and the release of dissolved matter from the sediment due to decomposition. Therefore, although the organic matter sedimentation is high, the relatively low sediment P values implies a high benthic flux owing to anaerobic conditions. The sediment Eh was low (−150 to −120 mV) in these areas, as described before, indicating reducing conditions. H2S is contained in the sediments of all stations in nHB as described in Section 3.2. As shown in Figure 8, peak phosphate concentration in the bottom layer of Station 31 is usually observed in September, which coincides with the occurrence of hypoxia. This phenomenon is the same at Station 22. The phosphate under anaerobic conditions originates from metal oxide-bound P, such as FeOOH≡PO4 in the reaction with H2S, in the sediment. The production of H2S in the sediments will continue if both sulfate and the organic matter are present in the sediments, even if it is refractory. This is the 'Legacy of Eutrophication' as previously mentioned [47–50].

An estimated 1.8 × <sup>10</sup><sup>5</sup> ton of P accumulated in the Hiroshima Bay sediment, equivalent to half the amount consumed in a year in Japan (3.5 × 105 ton) [71]. Even if P is liberated yearly from metal oxide-bound-P, the amount liberated from the anoxic nHB sediments is estimated to be 26 ton P yr−<sup>1</sup> from the average concentration in the bottom as 2 μmol L−1, assuming the bottom water depth containing the concentration is 3 m. If we assume the P accumulated in the sediment changed the form to those liberated yearly, it will take 6900 years until the end of release of all the accumulated P. Although this estimation is too simple and rough, it is certain that the annual hypoxia will last long even if a harsher nutrient reduction measure is taken.

Additional measures are required to remediate such deteriorated anoxic sediments. Two engineering methods have been implemented to remediate deteriorated anoxic sediments. One simply removes and dampens the deteriorated anoxic sediments in the other site. Because of the salt content, we cannot use the sediments in agriculture as fertilizer or anywhere else on land. Therefore, the most suitable place to dampen was the shallow coastal area along the coastline in terms of transportation cost. This kind of reclamation was effective during the rapid economic growth period because there was demand for land for industries and houses. However, such reclamation was restricted by the tentative law in 1973 because it damages coastal ecosystems. Another way is to remediate the sediment quality by adsorbing or oxidizing the reductive substances using materials that have such functions. Traditionally, natural sand was used for capping deteriorated coastal sediments. However, the collection of natural sand was restricted by the law because it may damage the ecosystems of original sites where the sand is collected regardless of whether it is land, river, or sea. Furthermore, natural sand has no function to reduce the reductive substances. Thus, we have developed several functional materials from industrial byproducts, which are certificated by the Ministry of the Environment of Japan [72], and hot-air dried oyster shell (HACOS), steel-making slag, and granulated coal ash are those in which the present

authors' group has been deeply involved [73–78]. Their use in the future will be beneficial due to cost-effectiveness and the aspect of creating a recycle-oriented society.

Third, the recovery of the accumulated P as a resource can be an alternative. Currently, P is recovered from sewage at treatment plants and used in fertilizers and other products [79]. Abelson [80] indicated that phosphorus is a resource facing shortages on the Earth; therefore, techniques should be developed in the future to recover P that is accumulated in coastal sediments.

The Ministry of the Environment of Japan observed the sediment quality in the Seto Inland Sea, Japan, including Hiroshima Bay [54–56]. There were 46 stations in the first observation in 1982, and 23 in the second and third observations, which did not differ drastically from our number of sampling stations (32) in October 2010, showing an even distribution of the stations. The results are summarized along with our data from October 2010 in Table 3. The sampling depth in the observation by the Ministry of the Environment was the same as those of the first layer in our observations (0–5 cm), with data in deeper layer >5 cm lacking. Decreasing trends over time in TOC and TN were found in a report by Komai [81]; however, the discussion is limited. In contrast, only TP values in the present study were slightly higher than the average TP value from 2003. A possible cause is the trapping/adsorption of phosphate by iron oxides and manganese oxides that are formed under aerobic conditions. While the sediment P retention capacity is small when the bottom condition is anaerobic because of the reductive dissolution of phosphate detached from iron/manganese oxides [82,83], it increases via trapping by iron/manganese oxides as oxygen is supplied [84]. A record flooding during July 10–16 of 2010 [85], 3 months before our observations, might have increased the sediment P content.

**Table 3.** Comparison of TOC, TN, and TP values between those reported by Ministry of the Environment, Japan [54–56] and the present study (2010). The sediments were collected from the same depth in all samples (0–5 cm) in Hiroshima Bay. A similar number of sampling stations were included in the observations (46 in 1982, 23 in 1993 and 2003, and 32 in 2010), and were evenly distributed throughout the bay.


To remediate the eutrophic conditions, the P and N reduction measure has been in place for 40 years for the Seto Inland Sea including Hiroshima Bay. Consequently, the water transparency has increased drastically [37]; however, the production of commercially important fish, bivalves, and cultured seaweed has almost collapsed in the Seto Inland Sea [41,51,52]. We previously reported this could be because of 'cultural oligotrophication' [41,86]. Thus, we conclude that the Seto Inland Sea, including Hiroshima Bay, has been experiencing conflict between the oligotrophication of the surface water and the legacy of eutrophication of the bottom sediments. The recent decrease in the oyster production of Hiroshima Bay is a serious issue for farmers. The line length farmers use to hang oysters is 10 m. Therefore, the oysters are not affected by the hypoxia, which usually forms in the bottom layer. However, the anoxic sediment conditions are unsuitable for benthic animals, including bivalves, as their rearing environment. They should be an important feed for both benthic and pelagic fishes; for example, red sea bream, which is a highly regarded fish in Japan, prefers to feed on shrimp dwelling on the bottom. Thus, the hypoxia and sediment which contains deadly poisonous H2S can cause the collapse of the total ecosystem in addition to the oligotrophication in the surface water.

To alleviate the oligotrophic conditions of the surface water and increase the bioproduction, the Japanese government revised the law in 2015. The measure was to increase the sewage effluent load to a level at which natural and farming seaweeds, bivalves, and major fish species can sufficiently grow [53]. In Hyogo Prefecture, in the eastern part of the Seto Inland Sea, they relaxed the sewage treatment effluent during the Nori growing season in the winter by suppressing the nitrification and denitrification processes at sewage treatment plants [87]. However, observations have proven the effect was spatially limited because of the diffusion feature of seawater [87]. Computer simulations, which combine physical and biological processes, also showed limited effects [88]. Another cause of the limitation of the effect was a permissible limit in the maximum concentrations of nitrogen and phosphorus of the effluents they can release. The nutrient concentration was still regulated by the policy in 2015. Then, the central government revised the law again to expand the target sites to any workplace which emits phosphorus and nitrogen. The revised bill has recently passed the House of Representatives on 3 June 2021 [89].

Since the oyster rafts are extensively located in nHB, the supply of nutrients by increasing effluents from workplaces may not impact the offshore rafts. Therefore, we are attempting to develop a fertilizer to enhance the growth of farmed oysters [90], which can also be used for other bivalves, Nori, and other farming organisms. The most important issue that we must address is the concentration and time we should increase the effluents from workplaces, how many fertilizers, and when it should be applied. As the environmental conditions and type of fisheries conducted are different in each area, each prefectural government must devise a practical and effective plan for individual sea areas. For example, Hiroshima Prefecture is assigned to plan for Hiroshima Bay. When increasing nutrient load to maintain the oyster growth, the prefectural government should consider the oyster growing season is October to March. However, nutrient increases must not induce harmful algal blooms. Simultaneously, we must consider the remediation of anoxic sediments. Researchers are expected to propose some acceptable scientific perspectives using a sophisticated simulation model, which includes integrated physical-biogeochemical processes in both the water column and sediments. Our group has previously established a model that can reproduce the bottom hypoxia and remediation effects of byproducts by applying them to the sediments [91]. The other models we have been developed can be used to measure the oligotrophication of the surface water and oyster production [43,66]. We must discuss the simulation outputs in detail as a future prospect with different stakeholders.

#### **5. Conclusions**

The C, N, S, and P showed different spatial distributions within Hiroshima Bay. C was influenced more by material loads from land, whereas N reduced by denitrification in river mouths and landlocked areas with reduced sediments. S was affected by oyster farming activities, which are intensively conducted in, for example, nHB and Etajima Bay, and sulfate reduction that may occur under reducing conditions. P was relatively low in reduced sediments of landlocked areas (northeast of nHB, Kure Bay and Etajima Bay), where P is released under reduced conditions, and distributed evenly in the vertical profile

compared to the other elements. P was somewhat high compared to C and N referring to the Redfield ratio, suggesting that it can accumulate in sediments. Compared to the governmental monitoring data on sediment quality, the P value of the present study is higher than that of 2003, whereas decreasing trends in C and N were found. Approximately 1.8 × <sup>10</sup><sup>5</sup> tons of P is accumulated in 0.28 km<sup>3</sup> of muddy sediment in Hiroshima Bay, with a maximum thickness of 0.5 m.

The Seto Inland Sea, including Hiroshima Bay, is facing conflict between the 'Legacy of Eutrophication' in the bottom sediments and the 'Cultural Oligotrophication' of the surface water. Particularly, the lack of feed phytoplankton for farming oyster is a serious challenge for farmers in the Hiroshima Bay. In contrast, the hypoxia observed every September is unlikely to recover forever. A major cause is the biodeposits supplied from farming oysters. The oxygen consumption in the bottom layer is attributed to the reductive substances such as H2S. Recovery of the sediment quality is important for animals which are feed for both benthic and pelagic fish. Therefore, we should apply materials that can adsorb the reductive substances to recover the total ecosystem, in addition to adding nutrients to the surface water. A science-based plan with several options with computeraided perspectives, including cost performance, is warranted for the local government to facilitate proper decision-making.

**Author Contributions:** Conceptualization, T.Y. and S.-i.O.; methodology, T.Y. and H.Y.; software, H.Y.; validation, S.A.; formal analysis, K.O.; investigation, T.Y., K.O., S.A. and H.Y.; resources, T.Y.; data curation, S.A.; writing—original draft preparation, T.Y.; writing—review and editing, S.A. and H.Y.; visualization, K.O. and H.Y.; supervision, T.Y.; project administration, T.Y. and S.-i.O.; funding acquisition, S.-i.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by JSPS Grant-in-Aid for Scientific Research (KAKENHI) (A), Grant Number JP21241011.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available because they are collected in an university student's graduation project. The corresponding author who was the supervisor of the student keeps all the data in the forms of CD and printed matter.

**Acknowledgments:** We thank the captain, Kazumitsu Nakaguchi and the crew of the R/T vessel Toyoshio Maru, Hiroshima University for their help during the research cruise.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Concept Paper* **Facing the Forecaster's Dilemma: Reflexivity in Ocean System Forecasting**

**Nicholas R. Record 1,\* and Andrew J. Pershing <sup>2</sup>**


**Abstract:** Unlike atmospheric weather forecasting, ocean forecasting is often reflexive; for many applications, the forecast and its dissemination can change the outcome, and is in this way, a part of the system. Reflexivity has implications for several ocean forecasting applications, such as fisheries management, endangered species management, toxic and invasive species management, and community science. The field of ocean system forecasting is experiencing rapid growth, and there is an opportunity to add the reflexivity dynamic to the conventional approach taken from weather forecasting. Social science has grappled with reflexivity for decades and can offer a valuable perspective. Ocean forecasting is often iterative, thus it can also offer opportunities to advance the general understanding of reflexive prediction. In this paper, we present a basic theoretical skeleton for considering iterative reflexivity in an ocean forecasting context. It is possible to explore the reflexive dynamics because the prediction is iterative. The central problem amounts to a tension between providing a reliably accurate forecast and affecting a desired outcome via the forecast. These two objectives are not always compatible. We map a review of the literature onto relevant ecological scales that contextualize the role of reflexivity across a range of applications, from biogeochemical (e.g., hypoxia and harmful algal blooms) to endangered species management. Formulating reflexivity mathematically provides one explicit mechanism for integrating natural and social sciences. In the context of the Anthropocene ocean, reflexivity helps us understand whether forecasts are meant to mitigate and control environmental changes, or to adapt and respond within a changing system. By thinking about reflexivity as part of the foundation of ocean system forecasting, we hope to avoid some of the unintended consequences that can derail forecasting programs.

**Keywords:** ocean forecasting; reflexivity; endangered species; fisheries; harmful algal blooms; coupled natural-human systems; Anthropocene ocean

#### **1. Introduction**

The convention of studying natural systems—atmosphere, biosphere, etc.—as separate from human systems, is beginning to change [1]. Earth System Science and Anthropocene Studies increasingly couple natural and social science. As these disciplines merge, and especially in the context of applications, viewing scientists as part of a system of study is a challenge.

In forecasting applications within the natural sciences, this challenge has generally meant that a forecast is assumed to have no direct bearing on the outcome of the event in question. For example, forecasting rain does not affect whether or not it actually rains. Forecasting the time of a lunar eclipse does not alter the time of the eclipse.

In social sciences, on the other hand, it is more common for scientists to view predictions as a dynamic part of the subject system. For example, the prediction of a stock market collapse can itself cause a market collapse if investors panic in response to the prediction. This could happen even if no collapse would have occurred in the absence of the prediction [2]. This phenomenon is known as "reflexivity" and is a common component

291

**Citation:** Record, N.R.; Pershing, A.J. Facing the Forecaster's Dilemma: Reflexivity in Ocean System Forecasting. *Oceans* **2021**, *2*, 738–751. https://doi.org/10.3390/ oceans2040042

Academic Editor: Diego Macías

Received: 6 July 2021 Accepted: 4 November 2021 Published: 12 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

of human systems (Figure 1). Reflexive prediction has been examined in many fields, such as economics, political science, and even studies of faith healing [3], and has a rich history of study in the social sciences [4,5]. Reflexivity itself is sometimes seen as the property that distinguishes the natural sciences from the social sciences [6,7].

**Figure 1.** (**A**) The conventional forecasting scheme, where a system informs a forecast, which informs some human response. (**B**) A reflexive forecasting scheme where the human response is part of the system dynamics.

Natural systems forecasting has deep roots in weather forecasting, which is generally non-reflexive. However, many natural systems do have reflexive dynamics. For example, the dissemination of epidemic forecasts can alter human responses, changing the dynamics of the epidemic itself. A dire epidemic forecast could prompt a severe lockdown, thereby stifling the epidemic. Yet without the prediction, the lockdown might have come too late, and the dire outcome might have come to pass. There is evidence that the COVID-19 pandemic has reflexive dynamics and that taking these dynamics into account alters forecasts and outcomes [8].

Ocean system forecasting differs from weather forecasting in that many societally important forecasts deal with reflexive systems. Fisheries management often depends on a prediction of the stock size in future years. In turn, yearly fisheries forecasts can alter both fishing and management behavior, changing the mortality dynamics of the fish stocks. Similarly, endangered species management often relies on forecasts from population viability analysis. Management actions based on these forecasts are aimed at changing the predicted population trajectories. Even predictions of the global ocean climate system depend strongly on the human response to climate predictions themselves, where one of the explicit goals of making projections is to inform policy choices that will change the human forcing of the climate system.

Despite the prevalence of reflexivity in natural systems, it is typically not proactively taken into account in natural systems forecasting. The effects of reflexivity often lead to unexpected or unintended consequences, thus many fishery and endangered species management efforts require periodic yearly or multiannual iterative reevaluation of the objectives and targets, responding to missed targets reactively. For example, the Common Fisheries Policy of the European Commission has the requirement that any multiannual plan will "provide for its revision after an initial ex-post evaluation, in particular to take account of changes in scientific advice" [9]. Yet with increasingly more real-time forecasts available at people's fingertips, it is more urgent to understand reflexivity on a foundational level and to develop approaches to incorporate this understanding into forecasts proactively.

The feedback dynamic of reflexivity makes the phenomenon somewhat paradoxical, frustrating many attempts at prediction. A classic example is the 1948 U.S. election between Thomas Dewey and Harry Truman. Forecasts of a Dewey victory were so confident that the Chicago Tribune published the now infamous "Dewey Wins" headline on the day after the election. While it is possible that the mistaken forecast could have been due to a statistical or methodological error, the evidence suggests that the perception of a likely Dewey victory altered voter behavior [10]. That is, the forecast itself, and its dissemination, influenced (i.e., reversed) the outcome of the election. If that is true, then forecasts of a

Truman victory would have led to the opposite outcome—a Dewey win. The election was close—would it have been impossible to correctly forecast this event precisely because of the reversing effect of the forecast? This case illustrates the paradox of trying to make an accurate prediction in a reflexive system.

The paradoxical nature of reflexivity has been explored in literature, ranging from Oedipus to Ebenezer Scrooge to Asimov's Foundation, and has the same self-referential characteristic employed in philosophical cornerstones, such as Russell's Paradox and Gödel's Incompleteness Theorem [7]. The apparent paradox has led some to conclude that prediction in reflexive systems is not possible. The "Law of Forecast Feedback" [11] argues that a reliable prediction is not possible in a reflexive system. This pessimism is understandable, particularly when it comes to forecasting single binary or low-frequency events, such as elections or market collapses. Although natural sciences have often omitted reflexivity, they may offer an opportunity to address this paradox. Many natural systems forecasting programs involve high-frequency iterative forecasting, where forecasts are made and evaluated on a short time scale. The iterative nature of these applications provides an opportunity to examine how reflexivity works, and whether there are patterns that emerge or strategies that can be employed to make prediction successful despite reflexivity.

This paper examines the consequences of an iterative forecasting system having a reflexive component. It builds from a first-principles framework for prediction in ecology, adding a reflexive term to the dynamics. In particular, we incorporate two main elements of reflexive prediction: first, that the outcome would have been different without dissemination of the forecast, and second, that the forecast was believed and acted on [6]. We do not explicitly treat the mode of forecast dissemination. In practice, the mode of forecast dissemination is a key part of its influence on human behavior. For the purpose of illustrating some foundational properties of reflexivity in forecasting, we do not expand on the modes of forecast dissemination and the wide range of potential responses, but we recognize it as another important component to forecast implementation. By mapping previous ocean forecasting efforts into a biparametric time–time space, we explore how the iterative nature of many ocean forecasting endeavors can inform our understanding of reflexivity in forecasting, and we chart possible ways forward.

#### **2. Theory**

A generalized formulation of an ecosystem forecast can be written in terms of component parts as [12]:

$$Y\_{t+1} = f(Y\_t, X\_t | \theta + a) + \epsilon\_t \tag{1}$$

where *Yt* is the state variable we are trying to forecast for time *t* + 1, *Xt* are environmental covariates, *θ* and *α* represent model parameter mean and error, and *<sup>t</sup>* represents process error. To analyze the effects of reflexive prediction in an iterative forecast, we will examine a simple example of this general formulation,

$$Y\_{t+1} = \beta\_1 Y\_t + \beta\_0 + \epsilon\_t \tag{2}$$

This is the basic case for discontinuous (discrete) forecasting, where the state variable *Y* at time *t* + 1 is a linear function of the previous time step—essentially a linear autocorrelative model. The simplified formulation allows us to explore basic properties of iterative reflexive prediction. The general idea can be extended to a more complex model. To account for reflexivity, we separate the actual system trajectory, *Yt*, from the disseminated forecast, *Zt*. As pointed out earlier, we don't explore different modes of dissemination here, but we note that different modes of dissemination can lead to different forms of reflexivity.

**Case 1: No reflexivity.** In the simple model system for *Y*, the best forecast equation would be

$$Z\_{t+1} = \beta\_1 \mathbf{Y}\_t + \beta\_0 \tag{3}$$

where the disseminated forecasted *Zt*+<sup>1</sup> is identical to *Yt*+<sup>1</sup> without the error term. A forecaster could use this equation to make reliable and unbiased predictions at each time step, with uncertainty described by plus whatever uncertainty exists around the parameter measurements. This case represents the conventional, non-reflexive point of view. The forecast has high accuracy, basically limited only by the magnitude of the error terms (Figure 2A).

**Figure 2.** (**A**) Simulation with no reflexivity. (**B**) Simulation with reflexivity. (**C**) Simulation with reflexivity including a response to forecast accuracy. Error has been set to zero to make the cyclicity apparent. (**D**) Simulation with reflexivity including a response to forecast accuracy that includes memory of the accuracy over the past five time steps.

**Case 2: Self-defeating reflexivity**. In a reflexive prediction system, the outcome depends on the prediction. One way to express this is to add a reflexivity term to the general forecast equation:

$$Y\_{t+1} = f(Y\_t, X\_t | \bar{\theta} + \mathfrak{a}) + \mathfrak{e}\_t + \mathfrak{g}(Z\_{t+1}) \tag{4}$$

where *g* is some function of the disseminated forecast. This function is analogous to the "internal decision model" [13]. Here, the outcome of the event at time *t* + 1 depends on what the forecast was for that time (i.e., *t* + 1). There are two types of reflexive prediction: self-fulfilling and self-defeating (also referred to as "bandwagon" and "underdog" [4]). In a self-fulfilling reflexive system, forecasting a particular outcome makes that outcome more likely (e.g., the market collapse example). In a self-defeating reflexive system, forecasting a particular outcome makes that outcome less likely (e.g., the Truman election). Here we take the self-defeating reflexive prediction as the illustrative case.

For self-defeating forecasts, *Y* could be an index of some effect, such as the magnitude of an epidemic or the mortality rate of an endangered species—something that stakeholders would generally want to minimize. Dissemination of the forecast causes an inverse response, decreasing the magnitude of *Y*. In the linear model example, we add a response term to the forecast equation:

$$Y\_{t+1} = \beta\_1 Y\_t + \beta\_0 + \epsilon\_t - \rho\_t(Z\_{t+1}) \tag{5}$$

where *ρ<sup>t</sup>* is an increasing function of the disseminated forecast *Zt*+1. The higher the forecast value *Zt*+<sup>1</sup> is, the stronger the counteracting response will be, leading to a smaller value of *Yt*+1. The equation for *Z* represents whatever the forecaster's strategy is for predicting the system. It could be the linear model that describes the non-reflexive *Yt*, or something different.

While the scientific process behind a forecast aspires to objectivity, it exists in a broader system with subjective goals. These goals may be expressed in the kinds of processes that are forecasted (i.e., through funding choices) or in how forecasts are disseminated (e.g., the difference between forecasts of a hurricane path and communication strategies designed to get people out of harm's way). For the sake of simplicity, we will not separate the objective goals of scientific accuracy from the societal goals embedded in the application and dissemination of the forecast. Thus, when we refer to the forecaster's goals, we are essentially talking about the goals of the entire forecasting program.

Suppose the forecasting program were put in place with the ultimate aim of minimizing the value of *Y* (i.e., stop the epidemic or eliminate the mortality rate of the endangered species). Under this reflexive scenario, the naive strategy would be to always provide the direct forecast. For this example, we formulate a response term: *ρ<sup>t</sup>* = *ρ*0tanh(*Zt*+1). In other words, a high forecast value for the next time step motivates a response that counters the expectation. Using the tanh functional form caps the magnitude of the response. The forecast would have low accuracy, but the desired outcome would be achieved. On the other hand, a high-accuracy forecast would not prompt the response that minimizes the negative effect *Y*.

The consequences are twofold. First, by responding to the forecast as a warning, the actual value of *Y* is driven down. This could be considered a case where the desired effect is achieved (Figure 2B). However, a second consequence is that the forecast is now never accurate. It always overshoots the actual by an interval equal to *ρ* (on average). For this scenario to actually work, the forecast users would have to never catch on to the fact that the forecast is always more dire than reality. To put this into real terms, it would be akin to forecasting a fishery collapse every year, and although none ever occurs, the fishery repeatedly reduces catches as though a collapse were always imminent. This contradiction between forecast accuracy and forecast utility (from the perspective of the desired societal outcome) is the central point to the Law of Forecast Feedback.

**Case 3: Iterative self-defeating reflexivity**. Realistically, people would lose trust in a consistently dire forecast due to its consistent lack of accuracy. This is where the iterative dynamic comes into play. To account for this, we introduce a scaling factor *τ* to the response term that represents the reliability of the forecast:

$$\mathcal{Y}\_{t+1} = \beta\_1 \mathcal{Y}\_t + \beta\_0 + \varepsilon\_t - \rho\_t \pi\_t \tag{6}$$

Here *τ* is an inverse function of the error in the forecast (i.e., <sup>|</sup>*Zt*−*Yt*<sup>|</sup> *Yt* ), where *τ* = 1 when accuracy is perfect and goes to zero for very low-accuracy forecasts. In other words, as the forecast becomes inaccurate, it also loses its influence. Now we have the simplest fully iteratively reflexive forecast model.

We can formulate this factor *τ* as *τ<sup>t</sup>* = *e* −*τ*0|*Zt*−*Yt*| *Yt* . For a perfect forecast, *τ<sup>t</sup>* = 1, yielding a full response to the forecast. For very high inaccuracy, *τ<sup>t</sup>* decays to zero, zeroing out the response term. The parameter *τ*<sup>0</sup> shapes how quickly (as a function of forecast inaccuracy) the response term goes to zero. A high *τ*<sup>0</sup> would mean that only a small amount of inaccuracy is needed for people to stop believing in and responding to the forecast. The

result is an oscillating pattern, where a reliable forecast is acted on, driving *Y* down, thus making the next forecast inaccurate, diminishing the response, and driving *Y* back up (Figure 2C). This is akin to the boom–bust reflexive dynamics seen in market systems [7].

**Case 4: Iterative + learning self-defeating reflexivity**. As a final note, there's no reason to assume that the response only depends on the previous time step. Depending on circumstances, it is possible that collective memory would evaluate the forecast reliability over multiple previous time steps. This can be added to the model using a number of time steps *m*, over which *τ* is computed and averaged. The result is a variably reliable forecast, with periodic lapses in accuracy (Figure 2D). From here, it is not difficult to imagine a wide range of periodic and quasi-periodic patterns that can occur depending on the form of *τ<sup>t</sup>* and other properties of these equations. All of the richness of dynamical systems modeling could appear in the formulation of reflexivity.

#### **3. The Forecaster's Dilemma**

The question for the forecaster now becomes: how to deal with these opposing forces? On the one hand, a theoretically reliable forecast can alter behavior, making the forecast unreliable. On the other hand, consistently unreliable forecasts are likely to be ignored. The problem for the forecaster can be framed as the tension between two goals:

**Goal 1: The accuracy directive.** Conventionally, forecasters have tried to make predictions that accurately describe a future event. This also corresponds with goals of science to improve our understanding of the natural world. When the event comes to pass, a comparison between the forecast and the event serves as the assessment. This amounts to minimizing ∑*<sup>t</sup>* |*Zt*−*Yt*| *Yt* .

**Goal 2: The influence directive.** The purpose of a forecast is usually to elicit some action. This generally corresponds with some practical societal goal. The *Y* variable represents a negative effect that the forecast is aspiring to diminish over time, so this amounts to minimizing ∑*<sup>t</sup> Yt* (This could also be framed as maximizing a positive effect, such as species recovery).

A forecaster in a reflexive system should consider whether it is possible to meet these two goals simultaneously, and if so, what is the best forecasting strategy i.e., the choice of function for *Z* that accomplishes both directives? The example provided here is convergent in a recursive sense. That is, one can iteratively plug *Yt*+<sup>1</sup> back into the equation as *Zt*+1, and the forecast for the next time step will converge on a value that is both accurate and minimizes the negative effect, basically toeing a line between the two cases. However, most real-world examples will probably be more complicated, with more dynamic and complex *g*(*Z*) functions.

#### **4. Solving the Forecaster's Dilemma**

Reflexivity is not just of academic interest. The coronavirus pandemic brought home the point that reflexivity in forecasts can have very real consequences. As people come to use and expect increasingly more real-time forecasting, the issue of reflexivity represents an emerging scientific challenge. With the field of ocean system forecasting growing rapidly and experiencing a "Cambrian Explosion" in applications [14], forecasters should consider the role of reflexivity when adopting a forecasting program. We suggest three components to guide consideration of reflexivity. In some cases, dealing with reflexivity will be a nonissue. When reflexivity is deeply embedded within a system, however, new scientific and pragmatic approaches will need to be developed.

#### *4.1. Step 1: Identify How Reflexivity Could Occur in the System*

Ocean systems are organized across a range of temporal, spatial, and biological scales. These scales are related such that to first order, organisms that are small tend to have lower trophic level and shorter generation times and vice versa due to the strong size structuring in ocean ecosystems [15]. Ocean system forecasts tend to align along similar temporal, spatial, and biological scales because of this natural hierarchal structuring (Figure 3).

Most ocean system forecasting programs reviewed here fall along the one-to-one line in Figure 3. This scale matching is probably not an accidental coincidence, as decision making and predictability often align with the dominant process scales in a system. At one end of the spectrum, long-term planning for long-lived protected or exploited species is informed by their comparatively longer reproductive schedules, which aligns with a longer forecasting range [16,17]. Similarly, plankton with rapid doubling times fall at the other end of the spectrum, with shorter forecast ranges [18]. However, there are some examples that deviate from this one-to-one line. In the upper left region are harmful planktonic taxa like *Vibrio* and *Alexandrium*, which are short-generation taxa for which longer forecasts have been made [19,20]. This can be motivated either by forecasting and monitoring limitations, or by a particular forecast use. In the lower right, there are examples of short-range forecasts for protected or endangered taxa, such as whales and sea turtles. These few short-range forecasts for higher trophic levels focus on forecasting location rather than abundance, making the generation length less relevant. This also represents a response to a particular forecast use.

**Figure 3.** Forecast ranges of ocean system forecasting research programs and their associated biological time scales sources from literature examples, extended from the list provided in [14]. Positions are estimated based on descriptions in the texts [16–18,20–53].

This schematic needs to be shifted slightly to understand the role of reflexivity. Most of these forecasting studies do not consider the time scale of human response. Considering the coupled natural human system angle offers a slightly different lens, where we consider the dominant response time scales of the entire system, including the ecosystem and the human system together (Figure 4). If the response time is much longer than the forecast range (lower right of Figure 4), for example, reflexivity will be minimal to non-existent. In this scenario, many iterative forecasts would be made before any response takes place, so the human response would not affect the forecast. There are also cases where the human response has no bearing on the forecast. For example, a jellyfish forecast [46] might guide recreational activities, but probably would not influence the jellyfish populations themselves. Similarly, a forecast of the abundance of a harmful algal species might lead

to a fishery closure, but the bloom would persist unaffected. In any of these scenarios, if forecasting is guiding monitoring efforts, there is the possibility of feedback even in non-reflexive systems that can lead to confirmation bias.

If the coupled system response time is similar to or shorter than the forecast range, then there is the potential for reflexivity. These cases would fall near the one-to-one line or upper-left triangle of Figure 4. Examples include short-range forecasts of harmful or toxic algal species, or mid- to long-range forecasts for protected or exploited species. If there is the potential for reflexivity, the next step is to try to measure and describe it.

**Figure 4.** Schematic of forecast and system process scales. Here, the system includes the coupled natural human processes.

## *4.2. Step 2: Determine Whether Reflexivity Is Self-Defeating, Self-Fulfilling, or a Combination of Both*

Once the potential for reflexivity has been characterized, the next step represents the empirical or data analysis side of the question. Detecting or measuring the human response can be challenging, especially if monitoring and modeling efforts have focused on the non-human parts of the system.

Consider, as an example, the management of endangered marine mammals. The Marine Mammal Protection act in the U.S. was amended in 1994 such that marine mammals are managed according to potential biological removal (PBR), which frames population viability as a forecast: i.e., "a certain probability of an event occurring in a given amount of time" [54,55]. The calculation of PBR includes a scaling factor *FR*, which is set between 0 and 1. If a forecast looks particularly dire, setting *FR* to zero sets human-caused mortality to zero (in theory) regardless of the value of other factors. *FR* offers flexibility—and reflexivity—to the human response via the PBR calculation [56], and these numbers inform management decisions designed to reduce (or not) human-caused mortality of marine mammals. This places these forecasts toward the upper left of Figure 4, with long-range forecasts but short system response times.

If there is a way to quantify human response to the forecast, such as PBR, one can then examine whether the dynamics of this response is coupled to the dynamics of the natural system, such as in a dynamical systems framework. Two-way causality can be explored using tools such as Granger causality or convergent cross mapping [57], though there are cases where these tests are flawed [58]. In principle, it may be possible to detect reflexivity in ecological time series without measurement of the human response using Takens' theorem [59].

The species with the longest running PBR time series is the North Atlantic right whale, for which there are 25 years of data [60]. For the first decade of the 2000s, PBR was reduced

to zero by setting *FR* to zero in an attempt to restore this highly endangered species. The population climbed until 2010, at which point PBR was raised. Shortly thereafter, mortality rates began to rise (Figure 5A), and the species plummeted into crisis mode again [61]. When comparing PBR to human-caused mortality for this species, the strongest lagged correlation is with PBR leading mortality by four years (Figure 5B). Granger analysis supports this direction of lagged causality (*p* < 0.05, but see [58]), consistent with the notion that changing the forecast for the species has the following effect of changing the actual population trajectory of that species.

In some sense, this direction of causality shouldn't be surprising: the goal of the management approach is to have a following effect on the population. But the right whale example is a case of self-defeating reflexivity, and self-defeating reflexivity can cut two ways. A dire forecast can motivate recovery efforts, thereby improving the forecasted outlook, as intended—but a more favorable forecast can have the opposite effect. In the right whale time series, increases in PBR were followed a few years later by increases in human-caused mortality, reversing the recovery trend that the right whale population had been showing. While there are likely multiple factors at play in the sudden reversal of the right whale population trajectory, including climate and oceanographic changes [61,62], the pattern is consistent with the dynamics we would expect in a self-defeating reflexive forecasting system, leading to an unintended consequence. Ideally, we would like reflexive feedback to take effect when the population is struggling, but not when it's doing well—in other words, a concave-down curve in Figure 5B rather than concave up. Accounting for this kind of reflexive dynamic more deliberately probably requires a more mechanistic understanding of the reflexive term *g*(*Z*).

**Figure 5.** (**A**) Time series of PBR for the North Atlantic right whale (red) and two estimates of mortality: documented human-caused mortalities (grey) and the annual population change from the Pace model, subtracting out new calves (blue). (**B**) Lagged relationship between potential biological removal (PBR) and mortality. Data aggregated from NOAA and the North Atlantic Right Whale Consortium reports [60].

## *4.3. Step 3: Incorporate Human Response into a Forecast Model*

If there is significant reflexivity in a forecasting system, with important consequences, the next step is to try to formulate that response mathematically and incorporate it into a model. This step represents an open area of scientific research and theory. A key question to answer here is: Can the accuracy and influence directives both be met?

Ocean forecasting up to now has largely followed the tradition of weather forecasting, combining mechanistic or processed based formulations, such as the advection-diffusionreaction equation, with statistical formulations. More recently, machine learning algorithms have been replacing the statistical components, and to some extent the mechanistic components as well, to derive predictive rules from data. This framework has so far been mostly an elaboration of the *f*(*Z*) term. The reflexive term *g*(*Z*) represents a largely unexplored

research opportunity where similar approaches could be used. In fisheries forecasting, stochastic models have been used to couple these components [63].

Coupled natural-human systems are complex. There probably is not an analog for the Navier-Stokes equations for the human part of the system. However, examining mathematical formulations in a theoretical context can help answer whether the accuracy and influence directives are at odds with each other. If they are, then it could be a sign that the forecast will do more harm than good—or at least not have the intended effect.

There is also the opportunity to think beyond dynamical systems and statistics. For example, because information is exchanged in at least two directions in reflexive systems, some have viewed forecasting in the context of game theory [64]. This view becomes particularly interesting when multiple agents, including forecasters and users, are exchanging wide ranging information about the system in question. On the ocean, forecast users can be capturing, processing, and exchanging real-time observations and knowledge, such as where a commercial fish species is found—knowledge that is not available to the forecaster but plays a role in human response [65]. Similarly, multiple forecasters might be using different knowledge and approaches, and exchanging some part of that information with each other. In this way, forecasting programs can be nested within networks of social-natural systems with complex information flow.

The challenge of reflexive systems forecasting highlights the need to be making more operational forecasts. The iterative property of these forecasting systems is key, providing the data and experience needed to build up an understanding of reflexive dynamics. The tendency is often to focus forecasting efforts on high-stakes problems, such as endangered species or health hazards, but lower stakes problems (e.g., nuisance species, ecotourism) could provide a safer arena for building up the datasets needed to analyze and understand reflexive dynamics in ocean systems in new ways.

## **5. Conclusions: Reflexivity in the Changing Ocean**

Reflexivity highlights the significant human dimension and associated challenges in emerging forecasting programs. Traditionally, whether forecasting the weather or some ecosystem process, the natural system is viewed in an objective sense, separate from the human observer. In the context of ocean systems, reflexivity is an emerging challenge that has bearing both to how we understand and interact with the ocean, and how we understand and make use of algorithms.

Regarding human interactions with the ocean, the "Anthropocene Ocean" is described as a socio-material space [66] where physical and biological systems are interlinked with social and scientific systems. In this context, there are two common frameworks useful for understanding ocean system forecasting. One framework is that of planetary security mitigating the risks of environmental damage due to human activities. In this context, forecasting would serve as an aid to monitoring and controlling environmental processes. Reflexivity is implicit in this construct, as the human response is the mechanism for influencing the environment. Here both the accuracy and influence directives are important for the forecast to be effective.

With the urgency around issues like climate change, there are practical limitations to the "measure and control" approach to dealing with the Anthropocene. An alternative emerging perspective is the idea of correlational and relational epistemologies, where management structures would sense and adapt to events in real-time [66,67], without a causal understanding or attempt to mitigate the dynamics. This perspective also relies on algorithmic and digital technologies, but in this context, forecasting serves as information connectivity and not as a means of system control. Reflexivity is not necessarily implicit in this perspective. When reflexivity is present, a forecast could potentially still be useful without both of the two directives met, if it is nested within a larger network of correlational human responses and feedbacks [67]. This view has drawbacks as well, and potential for unintended consequences. As each new forecasting system comes online, we should ask whether it is intended to understand (in a causal sense) and control aspects of the ocean

system, or whether it is intended as part of a network of adaptation tools. This distinction will help narrow the scope of reflexivity in the forecasted system.

The question of whether the goal is to predict and control the environment, or to adapt and respond to a changing environment, is at the core of many discussions about big data, algorithms, and artificial intelligence. As forecasting algorithms become more widespread and embedded in our social relationship with Earth systems, ocean science can take lessons from the growing field of algorithmic accountability. Across applications ranging from resume sorting to prison sentencing, algorithms are replacing human decision making. The proliferation of algorithms in this way has led to many unintended consequences [68]. Ocean system forecasting shares this risk of unintended consequences—something that has already occurred in a few ocean forecasting programs [69,70]. For reflexive forecasts, when the accuracy and influence directives are at odds with each other, there is high potential for unintended consequences. The field of algorithmic accountability is developing methodologies for addressing this, such as action plans for redress when unintended outcomes occur, which can be applied to ocean forecasting to help prevent unintended consequences or address them when they occur [71–73].

Despite the potentially confounding nature of reflexivity, the subject represents a rich area of scientific inquiry. The reflexive term in the forecasting equation—i.e., the *g*(*Z*)—captures an emerging challenge in natural systems forecasting. Many forecasting evaluation analyses choose not to treat the reflexive feedback dynamic [74], and others have simply ignored it [75]. Some leave the human response to the realm of policy, communications, or to forecast users, while others view this part of the equation as a focus for quantitative study and analysis in its own right [2]. The example developed here separates *f*(*Y*) and *g*(*Z*) into additive terms, but it is possible that they interact in more complex and nonlinear ways yet to be discovered. There is also another layer of complexity added to *g*(*Z*) dynamics when considering the mode of forecast dissemination. People respond differently depending on how a forecast is communicated, and if the system is reflexive, then communication choices can feed back on the natural system dynamics *f*(*Z*). By representing iterative system forecasting as a combination of two components, *f*(*Y*) and *g*(*Z*), we see a promising quantitative starting point for integrating natural sciences with social and behavioral sciences, as well as a pathway for using forecasts as a tool for navigating the complex interactions between humans and the ocean.

**Author Contributions:** Conceptualization, N.R.R. and A.J.P.; methodology, N.R.R.; software, N.R.R.; validation, N.R.R.; formal analysis, N.R.R.; investigation, N.R.R. and A.J.P.; resources, N.R.R. and A.J.P.; data curation, N.R.R.; writing—original draft preparation, N.R.R. and A.J.P.; writing—review and editing, N.R.R. and A.J.P.; visualization, N.R.R. and A.J.P.; supervision, N.R.R. and A.J.P.; project administration, N.R.R. and A.J.P.; funding acquisition, N.R.R. and A.J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** Support for this research came from institutional funds from the Tandy Center for Ocean Forecasting at Bigelow Laboratory for Ocean Sciences and the Otto Mønsteds Fond via Denmark Technical University (NRR) and from NSF grant OCE-1851866 (AJP).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Code and data used in this study are stored in the following open repository: https://github.com/SeascapeScience/forecasters-dilemma (accessed on 1 November 2021).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Review* **Glow on Sharks: State of the Art on Bioluminescence Research**

**Laurent Duchatelet 1,\*, Julien M. Claes 1, Jérôme Delroisse 2, Patrick Flammang <sup>2</sup> and Jérôme Mallefet 1,\***


**Abstract:** This review presents a synthesis of shark bioluminescence knowledge. Up to date, bioluminescent sharks are found only in Squaliformes, and specifically in Etmopteridae, Dalatiidae and Somniosidae families. The state-of-the-art knowledge about the evolution, ecological functions, histological structure, the associated squamation and physiological control of the photogenic organs of these elusive deep-sea sharks is presented. Special focus is given to their unique and singular hormonal luminescence control mechanism. In this context, the implication of the photophore-associated extraocular photoreception—which complements the visual adaptations of bioluminescent sharks to perceive residual downwelling light and luminescence in dim light environment—in the hormonally based luminescence control is depicted in detail. Similarities and differences between shark families are highlighted and support the hypothesis of an evolutionary unique ancestral appearance of luminescence in elasmobranchs. Finally, potential areas for future research on shark luminescence are presented.

**Keywords:** shark; luminescence; Etmopteridae; Dalatiidae; Somniosidae; photophore; hormonal control; counter-illumination

## **1. Introduction**

Bioluminescence is the ability of living organisms to produce visible light [1]. Mention of this phenomenon dates back to antiquity with the description of "cold light" by Aristotle in his book "De Anima". Many centuries later, Charles Darwin, on board the Beagle, witnessed and described light in water as "milky sea" in his logbook. The first studies demonstrating mechanisms underlying bioluminescence appeared in 1667 with Robert Boyle, who depicted the oxygen requirement for luminescence production.

Bioluminescence is the result of a spontaneous exergonic chemical reaction involving the oxidation of a luciferin catalyzed by a luciferase [2], which produces a transitory excited state that finally relaxes by emitting a photon with oxyluciferin as final product [3–5]. Luminous systems involve either luciferase and luciferin as separate components or a complex molecule called "photoprotein" comprising a preoxidized luciferin and a luciferasic activity [5]. Luciferins are often common to different taxa while luciferases are thought to be typically species-specific [6] but exceptions to this rule have been recently highlighted [7–9]. Similar luciferases, first described in species from same phyla (within a clade), are now found between phylogenetically distant species supporting the existence of multiphyletic distribution of these luciferases [7–9].

Rare in terms of species number, luminescence nevertheless arose in over 700 genera from 13 phyla covering all kingdoms, except plants and archaea, e.g., [6,10,11]. This phylogenetic diversity results from numerous independent evolutions of light-producing capability—in fact, it is estimated that bioluminescence arose independently more than 90 times during evolution [1,3,9,10,12,13], which suggests that luminescence is of paramount

**Citation:** Duchatelet, L.; Claes, J.M.; Delroisse, J.; Flammang, P.; Mallefet, J. Glow on Sharks: State of the Art on Bioluminescence Research. *Oceans* **2021**, *2*, 822–842. https://doi.org/ 10.3390/oceans2040047

Academic Editors: Antonio Bode and Eric Clua

Received: 1 September 2021 Accepted: 3 December 2021 Published: 17 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

ecological importance, but also that the acquisition of light emission capability is a relatively easy and quick process [1].

Eighty percent of luminous taxa inhabit marine environments, from coastal shallow waters to abyssal depths, and are mainly found in bacteria, protists, ctenophores, cnidarians, annelids, mollusks, chaetognaths, crustaceans, echinoderms, tunicates, and fishes [1,6,14,15]. Within the water column, the mesopelagic zone (i.e., 200–1000 m depth) has been estimated to host a high proportion of luminous organisms, e.g., 70% of mesopelagic bony fishes appear to be luminous [16]. Comparatively, terrestrial and freshwater luminous animals are only represented by earthworms, snails and limpets, fireflies, beetles, and some insect larvae, e.g., [17–23].

Although observations of shark luminescence have been reported for almost two centuries [24], the first research projects dedicated to shark luminescence were only initiated in 2005 and focused on the study of the ecological functions and physiological control of the photophores of a single species, the velvet belly lanternshark, *Etmopterus spinax* [25]. Since then, shark luminescence research has flourished, with detailed phylogenetical, ecological, histological, and physiological studies now available for numerous species, e.g., [26,27]. By synthesizing the findings of those studies, this article aims to provide not only a holistic view of current shark luminescence knowledge but also perspectives for future research.

#### **2. Luminous Shark Diversity**

Among cartilaginous fishes, only sharks have evolved the ability to emit light. Indeed, no report of bioluminescence exists for ratfishes (Chimaeriformes) apart from the mention of a luminous fluid on the rabbit fish, *Chimaera monstrosa*, in 1810 interpreted as coming from luminous bacteria due to the degradation of the harvested organism [28]. In addition, the anecdotal observation of putative photogenic organs in the deep-sea dark blind ray, *Benthobatis moresbyi* [29], has not been confirmed by more recent and detailed morphological studies [30,31]. Bioluminescence in sharks appears currently restricted to Squaliformes, where only three families (Dalatiidae, Somniosidae, and Etmopteridae) contain luminescent representatives (Figure 1). Indeed, although bioluminescence has once been suggested for the specific supralabial white band of the megamouth shark, *Megachasma pelagios*, recent work invalidated this assumption [32]. Although the luminescence reported for the species that belong to Etmopteridae and Dalatidae was initially thought to have been acquired independently [33,34], the recent discovery of luminescence from a somniosid, the velvet dogfish, *Zameus squamulosus* [35,36], now strongly suggests that the acquisition of luminescence capability in sharks represents a single evolutionary event, which occurred during the deep-sea radiation of Squaliformes at the end of the Cretaceous [35]. Although fossil studies estimate the Etmopteridae radiation around 90 million years ago [37], molecular data presents a separation of Etmopteridae from other Squaliformes during the Upper Cretaceous (i.e., 65–90 million years ago) [35,38]. The Dalatiidae family radiated later, during the Paleocene after the Cretaceous/Paleocene mass extinction, 65–105 million years ago, when they substituted extinct marine reptiles and fishes in the epipelagic fauna before reaching the deep sea [34,35,39].

Interestingly, the photogenic structures appear ubiquitous in Etmopteridae (four genera: *Trigonognathus*, *Aculeola*, *Centroscyllium* and *Etmopterus*; 52 species) and Dalatiidae (seven genera: *Dalatias*, *Isistius*, *Mollisquama*, *Euprotomicroides*, *Squaliolus*, *Euprotomicrus* and *Heteroscymnoides*; 10 species). Nevertheless, *Z. squamulosus* is the only somniosid shark known to possess photophores (Figure 1; Table S1). In parallel, luminescence has been observed in live in 15 species only, however, covering most clades, i.e., the blurred smooth lanternshark, *Etmopterus bigelowi* [40]; the southern lanternshark, *Etmopterus granulosus* [41]; the blackbelly lanternshark, *Etmopterus lucifer* [41]; the slendertail lanternshark, *Etmopterus molleri* [42,43]; the smooth lanternshark, *Etmopterus pusillus* [44]; *E. spinax* [45,46]; the splendid lanternshark, *Etmopterus splendidus* [47]; the green lanternshark, *Etmopterus virens* [40]; the smalleye pygmy shark, *Squaliolus aliae* [43,48]; the kitefin shark, *Dalatias licha* [41]; the taillight shark, *Euprotomicroides zantedeschia* [49]; the pygmy shark, *Euprotomicrus bispina-* *tus* [50]; the cookiecutter shark, *Isistius brasiliensis* [51,52]; *Z. squamulosus* [36]; the viper dogfish, *Trigonognathus kabeyai* (Mallefet, unpublished data); Figures 1 and 2; Table S1. In addition, expected luminous species are encountered in Etmopteridae and Dalatiidae, based on the presence of photophore structures and/or flank marks in the holotype description [33,35,36,38,40–49,51–90]; Table S1. This led to a total number of 62 luminous sharks, i.e., ~11% of currently described species (550 in total [91]), while in comparison, luminescence is thought to have appeared in only ~5% of bony fishes [12].

**Figure 1.** Shark luminescence distribution within Squaliformes families based on published phylogenies [35,38,53,92–96]. Circles inside the tree represent the luminous (blue), expected luminous (gray), and non-luminous (white) status of each represented species. Statuses are based on in vivo pictures, physiological studies (luminous), the presence of photophores or flank marks (expected luminous), and none of these criteria (non-luminous) (see Table S1). Circles outside the tree, scaled to the total number of species (number in brackets) within a given family [54,91], indicate the proportion of luminous (blue), expected luminous (gray), as well as non-luminous (white) species. For each family, total number of luminous/expected luminous/non-luminous species are given next to the outside circle. Blue star indicates the expected origin of luminescence in sharks.

**Figure 2.** Shark luminescent patterns. Squaliform cladogram with taxonomic grouping where photophores/flank markings (gray shade) or live luminescence (blue shade) have been observed [35,36,40,41,49–74,97,98] and spontaneous luminescence of representative shark species lateral view for *D. licha*, *S. aliae*, *E. bispinatus*, excreted bioluminescent fluid for *E. zantedeschia* (blue bubbles on the cladogram), and ventral view for the others]. D, Dalatiidae; E, Etmopteridae; FM, flank markings; S, Somniosidae. Indicative scale bars, 5 cm. Shark drawings are modified from [34,74]. Photographs courtesy by D. Perrine (*E. bispinatus*), T. Raczynski (*E. zantedeschia* [49]) and J. Mallefet (other species; [36,41,65], Mallefet, unpublished).

#### **3. Ecology of Shark Luminescence**

Deciphering the ecological functions of luminescence from elusive animals such as deep-sea sharks is extremely challenging. Indeed, field observations are scarce and unbiased lab experiments proved to be difficult to perform. As a matter of fact, most functions formulated over years regarding the function of shark luminescence remain undemonstrated, mainly due to the difficulties of observing and collecting these rare and elusive organisms to conduct ethological studies on the function of their luminescence. Fortunately, however, detailed analyses of photophore distribution (luminescent "patterns") as well as physical characteristics (intensity, wavelength and angular distribution) and kinetics of luminescence coupled to physical models for pelagic vision and molecular phylogenetic analyses, now allow us to draw a clearer picture of the adaptive benefits of luminescence in sharks.

In this context, counterillumination, i.e., a camouflage technique used by midwater organisms cloaking their silhouette from upward-looking organisms using a ventral glow mimicking downwelling sunlight, e.g., [99–103], is probably the primary function of shark luminescence, for both defensive and predatory purposes [41,65]. Indeed, shark photophores are predominantly situated on the ventral surface area (Figures 2 and 3) and produce a light whose color (wavelength) which is similar to that found in coastal (blue–green) and oceanic (blue) environment (Figure 3b; Table 1). Moreover, the intensity and angular distribution of light emission have been correlated to capture depth residual light parameters in three species: *E. spinax* [103], *E. splendidus* [65], and *S. aliae* [65].

**Figure 3.** Shark luminescence functions. (**a**) Composite pictures of *E. molleri* and *E. spinax* luminescence patterns. Spontaneous dorsal, lateral and ventral luminescence from *E. molleri*, with key photophores areas highlighted. Circular inserts represent spine-associated photophores (SAPs) of the dorsal fins of *E. spinax* [104]; these photophores are absent in *E. molleri* (which instead has spine base-associated photophores; [105]). (**b**) Angular distribution of dalatiid (*S. aliae*; mid body) and etmopterid (*E. spinax*; mid body and lateral flank marks) luminescence as well as downwelling sunlight; the perfect match observed for *E. spinax* is achieved thanks to a centripetal change in photophore orientation modified from [65,103]. (**c**) Simplified bioluminescent shark phylogeny providing an overview of demonstrated, experimentally supported and putative functions of shark luminescence [27,41,74,75,104–108], Claes and Duchatelet, unpublished observation. AP, associated photophores; FP, frontal photophores; LZ, luminous zone; OP, ocular photophores; SAP1, spine-associated photophores of 1st dorsal fin; SAP2, spine-associated photophores of 2nd dorsal fin; SBAP, spine base-associated photophores. Scale bar, 5 cm. Photographs by J. Mallefet.


**Table 1.** Shark luminescence color.

Interestingly, "water box" experiments suggested that *E. spinax*, contrary to other counter-illuminating animals, was not able to adapt rapidly the intensity of its luminescence by more than one order of magnitude in response to changing light regime [103]. Hence an "isolume follower" hypothesis was formulated, stating that sharks move up and down in the water column to remain cryptic [103]. This has been further supported by a large comparative study showing that daytime capture depth of luminescent sharks could be predicted from their ventral photophore cover (e.g., the percentage of the ventral surface area covered with photophores), which varies from 2% in the bareskin dogfish *Centroscyllium kamoharai* to 56% in the viper dogfish *T. kabeyai* [65]. In addition, recent studies show that some etmopterid species have a higher swimming speed and muscular enzymatic activities than their non-luminous counterparts living in the same deep environment [109,110]. These findings could be correlated with a greater physiological demand to perform vertical migration to remain camouflaged. The second most widespread and well-supported function of shark luminescence is intraspecific communication. Although Dalatiidae and Somniosidae present "simple" luminescent patterns where photophores follow a dorsoventral density gradient [33,34,36,41,52] (Figures 2 and 3), etmopterid sharks display complex luminescent photophore aggregations on the ventral area, but also on the flanks, the fins, the tail, around the eyes, the spiracles, the gills, and the epidermal tissue surrounding dorsal spines [34,104–106,111] (Figures 2 and 3). Since these patterns are species-specific and do not show sexual dimorphism (except photophores associated with primary sexual characters, e.g., male claspers, as explained below [112]), they are often used as taxonomic determination tools [34,54,63,68]. Interestingly the shape of lateral luminous areas (flank markings) appears to be clade-specific [38,65] (Figure 2) and their presence correlates with an increased speciation rate [75] and a moderate predation risk [65]. This, and the ability for those sharks to discern the flank marking shape (and other specific luminous zones) at a biologically meaningful distance (1–3 m, or three to ten body lengths, as indicated by visual modeling [111]) strongly support the idea that etmopterid flank markings represent an exaptation of counter-illuminating photophores to facilitate deep-sea communication [65]. Luminescence could also be used as a mating aid, which allows males to identify females from a distance (since photophores cover male claspers of all species [34]; Figure 3a) and visualize their cloaca and pectoral fins (brighter in etmopterid [34]; Figures 2 and 3a) in the darkness of the deep sea (sharks display internal fertilization during which the male stabilizes itself by biting the female's pectoral fin [97]). Sexual dimorphism in the luminescence time course of pharmacologically stimulated pelvic photophores from *E. spinax* further supports this hypothesis [107].

The last luminescence function for which experimental support is available is the aposematism, a mechanism by which an animal advertises potential predators that it is not worth attacking or eating [113,114]. Etmopterid sharks (contrary to dalatiid and somniosid species) have large sharp defensive spines associated with their dorsal fins [97]. In some species, photophores either placed on the edge of the dorsal fins (*E. spinax*; [104]) or around the base of the spines (*E. molleri*; [105]), allow these spines to be seen in the dark from a distance by potential predators, and hence potentially work as an aposematic signal as strongly suggested by experimental data [104,105] (Figure 3a–c).

In addition to these relatively well-established functions, researchers have formulated more speculative functional hypotheses regarding the luminescence of some species (Figure 3c). For instance, the luminescent liquid ejected by the pelvic pouch of *E. zantedeschia* probably works as a defensive "smokescreen" mechanism [49]. Given the morphological similarity of the pectoral glands, an identical function has been suggested for the pocket sharks, *Mollisquama* spp. [55,74]. The ocular and frontal photophores of *E. spinax* [103] and the ventral photophore of *D. licha* [41] could be used as a vision aid (Figure 3c). In addition, the tail of *D. licha*, *S. aliae*, *Z. squamulosus*, *T. kabeyai,* and *Etmopterus* species, which is more mobile and brighter than the rest of the ventral surface area, could work as a distracting lure, and hence be analogous to the caudal photophores of myctophid and tubeshoulder fishes [115–117]. Conversely, the idea that the dark, photophore-free collar of *I. brasiliensis* acts as a lure to attract bigger pelagic fishes or marine mammals on which it feeds [51], now seems at the very least dubious given that numerous common preys of this shark species are either filter feeders or top predators for which such a mechanism is useless [65], and that the closely related species, the largetooth cookiecutter shark, *I. plutodus*, which has a similar diet [97], lacks such a collar [70].

Etmopterid and dalatiid sharks display a set of visual features not found in nonbioluminescent sharks, which strongly supports the idea that their visual system coevolved with their ability to produce light. Aphakic gaps in *S. aliae* [111] and *E. bispinatus* Claes, personal observation] and translucent upper eyelid found in etmopterid sharks (genus *Trigonognathus* and *Etmopterus*; [111]) probably play a role in counterilluminating, e.g., facilitating perception of downwelling light intensity to ensure the perfect match, similarly to what was recently shown for deep-sea bony fishes [118]. On the other hand, the spectral absorbance (485–488 nm) of etmopterid rod photoreceptor as well as the retinal distribution of ganglion cells of those shark (which appear species-specific) appear finely tuned to detect their own bioluminescence, especially from flank markings [111].

#### **4. Photogenic Structures and Specialized Squamation of Bioluminescent Sharks**

Sharks display two types of photogenic structures: photophores (Figure 4), with internal luminescence present in all bioluminescent shark species, and secretory glands ejecting a bioluminescent fluid in the environment (external or secretory luminescence; Figure 4) present in *E. zantedeschia* (pelvic pouch) and in *Mollisquama* spp. (pectoral pockets).

Shark photophores are small, reaching a maximum diameter of ~50, ~100, and ~200 μm for Somniosidae, Dalatiidae, and Etmopteridae, respectively [34,36,41,48,52,65,76]. Comparatively, the photophores of bony fishes can reach up to 1 cm [119]. Given their high photophore density (sometimes over ten thousand per square centimeter; [65]), sharks are probably the luminous organisms with the highest number of photophores (a 52 cm TL male *E. spinax* has been estimated to bear about 440,000 photophores; [34]). The general structure of photophores is similar across families (Figure 4a,c,d).

**Figure 4.** Shark photogenic organs. (**a**) Photophores (modified from [32,34]); (**b**) secretory photogenic structures: lightproducing columnar epithelium from pouch of the taillight shark (modified from [34]; left) and putative light-producing cubic tissue from pocket shark's pectoral pockets (modified from [74]; right); (**c**) histological section of *I. brasiliensis* photophore; (**d**) histological section of *E. lucifer* photophore; (**e**) ultrastructure of *E. spinax* photophore; (**f**) ultrastructure of *E. spinax* photocyte; (**g**) autofluorescence from *E. lucifer* photophore; (**h**) autofluorescence from photocyte vesicle from *E. lucifer*; (**i**) autofluorescence from pocket shark's pocket epithelium. A, apical cells; B, basal cells; Ba, basal lamina; Bs, blood sinus; Ca, cartilage rod from pectoral fin. C1, cellular type 1; C2, cellular type 2; E, epidermis; G, granular area; I, inclusion; ILS, iris-like structure; L, lens cell; N, nucleus; P, photocyte; Ps, secretory photocyte; R, reticulated layer (reflector); S, pigmented sheath; V, photocyte vesicles; Va, vesicular area. Indicative scale bars in (**a**–**e**,**g**,**i**), 50 μm; in (**f**,**h**), 10 μm.

Etmopteridae photophores are composed of 6–14 emitting cells (photocytes) embedded in a cup-shaped pigmented cell sheath, covered by a reflective layer containing guanine crystals and capped by one or several lens cells [16,41,44,76,120]. A multilayer cell zone, called iris-like structure (ILS), is present between the lens cells and the photocytes and is used as a photophore shutter [43,76,77,120–122]. Renwart et al., 2014, described three different cell types constituting the ILS of *E. spinax* photophores: (i) the cellular type I, which contains fibrous material and was assumed to stabilize the lens cells; (ii) the cellular type II, with nucleus only visible via electron microscopy and whose function is unclear, and (iii) the pigmented cells, affiliated to melanophores, presenting pseudopodialike cellular projections [76] (Figure 4). Photocytes of *E. spinax* were ultrastructurally depicted as containing three distinct areas: the nucleus, the vesicular and the granular areas (Figure 4e,f). The granular area is assumed to be the location of the bioluminescence chemical reaction and associated microsources were named "glowons" [76,120]. Green autofluorescence is observed within photocytes (especially from their granular area) under blue/UV light exposure (Figure 4g,h), which is assumed to be due to the fluorescence properties of the bioluminescent substrate (as it is the case, e.g., for luciferin; [45,47,76,77]). Shark photophores appear to lack intra-organ innervation [46,76], with terminal epidermal nerves reaching only the surrounding of the photophore as demonstrated through acetylated-tubulin labelling [123].

312

Comparatively, Dalatiidae and Somniosidae harbor less complex and smaller light organs containing a single photocyte, a pigmented sheath and a group of small lens cells [33,34,36,41,52] (Figure 4a). Pigmented cells, attributed to ILS cells, surmounting the photocytes also act as a photophore shutter [41,77]. Photophore development has been studied in *E. spinax* and *S. aliae*, revealing a similar four-phase morphogenesis pattern: (i) apparition of pigmented cells; (ii) formation of the pigmented sheath; (iii) apparition of the protophotocyte (i.e., photocytes not able to produce light); (iv) maturation of the photocyte during which the photocyte acquires its luminescence competence (revealed by the presence of fluorescent vesicles in photocytes; [45,77], Duchatelet, personal observation).

The secretory epithelium from the pelvic pouch of *E. zantedeschia* displays a pseudostratified columnar epithelium with three distinct cell types: columnar cells with a large apical inclusion (which are probably secretory cells), flattened cells present between columnar cells and extending from the basal lamina towards the free surface of the epithelium, and superficial cells forming a thin cover above the apical part of columnar cells [61]. The epithelium of the putatively bioluminescent fluid-producing pectoral pockets of *Mollisquama mississippiensis* shows a surprisingly different structure, since it appears to be stratified, made of over 50 layers of cuboidal cells displaying green autofluorescence and probably involved in the production of a holocrine secretion [74] (Figure 4b,i).

Given that shark photophores are embedded in the shark epidermis, they compete for space with dermal denticles (placoid scales). As a consequence, bioluminescent sharks evolved specific squamation patterns, i.e., pavement-like, cross-, bristle-, hook-, and simple/alveolar leaf-shaped patterns [32,34–36,41,124,125]; Figure 5, which allows the accommodation of photophores between (pavement-like, cross-, bristle-, hook-shaped types) or below (leaf-shaped types) the placoid scales. Histology and light transmission analyses of leaf-shaped type denticles from *Z. squamulosus* recently revealed specific honeycomb structures allowing the transmission of at least 50% of the light produced [36]. On the specific rostral area, the kitefin shark, *D. licha*, also presents highly translucent leaf-shaped placoid scales but without honeycomb structure [41].

**Figure 5.** *Cont*.

**Figure 5.** Bioluminescent shark squamation. Placoid scale patterns found in bioluminescent sharks [32,34,36,41,124]. Photophore are represented with blue (open) and dark (closed) dots between placoid scales. Indicative scale bars, 50 μm.

#### **5. Control of Bioluminescence from Shark Photophores**

## *5.1. Hormonal Control*

Luminous sharks occupy a peculiar place among bioluminescent organisms regarding the control of their luminescence. Indeed, their photogenic organs are primarily controlled by hormones rather than by nerves, the general condition for intrinsic photogenic organs [42,43,47,48,108]. Although the majority of bony fishes have a nervous control of the light emitted (e.g., via adrenaline and noradrenaline, [126–129]), such classical neurotransmitters have been demonstrated to be inefficient to regulate the shark light emission [16,108]. Intraperitoneal injection of drugs such as adrenalin or acetylcholine in *I. brasiliensis* failed to trigger light emission [16]. Claes and Mallefet (2009) also attempted, without success, to induce light response in the etmopterid *E. spinax* after the application of neural agents (e.g., adrenalin, noradrenalin, serotonin, carbachol, and the classical depolarizing agent KCl) on isolated ventral photogenic skin patches, a technique which became the standard for pharmacological studies of organism luminescence control [108]. The absence of photophore innervations, and the inefficiency of classical neurotransmitters to trigger light emission, led to assumptions of a hormonal control of luminescence in sharks [28,108].

Using *E. spinax* as a model species, and hypothesizing that the shark light emission involves pigment motion within the ILS-melanophores, Claes and Mallefet (2009) have highlighted the implication of melatonin (MT), prolactin (PRL) and alpha melanocytestimulating hormone (α-MSH) in the light emission process [108] (Figure 6a). These three hormones are known to be involved in the physiological skin color modulation in elasmobranchs (i.e., sharks and rays): PRL, and α-MSH stimulate melanophore pigment dispersion and, thus, induce skin darkening, and MT regulate melanophore pigment aggregation, leading to a paler skin [130,131]. A few years later, the hormonal control was described in other species such as *E. splendidus* [47], *E. molleri* [42], *E. lucifer* [41] and *E. granulosus* [41], as well as in two Dalatiidae species, *S. aliae* [48] and *D. licha* [41]. Recently, the adrenocorticotropic hormone (ACTH), a melanocortin hormone also involved in vertebrate skin pigment motion within melanophores, was found to inhibit light emission induced by MT in both Etmopteridae and Dalatiidae [41,43] (Figure 6a). Conversely, the melanin-concentrating hormone (MCH), another hormone responsible of the pigment

granule aggregation in bony fishes [132–135] had no effect on shark luminescence [136]. Although the pharmacological control of *Z. squamulosus* photophores has not been investigated yet, these organs can be assumed to be under hormonal control as well, given their structural similarity with dalatiid photophores [35,36]. Future work on this subject will allow us to confirm this hypothesis.

**Figure 6.** Shark photophore control. (**a**) Overview of shark luminescence control agents mapped against experimentally investigated species over the last 15 years [25,27,41–43,47,48,77,107,108,121,122,137–139]. Letters in circles indicate the identified targets of the agent (ILS, iris-like structure; P, photocytes); (**b**) photophore-containing skin patch before hormonal stimulation (left) and at luminescence peak (right) from *E. spinax* [121]; (**c**) primary hormonal control pathway of shark photophore regulating ILS shutter function (modified from [122]).

Overall, the photophore responses to hormone application appear similar across investigated species, except for PRL (Figure 6a). In Etmopteridae, MT and PRL, both at a concentration of 10−<sup>6</sup> mol L−1, trigger light emission, while α-MSH at 10−<sup>6</sup> mol L−<sup>1</sup> and ACTH at 10−<sup>5</sup> mol L−<sup>1</sup> inhibit it [42,43,47,108] (Figure 6a). Both stimulating hormones present different light-emission curves when applied on photogenic ventral skin patches. Although MT triggers a slow increasing and long-lasting luminescence for up to one hour depending on the species, PRL application induces a fast–high response, decreasing rapidly after a few minutes. Nevertheless, it has been shown that both hormones act on the ILS

cells pigment movement to open the "iris" and allow the outward light emission [43,121] (Figure 6b). Conversely, α-MSH and ACTH inhibit the light emission with a decrease in the amount of light produced after both MT or PRL applications, and also act on the ILS cells pigment movement by occulting the photocyte, and hence preventing the light from being emitted outside [41,43,108,121,122]. In Dalatiidae, the hormonal control is similar, except for PRL, which was demonstrated to inhibit the light emission process in *S. aliae* [41,48]. The physiological effects of each hormone for each studied shark species are summarized in Figure 6a.

In silico transcriptome analyses and immune localization allowed us to confirm the presence and location of MT and α-MSH/ACTH receptors within the photophore from two Etmopteridae species, *E. spinax* and *E. molleri* [137] (Figure 6a). These G protein-coupled receptors (GPCR), known to be coupled with specific G proteins, are mainly localized within the cell membrane of photocytes and ILS cells [137]. Researchers also highlighted that PRL receptor mRNA sequences are absent from *E. spinax* ventral skin reference transcriptome [137], in agreement with previous works revealing the secondary loss of PRL receptor within the elasmobranch lineages during the growth hormone/prolactin protein family evolution [140]. Therefore, how PRL could cellularly have been perceived by the photophore to trigger luminescence remains enigmatic and assumptions have been made on the involvement of the closely related growth hormone and its specific receptor in the light emission control.

An integrative model of the hormonal control of shark photophore luminescence, highlighting the role of MT, α-MSH and ACTH on both the photocyte and the ILS, has been recently proposed [122] (Figure 6c). MT, through its perception by MT receptor, triggers the release of Gi proteins, inhibiting the adenylate cyclase activity [43,122]. Adenylate cyclase is known to directly act on cellular cAMP concentration [141–143]. Confirming results from Claes and Mallefet, 2009 [108], cAMP concentration assay after MT-induced light emission revealed a drastic decrease in cAMP levels [43,122]. In fish melanophore, a decrease in the cell cAMP concentration triggers the aggregation of melanin pigment toward the nucleus periphery [142,144,145]. Moreover, MT application was demonstrated to activate, through a calmodulin/calcineurin pathway, the cellular motor dynein, which carries pigmented granules towards melanophore-like nucleus periphery [122]. This MT pathway regulates the "opening" of the ILS cells and allows the emitted light to reach the outside of the photophore [122] (Figure 6). The ultrastructure analyses of the ILS structure before and after the MT-induced luminescence confirmed the pigment motion within the ILS-associated melanophores [120]. On the other hand, analysis of the α-MSH/ACTH pathway revealed how these hormones "close" the photophore and, hence, inhibit light emission. Through melanocortin hormone (α-MSH/ACTH) perception, its receptors release a Gs protein, activating the adenylate cyclase, and an increase in the intra-ILS cells cAMP concentrations [43,122]. The final step of this inhibiting pathway involves the cellular motor kinesin, which carries pigment granules on cytoskeleton towards the ILS melanophore pseudopodial-like projection, occulting the light produced [122] (Figure 6c). These two ILS-regulating pathways seem to be conserved across Etmopteridae and Dalatiidae [41,43,47,122] (Figure 6a).

In addition to the hormonal control, pharmacological studies also revealed that γaminobutyric acid (GABA) and nitric oxide (NO) can modulate the light emission in certain species [42,48,123,139] (Figure 6a). In the model species, *E. spinax*, GABA inhibits and NO modulates the MT/PRL-induced luminescence [123,139]. Conversely, NO application on *S. aliae* photogenic skin patches, after MT/PRL light emission triggering, had no effects [48] (Figure 6a).

Further research is needed to determine the intracellular interplay between the neuromodulators (NO and GABA) and the hormones as well as to better understand the intracellular events occurring in the photocytes upon hormonal stimulation.

## *5.2. Extraocular Photoreception and Pigment Motion Regulation*

Extraocular photoreception, the ability to detect and respond to light clues outside of the "eyes", has been suggested to be involved in the bioluminescence control in various invertebrate taxa. In the bobtail squid, *Euprymna scolopes*, extraocular photoreceptors and photocytes are colocalized within photophores. Expression of eye-like developing genes in the light "pocket" of *E. scolopes* underlines such a coupling between photoemission and photoreception mechanisms [146,147]. In the comb jelly, *Memniopsis leidyi*, and the burrowing ophiuroid, *Amphiura filiformis*, both photoreceptor and photocyte molecular markers (i.e., opsin and photoprotein/luciferase, respectively) are coexpressed within the same cells. These observations support the existence of a functional link between light perception and bioluminescence control for these species [8,148–150]. By analyzing the opsin and phototransduction genes expression within photophores, light organ photosensitivity was also suggested for a counterilluminating deep-sea shrimp, *Janicella spinicauda*, where it was assumed to play a role in the light emission regulation to ensure a match with the residual downwelling light [151].

As indicated by recent analyses of ventral skin and eye transcriptomes of *E. spinax,* photoreceptors and phototransduction actors were expressed in both tissues [152]. This species displays only two ocular opsins, one rhodopsin and one putatively associated peropsin, highlighting its monochromatic vision [152]. A specific extraocular photopigment, the encephalopsin (Opn3), was also detected in the skin transcriptome of *E. spinax* [152]. The colocation of this extraocular opsin (Es-Opn3) with the ventral skin photophores provides fuel for the putative functional coupling between light emission and light perception in luminous organisms [152]. More precisely, membranes of the ILS cells were shown to be the main expression site of the Es-Opn3, while no expression was found at the photocyte level [152]. Besides, the expression of this opsin was demonstrated to be concomitant with the *in utero* development of the photophores in *E. spinax* and *S. aliae* embryos, appearing with the transformation of protophotocytes to photocytes (i.e., when photocytes can produce luminescence; [77], Duchatelet, unpublished data), which supports the idea that the Es-Opn3 is used to detect embryonic luminescence [77]. The evaluation of the absorbance spectrum of this luminous shark extraocular opsin has added new evidence of the link between the two photobiological processes, since a clear overlap exists between the light emission spectrum of *E. spinax* and the Es-Opn3 photopigment light absorbance [122].

Going further in the phototransduction cascade, Duchatelet et al., 2020, demonstrated an impact of blue light exposure on the intracellular concentration of inositol triphosphate (IP3) from photogenic ventral skin [122]. This IP3 concentration modulation confirmed opsin activity and highlighted the first step of the phototransduction cascade occurring at the photophore level. Pharmacological analyses unveiled the next steps of the phototransduction pathway, with the implication of calcium, the Ca2+-dependent calmodulin, and the cytoplasmic motor dynein [122], clearly demonstrating the interconnection between these pathway steps and the pathway leading to the pigment granule aggregation in melanophores. Therefore, Duchatelet et al., 2020, proposed a model of light emission control in *E. spinax* based on the photoperception of its own luminescence that regulates ILS melanophore pigment aggregation and, thus, the aperture of the photophore, which regulates the light output [122] (Figure 6c). Interestingly, while hormones operate at both the photocyte, triggering light emission through an unknown biochemical reaction (i.e., unknown luciferin/luciferase system or photoprotein) and ILS (to regulate the amount of light transmitted as a camera diaphragm) levels, Es-Opn3 appears to be involved only in the pigment motion regulation of the ILS. This dual control is assumed to occur at least in the luminous species belonging to Dalatiidae and Etmopteridae (Duchatelet, unpublished data). The close link between photoemission and photoreception at the light-organ level, putatively regulating the amount of light emitted, seems to have emerged independently in phylogenetically distant luminous organisms: ctenophores, cephalopods, crustaceans, echinoderms, and luminous sharks [8,122,146–152]. It is a fascinating example of convergent evolution.

#### **6. Biochemistry of Shark Luminescence**

One of the remaining questions is what type of bioluminescent system is involved in the emission of light in sharks. Two categories of light systems are known to date: (i) luciferase–luciferin systems, and (ii) photoprotein systems. In the former, a substrate, the luciferin, is oxidized by an enzyme, the luciferase, in the presence of oxygen and often other cofactors [5]. In the latter, an enzyme complex as well as a preoxidized luciferin form together a complex protein, the so-called photoprotein, requiring only the contribution of a cofactor (often an ion) to produce light [5,153]. Various luciferin types have been molecularly described in the marine environments to date, such as the coelenterazine (the most phylogenetically widespread luciferin type in the oceanic environment) and its derivatives (dehydrocoelenterazine and enol-sulfate coelenterazine), the dinoflagellate luciferin, or the *Cypridina* luciferin (also called vargulin) [1,13,154]. Associated with these luciferins, different types of luciferases (nine different types) and photoprotein (three different types) have been molecularly depicted [1,5,9,13]. These substrates and enzymes, although initially considered to be species-specific [28], can sometimes be shared by phylogenetically very distant species [5–9]. This has led scientists to hypothesize that certain species can acquire the necessary components for the luminous reaction through their diet [1,5,155–159]. Recently, the golden sweeper fish, *Parapriacanthus ransonnetti*, has been shown to acquire not only its luciferin, but also its luciferase by feeding on luminous ostracods [160]. It is noted that the only species described as being able to synthesize de novo their luciferins are the midwater shrimp *Systellaspis debilis* [161], the ostracod *Metridia pacifica* [162], and the ctenophores *Bolinopsis infundibulum* and *M. leidyi* [163].

Attempts were made to identify the bioluminescent compound responsible for the light emission in *E. spinax* by analyzing the cross-reactivity of known luciferins with extract of shark skin assumed to contain the catalyst (i.e., luciferase). In parallel, the putative diet acquisition of a luciferin through the food chain was hypothesized. Only coelenterazine, the most commonly found luciferin in marine taxa [5,164–167], was found in the digestive tract of *E. spinax*, but none of the tested luciferins (coelenterazine, dinoflagellate/euphausiid luciferin, and *Cypridina* luciferin) has reacted with the shark photogenic skin extract [168]. Similar negative results have been obtained for the cross-reactivity with the respective luciferase: coelenterazine-dependent luciferase (*Renilla* luciferase), cypridinid luciferase, or euphausid luciferase [168]. Analysis of the transcriptomic sequence available for the photogenic skin failed to pinpoint any putative homologs of known luciferase/photoprotein sequences [52,152], leading to assume the involvement of an unknown light-emitting system in luminous sharks, including either an unknown light-emitting system (photoprotein) or a specific active or storage form of a known luciferin using a shark-specific luciferase. Research combining data from transcriptomics, proteomics, and bioinformatics modeling are ongoing, and could allow us to decipher this enigmatic bioluminescent system.

#### **7. Conclusions and Perspectives**

Bioluminescent sharks have fascinated humans for almost two centuries. Yet, dedicated research on these enigmatic deep-sea inhabitants involving spectrophotometry, luminometry, pharmacology, light/electron microscopy, biochemistry, molecular analyses, and transcriptomics started only 15 years ago. Results from those studies (over 50 publications in total) have been synthesized in the current review. Overall, findings suggest that luminescence acquisition has been a unique though pivotal evolutionary event for squaliform sharks, which greatly facilitated their radiation in deep-sea habitats and strongly shaped their visual system. From a primary function of camouflage (coopted from the hormonally controlled crypsis mechanism of shallow water elasmobranchs) in Dalatiidae, Somniosidae and basal Etmopteridae, shark bioluminescent patterns progressively became an intra- and interspecific communication tool in derived etmopterid sharks (e.g., *Etmopterus* species), an exaptation that considerably increased their speciation rate and has probably been facilitated by the increase in size and complexity of etmopterid photophores (which allows better orientation of light output, e.g., tangential to the body surface). The

secretory luminescence observed in *E. zantedeschia* (and putatively present in pocket sharks, *Mollisquama* sp.), which might have evolved from the invagination and modification of epidermal photophores, represents an example of the high selective pressure occurring in the darkness of the deep sea.

It clearly appears that the future of shark bioluminescence research will also be driven by new molecular data and techniques. The Next-Generation Sequencing methods recently allowed the emergence of transcriptomic studies in non-model organisms such as some selected luminous shark species (i.e., *E. spinax*, *I. brasiliensis* [52,152]). These studies paved the way for future transcriptomic, proteomic, and genomic studies on luminous sharks. Among an infinite number of fascinating questions, these studies could focus on the identification of the light-emitting molecular toolkit (luciferase, photoprotein, etc.) in luminous sharks.

Up to now, genomic resources for cartilaginous fish are still very limited (i.e., genomes available for the great white shark, *Carcharodon carcharias*, the whale shark, *Rhincodon typus*, the elephant shark, *Callorhinchus milii*, and the little skate, *Leucoraja erinacea*) and absent for luminous sharks. Genome size estimation in different Etmopterid species revealed that these species have among the largest genomes of all investigated Chondrichthyes (e.g., genome size of *E. spinax* might reach 16 Gbp, against 4.63 Gbp in *C. cacharias*, and 3 Gbp in *R. typus*) [169–173]. These impressive genome sizes possibly reveal lineage-specific genome expansion and large-scale alteration events such as gene/genome duplications, transposon insertions or events such as polyploidy [174], and confirm the challenging future of genomic studies on these species.

Further research on these sharks is planned to understand their ecology. In vivo ethological studies could lead to a better understanding of their behavior, notably the "isolumes follower" hypothesis by passive acoustic tagging of these sharks. Another approach could be to follow "in situ" reaction to stimuli mimicking the light emission of conspecifics.

Although recent research allowed to draw a clearer picture on the evolution, ecology and physiology of shark luminescence, gaps in our knowledge of these fascinating animals still exist. In particular, the interplay between hormonal control pathway and neuromodulation as well as the chemistry of shark luminescence remain to be determined. Work is underway to clarify these shadow areas.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/oceans2040047/s1, Table S1: Overview table of the luminous species within the luminous Squaliformes (Etmopteridae, Dalatiidae, Somniosidae) according to the literature.

**Author Contributions:** Conceptualization, L.D., J.M.C. and J.M.; software, L.D.; investigation, L.D., J.M.C., J.D. and J.M.; writing—original draft preparation, L.D., J.M.C. and J.D.; writing—review and editing, L.D., J.M.C., J.D., P.F. and J.M.; supervision, J.M.; project administration, L.D.; funding acquisition, P.F. and J.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Fonds National de la Recherche Scientifique (F.R.S.-FNRS, Belgium), grant number T.0169.20 awarded to the University of Louvain—UCLouvain Marine Biology laboratory and the University of Mons—UMONS Biology of Marine Organisms and Biomimetics laboratory. L.D., J.M.C., J.D., P.F. and J.M. are, respectively a postdoctoral researcher, scientific collaborator of the UCLouvain laboratory, postdoctoral researcher F.R.S.-FNRS, research director F.R.S.-FNRS, and research associate F.R.S.-FNRS.

**Acknowledgments:** This review is the contribution BRC # 381 of the Biodiversity Research Center (UCLouvain) from the Earth and Life Institute Biodiversity (ELIB) and the "Centre Interuniversitaire de Biologie Marine" (CIBIM).

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Oceans* Editorial Office E-mail: oceans@mdpi.com www.mdpi.com/journal/oceans

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com