Ten sites along the coast of Europe were selected for quantifying the diversity of benthic macrofaunal assemblages: four in the United Kingdom, four in France and two in Portugal (Figure 1). Sampling was carried out in the summer (spring tides of June and August 2017) following a standard protocol at each site. The sampling strategy aimed to maximize the number of species collected by sampling in a variety reef phases (prograding and retrograding, sensu Curd et al., 2019) within each site, as these are known to harbor different assemblages (Dubois et al., 2002; Jones et al., 2018). Reefs were sampled using eight PVC cores of 5 cm in diameter to a maximum depth of 15 cm. Since honeycomb worms occur within the first 15–20 cm of the reef (Gruet, 1986), only the living portion of the reef was sampled. At UK4, only veneer bio-constructions were available for sampling, and only five cores were collected because veneers were too scarce for further sampling. The contents of each core were preserved in 70% ethanol.

In the laboratory, cores were first weighed (wet weight, after removal of alcohol), then sieved on a 1 mm circular mesh. For a given core volume, the weight of the sediment provides means for comparing porosity or the void fraction of a sample. Macrofauna was then extracted from the sediments and enumerated. Individuals were identified to the lowest taxonomic level, most often to the species level. All species names were used according the World Register of Marine Species¹ and references used for taxonomic identification can be found in Supplementary Appendix 1. To ensure consistent taxonomic resolution across samples, the number of operators was limited (n = 4) and each uncertain identification was cross-verified by an expert in benthic taxonomy. However, due to the uncertainty regarding the morphological distinction at the species level between Mytilus edulis and Mytilus galloprovincialis, particularly at the juvenile stage (Jansen et al., 2007) and because hybridization occurs between the two (Daguin et al., 2001), all specimens sampled in the hybridization area (UK1, UK4, FR1, FR2, FR3, and FR4) were considered as Mytilus spp. (Wenne et al., 2020). All specimens are stored in the Laboratory of Coastal Benthic Ecology’s collections at Ifremer (Plouzané, France).

Given that the reef-building species S. alveolata affects resource availability for the associated community (Dubois et al., 2002), all analyses considered the density of this species as an explanatory variable. Hence, S. alveolata abundance was removed from the species (response) matrix. Although some studies exclude rare species (i.e., represented by a single individual in one or two samples) for the calculation of similarity (see Clarke and Warwick, 2001), they may represent a non-negligible portion of a functional group, and are likely to have an impact on ecosystem functioning (Leitão et al., 2016). Rare species were therefore retained in the analysis in order to avoid reducing the functional richness of the communities (Mouillot et al., 2013; Jain et al., 2014).

To obtain information on the environmental conditions at each site, water and air temperature data were recorded using iButton^® temperature loggers (accuracy ±0.5°C, hourly measurements) (Lima and Wethey, 2009), deployed between August and October 2016 for a period of approximately 1 year prior to the sampling campaign. In order to mimic temperatures similar to those experienced within the reefs, the loggers were coated with sand from the reefs and fixed onto rocky substrate at a constant shore level (corresponding approximately to the mid-tide level where the majority of reefs develop).

Multivariate analyses were used to test for differences in macrofaunal assemblages across four biogeographic provinces (as delimited by Dinter, 2001). Non-metric multidimensional scaling (nMDS) was used to plot sample stations on a two-dimensional ordination plane based on taxa composition dissimilarities and labeled with the corresponding collection site. nMDS was also run using species abundances averaged across all stations sampled in a given site. In addition, a hierarchical cluster analysis (HCA) was run on all samples using the Ward method (Ward, 1963). All analyses were carried out on the basis of Bray–Curtis similarity matrices. Abundance data were transformed by the log (x + 1) function to reduce the weight of the most abundant species.

In order to examine regional variation in α diversity, each site was coded as belonging to one of four biogeographic regions: Boreal, Boreal-Lusitanian, Lusitanian-Boreal, and Lusitanian (Figure 1; Dinter, 2001). Differences in community composition within and among regions were tested using PERMANOVA with a two-factor design (4999 residual permutations under a reduced model), with region as the fixed factor and site as the nested random factor. The weight of sediment was included as a co-variable in all analyses. Paired tests between regions were performed where the main effect was significant (P < 0.05). Prior to the PERMANOVA, differences in within-site multivariate dispersion were examined using the PERMDISP routine. When significant differences in assembly structure between regions were detected, a SIMPER analysis was performed to determine and rank the taxa responsible for the dissimilarities among sites and biogeographic regions. Variation of univariate assemblage metrics (i.e., abundance, species richness, the exponential of Shannon entropy (N1), and the inverse Simpson concentration (N2); Hill, 1973; Jost, 2006) were examined with permutational ANOVA, using the Euclidean distance in the PERMANOVA procedure (Anderson, 2017).

