Wave stress and biotic facilitation drive community composition in a marginal hardbottom ecosystem

Ecological patterns are inherently scale-dependent and driven by the interplay of abiotic gradients and biotic processes. Despite the fundamental importance of such gradients, there are many gaps in our understanding of how abiotic stress gradients interplay with biotic processes and how these collectively affect species distributions. Using a hierarchical design, we sampled two communities separated by depth along wave exposure and salinity gradients to elucidate how these two gradients affect species composition in habitats formed by the foundation species Mytilus trossulus and Fucus vesiculosus. Specifically, we looked at the impacts of regional salinity and temperature, local wave exposure, and site-dependent facilitation effects on the associated community composition. Wave exposure was the best predictor for species assembly structure, which was also affected byMytilus biomass and by salinity and water temperature. While the tested variables provided robust explanations for community structure and density, they did not provide conclusive explanations for variation in species richness or evenness. Mytilus biomass had a stronger effect on the associated community with increasing wave exposure at the deeper depth, but the patterns were less obvious at the shallower depth. The latter was also the case for Fucus. These findings comply partly with theoretical predictions suggesting stronger facilitation effects in physically harsh environments. Our results indicate that environmental drivers are the main structuring forces that affect species assembly structure, but also foundation species are important. Thus, predicting changes in species distributions and biodiversity requires the simultaneous consideration of environmental gradients, as well as the structure and composition of foundation species and the interplay between these factors. This work advances our understanding of the processes that modulate species distributions in a marginal marine area and broadens the knowledge of how biological and environmental factors interplay and have an influence on hard-bottom community structure in brackish water seas.

ents and biotic processes. Despite the fundamental importance of such gradients, there are many gaps in our understanding of how abiotic stress gradients interplay with biotic processes and how these collectively affect species distributions. Using a hierarchical design, we sampled two communities separated by depth along wave exposure and salinity gradients to elucidate how these two gradients affect species composition in habitats formed by the foundation species Mytilus trossulus and Fucus vesiculosus. Specifically, we looked at the impacts of regional salinity and temperature, local wave exposure, and site-dependent facilitation effects on the associated community composition. Wave exposure was the best predictor for species assembly structure, which was also affected by Mytilus biomass and by salinity and water temperature. While the tested variables provided robust explanations for community structure and density, they did not provide conclusive explanations for variation in species richness or evenness. Mytilus biomass had a stronger effect on the associated community with increasing wave exposure at the deeper depth, but the patterns were less obvious at the shallower depth. The latter was also the case for Fucus. These findings comply partly with theoretical predictions suggesting stronger facilitation effects in physically harsh environments. Our results indicate that environmental drivers are the main structuring forces that affect species assembly structure, but also foundation species are important. Thus, predicting changes in species distributions and biodiversity requires the simultaneous consideration of environmental gradients, as well as the structure and composition of foundation species and the interplay between these factors. This work advances our understanding of the processes that modulate species distributions in a marginal marine area and broadens the knowledge of how biological and environmental factors interplay and have an influence on hard-bottom community structure in brackish water seas.

