Biodiversity research on the Romanian shelf grew vigorously after WWII from the seminal work of Professor M. Bacescu and collaborators. Between the 1950s and 1970s, numerous species were described, including micro- and macrophytobenthos, Nemertini, Nematoda, Polychaeta, Mollusca, and Crustacea (Bacescu et al., 1971). The taxonomic and ecological contributions of these researchers remained underrecognized in a period of major environmental changes in the Black Sea (Gomoiu, 1992; Mee, 1992; Zaitsev, 1993; Bologa et al., 1995; Leppäkoski and Mihnea, 1996): growing coastal euthrophication (Mee and Topping, 1999; Cociaşu et al., 2008) and associated events of phytoplankton and macroalgal blooms (Bologa et al., 1995; Moncheva et al., 1995; Bodeanu, 2002), resulted in population declines or extinctions (Gomoiu, 1983, 1992), the introduction of new species (Gomoiu et al., 2002), and the modification of sediment regimes due to hydrographic disruption in tributary rivers (Panin et al., 1996; Panin and Jipa, 1998; Duţu et al., 2018). The combined impact of these stressors led to the restructuring of benthic assemblages with 60% of benthic fauna to be lost in three decades and local declines in population numbers and diversity (Petranu, 1997; Gomoiu, 1998).

The lag in marine benthic research in a period of such striking changes has been linked to the small number of specialists and limitations in resources due to the broader economic decline in the region (Gomoiu et al., 2003/2004). After 1990, the socio-economic recession in the riparian Black Sea countries and the introduction of national and international programs aiming at reducing river inputs and eutrophication, gave the Romanian coastline an opportunity to recover. However, despite the slow recovery observed in various marine ecosystems in the region (Mee et al., 2005; Friedrich et al., 2014; Mureşan et al., 2020; Teacă et al., 2020), it remained fragile and vulnerable to anthropogenic impacts, ecological accidents, and climate change.

Research on marine ecosystems, in general, and benthic communities, in particular, intensified in the 2000s with the application of the European Directives and the growth of research funding. Some of that research was directed toward macrozoobenthos in MPAs (Abaza and Godeanu, 2008; Abaza, 2010; Niţă et al., 2012), mapping of benthic habitats according to the Natura 2000 classification (Zaharia et al., 2012; Ungureanu et al., 2015), the study of benthic communities associated with rocky bottom substrates and macrophytic algae (Teacă et al., 2006a,b, 2010; Filimon et al., 2016), and the assessment of circalittoral habitats (Begun et al., 2010, 2018; Abaza et al., 2018a,2019; Mureşan and Teacă, 2019; Mureşan et al., 2019; Teacă et al., 2019, 2020). Today, the Romanian MPA network encompasses 7,457 km² of the Romanian shelf distributed among ten MPAs. Two MPAs in particular, the Cape Tuzla (ROSCI0273) and the Mangalia submerged Sulfur springs (ROSCI0094), are examined in this study and together represent 1.8% (107.3 km²) of the surface of all Romanian protected areas. Such area was reached only recently, when the area of both MPAs increased by 65 and 93%, respectively, following the approval of the MMAP Order no. 46/2016.

The growth of these protected areas also entailed a challenge now faced by many other MPAs: the knowledge on large portions of protected seafloor remains extremely limited, despite the need to apply European Union’s policies, such as the Marine Strategy Framework Directive, the Birds and Habitats Directives (Nature Directives hereafter) and the most recent European Green Deal (EC, 2019) and Biodiversity Strategy for 2030 (EC, 2020). As in other MPAs, the description, mapping and status of macrozoobenthic communities in the Cape Tuzla and the Mangalia Sulfur springs MPAs are critical and are examined in this study. We aimed to assess biodiversity, map benthic habitats according to the EUNIS classification, and to achieve a synoptic overview of the state of these communities. We attempted this by embracing the modern approach to management of marine habitats of European policies, such as the 2008/56/EC Directive, which requires European Member States to achieve a Good Environmental Status (GES) by 2020.

The Black Sea is an inland sea with a maximum depth of 2,245 m (Sorokin, 2002), connected in the northeast with the Sea of Azov and in the southwest with the Sea of Marmara and further south with the Aegean and Mediterranean Seas (Figure 1). The Romanian part of the Black Sea coastline spans 244 km and can be divided in a northern part, between Musura Bay and the Cape Midia Harbor jetties, and a southern part between the Cape Midia and the Romanian-Bulgarian border. The southern part includes the most diverse geomorphological features of the coastline (Panin and Jipa, 1998; Panin, 2005), and is characterized by the presence of almost continuous loess cliffs (12–20 m high). This area is also more stable in terms of geodynamic processes than the northern deltaic coast (Panin and Jipa, 1998), but is more exposed to anthropogenic stressors. The shoreline south of Constanta does receive a scarce input of sediments from the Danube system, so the sediments in this area are mainly of organic origin, derived from mollusk shells, and affected by continuous erosion caused by inappropriate hydrotechnical protective works.

