Seawater Acidiﬁcation A ﬀ ects Beta-Diversity of Benthic Communities at a Shallow Hydrothermal Vent in a Mediterranean Marine Protected Area (Underwater Archaeological Park of Baia, Naples, Italy)

: One of the most important pieces of climate change evidence is ocean acidification. Acidiﬁcation e ﬀ ects on marine organisms are widely studied, while very little is known regarding its e ﬀ ects on assemblages’ β -diversity. In this framework, shallow hydrothermal vents within a Marine Protected Area (MPA) represent natural ecosystems acting as laboratory set-ups where the continuous carbon dioxide emissions a ﬀ ect assemblages with consequences that can be reasonably comparable to the e ﬀ ects of global water acidiﬁcation. The aim of the present study is to test the impact of seawater acidiﬁcation on the β -diversity of soft-bottom assemblages in a shallow vent ﬁeld located in the Underwater Archeological Park of Baia MPA (Gulf of Naples, Mediterranean Sea). We investigated macro- and meiofauna communities of the ‘Secca delle fumose’ vent system in sites characterized by sulfurous (G) and carbon dioxide emissions (H) that are compared with control / inactive sites (CN and CS). Statistical analyses were performed on the most represented macrobenthic ( Mollusca , Polychaeta , and Crustacea ), and meiobenthic ( Nematoda ) taxa. Results show that the lowest synecological values are detected at H and, to a lesser extent, at G. Multivariate analyses show signiﬁcant di ﬀ erences between hydrothermal vents (G, H) and control / inactive sites; the highest small-scale heterogeneities (measure of β -diversity) are detected at sites H and G and are mainly a ﬀ ected by pH, TOC (Total Organic Carbon), and cations concentrations. Such ﬁndings are probably related to acidiﬁcation e ﬀ ects, since MPA excludes anthropic impacts. In particular, acidiﬁcation markedly a ﬀ ects β -diversity and an increase in heterogeneity among sample replicates coupled to a decrease in number of taxa is an indicator of redundancy loss and, thus, of resilience capacity. The survival is assured to either tolerant species or those opportunistic taxa that can ﬁnd good environmental conditions among gravels of sand.


Introduction
All biological systems have general properties that, coherently with environmental conditions, modify their structures in order to maximize resources exploitation tending to an equilibrium state (climax). In marine realms, benthic assemblages are heavily modified since sessile and sedentary species passively endure environmental shifts throughout their whole life [1,2]. Indeed, benthic species are often characterized by slow growth rates and low movement capacity that allow them to record environmental modifications taking place in a given area [3,4]. β-diversity is an assemblages gradient coherent with environmental variations coming from interactions between ecological features (e.g., species interactions, niche width, etc.) and edaphic abiotic characteristics [5] that need to be investigated in order to detect possibly recurring communities patterns.
Coastal marine ecosystems are under continuous pressure due to anthropic disturbances that determine climate changes effects (ocean warming, acidification, sea-level rise, and the increasing frequency and intensity of storms) [6]. Many coastal ecosystems are composed of species that are critical for supporting biodiversity, ecosystem function [7], and a suite of critical ecosystem services [8,9]. Nowadays, they are at risk of disappearing since climate changes are often faster than species/communities resilience capacity [10][11][12][13] and they have no time to modify themselves in order to reach a climax state. Investigations on these systems are often a challenge since it is difficult to disentangle the effects of climate change and anthropic impact on assemblage structure. In the case of Marine Protected Areas, we can exclude anthropic activities as a cause of assemblage modifications.
Small-scale heterogeneity is a measure of β-diversity [14,15], defined as the variations in the community composition among a set of sample units within a given spatial extent and without referring to any particular direction [16,17]. To this respect, how much a community is heterogeneous provide indications of its resilience capacity [18] since assemblages change their patterns in space according to the environmental drivers [19]. Assemblages particularly subjected to frequently disturbed environments might be more resilient than others [20] being composed of taxa that can cope with unstable environmental conditions (e.g., short life cycles, copious offspring, and high turnover rate) [21]. However, resilience capacity is also strictly related to species richness (α-diversity) and functional diversity that thus become measures of ecological stability [22]. In the natural systems, the higher the number of species within an ecosystem, the higher the likelihood that many ecological functions occur to stabilize the ecosystem [23]. In addition, where ecological functions of different species overlap (redundancy), these functions may persist also in the case of species removal; the loss of redundancy impairs the ecosystem ability to withstand disturbance, thus eroding its resilience [24].
Acidification is a reduction in the pH values of the ocean over an extended period, primarily caused by the uptake of carbon dioxide (CO 2 ) due to anthropogenic emissions (where it increases as an effect of climate changes) and to natural (mainly volcanic) activities. The latter is the case of hydrothermal vents in volcanic areas that within MPAs may be useful as natural laboratories to study microscale climate changes effects. Indeed, their gas emissions, when composed of CO 2 , might determine local seawater acidification affecting biological and ecological features [25]. Sea water acidification determines community modifications and loss of redundancy at α-diversity level [26][27][28]. However, to date many studies on acidification effects have considered target species and benthic community structures (i.e., [29][30][31][32] etc.) and, as far as we know, no attempt to assess the effects of this phenomenon on β-diversity has been made.
In the present study, we studied β-diversity as small-scale spatial heterogeneity of meio-and macrobenthic communities inhabiting the soft bottoms of a shallow hydrothermal vents area, and tested if we can consider β-diversity a measure of resilience capacity in acidified seawater environment. The investigation was carried out within the Underwater Archeological Park of Baia MPA located in the Gulf of Naples, where environmental physical parameters and community modifications are possibly less affected by anthropic activities.

