How Do Geological Structure and Biological Diversity Relate? Benthic Communities in Boulder Fields of the Southwestern Baltic Sea

Environmental factors shape the structure and functioning of benthic communities. In coastal zones of the southwestern Baltic Sea, boulder fields represent one of the most productive habitats, supporting diverse benthic communities that provide many ecosystem services. In this study, the influence of the geological characteristics of boulder fields on the biodiversity of associated hard-bottom communities was investigated at two different spatial scales (few kilometers and tens of kilometers). The analyses on overall richness (taxonomic and functional) and community composition revealed how: (i) locally the size of boulders and (ii) regionally site-specific factors like the boulder density distribution and the sediment distribution can act as environmental driving forces. The overall richness of assemblages was shown to increase with increasing surface area of boulders, by up to 60% for species and up to 40% for functional richness. At both investigated scales, differences in compositional variability (β diversity) of the communities were detected. Locally, smallest boulders hosted more variable communities (β diversity up to 2 times higher), while at the regional level, indications of a larger habitat heterogeneity featuring the highest β diversity were observed. This study exemplifies how geological habitat characteristics shape the biodiversity of boulder field communities. The obtained information could be considered in assessment strategies, in order to avoid misclassifications of habitats naturally limited in biodiversity, making a step forward to the desired objective of protecting, conserving, and managing boulder field communities in the study area and at other comparable sites.


Introduction
Coastal ecosystems belong to the global hotspots of biodiversity and productivity (Spalding et al. 2007). They are located at the interface between land, air, and marine waters, which all simultaneously influence resident communities (Cloern et al. 2016). Additionally, as a large share of the human population is concentrated along or near coasts, anthropogenic-driven impacts like coastal engineering measures, nutrient discharge, overexploitation of biological and geological resources, fishfarming, tourism, and the impact of climate change (among others) are increasing in these environments (Halpern et al. 2008). Many coastal ecosystems therefore experience habitat degradation and a loss of biodiversity, which slowly diminish their resilience (Lotze et al. 2006).
Conserving and restoring biodiversity requires the understanding of its driving processes and in this context, the separation of natural from human-mediated influences. For marine hard-bottom communities, the geological habitat characteristics represent important natural drivers, as they determine the heterogeneity and complexity of these abiotic environments (Le Hir and Hily 2005;Liversage et al. 2017). Both factors are thought to facilitate a higher biodiversity by increasing the availability of microhabitats and refuges from predation as well as changing recruitment patterns (Sebens Communicated by Paul A. Montagna 1991;Archambault and Bourget 1996;Kovalenko et al. 2012).
The coastal morphology of the southwestern Baltic Sea is the result of the last glaciation, the postglacial sea level rise, and current processes of erosion and accumulation, driven by waves and currents (Harff et al. 2011). The southern and southwestern Baltic Sea area is dominated by abrasion platforms, coarse-grained relict deposits (sand and gravel), sand and mud, with a general gradient from coarse to fine-grained sediments towards deeper waters (Niedermeyer et al. 2011). Moranic material left by the last glaciation provide boulders and stones which serve as the only natural hard-bottom substrate for benthic communities (Diesing and Schwarzer 2006;Schwarzer et al. 2014;Kaskela and Kotilainen 2017). In the southwestern Baltic Sea, natural hard-bottom substrates are mainly stones (6.3-20 cm) or boulders (20-630 cm) (Kolp 1966). Diverse assemblages of sessile organisms and mobile invertebrates reside here, which find food and shelter within the complex structures of stones and associated benthic flora (Le Hir and Hily 2005;Liversage et al. 2017). Furthermore, boulder assemblages serve as feeding areas for fish (e.g., Gadus morhua), birds, and marine mammals and are utilized by some fish species as spawning and nursery areas (Mikkelsen et al. 2013;Kristensen et al. 2017;Torn et al. 2017). Being a highly diverse community type within the Baltic Sea (Kautsky et al. 2017), hard-bottom assemblages associated with boulder fields provide various regulative (e.g., nutrient cycling), provisional (e.g., food resources), and cultural (e.g., recreational activities) ecosystem services (Rönnbäck et al. 2007).
Boulders in German and Danish coastal waters have been subject to intensive extraction down to − 20 m NN from about 1800-1974(Karez and Schories 2005Dahl et al. 2009). These boulders were mainly used for construction purposes, preferably as shore protection elements. Conservative estimates suggest that along the ca. 540 km long coast of Schleswig-Holstein (Germany), about 2.5 million boulders with a diameter from 60 to 153 cm have been removed, representing a surface area loss of 5.6 km 2 (Karez and Schories 2005;MELUR-SH 2012). As boulder assemblages are rather patchy in shallow environments of the Baltic Sea, the lost surface area implied considerable habitat degradation, even if the number of boulders is slowly increasing due to ongoing seafloor abrasion since the extraction has stopped (Bohling et al. 2009). European efforts to conserve marine habitats early addressed the high ecological value of stone and boulder accumulations. Consequently, those areas were included in the Habitats Directive (HD; habitat type 1170; European Commission 2007). The Water Framework Directive (WFD; Directive 2000/60/EC 2000), and the Marine Strategy Framework Directive (MSFD; Directive 2008/56/EC 2008) further extended the ambition of conserving characteristic habitats within their ecosystem-based approaches to reach a good ecological or environmental status. The principle of the ecosystem-based approaches is to link natural environmental factors with anthropogenic impacts in the area of concern (Borja et al. 2008;Long et al. 2015). Therefore, evaluation methodologies need to consider the correlation of biotic communities with the physical environment and apply this knowledge within robust and scientifically established ecological indicators (Borja et al. 2008). In contrast to other ecosystems like the North Sea (e.g., Michaelis et al. 2019), the relationship between geological characteristics and benthic hard-bottom communities has not been addressed sufficiently for the southwestern Baltic Sea. Thus, this study presents the results of a multidisciplinary approach joining geological investigations and biological assessment. We aimed to investigate the influence of boulder size (local factor) and the heterogeneity of the geological setting (boulder densities and sediment distribution; regional factors) on the biodiversity of hard-bottom communities in boulder assemblages.

Study Area
Three study regions of different exposure to wind and wave direction (Table 1) were selected along the shallow waters of the southwestern Baltic Sea coastline, consisting of Gelting Bay, Schönhagen, and Hohwacht Bay (Fig. 1). Gelting Bay is located at the eastern end of the outer Flensburg Fjord. This site stretches from Habernis in the northwest to the nature reserve of Geltinger Birk in the east. The geomorphology of Gelting Bay is unique in the southwestern Baltic Sea, as it is virtually unaffected by wind and waves except for strong winds from northwest. The cliff of Schönhagen extends from Damp in the south to Olpenitz in the north. This study region is directly exposed to wind and waves from eastern directions. Hohwacht Bay is located southwest of Fehmarn and extends from Todendorf in the west to Heiligenhafen in the east. The inner Hohwacht Bay is exposed to winds and waves from the west to the northeast.

Seafloor Mapping of the Study Regions
High-resolution hydroacoustic seafloor mapping techniques with sidescan sonar (SSS) were used to investigate the sedimentological conditions of the seafloor (Blondel 2009;von Rönn et al. 2019). A SSS (StarFish 452F, Tritech) was deployed from a rubber boat to detect stones and boulders as well as their spatial distribution. The SSS operates with a frequency of 450 kHz. Full coverage mapping was executed with a range of 50 m to each side and a line spacing between the profiles of 80 m, allowing 20% overlap. An average resolution of 0.1 m pixel −1 was achieved. Post-processing involving standard geometric and radiometric corrections (Blondel and Murton 1997;Lurton 2002) were done using SonarWiz 7.03 (Chesapeake Technology Inc.).
Ground truthing included underwater video observations as well as sediment sampling using a van Veen grab sampler to generate a sediment distribution map. Grain size distributions of sediment samples were obtained by mechanical sieving, using ¼ PHI intervals according to the mesh standard of the American Society for Testing and Materials (ASTM E11-09 2009). PHI is defined according to the following formula: where d is the grain diameter in millimeters. The grain size distribution for each sediment sample obtained by mechanical sieving is determined by the respective particle diameter (mm) and the corresponding class weight (%). The sediment composition was classified according to Folk (1954). For simplicity, stones and boulders were not explicitly differentiated by their size; hereafter, both are named "boulders." The SSS mosaics were used in combination with the sediment analysis and video recordings to create the sediment distribution map.

