The fieldwork was conducted in the Wadden Sea, the world’s largest system of tidal flats, gullies, and barrier islands. Its total tidal areas cover nearly 10,000 square kilometers along the coast of the Netherlands, Germany, and Denmark (Reise, 2005) with around 1,000 km² of soft-sediment intertidal flats at mean low water and a total surface area of around 2,500 km² at mean high water in the Dutch part of the Wadden Sea (Elias et al., 2012). The Wadden Sea is classified as a UNESCO World Heritage site (Unesco, 2014) and is protected as a Natura 2000 area under the European Bird- and Habitat directives, because of its high productivity and supports large numbers of invertebrates, fish, and birds with many long-distance migrants (van de Kam et al., 2004; van der Veer et al., 2015). Around 1990, the Dutch Wadden Sea experienced a drastic loss of nearly all (~4000 ha) intertidal mussel beds, primarily attributed to overfishing, storms, and recruitment failures (Dankers et al., 2001; Christianen et al., 2017). In 1991, the first Pacific oyster that attaches to mussel beds was observed outside cultivation plots (Reise, 1998). Consequently, mussels located within the interspace of oyster reefs exhibited reduced growth rates and body conditions due to interspecific competition (Dankers et al., 2001; Folmer et al., 2017). Moreover, climate change has impacted the Dutch Wadden Sea, leading to a temperature increase of 2.3°C over the last ~120 years, a decline in wind speed averaging -2% per decade, an increase in heavy rainfall events subsequently influencing sediment transport, and occurrence of extreme droughts in the years 2018, 2019 and 2020 (Hoekstra and Philippart, 2022).

To test for the effect of L. conchilega tubes on mussel settlement, we performed a manipulative field experiment on an intertidal flat south of the island Schiermonnikoog (53°27’22.0”N, 6°9’30.2 “E) from August 2020 until March 2022 ( Figure 1 and Appendix S1A ). This tidal flat was characterized by bare sandy sediment and a tidal emergence time of approximately 30% per day. L. conchilega was present in medium-dense aggregations (~800 individual m^-2), and at ~200-meter distance patches of mussels and oysters were observed. Within the L. conchilega aggregations, we created a randomized block design with two treatments in 3x3 m plots: ‘tubes present’ and ‘tubes removal’ (n=6). The plots were orientated in blocks of two, with four meters in between the plots and eight meters in between each block. No manipulation of the treatment ‘tubes present’ was necessary since the tubes were naturally present. A sled with a blade was pulled over the sediment surface (top 1 cm) to remove all tubes and associated fauna for the ‘tubes removal’ treatment. To reduce the modification of the sediment properties, we installed a sieve (Ø 1 mm) behind the blade that redeposited the sediment after cutting the tubes. After tube removal, we suffocated all fauna by digging in agricultural plastic sheets made from low-density polyethylene. After six days, we removed the plastic and removed all dead organic tissue on the sediment surface. Since it was impossible to selectively remove L. conchilega only, we had to remove all benthic fauna to prevent immediate new tube build-up. Subsequently, we performed measurements in August, October, and December 2020, June and September 2021 and March 2022. During the sampling rounds, we counted L. conchilega tubes with visible fringes and M. edulis abundance four times in the center of each plot with 20x20 cm frames. M. edulis was classified as a larva when its length measured less than one centimeter. We counted tubes with visible fringes because this indicates occupation by the worm itself (van Hoey et al., 2006). The fringes of the tubes are in general destroyed shortly after the worm abandons the tube or dies (van Hoey et al., 2006). In addition, we analyzed sediment properties according to the procedures described in section 2.2.2.

