Coupled effects of environment, space and ecological engineering on seafloor beta‐diversity

966 –––––––––––––––––––––––––––––––––––––––– © 2021 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Malin Pinsky Editor-in-Chief: Jens-Christian Svenning Accepted 31 January 2021 44: 966–974, 2021 doi: 10.1111/ecog.05440 44 966–974


Introduction
Human activities are transforming the biosphere and in the last 180 years ecosystem degradation and species extinctions are accelerating in an unprecedented way (Hoegh-Guldberg and Bruno 2010, Dirzo et al. 2014, McGill et al. 2015. Understanding the dynamics of ecological communities is essential for predicting the effects of homogenization and biodiversity loss in ecosystems (de Juan et al. 2013, Blowes et al. 2019 Eriksson and Hillebrand 2019). β-diversity is at the core of metacommunity theory as it describes variability in diversity across landscapes compared to the local site diversity (i.e. the ratio between α and ɣ diversity). The current view on metacommunity dynamics states that both environmental (i.e. niche) and spatial (i.e. neutral) factors act together, and the interdependency between these processes shape communities (Leibold et al. 2004, Brown et al. 2017. However, the effects of environmental heterogeneity on niche construction are often mixed with the biotic effects (Kraft et al. 2015, Thakur andWright 2017), and the role of species interactions are frequently neglected in community assembly models (Thakur and Wright 2017, Maynard et al. 2020). Here, we examine the relative roles (and their interdependencies) of environmental variability, spatial constraints and biotic interactions on the community assembly at the seafloor.
Biotic interactions have proven to be pivotal to determine species coexistence and niche construction in heterogeneous seascapes (Lohrer et al. 2013, Aller and Cochran 2019. Species can persist and modify the environment, thus affecting community assembly through biotic engineering, facilitation and competition (Thakur andWright 2017, Mod et al. 2020). Trait diversity and the degree of functional overlap within communities can tell us the importance of biotic interactions for the assembling process (Mouchet et al. 2013, Ulrich et al. 2018, Mod et al. 2020. For example, low diversity and high functional overlap might suggest a predominance of environmental filtering, on the other hand, low functional overlap could indicate strong competition and exclusion of species that display similar functions (Hewitt et al. 2008, Ingram and Shurin 2009, Mason et al. 2011). In addition, the effects of species competition tend to decrease with the increase in geographical distance while facilitation remains constant (Hewitt et al. 2008, Mod et al. 2020. Despite the proven importance of these biotic interactions, environmental drivers over multiple spatial scales and dispersal constraints are often considered the primary drivers of local community assembly in marine ecosystems (Kraan et al. 2015, Menegotto et al. 2019. However, the relative importance of environmental and spatial drivers in structuring communities in aquatic ecosystems also depend on the species dispersal ability and their degree of generalism (Heino 2013, Heino et al. 2015, Li et al. 2020. Generalists species and passive dispersers are more affected by spatial processes, e.g. spatial distance and connectivity among similar habitat types, whereas specialists are influenced by local environmental controls, e.g. small-scale variability in sediment grain size, organic matter and turbidity (Rodil et al. 2017). Hence, to fully understand the spatiotemporal variation in biodiversity, we must integrate biotic interactions into current methods for characterising community reassembly (e.g. β-diversity distance decay models), which are often based only on spatial and environmental distance.
It is widely recognized that β-diversity is a simple function of αand ɣ-diversity, although there are several ways to measure it (Kraft et al. 2011, Legendre andDe Cáceres 2013).
The additive partitioning of dissimilarities (Baselga 2012, Baselga andLeprieur 2015) is a multivariate distance-based method that assesses the contribution of β-diversity components. It breaks overall dissimilarity in samples into species turnover and nestedness, i.e. replacement and richness differences (Legendre 2014). Subsequently, these indexes can be combined with multivariate ordinations and variance partitioning to evaluate the relative contribution of environment, space and biotic interactions in community assembly, thereby providing a flexible way to investigate the influence of niche and neutral processes on community richness and composition (Ruhí et al. 2017, Boyé et al. 2019. Landscape heterogeneity is a key element in structuring biodiversity with biogenic habitats (e.g. shellfish or seagrass beds) and thus an important driver of environmental variability in aquatic ecosystems (de Juan et al. 2013, Boyé et al. 2019. Habitat filtering associated with biogenic structure, promotes β-diversity and variability in community structure across landscape (Hewitt et al. 2008, Lohrer et al. 2013, Boyé et al. 2017 leading to gradients of species replacement or richness differences in biodiversity (Legendre 2014, Boyé et al. 2017. Limited geographical connectivity and variability in species dispersal can also affect community structure at multiple scales (Lundquist et al. 2004, Moritz et al. 2013). Dispersal limitation is expected to increase with increasing geographical distance between locations (Heino et al. 2015). Therefore, coupled effects of both niche and neutral processes are likely to affect β-diversity, since spatially structured environmental gradients are significant drivers of variability in taxa replacement and richness differences.
Seafloor ecosystems are a good case study for showing these processes since they are environmentally and biotically heterogeneous, species-rich, multi-trophic, influenced by biotic interactions (such as facilitation and ecosystem engineering), highly connected and cover an area greater than all other habitats on earth combined (Solan et al. 2020). In this study, we combined information on invertebrate taxa occurrences from three different coastal soft-bottom seascapes, each one containing 400 observations sampled over 300 000 m 2 . We applied the additive partitioning of community dissimilarities to evaluate the relationship between β-diversity components with environmental properties, spatial constraints and biotic engineering in seafloor ecosystems. We used distance-based redundancy analysis (dbRDA) and variance partitioning to quantify the relative importance of these predictors for the variation in species richness and composition. We expected that: 1) variability in taxa replacement and richness differences would be associated with environmental variability (e.g. characterised by sediment properties), indicating environmental control on community assembly; 2) taxa replacement related to variations in the occurrence and abundances of ecosystem engineers, signal the important role of biotic interactions in community assembly; 3) largescale spatial constraints (variation among seascapes) will be more important than small-scale variability for β-diversity, due to the constraints of spatial connectivity in the community assembling process (species occurrences). Additionally, interactions between these factors would highlight the interdependencies of space, environment and biology in driving community assembly.

Dataset
The dataset consists of a total of 1197 sediment and macrofauna samples from 3 seascapes, containing 146 macroinvertebrate taxa and 13 environmental variables (all data available via Kraan et al. 2019Kraan et al. , 2020a. Samples were collected from three intertidal sandflats in Kaipara (174°17′S, 36°23′E), Manukau (174°41′S, 37°7′E) and Tauranga (175°56′S, 37°27′E), North Island, New Zealand ( Fig. 1) in the austral summer of 2012 (Kraan et al. 2015, Greenfield et al. 2016. Within each seascape, 400 cores for macrofauna (13 cm diam., 20 cm deep) were taken along 3 transects of 1 km length, established 100 m apart from each other. The sampling design was used to allow comparisons across different spatial scales, in which sampling points were spaced along three 1 km long transects repeating a sequences of inter-sample distances of 30 cm, 1, 5, 10, 30, 50, 100, 500 and 1000 m. The macrofauna cores were sieved on a 500 µm mesh sieve and macrofauna preserved in 70% isopropyl alcohol for later taxonomic identification (usually to species level). The percentage of seagrass Zostera mulleri, bare sand and shell hash cover was estimated through photographs (0.25 m 2 ) of each sampling point (Kraan et al. 2015, 2020a for detailed information). Mean sediment grain size and grain size fractions (silt < 63 μm, very fine 63-125 μm, fine 125-250 μm, medium 250-500 μm and coarse > 500 μm), as well as organic content (%), chlorophyll-a and phaeopigments (mg g −1 ) were estimated from a pool of three surface sediment cores (2 cm diam., 2 cm deep) at each sampling point. Grain size were measured using a Malvern Mastersizer, pigment concentrations were determined using a fluorometer, and loss on ignition was used to assess organic content (Kraan et al. 2015(Kraan et al. , 2020a. Missing data from the environmental factors (consisting of 2.4% of the data from the seagrass, bare sand and shell hash coverage) were estimated using the predictive mean matching method in the R package 'mice' (van Buuren and Groothuis-Oudshoorn 2011).
All three sites had a median grain size classified as fine sands (Kaipara = 213 µm, Tauranga = 197 µm, Manukau = 166 µm), yet the three locations differed in terms of environmental characteristics. The Manukau site is the muddiest (silt = 14%) with the highest concentration of chlorophyll-a (23 mg g −1 ) but also had the highest percentages of shell hash cover (16%). The Tauranga site has the highest proportion of coarse sediments (6%) and contains the highest coverage of seagrass (23%), whereas the Kaipara site is composed mostly of bare sands (84%) (Kraan et al. 2015). The ecosystem engineer bivalves, Macomona liliana and Austrovenus stutchburyi were abundant and widespread in the three locations. In Manukau an average of 38 ± 23 ind m −2 (mean ± SD) of M. liliana and 85 ± 115 ind m −2 of A. stutchburyi were found. Whereas in Tauranga, densities of M. liliana 31 ± 15 ind m −2 and A. stutchburyi 23 ± 46 ind m −2 were lower. In Kaipara, the average density of M. liliana was 38 ± 31 ind m −2 and of A. stutchburyi was 15 ± 38 ind m −2 .
To explore the role of biotic interactions, the abundances of M. liliana and A. stutchburyi and the percentage of seagrass Zostera mulleri and shell hash cover (which is influenced by infaunal bivalve populations) were also used as covariates in the subsequent analyses. We chose these ecosystem engineers because they are widely distributed on intertidal sandflats, have strong density gradients, with density-dependent effects on the ecosystem's biochemical and physical environment, and the recruitment of juveniles into patches (Turner et al. 1997, Sandwell et al. 2009, Kraan et al. 2020b).

Data analysis
We use additive partitioning of β-diversity derived from species occurrences and computed the total Sorensen dissimilarity between samples (βSor) and its two components: 1) turnover or taxa replacement (βSim), and 2) nestedness or richness differences (βNes) using the R package 'betapart' (Baselga et al. 2018). Biotic data were Hellinger-transformed (Legendre and Gallagher 2001) and environmental data were standardized. In order to account for the spatial effects on community assembly the geographical location of each sample was converted to UTM and then used to compute the weighted principal coordinates of neighbour matrices (PCNM, Borcard and Legendre 2002). PCNM scores were detrended and to avoid overparameterization, we performed a stepwise backward selection and the axes that showed significant explanatory power (PCNM axes 1-10) were used as predictors of the effects of spatial distance on beta-diversity within seascapes. The seascapes were converted to dummy factors and used as categorical descriptors of the largest spatial scale.
Subsequently, the relationships between βSor, βSim and βNes with environmental, spatial and biotic engineering drivers were investigated using distance-based redundancy analysis (dbRDA) followed by variance partitioning Anderson 1999, Peres-Neto et al. 2006). dbRDA is a multifactorial linear model that allows us to test the direct relations of multivariate response data with multiple explanatory factors. Through this technique, we can evaluate the amount of variation explained by single factors or groups of factors via partial canonical regression, i.e. variance partitioning (Legendre 2008, Dray et al. 2012. To avoid collinearity, highly correlated variables (p < 0.01) were removed and only the variables that most contributed to explain the constrained variability were included in the final dbRDA models after stepwise backward selection.
The spatial variation in β-diversity at regional scale (among seascapes) was evaluated through dbRDAs of the βSor, βSim and βNes based on the three seascapes together, being Kaipara and Manukau west coast, and Tauranga east coast. For simplicity, we choose to show only the dbRDA bi-plot based of overall Sorensen dissimilarity (βSor) among seascapes. Additionally, dbRDAs ordinations of βSor, βSim and βNes from each seascape were also conducted. The contribution of taxa replacement (βSim) and richness differences (βNes) components is shown as the relative change in the point size within the bi-plot ordination. Variance partitioning was interpreted hierarchically based on the influence of biotic engineers and their shared effects with environment and space. Then the influence of the environment and its shared effects with space, and later the effects of spatial constraints alone. The Venn diagrams showing the variance partitioning of all ordinations are presented. dbRDAs, variance partitioning and PCNM analyses were performed with the R package 'vegan' (Oksanen et al. 2019).

Results
The dbRDA ordination explained 44% of the variability in β-diversity at the regional scale (βSor), with the first two axes of the dbRDA accounting for 75.5% of the constrained variation. The first axis separated sites composed mostly of fine sand, lack of seagrass cover and high abundances of the bivalve Macomona liliana from those with coarse sand, high percentages of seagrass cover and higher content of labile food (i.e. high values of chl-a/phaeo ratio). The second axis summarized a gradient between sites with high percentage of medium sands to those characterized by very fine sediments and higher organic matter, chlorophyll-a and phytodetritus concentrations ( Fig. 2a-b). Large-scale effects (i.e. among seascape variability) were the most important drivers of β-diversity and shared effects between biotic engineering, environment and space had higher importance than factors in isolation (Fig. 3a). Environment effects summed (the sum of all components where the environment is involved except the ones shared with biotic engineering) explained 22% of the variance in βSor and environmental variability among seascapes alone contributed to 16% (Fig. 3a). Effects of biotic engineering summed (including shared effects with environment and space) explained another 16% of the variance in βSor (Fig. 3a). Figure 2. dbRDA ordination of the Sorensen dissimilarities (βSor) of macrobenthic communities. Bubble size represent variations in species turnover βSim (a) and richness differences βNes (b). Sites from distinct estuaries are shown in different colours, Kaipara (dark grey), Manukau (medium grey) and Tauranga (light grey); n = 1197. Environmental variables are show in solid red lines and biotic engineering variables in blue doted lines; Eigenvectors' abbreviations: percentages of medium sand (ms), fine sand (fs), very fine sand (vfs) and silt (silt). Organic matter (om), chlorophyll-a (chl-a) and phaeopigment (phaeo) concentrations. Chlorophyll-a/phaeopigment ratio (chl-a/phaeo). Percentages of seagrass (seagrass) and shell hash (shell) cover. Abundances of M. Liliana (Mac) and A. stutchburyi (Aus). The percentage of constrained variation explained by each of the first two RDA axes are also shown. Taxa replacement (βSim) was the component that most contributed to β-diversity at the regional scale with high replacement occurring within the Manukau (indicated by the larger medium grey bubbles in Fig. 2a) at sites characterized by higher percentages of shell-hash cover, fine sediments and higher organic matter, chlorophyll-a and phytodetritus concentrations (Fig. 2a). Lower replacement was observed in the Tauranga seascape (indicated by the smaller size of the light grey bubbles in Fig. 2a), associated with higher percentages of seagrass cover, and higher content of labile food, i.e. chlorophyll-a/phytodetritus ratio (Fig. 2a). The first two axes of the dbRDA ordination explained 52% of the variability on βSim (Fig. 3b). The sum of environmental effects accounted for 26%, being 20% related to between-seascape (large-scale) environmental variations. The sum of biotic engineering effects explained another 19% of the variance in βSim (Fig. 3b). The gradient of richness difference (βNes) had a low contribution to β-diversity and the first two axes of the dbRDA explained only 3% of the variability in βNes at the regional scale.
High values of richness differences (βNes) were observed in sites of Manukau and Kaipara seascapes associated with high percentages of fine sands and higher abundances of M. liliana (shown by the variable sized bubbles for Manukau and Kaipara in Fig. 2b). Despite this, βNes explained a minor fraction of variability in biodiversity across sites and variation in environment, space and engineering variables accounted for only 3% of the variability in βNes (Fig. 3c).
Kaipara and Manukau showed similar patterns of within seascape variation (Fig. 4). Environmental drivers explained a higher proportion of the variability in β-diversity at these seascapes. Environmental effects summed explained 15% of the variation in βSor at both Kaipara and Manukau (Fig. 4). Similarly, 19% of the variability in βSim was explained by the sum of environmental effects at Kaipara and 18% at Manukau (Fig. 4). In Kaipara, the sum of engineering effects explained another 11% of the variability in βSor and 14% in βSim (Fig. 4). In Manukau, biotic engineering effects summed explained 10% of the variance in βSor and 13% in βSim (Fig. 4). In Tauranga, located on the east coast, environmental drivers, biotic engineering and small-scale (withinseascape) spatial constraints explained similar proportions of the variance in both, βSor and βSim (Fig. 4).

Discussion
We found that engineering, environment and spatial factors in isolation explained a low amount of the variation in β-diversity. Biotic engineering is tightly coupled to spatial and environmental variation, and spatially structured abiotic and biotic effects were the main contributors of β-diversity (βSor) among seascapes. The β-diversity differences among seascapes were mainly explained by species turnover (βSim), whereas richness differences (βNes) among sites were less important. The abundance of the bivalve M. liliana had a positive effect on turnover, as observed by the dbRDA. This indicates that the spatial structure of these bivalves can modify the environment, indirectly structuring local communities with distinct composition from the background seascape.
More importantly, the various determinants of β-diversity did not simply operate at specific spatial scales. Therefore our initial expectations that 1) turnover (βSim) is influenced by environmental drivers; 2) variations in the distribution of ecosystem engineers contribute to β-diversity; and 3) large-scale spatial variation is more important to β-diversity than smallscale variability, could be accepted. Importantly, the specific role of each component cannot be disentangled from the others. Community structure within heterogeneous seascapes is dependent on direct and indirect effects that operate over multiple spatiotemporal scales. Environmental variation has a key role in this dynamic, however our results stress that spatial constraints (neutral processes) and biotic interactions must be considered to fully understand the mechanisms of community assembly. The current view on metacommunity dynamics states that both space and environment act together, and the interdependency between these processes shape communities (Brown et al. 2017  gradients were important drivers of species turnover, suggesting that spatial autocorrelation affect seafloor biodiversity. Environmental heterogeneity and community connectivity are coupled and influence communities over multiples scales (Thrush et al. 2005, Kraan et al. 2015, Menegotto et al. 2019. In our study, most of the β-diversity was captured by large-scale variability (among seascapes) whereas small-scale variability (within-seascapes) only explained a small portion of the variance in species composition. In addition, west coast seascapes (Kaipara and Manukau) showed more similarity in community composition compared to Tauranga, which is located in the east coast. This pattern could be related to differences in geomorphology and in the regime of waves and currents among coasts, but also to limited connectivity among the ecoregions (Ross et al. 2009). Therefore, neutral processes such as dispersal limitation and habitat connectivity are important drivers of seafloor β-diversity (Thrush et al. 2008). Hence, species life-history and dispersion potential, as well as disturbance-diversity relationships and priority effects, might be important for taxa turnover in marine ecosystems.
Environmental filtering is one of the most debated concepts in modern ecology (Kraft et al. 2011, Thakur and Wright 2017, since it can determine the prevailing mechanisms that drive community assembly, i.e. niche or neutral processes (Leibold et al. 2004, Rosindell et al. 2011, Brown et al. 2017. We showed that in marine soft-bottoms, environmental heterogeneity is a key driver of β-diversity, although, strongly interconnected with biotic engineering. Therefore, the mechanisms influencing species turnover are controlled by direct and indirect relationships of foundation species with environmental gradients (Boyé et al. 2019, Kraan et al. 2020b) and the resulting assembly pattern is taxa dependent (Boyé et al. 2017, Henseler et al. 2019. Hence, communities inhabiting heterogeneous landscapes are shaped by multiple factors, including natural and human-driven disturbances, trophic interactions, deterministic niche factors and stochasticity (Kraft et al. 2011, Mori et al. 2018, Limberger et al. 2019. For example, species dispersal is strongly dependent on hydrodynamic regimes and species life-history traits (Lundquist et al. 2004, Rodil et al. 2017 and can structure seafloor communities over multiple scales (Lundquist et al. 2004, Moritz et al. 2013). The effects of species competition are scale-dependent and are more significant on the small-scale (Mod et al. 2020). On the other hand, the effects of ecosystem engineering on community assembly are likely to propagate across different scales, which highlights the importance of seascape heterogeneity for maintaining β-diversity and ecosystem functions at larger scales (Hewitt et al. 2008, Mod et al. 2020. We found that biotic interactions significantly contributed to spatial turnover in species composition. Bivalves can colonize large extents of the seafloor, changing sediment biogeochemistry, food availability and habitat complexity (Lohrer et al. 2013, Thrush et al. 2017. Furthermore, they have a long legacy-effect in the environment, since calcium carbonate shells last longer than the life span of the organism and contribute to habitat heterogeneity (Hewitt et al. 2005). This has important implications for habitat complexity in a changing world (Hewitt et al. 2005, Blowes et al. 2019. For instance, regional and local decreases in population abundances of these habitat-forming species may boost biotic homogenization at the seafloor (Boyé et al. 2019). Moreover, differences between range distribution and abundance of engineer species may have implications for community resilience to disturbances and for biodiversity management (Greenfield et al. 2016, Gladstone-Gallagher et al. 2019. Macomona liliana is the most widespread species on Kaipara and Manukau seascapes, whereas Austrovenus stutchburyi is the most abundant at Manukau (Kraan et al. 2015). Hence, biotic engineers that are widespread across the seascape can be more resilient to small-extent disturbances, while species that occurred in smaller patches but in higher densities can be more sensitive.
Sites in Tauranga on the contrary were characterized by higher percentages of seagrass cover and lower abundances of bivalves. The benthic community was also more homogeneous, showing a lower β-diversity at site scale than the other two locations. This might suggest that seagrass meadows could provide a distinct habitat in the seascape, whereas bivalve engineers tends to promote small-scale heterogeneity in the sand habitat due to their bioturbation (Boström et al. 2010, Boyé et al. 2017, Kraan et al. 2020b). Moreover, seagrass meadows might limit dispersion of infaunal species by reducing hydrodynamics and bedload transport, and their physical structure restricts bulldozing species, which may result in a lower species turnover within the seagrass habitat (Boström et al. 2010, Henseler et al. 2019. Seagrass meadows provide substrate for both infauna and epifaunal invertebrates (Boyé et al. 2017, Henseler et al. 2019). These two components of benthic fauna can display distinct patterns of spatial variation, with epifauna in general showing lower β-diversity and high dominance of grazers than infauna, which suggest strong habitat filtering (Boyé et al. 2017).
Natural and human-driven disturbances are constantly changing the seafloor, shaping communities over multiple spatial extents and frequencies (Gladstone-Gallagher et al. 2019, Solan et al. 2020. For instance, stingrays can forage and disturb small patches of sediment at a high frequency, whereas trawling can disturb much larger extents of the seafloor at a different frequency. Biotic engineering operates at multiple spatiotemporal scales, but is frequently overlooked in metacommunity dynamics . We advocate that the interdependence between space, environment and species interaction must be recognized, even in ecosystems where competition is not the main regulator of biodiversity, such as the marine soft-bottom ecosystems, where food and space are not limiting factors. Ecosystem engineering and facilitation cascades can play a key role in shaping β-diversity, recognizing their importance for biodiversity and ecosystem function is thus essential for understanding the resilience of communities to multiple disturbances.