Spatial patterns and drivers of benthic community structure on the northern Adriatic biogenic reefs

The northern Adriatic Sea (NAS) hosts numerous biogenic subtidal reefs that are considered biodiversity hotspots. Several studies have already investigated the origin and biodiversity of these reefs. However, many of them are still unexplored and further knowledge is needed for their conservation. Here, the spatial variability, epibenthic community structure, and environmental features that characterize these habitats were investigated. Fifteen randomly selected reefs were sampled between 2013 and 2017, including some remote sites that have never been studied before. A fuzzy k-means clustering method and redundancy analysis were used to find similarities among sites in terms of epibenthic assemblages and to model relationships with abiotic variables. The results showed that these reefs are highly heterogeneous in terms of species composition and geomorphological features. The results were also consistent with previous studies and highlighted three main types of benthic assemblages defined by the dominance of different organisms, mainly reflecting the coastal-offshore gradient: nearshore reefs, generally dominated by stress-tolerant species; reefs at a middle distance from the coast, characterized by sponges, non-calcareous encrusting algae and ascidians; offshore reefs, dominated by reef builders. However, distance from the coast was not the only factor affecting species distribution, as other local factors and environmental characteristics also played a role. This kind of biogenic reefs in temperate seas are still poorly known. The present work contributed to shed further light on these habitats, by complementing the results of previous studies on their natural diversity, highlighting the specificity of the epibenthic communities of NAS reefs and the need to improve current, still inadequate, conservation measures.