Geological Characteristics of the Seafloor and Exposure of Sampling Stations
Within each region, three stations in areas with high boulder density (boulder fields) at a water depth of 5 m were selected for the biological sampling (see below), based on the analysis of the SSS mosaics ( Fig. 1). For each station, the number of boulders with a diameter > 25 cm within one 100 × 100 m cell was determined and the mean distance between the boulders in each cell was calculated. Accordingly, within the full coverage SSS mosaic, all boulders were counted manually and the distances between individual boulders were measured, using the software ArcGis (Esri

Biological Sampling
Hard-bottom communities were sampled in July 2017. At each sampling station, samples were taken by SCUBA divers at 5 m water depth. Prior to the sampling, the divers inspected the sampling location in a radius of approximately 30 m in order to visually estimate and classify the overall size range of boulders. After this, the divers collected samples of the benthic communities from three small (diameter < 30 cm), medium (30 cm < diameter < 50 cm), and large (50 cm < diameter < 130 cm) boulders along a coast-parallel transect (total N of samples = 178). The samples were obtained by scraping off the attached biota within a 10 × 10 cm sampling frame and transferring them directly into zipper bags. In addition, the length, width, and height of each sampled boulder were measured. As a result of the used sampling frame, the smallest boulder that was sampled had a diameter of ≈ 15 cm. In case the diameter of a boulder exceeded 30 cm, sampling was repeated along the longest axis of the boulder with 30 cm separating neighboring sampling frames. All collected material was fixed within 3 h in a buffered formaldehyde solution (final concentration of 4%).

Data Processing
The community composition and species richness were recorded for each boulder. Since larger boulders contained several samples, species lists for areas of single frames were cumulated to obtain the total species inventory of the respective boulder. The surface area of each boulder was calculated applying an approximation for the surface area of an ellipsoid: with p = 1.6075 as a constant and a, b, and c representing the length, width, and height of the boulder, respectively (Thomsen 2004). This shape resembles the appearance of boulders found in the study area. As shown in Eq. (2), the obtained value for the surface area was divided by two, since approximately half of the surface of a boulder would be available for settlement, while the other half is embedded in the sediment or underlies erosion processes at the seafloorwater interface.
Biological samples were analyzed for species occurrence to species level where possible (Table S1). For the majority of the samples, perennial macrophyte species (foundation species hereafter) dominated the surface of the boulder and, by providing substrate and shelter, facilitated a diverse community of sessile as well as mobile species (Jones and Thornber 2010). Therefore, for each sample, the foundation species was recorded. In addition, the functional diversity of the communities was examined. Therefore, the recorded species were categorized based on adult body size, growth form, trophic type and motility (Table 2). Consequently, each species was Fig. 1 Overview of the sampling stations in Gelting Bay (G1-3), Schönhagen (S1-3), and Hohwacht Bay (H1-3) along the Baltic coast of Schleswig-Holstein, Germany. Bathymetry data were downloaded from GeoSeaPortal (https://.geoseaportal.de/mapapps/?lang=de; accessed 13 November 2019) assigned to a functional group (four letters code; Table S1) by the combination of its functional traits, which represent their ecological role in the community (Wahl 2009). As an example, the filamentous red alga Ceramium virgatum was classified as XFAA and the polychaete Nereis pelagica as LMDB. Based on this classification, functional richness was calculated as the number of functional groups.

Statistical Analysis
Species and functional richness were modeled in relation to boulder surface area using a Generalized Additive Model (GAM). As GAMs are semi-parametric and data driven, no a priori assumptions had to be made on the expected shape of the species-area relationship (SAR), which is classically described by a power function (Preston 1960). The function gamm4 from the package gamm4 (Version 0.2-5; Wood and Scheipl 2017) was implemented to fit GAMs with Poisson distribution and log-link function (species richness) or Gaussian distribution and log-link function (functional richness). The boulder surface area was included as a smooth term using a penalized cubic regression spline restricted to up to three degrees of freedom. Each of the adjusted GAMs was compared to a mixed version of the model (GAMM), including the region (Gelting Bay, Schönhagen, Hohwacht Bay) as random factor. According to the Akaike information criterion corrected for small samples (AICc), the included random factor (region) only improved the performance of the model for functional richness and in consequence was kept only in this case. In order to evaluate the adequacy of all adjusted models, the plots of residuals were examined. To test for potential correlations of the predictors, the surface areas of boulders were compared among the study regions using an ANOVA and no difference was detected (F 2, 78 = 0.455; p = 0.636) (Fig. S1).
To compare the taxonomic and functional richness of small boulders cumulatively reaching a comparable surface area of a single large boulder, the surface area of the boulders was divided into four size classes (A ≤ 0.5 m 2 ; 0.5 m 2 < B ≤ 1.0 m 2 ; 1.0 m 2 < C ≤ 1.5 m 2 ; D > 1.5 m 2 ). Then, in an iterative process (N = 100), samples of small boulders (size class A) were randomly generated from the overall set of boulders sampled for this size class. The obtained random samples were filtered for those reaching the mean surface area of large boulders (size class D; 2.8 m 2 ). Consequently, there was no difference between the surface areas of small cumulated boulders and large boulders (t tests; taxonomic: p = 0.2389, N = 41; functional: p = 0.1931, N = 43). The mean cumulated richness (taxonomic and functional) of small boulders was compared against the mean richness of large boulders using a t test.
The taxonomic and functional structure of communities (presence/absence data) was visualized using non-metric multidimensional scaling (nMDS) plots based on Bray-Curtis dissimilarities. The assemblages were compared among sampled regions and boulder size classes (see definition above). The nMDS were generated by 2000 random iterations using the package vegan (version 2.5-4; Oksanen et al. 2019). The β diversity of the compared communities was calculated as mean distance to the spatial group median (multivariate group dispersion) (Anderson 2006;Anderson et al. 2011) using the betadisper function (vegan package). In a permutational test (9999 iterations), the homogeneity of multivariate dispersions was compared between groups using the permutest function (vegan package).
To statistically examine differences in taxonomic and functional community composition and to identify the species and functional groups responsible for observed community patterns (nMDS plots), multivariate Generalized Linear Models (ManyGLM; sensu Warton et al. 2015) were adjusted using the manyglm function from the mvabund package (Version 4.0.1; Wang et al. 2019). In the ManyGLM, a generalized linear model (GLM) was fitted for each species or functional group, and the log-likelihood ratios of the adjusted models were summed and used as a multivariate test statistic by randomization (Warton et al. 2015). The models included surface area of boulders and the sampling regions as explanatory variables. The binomial distribution and the logit link function were used. A model selection process based on the sum of AIC (AIC sum ) was performed to compare between models generated by running all potential additive and interactive combinations of the considered explanatory variables (including the null model, Table 3). The plot of residuals of the final models was inspected to ensure that the meanvariance relationship associated with the chosen distribution was appropriate. The anova.manyglm function was used to calculate univariate test statistics and p values for individual species and functional groups (999 iterations).
After identification of species and functional groups significantly contributing to the assemblage structure (based on the adjusted ManyGLM), their occurrences were modeled in relation to the surface area of boulders using GAMs with binomial distribution (logit-link function). The identity of the regions was included as random effect in the models (according Functional groups are defined as the four letter combination of the traits presented (Wahl 2009) to the AICc). To further examine the regional differences in community composition that were identified by the ManyGLM, the occurrence of individual functional groups was modeled as a function of the regions using GLMs with binomial distribution (logit-link function). Here, only the functional groups were modeled, since the number of species identified by the ManyGLM procedure was too large to produce a single model for each of them. The GLMs were followed by Tukey HSD tests, performed using the function glht from the multcomp package (Hothorn et al. 2008), to evaluate differences in functional group occurrence among the sampling regions.

Sediment Distribution
The backscatter mosaics combined with ground truthing by sediment analysis indicated different sedimentological zones. The SSS mosaic and sediment distribution maps display variable sediment compositions and boulder distribution in all study regions ( Fig. 2; Fig. S2). The highest densities of boulder accumulations always appeared within zones were coarsegrained sediments dominate the seafloor (lag sediments). Grain size distributions of sediment samples from these zones displayed a sediment composition primarily consisting of sand and gravel between the boulders (Fig. S2a). Individual boulders also occurred in zones where sandy sediments (fine to medium sands) characterize the seafloor ( Fig. 2; Fig. S2). All biological sampling stations are within zones of the seafloor that are dominated by lag sediments. The shallow water zone of Gelting Bay is characterized by extended zones of lag sediments on abrasion platforms, which are frequently interrupted by areas consisting of finer sediments (sand-sandy mud) ( Fig. 2; Fig. S2). Boulders occurred mainly within the lag sediments. The sampling stations (G1-3) are located on abrasion platforms, where patches composed of sandy sediment with extensions from 50 to 700 m are observed. A persistent abrasion platform is located offshore the cliff of Schönhagen, demarcated in the south and north by sandy environments. All three sampling stations (S1-3) are located within this abrasion platform. The seafloor of Hohwacht Bay is also characterized by zones of variable sediment composition. The lag sediments (abrasion platforms) with boulder accumulations alternate with areas of finer sediments (muddy sand-fine sand). Compared to Gelting Bay, these alternating sediment zones of sandy sediment in between lag sediments are continuous and have a larger extend. The three sampling stations (H1-3) are located on an abrasion platform, where patches composed of sandy sediment (extensions from 50 to 150 m) are observed within the dominating coarse-grained sediments.

Boulder Densities at the Sampling Stations
The highest average number of boulders per 0.01 km 2 was detected for the stations at Schönhagen (2557 ± 220 boulders/0.01 km 2 ), followed by stations in Hohwacht Bay (1220 ± 932 boulders/0.01 km 2 ), and Gelting Bay (431 ± 169 boulders/0.01 km 2 ). These results directly reflect the mean distances between boulders, which are smallest at Schönhagen with distances of 1.2 ± 0.1 m, followed by stations in Hohwacht Bay with distances of 1.6 ± 0.5 m, and Gelting Bay, where the distances were 2.3 ± 0.4 m.

Richness-Boulder Surface Area Relationship
The taxonomic and functional richness (cumulated per boulder) increased with boulder surface area (Table S2). Species richness initially increased strongly with surface area, but tended to saturate at a richness of approximately 45 species at surface areas > 3 m 2 (diameter of ≈ 80-130 cm) (Fig. 3a). The functional richness exhibited a similar initial increase and saturated at a richness of ≈ 19 functional groups at surface areas > 2 m 2 (diameter of ≈ 70-130 cm) (Fig. 3b).

Richness on Comparable Areas of Small and Large Boulders
The mean cumulated taxonomic and functional richness on small boulders (representing the surface area of a large boulder) were higher than the mean richness found on single, large boulders (p < 0.0001; Fig. 4). Taxonomic richness standardized for area was 1.7 times higher, while functional richness was 1.4 times higher on small boulders.

Community Structure
The nMDS plots highlighted region-specific differences in taxonomic composition (Fig. 5a). The communities of Gelting Bay and Schönhagen showed a lower variability around the spatial group median than the communities of Hohwacht Bay. Accordingly, the β diversity was found to be 1.5 times and 1.6 times higher for the communities of Hohwacht Bay compared to communities of Gelting Bay (p = 0.0002) or Schönhagen (p = 0.0001), respectively ( Fig.  5c; Table S3). Functionally, the communities were less different (Fig. 5b), but again a difference in β diversity between Hohwacht Bay and the remaining regions was observed ( Fig.  5c; Table S3). Here, the β diversity of the communities in Hohwacht Bay was 1.4 and 1.2 times higher compared to Gelting Bay (p = 0.0453) and Schönhagen (p = 0.0068), respectively.
The comparison of the community structure among boulder size classes revealed differences for both taxonomic (Fig.  5d-f) and functional (Fig. 5g-i) composition. In both cases, differences in the variability around the spatial group median were most pronounced between the smallest size class (A) and all remaining size classes (B-D). Noticeably, assemblages on smallest boulders always showed the highest within-group dispersion, a pattern confirmed by the calculated β diversity (Fig. 5j-l). The β diversity based on taxonomic (functional) composition was 1.5-2.1 times (1.3-2 times) higher on smallest boulders compared to the largest boulder size class (D) (Table S3). Only for the functional composition in Schönhagen, β diversities were similar among all boulder size classes (Fig. 5k). For taxonomic composition, the communities also showed differences in the overall composition of the assemblages, as the samples of size class A exhibited a low overlap with samples of the other size classes. In contrast, the functional composition differed only in terms of β diversity but not in the identity of functional groups.
The applied ManyGLMs statistically confirmed that both considered parameters (size of boulders and the regions) affect the overall taxonomic and functional composition of the communities, explaining 31% and 24% of the variability, respectively (Table 4). For boulder surface area, univariate tests revealed that the occurrence of eight species from the phyla Cnidaria (4), Porifera (2), Entoprocta (1), and Rhodophyta (1) were responsible for the observed differences in taxonomic community composition (Table S4). Regarding functional community composition, the functional groups XESA (100-1000 mm, encrusting, suspension feeder, attached) and MMDS (1-10 mm, massive, deposit feeder, swimming) were responsible for the detected differences related to boulder size    (H) (a, b), as well as in relation to boulder size classes (A ≤ 0.5 m 2 ; 0.5 m 2 < B ≤ 1.0 m 2 ; 1.0 m 2 < C ≤ 1.5 m 2 ; D > 1.5 m 2 ) within each region (d-i). Lines connect samples to their respective spatial group median. Stress values indicate for each nMDS how well dissimilarities are preserved. For each group visualized in the nMDS plots, the β diversity, based on the mean distance to the spatial group median, is presented (c, j-l). Different letters and symbols represent significant differences (p < 0.05) between compared groups based on permutational tests. The exact p values can be found in Table S3 ( Table S4). Region-specific differences in community composition were driven by the occurrence of 46 species and six functional groups (Tables S4 and S5).

Single Species and Trait Models
All species and functional group occurrences (identified by the ManyGLMs) were shown to increase with boulder surface area (Table S6). The probability of the considered species to occur in an assemblage generally increased with increasing boulder surface area. The trends for representatives of the phyla Entoprocta (Barentsia gracilis), Porifera (Halichondria (halichondria) panicea, Leucosolenia botryoides), and one hydrozoan (Opercularella lacerata) exhibited steep initial slopes, reaching highest occurrence probabilities already at surface areas ≥ 1 m 2 (Fig. 6a-c). Less extreme slopes were described for the hydrozoan Campanulina pumila and the red algae Delesseria sanguinea that reached highest probabilities of occurrence at surface areas ≥ 3 m 2 (Fig. 6b, d). Two species of the phylum Cnidaria (Bougainvillia muscus, Obelia longissima) only occurred on larger boulders (> 1 m 2 ) (Fig. 6b). The occurrence of functional groups increased with larger surface area of boulders as well (Fig. 6e, f). Here, the highest probability of occurrence is already attained at a surface area of 1 m 2 .
The post hoc comparisons of functional group occurrences between sampled regions revealed differences for four out of six functional groups previously identified by the ManyGLM (Fig. 7; Table S7). Large grazers (LMGB, LMGC; Fig. 7a) showed highest occurrences in Gelting Bay, while a medium grazer (MMGC, represented by Jaera (jaera) albifrons only; Fig. 7b) was mostly recorded in Hohwacht Bay. The largest differences were recorded for the functional group MMSA (1-10 mm, massive, suspension feeder, attached), where the average probability of occurrence was 60% lower in Hohwacht Bay compared to Schönhagen and Gelting Bay (Fig. 7b).

Discussion
The present study provides insights on how geological features and seafloor characteristics could shape the structure of hard-bottom communities in the southwestern Baltic Sea. Our analyses indicated that locally the size of the boulders and regionally site-specific factors like the boulder density distribution and the sediment distribution can act as environmental driving forces at the level of overall richness (taxonomic and functional), community structure, and individual species and functional group occurrences. In this context, we showed that assemblages found on larger boulders are richer in species and functional groups (α diversity; Fig. 3) and exhibit a lower β diversity compared to boulders of smaller surface area (Fig.  5j-l). The species and functional groups that were identified to be responsible for these differences in community composition were found to generally increase in occurrence on larger boulders (Fig. 6). In addition, region-specific differences in β diversity (Fig. 5c) as well as in the occurrence of certain species and functional groups were shown ( Fig. 7; Tables S4 and S5).
The observed relationship between species richness and the surface area of sampled boulders follows the widely accepted paradigm in ecology, predicting more species with increasing area (Arrhenius 1921;Preston 1960). By applying an approach with only a few assumptions on the final shape of the model (GAMM), the resulting SAR reached an upper asymptote in our study (for taxonomic and functional data; Fig. 3). Here, factors like the restricted pool of benthic species (associated to boulders) sampled, the comparatively species poor communities of the Baltic Sea, and most importantly the relatively narrow size range (diameters of ≈ 15-130 cm) of the investigated boulders likely determined the observed asymptotic behavior of the generated curves (Williamson et al. 2008;Tjørve and Turner 2009;Ojaveer et al. 2010). While the models for taxonomic and functional richness exhibited similar trends, their saturation point differed, with the latter showing saturation at a lower boulder area (> 2 m 2 ). Consequently, the further increase in species richness beyond this boulder size should contribute to functional redundancy within the communities (Cumming and Child 2009). In turn, functional groups on smaller boulders are mainly represented by single species, making them less resilient against environmental fluctuations (Yachi and Loreau 1999). For the presented SAR, it should be noted that the richness on largest boulders can be expected to be even higher, since the sampling effort scaled with the diameter of the stones, but not their area. For each parameter the summed deviances and the p value are given. The applied ManyGLMs were based on a binomial distribution Consequently, the sampling could have been incomplete for larger boulders. On the local scale, the performed community analyses provided further information whether the assemblages of smaller boulders, shown to have less species and functional groups, host unique communities or if they represent smaller subsets of the communities found on larger boulders. A considerably higher β diversity was observed in the smallest size class of boulders, especially in taxonomic composition (Fig. 5j-l). Species and functional groups responsible for this difference were all shown to increase in occurrence with increasing area (Fig. 6). As none of these species or functional groups solely occurred on smaller boulders, it can be concluded that the smaller boulders did not host unique assemblages, but simply random subsets of the overall pool of species recorded for larger boulders. This indicates that the structure of the Fig. 6 Generalized Additive Mixed Models (GAMMs) relating the occurrence of species (a-d) and functional groups (e, f) to boulder surface area. The mean effect of boulder surface area is presented as solid lines; shaded areas represent the 95% confidence intervals. The marks on the x-axis show the distribution of the data for the considered predictor. The estimated degrees of freedom (df) and the p values of the adjusted spline are given communities is mainly driven by random placement, i.e., the likelihood of a species to occur rises with increasing size of the boulders (Coleman 1981). In the specific case of boulder fields, the overturning of boulders and the abrasion and burial by mobile sediments could introduce further variability to the community structure. As smaller boulders might be more susceptible to these effects (Schrottke et al. 2006), the following recolonization would increase the variability of the assemblages. In comparison, larger boulders are rarely moved and are less affected by covering with sediment, consequently hosting more stable communities (Osman 1977;Sousa 1979). Other mechanisms like the diversity of microhabitats (Williams 1943(Williams , 1964 seem to be less important, since in those cases more similar communities on the smaller boulders, and/or species occurring exclusively on larger boulders, could be expected. In light of efforts targeting the conservation and restoration of boulder assemblages, these findings provide important hints on the habitat quality related to the occurring boulder sizes. As our results showed that single large boulders on average host more species (Fig. 3), they should receive a higher priority in conservation or restoration projects. However, the cumulated richness on a set of smaller boulders, which together represent the average area of a large boulder, was higher than on single large boulders (Fig. 4). This is likely caused by the differences among communities on different small boulders, i.e., their β diversity. The β diversity is not considered in the mean diversity of large boulders, which should caution against an assumption of larger diversity on (many) small boulders. Also, such an assumption is not taking into account the redundancy of functional groups, which was found to be larger with increasing area per individual stone. Therefore, the communities on larger boulders have the potential to maintain important functions during disturbances, a feature that is seen to be similarly important as biodiversity itself (Walker 1995). Støttrup et al. (2017) provide an example from practice, reporting the successful restoration of a boulder assemblage at the Danish coast by deploying large boulders (up to diameters of ≈ 160 cm). Surveys following the restoration showed that the created habitat was resistant to severe weather and, by increasing complexity, offered substrate for more diverse and productive communities.
In addition to differences in community structure caused by local drivers (boulder size), the community analyses detected differences related to the regions. While differences in β diversity were largest for the taxonomic composition among regions, the applied ManyGLMs confirmed differences in both, single species and functional groups. This is in line with the observed changes in the occurrence of certain functional groups among regions (Fig. 7). Generally, functional groups containing grazers (LMGB, LMGC, MMGC) were found to occur more in Gelting Bay, while medium-sized suspension feeders (MMSA) showed highest occurrences in Gelting Bay as well as Schönhagen. Grazers play a vital role in top-down control of epiphytic algae in marine systems, thus being able to mediate effects of eutrophication (Worm et al. 2000;Alsterberg et al. 2013). Suspension feeders promote benthic-pelagic coupling, as their filter feeding activity deposits, circulates, and regenerates nutrients (Kautsky and Evans 1987). Therefore, we further discuss the differences in geological characteristics between the investigated regions that might be the basis for the observed variability in functional group occurrences.
The geological characterization of the studied regions revealed grain size compositions typical for lag sediments of the southwestern Baltic Sea, primarily consisting of sand and gravel between boulder accumulations (Seibold et al. 1971;Diesing and Schwarzer 2006;Schwarzer et al. 2014). Comparing the grain size compositions of the surface samples between regions, similar distribution patterns are observed. Here, it should be noted that the sampling procedure of the grab sampler is limited to fine-grained sediments (e.g., Coggan et al. 2007), thus the range of grain sizes is likely Fig. 7 Generalized Linear Models (GLMs) comparing the occurrence of functional groups among sampled regions. Mean probability of occurrence (circles) and 95% confidence intervals (whiskers) are presented for each functional group and region. Letters indicate significant differences between regions according to a Tukey HSD test (see Table S7 for excact p values) underestimated. However, if the sedimentological information is combined with hydroacoustic, visual, or on-site investigations by SCUBA divers, the sedimentological composition of the seafloor can still be assessed in full extent (von Rönn et al. 2019). Accordingly, differences in the overall shape and sedimentological zonation of the abrasion platforms as well as the average number of boulders found per 0.01 km 2 were described for the three study regions. Schönhagen featured the boulder assemblages with the highest density of boulder distribution (and smallest distances between boulders). Gelting Bay and Hohwacht Bay exhibit a higher spatial heterogeneity. Here, the occurring distributions of boulders are separated from each other by larger areas of finer sediments. Additionally, the average distances between single boulders are up to 1.9 times larger. These subdivisions by finer sediments increase the proportion of habitat edges, which were shown to increase the variability in community composition (Hewitt et al. 2005). Edge effects in marine habitats are thought to result from hydrodynamic changes at transition zones, affecting larval dispersal by settlement shadows (at leeward edges), changing predation patterns, and increasing the provision of food (Tilman 1994;Bologna and Heck 2002;Zajac et al. 2003). For rocky reef communities, Roberts and Poore (2005) provide evidence for beneficial effects of an increased patchiness of habitats, shown to enhance the colonization of epifauna on a macroalga. Taken together, the mentioned effects of habitat subdivisions could have promoted a higher β diversity of the communities in this study, as seen for the communities of Hohwacht Bay (Fig. 5a) (Tilman 1994;Chesson 2000). The positive relationship of β diversity and habitat heterogeneity has been demonstrated in previous studies on benthic systems (Ellingsen and Gray 2002;Hewitt et al. 2005;Anderson et al. 2011). However, it should be noted that, at the regional level, other environmental factors like exposure or sedimentation effects during stormy conditions could additionally be responsible for the observed community patterns (Balata et al. 2007). In fact, differences in exposure could explain the lower β diversity of Gelting Bay compared to Hohwacht Bay, as these stations share similar geological and environmental settings. The average exposure of Hohwacht Bay (114,386 m 2 s −1 ; Table 2) is 3.6 times higher than for Gelting Bay (31,644 m 2 s −1 ), which additionally could promote variabilities in recruitment rates and survival of species, leading to a higher β diversity (Todd 1998). The potential overlap of effects related to habitat heterogeneity and other environmental factors (exposure, sedimentation) shaping the community structure on the regional scale, shows that the proposed assumption of the presented geological characteristics regionally driving patterns of β diversity remains speculative and needs further empirical evidences.
Throughout the history of human civilizations, coastal ecosystems have been a central target for resource exploitation and settlement (Lotze et al. 2006;Halpern et al. 2012). The resulting habitat degradation and losses in biodiversity led to various efforts during the last five decades, aiming to manage, conserve, and restore these environments (Boyes and Elliott 2014). For the majority of these efforts, the state of the considered habitat needs to be assessed before any measure can be applied. Assessing the condition of coastal habitats is challenging, as strong natural variabilities act in combination with anthropogenic impacts (Mieszkowska et al. 2014;Franz et al. 2019). In the specific case of the Baltic Sea, strong environmental gradients further exacerbate the definition of reference values. For example, Barboza et al. (2019) showed how fitness-related traits of the habitat-forming macroalga Fucus vesiculosus vary in response to the prevailing salinity gradient as well as other local environmental drivers (e.g., loads of nutrients), potentially affecting the functioning of associated communities. Thus, in hard-bottom habitats of the Baltic Sea, the analysis of the main processes shaping the structure of communities and functioning is especially important to disentangle natural drivers of biodiversity from factors caused by human activities. Torn et al. (2017) acknowledge this source of variability within their assessment system for different habitat types by including a so called "ecological zoning." This classification allowed the exclusion of depthand exposure-related biases within the calculation of a habitat quality index (Torn et al. 2017). Our study demonstrates that the geological characteristics of boulder assemblages, namely the size range of boulders, also have the potential to locally act as a natural driver of biodiversity. Therefore, attempts to evaluate boulder field habitats in the southwestern Baltic Sea should consider the size range of available substrate as a classification factor, e.g. in the framework of an ecological zoning. In this manner, misclassifications of habitats that are naturally limited in biodiversity as a result of abiotic conditions could be prevented. On the regional scale, we described tendencies of geological factors (e.g., habitat heterogeneity) shaping the structure of hard-bottom communities, but further evidences will be needed to clearly disentangle effects related to the geological setting from other regional factors. As an additional biological factor, the macrophytes identified as foundation species could be listed as typical taxa for communities at the sampled depth range, as their absence could strongly reduce habitat complexity. Since marine directives like the WFD and the MSFD require the assessment of boulder assemblages in the study region, our insights could contribute to develop or improve assessment systems in order to reach the desired objective of protecting, conserving, and managing boulder habitats in the southwestern Baltic Sea.

Conclusions
The study of hard-bottom communities in boulder fields of the southwestern Baltic Sea revealed how geological features and seafloor characteristics influence the richness and community structure at different spatial scales. Locally, a clear influence of boulder size was detected, at the level of overall richness (taxonomic and functional) and for community structure. Assemblages on smallest boulders exhibited the highest β diversity, a pattern most likely explained by random placement and/or earlier successional stages after disturbance. Region-specific differences in community composition were related to differences in boulder densities and sediment distribution, indicating varying degrees in habitat heterogeneity among the studied regions. However, at this scale, other factors (e.g., exposure and sedimentation effects) can be expected to influence the community composition as well. The results of this study highlight natural drivers of richness and community composition, which should find consideration in efforts to conserve or restore boulder field habitats in the study area.
Acknowledgments We thank Renate Schütt and Gesche Bock for their patient support during species identification, Eric Steen for the configuration of the hydroacoustic survey equipment, and Nils Steinfeld for his support during our research surveys. We also like to thank the research divers of Kiel University for their valuable field work. We wish to thank the crew of the research vessel "Littorina" for their excellent support.
Funding Open Access funding enabled and organized by Projekt DEAL. MF and GAvR acknowledge the financial support by the State Agency for Agriculture, Environment and Rural Areas of Schleswig-Holstein (LLUR; reference number: 0608.451426). FRB acknowledges the financial support by the German Academic Exchange Service through the project Doctoral Programmes in Germany 2015/16 (57129429). Information contained here has been derived from data that is made available under the European Marine Observation Data Network (EMODnet) Seabed Habitats project (http://www.emodnet-seabedhabitats.eu/), funded by the European Commission's Directorate-General for Maritime Affairs and Fisheries (DG MARE).

Compliance with Ethical Standards
Conflict of Interest The authors declare that they have no conflict of interest.
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/.