Structural and Functional Analyses of Motile Fauna Associated with Cystoseira brachycarpa along a Gradient of Ocean Acidiﬁcation in a CO 2 -Vent System off Panarea (Aeolian Islands, Italy)

: Ocean acidiﬁcation (OA), one of the main climate-change-related stressors linked to increasing CO 2 concentration in the atmosphere, is considered an important threat to marine biodiversity and habitats. Studies on CO 2 -vents systems, naturally acidiﬁed environments that mimic future ocean scenarios, help to explore the sensitivity of species and to understand how benthic communities rearrange their structure and functioning under the pressure of OA. We addressed this problem by studying the benthic invertebrates associated with a habitat-forming brown alga ( Cystoseira brachycarpa ) in the Bottaro crater vents system off Panarea island (Tyrrhenian Sea, Italy), by sampling along an OA gradient from the proximity of the main venting area (station B3, pH 7.9) to a control zone (B1 station, pH 8.1). Samples were collected in September 2016 and 2018. A total of 184 taxa and 23 different functional traits have been identiﬁed, considering feeding habit, motility, size, reproductive and developmental biology, and occurrence of calcareous structures. Invertebrates are distributed according to the distance from the high venting zone and low pH levels and results very consistent between the two investigated years. In the low-pH area (B3), 43% of the species are selected. The functional traits of the fauna mirror this zonation pattern, mainly changing the relative proportion of the number of individuals of the various functional guilds along the OA gradient. Invertebrates inhabiting the low-pH zone are mainly composed of weakly or non-calciﬁed species, with small size, burrower/tubicolous habit, omnivorous or suspension feeders, and with direct development and brooding habit. In the other stations, heavily calciﬁed forms, herbivore and herbivore/detritivore, and with medium (1–5 cm) and large (>5 cm) sizes prevail, showing indirect benthic and planktic development. The taxonomic analysis, coupled with functional aspects, increases our prediction of which traits could be potentially more advantageous for species to adapt to the hypothesized scenarios of OA, and identify present and future winner and/or loser organisms in the future ocean of the Anthropocene. ordered ANOSIM and SIMPER tests were performed to assess the signiﬁcance of the OA gradient models and highlight the species or the traits mostly contributing to the average dissimilarities between groups. K–W ANOVA was performed with Statistica TM 12 software, while multivariate analyses were performed with PRIMER PERMANOVA 7.0 software.