Introduction
Marine organisms generally exhibit complex spatial distribution patterns, often related to environmental gradients. Habitat types, species distribution, and composition are the result of various elements that include biotic and abiotic factors as well as historical events Brown et al. 2016;Silva et al. 2021). For example, large and rapid fluctuations in abiotic variables such as salinity and temperature can be critical for both vagile organisms, that move across spatial gradients, and sessile or sedentary species, that cannot escape such changes (Laprise and Dodson 1994;Pansch et al. 2018;Jones et al. 2020). These fluctuations can cause significant physiological stress to organisms, affecting their metabolism, reproduction, and ultimately their survival (Solan and Whiteley 2016;Pourmozaffar et al. 2020). Similar effects can result from increased sediment runoff, which reduces photosynthesis, recruitment, and growth rates, and increases burial or scouring (Airoldi 2000;Balata et al. 2005;Irving et al. 2009;Erftemeijer et al. 2012). Impacts can also be exacerbated by continued nutrient enrichment, which can lead to oxygen depletion and water quality degradation (Piazzi et al. 2012;Quimpo et al. 2020). These factors can then interact with other important biological and physical processes, such as larval and propagule supply, species competition, substrate types, and local oceanographic conditions (Rattray et al. 2016;Bandelj et al. 2020;Lurgi et al. 2020), resulting in complex responses that are often difficult to disentangle (Strain et al. 2014;Boyd and Brown 2015;Carrier-Belleau et al. 2021). In the long term, several studies have shown that this environmental variability plays an important role in regulating the distribution, abundance, and species composition of communities and shaping spatial patterns (Airoldi 1998;Mangialajo et al. 2008;Darling et al. 2013;Strain et al. 2014;Yang et al. 2015;Montefalcone et al. 2017). Therefore, the analysis of factors affecting species distribution and community composition is crucial for understanding environmental patterns and processes, adequately assessing ecological status, and for effective conservation measures (Crain and Bertness 2006;Pinedo and Ballesteros 2019).
Although the presence of some biogenic subtidal reefs on the incoherent substrates of the northern Adriatic Sea (NAS) has been known for at least two centuries (Olivi 1792), it is only in recent decades that the existence of a widespread mosaic of reefs has been confirmed (Stefanon and Zuppi 2000;Gordini et al. 2012). Locally referred to as "trezze" or "tegnùe," they are particularly important because they form hard substrates in an extensive soft-bottom area, suitable for the settlement of sessile species (Fava et al. 2016). The rich and diverse epibenthic communities they host are attractive to fishermen and divers and contribute to increasing the socio-economic value of the area (Tribot et al. 2016;Chimienti et al. 2017;Tonin 2018). However, several local stressors threaten their health and conservation: chemical pollution, changes in salinity and nutrient concentrations due to river runoff, direct physical damage from fishing, anchoring, or diving, reduction in water clarity due to suspended sediments, hypoxia/anoxia events, marine debris, smothering due to high sediment loads and mucilaginous aggregates, and potential oil spill due to the high marine traffic in the area (Pranovi et al. 2000;Melli et al. 2017;Moschino et al. 2019;Farella et al. 2021;Ponti et al. 2021). Climate change impacts leading to heat waves (Galli et al. 2019;Garrabou et al. 2022, and references therein), increased frequency and intensity of storms (Teixidó et al. 2013), and seawater acidification (Zunino et al. 2019(Zunino et al. , 2021 are additional sources of disturbance to these habitats. For these reasons, biogenic reefs are internationally protected (UNEP-MAP-RAC/SPA 2008), and some of them have been included in Natura 2000 sites.
The primary origin of these reefs can be linked to various processes of diagenesis (Stefanon 1970), methanogenesis (Gordini et al. 2012;Donda et al. 2015), and lithification of sandy fossil channel-levee systems (Tosi et al. 2017). Calcareous benthic organisms further shaped these reefs by forming calcareous bio-concretions with different geomorphologies on the primary hard substrate (Gordini et al. 2012;Turicchia et al. 2022). Benthic assemblages on these reefs differ markedly from those on the coralligenous banks elsewhere in the Mediterranean (Ingrosso et al. 2018), exhibiting greater spatial heterogeneity (at local and regional scales) than temporal variability (between years) in species composition and abundance (Ponti et al. 2011;Curiel et al. 2012;Falace et al. 2015). These differences are likely due to a combination of factors, including the extent of reef isolation, sedimentation, nutrient loading, and the large spatial variability of environmental properties of the region, which is largely influenced by local geology, hydrology, climate, and human use. Indeed, as a semi-enclosed basin (surrounded by extensive drainage areas) with a maximum depth of less than 70 m and an average depth of 30 m, the NAS is particularly influenced by freshwater inputs from rivers and lagoons, as well as by strong hydrological and meteorological variability (Aubry et al. 2004;Solidoro et al. 2009). Surface currents, which normally move counterclockwise along the coast (Poulain et al. 2001), change direction during strong winds (Bora or Burjan, the prevailing NE wind in the region) (Kuzmić et al. 2006) and during freshwater plume intrusion from the Po River (Jeffries and Lee 2007). The geomorphological characteristics of reefs (shape, height above the seafloor, depth, degree of isolation) also affect the structure of benthic communities. For example, species richness and diversity have been found to increase with depth and decrease with reef height and seafloor mud content (Ponti et al. 2011). Hydrodynamic conditions, and coastal freshwater, also play a critical role in regulating beta diversity and species composition, with species sensitive to sedimentation and high nutrient loading being displaced, especially on offshore reefs (Ponti et al. 2011;Falace et al. 2015;Bandelj et al. 2020). Benthic communities on biogenic reefs appear to be more interconnected than with the coast, depending on the pelagic larval duration (PLD) of species and small-scale variability of currents (Bandelj et al. 2020). However, some authors have shown that asexual reproduction is also a widespread strategy to colonize space and maintain species diversity in these communities (Fava et al. 2016).
Although many studies have documented various aspects of these biogenic reefs, more information is needed to better understand the role of geomorphology and the effects of environmental gradients on benthic community structure, especially given the large number of reefs that remain unexplored. In addition, most of the studies have considered the zoobenthic and phytobenthic communities separately, presented qualitative data, and have often been published in local journals. The objective of this study was to (1) investigate the spatial variability of both zoobenthic and phytobenthic assemblages; (2) increase knowledge of the structure of epibenthic assemblages; (3) unravel the key environmental variables and geomorphological features that influence the observed variability; (4) examine unstudied reefs to assess whether the observed dynamics and patterns are consistent with other reefs. The results of this study will be particularly useful for improving conservation efforts in the NAS, as current regulations appear insufficient to ensure preservation of species connectivity and benthic community representativeness (Bandelj et al. 2020).

Study area
The NAS, the northernmost part of the central Mediterranean ( Fig. 1), is characterized by sedimentary coasts on the northern and southwestern Italian sides of the basin, with numerous coastal lagoons (Grado-Marano, Caorle, and Venice) and estuaries (Isonzo-Soča, Tagliamento, Piave, Adige), including the large delta of the Po river. The eastern and northeastern coasts, on the other hand, are predominantly rocky with high and steep shores, bays, and deep fjords. The seafloor of the NAS consists of a large continental shelf characterized by distinct terrigenous sedimentary facies (Spagnoli et al. 2014). Several hundred biogenic reefs are scattered randomly on sandy and muddy bottoms in this area, at depths ranging from 7 to 30 m and 1 to 40 km from the coast. They are separated by distances ranging from a few meters to several tens of kilometers (Trincardi et al. 1994;Caressa et al. 2001;Gordini et al. 2012), have either large horizontal surfaces or high structures (up to 4 m high above the seafloor), and are either widespread or clustered ranging in size from 1 m 2 to 1 km 2 (Casellato and Stefanon 2008;Gordini et al. 2012). Most of them are in areas almost parallel to the coast between Chioggia and Grado, off the Venice and Grado-Marano lagoons, while others are more isolated and offshore ( Fig. 1).

Data collection
Fifteen biogenic reefs located between the Gulf of Trieste and the Gulf of Venice, including two remote previously unstudied reefs (West Bike and East Bike), were randomly selected considering different distances from the coast, depths, and morphologies ( Fig. 1; Table 1). Ten random 588 cm 2 areas (21 × 28 cm rectangles) were photographed on each reef in July-August 2013 or July-August 2017 ( Fig. 1; Table 1). Either a Canon PowerShot G12 or a Canon PowerShot G15 digital camera (with 10 and 12 megapixels, respectively; Canon, Tokyo, Japan) equipped with an aluminum underwater case, an S-TTL strobe (Inon D-2000), and a custom-made stainless-steel frame was used. Percent cover of sessile organisms, sediment, and mucilage was quantified by overlaying a grid of 400 equally sized squares using photoQuad image processing software (Trygonis and Sini 2012). Percent cover was based on the total readable area of each image, determined by subtracting dark and blurred zones or areas covered by vagile organisms (Ponti et al. 2018). Abundance of boring organisms was quantified based on their surface detection. Organisms were identified at the lowest possible taxonomic level (Table S1). Taxa were also assigned to different epibenthic categories: "Reef builders," "Bioeroders," "Non-calcareous (NC) encrusting algae," "Turf," "Sponges," "Cnidarians," "Ascidians," and "Other algae" (Table S1).
Geomorphological information on the reefs (i.e., height above the seafloor (m), depth (m), distance from the coast (km), and class of extent ("Very small": < 100 m 2 , "Small": 100-2500 m 2 , "Medium": 2500-10,000 m 2 , and "Large": 10,000-50,000 m 2 ) was obtained from available geophysical surveys conducted with a series of single and multibeam echosounders and side-scan sonars (Fortibuoni et al. 2020a, b;Gordini and Ciriaco 2020;Ponti 2020a, b), and visually confirmed by scientific divers during sampling. Sediment grain size reported as phi-scale (φ), was obtained from a detailed sedimentological map produced by interpolating the raw data (3-metric IDW) with a grid resolution of one kilometer. The inter-  (Bostock et al. 2018). In situ measurements of physicochemical variables at the sea surface and bottom of the sampling sites (annual median, fifth, and ninety-fifth percentiles of chlorophyll (µg/L, Med/P5/P95_CHL), seawater temperature (°C, Med/P5/P95_T), salinity (psu, Med/P5/P95_SAL), dissolved oxygen (mL/L, Med/P5/P95_OXY), ammonium (µmol/L, Med/P5/P95_NH4), nitrates (µmol/L, Med/P5/P95_NO3), phosphates (µmol/L, Med/P5/ P95_PO4), were obtained from the dataset reported in Solidoro et al. (2009), and from water sampling of the Regional Water Authorities (ARPA) of Veneto and Friuli Venezia Giulia Regions, using the same procedure described in Falace et al. (2015). The fifth and ninetyfifth percentiles of the physicochemical variables were considered instead of the absolute minimum and maximum values to avoid possible biases. Hydrodynamic data were obtained from a high-resolution numerical simulation of the NAS. In particular, the annual mean and maximum velocities (m/s, Mean_V, Max_V) and the mean kinetic energy (per unit mass) of the currents (m 2 /s 2 , Mean_KE) at the sea surface and at the seafloor were analyzed. The numerical simulation is based on the Massachusetts Institute of Technology general circulation model (MITgcm), a three-dimensional finite volume model for geophysical fluid dynamics (Marshall et al. 1997).

Data analysis
A Fuzzy k-Means (FkM) method was used to identify potential groups of sites with similar species composition and percent cover ). This non-hierarchical clustering method assigns a value to all objects based on the strength of their membership in different groups. In this case, a site with a strong membership value to a group is associated with that group; otherwise, it has a high degree of fuzziness and a higher degree of association with other groups (Bezdek 2013). Prior to clustering, a Hellinger transformation was applied to the data to overcome the "double zero" problem (Legendre and Gallagher 2001). Analyses were computed by increasing the fuzziness parameter from 1.1 to 2 in increments of 0.1 and setting the number of random initializations to 999 (Borcard et al. 2011). A fuzziness parameter of 1.2 was chosen because lower values produced almost a non-fuzzy clustering, while higher values were unable to unambiguously assign all sites to a group. FkM clustering was performed using the "cluster" package for R (Maechler et al. 2022). The best number of clusters was selected by calculating various cluster validity indices (i.e., partition coefficient, partition entropy, modified partition coefficient, silhouette, fuzzy silhouette, Xie, and Beni) using the "fclust" package in R (Ferraro et al. 2019).
For each site and group identified by FkM, we calculated: species richness (S, number of taxa); Hill's diversity index (N1), which can be interpreted as the effective number of species (i.e. the number of species in an equivalent community consisting of equally abundant species) derived from Shannon entropy (Cao and Hawkins 2019); the corresponding Hill's evenness index (N10); the percent cover of sessile organisms, sediment, and mucilage; and the percent cover and number of taxa of each epibenthic category. The height above the seafloor, depth, distance from the coast, class of extent, and sediment grain size for each group were also assessed using geomorphological information from the sites. In the text and graphs below, the mean values are given along with their standard errors (SE). To evaluate possible monotonic relationships, Spearman's rank correlation coefficients and corresponding p-values were calculated between geomorphological variables, mean values of diversity indices, and percent cover of epibenthic categories per site (α = 0.05, n = 15) using the "Hmisc" package for R (Harrell 2019). A two-way permutational analysis of variance (PERMANOVA; α = 0.05, (Anderson and Robinson 2001) based on Euclidean distances was applied to assess differences in mean values of geomorphological attributes, diversity indices, and epibenthic categories between groups (Gr, fixed) and among sites (Si) asymmetrically nested in Gr, as resulting from cluster analysis. P-values were obtained by permutation of residuals under a reduced model: when less than 1000 unique permutations were available, the asymptotic Monte Carlo p-value was used instead of the permutational one. Significant values for the Gr factor were examined by post-hoc pairwise tests. Analysis was performed using PRIMER 7 with the PERMANOVA + add-on package (Anderson et al. 2008).
The IndVal index (Dufrêne and Legendre 1997;De Cáceres and Legendre 2009) was applied to the discretized FkM results to identify characteristic taxa within site groups. This method combines the mean percent cover and frequency of occurrence of all taxa and determines their indicator values within all group combinations. Only taxa with a significant IndVal value (α = 0.05) after 999 permutations of samples were considered indicators. Maximum IndVal values of indicator taxa were then reported on all divisions of a dendrogram built on the results of the FkM to show taxa representative of each group with eurytopic and stenotopic characteristics. Eurytopic taxa have a large niche width, and their indicator value decreases when the sites, where they are abundant, are divided into several groups, while stenotopic taxa are indicators of only one group of sites, suggesting that they have a small niche width (Dufrêne and Legendre 1997). The IndVal was calculated using the "indicspecies" package for R (De Cáceres and Legendre 2009).
Finally, multivariate direct gradient analysis using redundancy analysis (RDA) (van den Wollenberg 1977) was applied to the FkM results to model the relationship with geomorphological, physicochemical and spatial descriptors . Spatial descriptors used in the present study were longitude and latitude, and the Moran Eigenvector Maps (MEMs) (Dray et al. 2006). The MEMs are a method to explicitly modelling the spatial signal by decomposing it into a set of orthogonal, i.e., mutually independent, spatial components (Borcard and Legendre 2002). The components thus represent the whole spatial signal, from the smallest scale (local autocorrelation) to the largest scale, encompassing the whole sampling area, and can be used as explanatory variables (Bauman et al. 2018a) along with environmental variables to understand the organization of benthic assemblages (Corte et al. 2018). While MEMs can represent any spatial signal, linear gradients can be modelled more efficiently using longitude and latitude (Borcard et al. 2011). Thus, the first step before deriving the MEMs consisted in checking with an RDA the existence of a significant linear gradient and successively detrending the dataset, i.e., removing the linear gradient as modelled by longitude and latitude. Based on this procedure, the MEMs we derived represented additional independent spatial components, not correlated with longitude and latitude. The MEMs were then derived from metric coordinates (WGS84 / UTM33N) by diagonalizing a doubly centered spatial weighting matrix (Bauman et al. 2018b). The spatial matrix is obtained by the product of a connectivity matrix and a weighting matrix. Different graphbased connection schemes (Delaunay triangulation, Gabriel's graph, relative neighborhood graph, and minimum spanning tree) were tested to build the connectivity matrix (Bauman et al. 2018b). Specific weighting functions were also tested to build the weighted matrix (the linear function, the concave-down function (with values from 1 to 10), and the concave-up function (with values from 0.1 to 1) (Bauman et al. 2018a, b)). Overall, 84 models (groups of MEMs) with a maximum of seven MEMs per group, were obtained considering different types of connection schemes and weighted functions. Only positively autocorrelated eigenvectors were retained for each model tested (Borcard et al. 2011). An RDA was run for every model on the residuals of the spatial model of longitude and latitude, and a forward selection with a double-stop criterion (the adjusted r 2 value of the model and a significant p-value, α = 0.05) (Blanchet et al. 2008) was performed to select the significant MEMs in each model. The model and the corresponding MEMs with the highest adjusted r 2 value were retained. A separate RDA was also performed for both the geomorphological and physicochemical datasets; the hydrodynamic variables were included in the physicochemical dataset (see results). Axes and models were tested for significance, and significant variables were selected by forward selection using a double-stop criterion (Blanchet et al. 2008). Highly correlated variables were removed after being tested for multicollinearity using Spearman's rank correlation coefficients (Borcard et al. 2011). The selected variables from the geomorphological, physicochemical and spatial datasets were then used as explanatory variables for the final RDA model. Variation partitioning (Borcard et al. 1992) was applied to the three groups of variables in the final RDA model to examine their mutual relationships and contribution to explained variance in the final model. Analyses were performed using the "vegan" (Oksanen et al. 2016), "adespatial" (Dray et al. 2023) and "packfor" (Dray et al. 2016) packages for R.
Cluster analysis showed that the reefs could be classified into four groups according to taxa composition and percent cover (Fig. 1, S2a- Table 1, and Table S2). All cluster validity indices were consistent with this result. The first group includes five sites off Grado: Biro, Gubana, La Bomba, Pescecane, Tettoia (Fig. S2a). The second group includes three sites: one closer to the first group (Colomba) and two off Chioggia and the Venice lagoon (TQS, MR08; Fig. S2a). The third group consists of five sites near Chioggia (P204, TM1, AL06, P208B, P213; Fig. S2b), and the fourth group includes two sites about 38 km from the coast (West Bike and East Bike; Fig. S2b). Most sites had a high membership value (> 0.90) and were clearly assigned to one group, whereas only P213 had a higher degree of fuzziness and showed some degree of membership to two other groups (0.12 for Group 1 and 0.16 for Group 2) (Table S2).
Group 4 was significantly less species-rich (S) and diversified (N1) than the other groups; Group 1 was more diverse (N1) than Group 2 (Fig. S3, Table S3). Although there were significant differences among sites, the groups did not differ from each other in terms of species evenness (N10). The differences between sites and groups were particularly evident when comparing epibenthic categories. Groups 3 and 4 were more similar than Groups 1 and 2 in terms of mean percent cover of epibenthic categories, especially for "Turf" (Fig. 2, Fig. S1). Mean percent cover of "Bioeroders" was significantly higher in Groups 1 and 2 than in the other groups (p < 0.05), while "Reef builders" predominated in Group 2 with a mean percent cover (67.9 ± 7.1%) that was statistically higher than that of the other groups (p < 0.01), which again did not differ from each other ( Table S4). Mean percent cover of "Ascidians" varied from 15.9 ± 4.2% in Group 1, 5.6 ± 1.5% in Group 2, and almost absent in Group 3 (1.4 ± 0.8%) and Group 4 (0.05 ± 0.05%) (Fig. 2); however, due to the high variability among sites within groups, statistically significant differences were found only between Groups 1 and 3 and between Groups 2 and 3 (Table S4). Cover of "Cnidarians," "Other algae," "Sponges," "NC encrusting algae," and "Turf" showed significant variability among sites within groups, but did not differ statistically among groups (Table S4). Except for "Cnidarians," the mean number of taxa by epibenthic category varied greatly among sites within groups (p < 0.01, Fig. S4 , Table S5); nevertheless, some trends are apparent (Figs. S4, S5, Table S5). In particular, "Reef builders" and "Bioeroders" were generally significantly more diverse in Group 2 than in the other groups. "Ascidians" were more diverse in Group 1 and 2 than in Groups 3 and 4. Similarly, "NC encrusting algae" (not found in Group 4) were more diverse in Group 1 and 2 than in Group 3. Finally, "Sponges" seems to be more diverse in Group 1, but no statistically significant differences can be detected, likely due to the high variability within groups.
A total of 20 indicator taxa were identified using the IndVal index (Fig. 3). The maximum IndVal values of the indicator taxa were high (> 0.70) and statistically significant in all combinations of the groups. The index showed that Groups 1 and 2 shared some eurytopic taxa, such as the sponges C. viridis and Aplysina cavernicola, the ascidian Aplidium tabarquensis, and the red alga Peyssonnelia rosa-marina. Two reef-builders, namely Lithophyllum incrustans and vermetids, a boring sponge (C. rhodensis) and an ascidian (Cystodytes del- lechiajei) distinguished Group 2 from Group 1; conversely, Group 1 was characterized by some sponges and ascidians. Groups 3 and 4 were both characterized by an encrusting sponge: Dictyonella incisa. However, Group 3 was characterized mainly by foliose algae such as Rhodophyllis sp. and Nitophyllum punctatum, by a turf-forming red alga (Antithamnion sp.) and by a cnidarian (Sarcodictyon catenatum), while Group 4 was characterized by two sponges (Polymastia mamillaris and Ulosa digitata) and by mucilage (Fig. 3).

Geomorphological and sedimentary features
The reefs studied also differed in terms of geomorphological features (Table 1; Fig. 4). Group 4 sites were significantly deeper (30 ± 0.3 m) (p < 0.01) and farther from the coast (37.9 ± 0.5 km) (p < 0.05) than those of the other groups (Table S6). On average, Groups 1 and 2 were also significantly farther from the coast (p < 0.05) than Group 3, which includes sites near the town of Chioggia (8.7 ± 1.5 km). The sediments on the seabed were finer near the sites of Group 3, consisting mainly of mud and sandy silt, but they were only significantly different from those of Group 1 (p < 0.01). No significant differences were found between groups in terms of average height above the seafloor (Table S6). In Groups 1, 2 and 3, the class of site extent was more heterogeneous than in Group 4, which was represented by only two reefs of medium extent. Group 1 reefs were predominantly small (60%). Group 2 reefs ranged from very small to medium, and Group 3 reefs ranged from small to large (Fig. 4). Mean percent sediment cover of reefs was significantly higher in Groups 3 and 4 than in Groups 1 and 2 (Fig. S6a, b, Table S6), although it differed significantly between sites (p < 0.001). High percent cover of mucilage was found only in West Bike (41.3 ± 3.8%) and East Bike (23.7 ± 3.1%) (Fig. S6c).

Correlations between epibenthic assemblages and abiotic features
Univariate rank correlations showed some monotonic relationships between mean characteristics of epibenthic assemblages and abiotic characteristics of the studied reefs. Species richness correlated positively with Hill's diversity index (r s = 0.78, p < 0.001) and negatively with depth and percent sediment cover (r s < -0.6, p < 0.01) (Fig. S7). Depth was also negatively correlated with Hill's diversity index (r s = -0.65, p < 0.01). Significant positive correlations were found between depth and distance from the coast (r s = 0.61, p < 0.05) and between percent sediment cover and sediment grain size on the phi scale (r s = 0.54, p < 0.05), indicating higher sediment deposition on reefs with finer sediments. However, percent sediment cover and sediment grain size were not correlated with depth or distance from the coast, suggesting different fluvial influences and sedimentation regimes in the area. Mean percent cover of "Turf" and "Other algae" were positively correlated with each other (r s = 0.85, p < 0.001), and the latter was also negatively correlated with distance from the coast (r s = -0.59, p < 0.05) (Fig. S8). Sediment grain size was negatively correlated with mean percent cover of "Ascidians" (r s =-0.66, p < 0.01), "Sponges" (r s = -0.55, p < 0.05), and "Bioeroders" (r s = -0.54, p < 0.05), while there is a positive correlation with that of "Turf" (r s = 0.52, p < 0.05), "Cnidarians" (r s = 0.57, p < 0.05) and "Other algae" (r s = 0.52, p < 0.05). Mean percent cover of "Ascidians," "Bioeroders," and "NC encrusting algae" were also negatively correlated with mean percent cover of sediment and positively correlated with each other (Fig. S8). Height above the seafloor was negatively correlated with mean percent cover of "Ascidians" (r s = -0.55, p < 0.05) (Fig. S8).

Effect of spatial and environmental descriptors on the epibenthic assemblage structure
The spatial analysis revealed a significant linear gradient of longitude and latitude (p = 0.001) explaining 54% of variance; thus, the dataset was detrended before deriving MEMs. After calculation of the adjusted r 2 and forward selection of the MEMs, only large-scale MEM2 and small-scale MEM5 of the model built with Gabriel's graph connection scheme and the concave-up function with a coefficient value equal to 0.1 were retained (Fig. S9). The best model had an adjusted r 2 of 0.62. MEM2 explained 29% of the variance (p = 0.025) and MEM5 33% of variance (p = 0.003). The same model explained 71% of the variance when it was run on the non-detrended dataset. MEMs consistently modelled a quite clear spatial pattern with sites closer to the coast of Chioggia and Grado having more similar values at large-scale (MEM2), while MEM5 mostly discriminated sites at a local scale in the southwestern part of the study area (i.e., sites TQS and MR08 having high values, and site P208B, P213, West Bike, and East Bike having low values on the MEM) (Fig. S9).
The RDA geomorphological model was run with eight variables (depth, distance from the coast, reef height, sediment grain size, class of extent (dummified) "Very Small," "Small," "Medium," and "Large"). The model had an adjusted r 2 of 0.55, with only the first axis being statistically significant (p = 0.003) and explaining 30% of the variance. After removing the most highly correlated variables and performing forward selection, only sediment grain size and distance from the coast were retained, and the adjusted r 2 was reduced to 0.41. The first and second axes explained 25% and 16% of the variance, respectively, and were both statistically significant (p = 0.001 and 0.01) with 999 permutations.
Because all hydrodynamic variables (annual Mean_V, Max_V, Mean_KE at the sea surface and at the seafloor) were highly correlated with each other (r s > 0.9), RDA was not performed for this subset of variables and only Mean_V at the sea surface was retained, because it was less correlated with the other variables, and included in the physicochemical dataset. The RDA model for the physicochemical dataset was built with forty-three variables (annual Med, P5, and P95 of CHL, T, SAL, NH4, NO3, PO4, OXY at the sea surface and bottom, and annual Mean_V at the surface). The model had an adjusted r 2 of 0.99; 44% of the variance was explained by the first axis (p = 0.02), while the second axis explained 34% of the variance after 999 permutations (p = 0.03). After reducing multicollinearity, the following variables were retained: Mean_V and P5_T on the surface and P5_CHL, Med_ NH4, and P95_PO4 on the bottom. The model had an adjusted r 2 of 0.65; both the first and second axes were statistically significant (p = 0.001 and 0.01) and explained 31% and 20% of the variance, respectively.
Variation partitioning performed on the final RDA model run with the variables selected in the three RDA subsets, showed that a high fraction of the variance was jointly explained by the three datasets (Table S7). The spatial and physicochemical datasets shared the highest joint contribution (r 2 adj.= 0.48). The unique contribution of each single dataset to the explained variance of the model was negligible (Table S7). According to these results, since the information brought by the three datasets was highly overlapped, we decided to run the final RDA model only with the geomorphological and physicochemical variables that shared the lowest fraction of variance (r 2 adj.= 0.09) and enables us to discuss the observed biotic patterns in terms of environmental factors. The adjusted r 2 of the final RDA model was 0.76; the first axis explained 37% (p = 0.001), whereas the second axis explained 23% of the variance after 999 permutations (p = 0.01) (Fig. 5). The RDA plot showed that P5_CHL at the bottom and P5_T at the surface tended to increase toward Group 1 sites; geomorphological variables such as Mean_V at the surface and sediment grain size, as well as P95_PO4 and Med_NH4 at the bottom, tended to increase toward Group 3 sites. Distance from the coast and, to a lesser extent, P5_T at the surface tend to increase toward Group 2 and 4 sites, while P95_PO4 and Med_NH4 at the bottom tend to decrease (Fig. 5).

Discussion
The present study has shed further light on the biogenic reefs of the NAS, by partially confirming previous findings and clarifying some ecological aspects related to the distribution of epibenthic assemblages. For the first time, reefs at about 38 km offshore were considered; very little information was previously available on these bioconstructions (Nicoletti et al. 2015).
The results showed that the biogenic reefs in the NAS are very heterogeneous in terms of epibenthic assemblages and show a strong correlation with some environmental variables (i.e., distance from the coast, nutrients, hydrodynamic conditions, sediment cover, and sediment grain size). The structure of the studied communities is consistent with previous reports and confirms the overall dominance of filter-feeding organisms (e.g., sponges, ascidians), turf-forming algae, and encrusting calcareous algae with a distinct distribution pattern (Molin et al. 2003Casellato and Stefanon 2008;Curiel and Molin 2010;Curiel et al. , 2012Curiel et al. , 2014Curiel et al. , 2017Ponti et al. 2011;Miotti et al. 2014;Falace et al. 2015;Fava et al. 2016;Nesto et al. 2020). The results also confirmed striking differences with assemblages growing on typical Mediterranean coralligenous reefs (Ballesteros 2006;Ponti et al. 2011;Falace et al. 2015;Ingrosso et al. 2018), as representative taxa such as hydrozoans, gorgonians, the green algae Flabellia petiolata and Halimeda tuna, and erect bryozoans were not recorded in our sampling, while encrusting bryozoans had a comparatively low percent cover (< 0.33 ± 0.24%). With few exceptions, Lithophyllum incrustans and Lithothamnion sp. are the most abundant calcareous algae on NAS biogenic reefs, while Lithophyllum stictiforme and Mesophyllum spp. are generally the predominant taxa on typical Mediterranean coralligenous reefs (Ballesteros 2006). The biogenic reefs of the NAS also show some differences with those of the southern Adriatic continental shelf in terms of the composition of the main taxa (Bracchi et al. 2017;Piazzi et al. 2019). In this respect, the present work confirms the specificity of the epibenthic assemblages inhabiting the reefs of the NAS (Ponti et al. 2011;Falace et al. 2015).
Cluster analysis provided clear evidence of the wide variability of reefs. Four distinct groups of sites were identified, and only one site (P213) showed some degree of membership with the other groups. The distribution pattern of the benthic communities seems to be mainly determined by specific environmental characteristics of the area, as shown by the RDA model results. Sites near Chioggia (Group 3) are characterized by more stressful environmental conditions. Here, the main influence of the rivers and the Venice Lagoon  Table 1 results in higher nutrient loading, stronger water flow, higher suspended sediment load and limited water transparency, affecting the biotic community structure and allowing the dominance of turf-forming algae. This is also confirmed by higher percent cover of sediment and stress-tolerant species associated with turf-forming algae, such as the encrusting sponge Dictyonella incisa and foliose erect algae (e.g., Rhodymenia sp.), as observed by Curiel et al. (2014). Conversely, sites relatively far from the coast (Group 2) with less energetic hydrodynamic conditions, lower nutrient loading, and less sedimentation appear to be characterized by more stable conditions, allowing the dominance of encrusting calcareous red algae. Indeed, several studies have shown that nutrient enrichment and sedimentation strongly alter coralligenous assemblages by reducing bio-builder taxa and increasing ephemeral algae (Piazzi et al. 2011(Piazzi et al. , 2012Montefalcone et al. 2017). Colonial ascidians such as Polycitor adriaticus, Aplidium conicum, and sponges such as Tedania (Tedania) anhelans, which are sensitive to fine sediments that may occlude their filtration systems (Naranjo et al. 1996;Molin et al. 2003;Turicchia et al. 2021), are also well represented at Group 2 sites. These results are consistent with previous observations by other authors who reported, for example, the presence of P. adriaticus and A. conicum on reefs with little or no hydrodynamic turbulence and particles in the water column and at greater distances from the coast (Ponti and Mastrototaro 2006;Ponti et al. 2011;Miotti et al. 2014;Falace et al. 2015). Interestingly, bioeroders are particularly abundant at Group 2 sites, confirming their association with reef builders and the presence of ongoing bioconstruction and bioerosion processes that are critical for maintaining reef habitats and their diversity (Turicchia et al. 2022). Group 1 reefs are instead dominated by sponges and ascidians, with a high percent cover of "Turf" (17.5 ± 6.7%) and can be considered as an intermediate state in terms of epibenthic assemblages between nearshore (Group 3) and relatively offshore reefs (Group 2). Such pattern was also confirmed by the spatial analysis (MEM) that showed a significant linear gradient modelled by longitude and latitude, explaining a high portion of variance and suggesting a strong influence of distance from the coast and river inputs on the presence and abundance of epibenthic assemblages. The high collinearity of the spatial components with the physicochemical variables gives further confirmation that environmental variables are strongly dependent on the spatial structure and shape the epibenthic assemblages on a coastal-offshore gradient. The residual spatial pattern was described by only two significant MEMs, probably due to the simple sampling design with a limited number of sites.
However, although our results confirm the habitat typologies and the role of the coastaloffshore gradient in shaping epibenthic assemblages, as previously shown by Falace et al. (2015), distance from the coast did not contribute to explain fully the observed species distribution patterns. In particular, compared to the results of Falace et al. (2015), a fourth group was detected by the clustering analysis comprising the two new study reefs, located farthest from the coast (West Bike and East Bike), and that have peculiar conditions. Based on the predictive model developed by Falace et al. (2015), it was expected that these two reefs would most likely host epibenthic assemblages dominated by reef builders. Conversely, the present study has shown that they do not have the same epibenthic assemblages as Group 2 sites, but are more similar to Group 3 sites, dominated by turf-forming algae and some sponges, while reef builders, ascidians, and NC encrusting algae are rare. The species richness and Hill's diversity index of Group 4 are significantly lower than those of the other groups. This condition can most likely be explained by the degree of isolation of West Bike and East Bike, which may limit larval recruitment. In addition, the high percent sediment cover measured on these reefs suggests that other local factors, such as human activities and possible seasonal variation in local water circulation, may limit the settlement and recruitment of benthic species, increasing sedimentation and by reducing light penetration at very deep sites (~ 30 m). Remarkably, West Bike was the only reef to host the sponge Polymastia boletiformis, which forms a large and dense population there (foreground in Fig. S2b). The presence of P. boletiformis, generally considered to be a North Atlantic species (Plotkin et al. 2011), represents an important finding, as there have been no records in the northern Adriatic since the first and only report (as P. robusta) in the late nineteenth century by von Lendenfeld (1896).
Illegal fishing, vessel traffic, diving, commercial fishing near reefs, and marine debris are among the main human activities that can affect benthic communities in the area (Falace et al. 2015;Melli et al. 2017;Tonin 2018;Moschino et al. 2019;Shabtay et al. 2019). For example, high densities of fishing nets covering benthic organisms on various reefs, carried by sea currents from nearby shellfish farms or lost/deposited by fishermen have been reported (Melli et al. 2017;Moschino et al. 2019). The occurrence of certain species may be thus the result of either environmental changes and/or direct human activities. However, partitioning the effects of these two factors is challenging, because information on small-scale environmental pressures and human impacts on species is not available. In addition, various biological and ecological processes (e.g., reproductive timing, larval and dispersal processes, settlement and recruitment processes, intra-and interspecific interactions) in combination with other processes could play an important role in maintaining species diversity at different spatial scales (Carson et al. 2010;Dawson et al. 2014;Fava et al. 2016;Costantini et al. 2018). Local hydrodynamic conditions, affecting both larval dispersal and their ability to find suitable substrates, may account for the large variability in reproductive success or mortality associated with pre-and post-settlement processes. Bandelj et al. (2020) demonstrated, for the study area, that species with different pelagic larval duration (PLD) follow different dispersal dynamics, and that passive dispersal by currents has a greater relative importance than geographic distance in predicting the beta diversity of biogenic reefs. The results of this study also offer some possible explanations for some of the observed patterns. For example, despite similar epibenthic assemblages, the Group 2 sites appear to be divided into two widely separated reef subgroups located in the northern and southern parts of the study area. Bandelj et al. (2020) showed that these sites are well connected at PLDs of 1 day and 3 days thanks to water circulation and the presence of other reefs in the middle, which were not considered in the present study. Therefore, different aspects and local factors need to be taken into account when developing management plans and conservation measures, as well as the monitoring of the characteristics and ecological status of biogenic reefs in the NAS.
Currently, there are only three Natura 2000 (N2K) sites in the study area: the "Trezze San Pietro and Bardelli" (Site code: IT3330009) in the north, the "Tegnùe di Chioggia" (Site code: IT3250047) in the south and the "Tegnùe di Porto Falconera" (Site code: IT3250048) in front of the town of Caorle, with a total area of about 5600 ha. These areas contain reefs that are strategic for maintaining cross-scale connectivity of biogenic reefs, but many other sites that harbor very rare species (e.g., P. boletiformis in West Bike) or act as source for species dispersal, according to Bandelj et al. 2020, are currently not protected. Moreover, recent projects (Interreg Italy-Croatia project ECOSS, and Interreg Italy-Croatia project CASCADE) reported the absence or inadequacy of management plans and suggested some monitoring activities to be carried out on a regular basis (Gianni et al. 2022). The existing protection measures restrict some activities (e.g., trawling, anchoring), but do not reduce the widespread anthropogenic pressures, such as high nutrient load, sedimentation, and old ghost nets; the latter affecting reefs in the Tegnùe di Chioggia N2K site (Melli et al. 2017;Moschino et al. 2019). Since physicochemical variables have a stronger influence on reef biodiversity, as shown in this study, a change in water quality due to increased nutrients and chemical pollution from rivers and coastal discharges, may result in changes of the benthic assemblages, reducing the presence of highly structured communities even on offshore reefs, with higher biodiversity and stress-sensitive species. Policies to conserve these habitats must, therefore, not disregard the proper management of human activities on land and the correct treatment of wastewater that reaches the sea. Moreover, specific measures to enhance NAS reefs conservation should focus on enforcing the activity regulations and improving surveillance of these areas supported by funding from national and local governments. An ecosystem-based approach that considers the development of joint management plans between several N2K sites, as conceived by the Habitat Directive (EC 1992;Bastmeijer 2018), and for neighboring areas, is strongly recommended. This is especially true due to the complexity of some territories and the multiplicity of human interests that require a more comprehensive management approach, especially regarding the inshore areas. A stronger cooperation among different stakeholders, which include dive centers, scientists, and fishermen, is also needed to share knowledge, best practices, and commitment to the NAS reef protection (Bertzky and Stoll-Kleemann 2009;Cvitanovic et al. 2014;Gianni et al. 2022).
In conclusion, the present work has highlighted the distinctiveness of the epibenthic assemblages inhabiting the NAS reefs, provided useful knowledge to refine the ecological quality assessment (Piazzi et al. 2023), and emphasized the need for better protection of these fragile habitats. Although the biogenic reefs investigated in this study probably represent only a small portion of the reefs in the NAS (Gordini et al. 2012), our results complement those of previous studies and represent the entire natural diversity of epibenthic assemblages in the area. However, further field studies are needed to better investigate the ecological dynamics and benthic communities of biogenic reefs that are very distant from the coast and have been poorly studied. In addition, a detailed and updated map of the distribution of the reefs, including the eastern part of the NAS, in combination with information on species genetic structure, inter-reef connectivity, and local human impacts are essential to achieve a reliable management of these marine habitats in the NAS.
Funding This research was supported by the Research Projects of National Interest (PRIN), funded by the Italian Ministry of University and Research: "Reef ReseArcH -Resistance and Resilience of Adriatic Mesophotic Biogenic Habitats to Human and Climate Change Threats" (Call 2015;Prot. 2015J922E4;. The manuscript was finalized during the participation of MP to the EU funded project "CoAStal and marine waters integrated monitoring systems for ecosystems proteCtion AnD managemEnt" (CASCADE; Interreg V-A Italy-Croatia CBC Programme 2014-2020, ID Number 10255941, from 01.01.2020 to 31.12.2022), and the participation of VB to the EU funded project "Ecological Observing System in the Adriatic Sea: oceanographic observations for biodiversity" (ECOSS; Interreg V-A Italy-Croatia CBC program 2014-2020, ID Number 10042301, from 01.01.2019 to 30.06.2021). FG and VB acknowledge financial support from PNRR MUR in the framework of the "National Biodiversity Future Center -NBFC" (CUP F83B22000050001). FC and MP also acknowledge financial support from PNRR MUR project ECS_00000033_ECOSISTER. This work was part of ET's Ph.D. thesis at the University of Bologna. Open access funding provided by Istituto Nazionale di Oceanografia e di Geofisica Sperimentale within the CRUI-CARE Agreement.

Data Availability
The datasets generated and analyzed during the current study are available from the corresponding authors on reasonable request.

Statements and Declarations
Conflict of interest The authors have no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.