Synergistic Effects of Rooted Aquatic Vegetation and Drift Wrack on Ecosystem Multifunctionality

Ecosystem multifunctionality is an increasingly popular concept used to approximate multifaceted ecosystem functioning, which in turn may help advance ecosystem-based management. However, while experimental studies have shown a positive effect of diversity on multifunctionality, observational studies from natural systems—particularly aquatic—are scarce. Here, we tested the relative importance of species richness and cover of rooted aquatic vegetation, as well as cover of the loose-lying form of the macroalgae bladderwrack (Fucus vesiculosus), for ecosystem multifunctionality in shallow bays along the western Baltic Sea coast. We estimated multifunctionality based on four indicators of functions that support ecosystem services: recruitment of large predatory fish, grazer biomass, inverted ‘nuisance’ algal biomass, and water clarity. Piecewise path analysis showed that multifunctionality was driven by high cover of rooted aquatic vegetation and bladderwrack, particularly when the two co-occurred. This synergistic effect was nearly three times as strong as a negative effect of land-derived nitrogen loading. Species richness of aquatic vegetation indirectly benefitted multifunctionality by increasing vegetation cover. Meanwhile, high bladderwrack cover tended to decrease vegetation species richness, indicating that bladderwrack has both positive and negative effects on multifunctionality. We conclude that managing for dense and diverse vegetation assemblages may mitigate effects of anthropogenic pressures (for example, eutrophication) and support healthy coastal ecosystems that provide a range of benefits. To balance the exploitation of coastal ecosystems and maintain their multiple processes and services, management therefore needs to go beyond estimation of vegetation cover and consider the diversity and functional types of aquatic vegetation.

tionality was driven by high cover of rooted aquatic vegetation and bladderwrack, particularly when the two co-occurred. This synergistic effect was nearly three times as strong as a negative effect of land-derived nitrogen loading. Species richness of aquatic vegetation indirectly benefitted multifunctionality by increasing vegetation cover. Meanwhile, high bladderwrack cover tended to decrease vegetation species richness, indicating that bladderwrack has both positive and negative effects on multifunctionality. We conclude that managing for dense and diverse vegetation assemblages may mitigate effects of anthropogenic pressures (for example, eutrophication) and support healthy coastal ecosystems that provide a range of benefits. To balance the exploitation of coastal ecosystems and maintain their multiple processes and services, management therefore needs to go beyond estimation of vegetation cover and consider the diversity and functional types of aquatic vegetation.

INTRODUCTION
The management of natural systems is challenged by the global erosion of biodiversity, which threatens ecosystem processes and services on a planetary scale (Steffen and others 2015). These changes have led to a demand for a more holistic management perspective grounded in integrated ecosystem assessments (Levin and others 2009). An increasingly common approach to quantify the integrated functioning of ecosystems is the estimation of 'ecosystem multifunctionality' (hereafter 'multifunctionality') or the simultaneous provisioning of multiple ecosystem processes (Hector and Bagchi 2007). Multifunctionality estimations can be based on any ecosystem processes but have mainly been used to quantify the simultaneous performance of multiple ecosystem functions (that is, processes that support ecosystem services) and/ or ecosystem services (Manning and others 2018;Garland and others 2020). Consequently, while multifunctionality is rarely considered in natural systems today, its specific objectives may offer important advances in ecosystem management (Manning and others 2018), where the overarching goal is to balance human use of natural resources to maintain ecosystem structure and processes and, indirectly, the services they support.
All species differ in their contribution to different ecosystem processes; therefore, species richness has often been found to correlate positively with multifunctionality (Hector and Bagchi 2007;Gamfeldt and others 2008;Zavaleta and others 2010;Isbell and others 2011;Maestre and others 2012;Byrnes and others 2014;Soliveres and others 2016). Several studies suggest that the effect of biodiversity (often estimated as species richness) on multifunctionality is different, and possibly stronger, than the effect on single functions (Hector and Bagchi 2007;Gamfeldt and others 2008;Meyer and others 2018; but see Gamfeldt and Roger 2017). Positive biodiversity effects on multifunctionality have been found in a large number of terrestrial and aquatic experiments, but with weaker effects among primary producers than among herbivores (Lefcheck and others 2015). Yet, there is limited knowledge of how pervasive the biodiversitymultifunctionality relationship is in natural systems, particularly in aquatic environments (Soliveres and others 2016).
Aquatic vegetation (that is, rooted plants and macroalgae) can be foundation species in shallow ecosystems and sustain a wide array of ecosystem functions that support ecosystem services, such as sediment stabilization (Madsen and others 2001;Scheffer 2004;Austin and others 2017), nutrient and carbon storage (McGlathery and others 2007;Wang and others 2016), and primary and secondary production (Duarte 2000). In many parts of the world, however, aquatic vegetation species are threatened by near-shore development, eutrophication, and cascading effects of overfishing (Lotze and others 2006;Orth and others 2006;Waycott and others 2009). Long-term efforts to reduce nutrient loads have led to partial or full recovery of aquatic plants both in some estuarine (for example, Lefcheck and others 2018) and shallow freshwater ecosystems (Sand-Jensen and others 2008), showing that concentrated management actions may reverse these negative trends. In Chesapeake Bay (USA), the recovery of aquatic plants was greatest in the less saline areas where more than a dozen aquatic plant species coexist and their species richness promoted increased plant cover (Lefcheck and others 2018). High species richness of aquatic vegetation has also been found to benefit seagrass restoration (Williams and others 2017), as well as the total cover and temporal stability of seaweed assemblages (Stachowicz and others 2008).
Drift algae are detached or naturally free-living macroalgae that are common in shallow coastal ecosystems (Valiela and others 1997). The term is often used to for ephemeral, filamentous taxa favored by eutrophication (Fletcher 1996), which are often seen as a nuisance because they (1) house lower abundance of fish in comparison with other types of vegetation (Deegan and others 2002), (2) compete for light and/or space with attached and habitat-forming vegetation like perennial plants or attached macroalgae, and (3) accumulate and decompose on beaches and in shallow areas, contributing to low-oxygen conditions (Hauxwell and others 2001;Berglund and others 2003). In northern Europe, common taxa that can occur as drift algae are fucoids (that is, Fucus spp.) commonly referred to as 'wrack.' One example is Fucus vesiculosus, a brown seaweed normally growing attached to rocks, shells or bedrock in wave-exposed areas. In the Baltic Sea (N Europe), F. vesiculosus can Aquatic Vegetation and Multifunctionality also be found growing loose-lying in sheltered sediment areas and there sometimes form dense mats. This particular growth form of F. vesiculosus was first described in 1901 (Svedelius 1901), and was in 2013 designated as a special biotope for protection (HELCOM 2013a). In the NE Atlantic similar types of drift wrack have been found to disturb seagrasses and seagrass restoration by shading, causing anoxic conditions and through uprooting and burial of shoots (Valdemarsen and others 2010;Moksnes and others 2018). Although the regular attached Fucus spp. provide important habitat for fish spawning, juvenile fish and their prey (Dempster and Kingsford 2004;Wikströ m and Kautsky 2007;Snickars and others 2010), knowledge of the functional role of the drift wrack is currently lacking.
The aim of this study was to disentangle the relative importance of species richness and cover of rooted aquatic vegetation, cover of drift wrack (Fucus vesiculosus, hereafter 'bladderwrack') and abiotic conditions, on ecosystem multifunctionality in a natural, coastal ecosystem. We conducted a field survey in 32 shallow coastal bays in the western Baltic Sea along natural (wave-exposure, bay openness and size) and anthropogenic gradients (land-derived nutrient load) known to strongly influence shallow bay ecosystems (Hansen and others 2008;Snickars and others 2009;Rosqvist and others 2010;Donadi and others 2017). In August and early September, right after ecosystem productivity peaks (Wijnbladh and Plantman 2006), we measured four indicators of key ecosystem functions known to be influenced by aquatic vegetation-recruitment of large predatory fish (as abundance of young-of-the-year perch Perca fluviatilis and pike Esox lucius), algal grazer biomass (biomass of invertebrate mesograzers), 'nuisance' algal biomass (inverted filamentous algal biomass per vegetation biomass), and 'water clarity' (inverted turbidity)-and used these to estimate multifunctionality. Previous studies show that (1) the abundance of juvenile pike increases with cover of (especially rooted) vegetation, while abundance of juvenile perch peaks at 40-80% cover (Hansen and others 2019); (2) the biomass of algal-feeding grazers increases with vegetation cover and in turn suppresses filamentous algae  (Agawin and Duarte 2002).
Based on theories of the influence of biodiversity on multifunctionality as well as system-specific drivers and responses (see above), we hypothesized that multifunctionality in shallow coastal bays is positively influenced by high diversity and cover of aquatic vegetation (which are in turn influenced by environmental factors like the size and topographic openness of the bays), and therefore may counteract a negative influence of nutrient loading (eutrophication). Moreover, because rooted vegetation and drift wrack (free-living Fucus vesiculosus) differ greatly in traits like morphology, their cooccurrence might benefit multifunctionality through a complementarity effect, evident as a synergistic interaction. To tease apart these direct and indirect drivers of multifunctionality we used piecewise path analysis, a form of structural equation modeling (SEM) (Shipley 2000;Grace and others 2010). Following standard procedure in SEM we outline the specified interactions ('paths') and their prior support in the Methods section. Finally, we also tested the robustness of the overall results for changes in the relative contribution of each of the four 'processes' to multifunctionality, using a weighted averages approach.

