Musical Chairs on Temperate Reefs: Species Turnover and Replacement Within Functional Groups Explain Regional Diversity Variation in Assemblages Associated With Honeycomb Worms

Introduction

Unraveling the processes that control how biodiversity is distributed over space and time are central objectives in macroecology and biogeography (Addo-Bediako et al., 2000; Ricklefs, 2004). While biodiversity was first observed to increase from the poles toward the tropics over 200 years ago (Hawkins, 2001), it is now recognized that latitude per se is not the main driver of spatial gradients in biodiversity, but rather a combination of variables and mechanisms that include ecological, evolutionary and historical processes (Hawkins and Diniz-Filho, 2004; D’Amen et al., 2017). For example, the rates of species diversification are thought to be greater in the tropics due to higher mutation rates in warmer regions (Rohde, 1992; Mittelbach et al., 2007). In addition, regions that have greater environmental stability over time may also host higher species diversity than areas that have suffered major environmental change, such as glaciations (Fine, 2015; but see also Fordham et al., 2019). Other variables shown to influence biodiversity include habitat or organism types, and organism properties, such as biomass, dispersal rates and physiology (Addo-Bediako et al., 2000; Hillebrand, 2004; Buckley et al., 2010; Rolland et al., 2015; Gaucherel et al., 2018). Given the multitude of factors that can affect diversity, understanding and predicting spatial and temporal patterns in biodiversity is challenging.

An important factor that influences biodiversity is whether a community is found on geogenic (of geological origin: sedimentary or rocky) or biogenic (of biological origin) substrate. Many communities co-exist in habitats that are altered by another living organism, and these are called foundation species or ecosystem engineers. Foundation species not only build habitat (Dayton, 1972) but also control the availability of resources for other organisms (Sarà, 1986; Jones et al., 1994). Through habitat modification, foundation species can alter the realized niche of species, at times facilitating niche expansion by buffering environmental conditions so that they continue to be favorable within the engineered habitat (Bulleri et al., 2016). Heterogeneity in the engineered habitat can promote facilitation of a greater number of species, while dominance of the engineering species or homogeneity in the habitat can enhance competition and limit the number of associated taxa (Schöb et al., 2012; Bulleri et al., 2016). Given that community composition can vary greatly within biogenic habitats across environmental gradients (Boström et al., 2006; Boyé et al., 2017), investigating within-habitat diversity is essential for guiding conservation actions (Airoldi et al., 2008) and enhances our understanding of biodiversity over broad geographical scales.

Although the northeast Atlantic is amongst the most studied marine regions on Earth (Hawkins et al., 2019), spatial structure in its coastal marine assemblages remains poorly understood. Variation in diversity over broad spatial scales may be related to species distributions, but even broad biogeographic delimitations continue to be contentious. For example, broadly accepted marine biogeographic frameworks consider two biogeographic provinces in the northeastern Atlantic: Boreal and Lusitanian (Briggs and Bowen, 2012), or Northern European Seas and Lusitania (Spalding et al., 2007), separated by “Forbes’ Line” (sensu Firth et al., 2021, after Forbes and Godwin-Austen, 1859). They differ, however, on boundary positions between provinces, and whether or not the southern province includes the Mediterranean Sea. Further biogeographic subdivision has been proposed for the northeastern Atlantic such that four provinces could be recognized: Boreal, Boreal-Lusitanian, Lusitanian-Boreal, and Lusitanian (Dinter, 2001), but these finer-scale subdivisions are less often considered or employed in macroecology. Lack of consensus on biogeographic delimitations are partially due to competing criteria used for setting boundaries but may also reflect incomplete distributional knowledge of many marine species, especially for poorly studied invertebrates. If the majority of species have restricted distributions, then species turnover might be expected to be higher over a given spatial scale. Community structure may therefore be related to some extent to biogeographic partitioning, such that communities may be more similar within rather than between biogeographic regions.

