Macrofaunal control of microbial community structure in continental margin sediments.

Through a process called "bioturbation," burrowing macrofauna have altered the seafloor habitat and modified global carbon cycling since the Cambrian. However, the impact of macrofauna on the community structure of microorganisms is poorly understood. Here, we show that microbial communities across bioturbated, but geochemically and sedimentologically divergent, continental margin sites are highly similar but differ clearly from those in nonbioturbated surface and underlying subsurface sediments. Solid- and solute-phase geochemical analyses combined with modeled bioturbation activities reveal that dissolved O2 introduction by burrow ventilation is the major driver of archaeal community structure. By contrast, solid-phase reworking, which regulates the distribution of fresh, algal organic matter, is the main control of bacterial community structure. In nonbioturbated surface sediments and in subsurface sediments, bacterial and archaeal communities are more divergent between locations and appear mainly driven by site-specific differences in organic carbon sources.

Organic matter (OM) is supplied to the ocean through in situ photosynthesis and by external input from terrestrial sources. Although only a small fraction of OM escapes mineralization in the water column and reaches the sediment surface, marine sediments are the largest OM sink on Earth over geologic timescales (1). Whether sedimentary OM is preserved or mineralized is in part controlled by microorganisms, which through a network of hydrolysis, fermentation, and respiration (i.e., terminal oxidation) reactions convert OM to its inorganic constituents. Despite the important role of microorganisms in determining the fate of sedimentary OM, the factors that control microbial community structure in marine sediments are not well understood. A vertical zonation of respiration pathways based on Gibbs energies has been proposed (2). However, such a clear zonation is absent from many surface (3) and subsurface sediments (4), and respiring organisms often only account for a small percentage of anaerobic microbial communities (5). Instead, genomic and cultivation-based data suggest that most anaerobic sediment microorganisms gain energy from the fermentative or acetogenic breakdown of OM (6–10). Additionally, correlations between certain bacterial taxa and input of fresh phytoplankton detritus in surface sediment (11), widespread occurrence of certain Thaumarchaeota in oxic sediments (12), and dominance of phylogenetically distinct microbial taxa in energy-depleted subsurface sediments (13) indicate that OM quality, presence of O2, and overall energy availability are important drivers of microbial community structure.

A further factor that influences microbial activity and community structure in marine surface sediment is bioturbation. Globally, >95% of the seafloor underlies oxygenated bottom water (14), with most of this seafloor being inhabited by macrofauna (15). These macrofauna interfere with the vertical zonation of respiration reactions, change the distribution of OM, and influence rates of microbial OM remineralization in surface sediments (15, 16) by two different modes of behavior. During “reworking,” macrofauna displace particles, including sediment grains and organic detritus, through locomotion, burrow construction, ingestion, or defecation. Reworking is often divided into “biodiffusional mixing,” which describes random particle mixing over short distances, as well as “nonlocal mixing,” which describes unidirectional particle transport over distances of centimeters to decimeters (17). During “burrow ventilation,” macrofauna flush their burrows with overlying water, meanwhile causing “bioirrigation,” which describes the enhanced solute transport through bulk sediments around the burrow (17). Both processes vary with the types of macrofauna that are present. Suspension feeders pump large amounts of oxygenated seawater through their burrows and can thereby induce redox fluctuations in burrows and sediment porewater (18). By contrast, deposit feeders mostly perform particle reworking, and thereby mix freshly deposited OM and oxidized solid-phase compounds from the seafloor into deeper layers and vice versa (19).

Despite its well-studied impact on sediment geochemistry, little is known about the importance of bioturbation in structuring microbial communities. Burrow wall microbial abundance and diversity are often elevated and resemble those at the sediment surface (20, 21). In addition, sulfate-reducing, denitrifying, and nitrifying bacteria are often enriched in burrow walls (22, 23). Less is known about how bioturbation influences the composition of microorganisms throughout the bioturbated layer, even though the geochemical influence of bioturbation extends beyond burrows (18). Recently, Chen et al. (24) reported a dominance of Bacteria over Archaea throughout bioturbated sediment, and equal abundances of both domains in underlying subsurface sediments. The same study also noted a clear change in bacterial and archaeal taxa from bioturbated to nonbioturbated layers. Chen et al. (24) hypothesized that introduction of electron acceptors (e.g., O2, NO3−) by burrow ventilation and labile OM (e.g., phytodetritus) by particle reworking were responsible for the observed patterns in microbial communities.

Here, we test this hypothesis by analyzing depth profiles of microbial communities at four continental margin sites across the Baltic Sea–North Sea transition. Sites differ in water depth, sedimentation rate, metal and OM content, microbial activity, and macrofaunal community structure, and range from nonbioturbated to ventilation and/or reworking-dominated. At each site, we quantify ventilation and reworking intensities by performing reaction-transport modeling on porewater solutes and particulate tracers, respectively, and determine redox state, total organic carbon content (TOC), distributions of fresh OM, and OM sources based on C-isotopic values (δ13C-TOC) and eukaryotic 18S rRNA gene sequences. We then analyze trends in microbial community structure in relation to modeled macrofaunal activities and measured environmental variables. We show that burrow ventilation mainly drives archaeal, whereas particle reworking predominantly drives bacterial community structure in surface sediment. We moreover observe that microbial communities in nonbioturbated surface sediments more closely resemble those in underlying subsurface layers than those from similar depths in bioturbated sediment.

Results

Sedimentological Characteristics.

Four sites (AU1, AU2, AU3, and AU4) were sampled along a water-depth gradient in the Skagerrak–Kattegat–Lillebælt region between the North Sea and Baltic Sea (586, 319, 43, and 37 m, respectively; ). AU1 is located in the Norwegian Trench, has a high solid-phase manganese content (), and receives sediments with low-reactivity OM through horizontal transport by currents (25). AU2 is located on the southern slope of the Norwegian Trench. Both AU1 and AU2 consist mainly of silty clay and have high rates of iron and manganese reduction in the top 0 to 10 cm (26). AU3 in the northern Kattegat is dominated by fine sands and silts, and also has high rates of iron cycling (26). AU4 in the Lillebælt region of the Baltic Sea is sulfidic with the exception of a 1-mm-thick oxidized surface layer. Like AU1 and AU2, AU4 is dominated by silty clay, but the sediment has a more muddy consistency. While porosities decrease and densities increase in the top 50 cm of all sites (), sediment grain sizes do not change notably with depth (for detailed core descriptions, see ). Sedimentation rates based on 137Cs profiles increase with decreasing water depth and distance to shore, and are 0.14, 0.27, 0.30, and 0.33 cm⋅y−1 at AU1-4, respectively (). Macrofaunal biomass increases from AU1 to AU3, whereas macrofauna are absent from AU4 due to seasonal bottom water hypoxia (26) ().