Introduction
As the ocean continues to take up carbon dioxide (CO 2 ), marine waters are progressively decreasing their mean pH, increasing pH temporal variability as well as altering the complex carbonate chemistry, an overall phenomenon known as ocean acidification (OA) [1,2]. Experimental studies assessing the potential impacts of ocean acidification on marine organisms have rapidly expanded and produced a wealth of empirical data over recent decades. Although it is difficult to predict the future of marine species, communities, and ecosystems, from the first reviews and meta-analyses on OA, based initially mainly on laboratory and mesocosm' studies [3][4][5][6], it was clear the overall negative effects of this stress-forcing factor on many basic life processes of marine organisms, likely except for photosynthesis. From the late 2000s, in the OA research scenario, naturally acidified systems, such as CO 2 vents, started to take a major role in assessing the effects of OA at different hierarchical scales of the biological complexity. Natural CO 2 vents or seep systems are subtidal areas affected by the emission of gasses of volcanic origin, and they are generally rich in CO 2 , although other gasses (also toxic) can be present (e.g., sulfur, methane). They can be found where there is or has been underwater volcanic activity, such as mid-oceanic ridges or island arcs, or intraplate magmatism [7]. One of the first vent systems studied to assess the effect of OA on the benthic biota was the Castello Aragonese islet off the Ischia island [8,9], while to date at least other 23 CO 2 -vent systems have been studied to some extent [10]. However, the reviews on the research conducted in these systems have been provided only for the Castello of Ischia [9], the boreal northern areas [10], and other systems [11,12]. Natural CO 2-vent sites generally provide a pH gradient and mimic near-future scenarios, of OA and variability of the pH and carbonate chemistry alterations, that are useful as a proxy to investigate ecological effects of ocean acidification. The effects of decreased pH can be assessed at increasing levels of organization, from the responses of individual species up through populations and communities to the whole ecosystems. As natural laboratories, CO 2 -vent sites incorporate a range of environmental factors, such as gradients of nutrients, currents, and species interactions, that cannot be replicated in the laboratory or mesocosms, with the caveat that some vent systems have confounding factors such as hydrogen sulfide and metals (e.g., Vulcano [13]).
Biodiversity, both expressed as taxonomic diversity of species and diversity of functions of such species, is an intrinsic property and value of the natural communities and largely contributes to the main services that ecosystems provide to human beings and society. From research to date conducted in various vent systems, OA is deeply altering habitat health [14], and although the studies at community levels are still relatively few, a general decrease in taxonomic biodiversity has been highlighted in temperate as well as in tropical systems [4,[15][16][17][18][19][20][21]. As regards the functional diversity, the available literature includes studies focused on specific groups (e.g., polychaetes [22]) or addressing only a specific trait, such as the type of reproduction and egg development (e.g., [23]), adult size (e.g., [24]), or trophic relationships and food-web structure (e.g., [5,18]). In this respect, a recent study was conducted at the Castello vent system off Ischia by Teixidó et al. [25], focusing on the functional traits approach of benthic organisms of the rocky shallow reef. The functional diversity of sessile benthic assemblages was assessed using 15 traits (combined for 72 species for 68 unique combinations), which describe complementary facets of species ecology. The research highlighted that functional diversity decreases with acidification and that functional loss is more pronounced than the corresponding decrease in taxonomic diversity. Moreover, most organisms accounted for a few functional entities (i.e., unique combination of functional traits) and low functional redundancy in acidified conditions. All these observations suggest that ocean acidification will result in a simplification of the community [5,6] and its trophic complexity [18], up to a threshold level where it could no longer sustain any ecosystem function [25].
The Panarea island, although being the smallest of the Aeolian Islands, has geological and geomorphological characteristics that make it, and its surrounding complex archipelago of small islets, the largest hydrothermal system of the whole Mediterranean Sea. The area consists of numerous and diversified hydrothermal submerged emissions of gas and hot waters of volcanic origin, including a recently described site with more than 200 deep chimneys [26]. These features make this zone one of the most suitable places in the Mediterranean to assess the effects of water acidification on benthic communities, as highlighted by various studies on different habitat and vent systems in the area [26][27][28][29][30]. Among the hydrothermal vents, the crater off the Bottaro islet (East of Panarea island) is particularly interesting, since in this shallow area a gradient of OA is established from the center of a crater, where most of the venting occurs, to the external part [19,28]. Moreover, the rim of the crater and the surrounding bottoms are colonized by a dense settlement of the canopy-forming brown algae Cystoseira brachycarpa (J. Agardh, 1896). This macroalga, which acts as an ecosystem engineer forming 'forests' in infralittoral and circalittoral areas, represents one of the most important marine habitats in the Mediterranean Sea, forming extended canopies comparable to land forests and providing refuge and subsistence for many organisms, including commercial species [31,32]. On the other hand, this habitat is very vulnerable to different anthropogenic disturbances, both direct (pollution, increased sedimentation, coastal urbanization, direct habitat transformation) and indirect (overfishing, invasive species, climate change), as well as destructive fishing (date mussel fishery) [33].
In the hydrothermal area off Panarea, C. brachycarpa is a fairly common and widespread seaweed and forms quite extended forests [34]. This macroalga not only provides space, shelter, and food but, thanks to the photosynthetic process, may buffer the negative effects of OA for the small associated invertebrate fauna. Indeed, some experimental studies revealed that large macroalgae can buffer the effects of OA on calcification processes [35,36], and influence plant-animal interactions [37], thus exerting various direct and indirect effects on this specific stressor.
In light of this, the present study, conducted in the frame of the Scientific Diving Summer School hosted in the ECCSELNatLab Italy infrastructure of Panarea island [38], is focused on the benthic motile invertebrates associated with the habitat-forming brown algae Cystoseira brachycarpa in the Bottaro crater, to assess differences in taxonomic and functional biodiversity and the zonation pattern along the identified OA gradient.

Study Area and Sampling Activities
The study area is located off the Panarea island (Aeolian Archipelago) at the Bottaro crater close to the Bottaro islet (8-10 m depth) (38 • 38 13.58 ' N; 15 • 6 33.95 ' E) ( Figure 1). The area was upset by a submarine explosion of low energy and degassing between 2nd and 3rd November 2002 (with emissions mainly of CO 2, with traces of H 2 S and other gaseous species) that formed a submarine crater [39]. Nowadays, most of the venting activity is concentrated within the crater characterized by a main depression at 11 m depth, 14 m wide, and 20 m long. The sampling of invertebrates was carried out in September 2016 and September 2018, at 3 stations located along a gradient of OA, at depths from 8 to 10 m: station B3 was located at the rim of the crater, and most affected by OA (7.8 pH), station B2 at 7 m distance (7.9 pH), and station B1 at 35 m from B3 in a control area with ambient pH values [28] (Figure 1). At each station, water samplings were conducted in September 2016, 2018, and 2019 in order to characterize the carbonate system parameters along the pH gradient. Precisely, water samples, in 2016 and 2019, were collected by scuba divers with Niskin bottles to determine pH (spectrophotometrically, pH T ), total alkalinity (AT, by potentiometric titration in an open cell), inorganic nutrient concentration, and carbonate system parameters (partial pressure of CO 2 , dissolved inorganic carbon, and calcite and aragonite saturation state), while in 2018 only pH was measured with a laboratory pH-meter (pH NBS ).
At the three stations, benthic invertebrates were collected by scraping 20 × 20 cm surface (4 replicates × station × years: 8 replicates per station in total); each plot was visually selected by having at least approx. 50% of the surface covered by Cystoseira brachycarpa, in order to avoid algal-bared samples since the alga shows quite strong patchiness in the area.
Samples were fixed in 4% sea-water formalin and then transferred to 70% alcohol for sorting of the algae and the associated invertebrates. In the laboratory, invertebrates were sorted into main taxonomic groups and then identified at the lowest possible level (family, genus, and species) ( Table S1). Coverage and biomass (dry weight) of C. brachycarpa were also measured in each plot. For the functional analysis, each species of the identified invertebrates was assigned to a trophic category considering the food type and the modality of feeding. Seven trophic categories were considered: herbivores, both feeding on micro-and macroalgae (many molluscs and polychaetes), herbivore/detritivore (many peracarid crustaceans), detritivore (many decapods), omnivore (some polychaetes), omnivore/suspension feeders (various polychaetes), carnivore (including a few parasites) (molluscs, polychaetes, decapods). The trophic definition was assumed according to current literature [40][41][42][43][44][45][46]. Invertebrate species were also classified from their motility and life habit distinguished as: burrower/tubicolous (various amphipods), crawler (various polychaetes, amphipods, gastropods), crawler/climber (various gastropods), sessile (a few polychaetes, bivalves), and swimmers (various peracarid crustaceans, mysids). The species were also divided into three coarse size ranges based on the adult maximum size recorded in the literature: smaller than 1 cm, between 1 and 5 cm, and larger than 5 cm, and into three different calcification rates: heavily calcified, weakly calcified, and non-calcified. Finally, main reproductive and developmental mode of the species were scored as: "direct", when a species incubates the broods and a juvenile directly comes out of the eggs (e.g., many polychaetes); "indirect planktic", when external fertilization and free spawning are involved and a larva with a short or long life in the plankton or near the bottom is present; and "indirect benthic", when the development includes a larval stage with some pelagic life periods, but the eggs are deposited in the benthic environment as capsules (e.g., many molluscs) or gelatinous masses (molluscs). By the combination of the above traits, 23 different functional groups were identified; matrix traits/stations were built for each sampling year and the pooled data of the two years (see Table S2).

Data Analysis
Differences among stations were assessed with non-parametric Kruskal-Wallis ANOVA for algal biomass, number of invertebrate species, abundance, and diversity index (H' Shannon log 2 ); for community analysis, a multivariate approach was performed on the abundance matrix of species/stations (community structure) and functional traits/stations (community functions). The multivariate analyses were performed using Bray-Curtis dissimilarities, cluster analyses (group average), and nMDS ordination (data square roottransformed). Finally, ordered ANOSIM and SIMPER tests were performed to assess the significance of the OA gradient models and highlight the species or the traits mostly contributing to the average dissimilarities between groups. K-W ANOVA was performed with Statistica TM 12 software, while multivariate analyses were performed with PRIMER PERMANOVA 7.0 software.

Physico-Chemical Variables
The bottom seawater temperature measured in September 2016 evidenced no anomalies compared to the surrounding seawater and showed similar values along the transect, varying from 23.1 • C (station B3) to 22.9 • C (stations B1 and B2). Salinity was 39.5 at all sampled stations, while a decreasing pH gradient was detected towards the more impacted stations with pH T 25 • C values of 7.990, 7.857, and 7.785 at B1, B2, and B3, respectively. The same pattern was observed for the saturation rates of calcite and aragonite, while inverse patterns (with highest values measured in B3) were recorded for TCO 2 and pCO 2 ( Table 1). In September 2018, the pH measured with the pH-meter (pH NBS 25 • C) revealed no differences among the stations, while in September 2019 the measured pH T 25 • C showed a similar trend as in 2016. Other chemical variables did not show significant differences among stations or evident patterns along the pH gradient (Table 1).

Taxonomic and Structural Analyses
During both years surveys, a total of 20,606 individual macrobenthic invertebrates were sampled, belonging to 184 taxa, 146 of which were identified at the species level (see Table S1). The benthic community was composed of crustaceans (77 taxa and 11,625 individuals), mainly amphipods (45, 4587) followed by isopods (12,887), decapods (11,306), and tanaids (4, 5844); molluscs (55,1305); polychaetes (46, 6629); pantopoda (9, 544), and echinoderms (2, 503). Forty species occurred with a single individual; station B1 hosted 127 taxa (69% of the total), B2 129 (70.1% of the total), and B3 104 (56.5% of the total). Stations B1 and B2 shared 88 taxa, B2 and B3 shared 81 taxa, and B1 and B3 shared 77 taxa. Station B1 showed 28 exclusive taxa, B2, 27, and B3, 17 taxa, respectively; all the exclusive taxa however occurred with low numbers of individuals (between one and nine), so the majority of the most abundant species occurred in all the three stations along the pH gradient. Among the most abundant species (occurring with more than 200 individuals) many are amphipods, such as Ampithoe ferox, A. helleri, A. ramondi, A. riedli, Apherusa vexatrix, Cressa cristata, Protohyale camptonyx, Pereionotus testudo, Stenothoe elachista, and Caprella acanthifera. Other abundant crustacean species are the isopod Janiridae and the tanaids of the Leptochelidae family, then polychaetes, such as Odontosyllis ctenostoma, Exogoninae spp., Amphiglena panarensis, and A. aeoliensis (two new species recently described [47]). Finally, there were the molluscs Granulina marginata, Rissoa variabilis, and Musculus costulatus, and the small echinoderm ophiuroid Amphipholis squamata. The trends along the OA gradient of the algal biomass, number of species, individuals, and Shannon diversity (H') for years 2016 and 2018 are shown in Figures 2 and 3; the pooled values for both years are displayed in Figure S1. The biomass of C. brachycarpa in both sampling years (Figure 2A, Figure 3A) shows statistical differences among the stations (ANOVA K-W: H = 6.73; p < 0.05 for 2016; H = 7.54; p < 0.05 for 2018), with the lowest values in both years recorded at station B3, characterized by the highest acidification conditions. A similar pattern is observed for the number of species showing significant differences along the gradient in both years (ANOVA K-W: H = 6.42; p < 0.05 for 2016, K-W: H = 7.73; p < 0.05 for 2018), with the lowest species levels in B3 ( Figure 2B, Figure 3B). In terms of abundance, there are no significant differences along the gradient, although in both years higher values occurred at both acidified stations B2 and B3 ( Figure 2C, Figure 3C). Finally, the diversity index H' clearly mirrors the trend of the species richness, with significant differences along the pH gradient (ANOVA K-W: H = 6.42; p < 0.05 for 2016; K-W: H = 7.73; p < 0.05 for 2018) and with much lower values at B3 ( Figure 2D, Figure 3D).  Regarding the multivariate analyses, both the cluster analyses ( Figure S2) and the nMDS (Figure 4) showed the samples of each sampling year, as well as the pooled years, clustered by stations and then distributed along the acidification gradient. This change in faunal composition, with increasing distance from the acidified zone, is supported by the ordered ANOSIM test (R = 0.69; p < 0.01 for 2016; R = 0.90; p < 0.001 for 2018; R = 0.294; p < 0.01 for the 2016 + 2018 pooled data). The SIMPER analysis was conducted for the 2016 and 2018 data pooled, since analysis of every single year gave very similar results. The taxa mostly contributing to the pooled nMDS pattern (Figure 4) were Leptocheliidae and Amphiglena panarensis, more abundant in the most acidified station (B3), A. aeoliensis, more abundant in the weakly acidified station (B2), and Cressa cristata, more abundant in B1 and completely absent in B3 (Table S3).

Analyses of Functional Traits
The distribution of the functional traits/groups along the gradient are shown in Figures 5 and 6. In each graph, the relative dominance (percentage) of the number of individuals and the number of species for each category were plotted against stations along the pH gradient. All the identified functional traits (23) are present at the three stations, and what changes is their relative proportions (Table S2). As regards the trait of calcification (Figure 5a), station B3 is characterized by numerical dominance of non-calcified forms and the lowest presence of heavily calcified organisms; in the other stations, the proportion of the various forms is almost equal. Concerning the number of species, the proportion of the various forms is almost equal among stations. Analyzing the size distribution (Figure 5b), station B3 is numerically dominated by small species (38%), while larger taxa prevail at both stations B2 and B1 (58% and 52%, respectively). As regards the number of species, the proportion of the different sizes was almost equal among stations.  Looking at the distribution of the main trophic guilds (Figure 6a), station B3 is characterized by a larger proportion of individuals belonging to omnivore and omnivore/suspension feeder forms; other guilds (detritivore, herbivore, herbivore/detritivore, and carnivore/parasite) show higher proportions at stations B2 and B1. The number of species instead shows a more homogenous distribution among stations. The traits related to motility (not shown) reveal a large numerical and species dominance of burrower/tubicolous forms (53% of the whole abundance), followed by crawler/climber (39%). The former guild is mainly represented by many amphipods and other peracarid crustaceans and sabellid polychaetes. Burrower/tubicolous are largely dominant at B3 (58%) and less represented at stations B2 and B1, where other guilds are also relatively well represented such as crawler/climber, sessile, and swimmers. Finally, the analyses of the reproductive traits (Figure 6b) reveal that station B3 is characterized by the numerical dominance of direct brooders (38%), which together with indirect benthic (20%) account for almost 60% of all reproductive types. Brooder and indirect benthic are also quite abundant at station B2, while at B1 the indirect planktic is the dominant guild (55%).
The multivariate analyses, based on the combination of functional traits (see Methods), are shown in Figure 7. The distribution pattern of stations along the OA gradient mirrors quite well the pattern observed for the structural analyses in both years, as well as for the pooled data (2016 + 2018) (Figure 4). Stations are distributed according to the pH gradient, and these differences in functional traits dominance with increasing distance from the acidified zone are also confirmed by an ordered ANOSIM test (R = 0.68; p < 0.02 for 2016; R = 0.64; p < 0.01 for 2018; R = 038; p < 0.01 for the 2016 + 2018 pooled data). Although the correlation coefficients of the ANOSIM test are lower than those obtained from the structural analyses, they are significant and still explain a larger proportion of the natural variability of the data.  (Table S4). The traits mostly contributing to the pooled nMDS pattern (Figure 7; Table S4) are the small size, the direct/brooder, the weakly calcified forms, omnivore, herbivore, and burrower/tubicolous as regards the motility.

Discussion
The main physico-chemical features measured at Bottaro crater confirmed that the different identified zones are characterized mainly by decreasing mean pH and increasing concentrations of pCO 2 and TCO 2 , moving from the control station (B1) to the crater rim (B3); even the control station, due to its relative proximity to vents, showed in some of the measures slightly acidified values. The intense CO 2 emissions inside the Bottaro crater also show direct consequences for aragonite and calcite saturation rates: aragonite saturation values increased rapidly with distance from the crater, showing low saturation conditions (Ω Ar < 2.5) at the crater rim (B3) compared to the supersaturated waters (Ω Ar > 2.5) of station B1. These values are comparable to those reported by Kroeker et al. (2011) [4] for the ambient and low-pH zones of the Ischia Castello vent system. Here, these authors assessed in situ three main pH zones (ambient, low, and extreme low), which differed in the mean and variability of seawater pH, and related factors, along a continuous gradient. Comparable values have been also observed for the calcite saturation state. These differences among the identified zones are mirrored by distinct benthic communities associated with the vegetated hard bottoms. Our results confirm that the organisms residing at stations B2 and B3 have been exposed to low pH and low aragonite and calcite saturation.
The implications for this are discussed below as regards the taxonomic and functional trait distribution of the benthic invertebrates.
The total species diversity observed (184 taxa) was unexpected and fairly high, especially considering the relative homogeneity of the habitat, albeit coupled with high algal patchiness, as well as the scarce presence of megabenthic species, as highlighted by other studies based only on visual census techniques [28]. Therefore, habitat-forming brown alga, C. brachycarpa, and its complex architecture and food and shelter provision, can play an important role in the biodiversity recorded in our samples.
In both the examined sampling years, the distribution of benthic populations is very consistent with the OA gradient, and changes in structure and composition of faunal assemblages, as well as relative dominance of the functional traits, are clearly expressed along the different stations and mainly related to the gradient. In station B3, the benthic diversity (number of species and H' index) was 43% lower compared to the total number of identified species (184 taxa), while at station B2 and B1 the observed species reductions compared to the total were of 30% and 29%, respectively. Species occurrence seems to be affected by both the harsher conditions, due to lower pH levels, coupled with the reduction in habitat complexity, due to lower coverage and biomass of Cystoseira. Species reduction/selection along the OA gradient have been observed also in the Ischia vents, where the higher variability and much lower pH values than that measured in the Bottaro crater determine a reduction of 67% in total biodiversity at extremely low-pH conditions (>7.4 units) and a reduction of 52% at low pH (7.8-7.9) [9,48]. However, our analyses evidence that, although algal biomass, number of species, and diversity index H' decrease towards the low-pH zone, the number of individuals shows different trends, with abundances comparable among stations and even higher in weak and low-pH zones compared to the control (station B1). This is mainly due to the redundancy of a few species, which were more resilient to OA, and species that, due to low competition with other taxa, can develop denser populations (e.g., A. panarensis, Lepocheliidae, Exgoninae; various amphipods; Table S1). A similar trend was observed in Ischia vents' hard-bottom communities, where Kroeker et al. [4] recorded a strong community rearrangement and compensation due to the redundancy of few species. In addition, some of these abundant species at the Ischia vent vegetated bottoms are similar or belong to similar genera as our samples [4,22,[49][50][51]. Furthermore, benthic invertebrates associated with vegetated habitats in the Vulcano vents (macroalgae and the seagrass Cymodocea nodosa) show a very similar pattern, with a community dominated just by an annelid polychaete (Platynereis cf massiliensis) and few amphipods [18]. Higher abundances were also observed in invertebrates associated with Posidonia oceanica meadows under OA conditions at the Castello vents off Ischia [51,52]. A similar pattern was observed in tropical reef vents, where few corals species still dominated in the low-pH conditions [16]. On the other hand, a different pattern was shown by the cryptic associated fauna (studied with the Autonomous Reef Monitoring Structures, ARMS), where other factors seem to hinder the direct influence of low pH, also for calcifying organisms, some of which (e.g., molluscs) seem more abundant under OA than in the control area [20]. Finally, invertebrates associated with hard bottoms in a Japan vent system showed a selection of 30% of intertidal taxa and 60% of subtidal taxa against high-CO 2 conditions (low pH), and most of the selected species under OA conditions are corals, bivalves, and gastropods [21].
In our case study, other than the high resilience of some well-adapted species and reduced competition, a pivotal role is exerted by the buffer capacity of the macrophytes against the too-strong negative effects of low pH on the associated fauna. Indeed, the patchy distribution of the motile invertebrate community seems to reflect the high patchiness of the algal cover observed across the pH gradient, due to the morphology of the bottom, which is an irregular mosaic of rocks and Posidonia patches. pH buffering by macroalgae and seagrasses, which may act as "chemical refugia" from climate change and OA [53], has been observed and hypothesized in various studies [54,55], especially with respect to small invertebrates [35][36][37]51,56].
The dispersion of sample replicates, shown in the nMDS models space, was quite similar among stations, and also consistent in the different sampling years. This aspect is likely related to the high algal cover patchiness, although patchiness observed under OA conditions in another vent system showed a decrease towards a more homogeneous and uniform community setting [5]. However, it is also worth mentioning the recording of two taxa new to science, the two sabellid polychaetes Amphiglena panareensis and A. aeoliensis [47]. These species have been to date recorded only in this area, while in other nearby vents, such as the Vulcano vents (located on the Vulcano island approx. 12 nautical miles from Panarea) and the vents off Ischia (Castello) and Pozzuoli (Secca delle Fumose; Gulf of Naples), other new species of the same genus were recorded [40]. This discovery is extremely interesting since these new species may be related to the unique conditions of these shallow vents, which may also represent evolutionary laboratories. However, further collections outside the vents could be necessary to demonstrate the vent-specific occurrence of these new taxa, since to date, in the available literature, no shallow hydrothermal vents exclusive species have been described.
Regarding functional diversity, the distribution of each considered functional trait shares a common feature, where the proportion of species belonging to each different trait group does not show relevant changes along the OA gradient. In other words, each trait/function identified contained a quite similar number of species along the OA gradient, and there is no loss of species redundancy. What strongly changes among traits in different stations is their relative abundance (dominance). Heavily calcified taxa (e.g., molluscs, echinoderms), for instance, show a clear decrease in abundance in the most acidified station, B3 (<20%), while non-calcified taxa (e.g., polychaetes) are more abundant (38%) in the control station B1. Our results on chemical characterization along the gradient, showing low pH and low aragonite and calcite saturation in stations B3 and B2, is consistent with this pattern. In general, reduction in calcifying organisms under more severe OA conditions is a common trait observed in various CO 2 -vent systems [11], and it has been recorded both in benthic populations of hard-bottoms and seagrasses as well as on softsubstrates [4,11,18,21,57], but also in plant-dominated habitats and their epiphytes [58,59].
Regarding the size, there is a clear pattern, with smaller species dominance (<1 cm) (38% of total abundance) in the proximity of the crater (B3), while at the other stations, larger-sized taxa (>5 cm) prevail (52% in B1 and 58% in B2). In general, invertebrates associated with macroalgae are composed of small, cryptic species that commonly find refuge and resources to live in the algal thalli and where they often reproduce. In our case, the large abundance of small organisms observed at station B3 is due to the redundancy of a few small species, such as Amphiglena spp. Leptochelidae, and Exogoninae. Size constraint in CO 2 vents has been observed in polychaetes [22,60] and molluscs [24]. However, in our case, small-sized species seem selected by the algal architecture and the cryptic needs, more than by size or growth-limitation due to low pH.
Small-sized invertebrates associated with macroalgae are also often herbivores and herbivore-detritivores (the so called mesoherbivores), or suspension feeders. Indeed, these trophic guilds are the most abundant in our samples, together with omnivores ( Figure 6a). The pattern of feeding guild distribution shows the numerical dominance of more generalistic habits (omnivores and suspension feeders) at the most acidified station (B3), while the other guilds prevail in the less acidified stations. Furthermore, this pattern seems somehow consistent with previous observations in other vent systems which highlight a general simplification and a prevalence of low-level food web consumers under more severe OA [18,21,61], and in general a functional trait oversimplification [25]. However, in these studies, functional simplification is coupled with a strong reduction in biodiversity and species redundancy within each trait that, as stated above, does not occur in our system where the number of species of each trophic guild is maintained quite stably across the stations, although with limited abundances.
Finally, the analysis of the main reproductive/developmental type clearly shows the prevalence of direct development and brooding habits in the most acidified zone, although indirect species which deposit eggs on the bottoms are also relatively abundant. Prevalence of brooders has been observed in invertebrates associated with other vents as the Ischia-Castello vents, especially in polychaetes [22,23], as well as at the Vulcano vents, where invertebrates are dominated by the brooding nereid polychaete Platynereis cf massiliensis [18]. This pattern has been related to the fact that brooding and bottom egg deposition can favor the success of offspring since they develop in a more sheltered environment where macrophyte photosynthesis can buffer the low pH [23]. However, since brooding is often a covariate trait with small dimensions, the dominance of brooders could also be related to the dominance of small-sized species.

Conclusions
Invertebrates associated with the habitat-forming macroalga Cystoseira brachycarpa, occurring at Bottaro crater off Panarea, were clearly distributed according to the distance from the high venting zone and the OA gradient. This macrofaunal distribution appears similar and very consistent between the two investigated years. In the low-pH area, around the crater rim, 43% of the species are selected, while the high algal patchiness is responsible for the variability of samples, and is comparable among stations. The faunal functional traits mirror this zonation pattern, where each of the considered traits changes in terms of relative numerical dominance along the OA gradient. Invertebrates inhabiting and more abundant in the low-pH zones can be considered winners, and to have more chance to cope with future marine conditions. These are mainly composed of small-sized (<1 cm), non-calcified or weakly calcified species (mainly polychaetes and amphipods), with burrower/tubicolous and omnivorous or omnivorous/suspension feeder habits, and with direct development and brooding habit. At the other stations, medium and large-sized, heavily calcified forms (e.g., molluscs), herbivores and herbivore-detritivores prevail, showing an indirect benthic and planktic development. These, according to studies conducted in other hydrothermal systems, appeared to be the most vulnerable (losers) to the acidified conditions foreseen for the future.
The taxonomical analysis, coupled with the functional aspects, increases our prediction of which traits could be potentially more advantageous for the species to adapt to the hypothesized OA scenarios, and identify actual and future winner and/or loser organisms in the future Anthropocene ocean.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/jmse10040451/s1, Figure S1: community structure variable of 2016-2018 data pooled; Figure S2: cluster analyses of taxonomic data of 2016 and 2018 samples; Table S1: list of species and their functional trait features; Table S2: matrix of functional trait/stations; Table S3: SIMPER analyses of structural data from 2016 and 2018 pooled data; Table S4: