The interplay of co-occurring ecosystem engineers shapes the structure of benthic communities – a mesocosm experiment

1 Introduction

Ecosystem engineers are essential organisms that transform the environment they inhabit by creating new habitats and modifying existing ones (Jones et al., 1994). Autogenic ecosystem engineers, also known as habitat-forming organisms, alter their environment through their own physical structures (Jones et al., 1994). On the other hand, allogenic ecosystem engineers modify the environment by transforming living or non-living materials from one physical state to another (Jones et al., 1994). These engineers play a crucial role in shaping benthic communities and in influencing the functioning of entire ecosystems, especially in marine sediments (Van Der Heide et al., 2012; Rossi et al., 2013). For instance, seagrass meadows provide habitats for a range of marine organisms (Gartner et al., 2013). Meanwhile, burrowing animals, such as polychaetes and bivalves, create oxygenated burrows that improve microbial activity and nutrient cycling (Rosenberg, 2001; Laverock et al., 2011). When ecosystem engineers are lost or decline in numbers, potential cascading effects on other species and ecosystem processes may ultimately reduce the resilience and functioning of marine ecosystems (Waycott et al., 2009; Rossi et al., 2013).

In marine sediments, the presence of habitat-forming organisms such as seagrasses and bivalves significantly affects the settlement of invertebrate larvae, ultimately shaping sediment communities (Gartner et al., 2013; Eklöf et al., 2015). The mechanisms that govern these effects are diverse, including physical and chemical changes. The habitat structure created by ecosystem engineers can alter local hydrodynamics, leading to increased particle retention, which in turn facilitates the passive settlement of invertebrate larvae (Boström and Bonsdorff, 2000). Moreover, the chemical changes induced by habitat-forming organisms act as chemical cues for pelagic larvae to settle (Welch et al., 1997). Through so-called facilitation cascades, ecosystem engineers can also have an indirect effect on recruitment by facilitating organisms that affect local recruitment before and after settling (Bishop et al., 2012). Therefore, the interactions between ecosystem engineers and their associated community can significantly influence the success of colonization by invertebrates and their subsequent development.

In the Baltic Sea, the eelgrass Zostera marina Linnaeus, 1753, the epifaunal mussel Mytilus spp. (Linnaeus, 1758) and the infaunal clam (Baltic clam) Macoma balthica (Linnaeus, 1758) are abundant ecosystem engineers which, naturally co-occur and interact (Reusch et al., 1994; Gagnon et al., 2020; Meysick et al., 2020; WoRMS Editorial Board, 2023). They play a vital role in shaping benthic communities and in the functioning of marine environments (Boström and Bonsdorff, 1997; Norling and Kautsky, 2008). Seagrass meadows offer critical habitats for diverse marine species, including shelter for epifauna and juvenile fish, and support for deposit feeders by stabilizing sedimentary bottoms and accumulating organic matter (Boström et al., 2006a). Furthermore, seagrass meadows serve as a substrate for epiphytic algae (Borowitzka et al., 2006) and provide a settling environment for bivalves (Boström and Bonsdorff, 2000). On the other hand, blue mussels (Mytilus spp.; a species complex including a hybrid between Mytilus trossulus A. Gould, 1850 and Mytilus edulis Linnaeus, 1758; Stuckas et al., 2017a, b; WoRMS Editorial Board, 2023) create complex structures by forming patches on soft bottoms, which modify sediment topography, provide habitat for invertebrates, and promote biogeochemical changes in sediment by accumulating organic matter through pseudo-feces and carbonates through shell fragments (Norling and Kautsky, 2008; van der Zee et al., 2012). Although there is limited knowledge about the direct effects of Baltic clams on the colonization and structuring of other invertebrates, they play a key role for coastal sediment dynamics and nutrient fluxes along the sediment-water interface through bioturbation and bioirrigation (Norkko et al., 2013; Meysick et al., 2022), suggesting their importance for the composition of invertebrate communities. Both, seagrasses and bivalves, are essential for maintaining the health and productivity of the benthic community, and their conservation is crucial for preserving the marine ecosystem’s functioning in the Baltic Sea.

Despite the recognized significance of ecosystem engineers in soft sediment habitats, our current understanding of how habitat structure, as well as pre-existing community structure and functioning, interact to determine benthic colonization in the Baltic Sea and other similar systems remains limited. This knowledge gap highlights the need for experimental studies assessing the role of ecosystem engineers and their interactions in shaping benthic communities in soft sediments. To address this, we used a mesocosm approach to investigate the effects of three ecosystem engineers and their combinations on invertebrate settlement. The treatments included (i) vegetated and unvegetated habitats, (ii) with or without adult epifaunal mussels, and (iii) with or without adult infaunal clams. Our mesocosm experiment was designed to test two working hypotheses: Firstly, that the community structure of newly settled invertebrates is habitat-dependent. Secondly, that the ecosystem engineering will promote biodiversity by creating distinct ecological niches. Furthermore, we assessed how the number of ecosystem engineers in a benthic habitat affects the (re)colonization of benthic invertebrates.

2 Materials and methods

2.1 Mesocosm settings

A mesocosm system was established at the Skärgårdscentrum in Korpoström, Finland (Figure 1A, Supplementary Figure S1), inspired by the Kiel Indoor (Pansch and Hiebenthal, 2019) and Outdoor (Wahl et al., 2015) Benthocosms. The system comprises 20 600-liter polyethylene tanks (Schwer-PALOXE 1130 x 930 x 600 mm) supplied with a constant flow-through of unfiltered seawater (GRUNDFOS SP 5A-12 R 5 pump installed at a depth of three meters, from the Archipelago Sea, around 20 m from the tank system). The system mimics natural conditions in nutrients and organic matter and daily fluctuation of abiotic factors such as water temperature, salinity, pH, and oxygen.

Figure 1 Experimental setting: (A) general view of the mesocosm system, encompassing 20 600-liter tanks, (B) the arrangement of the 80 experimental units (20-liter buckets, four in each tank and a header tank providing flow-through of water) during the present experiment, and (C) the treatments investigated during the present experiment, encompassing sand-only, and all crossed combinations of the bivalves Mytilus spp. (My), Macoma balthica (Ma), as well as the seagrass Zostera marina (Zo).

For this experiment, each tank contained four 20-liter experimental units (vinyl buckets PP Plastex^®; Figure 1B). Header tanks (10-liter buckets PP OrthexTM; Figure 1B) above the water level distributed water passively among the experimental units via 90 cm long silicon tubes (4 mm diameter, Versilic Silicon Tubing #760320), maintaining a flow rate of 15.7 L h^-1 (± 0.4), and replacing the water volume within each unit approximately every 1.3 hours, and about 18 times per day. The water within the experimental units remained isolated from the water bath of the main tank.

To mimic natural light conditions, two layers of 1 mm mesh-size mosquito netting were placed over the top of each mesocosm tank providing a light intensity (Li-Cor model Li-1400 data logger) of 440 UM (µmol photon m-2 s^-1) on average, representative of 3 m water depth in the area during summer.

2.2 Experimental communities

The experiment was designed to examine the impact of various combinations of dominant ecosystem engineers on the colonization of macrobenthic assemblages. Eight distinct treatments (Figure 1C) were established, encompassing: (1) a treatment without ecosystem engineers (Sand), (2) a monospecific treatment with infaunal clams (Macoma), (3) a monospecific treatment with epifaunal mussels (Mytilus), (4) a monospecific treatment with seagrass (Zostera), as well as all combinations of species presence/absence in a fully factorial design, i.e., (5) Mytilus + Macoma (My + Ma), (6) Zostera + Macoma (Zo + Ma), (7) Zostera + Mytilus (Zo + My); and (8) a treatment including all three ecosystem engineers (Zo + My + Ma).

A randomized (incomplete) block design was employed to arrange the treatments across the 20 mesocosms, with each mesocosm containing four treatment combinations (Figure 1C; Supplementary Figure S2). Consequently, each treatment was replicated ten times, resulting in a total of 80 independent experimental units.

Sediment, M. balthica, and Z. marina were collected from Ängsö Bay (60°06’31.3”N|21°42’42.8”E) as the bay features a patchy seagrass meadow, in which M. balthica and Mytilus spp. are abundant. The sediment in the area consists of organic-poor (<0.5%) fine-grained sand, while the water salinity measures 6 PSU (Gustafsson and Boström, 2014). Sediment was collected between the 1^st and 8^th of June 2021, sieved at a 1 mm mesh-size, homogenized, and air-dried for a minimum of 48 hours in the sun to eliminate any remaining macrofauna and larvae. This process, particularly effective in the subtidal and non-tidal conditions of the Baltic Sea, ensures the elimination of macrofauna and larvae, which are generally not adapted to prolonged desiccation. Seagrass shoots were obtained via SCUBA diving at water depths of approximately 2.3 to 3 m along the shoreline. The shoots were stored in cooling boxes filled with seawater on the boat and promptly transported to the field station, where they were placed in net-bags submerged in seawater for one day. Around 2000 M. balthica individuals were carefully hand-dug from the bare sand at a depth of about 0.2 to 0.5 m on the 8^th of June. On the same day, approximately 2000 Mytilus spp. were sampled from the pier structures at the field station (60°06’36”N 21°35’54”E). Any fouling organisms, such as barnacles and/or bryozoans, were gently removed from the mussel shells by hand or using a plastic spatula. Both bivalves were stored in the mesocosms filled with flow-through seawater until the experiment’s initiation.

Community assemblies were created at natural densities observed in the area, with 700 individuals of Mytilus spp. per m² (45 individuals per experimental unit), 716 individuals of M. balthica per m² (46 individuals per experimental unit), and 150 shoots of Z. marina per m² (10 individuals per experimental unit; Supplementary Figure S3).

On the 10^th of June, all 80 experimental units were filled with sediment to a height of 9 cm. For treatments involving Z. marina, ten shoots were randomly positioned within each experimental unit. Only terminal shoots of similar size (excluding side shoots) with a rhizome length of six internodes were used.

Following the planting of Z. marina shoots, a one-week acclimatization period was given before bivalves were introduced to the sediment surface of the respective experimental units. Prior to their addition, the bivalves were measured and divided into size classes that reflected a natural size distribution. For M. balthica, clams were measured and sorted into 7 size classes (0.6 cm to 2.0 cm in 0.2 cm steps) and the natural size distribution was mimicked by using the dominant 5 middle size classes (0.8 to 1.8 cm). For Mytilus, mussels were measured and sorted into 5 size classes (1.0 cm to 3.5 cm in 0.5 cm steps) and the natural size distribution was mimicked by using the first 3 size classes (1.0 to 2.5 cm). M. balthica and Mytilus spp. were added on the 18^th and 19^th of June, respectively. M. balthica individuals that had not burrowed into the sediment within two days were considered inactive or deceased and subsequently replaced by new individuals within the corresponding size class.

2.3 Sampling and laboratory procedures

The water temperature in the mesocosm system and at the pier near the water intake was continuously monitored using HOBO loggers (HOBO Pendant Data Logger Temperature 8K-64K, Onset Computer Corporation). These loggers recorded temperature at 15-minute intervals throughout the experiment. On the 12^th of June, 16 loggers were randomly distributed among the experimental units. Additionally, a multisensory WTW device (MULTI 3630 IDS, Wissenschaftlich-Technische Werkstätten) was used to measure oxygen (mg L^-1), pH (NBS scale), salinity (PSU), and temperature (°C) on a weekly basis in all experimental units, header buckets, and at the pier near the intake pump.

On the 17^th of August, midway through the experiment, we assessed the nutrient conditions in the mesocosms to examine the potential impacts of ecosystem engineers, such as bioturbation (by M. balthica), deposition (by both bivalves), and nutrient uptake (by Z. marina). Additionally, we assessed correlations between nutrient concentrations and the abundance of settling invertebrates, as well as the nutrient levels in both the water column and porewater.

For sampling the nutrients of the water column, seawater from the center of the experimental units was collected with a 60 mL syringe. Then about 30 mL of seawater was passed through a disposable syringe filter with a pore size of 0.45 μm (CHROMAFILR CA-45/25; MACHEREY-NAGEL GmbH & Co. KG; Düren, Germany) into a 50 mL falcon tube. The samples were kept in the dark and were immediately transferred to the freezer (-20°C).

The sediment pore water was sampled with Rhizon soil moisture samplers (10 cm; Rhizosphere Research Products; Wageningen, Netherlands) with an integrated 0.15 μm filter, attached to special vacuum bottles. The Rhizons were carefully inserted into the sediment at an approximate angle of 20°, allowing selective collection of pore water from the upper centimeters of the sediment. Approximately 35 mL of pore water was collected and then decanted into falcon tubes. Samples were kept in the dark and quickly brought to the freezer (-20°C).

Nitrite, Nitrate, Phosphate, and Silicate were analyzed by the autosampler SEAL QuAAtro Model 30 (Continuous Segmented Flow Analyzers) in the central lab of the GEOMAR, Kiel. Water column samples were analyzed without any dilution, while pore water samples were diluted 1:5 (2 mL sample volume and 8 mL artificial seawater) before analysis. Ammonium was analyzed after Holmes et al. (1999).

Since micro- and macroalgae growth can affect the settlement of invertebrates in hard and soft benthic bottoms (Österling and Pihl, 2001; Webster et al., 2015; Caronni et al., 2019), we quantified epiphytic growth in the experimental units. To determine epiphytic growth, PMMA (Plexiglass) slides of 9x9 cm were placed in each experimental unit on the 1^st of July, standing upright (70°) on the sediment. The slides were examined biweekly, by cleaning one side, while the other side of the panel was photographed in the lab under standardized light. Afterwards, this same side was cleaned, and photographed again to serve as a zero offset. Finally, the cleaned slide was placed back into the experimental unit. The pictures were analyzed using the software Image J (Schneider et al., 2012). Thereby the green color of the epiphytes in the pictures was converted with a grey coding (240 = white, 0 = black). The darkening was used as a proxy for the epiphyte cover on the slide. Hereby, a low darkening value corresponds to a low epiphyte cover on the slide, and a high darkening value corresponds to a high epiphyte cover on the slide.

Upon completion of the experiment on the 28^th of September, the sediment from each experimental unit was sieved using a 1 mm mesh-size. The collected residuals, along with the infaunal organisms, were transferred to labeled plastic vials, filled with 70% ethanol, and stored for subsequent analysis at -20°C. In the laboratory, each sample was examined under a stereomicroscope, and the individuals were separated and identified to the lowest possible taxonomic level. The number of individuals of each species in each sample was recorded. Dry weight biomass for each species was determined by weighing after drying for 48 hours at 60°C. The treatment (added) clams and mussels (>1cm) and the colonizing bivalves (maximal size 6 mm) collected at the end of the experiment, were distinguished by their very distinct sizes.

2.4 Data analysis

All statistical analyses were conducted using R 4.2.3 (R Core Team, 2021), and graphs were made with tools from the ggplot2 package (Wickham, 2016).

A set of community parameters, including total abundance, total biomass, species richness, and Shannon diversity index (both abundance- and biomass-based), were quantified for each experimental unit using the vegan package (Oksanen et al., 2020).

To explore the variation in water parameters such as pH, temperature, dissolved oxygen, and nutrient concentrations across different treatment levels, we employed non-parametric Permutational Multivariate Analysis of Variance (PERMANOVA; Anderson, 2001). This method was chosen as it does not require the assumption of multivariate normality and homogeneity of variances, conditions which our data did not satisfy. In our model, treatment was considered as a fixed effect factor. To account for random effects associated with experimental blocks, we incorporated these as a strata term within the model, thus restricting permutations to within-block levels. This approach ensures that the pseudo-F statistic accurately reflects the treatment effects while controlling for block variations. The PERMANOVA was executed on Euclidean distances and based on 999 permutations using the adonis2 function from the vegan package in R (Oksanen et al., 2020). We also used principal component analysis (PCA) as a complement analysis to assess the multivariate structure of the water parameters data.

Community parameters were modeled using Generalized Linear Mixed Models (GLMMs) from the lme4 package (Bates et al., 2015). All models treated Treatment as a fixed factor with eight levels, and potential block effects were accounted for by incorporating experimental blocks as a random factor. To select water parameters as covariates in these models, we first undertook a model reduction procedure to avoid overfitting. Initially, highly collinear variables were removed using a correlation matrix of Pearson correlations to identify pairs of variables with a Pearson correlation coefficient greater than |0.8|. The remaining variables were included as covariates in a full model. We then assessed the Variance Inflation Factors (VIF) of each covariate to further refine the model by removing collinear predictors. Subsequently, we performed a model selection routine based on the corrected Akaike Information Criterion (AICc) using tools from the MuMIn package (Bartoń, 2023). The final models considered only covariates that were present in the best model fits (ΔAICc < 3).

For modeling total abundance, total biomass, and Shannon indices, Gaussian error distributions were assumed. Species richness was modeled using a zero-truncated Poisson regression model with mixed effects, utilizing the glmmTMB package (Brooks et al., 2017). The importance of random effects was assessed using likelihood-ratio tests (LRTs). To maintain model parsimony and avoid overfitting, random effects that did not significantly improve model performance were dropped. To test for the significance of fixed effects, LRTs based on F-statistics were used for Gaussian models, and adjusted Chi-squared tests were used for mixed effects models. In instances where significant treatment effects were detected, pairwise comparisons between treatment levels were conducted using Tukey’s Honestly Significant Difference (HSD) tests from the emmeans package (Lenth et al., 2023). Model validation involved comparing the final models to a reduced model, excluding all predictors except for the intercept, through likelihood ratio tests.

We also compared the communities that settled in our mesocosms with five samples collected from Ängsö Bay during a pilot study in the summer of 2020. These samples were obtained through scuba diving and the use of sampling corers (15 cm diameter and 10 cm depth), and underwent the same protocols as those used for processing sediment communities from our experimental units. We normalized the abundances to match the area of our experimental units and employed Nonmetric Multidimensional Scaling (nMDS) based on Bray-Curtis distances of species abundances to contrast the field communities with those in the experimental units. This analysis was conducted utilizing tools from the vegan package (Oksanen et al., 2020).

3 Results

The PERMANOVA analysis of water parameters across treatments revealed a statistically significant effect of treatments on water conditions (pseudo-F = 2.344, df = 7, p < 0.05), accounting for approximately 20.6% of the total variation, with additional distinct block effects contributing another 24.3%. Further insights from PCA ordination (Supplementary Figure S5) identified epiphyte cover, pore water phosphate, and pore water ammonium as the most influential variables driving these differences. Notably, there was considerable variability within replicates of the same treatment level. However, higher levels of epiphyte cover and pore water ammonium were predominantly associated with the Mytilus, Zo + Ma, and Zo + My + Ma treatments. Conversely, elevated levels of pore water phosphate were generally linked to samples from treatments involving Z. marina (Zo + Ma, Zo + My, and Zo + My + Ma; Supplementary Figure S4, Supplementary Table S1).

After processing all 80 samples, we tabulated and identified a total of 34,131 settled macroinvertebrate individuals. Among these, we found 15 distinct species, in which four dominated more than 99% of the total abundance and biomass (Supplementary Tables S2, S3). The clam M. balthica stood out with 21,712 individuals, followed by the cockle Cerastoderma glaucum (Bruguière, 1789) (10,290 individuals) and the polychaete Hediste diversicolor (O.F. Müller, 1776) (1,369 individuals; Supplementary Tables S2, S3). Despite the numerical predominance of M. balthica, accounting for over 50% of individual counts in all treatments (Figure 2A, Supplementary Figure S5A), C. glaucum dominated in terms of biomass, representing more than 65% of the total dry weight across all treatments (Figure 2B, Supplementary Figure S5B). Remarkably, the four Zostera treatments exhibited reduced abundances of M. balthica and increased abundances of C. glaucum compared to the treatments without Zostera.

Figure 2 Relative abundance (A) and biomass (B) of the most common species that settled in each treatment. Stacked bars represent the sum of the total abundance or biomass of a species for each treatment. Treatments represent assembled communities encompassing sand-only, and all crossed combinations of the bivalves Mytilus spp. (My), Macoma balthica (Ma), as well as the seagrass Zostera marina (Zo).

Except for species richness, which mainly ranged between 4 and 5 species (Figure 3C) and no significant effect of predictors was detected, the results of the fitted GLMMs indicated that all community parameters exhibited significant variation across treatments (Table 1). Nonetheless, these models failed to identify any significant effect from the covariates integrated into the models on any of the community descriptors (Table 1). Overall, total abundance and biomass (Figures 3A, B) demonstrated similar trends across treatment levels, with the highest values generally observed in the Sand treatment and the lowest values in the Zo + My + Ma treatment. Total abundance ranged from 143 to 797 individuals per sample and was significantly lower in the Zo + My + Ma treatment (mean ± SD: 254 ± 58 ind. sample^-1) compared to most treatments but My + Ma and Zo + Ma (Figure 3A, Supplementary Table S4). The Zo + Ma treatment (320 ± 72 ind. sample^-1) was also significantly lower than the Sand (524 ± 192 ind. sample^-1) and Zostera (428 ± 166 ind. sample^-1) treatments. Total biomass ranged from 2.4 to 19.0 grams per sample, with the Zo + My + Ma treatment (5 ± 1.6 g sample^-1) being significantly lower than all treatments except My + Ma (Figure 3B, Supplementary Table S4). The total biomass in the My + Ma treatment (7.7 ± 2.5 g sample^-1) was also significantly lower than that in the Sand (13.6 ± 4.6 g sample^-1) and Macoma (12.6 ± 2.7 g sample^-1) treatments. Abundance-based Shannon diversity showed the highest values in the Zo + My + Ma treatment, being significantly higher than all treatments except Zo + Ma (Figure 3D, Supplementary Table S4). Conversely, the treatments lacking seagrass, such as Sand, Macoma, and Mytilus, displayed the lowest values (Figure 3D), being significantly lower than the Zostera, Zo + Ma, and Zo + My treatments (Supplementary Table S4). Notably, trends for biomass-based Shannon diversity differed from the abundance-based index, with the highest values observed in the Mytilus treatment (Figure 3E). Indeed, treatments with mussels exhibited higher values than other treatments, except for two comparisons: My + Ma vs Zo + Ma and Zo + Ma vs Zo + My (Supplementary Table S4). In general, treatments with a higher number of ecosystem engineers were associated with higher invertebrate diversity but negatively correlated with total abundance (Figure 4).

Figure 3 Variation in macrobenthic community parameters ((A) Total abundance, (B) Total dry weight, (C) Species richness, (D) Abundance-based Shannon diversity, (E) Biomass-based Shannon diversity) across treatment levels. Boxplots show the median (thick horizontal lines) separating the second (box below) and third (box above) quartiles with the data extension (vertical thin lines) and outliers (black dots). Treatments represent assembled communities encompassing sand-only, and all crossed combinations of the bivalves Mytilus spp. (My), Macoma balthica (Ma), as well as the seagrass Zostera marina (Zo).

Table 1 Summary of the likelihood-ratio tests checking for the effects of single exclusions of predictors used to model macroinvertebrate total abundance, total biomass, abundance-based Shannon, and biomass-based Shannon.

Figure 4 The relationship between total abundance, abundance-based Shannon diversity, and the number of ecosystem engineers. The fitted line represents a 2nd-order polynomial fit and its associated confidence interval (gray shade). Model fit coefficients and summary statistics are shown.

4 Discussion

Through our mesocosm experiment, we observed the distinct effects of three ecosystem engineers on different aspects of macroinvertebrate community structure. The outcomes were primarily influenced by the initial combination of ecosystem engineers present in the community. Particularly noteworthy were the effects of the presence of Mytilus spp. and the combined presence of all three ecosystem engineers, which had the most pronounced impact on community structure. Interestingly, the number of ecosystem engineers had a positive effect on the diversity of the settling communities. When all three ecosystem engineers were present together, a strong negative effect on total abundance and biomass was detected, suggesting possible synergistic interactions among ecosystem engineers on the assembly process of macrobenthic communities. These interactions likely prevent the dominance of C. glaucum and M. balthica, thereby favoring increased alpha diversity. The remarkable influence of ecosystem engineers, especially Z. marina and Mytilus spp., on the structure of benthic communities provides compelling evidence for their role in the (re)colonization of sedimentary benthic habitats.

4.1 Effects on macrobenthic diversity

In our ecological model system, despite its relatively low species diversity, we observed significant trends in the relationship between ecosystem engineers and macrozoobenthic colonizers, particularly when all ecosystem engineers were present. Here, we observed significantly lower total abundance and biomass of the colonizing benthic assemblages compared to other treatments. One plausible explanation for these observations is the intensification of competition for resources. The simultaneous presence of seagrasses, mussels, and clams likely leads to a scarcity of available substrate for settlement and likely intensifies the competition for food (Quinn, 1982). These altered environmental conditions reduce niche availability for new settlers, hindering their chances to establish and thrive (Quinn, 1982). Additionally, our data suggest a synergistic effect arising from the interactions among these ecosystem engineers, an effect that cannot be solely attributed to any individual species. The decreased abundance and biomass of the dominant bivalves C. glaucum and M. balthica, found under the presence of all three ecosystem engineers, primarily drove these effects. As a result, we generally observed a higher Shannon diversity, indicating a more species-diverse community. Similarly, other studies found elevated associated diversity in the presence of other co-occurring coastal ecosystem engineers such as seagrasses and pen clams (Zhang and Silliman, 2019), cordgrass and ribbed mussels (Angelini et al., 2015), or mangroves, fucoid algae, and oysters (Bishop et al., 2012). Therefore, our research indicates that in Baltic coastal ecosystems, the simultaneous presence of various ecosystem engineers can result in either synergistic interactions or elevated competition for space and nutrients. This dynamic leads to a notable reduction in the prevalence of specific bivalves, concurrently fostering an increase in overall species diversity.

The benthic system simulated in our treatments, which included all ecosystem engineers, likely represented a complex environment with increased microhabitat variety and a diverse array of chemical cues for pelagic larvae (Lecchini et al., 2010; Thiyagarajan, 2010; Cowan et al., 2016). The mutualistic relationships between seagrasses and infaunal bivalves exert a significant influence on sediment biogeochemistry. This influence is primarily manifested through seagrass root radial diffusion, which shapes sediment pore water chemistry, as well as the transformative effects of bivalve feeding activities, encompassing bioturbation and biodeposition (Meysick et al., 2020). Moreover, the presence of Mytilus spp. mussels in seagrass habitats leads to additional biogeochemical changes in the top layers of the sediment due to their biodeposition in the form of feces and pseudo-feces (Reusch et al., 1994; Mermillod-Blondin, 2011). These epifaunal bivalves also create gregarious patches that alter the microtopography of benthic surfaces (Borthagaray and Carranza, 2007; Norling and Kautsky, 2008). Consequently, this complex Mussel-Clam-Seagrass system, with enhanced habitat diversity and chemically rich cues, may support increased local diversity by limiting the dominance of single species (i.e., C. glaucum and M. balthica).

Mytilus spp. presence alone had a significant impact on biomass-based diversity, which can be attributed to facilitation and inhibition of species within mussel-dominated patches (Norling and Kautsky, 2008). One notable finding was the negative effect of Mytilus spp. presence on relative biomass and abundance of the cockle C. glaucum accompanied by a simultaneous positive effect on M. balthica compared to other simulated communities. This had a substantial influence on the computation of biomass-based Shannon diversity, which reached its highest values in the presence of Mytilus spp. alone. Although identifying the precise causal factors behind these patterns is challenging, they may be attributed to a combination of competition and facilitation processes among these bivalves (Miyamoto and Noda, 2004; Commito et al., 2005; Norling and Kautsky, 2008), leading to a more balanced distribution of biomasses among macrofaunal species in the presence of Mytilus spp.

The mussel Mytilus spp. is widely acknowledged as a dominant space competitor in the Baltic Sea and other regions, primarily due to its abundant recruitment and rapid growth in eutrophic environments (Norling and Kautsky, 2008). Possibly, Mytilus spp. has certain adverse effects on invertebrate settlement, primarily associated with mussel filtration activity and active substrate selection by larvae (Enderlein and Wahl, 2004). Consequently, it often dominates the community, leading to an overall decrease in local diversity (Reusch et al., 1994; Reusch and Chapman, 1997). However, Mytilus spp. also acts as a facilitator for numerous other invertebrate species, occasionally creating true hotspots of benthic diversity in areas that would otherwise be characterized by monotonous sandy bottoms (Norling and Kautsky, 2008). The association between the treatment with mussels alone and higher abundance and biomass of the clam M. balthica observed in our results possibly reflects the above-mentioned facilitation process. Although M. balthica may not directly increase structural complexity suitable for settling larvae, its role gains significance when interacting with other ecosystem engineers such as Mytilus spp. Clams such as M. balthica exhibit robust bioturbation effects, redistributing organic matter and releasing vital nutrients such as ammonium and phosphate from the sediment into the water column (Michaud et al., 2006). This interplay underscores the potential of combined engineer species to shape the ecosystem dynamics in soft benthic habitats in the Archipelago Sea.

4.2 Effects on the structure of settled communities

The structure of settled macrobenthic communities exhibited changes across the habitats simulated in our treatments in both, abundance-based and biomass-based analyses. The treatment consisting of Mytilus spp. alone showed a tendency to differ from all treatments encompassing the seagrass Z. marina. This indicates that the effects of mussels on settled macrobenthic communities differ somewhat from those observed in simulated seagrass meadows with or without adult infaunal bivalves (i.e., treatments with M. balthica). Furthermore, settled communities in treatments with seagrass tended to present similar species compositions, suggesting that the effects of Z. marina in structuring macrobenthic assemblages are strong enough to obscure the impact of co-occurring bivalve ecosystem engineers. Supporting this, Kindeberg et al. (2022) documented in Southern Sweden how Z. marina meadows can override the influences of abiotic sediment factors on macrobenthic assemblages. Similarly, González-Ortiz et al. (2014) observed in a field experiment in the Bay of Cadiz that the presence of seagrass significantly mitigates the impact of highly bioturbating fiddler crabs on sediments. These studies corroborate our findings that the effects of Z. marina are sufficiently dominant to mask the impacts of other co-occurring ecosystem engineers in shaping sediment fauna.

Among the tested treatments with ecosystem engineers, Mytilus spp. alone represents the most distinct habitat characterized by a system dominated by sessile bivalves forming stable patches, thereby covering and stabilizing the sediment surface (Reusch and Chapman, 1997). These conditions potentially favor bacteria with anaerobic metabolism (Chen et al., 2022). Indeed, our treatment with Mytilus spp. alone was characterized by higher levels of pore water Ammonium compared to Macoma and Sand treatments. The provision of such a chemically contrasting environment for microbes may impact the quality of biologically available organic matter for settling macrofauna, particularly deposit feeders (Miatta and Snelgrove, 2021). Indeed, Norling and Kautsky (2007) highlighted the functional role of live mussels in benthic communities, showing that community metabolism in mussel patches heavily depends on mussel biodeposition. This, in turn, can lead to differential settlement rates and contribute to contrasting community structures compared to benthic habitats dominated by other types of ecosystem engineers (Ysebaert et al., 2009). Additionally, Sokołowski et al. (2015) found that the quality of sedimentary food and macrophyte vegetation density in the Gulf of Gdańsk markedly affects macrofaunal communities. Therefore, the disparities we observed in the Mytilus spp. treatments relative to others may be the result of the unique chemical environment they create.

Although we cannot separate the physical effects of the presence of mussels and seagrass from their biogeochemical effects on soft sediments in structuring settling macrofaunal communities, our findings provide evidence that the variation in macrobenthic structure across benthic habitats, provided by mussel beds and seagrass meadows, is indeed correlated with chemical gradients in sediment nutrient concentrations. This underscores the functionally contrasting nature of the habitats examined in our mesocosm experiment and the implications that mussels and seagrasses have for the beta diversity across the Archipelago Sea, the Baltic Sea, and potentially in other systems of comparable habitat types and nutrient environments.

4.3 Ecological implications

Our findings provide compelling evidence and support previous findings that the effects of ecosystem engineers, such as seagrasses and bivalves, on recruiting benthic assemblages are intricately tied to the presence of other ecosystem engineers (González-Ortiz et al., 2014). This suggests that synergistic effects resulting from the interactions of functionally different ecosystem engineers play a significant role in shaping the soft benthos of the Archipelago Sea. Consequently, it can be inferred that the spatial variability and co-occurrence patterns of ecosystem engineers exert a critical influence on the alpha and beta diversities of benthic communities within the broader context of the Baltic Seascape.

Understanding the factors governing the spatial distribution of these ecosystem engineers is vital for unraveling the complex patterns observed in the region. The habitats created by the ecosystem engineers studied in this research significantly contribute to the heterogeneity of benthic communities (Norling and Kautsky, 2008). These habitats can be envisioned as functionally variable mosaics that not only interact with and influence one another but also have broader impacts on other components of the marine ecosystem (Meadows et al., 2012; van der Zee et al., 2012). As a result, the mosaic of benthic habitats shaped by these associations of ecosystem engineers assumes a pivotal role in the intricate food webs of the Archipelago Sea (Aguilar et al., 2017).

To effectively manage shallow marine benthic ecosystems, marine managers and spatial planners should understand and consider the complex interactions among different ecosystem engineers like seagrass, mussels, and clams. This knowledge is crucial for conservation and spatial planning, allowing for a better design of marine zones and resource allocation in a way that supports the natural dynamics of these ecosystems. Adopting a holistic approach, considering coexistence and interactions of key species, will enhance ecological integrity and lead to more effective preservation and sustainability of these unique environments.

4.4 The mesocosm settings

Our experimental settings encompass several factors that may have influenced our results, particularly regarding the scale of our mesocosm facility. These factors could partly account for the differences observed between the communities that formed within the mesocosm system and those found in natural seagrass meadows in, for example, Ängsö Bay (Supplementary Figure S6). When comparing our mesocosm-derived samples with those from other studies on benthic ecosystems in the Archipelago Sea, a notable contrast is the overwhelming predominance of only four species (H. diversicolor, C. glaucum, Hydrobia sp., and M. balthica) in our mesocosms, accompanied by a near-complete absence of Crustaceans. For example, the amphipod Monoporeia affinis (Lindström, 1855), identified as a dominant species in the benthic systems of the Archipelago Sea (Perus and Bonsdorff, 2004), failed to establish in our mesocosm setup. Likewise, other species such as the polychaetes Marenzelleria viridis (Verrill, 1873) and Pygospio elegans Claparède, 1863, common in the sediments of Ängsö Bay’s benthic systems (Boström et al., 2006b), were also absent in our mesocosms. These observations suggest that community dynamics within our experimental setups may differ significantly from those in natural environments in certain aspects.

As anticipated for a small semi-closed water habitat, the diurnal temperature variability within the mesocosms exceeded the ambient temperature recorded at the pier, although the general trend was similar to natural conditions (Figure S7). Such patterns are typical of closed or semi-closed environments, including tide pools and water ponds, which exhibit distinct ecological dynamics compared to more expansive bodies of water with greater thermal stability (Smith and Able, 2003; Smith et al., 2018; Vinagre et al., 2018).

The disparity in macrobenthic community structure when compared to natural benthic ecosystems can be attributed to the absence of vital post-settlement processes, which are fundamental in shaping benthic communities but are absent in our mesocosm setup. Notably, predation by fish and larger crustaceans, a critical factor influencing the composition of benthic communities in marine coastal ecosystems (Henseler et al., 2021), was entirely excluded from our mesocosm. Additionally, the recruitment process in our system might favor species that undergo indirect development and possess resilient larvae, which are more likely to survive being transported from the external water column into the mesocosm system.

Nonetheless, it is important to note that our experiment was conducted over a 101-day period during spring and summer months, which may have been insufficient for the benthic communities to reach a stable state within the experimental units (Boström and Bonsdorff, 2000). Additionally, it is crucial to recognize that our methodology to process sediment samples involved the use of 1 mm mesh-size sieves, effectively omitting meiofauna and smaller macrofauna elements from our analysis. This exclusion might be relevant since the meiofauna frequently accounts for a large portion of the biodiversity of invertebrate communities in soft sediments (Moreno et al., 2008). Yet, within these given boundaries, our experimental design effectively tested our specific hypotheses regarding the interactions of dominant ecosystem engineers in mediating macrobenthic colonization.

5 Conclusions

Our mesocosm experiment revealed that ecosystem engineers have a significant influence on the structure of settling macrobenthic assemblages. Seagrass, mussels, and clams possibly impact the benthic environment through direct effects on larval settlement, effects on new recruits, and indirect effects on sediment biogeochemistry. Interestingly, the combined effect of multiple ecosystem engineers differs from their individual effects, suggesting potential synergistic or antagonistic interactions of their effects on the benthic system. Understanding the spatial distribution and co-occurrence patterns of these engineers is essential for comprehending the complex patterns observed for the macrobenthic communities in the region. All in all, our mesocosm experiment provides valuable insights into the distinct effects of habitat-forming organisms and their contribution to the spatial structure of benthic assemblages in the Baltic Sea.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement

The manuscript presents research on animals that do not require ethical approval for their study.

Author contributions

JG: Conceptualization, Data curation, Formal Analysis, Methodology, Visualization, Writing – original draft, Writing – review & editing. SR: Conceptualization, Data curation, Investigation, Methodology, Project administration, Writing – original draft, Writing – review & editing. LK: Data curation, Investigation, Methodology, Writing – review & editing. LM: Conceptualization, Investigation, Methodology, Supervision, Validation, Writing – review & editing. CP: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was funded through the 2019-2020 BiodivERsA joint call for research proposals, under the BiodivClim ERA-Net COFUND programme, and with the Academy of Finland (grant nr.: 344743), and by Svenska Kulturfonden, as well as by The Sea at Åbo Akademi University.

Acknowledgments

This research gratefully acknowledges the contributions of C. Boström, N. Niemi, M. S. Machado Karekezi, A. Lindqvist, K. Gagnon, L.-R. Däschle, and L. Zimmermann in sampling and analyzing the infauna samples of the mesocosm experiment. Their expertise, dedication, and collaborative spirit were important to the success of this study. We also extend our appreciation to the entire research team and support staff involved in this project for their valuable contributions.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2024.1304442/full#supplementary-material

References