To measure the facilitation of L. conchilega abundance in the adjacent habitat zones of mixed mussel/oyster reefs (M. edulis mixed with Magallana gigas), we conducted surveys along six transects with two distinct characteristics: covering a mussel/oyster reef and a bare control at the north-east of the island Texel (53°8’41.5”N, 4°54’18.7”E), in August 2019, March 2020 and August 2020 ( Figure 1 and Appendix S1A ). The control site was chosen to be as near to the reef as possible and selected based on a similar habitat type: elevation measured with Real Time Kinematics dGPS (Trimble R8, GNSS system; –36 - –40 cm N.A.P. at 50 m distance from the reef), distance to the gully (~ 50 m during low tide), sediment type (sandy without dead shells). Hence, the distance between the reef and control sites was ~1500 m at Texel. The distance between the transects within one treatment and sample points was 50 meters, measured from the center of the reef (0 m distance from the reef, Appendix S1B ). The sampling station –50 was positioned 50 meters from the hydrodynamically exposed edge of the reef. The sampling stations 50, 100, 150, and 200 meters were located in the wake of the reef. The sample points were positioned perpendicular to the incoming tide ( Appendix S1B ). At each sampling point, we counted L. conchilega tubes with visible fringes four times in 20x20 cm frames and estimated mussel/oyster reef cover in 50x50 cm frames (Van Hoey et al., 2008). In addition, we analyzed sediment properties according to procedures described earlier in section 2.2.2 Sediment properties sampling procedures and measured elevation with RTK dGPS. Furthermore, we compared the stable mixed mussel/oyster reef at Texel with a young mussel bed at the south of the island Griend (53°14’20.6”N, 5°15’31.6”E). This mussel bed disappeared within one year after the start of the monitoring. The survey set-up replicated that of Texel, with the exception that at Griend the distance between the reef and nearby control was ~300 m. This adjustment aimed to maintain constant environmental parameters including elevation and sediment type.

To monitor the establishment of a mussel bed on L. conchilega aggregations on a temporal scale, we surveyed reef development south-east of the island Griend (53°14’34.6”N, 5°18’6.5”E; Figure 1 and Appendix S1C ). In August 2019, we measured the contours of dense L. conchilega aggregations with handheld GPS (Garmin GPSMap 66 Series). In April 2021, this reef was displaced by a young mussel bed. We measured the contours of this reef with drone imagery (DJI Mavic Pro 2) and handheld GPS and determined contours in qGIS (3.6.3-Noosa). In addition, we surveyed transects (n=9) that covered the young mussel bed and adjacent bare tidal flats in April and August 2021 ( Appendix S1C ). On these transects, we 1) counted L. conchilega tubes with intact fringes in 50x50 cm frames, and 2) measured elevation with RTK dGPS (Trimble R8, GNSS system) through sampling points (every 5 m).

We performed all statistical analyses in R (R Core Team, 2022). We validated model assumptions by 1) residuals versus fitted value plots to verify homogeneity of variances, 2) QQ-plots of the residuals to test for normality and 3) residuals versus each explanatory variable to check for independence. In addition, Shapiro-Wilks’s test (p > 0.05) and Bartlett’s test (p > 0.05) were used to test for normality and homogeneity of variance of the residuals, respectively.