Study Area
The study area is located in the north-west side of the Gulf of Naples (Italy, south Tyrrhenian Sea) within the Underwater Archeological Park of Baia MPA. Since 2002 the MPA protects the coastline consisting of beaches and tuff cliffs shaped by volcanic activity, evident for both bradyseism phenomena and hydrothermal vents [33]. The interaction between natural processes and human activities produced a natural environment characterized by a strong interdigitation between photophilous and sciaphilous populations (AP-C). Superficial muddy sands in sheltered-water habitats (SVMC), enriched with fragments of carbonate exoskeletons of benthic organisms, are common in the sampling area. Coarse sands and fine gravels mixed by the waves habitat (SGCF) [34,35] is also present below 15 m depth.
The sampled area, named "Secca delle Fumose" (Figure 1a), is characterized by two vent typologies, and environmental parameters featuring each sampled site are shown in Table 1 [36,37]. The Geyser (G) site is a vent characterized by water emission reaching~80 • C at the outlet, while at a distance of 20 cm from the center the sediments temperature is around 30 • C with a pH of 8. Rocky substrate surrounding the geyser are covered by yellow sulfurous deposits (Figure 1b), soft bottoms among the rocks are characterized by around 30% of TOC, and S 2− , SO 4 2− , Mg 2+ , Zn 2+ , and Pb 2+ are the most abundant ions detected in their interstitial water. The hydrothermal (H) vent site is characterized by carbon dioxide (CO 2 ) bubbles and its sediments show an almost homogenous temperature of~37 • C with a pH of 7.5. This site is characterized by a wide, white microbial mat (Figure 1c) on both rocky and sand bottoms, determining~35% of TOC and Na + , NO 3 − , and Pb 2+ are the most abundant ions detected in their interstitial water.
At the control/inactive sites located at the north (CN) and at the south (CS) of the study area, approximately 100 m from respectively G and H, water and gases emissions are not present, sediments temperature is around 22 • C with a pH of 8.1. Sediments are characterized by around 17.5% of TOC and the most abundant ions in interstitial water are NO 3 − , SO 4 2− , Zn 2+ , and Pb 2+ .