To examine diversity in marine ecosystems, it is important to consider how diversity is quantified and described. Biodiversity is a multifaceted concept that includes several components (Whittaker, 1972). In order to better understand what mechanisms influence biodiversity, it may be helpful to consider each of these different facets. While local patterns in diversity (α diversity; Whittaker, 1972) are most commonly assessed, regional differences in diversity due to variation in richness or species composition (β diversity; Whittaker, 1972; Airoldi et al., 2008) can also provide important insights into the mechanisms driving community structure (Hewitt et al., 2005; Anderson et al., 2011; Villéger et al., 2013). Focusing on β diversity is especially important in the context of global change, where ecological communities are subject to large environmental fluctuations and disturbances (Mori et al., 2018). Furthermore, it is now widely recognized that the integration of functional information based on species traits provides a better understanding of community functioning (Díaz and Cabido, 2001; Anderson et al., 2011; Pavoine and Bonsall, 2011; Münkemüller et al., 2012; Mouillot et al., 2014). Thus, comparing taxonomic with functional diversity (α and β) provides a better understanding of the ecological processes that shape community composition (Swenson et al., 2011; Villéger et al., 2013; Mori et al., 2018) and the impact of biodiversity loss on ecosystem functioning (Cadotte et al., 2011; Burley et al., 2016). For example, selective processes, such as environmental filtering, lead to homogenization of traits in communities, since only species with a specific set of traits could survive and develop under certain abiotic conditions. As a result, the loss of species with unique functional characteristics may have significant consequences for ecosystem functioning than the loss of species with characteristics that are more commonly expressed in the community (O’Connor and Crowe, 2005; Queirós et al., 2013). Nevertheless, the comparison of taxonomic and functional β diversity alone may not reveal the underlying ecological processes that structure communities (Baselga, 2010; Villéger et al., 2013; Legendre, 2014). This is in part because variation in species composition among sites (β diversity) is the resultant of two components: species turnover (i.e., replacement of species or functional strategies) and nestedness (i.e., dissimilarity associated with the loss of species or functional strategies, in which an assemblage is a strict subset of another). Partitioning β diversity into turnover and nestedness thus provides an additional facet for dissecting community assembly rules. In sum, a combination of tools and metrics, including taxonomic and functional α and β diversity (and their components) are essential for better understanding biodiversity in marine ecosystems.

The honeycomb worm, Sabellaria alveolata (Linnaeus, 1767), is a physical ecosystem engineer (Berke, 2010) commonly found along the European coast from northwest Ireland to northern Mauritania (Curd et al., 2020), where it builds biogenic structures of varying extent in the intertidal and shallow subtidal zones. Honeycomb worms build what are considered Europe’s largest biogenic reefs (Noernberg et al., 2010) and support a unique and rich assemblage of species (Dias and Paula, 2001; Dubois et al., 2006; Jones et al., 2018). Honeycomb worms play key functional roles in the ecosystems they support, by creating new three-dimensional habitat, which increases the physical complexity of the initial substrate, increases local biodiversity (Dubois et al., 2006; Jones et al., 2018), limits coastal erosion (Noernberg et al., 2010) and fashions biogenic structures (ranging from crusts and veneers to large reefs, hereafter “reefs” for simplicity) with high esthetic and recreational fishing value (Plicanti et al., 2016). Honeycomb worm reefs are broadly distributed across temperate Europe, however diversity investigations have only been carried out at local scales (Dias and Paula, 2001; Dubois et al., 2002, 2006; Schlund et al., 2016; Jones et al., 2018). It is therefore currently unknown how biodiversity supported by these reefs varies over its range.

The present study examined the patterns of diversity of benthic marine macrofauna associated with honeycomb worm reefs from sites spanning the entire European distribution of the species (but excluding North Africa), in order to address the following questions: (i) Does taxonomic and functional diversity of communities associated with honeycomb worms vary over broad geographical scales, and if so, what environmental drivers best explain this variation? (ii) Does community composition within honeycomb worm reefs vary with respect to currently described biogeographic provinces? (iii) Are there regional differences in taxonomic and functional β diversity in assemblages associated with honeycomb worm reefs? (iv) If so, are they mainly due to differences in species richness or in turnover? Finally, can differences between taxonomic and functional diversity help identify the ecological processes that affect biodiversity on honeycomb worm reefs?

Materials and Methods

Study Area and Sampling Methods

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.

Figure 1. (A) Map indicating the locations of the 10 study sites in the United Kingdom, France, and Portugal within the four biogeographic provinces defined by Dinter (2001). (Eight stations were sampled in each of the ten sites, except for UK4 where five stations were sampled.) (B) Abundance (±SD) of the associated macrofauna (in number of individuals.m^–2) and taxonomic indices (Richness, Hill diversity indices N1 and N2) per site.

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).

Data Processing and Environmental Variables

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).

Taxonomic Diversity

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).

Functional Diversity

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.

Table 1. Selected biological traits and modalities for analysis of biological traits (BTA).

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.

Partitioning of Taxonomic and Functional β Diversity

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.

Results

Taxonomic Diversity

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).

Table 2. Results of PERMANOVA to test for differences in benthic macrofaunal assemblages between bioregions (fixed) and sites (random, nested within bioregion).

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).

Figure 2. (A) nMDS representing the variability across sites and stations within sites. (B) Cluster constructed using the “Ward D2”’ method, showing the four groups of replicates. Height indicates the order in which the clusters were joined. (C) nMDS constructed from data averaged for each site, grouped by cluster groups. A–C were derived from a Bray–Curtis dissimilarity matrix constructed from log (x + 1) abundance data transformed.

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.

