Deposit-feeding worms control subsurface ecosystem functioning in intertidal sediment with strong physical forcing

Abstract Intertidal sands are global hotspots of terrestrial and marine carbon cycling with strong hydrodynamic forcing by waves and tides and high macrofaunal activity. Yet, the relative importance of hydrodynamics and macrofauna in controlling these ecosystems remains unclear. Here, we compare geochemical gradients and bacterial, archaeal, and eukaryotic gene sequences in intertidal sands dominated by subsurface deposit-feeding worms (Abarenicola pacifica) to adjacent worm-free areas. We show that hydrodynamic forcing controls organismal assemblages in surface sediments, while in deeper layers selective feeding by worms on fine, algae-rich particles strongly decreases the abundance and richness of all three domains. In these deeper layers, bacterial and eukaryotic network connectivity decreases, while percentages of clades involved in degradation of refractory organic matter, oxidative nitrogen, and sulfur cycling increase. Our findings reveal macrofaunal activity as the key driver of biological community structure and functioning, that in turn influence carbon cycling in intertidal sands below the mainly physically controlled surface layer.


Introduction
Burrowing macrofauna have impacted seafloor habitats and altered the marine oxygen and sulfur cycles since their first appearance during the Proterozoic-Cambrian transition (1). Today, these organisms are ubiquitous in seafloor sediments except under conditions of oxygen-depleted bottom water, toxic hydrogen sulfide concentrations, or high temperatures (2). These macrofauna influence sediment matrices through a process called "bioturba-Ventilation and reworking influence microbially mediated processes by modulating the redox state and gradients of electron acceptors and donors within sediment. By introducing the highenergy electron acceptors oxygen and nitrate from overlying water and creating redox oscillations, ventilation stimulates the microbial oxidation of organic carbon (OC) (4), coupled nitrificationdenitrification reactions (6), and oxidative removal of potentially inhibitory microbial end products, such as reduced metals, sulfide, and ammonium (7). In contrast, reworking can fuel microbial activities in deeper layers by introducing freshly deposited organic matter (8), nutrient-rich secretes, e.g. mucopolysaccharides (9), and metal oxides (10), and represents an important control of OC mineralization and burial in marine sediment (1-4, 11, 12).
Comparatively less is known about the interactions among sediment macro-, meio-, and microbiota, even though these interactions are crucial to understanding the flow of carbon through benthic food webs (12). Macrofauna alter the abundance, diversity, and community structure of micro-and meiobiota locally in burrow walls, feeding pockets, and feces (13,14). More recent studies on subtidal sediments show that the influence of macrofaunal activity is not restricted to burrows but controls microbial community structure throughout the entire bioturbated layer (15,16). This results in highly similar microbial communities across bioturbated coastal and continental shelf sediments that differ in macrofaunal species compositions, sedimentation rate, sediment geochemistry, and lithology (16). Hereby, the dominant lineages in bioturbated surface sediment differ from those that dominate nonbioturbated surface or subsurface marine sediment.
Compared to subtidal habitats, intertidal sediments experience stronger daily fluctuations in tidal and wave action, temperature, and water content (17). Despite these frequent perturbations, intertidal sediments have among the highest benthic microalgal productivities (10 3 to 10 4 mmol C m 2 y -1 ) and OC decomposition rates (10 3 to 10 5 mmol C m 2 y −1 ) across seafloor ecosystems (18)(19)(20), and can host significant populations of macrofauna (21). Both hydrodynamic forcing and bioturbation impact O 2 dynamics, OC remineralization, and nutrient fluxes in intertidal sediments (20,22). In addition, bioturbation changes food-web structures by promoting microbial grazing (13), symbioses (23), and meiofaunamicroorganism interactions (24). Despite the known impacts of hydrodynamic forcing and bioturbation, it remains largely unknown how both forces compare, or interact with each other, in intertidal ecosystems.
Here, we explore the impact of physical forcing and bioturbation on intertidal carbon cycling and community functioning based on sandy sediment of False Bay, Washington, United States. Study sites were dominated by lugworms, a group of polychaetes with worldwide distributions in coastal sediments (25) that account for up to 30% of benthic biomass in intertidal sediments (21). Past research on Abarenicola pacifica, the dominant lugworm species in False Bay, led to the concept of "microbial gardening," which originally described the stimulation of microbial growth by A. pacifica in subsurface sediments (13).
We investigate the importance of hydrodynamic forcing and macrofaunal activity based on adjacent plots with and without lugworms that share similar hydrodynamic and geochemical conditions. To this end, we integrate analyses of sedimentary features (grain size and biogenic structures), porewater geochemical gradients (electron acceptors, metabolic end products, and redox potential), OC (content, isotopic compositions, carbon-to-nitrogen ratios (C:N), and chemical reactivity), ribosomal RNA gene sequences of Bacteria, Archaea, and Eukarya (abundance, phylogenetic richness, and community structure), and sulfur (S)-and ni-trogen (N)-cycling catabolic genes with modeled rates of hydrodynamics and macrofauna-driven porewater and particle mixing. Based on our findings, we propose that lugworms, and possibly other burrowing macrofauna, drive intertidal sand ecosystems below the physically controlled surface layer.

Impact of A. pacifica on sediment particle distributions
Previous studies indicate that A. pacifica forms "J-shaped" burrows in sandy sediments (Fig. 1A), at the end of which it selectively ingests fine particles, and thereby causes a constant downward movement of sediment (13,25). The digested fine particles are expelled back to the sediment surface and homogenized with surrounding sediment during tidal inundation. Matching this feeding behavior, we observed on average coarser grain sizes in A. pacifica-inhabited (bioturbated) sediments compared to A. pacificafree (nonbioturbated) sediments (P < 0.05, pairwise t test), with maximum grain sizes at ∼15 to 25 cm (Fig. 1A). A further consequence of particle preference and feeding-related downward movement of sediment by A. pacifica was that sediments below burrows (∼20 to 25 cm) were highly enriched in woody debris, and to a lesser degree macroalgae (Ulvophyceae) detritus and bivalve shells. These biogenic macrostructures were largely absent in nonbioturbated sediments. Based on observed living depths of A. pacifica and distributions of macrostructures, measured sediment grain size distributions, and modeled coefficients of physical and biological mixing (details see later), we divide sediments with worms into (1) physically and biologically impacted layers (PBL; 0 to 12.5 cm) that overlie the main living depth of A. pacifica; (2) biologically impacted layers (BL; 12.5 to 25 cm), where lugworms perform most of their feeding and reworking; horizontal burrows, which correspond to the main living depths of A. pacifica, were typically located in the upper half of the BL; and (3) undisturbed layers (UL; > 25 cm). Nonbioturbated sediments were divided into physically impacted layers (PL; 0 to 12.5 cm), where significant sediment mixing and porewater exchanges due to tide-related water movements were detectable, and UL below (> 12.5 cm).