The two MPAs assessed in this study, the Cape Tuzla and the Mangalia Sulfur springs are located in the Romanian southern part extending down to 50 m depth (Figure 1). They occupy 107.3 km² and comprise loess cliffs which are constituted by discontinuous Quaternary red clays and Sarmatian limestone. Their erosion has created a heterogeneous subtidal bottom with extensive blocks of rock and boulder fields going down to 25 m in depth (Panin, 2005). These MPAs were officially designated under the EC Habitats Directive in 2009 (Decision 2009/92/EEC) to protect these unique and rare marine ecosystems, the submerged sulfur springs, extensive canopies of Cystoseira and Zostera meadows, and their associated natural resources from multiple socio-economic impacts. The MPAs support a rich biodiversity and have been recognized as biodiversity hotspots and important nursery grounds for several protected fish species through Habitat Directives and have been also listed in the Annex II of the Directive 92/43/EEC. They have been also included in the national list of the Natura 2000 sites, hosting four distinct habitats (Table 1).

During the seven field campaigns conducted between 2012 and 2018, 75 stations were sampled in total—47 stations in the Cape Tuzla and 28 stations in the Mangalia Sulfur springs MPAs (Figure 2). Other 20 samples were also collected for visual inspection and habitat identification on board of the research vessel R/V “Mare Nigrum.” The W-SE bathymetry gradient of the study areas ranged between 11.1 and 50.4 m deep. Folk 5/7 was used for sediment classification of the marine habitat delineation according to the EMODnet Seabed Habitats substrate and National classification schemes (Populus et al., 2017; Vasquez et al., 2020). Habitat distribution maps were then developed using the combination of sediment type, depth and point (station) data. Such habitat distribution was subsequently confirmed (in 2021) using remote sensing techniques such as multibeam and Side-Scan Sonar (Popa, Pers. Comm.). A single sample was collected at each station using van Veen grabs with grasping areas of 0.135 m² and 0.125 m². Samples were sieved through a 0.5 mm mesh (Bacescu et al., 1971; Todorova and Konsulova, 2005) and the retained organisms were preserved in 4% neutralized formalin-seawater solution. Density and biomass (wet weight) were estimated as ind. m^––2 and g m^––2. Species nomenclature was checked using the World Register of Marine Species¹ and all the specimens collected were deposited at the National Research-Development Institute for Marine Geology and Geoecology—GeoEcoMar.

Macrozoobenthic community structure was analyzed using species composition (S), their density (A), dominance (D), frequency (F), diversity (H’), and biomass (B). Diversity was calculated using the Shannon-Wiener index (H’) on a log₂ base (Shannon and Weaver, 1963). Spatial analyses of that data were conducted with Ocean Data View, version 5.4.0 (Schlitzer, 2018). Univariate and multivariate analyses were performed using the PRIMER 7 and PERMANOVA + software package, version 7.0.17. Multidimensional Scaling (MDS) and cluster analyses, in particular, were conducted with density and biomass, but the biomass offered a better visualization of the dissimilarity between stations and is the one presented here. The AZTI Marine Biotic Index, AMBI (Borja et al., 2000), and the multivariate AMBI, also called M-AMBI (Muxika et al., 2007), were calculated using the software (version 5.0),² and the assignment of the identified species to the five Ecological Groups (EG) was based on the species list (June 2017) in the AMBI index software (Borja and Muxika, 2005). Subsequently, M-AMBI*(n) was calculated using the arithmetic mean of the minimum and maximum normalized AMBI, H’ and S (Sigovini et al., 2013). The use of AMBI and related indexes follows the recommendations given by National Authority Reports (Abaza et al., 2018b) and the intercalibration initiative between Romania and Bulgaria (Todorova et al., 2018). For the infralittoral and circalittoral sedimentary habitats, a threshold value of 0.68 for M-AMBI*(n) was used to discriminate between “Good Environmental Status” (GES) and non-GES (Abaza et al., 2018b).