Study System
The Baltic Sea is a large (377,000 km 2 ) brackish marginal sea in northern Europe, with some of the largest archipelagos in the world with thousands of islands and shallow, wave-sheltered bays (Figure 1). The bays often house abundant and diverse assemblages of aquatic vegetation and macroinvertebrates of both marine and freshwater origin (Hansen and others 2008;Rosqvist and others 2010). The bays are also important spawning and nursery grounds for coastal fish due to rapid warming of the water during spring, where several fish species lay their eggs and use the aquatic vegetation for foraging and shelter (Snickars and others 2009;Hansen and others 2019). However, the Baltic Sea and its archipelagos are highly affected by anthropogenic activities such as eutrophication (Elmgren 1989;Gustafsson and others 2012) and coastal exploitation (Sundblad and Bergströ m 2014; Hansen and others 2019), which is why these vegetated habitats are threatened and targets for conservation (HELCOM 2013b).
The vegetation community in the shallow wavesheltered (soft-bottom) bays of the Baltic Sea consists of a total species pool of ca 15-40 plants and perennial, non-filamentous macroalgae (depending on region, Hansen 2016), of which pondweeds (Potamogeton/Stuckenia), milfoils (Myriophyllum), stoneworts (Chara) and (in some areas) dense mats of loose-lying Fucus vesiculosus dominate. Moreover, many bays harbor high abundances of filamentous algae, mainly as attached epiphytes or loose-lying, with Cladophora glomerata dominating (Hansen and others 2008).