Impact of lugworm bioturbation on porewater geochemical gradients
Microsensor profiles reveal shallower oxygen (O 2 ) penetration, lower hydrogen sulfide (H 2 S) accumulation, and more variable but higher pH and redox potentials in bioturbated compared to nonbioturbated sediments (Fig. 1B). While dissolved O 2 concentration peaks near the sediment-water interface likely reflect microphytobenthic photosynthesis, O 2 penetrations are shallow (< 0.5 cm) in both types of sediments. H 2 S concentrations in the top 4 cm remain low (< 5 μM) in bioturbated sediments but increase to 20 to 100 μM in nonbioturbated sediments. pH values decrease with depth and are slightly higher at > 1 cm in bioturbated (7.6 ± 0.2) than in nonbioturbated sediments (7.3 ± 0.1). Redox potentials decrease from 400 to 600 mV in overlying water to 333 ± 84 mV and 86 ± 56 mV at > 1 cm in bioturbated and nonbioturbated sediments, respectively, indicating more oxidizing conditions in bioturbated sediments.
Bioturbated sediments generally have higher concentrations of electron acceptors (nitrate (NO 3 -), sulfate (SO 4 2-)) and end products of microbial metal cycling (Fe 2+ , Mn 2+ ), but lower accumulations of other reduced end products (ammonium (NH 4 + ), hydrogen sulfide (H 2 S)) compared to nonbioturbated sediments  Mn 2+ concentrations also remain in the low micromolar range, but increase in deeper layers, especially in the BL. DIC concentrations fluctuate between ∼1 to 3 mM in both types of sediments.
Porewater exchange intensities (rightmost panel in Fig. 1C) are estimated from model simulations matching the DIC concentration profiles in manipulation experiments with and without worms (see "Methods"), and indicate physically induced porewater exchanges (α P ) in the top ∼10 cm of sediment. Biologically induced porewater mixing (α B ) indicates active pumping of seawater into worm burrow that causes bioirrigation to ∼10 to 20 cm in bioturbated sediments, but not in nonbioturbated sediments.