INTRODUCTION
Foundation species (FS), such as many trees, shrubs, corals, bivalves, macroalgae, and seagrasses, are species that often have a disproportionately large influence on community structure and ecosystem function (Dayton 1975). They are vital for many systems as they have a propensity to increase diversity, biomass, or ecosystem stability through, for example, enlarging niche space (Bulleri et al. 2016), enhancing habitat complexity (Kostylev et al. 2005), reducing physical stress (Bertness and Callaway 1994), supporting gas exchange (Attard et al. 2019), and increasing resource availability (Norkko et al. 2006, Norling and. High diversity is important, as rich communities are better equipped with functions that buffer against future changes in the environment, caused by stressors such as climate change, species introductions, or eutrophication. The community effects of FS depend on the environmental settings and are often conditioned by their individual-or population-level properties (density, age, size). Net effects of FS may therefore vary from positive, neutral to negative (e.g., reviewed in Bateman and Bishop 2017). The stress-gradient hypothesis (SGH; Bertness and Callaway 1994) suggests that the strength of facilitation and competition should be inversely related. At the highstress end, the effect of facilitation on species composition should be strong as facilitation alleviates environmental stress allowing a higher diversity or abundance of the associated community. Yet, mainly due to competition effects, where FS may monopolize resources, interactions may be neutral or even negative in physically low-stress environments (Bruno et al. 2003, Crain andBertness 2006). The influence of FS should increase with the degree to which they alleviate stressful conditions. Diversity may therefore peak in regions where habitat-forming species are most abundant and where physical structure is a limiting resource (Kimbro and Grosholz 2006). Thus, many factors affect the facilitating mechanisms of habitat-forming species and the challenge is to identify circumstances under which these habitat-forming species alleviate stress and provide facilitating services to the community (Crain and Bertness 2006). Two decades ago, Bertness et al. (1999) called for a more diversified research related to biogenic habitat provision. Even if there has been a renaissance in studies related to positive interactions since then, facilitation as a mechanism remains understudied in many systems (He andBertness 2014, Wright andGribben 2017). Moreover, while most studies have focused on pairwise species interactions (Brooker et al. 2008), there has been a call for studies on how facilitation affects whole communities (Soliveres et al. 2015, Wright andGribben 2017) and how cooccurring or adjacent habitat-forming species in conjunction determine the species assemblage (Yakovis et al. 2008, Bishop et al. 2013. Apart from facilitation, many environmental drivers affect species distributions that ultimately set the local and regional species pool. These drivers occur in spatial hierarchies, they are connected, and they integrate the interplay of ecological processes that form the patterns we observe (Hewitt et al. 2017). In the Baltic Sea, for example, sea salinity forms an upper-level driver that ultimately determines how deep into the Baltic Sea marine species penetrate (Remane 1934). Another upper-level regional driver is temperature that may determine how south cold-water species are established. A local driver is formed by, for example, wave exposure that may affect focal response variables, such as species interactions. Ecological interactions, however, often cross levels of organization or scales, and the interplay of processes at different scales may be much more complicated than first anticipated (Soranno et al. 2014, Hewitt et al. 2017. Consequently, the inclusion of multiple scales in ecological studies is important in order to test predictions about the influence of processes on ecological patterns, as the effects of these processes may be highly context-dependent (Hewitt et al. 2017). Limiting the scope and focusing on only on a narrow set of scales may also lead to false conclusions (Hewitt et al. 2017).
Here, we examine the community structure along wave exposure gradients in three nearby regions that are characterized by slightly different salinity conditions. In these regions, two habitat-forming species live at the limit of their lower salinity tolerance (Vuorinen et al. 2015). Specific interest lies in the analysis of changes in ❖ www.esajournals.org benthic species composition in rocky shore habitats dominated by the brown macroalga Fucus vesiculosus and the blue mussel Mytilus trossulus (henceforth Fucus and Mytilus). Both Fucus and Mytilus are dominant space occupiers in the area and provide several important ecosystem services including supporting and regulating services. Mytilus has a central role in the Baltic Sea food web. Blue mussels circulate nutrients and enhance biodiversity both directly, via substrate provision, and indirectly, by affecting near-bed hydrodynamics, affecting sedimentation rates, and clearing the water (Norling and Kautsky 2008 with references therein). Fucus spp. is the only perennial canopy-forming macroalga in the northern Baltic Sea and provides living space for a diverse assemblage of species. It also functions as a spawning, breeding, and foraging ground for many fish (Saarinen et al. 2018 with references therein).
While studies generally have looked at single systems along gradients of stress or at onedimensional gradients of disturbance in the same direction with a non-nested nature, we focus simultaneously on two environmental gradients and two FS and ask if their facilitating functions show responses to changes in the environment. In this rocky shore system, the importance of mussels and seaweeds as FS has experimentally been shown at patch scales (Koivisto and Westerbom 2010). However, it has largely remained unclear how different environmental gradients affect their facilitating functions. Nor have studies of scale-dependent biodiversity been carried out in a hierarchically balanced way, covering entire wave exposure gradients.
The purpose of this study is to examine and define the sublittoral hard-bottom community in terms of its structure and composition along two nested stress gradients-a regional and physiological salinity gradient and a local and physical wave exposure gradient. Our specific hypotheses were as follows: (1) Community structure and composition differ along environmental gradients and show responses that relate to the quantity of the two habitat-forming FS. Salinity is a highly important driver in this ecosystem with large effects on the resident Mytilus population, which shows evident structural changes when moving from west to east toward decreasing salinity (Westerbom et al. 2002(Westerbom et al. , 2019. Based on prior experimental work (Koivisto and Westerbom 2010), and knowledge of the regional differences in Mytilus population structure and biomass (Westerbom et al. 2002(Westerbom et al. , 2019, we expected to see the most abundant community in the west, where the biomass of Mytilus is highest and the bed structure is the most complex. Similarly, we expected to see a less dense community in the east, where the biomass of Mytilus is the lowest within the region. Since differences in Fucus biomass were not expected among regions, as Fucus is less responsive to salinity changes, we did not expect to see any region-wide effects of Fucus. Similarly, as the structure and biomass of the habitat-forming species vary among wave exposures (Westerbom and Jattu 2006), these changes in structure and composition were expected to be seen also in the local resident communities. According to our second hypothesis (2) the facilitating effects of the FS (measured as high diversity and/or abundance of associated fauna) are most obvious at sites with high wave exposure, thus supporting the stress-gradient hypothesis (SGH; Bertness and Callaway 1994).

Study area and species
The Baltic Sea is among the largest brackish water areas in the world offering a physically demanding, low-saline, temporally ice-covered, non-tidal, and evolutionarily young environment. The species pool is small, but species often occur in great density. The coastline of the study area, the western Gulf of Finland, is dominated by sublittoral rocks and reefs that typically are covered by dense blue mussel beds extending from the water surface down to more than 30 m depth with densities sometimes exceeding 10 5 individuals/m 2 (Westerbom et al. 2008). Fucus is the largest perennial alga in the study area, and it is typically belt-forming at shallow depths or organized in patches maintained by a combination of availability of light, ice abrasion, herbivory, and competition (Korpinen et al. 2007, Kraufvelin et al. 2007. Gradients in salinity, and to a lesser extent temperature, characterize the Gulf of Finland (e.g., Engstr€ om-€ Ost et al. 2015, Westerbom et al. 2019). Sea temperature increases slightly toward the east, whereas salinity declines from west to east. Even if the decline ❖ www.esajournals.org in salinity is subtle within the focal area (from 6.0 to 5.6), it is reflected in the Mytilus population as higher stability in the west where also bigger mussels are found compared to the east, where populations are much more dynamic (Westerbom et al. 2019). In addition to salinity, wave exposure is a major structuring force affecting local distribution patterns (Westerbom and Jattu 2006). Fucus populations have not been systematically surveyed in a comparable way in this region, but the coverage of Fucus is approximately similar among regions. Within regions, the coverage of Fucus is typically higher at intermediately exposed sites compared to very sheltered or to very wave-exposed sites.

Sampling
Using identical sampling methods, we sampled mussel bed communities and macroalgal communities in three regions along a 60-km stretch of coastline in SW Finland (approximately 59°45 0 N, 22°40 0 -59°47 0 N, 23°42 0 ). These regions were Hanko in the west, Tv€ arminne in the center, and Jussar€ o in the east (Fig. 1). We collected samples from 36 stations (12 per region) including a range of wave exposures where stands of Fucus and beds of Mytilus occur. To estimate wave exposure, we used the Baardseth wave exposure index (BI; see Westerbom and Jattu 2006), and consistent with recommendations of He and Bertness (2014), we grouped the stations into four exposure classes (sheltered [s], BI = 0-2; moderately sheltered [ms], BI = 7-9; moderately exposed [me], BI = 13-15; and exposed [e], BI = 18-22). Each exposure level included three stations per region (Fig. 1). Because shore slope affects the distribution of habitat-forming mussels (Westerbom et al. 2008), we chose only islands with a smooth rocky slope. We executed the sampling hierarchically, with samples nested within stations and stations nested within exposure levels and regions.
Within each station, we took ten samples, five from each of the two depths, 3 and 8 m, at distances of 10-20 m within depths (total number of samples, N = 360). At each station in August-September, we collected benthic samples from slopes opening toward the highest wave exposure, using a 20 9 20 cm frame with a 0.1-mm net bag. Samples were systematically dispersed in time within regions and exposures to avoid temporal biases in sampling. The method is standard in the Baltic Sea and used in, for example, national monitoring programs of the zoo-and phytobenthos (Kautsky and Van der Maarel 1990, Westerbom et al. 2002, R aberg and Kautsky 2007, Wikstr€ om and Kautsky 2007. Then, we processed the samples in the laboratory and counted the species according to common practices (Koivisto and Westerbom 2010, D ıaz et al. Stations shown with letters, where s denotes sheltered, ms moderately sheltered, me moderately exposed, and e exposed. ❖ www.esajournals.org 2015). Biomasses for Mytilus were estimated according to Westerbom et al. (2002). Fucus and other algae in samples were dried at 60°C for 5 d and then weighed.
In addition to sampling, we lowered a 100-m transect line perpendicularly from the shoreline and stretched it out to a minimum depth of 15 m. Once in place, a scuba diver estimated the coverage of flora and fauna over an area of 4 m 2 on both sides of the transect for every depth meter. Then, we extracted transect data from 3 and 8 m depths for use in multivariate models as measures of seascape structure. The coverage of Fucus, perennial algae, and red algae was used as a variable for seascape structure. To account for spatial variability and the lack of environmental monitoring data from each of the 36 stations, we extracted environmental data from coastal data models that are based on Finnish governmental in situ monitoring or satellite data and hosted by the Finnish Environment Institute. These data included sea surface salinity, water transparency (i.e., Secchi depth), phytoplankton chlorophyll a (Chl-a), and sea surface temperature.

Data analysis
We applied parametric univariate techniques to the facilitator variables (Mytilus biomass and abundance, Fucus biomass) and to summary variables of the associated macrofauna community (total abundance, Margalef species richness, Pielou's evenness) to test for differences among regions, wave exposure levels, and sampling stations. Due to the broad difference between mussel and seaweed habitat structure and the lack of Fucus at 8 m depth, we ran the analyses separately for 8 m depth (only facilitation by Mytilus) and for 3 m depth (Fucus and Mytilus facilitation).
For univariate analyses of the facilitator variables and summary variables for the associated macrofauna, we used three-way nested mixedmodel ANOVAs (Underwood 1997) following the model: where R i denotes level i of the region (fixed), E j denotes level j of exposure level (random), and S (RE) k(ij) denotes level k of station (random) nested in exposure level j and region i. RE ij represents an interaction component of the ith level of factor region and the jth level of factor exposure level, and e l(ijk) indicates the difference between an individual sample at station k of combination i and j of region and exposure level, respectively (Underwood 1997). Wave exposure was treated as random factor because the used wave exposure values represent a continuum of exposure levels rather than clearly separated classes (see Underwood 1997). We used a posteriori Student-Neuman-Keuls (SNK) tests to reveal the occurrence and direction of pairwise differences.
We studied facilitation effects with generalized additive mixed models (GAMMs) since they open the possibility to examine non-linear relationships (Wood 2006). We selected one model for each of the two depths to describe the variation in abundance of associated individuals between the two setups of habitat-forming FS. Initially, the factors used in each model were region (three levels) and wave exposure (four levels), both representing gradients of stress and disturbance, respectively. Since region had no effect, the final model excluded region as a factor. We considered the station factor as a random effect following the guidelines listed by Zuur et al. (2010). In the first model (8 m), Mytilus was the only foundation species, while in the second model (3 m), the factor for FS (BiomassFS) contained two levels, the biomasses of Fucus and Mytilus. On the continuous variable (abundance), we applied a smoother to counteract for residual heteroscedasticity due to non-linearity in the facilitation relationships. In the 3 m model, we used one smoother for each of the two FS. The dependent variable for associated fauna, that is, abundance of associated individuals, was discrete and over-dispersed showing a negative binomial distribution, and no transformation was applied to this variable. We applied the model selection following the criterion established by Zuur et al. (2010) on which outliers were identified and eliminated. We examined collinearity by plotting covariates together, while we detected variance inflation by the VIF index (Zuur et al. 2010). Collinearity was not an issue in these models. In both models, we used a cubic spline smoother (Zuur 2012) and we controlled overfitting using the integrated square second derivative cubic spline penalty (Wood 2006, Zuur 2012. Then, we used backward AIC selection (AIC, Akaike information criterion) to find the ❖ www.esajournals.org optimal models (Zuur et al. 2010). Finally, we applied model validation to verify that there were no violations of heteroscedasticity. Specifically, we plotted the relationship between residuals and fitted values. Afterward, we plotted residuals against each co-variable to investigate possible model misfit. The models we used were as follows: To visualize and to test for differences in community structure, we ran nonparametric multivariate techniques (Anderson et al. 2008) on the entire macrofauna data set applying the same model structure as for the three-way nested mixed-model ANOVAs above. Using three-way nested PERMANOVA, separately for the two study depths, we thus tested the data for differences in community structure among the three regions, four exposure levels, and 36 sampling stations. To even out the relative influence from very dominant and more rare species in the multivariate analyses, we square-root-transformed the data before running Bray-Curtis similarity analyses. Furthermore, we analyzed the association between the observed species composition pattern and the registered environmental variables using distance-based linear modeling, DISTLM, separately for each depth and using normalized environmental variables. Pearson's correlation coefficient indicated low correlation between environmental variables. We ran the DISTLM analyses using the "Best" option within the PRIMER software, which examines all possible combinations of predictor variables within the PRIMER software Gorley 2006, Anderson et al. 2008). Within this process, we initially obtained the best models for a series of alternative models including an increasing number of explanatory variables, going from a model including only one variable and ending with a model including all potential predictor variables. Then, we identified the most parsimonious model among these by AIC and we visualized the final model by distance-based redundancy analysis, dbRDA, which is an ordination technique constrained to find linear combinations of predictor variables that explain the greatest variation in the data cloud (Anderson et al. 2008). We assessed the relative influences of each response variable and environmental variable based on the length of their overlaid vectors in the resulting ordination plot. In addition, we identified the individual relationship of each environmental variable to the observed pattern in the species data set based on marginal F tests.
We ran the initial univariate analyses in GMAV5, the GAMM analyses in R version 3.4 package mgcv, and the multivariate statistical analyses in PRIMER 7.0 and PERMANOVA+ for PRIMER. Before running the parametric univariate tests, we checked the normality by Kolmogorov-Smirnov's test and homogeneity of variances by Cochran's C test. To homogenize variances, we used a square-root or logarithmic transformation.

Distribution of Mytilus and Fucus along gradients and their effects on fauna
At 8 m depth, there was a significant interaction effect (P = 0.01) between region and exposure for Mytilus biomass ( Fig. 2A; Appendix S1: Table S1a). Hanko showed higher Mytilus biomass than Tv€ arminne and Jussar€ o in the three lowest exposure levels, while for the highest exposure level, Tv€ arminne showed the highest biomass. At Hanko and Tv€ arminne, there were also significant withinregion differences ( Fig. 2A), whereas at Jussar€ o, no exposure differences were seen. Mytilus abundance showed no region-wide effects, although there were highly significant exposure level effects (P < 0.001) with abundances generally increasing with increased level of wave exposure ( Fig. 2A; Appendix S1: Table S1A).
At 3 m depth, a significant interaction (P < 0.01) between region and exposure was seen for Mytilus biomass ( Fig. 2B; Appendix S1: Table S1b) arising from the very high biomasses at sheltered shores in the Hanko region and no exposure differences within the other two regions (Fig. 2B). Mytilus abundance showed Fig. 2. Distribution of the predictor variables Mytilus biomass, Mytilus abundance, and Fucus biomass and the response variables richness, evenness, and total abundance of associated macrofauna for (A) 8 m depth and (B) 3 m depth. Note differences in scale among y-axes. Bars are means + SE. S stands for sheltered, MS for moderately sheltered, ME for moderately exposed, and E for exposed. neither region nor exposure effects. Fucus biomass at 3 m depth showed no regional differences, but showed wave exposure effects (P = 0.013), with the medium-sheltered region (ms) presenting higher biomasses than the other three exposure levels averaged over the regions (Fig. 2B; Appendix S1: Table S1b).
In total, we identified 52 invertebrate taxa or taxon groups associated with Mytilus or Fucus. At 8 m depth, the total number of individuals of associated fauna showed a strong response to the interaction between region and exposure (P < 0.001) that largely followed the results of Mytilus biomass (Fig 2A; Appendix S1: Table S1c) and varied with exposure (P = 0.005). Species richness showed a small effect of region at 8 m (P = 0.02), due to higher richness at Jus-sar€ o compared to Tv€ arminne and Hanko. No differences in evenness were seen among regions, nor among exposure levels (Appendix S1: Table S1c). At 3 m depth, the total abundance of associated macrofauna showed a weak exposure level effect that was due to constantly higher abundances at the most sheltered stations compared to the other wave exposure levels (P = 0.02; Fig. 2B; Appendix S1: Table S1d). No regional effect on richness was seen, but richness was higher at the two lower exposure levels than at the other two levels (P = 0.002). Evenness was slightly lower at Hanko compared to the other regions (P = 0.022; Appendix S1: Table S1d). Differences in evenness were, however, very small within and among regions (Fig. 2B).

Facilitation effects along stress gradients
Wave exposure had a clear effect on the relationship between the biomass of Mytilus and the abundance of associated species. This relationship was non-linear at both depths, showing that the facilitating function of Mytilus increased at higher Mytilus biomasses (P ≤ 0.001 in both cases). At 8 m, the facilitating effect increased sigmoidally with increasing biomass of Mytilus (Appendix S2: Table S1a). The facilitating effect was strongest at the most exposed stations and approached an asymptote at all exposures, indicating that Mytilus had reached a biomass that maximized its facilitating potential ( Fig. 3A; Appendix S2: Table S1a). At 3 m depth, the effect of the facilitating function of Mytilus resembled the exponential function with stronger positive effects at sheltered and exposed stations than at intermediate stations ( Fig. 3B; Appendix S2: Table S1b). The smoother on Fucus showed that while there was a positive linear relationship between the biomass of Fucus and the density of the associated community, the importance of Fucus appeared not to change along the wave exposure gradient (Fig. 3B; Appendix S2: Table S1b).

Associated fauna: multivariate responses
The community structure of associated macrofauna at 8 m depth showed a significant interaction effect between region and exposure (P < 0.01). There were also highly significant differences among all stations nested in region and exposure (P < 0.001; Appendix S3: Table S1). At 3 m depth, the structure of the associated macrofauna community differed significantly among exposure levels as well as among stations (P < 0.001), but no differences occurred among regions. Pairwise tests revealed that the differences between exposure levels occurred between all pairs, except for the two intermediate (ms and me) and the two outermost levels (me and e). At 8 m depth (with only Mytilus as facilitator), we identified the most parsimonious DISTLM model for a combination of nine environmental variables with wave exposure (21%), Mytilus abundance (20%), temperature (14%), salinity (13%), and Mytilus biomass (13%) as the strongest predictors ( Fig. 4; Appendix S4: Table S1). Samples from different exposure levels grouped distinctly and systematically along the dbRDA1 axis in the ordination (Fig. 4). Mytilus abundance, wave exposure, and salinity were strong predictors for positive values along the dbRDA1, while temperature predicted negative values. The samples from different regions were largely grouped along the dbRDA2 from top (mainly Hanko) to bottom (mainly Jussar€ o). Mytilus biomass was a strong predictor of positive values along the dbRDA2, while the transect coverage of red algae was a strong predictor of negative values. The taxa showing the strongest influence on the observed patterns were amphipods of the genus Gammarus that demonstrated a positive relationship with wave exposure and the bivalve Mya arenaria, which correlated negatively with wave exposure (Fig. 4). The second axis was mainly influenced by the gastropod Hydrobia, which Results from 3 m depth for both Mytilus and Fucus. S stands for sheltered, MS for moderately sheltered, ME for moderately exposed, and E for exposed. Fig. 4. Distance-based redundancy analysis (dbRDA) for square-root-transformed species abundance data from 8 m depth with the most important discriminating environmental/background variables and species as vector overlays. As shown, different regions (along the y-axis) and different exposures (along the x-axis) appear in different ends of the cluster, largely paralleling differences in Mytilus bed structure among regions and exposures. In the legend, H stands for Hanko, T for Tv€ arminne, and J for Jussar€ o, while s stands for sheltered, ms for moderately sheltered, me for moderately exposed, and e for exposed. showed a positive relationship with Mytilus biomass.
We identified the most parsimonious model at 3 m depth (with both Mytilus and Fucus as facilitators) for a combination of 12 environmental variables. Marginal tests identified wave exposure (22%) as the single variable that was most strongly attributed to changes in the associated macrofauna community, followed by temperature (16%) and salinity (7%; Appendix S5: Table S1). Fucus and Mytilus, however, explained only a fraction of the variation (2% and 4%, respectively). As Fig. 5 shows, samples from different exposure levels were grouped along the dbRDA1 with samples from sheltered stations to the right and samples from increasingly more exposed stations further to the left, the latter also correlating positively with salinity and negatively with temperature. With regard to dbRDA2, this axis encompassed a considerably smaller share of the explained variation with biomass and abundance of Mytilus as predictors of Fig. 5. Distance-based redundancy analysis (dbRDA) for square-root-transformed species abundance data from 3 m depth with the most important discriminating environmental/background variables and species as vector overlays. As shown, there is a high similarity among regions (y-axis), but clear differences among wave exposure levels (x-axis). In the legend, H stands for Hanko, T for Tv€ arminne, and J for Jussar€ o, while s stands for sheltered, ms for moderately sheltered, me for moderately exposed, and e for exposed. positive values and transect coverage of red algae as a predictor of negative values. The taxa showing the strongest influence on the observed patterns were Gammarus spp. that showed a positive relationship with wave exposure, and three species of gastropods together with Macoma balthica showing negative relationships with wave exposure (Fig. 5).

DISCUSSION
A peculiar feature of the rocky sublittoral in the Baltic Sea is that species diversity on local and landscape scales shows relatively small variations. Quality differences within and between rocky habitats in the Baltic Sea are typically expressed in dominance or abundance shifts among species, rather than in richness. The small difference in terms of biodiversity in this study, within and among regions and exposure levels, is likely a result of this omnipresent nature of species and an overall small species pool. The large-scale environmental variables salinity and temperature appear to cause region-wide effects, displayed as distinct groupings of samples from the different regions in the sample cloud of the dbRDA, with regions falling in a logical longitudinal order (Fig. 4). However, this distinctive pattern only characterizes samples from 8 m depth, whereas samples from 3 m depth do not reflect any distinct patterns driven by gradients in salinity or temperature. The effects of Mytilus are also clearly displayed in the dbRDA cloud (Fig. 4), where samples with lowest Mytilus biomass at Hanko and especially samples with highest biomass at Tv€ arminne (Fig. 2) are grouped apart from the main regional cloud, characterized by significantly more (Hanko) or significantly less (Tv€ arminne) Mytilus biomass. These findings from 8 m depth support Hypothesis 1predicting community differences driven by the interplay of regional-scale and local-scale gradients. Samples from 3 m depth show patterns that are driven by predominantly local drivers, suggesting that the local wave exposure gradient completely overrides the imprint of regional factors.
This paper demonstrates that faunal composition and faunal abundance in our study area show distinct differences that are (1) influenced by large-scale abiotic factors including salinity and temperature (supporting Hypothesis 1), (2) locally modified by wave exposure (supporting Hypothesis 2), and (3) driven by the facilitating function of Mytilus (supporting Hypothesis 2) while being less affected by Fucus (supporting weakly Hypothesis 1). With regard to the stressgradient hypothesis (Bertness and Callaway 1994), Mytilus at 8 m depth clearly affects the associated community and increases the total abundance of associated fauna with stronger effects at exposed stations compared to sheltered ones (supporting Hypothesis 2; Fig. 3A). However, at the higher wave exposure stress at 3 m compared to 8 m depth, the pattern is less obvious and indicates that the facilitation function, in fact, may be strongest at the most sheltered stations (see also D ıaz et al. 2015). The challenge in this context is therefore to develop an understanding of how the environment modulates the community structure differently at the two depths and how wave exposure changes the interaction between the invertebrate community and the two habitat-forming species. The difference in community response between depths and among exposures may arise from two interrelated factors: (1) those of how increasing physical stress changes the intensity or function of the facilitators affecting their facilitating value, and (2) those arising from depth-dependent differences in the dominance and heterogeneity of different habitats. We address these two possibilities below.
Many studies suggest a decline in facilitation toward the most severe conditions, as facilitators are less successful in ameliorating stress in very harsh environments (Kawai and Tokeshi 2007; but see He andBertness 2014, Qi et al. 2018). Therefore, while the facilitation of mussels at 8 m depth increases with exposure throughout the wave exposure gradient (Fig. 3A), mussels are unable to ameliorate stress at the more stressful physical conditions that are prevailing at 3 m depth at the exposed stations (Fig. 3B). This effect may arise from the interacting effects of mussel bed structure and wave force strength that differs among stations and between depths within stations. Such effects may result from environmental effects on the facilitating traits of the FS (Irving and Bertness 2009, Bishop et al. 2013, Bateman and Bishop 2017. At 3 m depth, mussel beds are mainly composed of many small individuals that are tightly packed and have a less complex bed structure but more patchy distributions (D ıaz et al. 2015) than at 8 m (see the ratio between abundance and biomass at different exposures in Fig. 2A, B). This difference in mussel bed structure is driven by wave-induced dislodgement that increases isometrically with increasing size of mussels and affects beds at 3 m depth significantly more than at 8 m depth (Westerbom and Jattu 2006). Such loss of architectural complexity provided by habitat-forming mussels and oysters has been shown to reduce biodiversity, species richness, or evenness (Kostylev 1996, Gutierrez et al. 2003, Kimbro and Grosholz 2006, Koivisto and Westerbom 2010. Results here suggest that the bed structure of Mytilus at 8 m depth offers a more complex three-dimensional architecture throughout the gradient, whereas a similar structure at 3 m is only found at sheltered stations where the wave force is not sufficiently strong to remove the largest individuals (Westerbom and Jattu 2006).
An alternative explanation may arise from the difference in physical structure between depths, as bed surface area may be a stronger predictor for individual abundance and species richness than environmental stress (Kostylev et al. 2005, Kimbro andGrosholz 2006). Large mussels increase the available surface area for fouling and hiding organisms, and this surface area may be largest at low-disturbance stations where mussels are larger. Additionally, the degree of shell packing and the size of mussels affect refuge provision (Kostylev et al. 1997, Commito andRusignuolo 2000) and affect further the enrichment of organic material and nutrients, which typically are limiting resources on wave-beaten shores (Tokeshi and Romero 1995, Norling and Kautsky 2008, Kotta et al. 2009). FS may have a particularly large influence on associated communities where they increase the availability of limiting resources (Bateman and Bishop 2017). At highly wave-beaten shores at 3 m, the availability of these mussel matrix resources wanes directly through increased flushing and indirectly through packing the mussel beds and thereby reducing their sediment storage potential.
A third explanation may arise from depthdependent differences in the dominance structure of habitats. Mussels below the algal belt are primary space holders that almost solely provide habitats and resources for secondary species. In the windswept Patagonian intertidal, Silliman et al. (2011) showed that competition effects by mussels become unimportant as mussels are the sole dominant space holders capable of buffering against lethal conditions, and therefore, the whole community is completely dependent on mussels for its persistence. Similarly, Norling and Kautsky (2008) and Koivisto and Westerbom (2010) showed experimentally that mussel beds host substantially more species and individuals than non-mussel habitats in the Baltic Sea. McAfee et al. (2016) showed among Australian oysters that the effects of facilitators may be dependent on the extent to which facilitators replicate the function of other co-occurring FS. Hence, at 8 m depth where there are no alternative habitat-forming biogenic structures, species richness and individual density may be particularly dependent on the availability of more homogeneous mussel beds (D ıaz et al. 2015). However, at shallow depths, mussel beds are not the sole biogenic habitat providers (Kraufvelin and Salovius 2004), likely deflating the importance of one habitat provider (Saarinen et al. 2018). In a similar way, two habitat-forming mussel species coexist and facilitate each other in the mid-zone of rocky shores on the South African south coast, while on the lower and upper zones, either species dominates due to contrasting species-specific tolerances to wave and heat stress (Erlandsson et al. 2011).
A corresponding mechanism driven by the amelioration of wave exposure was not found for Fucus. Fucus explained only a fraction of the community structure (Appendix S1: Table 4). The community effects of canopy-forming species have been widely studied in the intertidal where canopies provide shelter from temperature stress (Watt and Scrosati 2013). Under those conditions, fucoids have shown positive effects on understory organisms in physiologically harsh and physically benign environments by providing shelter from heat, whereas negative effects have been shown in physiologically benign and physically harsh environments as algal fronds smash and scour their immediate surrounding due to waves (Beermann et al. 2013 with references therein). In the shallow subtidal, disturbance caused by the thallus movement can be very pronounced and accentuate with higher wave exposure. Kiirikki (1996), for example, showed that Fucus controls the volume of filamentous algae on exposed shores, but on sheltered shores, its effect was negligible. Algal smashing and scouring may therefore reduce total individual abundance directly, but also indirectly by reducing Mytilus coverage or the size of mussel individuals in the patches below the Fucus (Appendix S6: Fig. S1). One foundation species may therefore affect the traits and the facilitating function of the other (Bishop et al. 2013). In addition, many other plausible explanations may clarify the observed patterns. For instance, when multiple stressors work in the same direction, the predictability may suffer, as species may be tolerant to one stressor and intolerant to another (Qi et al. 2018). Moreover, facilitation as a mechanism, as other biotic interactions, is primarily species-specific and responses may vary from negative to positive among species due to their different ecological requirements (Soliveres et al. 2015).
The species community and diversity patterns observed in this study are much in line with previous research efforts in the area. These show weak effects of Mytilus on species richness, but large effects on invertebrate abundance and community structure (Norling and Kautsky 2008, Koivisto and Westerbom 2010, Koivisto et al. 2011, as well as greater facilitation effects at 8 m than at shallower depths (D ıaz et al. 2015). The overall weak facilitator effect of Fucus across different wave exposures is surprising, as Fucus traditionally has been considered as an important FS in the Baltic Sea (R aberg and Kautsky 2007, Wikstr€ om and and has experimentally been shown to increase diversity and individual density (Koivisto and Westerbom 2010). Nevertheless, previous studies have been rather inconclusive when comparing Fucus with ephemeral algal habitats. Some have shown higher density of the associated community in Fucus habitats, but similar diversity Salovius 2004, Wikstr€ om and, while others have shown lower density but higher diversity (Saarinen et al. 2018). These inconclusive results suggest that the positive community effects of Fucus are highly contextdependent.
Even if the influence of wave exposure on Baltic Sea communities has been recognized for long (Kautsky and Van der Maarel 1990) and is also supported by an ample global literature (McQuaid and Branch 1984, Menge and Sutherland 1987, Harley and Helmuth 2003, many studies in the Baltic Sea still overlook the impact of wave exposure as a determinant of community composition. We show that apart from influencing on community composition, wave exposure seems to affect the interactions between species. Therefore, studies focusing on the habitat-forming function of aquatic organisms need to consider the influence of wave exposure that affects the community not only directly, but also indirectly through changing the facilitating traits of the habitat-forming species. Failure to integrate site-specific interaction reduces our ability to predict community responses to, for example, a changing climate and falls short in identifying the generality of the mechanisms behind the factors that affect the structure of communities (Gilman et al. 2010). Overall, we need a better understanding of the mechanisms that affect the relationships between foundation species, the environment within which they interact, and the shared effects these variables impose on the associated community. The challenge for scientists is to develop accurate predictive models about how these factors interplay, especially as global seas are rapidly changing due to a shifting climate. As other seas, the Baltic Sea shows rising temperatures, declining salinity conditions, and possibly an altered wave climate. Changes in these factors are predicted to shift the distribution and structure of the focal foundation species of this study (Vuorinen et al. 2015, Westerbom et al. 2019. Such structural changes among FS together with an altered sea climate (Westerbom et al. 2019) will likely affect the resident community (see also Westerbom 2010, D ıaz et al. 2015). Consideration of the structure and impact of FS along spatial and temporal scales is therefore critical for understanding the underlying mechanisms affecting species distribution and for finding management solutions to mitigate biodiversity loss (Hewitt et al. 2017). Using a hierarchical sampling design, we have gained insights to the ways positive species interactions vary with environmental conditions and have empirically shown that facilitation contributes to community structure in a sparsely studied community at the range border. However, using a mensurative approach, this study was not able to provide any mechanistic explanations for the processes underlying the observed patterns. We underline the need for manipulative experiments that can test mechanistic relationships to further our understanding of species distribution patterns that are vital for finding efficient management solutions and conservation strategies.