To qualify the link between the environment and macrobenthic community structure at each site, a distance-based redundancy analysis (dbRDA; Legendre and Andersson, 1999) was carried out. In addition to the variables recorded in the field (density of the engineering species S. alveolata, weight of sediment, maximum water and air temperatures), mean monthly values spanning 2000–2014 were obtained for 30 variables from BioORACLE² (Tyberghein et al., 2012; Assis et al., 2018). In order to eliminate multicollinearity among these environmental variables, Spearman rank correlations were calculated for all pairs of variables. Pairs with a Spearman correlation coefficient >0.7 were considered highly correlated. Only the following uncorrelated variables were kept in the analysis, in addition to the variables recorded in the field: mean surface water temperature (°C), mean chlorophyll a concentration (mg.m^–3) and maximum current speed (m.s^–1) (Dormann et al., 2013). Prior to analysis, all values for the environmental variables were standardized, then a DISTLM routine was used to obtain the most parsimonious model using a stepwise selection procedure and adjusted R² selection criterion (McArdle and Anderson, 2001).

To characterize the functional diversity at each site, a biological trait analysis (BTA) was conducted (Statzner et al., 1994). Eight biological traits (divided into 32 modalities) were selected (Table 1), providing information linked to the ecological functions performed by the associated macrofauna. The selected traits provide information on: (i) resource use and availability (by the trophic group of species, e.g., Thrush et al., 2006); (ii) secondary production and the amount of energy and organic matter (OM) produced based on the life cycle of the organisms (including longevity, maximum size, and mode of reproduction, e.g., Cusson and Bourget, 2005; Thrush et al., 2006); and (iii) the behavior of the species in general (i.e., how these species occupy the environment and contribute to biogeochemical fluxes through habitat, movement, and bioturbation activity at different bathymetric levels, e.g., Solan et al., 2004; Thrush et al., 2006; Queirós et al., 2013). Species were scored for each trait modality based on their affinity using a fuzzy coding approach (Chevenet et al., 1994), where multiple modalities can be attributed to a species if appropriate, and allowed for the incorporation of intraspecific variability in trait expression. The information concerning polychaetes was derived primarily from Fauchald and Jumars (1979) and Jumars et al. (2015). Information on other taxonomic groups was obtained either from databases³ of biological traits, publications (Caine, 1977; Leblanc et al., 2011; Rumbold et al., 2012; Jones et al., 2018) and publications listed in Supplementary Appendix 1.

Ordination of the functional trait data was done using a Fuzzy coded multiple Correspondence Analysis (FCA) (Chevenet et al., 1994). Then, a hierarchical clustering analysis based on the Ward algorithm (Ward, 1963) was carried out using Euclidean distances (Usseglio-Polatera et al., 2000) to define homogeneous functional groups comprising species with similar biological trait associations. The frequencies of the modalities of each trait were calculated in order to visualize the biological profiles of identified functional groups.

Regional differences in diversity (β diversity) were estimated from presence-absence data using Sørensen’s (1948) dissimilarity. For each pair of cores, taxonomic β diversity and its two components, turnover and nestedness, were computed using the Baselga partitioning scheme (Baselga, 2017; Schmera et al., 2020). Functional β diversity was computed based on the multidimensional functional space from the Fuzzy Correspondence Analysis, where axes were synthetic components summarizing functional traits (Villéger et al., 2010). The first four axes were used for calculating Sørensen dissimilarity according to Villéger et al.’s (2013) equation for all pairwise comparisons between samples (1) belonging to the same region (within bioregion), or (2) belonging to different regions (among bioregion). Correlations between taxonomic and functional β diversity as well as between their respective components were tested using Mantel permutational tests (Villéger et al., 2013).