A total of 119 taxa were identified from both MPAs along the 7-year study period (Supplementary Material 1). Diversity was dominated by polychaetes (28.6%), followed by mollusks (25.2%), crustaceans (25.2%), and to a lesser extent by nemerteans (5%), cnidarians (5%), and “others,” encompassing a combined 11% of all other taxa identified. 95 taxa were found in the Cape Tuzla MPA, dominated by polychaetes, molluscs, and crustaceans (26, 25, and 24%, respectively). Fairly similar results were gathered in the Mangalia Sulfur springs MPA, where 98 taxa belonged to twelve groups and were dominated by the same taxa (28.6, 26.5, and 23.5%, respectively). 72 taxa (60.5%) were found in both MPAs, while 23 and 26 taxa were collected exclusively in the Cape Tuzla and the Mangalia Sulfur springs MPAs, respectively (Supplementary Material 1). Out of the species collected, five are currently listed in the Red Book of the Invertebrates from Romania (Murariu and Maican, 2021): Upogebia pusilla, Carinina heterosoma, Liocarcinus navigator, Gastrosaccus sanctus, and Brachynotus sexdentatus.

Our point-sample study covered an overall area of 92% of the Cape Tuzla MPA and 90% of the Mangalia Sulfur springs MPA. Seven large benthic habitats were delineated (see below) based on the broad approach of manually drawing on a sampling point grid, following community structure in terms of presence of ecosystem engineers, their coverage, and the type of substrate. Cluster and multidimensional scaling plot (MDS) using biomass data (see section “Materials and Methods”) showed that the dissimilarity between benthic assemblages in different stations ranged from 83 to 96% (Figure 3), whereas similarity values ranged between 34 and 56% (Figure 4) (ANOSIM global R = 0.552, p = 0.1%, permutations = 999). Besides the clear spatial ordination based on bootstrap biomass, the results showed resemblance between habitats with similar type of substrate (as in the case of Infralittoral and Circalittoral coarse sediments). Moreover, according to the SIMPER analysis (see Supplementary Table 2), the species which contributed the most to the habitat differentiation were ecosystem engineers, including Mytilus galloprovincialis, Chamelea gallina, and Upogebia pusilla.

Following the depth gradient, three infralittoral and four circalittoral habitats were distinguished (Table 2 and Figure 5). For both MPAs, the total number of species in the infralittoral and circalittoral added up to 62 and 103, respectively. Benthic communities were dominated by polychaetes (23 and 28 species in infra- and circalittoral, respectively) and crustaceans (18 and 24 species, respectively). Among polychaetes, 17 species were common to habitats present in both biozones, whereas 85% of the abundance and 75% of the biomass were accounted by the detritivorous species Heteromastus filiformis and Prionospio maciolekae and the predatory species Nephtys hombergii and Micronephthys longicornis (Supplementary Material 1). Mollusks reach the highest diversity in circalittoral habitats (27 in comparison to 11 species in the infralittoral), due mainly to soft-bottom engineering species such as Acanthocardia paucicostata, Papilicardium papillosum, Spisula subtruncata, Abra alba, Abra nitida, Pitar rudis, and Modiolula phaseolina.

The infralittoral rocks with photophilic algae and/or associated Mytilus galloprovincialis reefs were treated as a major geological unit and ranged between 1 and 20 m deep (based on Bacescu et al., 1971; Teacă et al., 2010; Teacă et al., pers. obs.), with a larger breadth in the Cape Tuzla MPA than in the Mangalia Sulfur springs MPA. The seabed in the former MPA was a rugged landscape marked by boulders of different sizes, whereas the seabed in the latter was rather flat in its upper part and more irregular toward deeper seafloor. Nearby and as a continuation of these reefs, the following benthic habitats were described (Table 2):

The AMBI values for the 75 stations ranged between 0.72 and 5.14, indicating that macrobenthic communities from 25 stations (33.3%) achieved GES and that from 50 stations (66.6%) did not achieved GES. M-AMBI*(n) values ranged between 0.32 and 1.00. Based on this index, GES was assigned to a total of 50 stations and non-GES to 25 stations (Figure 8). M-AMBI*(n) values were generally lower in offshore (circalittoral) than in inshore (infralittoral) bottoms. The circalittoral mud habitat achieved a 54% GES, while other circalittoral habitats (sandy mud, shelly coarse sediment, or mixed sediment) reached between 70 and 100% of GES. Infralittoral shelly coarse sand and shellbed with various infauna and infralittoral fine and medium sand dominated by Chamelea gallina habitats achieved 88 and 60% of GES, respectively. In general, relatively high local diversity (20–35 species in average per habitat) translated into high values of H’ and M-AMBI*(n) index. Only the infralittoral muddy sand habitat accounted for a 100% non-GES due to the small number of samples available in this habitat (Figure 8).