Field Survey
The field survey was performed in August and early September 2014 in 32 shallow bays along 360 km of the central Swedish Baltic Sea coast ( Figure 2). The bays were chosen to form gradients in topographic openness and land-derived nutrient load, and were situated more than 10 km apart (water distance) to ensure that fish communities were independent based on connectivity and migration (Saulamo and Neuman 2002; Berkströ m and others 2019).
The water surface area (a, m 2 ) of each bay was estimated from satellite images in Google Earth Pro (Version 7.1.5.1557), and ranged from 0.4 to 25 ha (mean: 5.5 ha). Topographic openness (Ea, dimensionless) was estimated as Ea = 100 * At/a, where At is the smallest cross-sectional area (m 2 ) connected to the sea (calculated from field mea-surements of inlet depth and satellite image-derived width) (Persson and others 1994). Modeled data of land-derived nutrient load (kg year -1 ) from watersheds surrounding each bay were derived from the Swedish Meteorological and Hydrological Institute's database 'SVAR' [version 2012_2] and 'S-HYPE' model [version 2012_1_2_1] (Arheimer and others 2012).
Within each bay, sampling was conducted at 6-9 stations, with higher replication in larger bays. To gain a representative sampling of the whole bay area, the stations were distributed at least 30 m apart across the bay at 0.5-3.5 m depth. The stations were marked with weighted buoys. We first sampled water for turbidity. Water was sampled at 0.5 m depth, stored dark and cold for ca. 6 h until turbidity was measured using a handheld turbidimeter (AquafluorÒ, Turner Designs, USA). The average turbidity (NTU) was calculated from 3 measurements per station. After the water sampling we took a break on land to reduce disturbance on fish. Second, we sampled young-of-theyear (YOY) fish using small underwater detonations (non-electric 1 g charge igniting a 10 g primer), which kills or stuns fish within a ca. 5 m radius (ca. 80 m 2 ) (Snickars and others 2007; Sundblad and others 2014). The point of the detonation was directly marked with a weighted buoy, and floating fish were then netted and sunken fish collected by a snorkeler. YOY perch (Perca fluviatilis) and pike (Esox lucius) were identified and counted, Third, the cover of benthic vegetation was surveyed within 5 m of the buoy marking the detonation. The cumulative % cover of all macroscopic vegetation species was visually estimated by a snorkeler within five randomly placed 0.5 9 0.5 m frames (0.25 m 2 , Figure 3) (for a more detailed description, see Austin and others 2017). The rooted vegetation was dominated by the angiosperms Stuckenia pectinata, Potamogeton perfoliatus and Myriophyllum spicatum (Table S1). We here also classified Characeans (Chara baltica being the most common) as rooted, because of their similar growth form and root-like rhizoids. Adding bladderwrack, this classification captured 23 of 27 habitat-forming species occurring in the 32 bays (Table S1). Nearly all bladderwrack was loose-lying (for example, 80% of observations were in areas with pure sediment bottom, and 95% in areas with sedimentdominated bottoms). Species richness of the rooted vegetation was estimated as a-diversity, that is, the mean rooted species richness per station for each bay.
Finally we sampled vegetation including filamentous algae (mainly epiphytic) and associated macroinvertebrates within a randomly placed 0.2 9 0.2 m frame (0.04 m 2 ), with a 1 mm-mesh bag attached. In the laboratory, 3-8 samples from each bay were analyzed. First, filamentous algae and macroinvertebrates were separated from the vegetation. The mean dry biomass of filamentous algae and vegetation per station were determined after drying at 60°C for 48 h. The macroinvertebrates were identified, their body length measured and their biomass estimated as gram ash-free dry mass per sample (0.04m 2 ) using taxon-specific correlations between length and ash-free dry mass (Eklö f and others 2017).

Estimation of Ecosystem Multifunctionality
After averaging all station-level variables per bay (N = 32) we estimated ecosystem multifunctionality based on four indicators of ecosystem functions, that is, processes that directly and indirectly support ecosystem services in shallow coastal areas of the Baltic Sea (HELCOM 2010). These were: (1) large predatory fish recruitment quantified as the abundance of young-of-the-year (YOY) perch and pike, (2) grazer biomass (see Donadi and others 2017 Appendix D), (3) biomass of filamentous 'nuisance' algae quantified as the inverted biomass of filamentous algae per biomass of vegetation (which yields a ratio between 0-1), and (4) water clarity, quantified as the reflected turbidity (-NTU + max(NTU)). Biomass of the filamentous algae and water turbidity were inverted and reflected, respectively, because high values of 'multifunctionality' should express the desirable level of the variable. The four indicators were standardized by dividing each variable by the average of its three highest values (yielding an overall multifunctionality range of 0.25-0.99). A correlation matrix of the average multifunctionality and the four indicators is shown in Figure S1. Although multifunctionality can be estimated in several ways, we chose the 'averaging' approach (where multifunctionality is the mean of the included functions) because it is simple, intuitive and therefore easy to communicate. We refrained from using 'threshold' approaches (Byrnes and others 2014) because they are very sensitive to standardization method and the number of functions included (Gamfeldt and Roger 2017).