Figure 3. Result of the distance-based redundancy analysis (dbRDA) using the Bray–Curtis dissimilarity matrix computed on the log-transformed abundance data and the seven selected environmental variables. Field data: “Water” and “Air (T°C max)” = maximum air and water temperature, “Sediment” = amount of sediment in g.m^–2, “Density SA” = density of S. alveolata (ind.m^–2). Data extracted from BioORACLE: “Chla” = average chlorophyll level (mg.m^–3), “Water (T°C avg.)” and “Current (max)” = average maximum current velocity (m.s^–1). The first two axes capture 22.9% of the variation explained by these seven variables.

Table 3. DISTLM marginal test results for each environmental predictor variable selected for the most parsimonious model for macrofaunal assemblages.

Functional Diversity

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).

Figure 4. (A) Cluster constructed using the “Ward D2”’ method performed from the Fuzzy Correspondence Analysis (FCA) results. (B) Projection of modalities and taxa by biological traits on the first two axes of the FCA. (C) Histogram of eigenvalues. The first two axes (in black) were kept in the FCA.

Taxonomic and Functional β Diversity

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.

Figure 5. Triangular plots illustrating the geographical pattern of (A,B) taxonomic and (C,D) functional β diversity Sorensen dissimilarity between the species composition (presence/absence data) of the 10 study sites was used to quantify their similarity, and the two components of their beta diversity nestedness (i.e., influenced by the difference in number of species between the two communities) and turnover (i.e., species replacement between two communities). Contributions were calculated for comparisons between samples belonging either to (A,C) the same bioregion (within bioregion) or to (B,D) different bioregions (among bioregion). Red lines indicate the centroid value for each graph with its associated mean values for the three components of the Sorensen dissimilarity.

Table 4. Taxonomic and functional β diversity and its two components; mean ± standard deviation.

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).

Discussion

Influence of Honeycomb Worm Reefs on Local Diversity

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.

Relationship Between Taxonomic and Functional β Diversity

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).

Biogeographic Regions

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.

Conclusion

Our results highlight the importance of considering various aspects of diversity in order to have a more comprehensive understanding of the ecological processes that shape marine communities, which will ultimately better inform conservation strategies (Meynard et al., 2011; Villéger et al., 2013; Loiseau et al., 2017). In the case of honeycomb worm reefs, the environmental filtering that excludes some functional groups from the reefs would not have been detected if only patterns in taxonomic diversity had been examined. Similarly, the high level of taxonomic turnover observed among and within biogeographic regions would have been overlooked if functional diversity had been considered alone. Examining multiple facets of community diversity has enhanced our understanding of the factors that shape assemblages associated with S. alveolata. Preserving taxonomic diversity is and will continue to be valuable for maintaining regional levels of species diversity (De Juan and Hewitt, 2011). However, while many conservation programs prioritize the conservation of local taxonomic community diversity (Socolar et al., 2016), considering the functional complementarity of communities across broader spatial scales in addition may prove to be more efficient in maintaining healthy ecosystems (Mori et al., 2018). The results presented here show that honeycomb worm reefs support high taxonomic diversity, but functional diversity is more limited, with some key functional groups (such as grazers) being absent from the community. Adjacent habitats may therefore host very different sets of species with equally important ecological roles. Benthic homogenization and loss of complexity of the sea floor may reduce overall functional diversity in the marine benthos (Airoldi et al., 2008). Protecting a diversity of benthic habitats may therefore be necessary for ensuring good ecosystem functioning in marine communities.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. The functional trait dataset for this study can be found in SEANOE (https://www.seanoe.org/; https://doi.org/10.17882/79817). Further inquiries can be directed to the corresponding author/s.

Author Contributions

SD, LF, and FN conceived the study. AM led analysis and drafting of the manuscript. LB, AC, AD, SD, LF, FL, CM, FN, and RS carried out the field work. CP, AM, SD, and ND conducted species and trait identification. AB, CP, MM, and SD contributed to data analysis. All authors contributed to editing the final manuscript.

Funding

This work was supported by the Total Foundation (Grant No. 1512 215 588/F, 2015). AM and AC were funded by Ph.D. grants from Ifremer, Region Bretagne, and the ISblue project (the Interdisciplinary graduate school for the blue planet, ANR-17-EURE-0015, co-funded by the “Investissements d’Avenir” program).

Conflict of Interest

LB was employed by the company Fugro GB Marine Limited.

The remaining 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.

Acknowledgments

The authors would like to thank Céline Cordier, Mickael Vasquez, and Anthony Knights for assistance with field collections. The authors would also like to thank Maïwenn Lescop, Gautier Sourimant, and Jean-Dominique Gaffet for their help with sorting through reef samples and measuring worms. This study forms part of the REEHAB (REEf HABitat) project (www.honeycombworms.org).

Supplementary Material

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

Footnotes

1. ^ http://www.marinespecies.org/index.php (accessed January 15, 2021)
2. ^ https://www.bio-oracle.org/
3. ^ www.marlin.ac.uk/biotic

References