The proportion of ecological groups (EG) to which benthic species were allocated, showed the dominance of species “sensitive” and “indifferent” to organic pollution in infralittoral habitats (Chamelea gallina, Mytilaster lineatus, Parvicardium exiguum, Apherusa bispinosa, Melita palmata, Microdeutopus versiculatus, Iphinoe tenella, Diogenes pugilator, Upogebia pusilla, and Liocarcinus vernalis). In these habitats, 7–57% of the species were considered sensitive and 13–75% of the species were considered “indifferent.” In contrast, first-order “opportunistic” species (Heteromastus filiformis, Oligochaeta) were better represented in the circalittoral mud habitat (15–69%).

Mapping provides basic information about the distribution, pattern, extent and area of habitats of conservation importance, and through that supports spatial planning and ecosystem-based management of MPAs (Todorova et al., 2015). Our study encompassed habitats associated with the first two (out of five) biological benthic zones that have been identified in the Black Sea: infralittoral, circalittoral, offshore circalittoral, bathyal, and abyssal. The two first zones, infralittoral, and circalittoral bottoms, represent almost 14% of the entire Black Sea seafloor, and ∼2.3% of the Romanian coastline. The soft-bottoms associated with infralittoral and circalittoral habitats are affected by wave action (2.9 N/m²), whilst the subsequent zones, offshore circalittoral, bathyal and abyssal bottoms, only start at the sea-bottom 9.7°C isotherm (at about 50 m depth) (Populus et al., 2017; Vasquez et al., 2020).

The seven habitats described in this study span only 107.3 km² of the infralittoral and circalittoral zones, a relatively small surface when compared with the extent of lifeless areas like the anoxic bathyal and abyssal (Populus et al., 2017). Nevertheless, the ecological roles of these habitats are deemed essential for the overall health of the Black Sea. Important benthic species and communities associated with these habitats provide goods (such as food provision), cultural services (e.g., recreational diving in biogenic reefs) and have key regulatory functions including sediment bio-irrigation, reworking and mobilization, in addition to serving as carbon and nutrient sinks (see Lenihan and Peterson, 1998; Aller and Cochran, 2019).

Hard bottoms constituted a continuous geological unit stretching along both MPAs, and which probed to be associated with several types of benthic habitats. The main habitat linked to the rocky facies is the Infralittoral rock with photophilic algae and/or Mytilus galloprovincialis reefs, ranging between 1 and 20 m deep (Bacescu et al., 1971; Teacă et al., 2010). In the Mangalia Sulfur springs MPA, the rocky facies reach slightly deeper bottoms in the circalittoral, causing the partial overlap with the Infralittoral and Circalittoral coarse habitats, whose precise limits remain to be defined at a later stage using remote sensing. It is also noteworthy to mention that some species that inhabit shallow coarse sediments in the infralittoral and circalittoral bottoms (e.g., H. filiformis and U. pussila), exhibit a large biotope plasticity, with no apparent preferences for infralittoral sand or circalittoral muddy sands. However, such plasticity warrants further research attempting to understand the dynamics of species colonizing or migrating among distinct habitats.

In spite of their small surface (0.75–2.14 km²), the infralittoral as well as the circalittoral mud, coarse and mixed habitats, sustained a significant pool of benthic diversity. The infralittoral habitats (down to 20 m depth) had been already the subject of prior mapping research for the ‘‘Management Plan of the Tuzla and Mangalia Site of Community Importance,’’ as defined according to the Natura 2000 classification.³ However, no data was available before our study on macrozoobenthic diversity, their distribution, and benthic habitat assessment for these MPAs. Highlighting the contribution of this study also entails acknowledging its limitations: due to feasibility, not every habitat could be sampled with the same intensity. Still, this is an issue that should promote further studies on habitats still understudied, for which this study provides the first available information.

Over 80% of the species included in the Red Book of the Invertebrates from Romania are found in the infralittoral habitats (< 20 m) of the southern Romanian Black Sea coast (Murariu and Maican, 2021). Some of these species, including Dysidea fragilis, Halichondria panicea, Petricola lithophaga, Gastrana fragilis, Flexopecten glaber ponticus, Donax trunculus, Pholas dactylus, Palaemon adspersus, Eriphia verrucosa, had been already mentioned in various studies in the region (Begun et al., 2006; Teacă et al., 2006a,b, 2010; Filimon, 2020; Surugiu et al., 2021; Teacă, pers. communication). However, our study collected only a fraction of the species (5) included in the Red Book of the Invertebrates from Romania, likely due to the rather small surface of the investigated infralittoral bottoms.