Biogeochemical Profiles.

Vertical porewater profiles of oxygen (O2), sulfate (SO42−), methane (CH4), and dissolved inorganic carbon (DIC) suggest that microbial sulfate reduction, methanogenesis, and total respiration (DIC production) increase from AU1 to AU4 (Fig. 1). Accordingly, DIC profiles steepen, while both O2 penetration depth (AU1, 15 to 18 mm; AU2, 8 mm; AU3, 4 mm; AU4, <1 mm; from ref. 26) and the depth of the sulfate–methane transition (SMT), below which SO42− is depleted and CH4 accumulates, become shallower. Despite this general trend in microbial activity and despite ongoing sulfate reduction, SO42− and DIC concentrations are nearly constant in the top 60, 10, and 25 cm of AU1, AU2, and AU3, respectively.

Fig. 1.

Geochemical characterizations and indicators of bioturbation activity at AU1 to AU4. (A) Depth profiles of SO42−, CH4, DIC, TOC, δ13C-TOC, and C:N ratios across the entire cored intervals. Ventilation-dominated layers (VLs), reworking-dominated layers (RLs), and sulfate–methane transitions (SMTs) are indicated by horizontal dashed, solid, and dotted lines, respectively. (B) Bioturbation activities indicated by profiles of chl a, modeled rates of biodiffusional mixing (DB), nonlocal mixing (Na), and bioirrigation (α), as well as % Fe(III) in the top 0 to 40 cm. (C) Composition of macrofaunal communities, based on counts of individuals.

TOC trends vary between stations and are not correlated with water, sediment depth, or microbial activity (Fig. 1). The highest TOC is at AU4, where values increase with depth from 4.7 to ∼6% in the upper 50 cm. At AU3, which has the lowest TOC, values decrease from 0.8% at the seafloor to 0.5% at 50 cm—showing a local peak at ∼30 cm—and gradually increase back to 0.8% below. TOC decreases gradually from the seafloor to the core bottom at AU2 and AU1, from 2.1 to 1.7% at AU2 and from 1.9 to 1.3% at AU1.

δ13C-TOC and carbon-to-nitrogen (C:N) ratios at AU1 to AU3 fall within typical ranges of marine phytoplankton (27). δ13C-TOC is relatively constant with depth at AU1 (−22.8 ± 0.2‰) and AU2 (−22.5 ± 0.2‰), while at AU3 values decrease from −23.0‰ to −24.1‰ in the upper 50 cm and stabilize below. AU4 has δ13C-TOC similar to the other stations from 0 to 20 cm (−22.7 ± 0.3‰), but increases to approximately −20‰ at 50 cm and below suggest a significant additional TOC source, e.g., seagrass (28). C:N values vary minimally at AU1, AU2, and AU4 (mostly 8 to 9). At AU3, C:N ratios increase from ∼8 to 10 in the upper 40 cm, and further to ∼12 in the deepest sample.

Macrofaunal Activity and Community Composition.

Chlorophyll a (chl a) depth profiles reflect the input of fresh phytodetritus by sedimentation and translocation of this phytodetritus to greater depths by reworking (29) (Fig. 1). The contribution of chl a to TOC is negatively correlated with water depth (), indicating a decrease in labile OM deposition with increasing water depth. At AU1, chl a decreases threefold from 3 cm to underlying layers. At AU2 and AU3, chl a shows bimodal distributions, with surface peaks at ∼7 and ∼3 cm and subsurface peaks at ∼27 and ∼30 cm, respectively. The steepest chl a decrease is at the nonbioturbated AU4, where chl a values drop over 10-fold in the upper 5 cm and show only minor increases from 10 to 20 cm. Ratios of chl a/(chl a + phaeopigments), which are another proxy for fresh phytodetritus (30), confirm the trends observed in chl a ().

Macrofaunal reworking was quantified by modeling solid-phase biodiffusional mixing (DB) and nonlocal mixing (Na), whereas published solute-phase bioirrigation coefficients (α) on the same cores were used as proxies for ventilation rates (26). DB, Na, and α generally decrease with water depth, matching the decreasing macrofaunal biomass from AU3 to AU1 (Fig. 1 and ). AU1 is dominated by the stationary tube-dwelling polychaete Spiochaetopterus typicus, which feeds on suspended particles and surface deposits, and the surficial burrow-building, mainly ventilating chemosymbiotic bivalve Thyasira equalis (31) (Fig. 1C). DB and α values are relatively high in the upper 5 cm but decrease below detection below 10 cm. The macrofaunal community at AU2 is dominated by Thyasira equalis, the subsurface deposit-feeding bivalve Yoldiella lucida, surface deposit-feeding and ventilating Abra nitida, in addition to several deep-dwelling, deposit-feeding worms (Melinna cristata, Galathowenia oculata, Maldanidae sp., and Neoleanira tetragona). This functionally diverse assemblage matches the high DB and α values from 0 to ∼13 cm. However, the continuously high excess 210Pb and Na values to 40 cm are difficult to explain with these fauna. One possibility is that the thaliassinidean shrimp, Calocaris macandreae, which is common in the area and can dig intensively to 60 cm or deeper (32), was missed due to limited depth penetration of our box corer. Matching the high DB and α values throughout the upper 10 cm, AU3 is dominated by dense populations of the brittle star Amphiura filiformis (>1,000 individuals [ind.]⋅m−2) and the bivalve Kurtiella bidentata. Both are biodiffusors that typically live in <10-cm depth and feed on fresh phytodetritus (33). In addition, specimen of the ghost shrimp Callianassa tyrrhena, which forms deep, gallery-like burrows, were found in deeper layers (16 ± 4 ind.⋅m−2), and explain the clear subsurface peaks of 210Pb and Na at ∼20 cm. Matching the absence of macrofauna, DB, Na, or α values indicate no significant bioturbation activity at AU4.