All index calculations and statistical analyses were performed using the “betapart” (Baselga, 2012), “ade4” (Dray and Dufour, 2007) and “vegan” (Oksanen et al., 2020) packages in R 3.4.4 (R Development Core Team, 2008) as well as using PRIMER 7 (Clarke and Gorley, 2015) for PERMANOVA, PERMDISP, SIMPER and DISTLM analyses.

A total of 129 taxa were observed in association with honeycomb worm reefs across the 10 sampled sites (77 stations, Supplementary Table 1). Taxon richness varied from two taxa (core from UK3, FR2, and PO2) to 34 taxa per station (core from FR2). In all sites except FR3, S. alveolata was the dominant species, with densities ranging from 6,450 ind.m^–2 at UK2 to 80,000 ind.m^–2 at FR2 (Supplementary Figure 1A). In all sites except FR3, S. alveolata was the dominant species, with densities ranging from 6,450 ind.m^–2 at UK2 to 80,000 ind.m^–2 at FR2 (Supplementary Figure 1A). The highest densities were observed in UK4 (130,000 ind.m^–2) but this most likely corresponds to a recent recruitment event, as individuals were on average much smaller (diameter of the opercular crown less than 2 mm; Gruet, 1986). The ratio of individuals of S. alveolata to individuals of associated macrofauna showed a clear dominance of the engineering species at UK4 and UK3 (94 and 87%), while at other sites, this ratio varied from 43 to 68% (Supplementary Figure 1B). It was also at these two sites that the number of individuals of the associated macrofauna were the lowest, reaching up to 5,820 ind.m^–2 at UK3 (with an average of 3,507 ± 2,313 ind.m^–2) and 7,500 ind.m^–2 283 at UK4 (with an average of 6,194 ± 1,304 ind.m^–2; Supplementary Figure 1A). The fauna associated with honeycomb worm reefs were primarily of annelids and arthropods and to a lesser extent, mollusks and nematodes (Supplementary Figure 2). However, community composition varied significantly among sites. Annelids were dominant in UK2, UK3, FR4, and PO2, where they represented between 36 and 56% of individuals. At UK1 and PO4, arthropods dominated the community, representing 47 and 51% of individuals, respectively. For the other sites, mollusks were dominant (53% of individuals in UK4 and 39% in FR2) as well as nematodes (36% of individuals in FR1 and 39% in FR3). Sites FR1 and FR3 had a higher abundance of nematodes and mollusks compared to other sites.

Taxon richness and abundance did not show any significant variation between sites within the same region but showed significant differences between regions (Table 2). As for N1 and N2 indices, honeycomb worm reef communities were characterized by low values (Figure 1). For both metrics, there was no significant effect of site or region. Sediment weight had a significant effect, but only on the N2 index (Table 2). The PERMANOVA also detected a significant effect of sediment weight on taxon richness but not on abundance (Table 2).

Hierarchical cluster analysis based on species composition and abundances per station defined four groups of assemblages (Figure 2B). Group I included stations sampled in UK1 and UK4, plus two stations from UK3 and one station from FR1. Group II included all stations in UK2 and the remaining stations from UK3. Group III included all but one station in France (FR1, FR2, FR3, and FR4) and Group IV included stations in Portugal (PO1 and PO2). The nMDS indicated some degree of partitioning among sites, with two of the British sites (UK1 and UK4) being distinct from two of the southern sites (nMDS1 and nMDS2; Figure 2A). However, much overlap was observed in community composition among sites found in the center of the distribution, particularly FR1, FR3, and UK3, which exhibited considerable variability among stations within each site (F = 7.98, P < 0.001; PERMDISP; Figure 2A). PERMANOVA detected significant variability between regions and between sites within regions, as well as a significant effect of the sediment covariate (Table 2). Regional differences were driven by differences in intra-regional variability (F = 18.08, P < 0.001; PERMDISP; Figure 2A) but also by changes in mean species composition (Figures 2A,C). Paired tests between regions showed that assemblages from the Lusitanian Province (PO1 and PO2) were distinct from those from the Boreal province (UK1, UK2, and UK3) and the Lusitanian-Boreal province (FR1, FR2, FR3, and FR4) (Supplementary Table 2). SIMPER analysis indicated that the differences observed between regions were mainly due to a higher abundance of the polychaete Syllis armillaris (Müller, 1776) and the mussel M. galloprovincialis (Lamarck, 1819) in the Lusitanian Province compared to the Lusitanian-Boreal and Boreal Provinces (Supplementary Table 2).