For the removal experiment on Schiermonnikoog, M. edulis and L. conchilega tube counts were analyzed using generalized linear mixed models (GLMM) assuming a negative binomial, zero-inflated error distribution (absence was affected by the removal treatment; r-package: ‘glmmTMB’ (Brooks et al., 2017). The zero-inflation model (ZI) estimates the probability of an extra zero such that a positive estimate value indicates a higher chance of absence (e.g. ZI: tubes removal >0 means more absence). This is the opposite of the negative binomial model (NB) where a positive contrast indicates a higher abundance (e.g., NB: tubes removal >0 means higher abundances in this treatment) (Brooks et al., 2017). The sampling date and plot nested in blocks were included as random effects. For the transects on Texel, we summarized the counts of four pseudo-replicate L. conchilega tube counts per sampling station and analyzed these numbers by generalized linear models (GLM) with negative binomial distribution (glm.nb function from the R package MASS (Venables and Ripley, 2002). Tukey’s post-hoc comparisons were used to test for significant differences (p < 0.05) between treatments within distance from the reef (R package ‘emmeans’; (Lenth, 2019). Sampling dates (August 2019, March 2020, and August 2020) were used as replicates since there was no variation in means among the sampling time points. To test for differences between the treatments within distance from the reef, sediment D50 and silt percentage were analyzed with GLM with a negative binomial distribution, organic matter content was log-transformed and analyzed with linear models, and Tukey’s post-hoc comparisons were applied. To monitor biogenic reef development, we spatially visualized L. conchilega tube counts by ordinary kriging r-package: ‘automap’; Hiemstra, 2022). The function ‘autoKrige’ fits a variogram model to the given data set and returns the results of the interpolation: prediction, variance, and standard deviation. Consequently, the interpolated prediction was used to visualize L. conchilega tube counts.

To investigate the interlinkage between M. edulis and L. conchilega, we examined the occurrence frequency of both species at individual sample locations within the SIBES dataset. We found that both foundation species were positively associated with each other because the observed value of sample locations with the presence of both M. edulis and L. conchilega was 2.45 times higher than the expected value (X² (1, N=42409) = 831.45, p < 0.001, Figure 2 and Appendix S2 ). In addition, MaxEnt species distribution modeling revealed that 125 km² and 81 km² of the total intertidal area of the Dutch Wadden Sea was highly suitable (occurrence probability of > 0.6) for M. edulis and L. conchilega, respectively ( Figure 3 ). Only 37 km² of the ecological niches between both foundation species overlapped ( Figure 3 ), which means that they share 30% and 46% of their niche with the other foundation species for M. edulis and L. conchilega respectively. The highly suitable areas for M. edulis were positioned closer to the mainland (training data AUC = 0.705), while these areas for L. conchilega are found closer to the barrier islands and tidal inlets (training data AUC = 0.734) ( Figure 3 ).

To investigate the role of L. conchilega tubes as settlement substrate for M. edulis spat within nested assemblages, we experimentally removed L. conchilega tubes and measured the effect on M. edulis. High L. conchilega tube densities of 1542±160 m^-2 (mean±SE) facilitated the settlement of M. edulis larvae (length: ~ 1 cm) (Figure 4A). M. edulis counts were 400 times higher in the ‘tubes present’ treatment compared to the ‘tubes removal’ treatment with 177±36 individual tube counts (mean±SE) ( Figure 4B and Table 1 ). This facilitation within nested assemblage was visible until the end of the experiment in April 2022, almost one year after the first M. edulis settlement. Sediment properties did not differ between treatments, indicating that there was no significant attenuation of hydrodynamics by L. conchilega tubes ( Appendix S3 ).

To explore the facilitation of L. conchilega by mixed mussel/oyster reefs (M. edulis mixed with M. gigas), we surveyed L. conchilega tube densities and reef cover. No tubes 0±0 (mean±SE) were counted at sampling stations -50 m and 0 m mussel bed. The highest tube densities of 1606±122 m^-2 (mean±SE) were found at 50 meters in the wake of the reef, and the lowest tube densities of 128±38 m^-2 (mean±SE) were found at 200 m in the wake of the reef. The tube counts were on average 24 and 22 times higher at 50 (β = 3.17, SE ± 0.66, p < 0.001) and 150 m distances (β = 3.08, SE ± 0.67, p < 0.001), respectively, in the wake of the reef compared to the nearby control ( Figure 5 and Table 1 ). The mussel/oyster reef area extent covered the sampling stations -50 and 0 m with an average reef cover of 45% at sampling station -50 and 75% at sampling station 0 ( Figure 5A ). At 50 meters in the wake of the reef, sediment silt percentage was higher (+17% silt), median grain size was lower (D50: –32 μm) and organic matter content was higher (+1% LOI) compared to the nearby control ( Figure 5B and Appendix S4 ). The sediment surface around the reef center was elevated by ~ 0.35 meters (0 meters; Figure 5C ). To compare these counts with a young mussel bed located south of the island Griend that disappeared within one year after the start of the monitoring, we observed the highest L. conchilega tube densities 3956±631 m^-2 at 50 meters in the wake of the reef compared to the control with 171±106 tubes m^-2. One year later, however, this mussel bed had disappeared, and L. conchilega tube densities were reduced to almost zero counts ( Appendices S1 , S5 ).

To study the interlinkage between M. edulis and L. conchilega, we monitored biogenic reef development ( Figure 6 ). In 2019, the monitored reef consisted of almost exclusively L. conchilega tubes with 2201±121 individual tubes m^-2 (mean±SE, n=20) and covered a surface area of 383,507 m² ( Figure 6 ). In 2021, the reef was instead dominantly shaped by M. edulis, and L. conchilega tubes were mostly present in the adjacent habitat zone (50 – 150 m). The total surface of the biogenic reef dominated by M. edulis was 33,047 m² and the highest L. conchilega tube densities (660 individual tubes m^-2) were found within 50 m in the wake of the reef in 2021. Hence, biogenic reef development of M. edulis and L. conchilega indicated spatiotemporal dynamic habitat heterogeneity that consequently is important to promote biodiversity (van der Ouderaa et al., 2021).

Our results showed a positive association between the occurrence of M. edulis and L. conchilega. In addition, M. edulis and L. conchilega share 30% and 46% of their ecological niche, respectively. These combined results suggest that interspecific facilitation may contribute to the enlargement of the ecological niche into the realized (and observed) ecological niche. Enlargement of the niche by facilitation has been described previously for mussel beds that stabilize the substratum and increase salt marsh production via the deposition of nutrients into the sediment (Bertness and Leonard, 1997; Stachowicz, 2001). Analyses of the large SIBES dataset (sampling stations: n=42409) establish a robust basis for the assumption that M. edulis and L. conchilega are positively correlated, potentially enhancing each other’s fundamental niche. For example, interspecific facilitation can ameliorate the environmental stressors that otherwise would limit species occurrences, such as providing settlement substrate or reducing hydrodynamic forces (Callaway, 2003; Walles et al., 2015).

Previous studies found that M. edulis settlement prefers structured habitat surfaces over bare substrata (Dean, 1981; Bourget et al., 1994) and that the structure of imitation tubes can provide settlement substrate for M. edulis larvae (Callaway, 2003). However, imitation tubes differed from natural tubes of L. conchilega because they were more rigid, more solid, not hollow, more persistent over time, and had no tentacle fringes (Callaway, 2003). This raised the question of whether natural L. conchilega tubes can also facilitate M. edulis larvae settlement. In this study, we showed that L. conchilega can indeed provide a suitable settlement substrate for M. edulis spat in natural conditions. The mechanisms behind this facilitation may be explained by two factors: the structure of the tubes reduces the hydrodynamic forces (Borsje et al., 2014) and/or the tubes provide a complex attachment substrate (Callaway, 2003). M. edulis larvae may use different factors to find a suitable substratum for settlement, e.g., physical [flow conditions, light, or gravity; (Bayne, 1964; Pernet et al., 2003)] or chemical [metabolites of algae or invertebrates; (Dobretsov, 1999)] cues. Our findings suggest that L. conchilega tubes provided settlement substrate for M. edulis spat, as we did not find indications of altered hydrodynamics based on sediment properties in the tubes removal plots compared to the tubes present plots (Le Hir et al., 2000). However, an additional effect of chemical metabolites produced by L. conchilega attracting M. edulis cannot be excluded since we have not been able to measure such cues.

By attenuating hydrodynamic forces from waves and currents, bivalve reefs typically have a much larger-scale impact on their surrounding physical landscape than the area they actually cover (Walles et al., 2015). This may facilitate other species such as seagrasses, salt marsh plants, and cockles in adjacent assemblages i.e., long-distance, cross-habitat facilitation (Wall et al., 2008; Donadi et al., 2013; van de Koppel et al., 2015; Walles et al., 2015). These reefs create habitat for other species by e.g., trapping water and increasing tidal emergence time (Nieuwhof et al., 2018; van der Ouderaa et al., 2021) or wave attenuation (Cheong et al., 2013). Walles et al. (2015) found that the spatial extent of habitat modification by bivalve reefs is of the same order of magnitude as the reef area (length, width, height). In our study area, the sediment surface of the bivalve reef was around 100 meters wide and on average 0.35 meters higher than the surrounding tidal flat. A similar elevation of the bare tidal flat to the height of the reef was reached at around 110 – 140 meters in the wake of the reef. Hence, the modified area leeward of the reef was expected to be around 100 – 140 meters. In addition, we found that the bivalve reefs increased silt percentage and reduced median grain size up to 50 meters in their wake. Both measurements are indicators of reduced hydrodynamic forces (Le Hir et al., 2000). However, mixed mussel/oyster reefs facilitated L. conchilega occurrence at 50 and 150 meters in their wake, but not at 100 meters. This might be because of the high patchy distribution of L. conchilega (Callaway et al., 2010). In general, our results imply that bivalve reefs create suitable habitats for L. conchilega in adjacent assemblage by mitigating hydrodynamic forces.

We observed that L. conchilega tubes facilitated M. edulis settlement and consequently reef establishment south-east of the island Griend. In addition, we observed that the M. edulis reef facilitated L. conchilega densities in the wake of the reef, which is hydrodynamically less exposed. This suggests that interspecific facilitation may promote the co-occurrence of both species by creating temporal shifts between both species. Our findings find support in the observations made by (Hertweck, 1995), who found M. edulis spatfall attached to L. conchilega tubes in dense aggregations on the tidal flats south of the island Spiekeroog, German Wadden Sea. In addition, L. conchilega is frequently found surrounding mussel beds (Hertweck, 1995; Bungenstock et al., 2021). Although the interplay between both species appears evident, these cause-effect relationships must be interpreted with caution since they rely on observational data (Hertweck, 1995; Bungenstock et al., 2021). Moreover, additional research approaches, such as experiments and species distribution modeling, further support the interdependency of the two species. This extensive effort was essential for testing the interdependency across various local contexts, long-term temporal scales, and large spatial scales, ensuring objectivity in statistical analyses. Yet, the question remains which environmental factors promote the dominance of a certain foundation species? For example, it is known that cold winters limit the growth of L. conchilega populations (Strasser and Pieloth, 2001). On the other hand, L. conchilega tubes are more resistant to storms or winter ice because they are anchored deeper (~30 cm) into the sediment than M. edulis beds (Alves et al., 2017). Furthermore, M. edulis beds are attached to the sediment surface and winter ice can dislodge the mussels from the substratum (Bertness and Grosholz, 1985; Bungenstock et al., 2021). Therefore, this study emphasizes the importance of supporting optimal growth strategies for biogenic reef development by biological interactions.

Long-term reef development of M. edulis and L. conchilega on tidal flats is a dynamic process shaped by complex interactions between environmental factors and biological interactions. M. edulis spat aggregate densely on substrates, such as rocks, and shells (Wehrmann, 2003). Over time, the accumulation of mussel beds alters sedimentation patterns, stabilizes substrates, and enhances the survival of the reef (Widdows et al., 1998). In this study, we described that L. conchilega tubes function as settlement substrates for M. edulis larvae. Eventually, the M. edulis spat may develop into the physical habitat structure of a M. edulis bed, as evidenced by Callaway (2003) findings from a six-year-long experiment. The settlement of L. conchilega larvae relies on the hydrodynamic conditions and the substrate to settle on (D’Hurlaborde et al., 2022). Whether L. conchilega aggregations will develop into a biogenic reef, depends on their resistance to disturbance, recovery potential, and their large-scale persistence (Rabaut et al., 2009; Callaway et al., 2010). We described that M. edulis beds create favorable conditions for L. conchilega reef development in their wake. Although we were unable to assess this long-distance facilitation over an extended duration, similar long-distance facilitation of the edible cockle Cerastoderma edule by M. edulis beds was reported in 2009 and 2011 (Donadi et al., 2013). Hence, the influence of interspecific facilitation seems to contribute to the long-term development of biogenic reefs.

The existence of biogenic reefs such as M. edulis and L. conchilega relies on the establishment and survival of the reefs. This existence can be enhanced by the optimization of the environmental variables as settlement substrate and/or reduced hydrodynamics. For instance, the protection of areas with relatively stable L. conchilega aggregations can improve reef settlement of M. edulis. Re-establishment of mussel beds after storm events may be encouraged when L. conchilega aggregations are present. For active restoration, it might be necessary to temporarily mimic the physical structure of biogenic reefs with artificial structures to improve the local environmental conditions for reef establishment (Temmink et al., 2022; Nauta et al., 2023). However, sufficient structure and size of the biogenic reef is needed to achieve the restoration of landscape-scale connectivity (Gillis et al., 2017). Biogenic reef structure and size of L. conchilega aggregations can be classified by using a ‘reefiness’ score with highest ‘reefiness’ scores for densities of >1,500 individual m^-2, a total area extent of >100,000 m², and a relative height of >9 cm (Hendrick and Foster-Smith, 2006; Rabaut et al., 2009). Based on this ‘reefiness’ classification, the L. conchilega aggregation in this study can be classified as high-scoring biogenic reefs. Restoring landscape-scale connectivity means that the one biogenic reef-building species needs to be restored to a sufficient physical scale to generate positive effects on the other (foundation) species. The methodology applied in this study, incorporating field surveys, experiments, and species distribution models, could serve as a recipe for further understanding and co-restoration of interspecific dependency in the Dutch Wadden Sea. Moreover, biogenic reefs including M. edulis and L. conchilega provide numerous ecosystem services such as enhanced biodiversity and fisheries, sediment stabilization, and reduced hydrodynamics (Rabaut et al., 2007, 2010; De Smet et al., 2015; Alves et al., 2017; zu Ermgassen et al., 2020). Therefore, we suggest that possibilities for co-restoring foundation species with other foundation species need to be further explored.