To facilitate the comparison of geochemical and microbial community data, we define the upper layers of sediment, where bioirrigation rates (α) are high, as ventilation-dominated layers (VLs hereafter). Deeper layers, where bioirrigation is largely reduced but significant reworking, in particular nonlocal mixing, is still present, are termed reworking-dominated layers (RLs hereafter). The inferred zonation of dominant bioturbation activities generally matches depth profiles of chl a and %Fe(III) [percentage of Fe(III) of reactive solid-phase iron (Fig. 1); data from ref. 26]. %Fe(III), an indicator of sedimentary redox conditions, decreases in a similar fashion to α values throughout the VLs and remains at background values below and throughout the nonbioturbated site AU4. However, to avoid definitional biases, all correlation analyses between environmental variables and microbial community data performed later in this study are independent of the proposed divisions into VLs and RLs.

Microbial Abundance and Richness.

16S rRNA gene abundances of Bacteria show similar depth profiles across sites and decrease by one to two orders of magnitude from core tops to bottoms (Fig. 2). Except for a decrease in the top 4 cm of AU1 and AU3, archaeal gene abundances are relatively stable in the top 200 cm of AU1, 40 cm of AU2, 35 cm of AU3, and 20 cm of AU4, and only decrease with depth below.

Fig. 2.

Depth profiles of microbial communities at AU1 to AU4. (A) 16S rRNA gene abundances of Bacteria and Archaea (red triangles, Bacteria; blue diamonds, Archaea); (B) Bacteria-to-Archaea ratios (BARs); (C) Chao1 richness of Bacteria and Archaea; (D and E) community compositions of Bacteria and Archaea. The VLs, RLs, and SMTs are indicated by horizontal dashed, solid, and dotted lines, respectively. The black arrows indicate local peaks in chl a or excess 210Pb from Fig. 1.

Bacteria-to-Archaea ratios (BARs) decrease with depth, showing the highest values in the VLs of AU1 to AU3 and at the seafloor of AU4 (Fig. 2). At AU1 to AU3, BARs increase in surface sediments and then show a clear drop in the RLs. Local peaks in VLs and RLs generally match local peaks in chl a or excess 210Pb. Below the RLs, BARs continue to decrease to ∼100 cm, below which values are relatively stable. AU4 shows a different trend: Except for the peak value at the seafloor, BARs are lowest above the SMT, and increase gradually below the SMT.

Chao 1 richness is higher in Bacteria than in Archaea except at the sediment surface of AU4 (Fig. 2). At AU1 to AU3, bacterial and archaeal richness increase throughout the VLs and reach their highest values in the RLs before decreasing below. At AU4, bacterial richness stays relatively constant, whereas archaeal richness decreases strongly in the top 30 cm.

Depth Profiles of Microbial Communities.

Microbial community structure changes consistently from bioturbated to nonbioturbated sediments, and between the VLs and RLs of bioturbated sites (Fig. 2 ). In the VLs of AU1 to AU3, Bacteria belonging to Bacteroidetes, Acidobacteria, Alphaproteobacteria, Deltaproteobacteria, Gammaproteobacteria, and Planctomycetes are dominant. In the RLs, percentages of Gammaproteobacteria and Bacteroidetes decrease, while fractions of Deltaproteobacteria, Planctomycetes, and Chloroflexi increase (Fig. 2). The nonbioturbated surface sediment of AU4 (0 to 20 cm) is also rich in Bacteroidetes, Gammaproteobacteria, Deltaproteobacteria, and Planctomycetes, but here these groups coexist with a significant percentage of Chloroflexi. The nonbioturbated layers of AU1 to AU3 and layers below 30 cm at AU4 share Deltaproteobacteria, Planctomycetes, and Chloroflexi as dominant bacterial fractions. Furthermore, Atribacteria clearly increase below the SMTs, where methanogenesis dominates the terminal OM degradation. Additional, site-specific characteristics in subsurface layers include high percentages of Aminicenantes at AU1 and AU3, and high percentages of Spirochaetae at AU2.

Among the Archaea, the Marine Group I class (Thaumarchaeota) dominates the bioturbated layers of AU1 to AU3 (Fig. 2). In the VLs, two phylotypes, Nitrosopumilus sp. 1 [98.7% similarity to Nitrosopumilus cobalaminigenes (34)] and Nitrosopumilus sp. 2 (95.8% similarity to N. cobalaminigenes), account for 33 ± 14% and 30 ± 8% of total archaeal abundance, respectively. While Nitrosopumilus sp. 1 is virtually absent from nonbioturbated samples, Nitrosopumilus sp. 2 occurs locally at high relative abundances in nonbioturbated subsurface layers. Percentages of Woesearchaeota are high in bioturbated layers (15 ± 6%), but also in nonbioturbated layers of AU1 and AU4. Relative abundances of Bathyarchaeota (mainly Group C3) are high everywhere except in VLs of AU1 and AU3 and the surface layer of AU2. Euryarchaeota (mainly Marine Benthic Group D [MBG-D] of Thermoplasmata), Lokiarchaeota (mainly gamma subgroup), and Thorarchaeota are rare or absent within VLs, but increase and are abundant in RLs of AU1 to AU3 and in nonbioturbated sediments. Methanomicrobia are rare at all stations, except below the SMT of AU3 and AU4 (>99% ANME-1b).

Microbial Community Assembly and Its Potential Drivers Across the Four Sites.

A principal coordinates analysis (PCoA) (also known as multidimensional scaling analysis) indicates significant differences in microbial community compositions between VLs, RLs, and nonbioturbated sediments (Fig. 3). These differences are robust across different distance matrices (e.g., unweighted Unifrac, Bray-Curtis; ) and taxonomic levels (e.g., phylum, class, order, and zero-radius operational taxonomic unit [ZOTU]; ). Burrow wall communities (AU3 only) cluster with those in adjacent bulk sediment, and thus follow the division between VLs and RLs (Fig. 3 and SI Appendix, Fig. S8). By contrast, microbial communities in nonbioturbated samples are separated by sites, and show strong depth-related changes, especially in Bacteria. The only nonbioturbated sample that clusters with bioturbated samples is the bacterial community at the seafloor (0 to 2 cm) of AU4.