Distance-based redundancy analysis analyses indicated some degree of partitioning between regions. The first two axes represent 12.4 and 10.5% of the explained total variance, respectively (Figure 3). The species assemblages were structured along two gradients. The first was driven by the mean chlorophyll a and mean water temperature variables, which were negatively correlated, highlighting the differences between the northern and the southern assemblages. The second was driven by the amount of sediment in the cores and the maximal temperature of the air, which separated the middle range sites from the southern and northern sites, with higher values in the middle range sites (Figure 3). Note that the amount of sediment in the cores was also negatively correlated with maximal current velocity, which was higher in the northern sites. The DISTLM routine was used to determine links between environmental predictor variables and variability in assemblage structure (Table 3). Marginal tests showed that mean water temperature and the amount of sediment in the cores were, individually, the most important predictor variables. Surprisingly, the density of S. alveolata did not appear to be a structuring factor for these intertidal communities. The stepwise selection procedure indicated that the most parsimonious model included all environmental variables, explaining 0.38% (adjusted R² = 0.31) of the total observed variation in this assemblage structure.

Changes in community structure were also analyzed in terms of functional diversity. Cluster analysis carried out on the FCA axes revealed five main groups of taxa with distinct trait combinations (Figure 4A). The most clearly delineated group (Group 1) was composed mostly of intertidal, suspension-feeding, small, long-lived organisms that live mostly fixed or in tubes and release their gametes into the water column (Figure 4B). This group was represented by 14 taxa, the majority of which were bivalve mollusks, sabellidae polychaetes, and barnacles. The four other groups were largely composed of infaunal taxa. Group 2 was composed of very small, free-living and tube-dwelling, short-lived, sometimes annual, organisms. Most species in this group lay or incubate eggs and have no larval phase. Their contribution to sediment reworking was mainly at the sediment surface. This group was composed of 47 taxa that included amphipods and isopods but also pycnogonids and nematodes. Group 3 comprised large, average-lived organisms, that freely release their gametes into the water column, with feeding modes mostly associated to scavenging and sub-surface deposit-feeding, living free or in burrows, and participating to sediment reworking, either as biodiffusors or through vertical sediment transport. This group included 23 taxa, belonging to the polychaete annelids Eunicidae, Lumbrinereidae, and Oenonidae (formerly part of Eunicidae) and Terebellidae. Group 4 and Group 5 comprised species with heterogeneous and intermediate trait characteristics relative to the more functionally homogeneous Groups 2 and 3, the former being represented by 14 taxa, mainly polychaete annelids belonging to the Spionidae, Capitellidae, and Cirratulidae families, while the latter included 23 taxa, most of them being polychaetes belonging to different families of the order Phyllodocida (Phyllodocidae, Nereidae, Syllidae, Glyceridae, and Polynoidae). It also included decapod crustaceans, gastropods, and oligochaetes (all groups are detailed in Supplementary Figure 3). The relative frequencies of each group did not show a significant difference in the proportion of each functional group with latitude (Supplementary Figure 4).

Taxonomic β diversity values for macrofauna associated with honeycomb worm reefs showed greater similarity on average within regions (19–51%; Figure 5A and Table 4) compared to among regions (9–33%; Figure 5B and Table 4). However, levels of similarity within regions remained low, indicating important heterogeneity across sites of a given region. On average, when considering pairs of assemblages within regions, 60% of the species were found in only one assemblage: 50% of them changed in terms of species identity (turnover) and 10% were unique to the richest assemblage (nestedness) (Figure 5A – within region, Table 4). For pairwise comparisons among different regions, differences were even more pronounced, with an average of 80% of species being found in only one assemblage, with 70% due to species turnover and 10% linked to nestedness (Figure 5B – among bioregion, Table 4; for all pairwise comparisons among regions, see Supplementary Figure 5 and Table 4). The contributions of nestedness to β diversity were on average similar within and among regions. Overall, variation in species composition within and between bioregions were primarily due to changes in species identity.