Path Analysis
To assess the direct and indirect drivers of multifunctionality we used piecewise path analysis (Shipley 2000  vegetation, and (3) increases the cover of drift wrack (Rosqvist and others 2010). These structural changes are a result of increased wave exposure and water exchange leading to decreased sedimentation and accumulation of fine organic particles, resulting in less nutrient-rich and more heterogeneous and coarser substrates, with increasing openness. In addition, open bays are more accessible to both drifting filamentous algae and bladderwrack than enclosed bays (Berglund and others 2003;Rosqvist and others 2010). We also hypothesized that (4) bay surface area increases the species richness of vegetation (Rosenzweig 1995), and that (5) land-derived nitrogen load has a direct negative effect on multifunctionality by increasing growth of phytoplankton (increasing turbidity) and ephemeral filamentous algae (for example, Schramm and Nienhuis 1996). We initially included a direct effect of land-derived nitrogen load on bladderwrack and rooted vegetation cover (negative eutrophication effect), but the relationship was not significant in any of the models (see also Austin and others 2017) and was therefore left out to decrease model complexity given the relatively low sample size. We focused on nitrogen (N) rather than phosphorous (P), as the ratio between the two elements suggested that most of the bays were N-limited at the time of sampling (based on measured nutrient concentrations in the bays and the Redfield ratio; see Austin and others 2017). However, the modeled load of nitrogen and phosphorous were highly correlated in our data (Spearman rho = 0.96, p < 0.0001). Further, we hypothesized that the species richness of rooted vegetation (6) (7) is negatively affected by increased cover of loose-lying bladderwrack due to shading and competition for space (Valdemarsen and others 2010; Moksnes and others 2018). We initially also included a direct effect of bladderwrack on rooted vegetation cover, but this path was not significant in any of the models and was therefore removed. However, the upper limit of rooted vegetation cover seemed to be limited by high bladderwrack cover ( Figure S2). Similar to many previous multifunctionality studies we initially also included a direct effect of species richness on multifunctionality (for example, Hector and Bagchi 2007), but the effect was not significant in any of the models and was therefore removed. Instead, we hypothesized that there was an indirect effect of vegetation species richness on multifunctionality mediated by increased vegetation cover (for example, Lefcheck and others 2018), a variable well known to affect many ecosystem processes. Finally, we hypothesized that (8) the cover of bladderwrack and rooted vegetation increases the average multifunctionality due to their habitat-forming characteristics, contribution to secondary production, ability to compete with filamentous algae and phytoplankton for nutrients and to increase sedimentation and decrease re-suspension of sediment ( Additionally, rooted vegetation species such as Myriophyllum spicatum, Potamogeton perfoliatus and Stuckenia pectinata can grow much taller (2-3 m, Figure 5) than bladderwrack ( £ 0.5 m) and form complex canopies above it. This led us to hypothesize that the effect of the two vegetation types on multifunctionality might even be synergistic, rather than simply additive. We tested if the effects of bladderwrack and rooted vegetation cover on multifunctionality were additive or synergistic by comparing two separate SEMs: one including the additive effects of bladderwrack and rooted vegetation (model 1, Figure 4) and one also including their interaction (model 2, Figure 4). The models thus included three exogenous variables: (1) topographic openness (log 10 -transformed to reduce skewness), (2) bay surface area (log-transformed to reduce skewness) and (3) land-derived nitrogen load (log 10 -transformed to reduce skewness). The endogenous variables were (1) bladderwrack cover (average per bay), (2) rooted vegetation cover (average per bay), (3) species richness of rooted vegetation (a-diversity, the mean richness per station per bay), and (4) average multifunctionality (per bay).
To better understand the underlying processes driving the multifunctionality, and to test the robustness of the results to changes in individual functions, we also estimated the drivers of weighted multifunctionality. We doubled the individual contribution of each function to the multifunctionality metric (from 25 to 50%, using weighted means) in separate path models. Models including weighted multifunctionality were then compared to the two original path models in terms of their model fit (Akaike's Information Criterion, AIC) and effect sizes (that is, standardized coefficients) of the predictors of multifunctionality.
The distributions of all variables were examined graphically before model fitting. When normality assumptions were violated, a negative binomial error distribution with a log-link function was applied. Predictors were tested for multicollinearity by calculating the variance inflation factor (VIF < 2) and the correlation coefficients for all predictor combinations (< 0.63).
Model fit was assessed using the test of directed separation that produces a v 2 -distributed Fisher's C statistic where models are indicated to fit the data when p values are larger than 0.05 (Shipley 2009). We compared the models with and without the interaction between bladderwrack and rooted vegetation using AIC. Competing models with DAIC < 4 were considered to fit the data equally well (Burnham and others 2011;Shipley 2013). To compare the relative strength of the significant paths, standardized coefficients were calculated using the latent-theoretical method developed for binomial response models (Grace and others 2018), which was adjusted for negative binomial distribution with log-link function by log-transforming the denominator (sd(ln(Y + 1)) (JB. Grace, pers. comm. in April 2020). Models were visually validated by plotting the residuals against the fitted values.

Results
Multifunctionality was best described by a combination of direct and indirect effects of bay morphology, species richness of vegetation, vegetation cover and nutrient loading. The model including the interaction between bladderwrack and rooted vegetation cover (model 2) fitted the data as well as the additive model as shown by the test of directed separation (model 1, see Table 1), but explained a much higher proportion of the variability in multifunctionality (R 2 = 0.56 vs. 0.30, Table 1). In model 2, there was a positive effect of the interaction between bladderwrack and rooted vegetation cover on multifunctionality (p < 0.001). The interaction was a result of a synergistic effect (Figure 6): multifunctionality increased more steeply with higher cover of both rooted vegetation and bladderwrack together, than by rooted vegetation alone (Figure 7 and Figure S3). For example, when rooted vegetation cover was 30%, multifunctionality was 0.741 when bladderwrack cover was 10%, but 0.995 when bladderwrack cover was 40%. The interaction was the second strongest effect in the two models (standardized coefficient: 1.18), following a positive effect of bay area on species richness of rooted vegetation (standardized coefficient: 1.45). There was also a marginally significant negative effect of land-derived nitrogen load on multifunctionality ( Figure 6, p = 0.055). Bladderwrack cover increased with topographic openness, while rooted vegetation cover increased with species richness of rooted vegetation. There was a trend toward a negative effect of bladderwrack cover on species richness of rooted vegetation ( Figure 6, p = 0.072), and a trend toward a negative effect of topographic openness on rooted vegetation cover (Figure 6, p = 0.073), but there was no statistically significant effect of topographic openness on species richness of rooted vegetation.
In the additive model (model 1 effects) multifunctionality increased with cover of bladderwrack and rooted vegetation (Table 1, Figure S4 and Figure S5). Here, the negative effect of land-derived nitrogen load on multifunctionality was also marginally significant (p = 0.051) and comparable in strength to the positive effects of bladderwrack and rooted vegetation cover (Table 1, Figure S4 and Figure S5). The remaining paths did not differ from model 2.
The weighted multifunctionality models largely confirmed the results of the original (unweighted) model. Accordingly, also in the weighted analyses a higher proportion of the multifunctionality was explained in models including the interaction, than in the additive models (Table 1), and models did not differ in terms of AIC (D AIC < 4 units). In the weighted models with higher influence of large predatory fish recruitment, or 'nuisance' algal biomass, the interactive effect of bladderwrack and rooted vegetation cover on multifunctionality was slightly stronger than in the original model (standardized coefficient: 1.27 and 1.30 vs. 1.18). When instead the influence of water clarity was increased, the additive effects of bladderwrack and rooted vegetation cover on multifunctionality in model 1 increased in strength. Finally, when the influence of water clarity or grazer biomass were increased, the negative effect of land-derived nitrogen load on multifunctionality became significant (compared to marginally significant), both in model 1 and 2. For a more detailed description of the results from the weighted multifunctionality models see Text S1.

Discussion
The aim of this study was to assess the relative importance of species richness and cover of rooted aquatic vegetation, together with cover of drift wrack, for ecosystem multifunctionality in temperate shallow coastal bays. We found that high cover of rooted vegetation and high cover of bladderwrack (mainly loose-lying and drifting) had a synergistic effect on multifunctionality. This synergy appeared especially important for large predatory fish recruitment and the filamentous algal biomass, whereas water clarity was more influenced by the additive effects of bladderwrack and rooted vegetation cover. The rooted vegetation cover was driven by the species richness of vegetation, while bladderwrack cover was driven by bay openness (wave exposure and water exchange). In contrast to the positive vegetation effects on multifunctionality, there was a direct negative effect of land-derived nitrogen on multifunctionality; an effect that increased in strength when water clarity or grazer biomass was given more importance to multifunctionality.
Our results suggest that the species richness of the vegetation had an indirect positive effect on multifunctionality, by increasing the cumulative  Aquatic Vegetation and Multifunctionality vegetation cover. Plant species richness has been experimentally shown to increase plant cover and biomass in hundreds of experiments (reviewed by Cardinale and others 2011). The positive effect of species richness on rooted vegetation cover in our study, combined with a lack of a statistically significant direct effect of richness on multifunctionality, indicates that rare species are important for achieving high cover, even though they do not contribute directly to multifunctionality (see also Gamfeldt and Roger 2017). Plant cover is a key composite plant trait often linked to ecosystem structure and function in coastal ecosystems (for example, Pillay and others 2010). Species richness has often been found to predict multifunctionality (Lefcheck and others 2015), but had no direct effect in this study (see also review by van der Plas 2019). The lack of a direct diversity effect on multifunctionality could be caused by a few dominant species sustaining high multifunctionality, whereas a high vegetation diversity allows more complete use of the substrate and therefore increases vegetation cover. The three most common rooted vegetation species in our study were Stuckenia pectinata, Potamogeton perfoliatus, and Myriophyllum spicatum (occurring in 29, 24 and 22 of the 32 bays, respectively, Table S1). They are all dominant, canopy-forming and competitive species with shoots reaching up to 3 m from the seabed (unpublished data), forming structurally complex underwater 'jungles.' Viewed together, these results suggest that both vegetation diversity and identity benefit the functioning of these shallow, coastal ecosystems. High cover of loose-lying and drifting bladderwrack also increased multifunctionality, particularly when combined with high cover of rooted vegetation (interaction effect; see model 2, Figure 6). This synergistic effect suggests that wrack and rooted vegetation complement each other functionally, potentially because high habitat volume, surface area, complexity and height cannot be achieved by high wrack or high rooted vegetation cover alone. But when the two vegetation types co-occur, a more complex and (locally) variable habitat is created, allowing different species and life stages to find refuge, spawn or feed. The importance of the synergistic effect on associated consumers is supported by our findings that the synergy increased in strength when increasing the contribution of either large predatory fish recruitment and biomass of filamentous algae to multifunctionality (Table 1). Further, when increasing the contribution of grazer biomass, the additive effects of bladderwrack and rooted vege- Table 1. tation cover on multifunctionality disappeared, while their synergistic (interactive) effect remained. This reflects the importance of the highly complex and variable habitat also for grazers. In contrast, water clarity was mainly explained by the additive effects of high rooted vegetation and bladderwrack cover (Table 1). This is in line with previous studies showing the importance of total vegetation cover for increased competition with phytoplankton, increased sedimentation and decreased re-suspension (Austin and others 2017). In summary, this plant-wrack synergy effect constitutes a form of higher-order, macrophyte functional diversity effect that complements the aforementioned indirect positive influence of rooted vegetation diversity on multifunctionality.
Besides the direct positive effect of bladderwrack cover on multifunctionality, there was also a trend toward a negative effect of bladderwrack cover on species richness and cover of rooted vegetation, which would indirectly reduce multifunctionality. We speculate that while drift wrack and rooted vegetation complement each other functionally (see above), they can also compete for limiting resources such as light and space. Previous studies suggest there may be a threshold in rooted vege- Figure 6. Path diagram of model 2 with the interaction between bladderwrack and rooted vegetation cover. The width of the arrows is proportional to the standardized coefficients (also shown in numbers next to arrows). Significance level was set at alpha = 0.05. Nonsignificant relationships are indicated with dashed arrows. The positive effect of rooted vegetation cover increases with increased cover of bladderwrack and vice versa. This interactive effect is also shown in Figure S3, together with the data points. tation density, beyond which the rooted vegetation is dense or tall enough to survive and continue to grow in the presence of drifting wrack (Moksnes and others 2018). However, small and/or sensitive species (such as Chara aspera, Chara globularis, and Ruppia cirrhosa), or sparsely growing meadows are likely to be more severely affected by drift wrack. This could explain the near-significant trend of a decrease in rooted vegetation species richness with increasing bladderwrack cover, instead of a direct effect of bladderwrack on rooted vegetation cover. Drift wrack has previously been found to hinder seagrass restoration by outcompeting underlying seagrass shoots, but if the seagrass instead grows densely, the wrack stays at the edges of the meadow (Valdemarsen and others 2010; Moksnes and others 2018). Our findings indicate similar competitive interactions in shallow Baltic Sea bays, which probably manifest in spring when the rooted vegetation is generally sparser and smaller. However, the direct positive influence of drift wrack on multifunctionality suggests that drift wrack may compensate for its indirect negative effect on multifunctionality, by in itself forming an important habitat. Consequently, the relationship between drift wrack and ecosystem processes and services can be much more complex than the often-observed negative impacts of more common mats of filamentous drift algae (Hauxwell and others 2001). Additionally, the abundance of drift wrack in our study is more likely driven by water movement alone (reflected by topographic openness), in contrast to the filamentous algae that are influenced by both wave action and high nutrient loading (Fletcher 1996;Berglund and others 2003). Hence, reducing eutrophication should in theory increase the amount of drift bladderwrack through increased light penetration, as opposed to reducing drifting nuisance algae. However, the persistence of drift wrack beds, particularly in relation to rooted vegetation, must be studied in greater detail.
The direct negative effect of land-derived nutrient load on multifunctionality suggests that eutrophication has negative effects not just on single ecosystem components, but on the entire ecosystem functioning. Simultaneously, this finding indicates that measures that reduce eutrophication may increase multifunctionality. Regulation of fisheries may have similar effects on fish communities as eutrophication mitigation measures, as shown in the Baltic Sea where fishing targeting large predatory fish has been found to generate eutrophication-like effects (Eriksson and others 2011; Bergströ m and others 2019). Moreover, protecting rooted vegetation against effects of, for example, dredging and recreational boating (Short and Wyllie-Eciieverria 1996; Hansen and others 2019) is important to maintain or increase multifunctionality, since high species richness and cover of vegetation counteracts the negative effect of nitrogen load on multifunctionality.
In conclusion, this study links ecological theory with management applications by illustrating how the diversity and abundance of habitat-forming aquatic vegetation jointly affect ecosystem multifunctionality. We estimated multifunctionality based on four functions tightly linked to ecosystem services such as food provisioning, recreation, maintenance of biodiversity and habitats, and biological regulation in coastal areas of the Baltic Sea (HELCOM 2010). These ecosystems have been classified as having bad or moderate status in the study area, where excess nutrient loading and commercial fishing are anthropogenic drivers with large negative net effects on habitats and primary production others 2015, 2020). Other anthropogenic pressures in these shallow coastal ecosystems include recreational boating, mooring structures and shoreline development Sagerman and others 2020). To balance the human use and exploitation of coastal ecosystems and the maintenance of ecosystem services, we need a better understanding of the contribution of coastal habitats to ecosystem services in environmental decision making, where the use of an ecosystem-based approach is increasing (Kosenius and Ollikainen 2015). Many people that live near the coast value habitatforming vegetation, conservation of pristine areas and large stocks of predatory fish (Kosenius and Ollikainen 2015), whereas others find dense and high submerged vegetation to be a nuisance that clogs waterways (van Nes and others 2002, Sandströ m and others 2005). Our results indirectly suggest that protecting high vegetation diversity may counteract the negative effects on multifunctionality imposed by anthropogenic pressures such as increasing nutrient load. The positive effects of the diversity-driven rooted vegetation cover and drift wrack on ecosystem multifunctionality were of equal or greater strength (additive vs. synergistic effects) than the negative effect of nutrient load, indicating the importance of aquatic vegetation for healthy and highly valuable coastal ecosystems.

FUNDING
Open access funding was provided by Stockholm University. This study is a product of project Plant-Fish and was funded by Formas (grant 2013-1074), HM Carl XVI Gustaf's Foundation for Science and Education (2014-0002), the Baltic Sea 2020 foundation and the Baltic Sea Centre (Askö grants).

DATA ACCESSIBILITY
The data used for the manuscript are included as electronic supplemental material including complete metadata (Appendix 1).

Co mpliance with E thical St andards
Et hi cal Approval The fish sampling procedures were judged and approved by the Ethical Board on Animal Experiments of the County Court of Uppsala, Sweden, permit C 139/13 (to the Department of Aquatic Resources, Swedish University of Agricultural Sciences).
Conflict of interest The authors declare that they have no conflict of interest.

OP E N A C C ES S
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 h ttp://creativecommons.org/licenses/by/4.0/.