Impact of reworking on distributions of OC pools and sources
Depth profiles of OC quantity (total organic carbon; TOC), source (C:N, δ 13 C-TOC), and chemical reactivity (chlorophyll a (chl a), pheopigments, freshness index (ratio of chl a to the sum of chl a and pheopigments) indicate that A. pacifica modifies the distributions of OC fractions (Fig. 2). This is confirmed by copy numbers of ribulose-1,5-bisphosphate carboxylase genes (rbcL) of Ochrophyta (mainly diatoms based on 18S rRNA gene sequencing) and vascular plants. TOC and C:N in bioturbated sediments are constantly low in the PBL (0.13 ± 0.04% and 8 ± 2, respectively), increase to ∼0.6% and ∼30 at the bottom of the BL, and decrease to ∼0.1% to 0.2% and ∼10 to 20 below in the UL, respectively ( Fig. 2: Bulk OC). δ 13 C-TOC generally decreases from ∼ −19‰ at 1 cm to ∼ −27‰ at the bottom of the BL, with a slight local increase (∼ −21‰) at ∼16 cm. The maxima in TOC and C:N and minima in δ 13 C-TOC at the bottom of the BL correspond to high-sand layers below A. pacifica burrows. In comparison, TOC, C:N, and δ 13 C-TOC in nonbioturbated sediments are more variable and show no consistent trends related to depth.
Chl a values decrease with depth, with slightly slower depth attenuation and consequently higher chl a at 12.5 to 25 cm in bioturbated sediments. By contrast, pheopigments (chl a degradation products) are markedly lower in bioturbated sediments compared to nonbioturbated sediments (0.5 ± 0.7 vs. 2.0 ± 1.2 μM, P < 0.001). The freshness index decreases more gradually and is higher (P < 0.001) in the top 25 cm of bioturbated sediment compared to nonbioturbated sediments. rbcL abundances of Ochrophyta and vascular plants in bioturbated sediments are stable around 10 6 and 10 4 copies g -1 from 0 to 12.5 cm, respectively, and drop to ≤ 10 2 genes g -1 in the BL. rbcL abundances of Ochrophyta re-main low, whereas those of vascular plants increase again by 1 to 2 orders of magnitude in the bottom of the BL or UL. In nonbioturbated sediments, Ochrophyta and vascular-plant genes are slightly higher (∼10 6 to 10 7 and ∼10 4 to 10 5 genes g -1 , respectively) from 0 to 12.5 cm and decrease less with depth, remaining at ∼10 4 and ∼10 4 to 10 5 genes g -1 , respectively.
Coefficients of physically (D P ) and biologically (D B ) induced particle mixing (see "Methods" for details) suggest a strong decrease in D P from surface to subsurface sediments, and an increase in D B in subsurface layers with A. pacifica ( Fig. 2: Reworking). Matching the latter, feeding intensity (F B ; see "Methods") increases significantly in the BL.

Influence of lugworm bioturbation on gene abundances
Abundance profiles of bacterial and archaeal 16S rRNA and eukaryotic 18S rRNA genes, and of functional genes involved in microbial S-and N-cycling (dsrB: ß subunit of dissimilatory sulfite reductase; soxB: ß subunit of thiosulfohydrolase; narG: respiratory nitrate reductase 1 α chain; and amoA: α subunit of ammonia monooxygenase) show a clear impact of lugworms. This impact is greatest in the BL where lugworms live and feed (Fig. 3A). Within the PBL, bacterial gene abundances decrease from ∼10 9 to ∼10 8 genes g -1 , archaeal abundances increase from ∼10 6 to ∼10 7 genes g -1 , and eukaryotic abundances decrease from ∼10 8 to ∼10 6 genes g -1 . Similar trends hold in the PL, yet with, on average, ∼3-fold higher bacterial and archaeal and ∼3-fold lower eukaryotic gene abundances. In the BL, prokaryotic and eukaryotic gene abundances decline by ∼2 and ∼1 orders of magnitude, respectively, before increasing again in the UL. rRNA gene abundances decrease less over the interval from 12.5 to 25 cm in nonbioturbated sediments and show a more continual decrease to the deepest sam- gene copies (as log10 gene copies g -1 wet sediment, upper row), as well as Bacteria-to-Archaea (BAR), Archaea-to-Eukarya (AER), and Bacteria-to-Eukarya (BER) rRNA gene ratios (lower row). (B) Depth profiles of S-and N-cycling catabolic gene copies (sulfate reduction: dsrB; sulfide oxidation: soxB; nitrate reduction: narG; and ammonium oxidation: amoA; as log10 gene copies g -1 wet sediment, upper row), and percentages of these genes in %, calculated by dividing dsrB, soxB, narG, and amoA gene copies by corresponding total 16S rRNA gene copy numbers (lower row). All values represent averages from three plots that were sampled at different time points (error bars denote standard deviations). Significance level of Welch's t test: * * * P < 0.001, * * P < 0.01, * P < 0.05, 0.05 < P < 0.1. (Abbreviations: PBL = physically and biologically impacted layer, BL = biologically impacted layer (both "Bioturbated" sediments); PL = physically impacted layer, and UL = undisturbed layer (both "non-bioturbated" sediments.) ple (30 cm). Notably, the equally nonbioturbated deepest samples have similar gene abundances in both types of sediments. Gene abundance ratios between domains, calculated as a proxy for relative abundance changes (Fig. 3A), vary with depth in both bioturbated and nonbioturbated sediments. Bacteria:Archaea Ratios (BARs) decrease from ∼10 2 to ∼10 1 , with slightly higher values in the BL than at comparable depths in nonbioturbated sediments. Eukarya:Archaea Ratios (EAR) and Eukarya:Bacteria Ratios (EBR) decrease from ∼10 1 to ∼10 -2 and 10 -2 to ∼10 -3 , respectively, and are on average 5-and 8-fold higher in the PBL than in the PL, respectively. Notably, EARs and EBRs in bioturbated sediments have localized subsurface peaks (∼16 to 25 cm) due to increases in eukaryotic gene copy numbers in the lower part of the BL.
Genes involved in reductive (dsrB) and oxidative (soxB) S cycling are generally 1 to 3 orders of magnitude more abundant than genes involved in reductive (narG) or oxidative (amoA) N cycling (Fig. 3B). Depth trends are similar to those in 16S rRNA gene copy numbers. In bioturbated sediments, gene abundances are highest in the PBL, followed by a strong decline in the BL, and a slight increase in the UL. In nonbioturbated sediments, gene abundances decrease less with depth. In both surface sediments (0 to 12.5 cm) and deeper layers (12.5 to 25 cm), dsrB, soxB, and narG copy numbers are lower in bioturbated than in nonbioturbated sediments. By comparison, amoA gene copies are higher in surface layers of bioturbated sediments, and do not differ significantly between bioturbated and nonbioturbated sediments from 12.5 to 25 cm. Ratios of S-and N-cycling to total 16S rRNA gene copies suggest higher variability in percentages of S-and N-cycling microorganisms in bioturbated sediments (Fig. 3B). Overall, ratios are lower (dsrB), slightly higher (soxB), or clearly higher (narG and amoA) in bioturbated compared to nonbioturbated sediments.

Impact of lugworm bioturbation on benthic community structure
The impact of bioturbation on relative abundances of microbial and eukaryotic taxa is strongest in the BL, where worms mainly live and feed (Fig. 4A). In the following text, we present major groups at the class and phylum level with mention of their dominant subgroups in parentheses.
Microbial communities in surface layers of both bioturbated and nonbioturbated sediments are dominated by Flavobacteriia (Flavobacteriaceae), Acidimicrobiia (OM1_clade), Alphaproteobacteria (Rhodobacteriaceae), Gammaproteobacteria (Woeseiales and Halieaceae), and Woesearchaeota (unclassified). These groups, which we refer to as "surface lineages" from here on, gradually decrease with depth, whereas "subsurface lineages" Latescibacteria (unclassified), Chloroflexi (Anaerolineae and Dehalococcoidia), and archaeal Thermoplasmata (Marine Benthic Group D) show the opposite trend. Other groups, such as bacterial Deltaproteobacteria (Desulfobulbaceae and Desulfobacteracea) and Planctomycetes (Planctomycetaceae), or archaeal Bathy-, Thor-, and Lokiarchaeota show no systematic changes with depth. Sediments of the BL show distinct differences compared to nonbioturbated sediments. Relative to the same sediment depths in nonbioturbated sediments, the BL has lower percentages of most "surface lineages" and "subsurface lineages," but higher percentages of Gammaproteobacteria (Thiotrichaceae, Piscirickettsiaceae, and Ectothiorhodospiraceae) and several bathyarchaeotal subgroups (Group C3, MCG-8, and MCG-17). Eukaryotic communities in surface sediments are dominated by Metazoa and Ochrophyta, with higher percentages of Metazoa (Arthropoda and Annelida) and lower percentages of Ochrophyta (benthic diatoms Gedaniella sp.) in bioturbated compared to nonbioturbated sediments. In bioturbated sediments, Dinoflagellata (Dinophyceae), and Apicomplexa (Eugregarinorida, obligatory invertebrate parasites) increase in deeper layers of the PBL (7 to 12.5 cm), while Streptophyta (seagrass Zostera marina, which forms a meadow at the subtidal end of False Bay) and Chlorophyta (macroalgal Ulvophyceae, unicellular Chlorodendrophyceae) become abundant in the BL. Similar trends occur in nonbioturbated sediments, except that Ochrophyta and to a lesser extent Apicomplexa dominate over a larger interval (1 to 16 cm), whereas Strepto-and Chlorophyta only become dominant below 16 cm. Notably, relative abundances of Metazoa (meiofaunal worms Nematoda and Gnathostomulida) are significantly elevated in the BL while those of Ochrophyta are significantly lower compared to the same depths in nonbioturbated sediments (both P < 0.01).

Assembling mechanisms of benthic community
Richness analyses confirm the strong impact of lugworms on intertidal sediment communities (Fig. 4B). While bacterial and archaeal richness are comparable in the PBL and PL, they are clearly lower in the BL than in the same depths of nonbioturbated sediments. In contrast, eukaryotic richness is higher in the PBL than in the PL but similar in subsurface layers of both bioturbated and nonbioturbated sites.
Principal Coordinates Analysis (PCoA; weighted Unifrac distance that considers both phylogenetic distances between taxa and their read percentages) based on zero-radius operational taxonomic units (ZOTUs) confirms the major role of lugworms in structuring benthic communities (Fig. 4C). While Bacteria and Archaea share similar communities between the PBL and PL (both PERMANOVA P > 0.05), their communities in the BL differ significantly from those at the same depths in nonbioturbated sediments (P < 0.01). The same result is obtained with an unweighted Unifrac algorithm that only considers the phylogenetic distances between communities ( Figure S1). By contrast, eukaryotic communities only show significant differences between bioturbated and nonbioturbated sediments when an unweighted Unifrac distance algorithm is used, indicating that bioturbation more strongly affects low-abundance eukaryotic taxa ( Figure S1). These patterns of eukaryotic community assembly do not change when Metazoan 18S sequences are excluded ( Figure S2).
Canonical Analysis of Principal Coordinates (CAP) was used to identify potential drivers of benthic communities (Fig. 4C). The included variables explain 76% of bacterial, 71% of archaeal, and 60% of eukaryotic community variations. The first axis (CAP1) generally differentiates bacterial and archaeal communities of surface (0 to 12.5 cm) and subsurface (12.5 to 30 cm) sediments. Surface communities are correlated with physically induced porewater exchange (α P ) and particle mixing (D P ), indicators of reactive OC sources (chl a, freshness index, Ochrophyta, and vascular plant rbcL copy numbers), and concentrations of NO 3 -
Additional, interdomain correlation analyses on bioturbated sediments show that surface sedimentary clusters of Bacteria (B1 and B2) and Archaea (A1) are positively correlated, while the eukaryotic cluster E1 only correlates with bacterial cluster B2 (all distance correlation R ≥ 0.6, P < 0.01; Fig. 5B). The prokaryotic subsurface clusters B3, B4, and A2 are also positively correlated with each other, and generally negatively correlated with surface clusters. By comparison, nonbioturbated sediments only show correlations between bacterial and archaeal surface (C1 and C3) and subsurface clusters (C2 and C4, Figure S4).

Discussion
Despite being subjected to strong physical forcing by wave action, tides, and daily fluctuations in temperatures, intertidal sediments are rich in sediment macrofauna and represent global hotspots  . Leftmost panels show differences in porewater and particle mixing between nonbioturbated and bioturbated plots. In nonbioturbated sediments, physical forcing dominates, causing porewater and particle mixing to be mainly limited to surface sediments (PL). In bioturbated sediments, worm activity in subsurface sediments profoundly impacts porewater and particle mixing throughout. This results in strong differences in redox gradients, whereby lugworms increase the depth of oxidizing conditions (light gray background color) and consequent onset of reducing conditions (black color). The biggest impact of worms is in the BL, where worms spend most of their time in the horizontal section of their burrows. High rates of sediment ingestion, whereby worms selectively feed on labile OC-rich fine particles, lowers the pool size of labile OC and effectively drives down bacterial, archaeal, and eukaryotic populations. Simultaneously, worm feeding avoidance of (mainly terrestrial) macrostructures increases the burial and pool size of refractory OC. In sediments overlying burrows (PBL), lugworm activity strongly promotes intra-and interdomain bacterial and eukaryotic network connectivity, while the opposite is the case in the BL where worms feed. By contrast, lugworm activity only has a minor impact on intradomain archaeal networks, or archaeal networks with other domains. Finally, bioirrigation by lugworms alters the distribution of oxidative and reductive N-(and to a lesser extent S-) cycling, enhancing the relative importance of ammonium oxidation (nitrification) and nitrate reduction (denitrification) throughout bioturbated sediment. (Abbreviations: PBL = physically and biologically impacted layer, PL = physically impacted layer, BL = biologically impacted layer, and UL = undisturbed layer.) of terrestrial and marine OC cycling (12,(17)(18)(19)(20)(21)(22). While it remains unclear how physical forcing and macrofaunal populations compare in their impact on intertidal ecosystem functioning, a diminished importance of bioturbation has been proposed for sediments where hydrological processes dominate (5,26). Our study shows that, even in intertidal sands that experience strong physical forcing, benthic macrofauna-in this case the lugworm A. pacifica-can have a major impact on the transport and cycling of OC, and on the abundances, community structures, and networks of sediment biota.
While modeled rates of physical and biological processes suggest dominance of physical over biological processes in the top 5 to 10 cm of sediments (Figs 1 and 2), a strong impact of bioturbation is evident in deeper sediments that are sheltered from physical forcing (Figs 3 and 4). Here, lugworms significantly alter the porewater chemistry through ventilation. This introduces seawater-derived electron acceptors, e.g. O 2 and NO 3 -, while removing end products of anaerobic mineralization, e.g. NH 4 + and H 2 S. In addition, selective feeding on fine, algal organic matterrich particles in deeper layers induces a downward movement of sediment, including fresh algal organic matter, and strongly enriches coarse sediment grains and organic macrostructures below burrows. This selective feeding greatly lowers the abundance and richness of, presumably fine particle-associated, Bacteria and Archaea, and negatively affects all but a subset of microbial taxa involved in oxidative N-and S-cycling, macrofaunal symbioses, and degradation of organic macrostructures. Lugworm feeding also negatively impacts eukaryotic abundances and transforms organismal networks, mainly by enhancing network connectivity among "surface lineages" and reducing correlations among "subsurface lineages" (Fig. 5). This strong impact on organismal networks suggests that lugworms fundamentally alter the functioning of intertidal ecosystems. We synthesize the most important findings of our study as described above in a conceptual plot (Fig. 6) and provide detailed explanations in the following sections.

Physical processes mainly control surface sedimentary communities
Bacterial and archaeal communities in surface sediments (< 10 cm) are distinct from those in other sediment layers (Fig. 4), but do not change significantly in the presence of A. pacifica. Though physical and biological processes both enhance porewater exchange and sediment mixing, our data indicate that physical ventilation (α P ) and particle mixing (D P ) dominate in surface sediments (Figs 1 and 2). These results are in line with those of a tracerbased study conducted in a nearby but lower part of the False Bay intertidal area, that reported physical porewater advection throughout the upper 10 cm of sediment (27). α P and D P are the variables that explain the highest percentages of bacterial (∼30% and ∼28%) and archaeal (∼33% and ∼30%) community variation in surface sediments (Fig. 4C), whereas α B and D B appear more important in the BL, where both reach their peak values. Physical factors appear to control microbial communities in surface sediments through their strong impacts on sediment geochemistry. α P and D P are significantly positively correlated with NO 3 concentrations, chl a content, OC freshness, and rbcL copy numbers (especially Ochrophyta), and significantly negatively correlated with C:N, δ 13 C-TOC, H 2 S, and NH 4 + ( Figure S5A). All of these variables correlate significantly with the bacterial and archaeal community structure in surface sediments (CAP, Fig. 4C), indicating that physical processes drive microbial communities in surface sediments through increased input of seawater-derived electron acceptors, removal of potentially inhibitory or toxic (anaerobic) mineralization end products, and transport of labile, e.g. benthic Ochrophyta (diatom)-derived, OC from the sediment surface to deeper layers. While the geochemical agents are similar to those previously proposed to drive microbial community structure in continental shelf surface sediments (16), it is physical forcing rather than bioturbation that controls these variables in intertidal surface sands of False Bay. Matching the fluctuating physical and redox conditions in sandy, advection-controlled surface sediments (28,29), phylogenetic analyses indicate high metabolic flexibility and/or physiological resilience among "surface lineages." Dominant "surface lineages," such as Flavobacteriaceae (Flavobacteriia) and Rhodobacteraceae (Alphaproteobacteria), are widespread across diverse habitat types (seawater, sediment, soil, and biofilms) and include (facultatively) aerobic and anaerobic chemoorgano-and photo(organo)trophic members (30,31). Other dominant "surface lineages" can switch between aerobic respiration and fermentation during redox oscillations (Woeseiales (Gammaprotebacteria) (29)), scavenge oxygen for anaerobic growth (Desulfobulbaceae (Deltaproteobacteria) (32)), and occur widely in oxic and anoxic surface sediments (Woesearchaeota (33)).
Despite the likely dominance of physical over biological processes, bioturbated surface sediments have a higher eukaryotic abundance and richness, and different eukaryotic community structure than nonbioturbated sediments (Figs 3 and 4). Of all sequenced eukaryotic ZOTUs, 44% are shared between the PBL and PL, while the rest are unique to the PBL (35%) or PL (21%). Dominant unique taxa in bioturbated sediments include many meiofauna ( Figure S6), such as parasitic worms (Chromadorea within Nematoda (34)) and alveolates (Actinocephalidae within Gregarinomorphea (35)), and worms that feed on microbial populations in suboxic, sulfur-rich environments, often near polychaete burrows (Haplognathia within Gnathostomulida (36)). In addition, relative abundances of surface deposit-feeding polychaetes (Boccardiella spp. (37)) were clearly elevated in the PBL and explain the high percentage of Metazoa at 5 cm sediment depth in bioturbated sediment ( Fig. 4A; Figure S7).
Enhanced sediment permeability leading to increased input of oxygenated seawater and removal of toxic metabolites, e.g. H 2 S, in bioturbated surface sediments may promote the thriving of these meiofauna and other eukaryotic organisms (22,36), as could higher energy supply due to enhanced input of labile OC. The latter explanation is supported by higher freshness values and lower pheopigment contents in the PBL compared to the PL (Fig. 2B). Higher input and turnover of labile OC likely results from subsurface "conveyor feeding" by A. pacifica, whereby the microalgae-rich surface sedimentary layer is perpetually "subducted" into deeper layers via worm feeding funnels (Fig. 6). The subsequent return of ingested sediment to the surface as nutrient-rich feces, and the upward advective flow of nutrient-rich porewater from depth, could foster rapid regrowth of microalgae and lead to enhanced biomass turnover in worm-inhabited sediments (38).
Yet, the reason why higher input and turnover of labile OC, and thus bioavailable energy, in the PBL would increase eukaryotic biomass but not that of Bacteria and Archaea (Fig. 3A), is unclear. Past studies have shown that worm feces are microbially recolonized from surrounding sediments within hours (14) and support high growth rates that lead to recovery of microbial communities within 24 hours (39). We speculate that motile, grazing meiofauna (e.g. Gnathostomulida) and protists (e.g. Ciliphora) actively move through the interstitial space of sandy sediments (36,40) and maintain their vertical position even when exposed to a downward movement of surrounding sediment particles. This may enable these meiofauna to exploit the enhanced supply of microbial biomass in the PBL while at the same time avoiding predation by A. pacifica in the underlying BL. The same would not apply to most Bacteria or Archaea, which are particle-attached and grazed upon by motile meiofauna during transport through the PBL and subsequently by A. pacifica in the BL. Exceptions might include free-living prokaryotes, e.g. ammonium-oxidizing archaea (41), which occur at an order of magnitude higher abundances in the PBL than the PL ( Figure S8). These archaea could evade ingestion by particle-feeding A. pacifica, and benefit from the elevated supply of O 2 from bioirrigation in sediments above A. pacifica burrows, as indicated by the significant positive correlation between amoA copies and α B ( Figure S5B). In addition, grazing by unicellular eukaryotes such as ciliates and flagellates may lower the microbial biomass in this layer.

Lugworm activity drives subsurface communities
Worm feeding (F B ) and sediment mixing (D B ) have a strong impact on the community structure of Bacteria, Archaea, and Eukarya in subsurface sediments (Fig. 4C), and both negatively affect prokaryotic gene copy numbers, and prokaryotic and eukaryotic richness (Figs 3 and 4B; Figure S5B).
While D B values are positively correlated with freshness index, in line with increased input and cycling of labile OC in bioturbated sediments, F B values have positive correlations with C:N ( Fig. 2; Figure S5A). Pheopigment contents and rbcL copy numbers are furthermore more strongly (negatively) impacted by F B than by D B (Fig. 2; Figure S5A). These correlations are in line with the observation that A. pacifica selectively feeds on microalgae-and prokaryote-rich fine particles (13). This leads to the gradual enrichment of noningestible coarse sands (Fig. 1A) and macrostructures, such as woody debris, seaweed detritus, and bivalve shells below lugworm burrows ("biogenic bedding"), as observed for other ecosystems inhabited by deposit-feeding fauna (4,42). This selective feeding has a strong impact on bulk OC compositions, as evidenced by C:N and δ 13 C-TOC. These range from fresh microalgae-dominated in surface sediments (C:N: 7 to 10; δ 13 C-TOC: −23 to −18‰) to typical of terrestrial C3 vascular plant origin (C:N: ∼30; δ 13 C-TOC: ∼−27‰; (43)) below lugworm burrows (Fig. 2). The low rbcL copy numbers of vascular plants (Fig. 2B) are consistent with this vascular plant matter being dominated by wood and/or highly degraded detritus. Previous studies have suggested that lugworms, through downward mixing of fresh detritus, bioirrigation, and release of nutrient-rich secretes, stimulate microbial and meiofaunal growth in subsurface burrows ("microbial gardening" (13,44)). By contrast, our study suggests that these potentially positive effects are largely offset by selective feedinginduced negative impacts of A. pacifica on microbial and meiofaunal biomass.
In addition, subsurface bioirrigation and commensal interactions could play a role in structuring subsurface communities. Percentages of several gammaproteobacterial families (Thiotrichaceae, Piscirickettsiaceae, and unclassified Ca. Thiobios) increase in the BL (P < 0.05). Members of Thiotrichaceae, Piscirickettsiaceae, and Ca. Thiobios include sulfide oxidizers that use O 2 (49,50) or NO 3 as electron acceptors (51). These groups may benefit from input of O 2 and NO 3 by subsurface bioirrigation, which would match the higher relative percentages of narG and soxB genes in the BL (Fig. 3B). Intriguingly, members of these families have also been found in associations with marine fauna (Thiotrichaceae: echinoids (52); Piscirickettsiaceae: fish (49); and Ca. Thiobios: ciliates (50)). In contrast, relative abundances of "subsurface lineages" such as bacterial Chloroflexi (mainly Anaerolineaceae (53) and Dehaloccocoidia (54)) and archaeal Euryarchaeota (mainly Marine Benthic Group D within Thermoplasmata (55)), that include anaerobic fermenters and acetogens, decrease in the BL, matching the reported negative impact of bioirrigation on these anaerobic lineages (16). The metazoan community is dominated by free-living Nematoda (Trefusidae) and Gnathostomulida (again Haplognathia), which both increase in diversity and relative abundances in the BL (Figure S7). Members of both groups are frequently enriched in suboxic sediments but differ in diets (Trefusidae: selective deposit feeders (56); Haplognathia: microbial grazers (36)). Both groups may live by scavenging on ventilation-induced O 2 , food scraps, and metabolites released by macrofauna, and evade worm predation by means of their free-living lifestyles and/or by inhabiting coarse particles that are not ingested by A. pacifica (36,56,57).

Lugworms modify sediment community networks
Our research shows that lugworms strongly alter community networks by increasing positive correlations among surface taxa (B2, A1, and E1) while reducing those among subsurface taxa (B3, B4, and A2; Fig. 5). In addition, negative correlations between surface and subsurface bacterial lineages (C1 vs. C2) are greatly reduced in the presence of lugworms.
Worm bioturbation supports unique bacterial and eukaryotic clusters (B2 and E1) among "surface lineages," which are, moreover, significantly correlated with each other (Fig. 5). The bacterial cluster B2 contains putative aerobic degraders of algal polysaccharides (Planctomycetacia, mainly Pirellula and Rhodopirellula (58)) and cellulose (Bacteroidetes VC2.1_Bac22 (59)), and metazoan symbionts (PAUC43f_marine_benthic_group (60)). The eukaryotic cluster E1 includes microalgae (Ochrophyta and Prasino-Clade-V), plant pathogens (Endomyxa-Phytomyxea (61)), bacterivorous protists (Filosa Granofilosea (62)), and ciliates that are known to host denitrifying endosymbionts (Plagiopylea (63)). While the metabolisms of these groups are typical of intertidal surface sediments, they are diverse and do not give clear indications that members of B2 and E1 enter biological consortia or codependencies in the presence of bioturbation. Instead these correlations could also be driven by factors that co-vary in lugworm-inhabited sediments, e.g. higher vertical sediment transport, increased input of labile OC from above and of O 2 from below, higher meiofaunal, and metazoan population size. These same variables, especially the higher rates of vertical sediment transport and bioirrigation, may also reduce the vertical zonation of surface and subsurface bacterial lineages and explain the absence of negative correlations between surface and subsurface bacterial lineages in the presence of lugworms.
Interestingly, presence of lugworms reduces correlations in subsurface sediments compared to A. pacifica-free sediments. This is most apparent in the bacterial cluster C2, which is broken down into two separate smaller clusters (B3 and B4; Fig. 5A). Major lineages of B3 include putative sulfate reducers (Zixibacteria (64)), widespread but poorly known members of the phylum Bacteroidetes (BD2-2), most members of which are primary fermenters of carbohydrates and proteins (65), and members of Latescibacteria, which genomic analyses suggest are primary fermenters of algal polysaccharides and glycoproteins (66). B4 also comprises groups that are linked to the anaerobic degradation of carbohydrates and proteins (Aminicenantes (67); Anaerolineae, mainly Anaerolineaceae (53); Dehalococcoidia (54); and Deferribacteres, mainly Caldithrix (68)). Yet, while B3 member are mostly found in marine sediments and seawater, B4 members are also widespread in activated sludge, terrestrial aquifers, and terrestrial deep subsurface environments (see references above), where the content of algal detritus is low. It is, thus possible that differences in primary OC sources drive the observed patterns. The near absence of B3 from the terrestrial OC-dominated deep part of the BL matches the fact that most of its members are linked to algal polysaccharide degradation ( Figure S9). By the same token, cluster B4 occurs at high relative abundance in the terrestrial OC-dominated deep layer, and could include key degraders of vascular plant detritus. Thus, the separation of B3 and B4 in bioturbated sediment may be the outcome of feeding-related partitioning of OC pools by A. pacifica.

Conclusions
Even though hydrodynamic forcing and macrofaunal bioturbation are the main sources of disturbance in coastal sediment (12,(15)(16)(17)(18)(20)(21)(22), the relative impact of both processes on biological communities was previously not known. We show that intertidal sand ecosystems dominated by lugworms can be divided into distinct hydrodynamically and biologically controlled layers. Hydrodynamic forcing, by regulating OC and electron acceptor inputs, is the main driver of biological communities in surface sediments. In contrast, communities in deeper layers, that are sheltered from hydrodynamic forcing, are largely controlled by lugworm activities. Though lugworms clearly modify OC and electron acceptor inputs in these deeper sediments, their most important influence on organismal communities is through selective feeding on fine particles. This selective feeding effectively depletes labile microalgal organic matter, drastically lowers abundances of Bacteria, Archaea, and Eukarya, strongly weakens biological networks, and promotes the burial of refractory OC. The clear division into hydrodynamically and bioturbation-controlled vertical zones provides a basis for integrating physical and biological factors into ecological and biogeochemical models of globally distributed intertidal sand ecosystems.

Study area
The study area is located in the northwest part of the False Bay Biological Preserve on San Juan Island, United States ( Figure S10). False Bay is an intertidal embayment (∼1 km 2 ), that is covered by coarse sand at the mouth and finer silty and muddy sand along the margins, and is subject to mixed semidiurnal tides (69,70). A total of two small streams enter the bay from north and west, respectively, but both were nearly dried up during the sampling period. Abarenicola pacifica dominates several upper intertidal areas of False Bay, one of which, and an adjacent nonbioturbated area were the focus of this study ( Figure S10). Throughout the study period in July/August 2017, sediment porewater in the study areas had chloride (Cl -) concentrations (470 to 510 mM) and salinities (∼30 to 33 psu; Figure S11) in the range of seawater from the Strait of Juan de Fuca (30 to 34 psu) (71), indicating minimal groundwater seepage. Temperatures during the day varied from 17 to 27 • C at the sediment surface and decreased to steady values of 15 to 17 • C below ∼15 cm ( Figure S12). Sediments were exposed for 6 to 8 hours during each low tide. For the first 2 to 3 hours of exposure, sediments remained water-saturated to the sediment surface.

Study design
To investigate the relative importance of hydrodynamic forcing and bioturbation in driving intertidal sediment geochemistry and biological communities, we performed field investigations on natural bioturbated and nonbioturbated areas. In addition, we used manipulation experiments within the bioturbated area to quantify rates of physical and biological porewater mixing. For an overview of the study design see Table 1 Investigations on natural sediments A "bioturbated" area dominated by A. pacifica (30 to 40 individuals m -2 ) was compared to an adjacent macrofauna-free "nonbioturbated" area ( Figure S10; Table 1). Both areas had similar hydrological conditions, lithologies, and OC sources. The "Bioturbated" area was sampled on August 5 (T1), 19 (T2), and 27 (T3), and September 7 (T4), while the nonbioturbated area was only sam-pled at T2, T3, and T4. Sampling was initiated around the time of tidal exposure (T1: ∼2 hours after exposure; T2 and T3: 1 hour before exposure with 10 to 20 cm of overlying water; and T4: ∼1 hour after exposure). Microsensor measurements and subsequent porewater and sediment sampling were completed within 1 hour. All data presented in this study except modeled porewater mixing rates (next paragraph) are based on these investigations on natural sediments.

Manipulation experiments
Manipulation experiments to model rates of physical porewater exchanges (α P ) and bioirrigation (α B ) were conducted (Table 1; model output shown in Fig. 1C). A 2 × 3 m 2 sediment plot in the bioturbated area was marked by wooden posts, which held up a 50-cm tall barrier made of thin garden netting (nylon, 5-mm mesh size) of which 20 cm were belowground to prevent recolonization by burrowing macrofauna. Macrofauna were removed from the entire plot by sieving the top 30 cm of sediment through a metal sieve with a 1-mm mesh size. Sieved sediments were homogenized and returned to the original plot. After 3 days of preincubation, the manipulated plot had recovered its original sediment compaction and surface features (e.g. ripples). The plot was then divided into two large subplots by inserting the same nylon net barrier to ∼20 cm depth through the middle. One half was kept macrofauna-free ("defaunated") to mimic natural nonbioturbated sediment, while A. pacifica were added at natural densities (35 individuals m -2 ) to the other half ("refaunated"). Over the next ∼4 weeks, porewater DIC and SO 4 2concentration profiles were measured and compared between defaunated and refaunated treatments to determine porewater exhanges by physical mixing and bioirrigation. Prior to T4, the external nylon net was damaged and plots invaded by small crabs. Thus, all modeling was exclusively based only on T1 and T2 data.

Microsensor profiling
In situ depth profiles of dissolved O 2 , H 2 S, pH, and redox potential were determined prior to sediment coring using a field micropro-

Sampling
Sediment cores (∼40 cm) were taken from the same location where microsensor profiling was performed using 90-mm diameter plexiglass core liners. Porewater was sampled using rhizons (0.15 μm pore size, Rhizosphere, Netherlands) that were inserted through predrilled holes in the liner. The initial dead spaces of 20-mL syringes, three-way stop cocks, and rhizons were flushed two successive times, each time with 1 to 1.5 mL of porewater, to remove atmospheric gasses and prevent oxidation of redoxsensitive ions with O 2 . After this flushing, ∼10 mL of sediment porewater were extracted from each of the 2-cm depth intervals and fixed immediately with acid or base for cation or anion analyses, respectively (for details see (72)). While in the field, all porewater samples were stored inside coolers with −20 • C ice packs. After porewater sampling, sediments were extruded and sliced into 1 to 2 cm thick vertical sections. From the center of each  I-III  I-III  I-III  I-III   Non-bioturbated  n. a.  I-III  I-III  I-III   Manipulation  experiments   Defaunated  I-II  I-II  n. a.  I-II   Refaunated  I-II  I-II  n. a.  I-II slice, ∼4 cm 3 of sediment were sampled for microbiological and solid-phase geochemical analyses using sterile, cut-off 5 mL syringes, and immediately frozen in coolers with dry ice (Note: due to an error, we subsequently analyzed T1, T2, and T4 samples for "Bioturbated" and T2, T3, and T4 samples for "Non-bioturbated" sediments. This does not impact our conclusions as geochemical and microbiological differences between "bioturbated" and "nonbioturbated" sediments clearly exceed temporal variations within each category over the study period.).  (72)). Standards consisted of solutions of MilliQ water containing analytical-grade ammonium chloride, zinc sulphide, and sodium salts of sulfate, chloride, nitrate, and nitrite, and an ICPmultielement standard solution (all Sigma-Aldrich, Switzerland).

Sediment geochemical measurements
Chl a and its degradation products (pheo-pigments) were extracted from ∼1.5 g wet sediments using acetone and quantified spectrophotometrically using an acidification protocol ((76), Cary 50, UV-Vis, Varian). Chl a freshness indices were calculated based on ratios of chl a to the sum of chl a and pheopigments (77). TOC, total nitrogen (TN), and their isotopic compositions (δ 13 C-TOC and δ 15 N-TN) were measured on decarbonized samples by elemental analyzer/isotope ratio mass spectrometry (16). Ratios of TOC to TN (C:N) served as indicators of OM source and/or quality (43). Grain size distributions were determined by dispersing and mixing ∼0.5 cm3 of sediments in sodium monophosphate prior to analysis by a laser diffraction particle size analyzer (LS13 320 with Aqueous Liquids Module and autosampler, Beckman Coulter, Indianapolis, USA) using a Polarization Intensity Differential Scattering detector (PIDS, for particles 0.04 to 0.4 μm) and a light scattering detector (according to ISO 13320-1, for particles 0.4 to 2000 μm) (78). Fractions of the volume distribution are grouped into the categories clay (0.04 to 2 μm), silt (2 to 63 μm), and sand (63 to 2000 μm).

DNA and RNA extraction
Total DNA and RNA was simultaneously extracted from ∼0.2 g of wet sediment following lysis protocol I of Lever et al. (79). This protocol involves chemical (phenol-chloroform-isoamylalcohol, lysis solution I) and mechanical cell lysis (bead-beating: 0.1-mm Zirconium beads), followed by 2 × washing with chloroform:isoamyl alcohol (24:1), precipitation with a mixture of linear polyacrylamide, sodium chloride, and isopropanol, and purification with the CleanAll DNA/RNA Clean-Up and Concentration Micro Kit (Norgen Biotek Corp., Canada).

Quantitative PCR (qPCR)
Bacterial and archaeal 16S rRNA genes, eukaryotic 18S rRNA genes, Ochrophyta rbcL genes, and functional genes encoding the ammonia monooxygenase subunit A of ammonia-oxidizing Archaea (AOA) and Bacteria (AOB) (amoA), respiratory nitrate reductase (narG) of denitrification, dissimilatory sulfite reductase (dsrB) of sulfate reduction, and thiosulfohydrolase (soxB) of sulfide oxidation were quantified by SYBR-Green based qPCR assays using a LightCycler 480 II (Roche Life Science, Switzerland). Published primer pairs and a newly designed primer mixture for soxB were applied (Tables S1 and S2). The new soxB primer mixture was designed based on a soxB database that was constructed in ARB (http://www.arb-home.de). This database consisted of 340 longread (≥ 900 bp) soxB sequences that were mostly from wholegenome sequencing studies. Primers were designed based on best matches with manually optimized sequence alignments. For further details, see the caption of Table S1.

Sequence and bioinformatic analyses
Amplicon libraries were prepared according to (16) using published 16S and 18S rRNA gene primers (Table S1) and by sequencing via the MiSeq platform (Illumina Inc., USA). Raw sequence reads were processed following a bioinformatic pipeline in Deng et al. (16).

Quantification of porewater and sediment particle mixing
Physical and biological porewater mixing was determined by reaction-transport modeling of DIC and SO 4 2concentration profiles in defaunated and refaunated sediments. Herein, we assumed that defaunated treatment was subject only to physical mixing, while refaunated treatment was subject to both physical mixing and lugworm bioirrigation. Both types of porewater exchanges were described as nonlocal exchanges with depth-dependent physical mixing frequencies α P (x) and bioirrigation coefficients α B (x) (for further details, see Supplementary Text). Solid-phase mixing was determined by adjusting diffusive mixing coefficients (D P and D B ) to match simulated with measured chl a and pheopigment profiles in natural sediments (see Supplementary Text for further details).

Feeding intensity
Feeding intensity (F B ), which we define as the selective feeding of lugworms on fine sediment particles, enriches coarse grain sizes within and below the living depths of lugworms (Fig. 1A) (4, 13). We provide a proxy for F B based on differences in sand fractions between bioturbated ( f sand )sand and nonbioturbated sediments ( f sand ) in natural sediments: (1)

Multivariate statistical analysis
All statistical analyses were performed in R (80). Microbial richness was calculated based on rarefied datasets of observed ZOTU numbers. Ordination analyses (PcoA and CAP) with weighted Unifrac distances were performed using the "phyloseq" package (81). Permutational multivariate analysis of variance (PER-MANOVA, permutations = 999) and statistical tests (two-tailed pairwise t test and Welch's t test) were performed using the "vegan" and "stats" packages, respectively (82). Community networks were constructed using pairwise Spearman correlations between taxa, which were based on all samples and time points recovered from the "bioturbated" and "non-bioturbated" areas. The networks were visualized using the "corrplot" package with the hierarchical clustering method (83). Interdomain correlation strength was quantified by calculating the coefficients of distance correlations (dcor R) using the "energy" package (84).