Functional β diversity values for macrofauna associated with honeycomb worm reefs showed a comparable range in similarity within regions (38–88%; Figure 5A and Table 4) and among regions (34–84%; Figure 5B and Table 4). This similarity within and among regions, indicates high levels of overlap in functional space. On average, two assemblages shared 40% of their functional space, while functional β diversity was mostly driven by nestedness (i.e., by difference in the volume of the functional space filled by the assemblages; 24%) rather than by turnover (i.e., functional spaces not shared by the two assemblages; 15%) (Figures 5C,D and Table 4). The contributions of nestedness to functional β diversity were similar within and between regions (for all pairwise comparisons among bioregions, see Supplementary Figure 6).

Honeycomb worm reefs host diverse invertebrate assemblages. Here we examined how multiple facets of diversity, including taxonomic and functional α and β diversity, vary over much of the Atlantic coast of Europe. In terms of local levels of diversity, no significant differences were observed in Hill diversity indices (including richness) over the 10 study sites. Only the abundances of macrofauna were relatively higher in the southern sites compared to the northern sites. Our results are in agreement with a growing number of examples that show that there are many exceptions to the latitudinal diversity gradient described by Brown and Lomolino (1998) and Gaston and Chown (1999). Recent investigations have shown little or no relationship of diversity with latitude for the European marine benthos (Renaud et al., 2009; Hummel et al., 2017), particularly for soft sediment communities (Kendall and Aschan, 1993; Wilson et al., 1993; Kendall, 1996). Latitude is not a unidimensional environmental variable but a proxy for a number of primary environmental factors that interact and correlate with each other (Hawkins, 2003). For honeycomb worms, it appears that biotic and abiotic factors associated with the reef environment have contributed to maintaining constant levels of diversity over broad geographical scales, as discussed further below.

The assemblages sampled in our study showed a high diversity of macrofaunal organisms and are typical of the honeycomb worm reef assemblages reported in previous studies (Gruet, 1986; Dias and Paula, 2001; Dubois et al., 2002, 2006; Schlund et al., 2016). Mean species richness was comparable, but notably lower (7–16 species), than most previous reports for honeycomb worm reef richness (7–33 species; Holt et al., 1998; Dubois et al., 2002, 2006; Schlund et al., 2016; Jones et al., 2018). These values are also comparable to the values reported for epifauna from macroinvertebrate assemblages associated with Laminaria hyperborea (Gunnerus) Foslie, 1884 kelp holdfasts (10–30 species; Teagle et al., 2018) but lower than infaunal assemblages reported for the eelgrass Zostera marina Linnaeus, 1753 (10–60 species; Blanchet et al., 2004; Boström et al., 2006) in the northeast Atlantic. In terms of richness, maerl beds (composed of at least two species: Lithothamnion corallioides and Phymatolithon calcareum; Riosmena-Rodríguez et al., 2017) hosted, by far, the richest assemblages sampled thus far, with an average richness of 53 species (32–73 species; Boyé et al., 2019). Our results confirm that habitats engineered by honeycomb worms, as with other bioengineered habitats, host higher species richness than geogenic habitats (8–12 species; Jones et al., 2018; Boyé et al., 2019).

While species richness in assemblages associated with honeycomb worms remained within a narrow range throughout the coast of Europe, faunal composition did vary among sites. Two environmental variables were found to significantly structure assemblages: mean annual water temperature and the quantity of sediment in the cores. Mean annual temperature distinguished the United Kingdom sites and France sites from the Portugal sites and one site (FR4) in the Bay of Biscay. Thermal regimes affecting faunal composition appear to change within the Bay of Biscay, south of the Brittany peninsula, consistent with higher sea surface temperatures in the Bay of Biscay than in the surrounding areas connected by the Gulf Stream (Jenkins et al., 2008; Hummel et al., 2017). Sediment content was a key variable that structured communities, with sites that had higher sediment content, typically the France sites, being distinct from sites with lower sediment content, namely the United Kingdom and Portugal sites. Sediment content was negatively correlated with current velocity, such that sites that had lower hydrodynamics accumulated more sediment, while sites with higher hydrodynamics had higher porosity within the reef. In areas with high current velocities, the reef-building activity of S. alveolata is challenged by wave erosion, generating higher porosity within reefs. Conversely, low current velocities allow reefs to grow homogeneously but also allows unconsolidated particles to settle within fissures in the reefs. Previous studies have reported that within the same site, dense sections of reef, where S. alveolata is in an active growth phase (prograding reef, sensu Curd et al., 2019) tend to host assemblages with lower abundances and diversity than parts of the reef that are more fragmented (retrograding reef) (Dubois et al., 2002; Jones et al., 2018).