To the best of our knowledge, our study provides the most detailed analysis of benthic diversity at Romanian littoral MPAs to date. This study also documents some of the highest diversity levels (a grand total of 119 taxa), considering that previous studies (i.e., Abaza et al., 2019) have reported slightly lower diversities (113 taxa) from an area encompassing all the Romanian benthic habitats. Our study is comparable with the studies conducted at the Bulgarian littoral, where fine-scale surveys mapped benthic habitats within the framework of the MSFD (Todorova et al., 2015, 2020; Nikolopoulou et al., 2021). Similar to our results, those studies revealed a good compliance between delineation of distinct benthic habitats and the presence of relevant engineering species contributing to biomass as well as both density and diversity. Such approach has served as a common assessment scheme across nations for the Black Sea EUSeaMap habitats mapping (Populus et al., 2017; Vasquez et al., 2020). It follows that several of the habitats identified within the two Romanian MPAs studied here, correspond with those described from the nearby Bulgarian coast (i.e., infralittoral sand dominated by Chamelea gallina, circalittoral mud with Spisula subtruncata, etc.). Such correspondence suggests the potential for upscaling the current study to assess regional interconnectivity among similar benthic habitats across the Black Sea region.

The Cape Tuzla and the Mangalia Sulfur springs MPAs are wide open coastal areas exposed to a variety of anthropogenic activities (Begun et al., 2018). Shipping, tourism, and pelagic and bottom fishing can individually or in combination cause environmental disruptions and likely have an impact on benthic habitats (Begun et al., 2018). In fact, due to their small spatial extent and/or their vulnerability to environmental stressors, several species inhabiting both the infralittoral and circalittoral habitats have received a protection status by being included in the afore-mentioned Red Book of the Invertebrates from Romania (Murariu and Maican, 2021). Invasive species may also represent a threat, depending on the nature of their interactions with resident (native) species (see Gurevitch and Padilla, 2004; Quijón et al., 2017). In this area, in particular, several invaders have become established and do represent a serious threat (Dipolydora quadrilobata, Polydora cornuta, Streblospio gynobranchiata, Corambe obscura, Anadara kagoshimensis, Mya arenaria, and Amphibalanus improvisus) to both vulnerable species and the diversity of this coastline. Among them, the latest invader with rapid spreading behavior detected in Romanian benthic habitats is the spionid polychaete Streblospio gynobranchiata (Teacă et al., 2021).

Even though those stressors amount to a certain risk for the diversity in the Romanian Black Sea, these coastal habitats seem more severely affected by eutrophication and industrial pollution (Todorova et al., 2020). Despite that, our study indicates that 88 and 60% of the habitats and associated communities in the infralittoral of the Cape Tuzla and the Mangalia Sulfur springs MPAs, respectively, are achieving good ecological status (GES). Meanwhile, only ∼54% of the habitats and communities associated with circalittoral bottoms are achieving GES. The latter result contrasts with the GES of other habitats from the circalittoral (namely from mud and sandy mud, shelly coarse sediment and mixed sediment) which are achieving between 70 and 100% GES due to their higher biodiversity. These differences suggest that benthic habitats are for the most part resilient to the anthropogenic stressors present in this area, but also imply that there are particular habitats that do require further attention. We therefore call for further research to be conducted on specific sources of stress (first and foremost organic and industrial pollution) and the response of sensitive species to short- and long-term exposure to these stressors.

Recognizing the importance of mapping of benthic habitats in MPAs such as those studied here in the southern Romanian Black Sea, would be a key step toward better managing these areas. The unique biological value and aesthetic quality of the areas currently protected by the MPAs cannot oppose the reasonable use of these waters. Habitat mapping constructed from the combination of seafloor characteristics and main benthic community features represents a useful tool and therefore should play a critical role on the management of these systems. Seafloor structure and the species that colonize these habitats, even at relatively small spatial scales (e.g., Surugiu et al., 2021), not only maintain the function of these systems, but also act as integrative indicators of their quality and health. It follows that seafloor structure and function should be an integral part of the framework to integrate research and management programs (Sullivan et al., 2017). Such approach would alleviate some of the concerns that the larger temporal and spatial scope of our study makes evident, and that should be effectively communicated to coastal managers.