Fig. 3.

Ordination analyses of microbial community composition and its environmental drivers. (A) PCoA plots of bacterial and archaeal communities based on weighted Unifrac distances. Clusters consisting of samples from VLs, RLs, and nonbioturbated layers are shaded in deep blue, light blue, and red, respectively. Intercluster dissimilarities analyzed by PERMANOVA were highly significant (P < 0.001) in all cases. (B) Relationships between microbial communities and environmental variables examined by CAP. The asterisks indicate significance levels based on PERMANOVA (***P < 0.001; **P < 0.01).

A canonical analysis of principal coordinates (CAP) shows that different environmental variables are correlated with bioturbated and nonbioturbated communities (Fig. 3). Bacterial and archaeal communities in bioturbated sediments are significantly correlated with the modeled bioturbation variables (α, DB, Na), and with SO42− and %Fe(III). Furthermore, bacterial, but not archaeal, communities show significant correlations with chl a, and neither are significantly correlated with C:N ratios. In nonbioturbated sediments, bacterial and archaeal communities correlate with CH4, TOC, and δ13C-TOC. The above variables explain 68% of bacterial and 75% of archaeal total community variation. Most of this variation is accounted for by the first two axes displayed in Fig. 3, which explain 50.3% of bacterial and 54.1% of archaeal variation.

Eukaryotic 18S Gene Abundance and Community Structure.

The observed correlations with TOC, δ13C-TOC, and chl a suggest an important role for OM content, sources, and degradation state in driving microbial community structure. Here, we explore eukaryotic 18S rRNA sequences for insights into OM sources, inputs of OM by ventilation and reworking, and distributions of living eukaryotes (Fig. 4).

Fig. 4.

Depth profiles of 18S rRNA gene abundances and community composition. Photosynthetic lineages, i.e., Bacillariophyta, Dinoflagellata, Chlorophyta, and Streptophyta, are indicated by green colors.

18S rRNA gene copies decrease in the upper meter, with depth profiles at each site matching the distributions of modeled macrofaunal activity, but are similar across sites below 100 cm. At AU4, gene copy numbers drop fivefold in the top 3 cm, but only decrease slightly below. At AU1 to AU3, 18S gene copy numbers decrease by two orders of magnitude throughout the bioturbated zone, with the steepest decreases occurring within the VLs. Below the RLs, gene copy numbers decrease gradually at AU2 and AU3, and increase slightly at AU1.

18S rRNA gene sequences suggest dominant contributions of photosynthetic organisms, with the highest contributions at AU4 (>90% in all but surface sediments). At all sites, fractions of diatoms belonging to diverse Bacillariophyta decrease steeply with depth and are correlated with chl a distributions (Spearman r = 0.76, P < 0.001; ). 18S genes of the diatom Chaetoceros spp. (shown separately) are an exception. Sequences of these diatoms show less depth attenuation compared to other diatoms, and even dominate the subsurface of AU4. Dinoflagellates (Dinoflagellata, mainly Dinophyceae) have low relative abundances in the VLs but increase in the RLs and some subsurface layers. AU3 and AU4 show elevated contributions of green algae (Chlorophyta, mostly unicellular Trebouxiophyceae) and plants (Streptophyta, mainly Zostera marina, minor contributions of terrestrial angiosperms) in deeper layers, confirming the high DNA preservation potential of both phyla in marine sediments (35).

Among the nonphotosynthetic eukaryotes, protistan bacterivores, dominated by Excavata (mainly Euglenida and Kinetoplastida) and benthic Foraminifera (mainly Globothalamea and Monothalamida), are most abundant in the VLs of AU1 to AU3. These distributions are consistent with the O2 or nitrate requirements of these groups and imply that both groups may live in the VLs of these stations. By contrast, uncharacterized marine stramenopiles (MAST) have increased percentages in the nonventilated horizons of AU1, uncharacterized Cercozoa occur at high proportions throughout AU1 to AU3, and Fungi account for high fractions in nonbioturbated layers of AU1 and AU3. At AU1, metazoan 18S genes of worms (Annelida, mainly Polychaeta) and flat worms (Platyhelminthes) dominate the VLs, whereas tunicates (mainly Ascidiacea) dominate in deeper layers. The other sites have lower metazoan 18S gene contributions dominated by Nematoda (AU2 to AU4), Annelida (AU2), and Arthropoda (surface sample of AU4; mainly Copepoda) in surface sediments and by Urochordata, Annelida, Cnidaria, and Arthropoda in subsurface layers. Bivalves and brittle stars were only detected at very low (<0.1%) relative abundances, despite being dominant members of the macrofaunal communities at AU1 to AU3.

Drivers of Bacterial vs. Archaeal Community Structure in Surface and Subsurface Sediment.

Variation partitioning analyses coupled with redundancy analyses (VPA-RDA) and Mantel tests were incorporated into a path model to analyze the most important drivers of microbial community structure in surface (0 to 40 cm) vs. subsurface sediments (>40 cm; Fig. 5). Hereby the relative importance of all environmental variables tested in Fig. 3 was determined.

Fig. 5.

Conceptual path model to investigate the relationships between bioturbation activity, sediment biogeochemistry, and microbial and eukaryotic community structure: (A) surface sediments (0 to 40 cm) and (B) subsurface sediments (>40 cm). The red and blue colors distinguish solid-phase (DB, Na) and solute-phase (α) bioturbation parameters. All calculations are based on PLS regression, VPA-RDA, and Mantel test. Widths of arrows and wave lines are proportional to the explanatory power of VPA-RDA variables and Mantel coefficients, respectively. ***P < 0.001; **P < 0.01; *P < 0.5. Note: Surface and subsurface sediments are divided at 40 cm based on Fig. 1. Changing this division to other depths does not significantly alter the statistical outcomes, however ().