Sediment content has been found to be an important variable structuring communities at local scales. In the Bay of Mont-Saint-Michel, for example, higher sediment content in retrograding reefs explained the presence of many species typically belonging to muddy sandy bottom communities (Dubois et al., 2002, 2006). At the regional level, differences in hydrodynamic and sediment accumulation regimes may therefore lead to differences in reef density, which in turn affect assemblage compositions over the Atlantic coast of Europe. Unlike other engineering species such as haploops (gregarious tube-dwelling amphipods) (Rigolet et al., 2014), the density of S. alveolata was not the main factor structuring communities. Unlike haploops, the reef structures developed by S. alveolata persist after the death of the individuals, such that characteristics of the reef (here sediment content) better explain variation in communities than the density of the engineering species. Our results are consistent with other studies that have shown that reef structure is more important for explaining community composition than density of the engineer, such as in habitats built by Owenia fusiformis Delle Chiaje, 1841 (Fager, 1964); Spiochaetopterus bergensis Gitay, 1969 (Munksby et al., 2002; Hastings et al., 2007), and Lanice conchilega (Pallas, 1766) (Zühlke et al., 1998; Zühlke, 2001; Callaway, 2006; Rabaut et al., 2007; Van Hoey et al., 2008; De Smet et al., 2015).

We found five functional groups in association with honeycomb worm reef formations. Our results show that changes in taxonomic composition did not result in changes in ecological role, but rather that the same functional groups were found in association with honeycomb worm reefs throughout the coasts of Europe. Foundation species greatly influence the structure and functioning of species assemblages (Bruno et al., 2003). However, the effects of foundation species on biodiversity are not necessarily positive for all species, providing resources for some but excluding others (Rigolet et al., 2014). Through tube building activity, honeycomb worms transform unconsolidated sediment into a complex three-dimensional structure with properties that differ from both rocky shores or bare sediment. The reefs attract mostly soft sediment infauna (other polychaetes) and provide pockets of soft sediment for burrowers, but exclude taxa that require rocky substrate to settle upon or that compete for space with S. alveolata, such as barnacles and mussels (Holt et al., 1998; Dubois et al., 2002, 2006). Brown and red macroalgae cover on honeycomb worm reefs is reduced compared to rocky shores (Dubois et al., 2006), hence, excluding a large set of herbivores, and favoring deposit-feeders, as can be seen in the functional groups recovered here. Honeycomb worm reefs may therefore act as a biological filter for a given local pool of organisms. Similar results have been reported in different bioengineered habitats such as the communities associated with haploops in the bay of Concarneau, France (Rigolet et al., 2014). Contrary to adjacent sandy and muddy bottom communities, the establishment of haploops communities excluded or limited the colonization of other burrowers and tube-dwelling suspension-feeders, but attracted small mobile predators which possibly predate on haploops or other small associated organisms. The biological filtering that occurs in honeycomb worm and haploops habitats may therefore be applicable to other bioengineered habitats.

Our results show that while the taxonomic composition of the fauna associated with honeycomb worm reefs varied over broad geographical scales, species were replaced by another with the same functional role, such that only few differences in functional groups occurred across the reefs of the northeastern Atlantic. Despite high species turnover observed between biogeographic regions (70% on average), functional turnover was only 15% on average. As a result, high functional similarity was observed between regions and most functional changes across regions were due to one assemblage being a subset of another (24% of functional nestedness). For instance, the functional role performed by the isopod Lekanesphaera levii (Argano and Ponticelli, 1981) sampled in the northern reefs is the same as another isopod, Dynamene bidentata (Adams, 1800) in the southern reefs. This is also the case for two species of Phyllodocidae: Eulalia clavigera (Audouin and Milne Edwards, 1833) and Eulalia ornata (Saint-Joseph, 1888). E. ornata is more abundant at Boreal and Boreal-Lusitanian sites but tends to be replaced by E. clavigera at Lusitanian-Boreal and Lusitanian sites.