Sampling Campaign and Data Analysis
A macro-and meiofauna sampling campaign was performed in October 2016 in H, G, CN, and CS sites. All samples were collected in the same habitat: SVMC enriched with fragments of carbonate exoskeletons of benthic organisms.
Macrofauna samples were obtained at each site by scuba-diving through an air-lift pump equipped with a 0.5 mm nylon mesh size bag within a 50 cm × 50 cm frame, reaching a depth of 10 cm in the sediment (for further information on methodology see [36]).
Meiofauna samples were collected at each site by scuba-diving through cylindrical corers with a radius of 2.75 cm and reaching a depth of 10 cm in the sediment (for further information on methodology see [37]).
Mollusca (M), Polychaeta (P), Crustacea (C) (hereafter called macrofauna), and Nematoda (N) (hereafter called meiofauna) assemblages were the main taxa considered in the present study, being representative of more than the 80% of the total sampled organisms [36,37]. They were analyzed separately in order to describe each assemblage and to assess their β-diversity (measured as small-scale heterogeneity among sample replicates).
Descriptive analysis takes into account synecological indices: α-diversity, measured as species/taxa richness (α M, P, C = n • of species on 25 dm 3 ; α N = n • of taxa on 0. Distance-based permutation univariate analysis of variances (PERMANOVA) [41] was based on synecological indices' Euclidian distances [42] and each term in the analysis was tested by 4999 random permutations. Post hoc pair-wise comparisons using the PERMANOVA t-statistic and 4999 permutations were also conducted if necessary.
PERMANOVA multivariate analyses were carried out on abundances, in order to test for differences in assemblage structure due to environmental characteristics of each sampled site [experimental design involved one factor Site (Si, fixed, four levels) with n = 3]. Prior to the analyses data were log(x+1) transformed to reduce the weight of scanty taxa. All multivariate analyses were based on Bray-Curtis similarity, and each term in the analysis was tested by 4999 random permutations [43,44]. Post hoc pair-wise comparisons using the PERMANOVA t-statistic and 4999 permutations were also conducted if necessary. Multivariate patterns were visualized through Canonical Analyses of Principal coordinates (CAP) [45] of sites' elements. Permutation tests of multivariate dispersion (PERMDISP) [46], a measure of β-diversity based on Bray-Curtis similarity matrices, were also performed for the term Si to investigate possible effects of environmental characteristics on small-scale heterogeneities among assemblages (heterogeneities are measured as mean ± SD distance of sites replicates from the centroid of the cluster they form). All analyses were performed using r-project software [47] with BiodiversityR package [48].
Environmental data measured at each site and reported in Table 1 [36], are here related to β-diversity measures of macro-and meiofauna assemblages. In particular, a distance-based linear modelling (DistLM) [49] using step-wise as selection procedure and adjusted R 2 (hereafter called only R 2 ) as selection criterion were performed in order to find correlations among the distances from the centroids in each site and environmental data. Relations between assemblage heterogeneities and environmental variables that significantly affect β-diversity were visualized through distance-based redundancy analysis (dbRDA) [50].
Finally, environmental parameters that significantly affect β-diversity were related to synecological indices through a linear regression model in order to assess their effects on assemblages structures. For each pair composed by environmental parameters (independent variable)/synecological index (dependent variable) the R 2 was calculated to find possible correlations among synecological indices and environmental parameters.

Mollusca Assemblages
The list of Mollusca collected in the four sites is shown in Supplementary Material (A). The highest numbers of Mollusca species and individuals were detected at CS (α CS = 28.33 ± 4.76 species/25 dm 3 ; A CS = 106.67 ± 52.50 individuals/25 dm 3 ) while the lowest were found at H (α H = 1.33 ± 2.15 species/25 dm 3 ; A H = 3.33 ± 2.89 individuals/25 dm 3 ). α-diversity showed significant differences (p < 0.001) between hydrothermal vents sites and control/inactive sites, while abundances showed significant differences (p < 0.05) between H and the other sites (Figure 2i). The diversity index showed significant differences between the hydrothermal vents sites and control/inactive sites (p < 0.01) while evenness showed a significant difference between CN and G (p < 0.05) (Figure 2ii).
PERMANOVA analyses show significant differences between sites (p < 0.01); in particular the pair-wise test highlights significant differences between: H and control/inactive sites (p H,CN = 0.043; p H,CS = 0.031), and G and CN (p = 0.043). CAP analysis shows four different clusters separating the four levels composing the factor Site ( Figure 2iii). PERMDISP does not show significant differences between small-scale heterogeneities of Mollusca assemblages but H is characterized by the greatest mean distance (d = 46.121 ± 9.893), while the lowest value belongs to the CS site (d = 23.551 ± 1.594) (Figure 2iv).
In Figure 3 correlations between environmental parameters and synecological indices that significantly affect Mollusca heterogeneities are shown. In particular, pH is always directly related to synecological indices showing the highest correlation with H' (R 2 = 0.75); TOC is always inversely related to synecological indices showing the highest correlation with α (R 2 = 0.62) and J (R 2 = 0.6); Na+ is always directly related to synecological indices showing the highest correlation with α (R 2 = 0.87) and H' (R 2 = 0.92).

Polychaeta Assemblages
The list of Polychaeta collected in the four sites is shown in Supplementary Material  . α-diversity showed significant differences (p < 0.01) between hydrothermal vents sites and control/inactive sites except between G and CS; while abundances showed significant differences (p < 0.001) between all sites except between those characterized by hydrothermal vents (Figure 4i). PERMANOVA analyses show significant differences between sites (p < 0.01) although pair-wise test highlights significant differences only between the CN and CS sites (p = 0.012). CAP analysis shows four different clusters separating the four levels composing the factor Site ( Figure 4iii). PERMDISP does not show significant differences among small-scale heterogeneities of Polychaeta assemblages but site H is characterized by the greatest mean distance (d = 57.735 ± 0.1), while the lowest value was recorded at site CN (d = 22.527 ± 1.568) (Figure 4iv).
In Figure 5, correlations between environmental parameters and synecological indices that significantly affect Polychaeta heterogeneities are shown. In particular K + is always directly related to synecological indices although with low value of R 2 for which the highest is related to H' (R 2 = 0.48); Mg 2+ is directly related to α, A, and H' showing almost the same R 2 (about 0.6) and inversely very weakly related to J (R 2 = 0.037); Zn 2+ is directly related to α, A, and H', showing the highest correlation with α (R 2 = 0.63), while a very low inverse relation is detected with J (R 2 = 0.069).

Crustacea Assemblages
The list of Crustacea species collected in the four sites is shown in Supplementary Material (C). The highest number of Crustacea species was recorded at CN (α CN = 9.66 ± 1.52 species/25 dm 3 ) and the highest number of individuals was found at CS (A CN = 28 ± 19.67 individuals/25 dm 3 ), while both the lowest number of species and individuals were detected at H (α H = 0.67 ± 0.58 species/25 dm 3 ; A H = 0.67 ± 0.57 individuals/25 dm 3 ). α-diversity showed significant differences (p < 0.05) between H and the control/inactive sites, while abundances exhibited significant differences (p < 0.05) only between the H and CN sites (Figure 6i).
Diversity and evenness indices were not calculated in the H site, since only one single species with one individual was recorded in two replicates. The highest value of diversity index was detected in CN (H' CN = 2.01 ± 0.19) and the lowest in G (H' G = 2.01 ± 0.19). The highest value of evenness was detected in G (J G = 0.96 ± 0.02) and the lowest in CN (J CN = 0.88 ± 0.03). Neither diversity and evenness indices sowed significant differences between sites (Figure 6ii).
PERMANOVA analyses were performed only among G and control/inactive sites since the low number of Crustacea species and individuals at site H did not allow the use of Bray-Curtis similarity.
The G and control/inactive sites show significant differences (p < 0.01; in particular p G,CN = 0.043; p G,CS = 0.033) as confirmed by CAP analysis that shows two clusters: one polarized on the positive area of CAP2 and composed by G replicates, and one polarized on the negative area of CAP2 and composed by control/inactive site replicates (Figure 6iii). PERMDISP does not show significant differences between small-scale heterogeneities of Crustacea assemblages, but the G site is characterized by the greatest mean distance (d = 44.828 ± 2.108), while the lowest value belongs to the CS site (Figure 6iv).
In Figure 7, correlations between environmental parameters and synecological indices that significantly affect Crustacea heterogeneities are shown. In particular, pH is directly related to α, A, and H', showing almost the same R 2 (about 0.5), and inversely related to J (R 2 = 0.51); T is inversely related to α, A, and H', showing the highest correlation with H' (R 2 = 0.73), and directly related to J (R 2 = 0.23); Na + is directly related to α, A, and H', showing almost the same R 2 (about 0.69), and inversely related to J (R 2 = 0.038).   . α-diversity showed significant differences (p < 0.05) between H and the other sites, while abundances did not show significant differences (Figure 8i).
PERMANOVA analyses show significant differences between sites (p < 0.01)-in particular, the pair-wise test highlights significant differences between H and other sites (p H,CN = 0.0054, p H,CS = 0.007, p H,G = 0.003). CAP analysis shows three clusters: one polarized on the negative area of both CAP1 and CAP2, composed by H replicates; one polarized on the positive area of CAP1 and the negative area of CAP2, composed by G replicates; and one polarized on positive area of both CAP1 and CAP2, composed by control/inactive site replicates (Figure 8iii). PERMDISP does not show significant differences between small-scale heterogeneities of Nematoda assemblages whose distances are similar among the four sites (Figure 8iv). DistLM shows that Nematoda heterogeneities at sites H and G are significantly and inversely correlated to pH (p = 0.028, R 2 = 0.66) and to the % of mud (p = 0.008, R 2 = 0.54), and directly correlated to Pb + (p = 0.011, R 2 = 0.737) and temperature (p = 0.032, R 2 = 0.79), which explain, altogether, about 65% of the total variation (Figure 8v).
In Figure 9, correlations between environmental parameters and synecological indices that significantly affect Nematoda heterogeneities are shown. In particular pH is always directly related to synecological indices showing the highest correlations with α (R 2 = 0.71), H' (R 2 = 0.80) and J (R 2 = 0.79); T is always inversely related to synecological indices showing the highest correlations with α (R 2 = 0.76), H' (R 2 = 0.80), and J (R 2 = 0.75); Pb 2+ is inversely related to to synecological indices showing the highest correlations with J (R 2 = 0.49). Mud always shows a very weakly correlation with synecological indices.

Effects of Acidification on β-Diversity
The present study aims at detecting the possible effects of hydrothermal vents, located within an MPA, on β-diversity of soft-bottom meiobenthic and macrobenthic assemblages. Faunistic results from the sediment samples collected at each site were divided in four main sub-assemblages that represent the vast majority of macrobenthic communities (Mollusca, Polychaeta, and Crustacea populations) and meiobenthic communities (Nematoda populations). Coherently with previous studies [36,37] the lowest α-diversity and abundance were detected at site H, followed by the G site. The control/inactive sites were characterized by soft-bottom community values comparable with other investigations in similar habitats (i.e., [51][52][53][54][55]). In our study, PERMANOVA highlighted significant differences between the hydrothermal vents and the CN and CS sites, while PERMDISP did not show significant differences between small-scale heterogeneities (differences between same-site replicates) [14]. In classical BACI (Before and After Controls Impacts) studies [56] high small-scale heterogeneities are related to lowly disturbed areas since an impact often tends to homogenize community structures, lowering redundancy [57,58] and resilience capacity [18,59]. However, this is not always true since it may depend on what kind of species or taxa compose the assemblages. Indeed, a recent paper [60] shows evidence that good quality of coralligenous assemblages results in more homogeneous assemblages since the dominance of sessile species homogenizes community structures over a long timescale, not allowing the easy settlement of minor taxa. In the present investigation, site H showed the highest heterogeneity values although composed by the lowest number of taxa at both macro and meiobenthic level. This can be explained considering that the low number of taxa and individuals at H are represented by tolerant species, such as Tritia cuvierii (Payraudeau, 1826) and the sibling Capitella capitata (Fabricius, 1780), that can survive in such extreme conditions as those created by the presence of CO 2 [36,61]. Most macrobenthic and meiofaunal species recorded at site H are also present at site G but not at the control/inactive sites, confirming that hydrothermal vent conditions are the drivers shaping the structure of these benthic assemblages. Although meiobenthic communities are generally more structured than macrobenthic ones, site H was characterized by some dominant Nematoda genera such as Daptonema (Cobb, 1920) and Oncholaimus (Dujardin, 1845) [37] that, as previously shown, exhibit high tolerance to CO 2 [62,63].
In the present investigation some abiotic parameters, such as: pH, temperature, cation concentrations, etc., affected β-diversity, (assessed through small-scale heterogeneity measures among replicates). In particular, distLM show that temperature is significantly related to crustacean and nematodean assemblages, while acidification (measured as variation of either pH or ions concentrations) affects both macro-and meio-assemblages. Such findings are coherent with previous studies that show how the increase in temperature and the decrease in pH synergistically act on assemblages' structures [28,64]. Indeed, on one hand, low α-diversity and abundances of macrobenthic assemblages are related to the temperature rising, that eliminates low-resilient species, primarily crustacean, determining: reduction in growth and survival rates [65], and thus destructuring the community and the evenness among species. On the other hand, high heterogeneity is strictly related to the acidified environments [66,67], highlighted by low values of pH and by the presence of some cations resulting from seawater and sediments (Na + , Mg + and K + ) [68]. It is known, as previously shown by Linares et al. (2015) [32] that even a small decrease of pH may lead to dramatic shifts in highly diverse and structurally complex benthic systems. Acidification affects sensitive macrobenthic species that completely disappear, decreasing redundancy and generating available niches for other more tolerant species. This is particularly evident with Mollusca α-diversity and abundances, which are correlated with TOC at the G and H sites, as shown by the disappearance, at these sites, of Hiatella arctica (Linnaeus, 1767) (abundant in control/inactive sites) that, as also shown in other studies (i.e., [69]) cannot survive since acidification affects their shell texture [70], calcification [31] and growth rate [30].

Assemblages' Strategies to Mitigate Acidification Effects
While previous studies showed that seawater acidification causes important changes at species/taxa (α-diversity) level (i.e., [25,71]), to our knowledge this study may be the first one highlighting that β-diversity might be directly related to acidification. The protection level of this MPA excludes marked anthropic impacts, thus confirming that acidification, although increasing the heterogeneity in H and G sites, may lower β-diversity coherently with the pH gradient that, in turn, generates a loss of redundancy. The highest heterogeneity present at H and, to a lesser extent, at G is due to microhabitats available among sediment particles. Acidification might inhibit interstitial species' existence; however, they may be able to find micro-areas where they can live if the small-scale environmental conditions are not so extreme [72]. The concept of microhabitat may change when we consider macro-or meiobenthic assemblages. Macrobenthic species, due to their larger size, tend to be borrowers (endobenthic) rather than interstitial (mesobenthic), like most meiofauna; consequently, they do not occupy interstitial spaces. This causes high dispersions of macrofauna and thus great heterogeneity also for poorly structured assemblages. Conversely, meiofauna, and Nematoda in particular, find more available microhabitats in the interstitial space, due to their smaller size, giving rise to well-structured and diversified assemblages [73][74][75] and, thus, increasing synecological values and maintaining unchanged β-diversity.
The present investigation shows that more heterogeneous assemblages may be related to hydrothermal vents where the acidification effects are more evident, though reducing number of species and destructuring the communities. Thus, in these conditions high heterogeneity is representative of low β-diversity and, consequently, low resilience capacity due to the loss of redundancy at the H and G sites, since benthic organisms are forced to use only the microhabitat where extreme conditions are not present. In particular, this is more evident for macrobenthic than meiobenthic assemblages since for the former the likelihood to find survival habitats among grains of sand and mud is very scanty.

Conclusions
According to the results we obtained in this study, we can maintain that acidification affects not only benthic community structures but also their β-diversity and, in turn, their resilience capacity.
These effects, here related to a particular shallow, acidified, volcanic environment at SdF may be extended to other acidified environments due, for instance, to climate changes. Therefore, it may be useful to measure community resilience capacity of soft-bottom assemblages in order to value, at the same time, the severity of climate change effects. In this study, high heterogeneities among replicates, linked to low synecological indices, may be indicators of low resilience capacity due to environmental changes that occurred so quickly that communities could not recover to a more stable climax stage. Finally, as previously shown in other studies, high heterogeneity is not always an indicator of high β-diversity and its values should be considered in the context of the study area. Thus, in extreme marine environments, where several factors exert limiting effects on many species, the heterogeneity, considered as new opportunities for some species but also as negative conditions for other ones, may be associated with lower biodiversity values.