According to VPA-RDA, the variables that provide the most explanatory power in surface sediment account for 50% and 58% of bacterial and archaeal community variation, respectively (Fig. ). Notably, ventilation and reworking differ in their impacts on bacterial and archaeal communities. Reworking, based on nonlocal mixing (Na) and biodiffusional mixing (DB), appears to drive bacterial community structure by controlling OM compositional variables (chl a, δ13C-TOC, and TOC; r2 = 0.50, P < 0.05, partial least-squares [PLS] regression), whereas SO42− concentrations and bioirrigation rates (α) are less important. By contrast, archaeal community variation is mainly explained by bioirrigation. Mantle tests, moreover, show that eukaryotic 18S rRNA gene sequences are highly correlated with bacterial and archaeal community structure and explain a higher fraction of bacterial (r = 0.83) than archaeal community variation (r = 0.50).

In subsurface sediments, VPA-RDA suggest that TOC and δ13C-TOC, and to a lesser extent SO42− and CH4, are the main drivers of microbial community structure (Fig. 5). Measured variables explain 56% of bacterial but only 29% of archaeal community variation. According to Mantle tests, eukaryotic 18S rRNA gene sequences are again highly correlated with and explain comparable fractions of bacterial (r = 0.40) and archaeal community variation (r = 0.46).

Drivers of Microbial Taxa Composition from Bioturbated to Nonbioturbated Sediment.

To identify the most important factors that drive the relative abundances of specific microbial taxa, we examined pairwise correlations between dominant bacterial and archaeal lineages and 1) modeled bioturbation activity (α, DB, Na), 2) geochemical variables [%Fe(III), chl a, SO42−, CH4, δ13C-TOC, TOC], and 3) dominant eukaryotic taxa (Fig. 6).

Fig. 6.

Pairwise Spearman correlations of relative abundances of (A) bacterial and (B) archaeal lineages with environmental variables and relative abundances of dominant eukaryotic groups. Only significant (P < 0.05) correlations are shown. The asterisks mark which environmental variable shows the strongest correlation with each given microbial lineage. Note: Correlations of α, %Fe(III), DB, Na, and chl a with microbial lineages were only calculated for surface sediment where data exist for both.

Based on observed correlation patterns, which are supported by clustering patterns in a co-occurrence network of microbial lineages (), we propose the following categories of microorganisms: Type I “bioturbation lineages” are dominated by bacterial Bacteroidetes, Verrucomicrobia, Acidobacteria, and Gammaproteobacteria, and thaumarchaeal Nitrosopumilus sp. 1 and sp. 2. These lineages generally show significant correlations with α, DB, Na, chl a, SO42−, δ13C-TOC, and distributions of Bacillariophyta others, Foraminifera, and Excavata. Type II “ubiquitous lineages” are abundant in surface and subsurface sediments and in many cases only show weak or inconsistent trends in relation to bioturbation activity, geochemical variables, or 18S genes. This category is dominated by bacterial Delta- and Alphaproteobacteria and Planctomycetes, as well as archaeal Woesearchaeota, Diapherotrites, and Lokiarchaeota (beta subgroup). Type III “subsurface lineages” correlate most strongly with CH4 concentrations, δ13C-TOC, and/or with 18S genes of Streptophyta, Chlorophyta, Fungi, and MAST. These lineages are dominated by bacterial Chloroflexi (all but Caldilinea), Aminicenantes, Spirochaetae, and Atribacteria, as well as archaeal Bathy-, Odin-, Thor-, Eury-, and Aenigmarchaeota.

Notably, many bacterial but only few archaeal Type III lineages are correlated with δ13C-TOC. Also, while Type I Bacteria are typically negatively correlated with δ13C-TOC, Type II and Type III Bacteria are mostly positively correlated with δ13C-TOC. Furthermore, within the Lokiarchaeota, which we classify as Type II above, the dominant gamma subgroup shows correlation patterns typical of Type III lineages. Finally, several variables are only weakly or rarely correlated with dominant bacterial and archaeal taxa. These include TOC and relative abundances of Chaetoceros spp. of Bacillariophyta, Dinoflagellata, and Metazoa.

Most of the Type I “bioturbation lineages” can be further subdivided based on correlations with measured variables: Type Ia members, consisting of dominant Archaea (Nitrosopumilus sp. 1 and sp. 2) and several less dominant bacterial groups, e.g., Opitutae and Verrucomicrobiae (both Verrucomicrobia), Oceanospirillales, BD7-8 marine group, and Alteromonadales (all Gammaproteobacteria), show strongest correlations with bioirrigation rates (α) and contributions of Foraminifera and Excavata. By contrast, Type Ib members, which include Cyanobacteria, most Bacteroidetes, Verrucomicrobia, Acidobacteria, and Gammaproteobacteria, as well as several deltaproteobacterial lineages (Myxococcales, Sva0485), mainly correlate with reworking-related variables, i.e., chl a, DB, and Na, and/or with 18S gene percentages of Bacillariophyta others and Excavata. This subdivision of “bioturbation lineages” into Type Ia and Type Ib lineages is also visible in the co-occurrence network ().

Discussion

Our study indicates that bioturbation controls microbial community structure in surface sediments across continental margin sites that differ greatly in water depth, sedimentation rates, input of fresh OM, TOC content, redox conditions, microbial activity, and macrofauna (Figs. 1–3). Independent of macrofaunal abundance, community structure, or dominant mode of behavior (ventilation, reworking), bacterial and archaeal communities inhabiting bioturbated surface sediment are highly similar to each other, and differ significantly from those in nonbioturbated surface or subsurface sediments (Figs. 1 and 3). While bioturbation generally favors Bacteria over Archaea (BARs; Fig. 2), both bacterial and archaeal richness are lowest in sediments with high ventilation rates and highest in reworking-dominated sediments (Fig. 2C). Ventilation, moreover, selects for different lineages (Type Ia) than reworking (Type Ib), and is most likely the main driver of archaeal community structure in surface sediment through its impact on O2 availability and redox conditions (Figs. 5 and 6). By contrast, reworking appears to drive bacterial communities through its influence on fresh OM distributions.