Biogenic habitats tend to show little variation in functional groups over broad spatial scales with high levels of redundancy within each group (Hewitt et al., 2008; Barnes and Hamylton, 2015), although this depends somewhat on the foundation species (Boyé et al., 2019). High species turnover accompanied by constant taxonomic richness has also been observed in eelgrass assemblages (Boyé et al., 2017). Similarly, variation in taxonomic composition associated with eelgrass beds, mangroves, maerl beds, and coral reefs did not result in differences in functional trait composition across approximately 500 km of coastline from both sides of the Atlantic and from the Caribbean and coral seas (Hemingson and Bellwood, 2018; Boyé et al., 2019). Our results support previous work that has shown that biogenic habitats are important in structuring benthic assemblages at the regional scale. And as with other biogenic habitats, communities supported by honeycomb worms show high functional redundancy that is thought to provide spatial insurance for benthic ecosystem functioning at local and broad spatial scales (Boyé et al., 2019). In addition, high functional redundancy may indicate that honeycomb worm reefs can offer similar niche properties to its associated assemblages across varying environmental conditions, as has been found for North Atlantic eelgrass, Caribbean mangroves or Indo-Pacific coral reefs (Cornell and Lawton, 1992; Boyé et al., 2017; Hemingson and Bellwood, 2018; Storch and Okie, 2019). Indeed, the reefs themselves provide protection from the physical elements, such as wind, waves, sun exposure, and desiccation, which may mean that they also act as environmental filters, buffering extremes in temperature or other environmental variables as is observed for many foundation species (Bertness and Callaway, 1994; Bruno et al., 2003; Bouma et al., 2009). As such, they may override the effect of large environmental gradients such as latitudinal temperature gradients (Jurgens and Gaylord, 2018), with important consequences for the spatial and temporal variation of their associated communities (Bulleri et al., 2018; Boyé et al., 2019).

Taxonomic differences were observed in macrobenthic faunal assemblages associated with honeycomb worm reefs along the coast of Europe. Hierarchical clustering analyses showed that assemblages were grouped together into four main clusters, two of which were found in the United Kingdom, a third represented by all assemblages in France, and a fourth represented the assemblages from Portugal. These clusters suggest there may be greater regional differences in macrobenthic communities than is currently recognized in biogeographic frameworks which consider only two main provinces in the northeastern Atlantic (Spalding et al., 2007; Briggs and Bowen, 2012). The assemblages observed in our study suggest important differences in species composition within both Lusitanian (Portugal vs. France sites) and the Boreal/Northern European Seas provinces (UK1 + UK4 vs. UK2 + UK3). While our results show good agreement with assemblage differences that could correspond to differences in species distributions associated with Lusitanian-Boreal and Lusitanian province subdivisions (defined by Dinter, 2001), it is less clear whether the northern assemblages found in honeycomb worms support a Boreal and Boreal-Lusitanian subdivision. Nevertheless, the differentiation observed here does suggest a higher degree of partitioning of species within the northeastern Atlantic that may be in part related to biogeography, with many species having restricted distributions. Our findings for the four United Kingdom study sites are in agreement with patterns of community composition associated with another important temperate foundation species, the kelp L. hyperborea (Teagle et al., 2018). While communities associated with holdfasts were fairly variable throughout the study area, six sites collected in northern and southern Scotland (Boreal province) were distinct from six sites in northwestern Wales and Southern England, supporting partitioning within the Boreal/Northern European Seas province.

In honeycomb worm assemblages, differences in Boreal and Boreal-Lusitanian regions were mainly driven by the occurrence or absence of M. edulis and potential hybrids of M. edulis × M. galloprovincialis in the two studied regions. Given the morphological similarities among the three taxa, identification to the species level was uncertain in UK1 and UK4, with most identifications being kept to the genus level. This may partly explain the similarity among the Scottish and southern England sites in our dataset. Molecular identification was not attempted here, but could help resolve some of these taxonomic uncertainties in future work, as has been done for Mytilus spp. (Wenne et al., 2020). Similarly, nematodes were identified only to the phylum level, but did contribute significantly to the observed differences in assemblages between the Lusitanian and Lusitanian-Boreal regions. Identification to the species level would likely further differentiate these two regions (Bhadury et al., 2006). Overall, our results indicate strong differences in community composition that are related to taxonomic turnover, which in turn indicates that species ranges may better correspond with finer-scale biogeographic partitioning within the northeastern Atlantic, as proposed by Dinter (2001). Future studies that take into consideration more extensive species inventories across phyla and habitats may help refine our understanding of marine biogeography in the northeastern Atlantic, and resolve some of the current discord in various frameworks.