In subsurface sediments, where bioturbation is absent, bacterial community structure is mainly correlated with organic carbon content and sources (TOC, δ13C-TOC) and methane concentrations, whereas the main drivers of archaeal community structure remain unclear (Figs. 5 and 6). Stronger site-specific community clustering in Archaea suggests that site-related variables may play a more important role in structuring archaeal than bacterial communities (Fig. 3). Type III “subsurface lineages” of Bacteria and Archaea show strong correlations with 18S genes of certain eukaryotic groups (Fig. 6). Archaea generally show stronger correlations with nonphotosynthetic MAST, Fungi, and Cercozoa, whereas Bacteria show stronger correlations with photosynthetic Chlorophyta and Streptophyta. The implications of these correlations are unclear, though it is worth mentioning that bacterial communities in surface sediments also show strong correlations with biomarkers of photosynthetic organisms, i.e., chl a, 18S genes of diatoms (Bacillariophyta others; Figs. 5 and 6). Yet, the existence of Type II “ubiquitous lineages” indicates that not all microbial lineages are controlled by bioturbation or OM sources. Correlations of several Type II members with SO42− or CH4 suggest possible involvement in respiration reactions or syntrophic associations with microorganisms involved in respiration reactions.

Impact of Ventilation on Microbial Community Structure.

It appears that strong ventilation activity reduces microbial richness and selects for distinct “Type Ia” aerobic and aerotolerant microorganisms that are highly similar across sites (Figs. 1 and 2). The impact of ventilation is not restricted to burrow walls but extends throughout the entire VL, as indicated by only minor differences in community structure between burrow walls and bulk sediments (Fig. 3 and ). Possibly the high rates of porewater exchange and sediment mixing cause sediments throughout the entire VL to come into contact with O2 or experience strong redox fluctuations on a frequent basis. These frequent chemical disturbances may exclude many anaerobes and select for a low diversity of redox-resilient microorganisms.

Ventilation is the main driver of archaeal community variation in surface sediments (Fig. 5), and in particular favors two phylotypes of the aerobic nitrifying thaumarchaeotal genus Nitrosopumilus (Fig. 2). The relative abundance of one of these phylotypes, Nitrosopumilus sp.1, is even highly correlated with ventilation rates (Spearman r = 0.86, P < 0.001; ). While ventilation is only a minor driver of bacterial community variation in surface sediments, relative abundances of several verrucomicrobial classes (Opitutae, Verrucomicrobiae) and gammaproteobacterial orders (Oceanospirillales, BD7-8, Alteromonadales) are significantly correlated with α. Previous marine isolates of Verrucomicrobiae and Oceanospirillaceae (dominant family of Oceanospirillales in our study) were aerobes with broad substrate spectra (including monosaccharides, amino acids, alcohols, organic acids), and some Oceanospirillaceae additionally perform nitrate reduction (36, 37). Opitutae, Psychromonadaceae (dominant family of Alteromonadales), and BD7-8 subgroup commonly live in symbiotic relations with marine sediment invertebrates and are facultatively anaerobic carbohydrate degraders (Opitutae, Psychromonadaceae) or putatively chemoautotrophic, aerobic sulfide oxidizers (BD7-8) (38–40).

Besides correlating with α, many Type Ia taxa also strongly correlate with 18S gene percentages of labile algae (Bacillariophyta others) and protists (Foraminifera, Excavata). Correlations with benthic Foraminifera and Excavata match the aerobic or facultatively anaerobic metabolism of Type Ia microorganisms and these protists. By contrast, correlations with carbohydrate- and amino acid-rich diatoms [Bacillariophyta others (41)] are consistent with the use of carbohydrates and amino acids as energy substrates by several Type Ia taxa.

While all “Type Ib” lineages and many Type II bacteria (e.g., Alpha- and Deltaproteobacteria) are also present in high read percentages throughout the VLs, high ventilation rates seem to select against Type III “subsurface lineages.” These lineages may not tolerate O2 exposure and/or the strong redox fluctuations that are the norm in bioturbated surface sediments (16). Several Type III groups, including Bathyarchaeota, Lokiarchaeota (gamma subgroup), Thorarchaeota, Euryarchaeota (MBG-D), and Chloroflexi (Anaerolineae and Dehalococcoides) already increase in percentages in deeper parts of the VLs and the underlying RL, and are abundant in nonbioturbated surface sediment of AU4.

Impact of Reworking on Microbial Community Structure.

Our study indicates that reworking strongly impacts bacterial but not archaeal community structure by mediating the transport of fresh OM into sediment (Figs. 3 and 5). Matching this trend, bacterial but not archaeal gene abundances are significantly correlated with chl a content (Spearman r = 0.86, P < 0.001) and eukaryotic gene abundances (Spearman r = 0.92, P < 0.001) in surface sediments of bioturbated sites ().

Reworking-related Type Ib bacterial taxa show strong correlations with chl a content and 18S gene percentages of diatoms (Bacillariophyta others; Fig. 6). These lineages comprise mostly aerobic, microaerobic, and facultatively anaerobic Bacteria that degrade algal polysaccharides, e.g., Flavobacteriia (42) of Bacteroidetes, Sandaracinaceae of Myxococcales (43), and organo- and photoorganotrophic Halieaceae [main family of gammaproteobacterial Cellvibrionales (44)]. Many Type Ib lineages, e.g., Flavobacteriia, Sphingobacteriia (both Bacteroidetes), and Sav0071 (Gammaproteobacteria), are, furthermore, distributed independent of bioirrigation rates (α) or redox status [%Fe(III)] (Fig. 6). In addition, several metabolically diverse lineages correlate significantly with nonlocal mixing (Na). The metabolisms of these lineages include nitrification [e.g., Nitrospinae (45), Nitrosomonadales of Betaproteobacteria (46)], aerobic methylotrophy [E01-9C-26 of Gammaproteobacteria (47)], aerobic organotrophy [Arenicella of Gammaproteobacteria (48)], anaerobic organotrophy [Bacteroidia of Bacteroidetes (49)], and sulfate reduction [Sva0485 of Deltaproteobacteria (50)].

The diversity of aerobic, microaerobic, and anaerobic metabolisms among Type Ib lineages, and the fact that many Type III lineages already emerge in the RLs, explains the high diversity in sediment layers where reworking is the dominant macrofaunal activity. We suggest that the decreased O2 input but continued significant input of fresh OM results in a spatially and temporally heterogeneous mosaic of oxic to anoxic conditions and OM ranging from fresh detritus to diagenetically altered, low-reactivity compounds. The richness in chemical microenvironments results in a high diversity of metabolic niches that support the coexistence of physiologically diverse Type I bioturbated, Type III subsurface, and Type II ubiquitous lineages within the same macroenvironment.

Changing Communities from Surface to Subsurface Layers.

The shift toward Type III subsurface lineages in nonbioturbated sediment is accompanied by a decrease in microbial richness as Type I lineages, which appear linked to O2 and fresh OM input, decrease in relative abundance (Figs. 2 and 6). At AU4, where O2 is absent, Type III Archaea (gamma subgroup of Lokiarchaeota, Bathy- and Thorarchaeota, Thermoplasmata) and Bacteria are present at the seafloor (Anaerolineae, Dehalococcoidia) or emerge at depths where labile phytodetritus (chl a, Bacillariophyta others) is depleted (Atribacteria) (Figs. 2 and 4).

Despite these consistent changes from bioturbated to nonbioturbated sediments, the reasons for the high site dependency of microbial community structure in nonbioturbated sediments remain unclear. Percentages of many Type II taxa are significantly correlated with SO42− concentrations, whereas percentages of numerous Type III taxa are significantly correlated with CH4 concentrations (Fig. 6). The correlations of Methanomicrobia with CH4 concentrations and of Desulfobacterales with SO42− concentrations match the known CH4-cycling and SO42−-reducing metabolisms of these taxa, respectively. In other cases, syntrophic partnerships with sulfate reducers or methanogens offer a potential explanation.

OM composition may also be important, given that many subsurface microorganisms gain energy by fermentative or acetogenic breakdown of OM (6–10). Consistent with this notion, a major fraction of subsurface bacterial community variation can be explained with TOC and δ13C-TOC (Fig. 5). Moreover, δ13C-TOC and percentages of major 18S lineages correlate significantly with major Type III groups (Fig. 6). Dominant bacterial lineages degrade aromatic compounds and fatty acids [Dehalococcoidia (8)], or ferment carbohydrates or proteins [Anaerolineae (51), Omnitrophica (52), Aminicenantes (53), Spirochaetes (54)]. Atribacteria have been linked to syntrophic growth on organic acids with methanogenic partners (55) and to fermentation of sugars and organic acids (9). Dominant Archaea have been linked to protein fermentation [MBG-D (7), Thorarchaeota (56)], fermentative or acetogenic breakdown of proteins, carbohydrates, and lignin [Bathyarchaeota (10)], fermentation of unknown substrates [Aenigmarchaeota (52)], and syntrophic growth on amino acids [gamma subgroup of Lokiarchaeota (57)].

Conclusions and Further Implications.

Our study indicates that macrofaunal bioturbation controls microbial community assembly at all phylogenetic levels and promotes dominance of Bacteria over Archaea in surface sediments from diverse continental margin settings. Burrow ventilation appears to control archaeal community structure mainly by regulating O2 input. By contrast, reworking appears to be the main driver of bacterial community structure, presumably by controlling the distribution of reactive OM. The reasons for bacterial dominance over Archaea remain uncertain, but may be related to higher metabolic versatility among surface sedimentary Bacteria. Certain Bacteria have been shown to thrive in natural and engineered habitats under oxic–anoxic fluctuations by means of their ability to rapidly shift between aerobic respiration and fermentation (58, 59). Similar adaptations may confer fitness advantages to certain groups of Bacteria in bioturbated sediments.

While literature data suggest that the energy substrates of major bacterial and archaeal lineages overlap at the compound level, relative abundances of many bacterial classes are strongly correlated with 18S percentages of photosynthetic eukaryotes. By contrast, archaeal lineages, in particular Type III subsurface lineages, tend to be more strongly correlated with nonphotosynthetic MAST, Fungi, and Cercozoa. Future investigations will reveal whether these trends reflect general differences in detrital food sources among Bacteria and Archaea.

Materials and Methods

Sampling.

All samples were taken during a cruise of the R/V Aurora in August to September 2014. The top 50 cm of sediment were collected using a Rumohr corer, which is a lightweight gravity corer that enables access to nearly undisturbed surface sediments (60). All deeper layers were sampled using a gravity corer (down to ∼500 cm). Sediments for faunal analysis were sampled using a box corer (40 cm × 40 cm × 60 cm). For detailed sampling scheme, see .

Macrofaunal Analysis.

Faunal sampling, identification, and quantifications were as outlined by Kristensen et al. (26). Sediments from two to three box cores per station were sieved through 1-mm mesh on board, and the retrieved material preserved in 4% formaldehyde. Fauna were separated from debris in the laboratory and stored in 70% ethanol for later counting, identification to lowest possible taxon, and classification into functional groups based on faunal mobility and feeding behavior. Macrofaunal biomass was estimated based on abundance data and known relationships between water depth and average individual biomass from the Kattegat–Skagerrak area (26, 61, 62).

DNA Extraction.

DNA was extracted from unhomogenized ∼0.2 g of wet sediment following lysis protocol II of the modular method of Lever et al. (63), which combines chemical (lysis solution I) and mechanical cell lysis (bead-beating: 0.1-mm Zirconium beads), 2× washing with chloroform:isoamyl alcohol (24:1), and precipitation with a mixture of linear polyacrylamide, sodium chloride, and ethanol. DNA was then purified according to protocol A of the CleanAll DNA/RNA Clean-Up and Concentration Micro Kit (Norgen Biotek Corporation). To minimize adsorptive losses of DNA during the extraction, samples were homogenized with 0.1 mL of 10 mM sodium hexametaphosphate solution prior to cell lysis.

Quantitative PCR.

Abundances of bacterial and archaeal 16S rRNA genes and eukaryotic 18S rRNA genes in DNA extracts were quantified on a LightCycler 480 II (Roche Life Science) by SYBR-Green I-based quantitative PCR (qPCR). The primer pairs for Bacteria, Archaea, and Eukarya were Bac908F_mod (5′-AACTCAAAKGAATTGACGGG-3′) (63)/Bac1075R (5′-CACGAGCTGACGACARCC-3′) (64), Arch915F_mod (5′-AAT​TGG​CGG​GGG​AGC​AC-3′) (65)/Arch1059R (5′-GCCATGCACCWCCTCT-3′) (66), and All18SF_mod1 (5′-TGCATGGCCGTTCTTAGT-3′)/All18SR_mod1 (5′-CTA​AGG​GCA​TCA​CAG​ACC-3′) (ref. 35, modified from ref. 67), respectively. Plasmids of 16S rRNA genes from Holophaga foetida (Acidobacteria) and Thermoplasma acidophilum (Euryarchaeota) and of 18S rRNA genes from Tubifex (Oligochaeta) were applied as bacterial, archaeal, and eukaryotic standards, respectively (further information in ).

Sequencing and Bioinformatic Analyses.

Sequence libraries were prepared according to a standard workflow (). 16S rRNA gene amplicons obtained with the primer pairs S-D-Bact-0341-b-S-17 (5′-CCTACGGGNGGCWGCAG-3′)/S-D-Bact-0785-a-A-21 (5′-GACTACHVGGGTATCTAATCC-3′) (68) for Bacteria and S-D-Arch-0519-a-A-19 (5′-C AGCMGCCGCGGTAAHACC-3′) (69)/967Rmod (5′-GTG​CTC​CCC​CGC​CAA​TT-3′) (65) for Archaea were sequenced via the MiSeq platform (Illumina). Raw sequence reads were quality-checked by FastQC (www.bioinformatics.babraham.ac.uk/projects/fastqc), read ends were trimmed using seqtk (https://github.com/lh3/seqtk), paired end reads were merged into amplicons by flash (70), primer sites were trimmed by usearch (71), and quality filtering was done by prinsEq (72). ZOTUs were generated using the UNOISE3 algorithm (73) and clustered using a 97% identity threshold. Bacterial 16S genes and eukaryotic 18S genes were taxonomically classified using the SILVA 16S database [release 128 (74); confidence threshold, 0.7] and Protist Ribosomal Reference database [PR2 v.19 (75); confidence threshold, 0.9], respectively. Archaeal 16S genes were assigned in ARB (www.arb-home.de) using neighbor-joining phylogenetic trees that were based on a manually optimized SILVA database that was expanded to include 16S gene sequences from whole-genome sequencing studies.

Geochemical Analyses.

Sediment chl a and pheophytin from ∼1.5 g wet sediment were extracted in 90% acetone and quantified spectrophotometrically following an acidification method (76). To measure TOC, total nitrogen (TN), δ13C-TOC, and δ15N-TN isotopic compositions, sediment samples were dried, decarbonized, and homogenized. Analyses were done on 5-30 mg of homogenized sediments by elemental analyzer/isotope ratio mass spectrometry, as described previously (77). Porosity and wet density of sediments were calculated from the wet and dry weights of 2-cm3 sediments. CH4 concentrations were measured by gas chromatography as published previously (78). Porewater concentrations of SO42− and DIC were measured by ion chromatography and gas chromatography and were published previously (79), as were bioavailable Fe(II) and Fe(III) contents, which were extracted using 0.5 M HCl and measured spectrophotometrically (26).

Radionuclides.

Two HPGe gamma spectrometry systems (ORTEC GMX-120265 and GWL-100230) were used for radionuclides analysis (210Pb, 46.52 keV; 214Pb, 351.99 keV; and 137Cs, 661.62 keV) (80). Efficiencies of gamma detectors were calibrated using IAEA reference materials (RGTh and RGU) and an in-house secondary standard (“Rock-Falling Mountain Soils”; RMC-AEC Taiwan). Excess 210Pb was calculated by subtracting 214Pb (used as an index of 226Ra, i.e., supported 210Pb) from total 210Pb. The activities of radionuclides were decay-corrected to the sampling date. All radionuclide data were calculated on a salt-free dry-weight basis.

Transport–Reaction Models.

Biological mixing rates of sediment particles were modeled based on the steady-state conservation equation for radiotracers (210Pb) (81) with consideration of biodiffusional and nonlocal mixing effects:

where x is depth; , the concentration of excess 210Pb; , the solid volume fraction; , the constant flux of solid sediment to the seabed; , the solid-phase density; and , the excess 210Pb concentration at; is the 210Pb decay constant. Fixed concentration and zero concentration gradients were imposed at the sediment–water interface and at 50-cm depth, respectively. Biodiffusional mixing and nonlocal mixing were described by the following:

Parameters that describe the mixing coefficients at , mixing depth (, ), and spread (, ) were adjusted to reproduce the observed excess 210Pb profiles. The bioirrigation coefficients (α) were calculated by Kristensen et al. (26), by fitting a double-exponential nonlocal exchange function (similar to Eq. ) to the measured DIC profile.

Multivariate Statistics.

Statistical analyses were performed in R (www.R-project.org) and based on ZOTU-level phylogenetic assignments unless stated otherwise. Chao1 richness, and PCoA and CAP coordinates using Unifrac and Bray–Curtis distances were calculated using the “phyloseq” package (82). In the CAP analysis, α, DB, Na, and %Fe(III) were assumed to be zero at >40 cm. Using the “vegan” package (83), we performed 1) a PERMANOVA analysis (permutations = 999) to examine the (dis)similarity between different groups of samples; 2) VPA-RDA to determine how environmental factors individually or jointly contribute to microbial community variation; and 3) Mantel tests to examine correlations between microbial and eukaryotic communities. Using PLS regressions in the R package “plsm,” we examined impacts of reworking activities on variations of OM distributions (84). A cooccurrence network of bacterial and archaeal lineages was constructed using Cytoscape (https://cytoscape.org/).

Data Availability.

Gene sequences are available online through the National Center for Biotechnology Information, https://www.ncbi.nlm.nih.gov (bacterial and archaeal 16S, BioProject PRJNA565394; eukaryotic 18S: Nucleotide accession nos. MN487107-MN488526). All geochemical, qPCR, and modeled bioturbation rate data are included in Dataset S1.