Single-celled bioturbators: benthic foraminifera mediate oxygen penetration and prokaryotic diversity in intertidal sediment

Bioturbation processes influence particulate (sediment reworking) and dissolved (bioirrigation) fluxes at the sediment–water interface. Recent works showed that benthic foraminifera largely contribute to sediment reworking in intertidal mudflats, yet their role in bioirrigation processes remains unknown. In a laboratory experiment, we showed that foraminifera motion behaviour increased the oxygen penetration depth and decreased the total organic content. Their activity in the top 5 mm of the sediment also affected prokaryotic community structure. Indeed, in bioturbated sediment, bacterial richness was reduced, and sulfate-reducing taxa abundance in deeper layers was also reduced, probably inhibited by the larger oxygen penetration depth. Since foraminifera can modify both particulate and dissolved fluxes, their role as bioturbators can no longer be neglected. They are further able to mediate the prokaryotic community, suggesting that they play a major role in the benthic ecosystem functioning and may be the first described single-celled eukaryotic ecosystem engineers.

1 Introduction

Intertidal mudflats are among the most productive ecosystems on Earth (Heip et al., 1995). Given their natural features, they are zones of prime importance for organic matter (OM) accumulation (Jickells and Rae, 1997), which can sequester more than 200 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ (Chmura et al., 2003). Mudflat sediments usually host intense biological activity, and OM is rapidly mineralized (Mayor et al., 2018) via a series of diagenetic reactions from oxygen respiration to methane production (Froelich et al., 1979). In such cohesive environments, dissolved oxygen (O₂) is usually only available in the top millimetres of the sediment, and transport of solutes is assured by molecular diffusion (Aller, 1988).

Burrow-dwelling macro-invertebrates (organisms larger than 500 µm) greatly influence intertidal mudflat functioning through bioturbation (Meysman et al., 2006) – a process which combines sediment reworking (i.e. transport of particles) and burrow ventilation (which causes bioirrigation: the transport of water and solutes; see the review in Kristensen et al., 2012). The effect of bioturbation by macro-invertebrates on the benthic ecosystem functioning is mediated by complex interactions with meiofaunal organisms (organisms smaller than 500 µm; Piot et al., 2014; Lacoste et al., 2018; Schratzberger and Ingels, 2018). Indeed, meiofauna may also contribute significantly to sediment reworking (Bradshaw et al., 2006) and bioirrigation (Cullen, 1973; Aller and Aller, 1992). Noticeably, meiofauna was reported to improve sediment oxygenation and sulfide removal (Bonaglia et al., 2020), to affect the nitrogen cycle by stimulating nitrate reduction (Prast et al., 2007; Bonaglia et al., 2014), and to enhance OM mineralization (Rysgaard et al., 2000; Nascimento et al., 2012). Meiofaunal bioturbation can further lead to changes in the abundances of all and specific groups of bacteria in sediments (Prast et al., 2007; Lacoste et al., 2018; Bonaglia et al., 2020), but these studies did not evaluate its effect on the whole bacterial and archaeal community structures. Bioturbation by macro-invertebrates may significantly impact bacterial community structure by modifying biogeochemical gradients and by modifying the availability and quality of OM (e.g. mucus production) in sediments (Papaspyrou et al., 2006; Cuny et al., 2007). For example, Laverock et al. (2010) demonstrated that bacterial communities from irrigated burrows of the ghost shrimp (Upogebia deltaura and Callianassa subterranea) were more diverse than bacterial communities from non-bioturbated sediments. In this context, it can be expected that bioirrigation by meiofauna would similarly increase oxygen availability in sediments overall, hence favouring aerobic prokaryotes over strictly anaerobic species in sediments but also increasing the sediment heterogeneity, enhancing microbial diversity.

In spite of their role in benthic ecosystem functioning (Moodley et al., 2000; Geslin et al., 2011), the role of foraminifera as bioturbators remains a fairly untapped question, with only a few pioneer studies looking at how their displacement may affect the sediment reworking process (Severin et al., 1982; Hemleben and Kitazato, 1995; Groß, 2000). Noticeably, their ability to move in the sediment column affects the surface sediment cohesiveness (Cedhagen et al., 2021) and contributes to the horizontal and vertical transport of sediment particles (Groß, 2002; Deldicq et al., 2020, 2021, 2023). Consequently, foraminifera are assumed to affect sediment porosity and allow for “good sediment ventilation” (Hemleben and Kitazato, 1995; Groß, 2002). Supporting this assumption, foraminiferal activity was shown to affect dissolved cadmium concentrations in the pore water and overlaying water (Green and Chandler, 1994), suggesting that foraminifera influence the water and solute exchanges at the sediment–water interface. However, studies based on two-dimensional oxygen measurements did not report a positive effect of foraminifera on dissolved oxygen concentrations in sediments, as their aerobic respiration produced a decrease of oxygen concentration in foraminiferal burrows (Oguri et al., 2006; Heinz and Geslin, 2012).

In this context, it appears critical to further describe the role of foraminifera in bioirrigation processes and quantify their contribution to solute fluxes at the sediment–water interface. To do so, the impact of foraminiferal displacements in the sediment matrix was assessed on (1) the oxygen vertical distribution in homogenized sediment; (2) the subsequent oxygen fluxes at the sediment–water interface; (3) the resulting influences on OM content (total organic carbon and total nitrogen); and (4) the prokaryotic (archaea and bacteria) community structure to ultimately determine their role in bioirrigation processes, OM mineralization, and the microbenthic communities.

2 Material and methods

2.1 Sediment and living foraminifera collection

Surface sediment (top 10 mm) from Authie Bay (northern France, English Channel; 50^∘22^′20^′′ N, 1^∘35^′45^′′ E) was collected in January 2018 and kept frozen in the dark at −20 ^∘C to kill any potential bioturbators before being used in the experimental cores.

Living foraminifera were extracted from surface sediment (top 10 mm, sieved over a 125 µm mesh) collected in the Boulogne-sur-Mer harbour (50^∘43^′04^′′ N, 1^∘34^′26^′′ E) in November 2019. Only active individuals (i.e. leaving a displacement track on a thin layer of sediment) were selected for the experiment.

Figure 1 (a) Schematic representation of experimented control cores and cores with foraminifera (side view) at the beginning (day 0) and the end of the experiment (day 85), with cores sampled for measurements of microporosity (grey), prokaryotic diversity and sediment TOC and TN (brown), foraminifera survival (green), and O₂ microprofiling (blue). (b) Location of the five microprofiling zones (top view of the cores) and (c) picture of the cores placed in the aquaria during oxygen microprofiling.

2.2 Experimental design

A total of 17 cores (45 mm height and 10 mm × 10 mm square section, Fig. 1a) were filled with homogenized thawed Authie Bay sediment (sediment was defrosted and stirred in a glass beaker before being transferred in the cores), placed in an air-bubbled 7 L aquarium (closed system filled with 35 PSU unfiltered English Channel seawater), and left for 14 d prior to adding foraminifera to give enough equilibration time to establish redox fronts and microbial processes in the sediment column. The experiment was carried out for 85 d in the dark (with a photosynthetic active radiation < 0.7 $\mathrm{\mu }\mathrm{mol}\mathrm{photons}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$; SA-190 quantum sensor, LI-COR) in a temperature-controlled room (at 18 ± 1 ^∘C).

Oxygen microprofiles were realized in control cores (n=6, without any foraminifera) and cores with foraminifera (n=6, abundance = 30 individuals cm⁻²). Foraminiferal species composition (78 % Haynesina germanica, 10 % Ammonia tepida, 8 % Quinqueloculina seminulum, and 4 % Cribroelphidium excavatum per core) and abundance selected for the experiment were chosen based on their natural densities and species composition in local mudflats (Francescangeli et al., 2020). From these 12 cores, 3 control cores and 3 cores with foraminifera were randomly selected at the end of the experiment to evaluate the influence of foraminifera on organic matter (OM) content and microbial community structures at two sediment depths (0–5 and 5–10 mm). The remaining three cores with foraminifera were used to determine the foraminiferal survival rate. Eight cores containing no foraminifera were dedicated solely to microporosity measurements at the beginning (n=4) and at the end of the experiment (n=4).

2.3 Foraminifera survival

At the end of the experiment, three cores with foraminifera were placed in a 1 µmol L⁻¹ CellHunt Green CMFDA solution (5-chloromethylfluorescein diacetate, Setareh Biotech) for 24 h, fixed with 70 % ethanol, and sieved over a 125 µm mesh (Choquel et al., 2021; Langlet et al., 2013). Foraminifera exhibiting a bright fluorescence under an epifluorescence stereomicroscope (Olympus SZX16 with a fluorescent light source Olympus KL1600pE-300) at 492 nm excitation and 517 nm emission wavelength (Langlet et al., 2014) were picked and identified to determine the foraminiferal survival rate.

2.4 Organic matter measurements

The total organic carbon (TOC) and total nitrogen (TN) contents of sediment samples were measured in two subsamples following the capsule method (Brodie et al., 2011). They were determined by high-temperature combustion of pre-acidified (HCl, 2N) dry samples (60 ^∘C, 48 h) and subsequent measurement of CO₂ and N₂ by thermal conductometry using an elemental analyser (FlashEA, Thermo Electron Corporation). The average differences between two subsamples were 0.06 % and 0.007 % for TOC and TN, respectively. Since the sediment was homogenized before the experiment, we assume that initial OM content was the same in control and cores with foraminifera.

2.5 Microporosity measurements

At the beginning and at the end of the experiment, four sediment cores were frozen at −20 ^∘C and sliced with a razor blade from 0 to 10 mm depth with a 1 mm vertical resolution to measure water content. For each slice of sediment, we measured on a precision microbalance (Sartorius R160P) the humid (m_h) and dry (m_d) masses (before and after drying at 40 ^∘C for 48 h) to determine water mass (m_w such as $m_{\mathrm{w}}=m_{\mathrm{h}}-m_{\mathrm{d}}$) and calculate the sediment microporosity (Φ) with ϱ_w = 1.035 and ϱ_s = 2.65 the density of water and sediment, respectively (Berner, 1980). Microporosity vertical distribution was modelled following an exponential decrease with depth (Fig. 2 in the Supplement). To estimate microporosity at each sampling time, we assumed that it was decreasing linearly with time.

2.6 Pore-water dissolved oxygen distribution

2.6.1 Sampling strategy

At each measurement time (from 1 d before adding foraminifera to 85 d after introduction of the living foraminifera), two cores containing foraminifera and two control cores were randomly selected, and three oxygen microprofiles were realized in each core. Each core was subdivided into five zones (Fig. 1b) sampled at five different times to ensure that microprofiling was not realized twice in the same area (Fig. 1c). All measuring cores and zones were selected randomly to minimize any potential effect of microtopography and core-specific responses (Table 1 in the Supplement).

2.6.2 Oxygen microprofiling

At each sampling time, a 50 µm diameter tip Clark type microelectrode (Revsbech, 1989) (Unisense, Denmark) was two-point-calibrated using the overlying water in the air-bubbled aquarium as a 100 % saturation reference and the signal at 10 mm depth in the experiment sediment as an anoxic reference. Oxygen concentration at 100 % saturation in 18 ^∘C and 35 PSU seawater was 239.7 µmol L⁻¹. The microsensor was placed on a motorized micromanipulator (Unisense, Denmark), and vertical profiles were realized from about 2 mm above the sediment–water interface down to the anoxic zone of the sediment with a 150 µm vertical resolution. Three microprofiles were realized in each selected zone, and the distance between two replicate profiles ranged from about 1 to 2 mm.

2.6.3 Oxygen profile interpretation

The oxygen penetration depth (OPD) was selected as the shallowest point with a dissolved oxygen concentration lower than 1 µmol L⁻¹ (Bonaglia et al., 2020).

We computed diffusive oxygen uptake (DOU) following Berg et al. (1998), using their Eqs. (1)–(10). We minimized the cost function, which includes data from the three replicates, using the L-BFGS-B algorithm (Byrd et al., 1995) with bounds to ensure that production remained negative. Berg et al. (1998) employed the stepwise regression algorithm that results in piecewise constant “production zones” (their Eq. 11) to limit the complexity of the model. Instead, we regularized the total variation (i.e. the sum of the absolute first-order derivative) using the elastic net algorithm (Rudin et al., 1992). Like the number of zones in Berg et al. (1998), the regularization intensity is a hyperparameter that controls the complexity (i.e. smoothness) of the optimized profile. We provide the algorithm, data, and Jupyter Notebook to reproduce our analysis (see the Supplement).

We imposed nil oxygen concentration and nil DOU in the sediment at the bottom of the calculation zone (Bonaglia et al., 2014). The sediment diffusion coefficient (Ds) was calculated using the microporosity (Φ) measurements (Ds = D₀⋅Φ²; Ullman and Aller 1982) and a D₀ coefficient of 1.854 × 10⁻⁵ cm² s⁻¹ (oxygen diffusion coefficient at 18 ^∘C and 35 PSU).

2.7 Prokaryotic diversity

At the end of the experiment, three sediment cores with foraminifera and three control cores were frozen at −20 ^∘C and sliced with a sterile razor blade in two 5 mm depth intervals (0–5 and 5–10 mm). For each sample, DNA was extracted from 0.25 g of wet sediment using the ZymoBIOMICS™ DNA Miniprep Kit (Zymo Research), according to the manufacturer's instructions. The quantity and the quality of extracted DNA were quantified and controlled using PicoGreen and a capillary electrophoresis (QIAxcel), respectively. The V3–V5 hypervariable regions of the 16S gene were amplified to target the bacterial community and archaeal community and to evaluate the respective abundances of archaea and bacteria in sediments. Amplifications were done using the following primer pairs: 357F_ILMN (5^′-CCTACGGGAGGCAGCAG-3^′) and 926R_ILMN (5^′-CCGYCAATTYMTTTRAGTTT-3^′) for bacteria, 519F_ILMN (5^′-CAGCMGCCGCGGTAA-3^′) and 915R_ILMN (5^′-GTGCTCCCCCGCCAATTCCT-3^′) for archaea, and 515F_ILMN (5^′-GTGYCAGCMGCCGCGGTA-3^′) and 909R_ILMN (5^′-CCCCGYCAATTCMTTTRAGT-3^′) for relative abundances of archaea and bacteria. The first PCR (PCR1) was performed with 35 cycles at 50 ^∘C for bacteria and at 58 ^∘C for archaea and relative abundances. Each PCR1 was performed in a 25 µL reaction volume, using “5x HOT BIOAmp^® Blend Master Mix” DNA Polymerase, 2 µL of DNA template, 0.24 µmol L⁻¹ reverse and forward primers, MgCl₂ at 12.5 mmol L⁻¹, bovine serum albumin at 20 mg mL⁻¹, “10x GC-rich Enhancer”, and nuclease-free water. Thermal cycles were as follows: 95 ^∘C for 3 min (95 ^∘C for 30 s, 55 ^∘C for 30 s, 72 ^∘C for 1 min) 25 times and 72 ^∘C for 5 min. The PCR was replicated three times for the 12 samples and 2 controls (extraction and PCR controls) for each couple of primers. Amplification replicates were then pooled and purified using Agencourt AMPure XP beads. A second PCR (using PCR1 as a DNA template) with 15 cycles for bacteria and archaea and 10 cycles for relative abundances was performed for sample indexing (indexes+P5/P7). PCR2 products were also purified with AMPure beads. Then, DNA was quantified using the QuantiFluor dsDNA kit (Thermo Fisher). All samples were pooled in equimolar proportions and sequenced on an Illumina MiSeq platform with 5 % PhiX (Flow Cell V3, paired-end 2 × 300 bp) by Biofidal (Vaulx-en-Velin, France).

Bioinformatic processing of the merged 2 × 300 bp paired-end reads followed sequential steps: (1) dereplication and filtering (keeping only 300 to 500 bp long reads containing a valid mismatch-free tag and no ambiguous base), (2) clustering into operational taxonomic units (OTUs) with SWARM (Mahé et al., 2014) (two-step procedure: local clustering threshold d=1 and then d=3), (3) removal of chimera, (4) removal of OTUs detected in only one out of three replicates from same condition, (5) abundance normalization (by rarefaction, i.e. subsampling at 33 885 reads per sample for bacteria, 33 834 reads per sample for archaea, and 15 645 reads per sample for respective abundances, to correct for variability in sequencing depths among samples), and (6) taxonomic affiliation against the 16S SILVA database release 138 (Quast et al., 2013) based on NCBI blastn+ (Altschul et al., 1990) and allowing for multiple affiliations. These different steps were performed with FROGS (Find, Rapidly, OTUs with Galaxy Solution; Escudié et al., 2018) on the Galaxy web platform (Afgan et al., 2018) of the Pôle Rhône-Alpes de Bioinformatique. The OTU abundance tables and taxonomic assignments produced at this stage were then analysed using the vegan R package (Oksanen et al., 2020) to calculate alpha diversity indexes (OTU richness and Shannon diversity index).

2.8 Statistical analysis

Since oxygen microprofiles were measured several times in the same core, we chose to analyse the effect of foraminiferal bioturbation using linear mixed-effect models (Pinheiro and Bates, 2000) with “core” as a random effect in all models. Oxygen penetration depth (OPD) and dissolved oxygen uptake (DOU) were set as response variables, while experimental duration (time), treatment (control or with foraminifera), and time–treatment interaction were selected as fixed effects. Preliminary segmented analysis showed a shift in oxygen conditions between −1 and 1 d, 1 and 9 d, 9 and 22 d, 22 and 55 d, and 55 and 85 d; hence, modelling was performed on data acquired from these five time intervals separately. Due to the peculiar shape of the oxygen distribution profiles, data acquired on day 5 (zones J4, K2, D2, and F2) in both controls and cores with foraminifera were removed from the analysis (see Fig. 1 in the Supplement).

The influence of the sediment layer and treatment on sedimentary bacterial (or archaeal) community structure was visualized using a non-metric multidimensional scaling (NMDS) performed with data of OTU abundances obtained from the different cores. Differences in bacterial (or archaeal) community structures between sediment layers and treatments were tested using permutational multivariate analyses of variance (PERMANOVA; Anderson 2001). Statistical tests were based on 999 permutations of the Bray–Curtis matrix.

To determine whether the experiment affected strictly anaerobic micro-organisms, supplementary analyses were performed on bacterial taxa involved in sulfate reduction and archaeal taxa involved in methane production. Three sulfate-reducing bacterial orders (Desulfobacterales, Desulfovibrionales, and Syntrophobacterales) were selected based on the literature (Wasmund et al., 2017). Their relative abundances (proportion of reads) in bacterial communities were determined for each sample. The same procedure was applied to the relative proportion of methanogens from three orders of archaea (Methanobacteriales, Methanosarcinales, and Methanomicrobiales). Relative abundances of sulfate reducers and methanogens were logit-transformed to normalize their distributions.

The influence of sediment depth (0–5 and 5–10 mm) and treatment (control or with foraminifera) on TOC and TN content was tested using linear mixed-effect models, with sediment layer, treatment, and the sediment layer × treatment as fixed effects and the core identification as a random effect.

Differences in bacterial and archaeal diversity indexes (OTU richness and Shannon diversity index), sulfate-reducing bacteria, and methanogenic archaea were tested using a two-way ANOVA (ANOVA2), with sediment layer and treatment as the main effects. For all variables, the normality and the homoscedasticity of the residues were tested using the Shapiro–Wilk test and Levene's test, respectively. Bacterial and archaeal richness data were log-transformed before statistical analyses using a two-way ANOVA to meet these assumptions. An additional Pearson test was carried out to quantify the correlation between bacterial richness and TOC. Data analysis was carried out in R v.3.5.3 using segmented, nlme, ade4, and vegan packages (Pinheiro and Bates, 2000; Dray and Dufour, 2007; R Core Team, 2019; Oksanen et al., 2020; Muggeo, 2008).

3 Results

3.1 Foraminiferal activity observations

Non-quantitative observations showed sediment displacement at the sediment surface as well as burrow formation on the sides of sediment cores down to about 7 mm depth. Newly formed burrows were frequently observed during the first 3 weeks of the experiment, but no new burrows were found after 3 weeks. Investigation of the CellHunt Green-labelled sediment at the end of the experiment showed 19, 22, and 26 living foraminifera corresponding to a survival rate of 63 %, 73 %, and 87 % in the three tested cores.

Figure 2 Percentages of total organic carbon and total nitrogen per sediment dry mass for control (open black circles) and foraminifera (open red circles) treatments in two sediment layers sampled at the end of the experiment (85 d) in three replicate cores.

3.2 Sediment organic carbon and total nitrogen content

At the end of the experiment, total organic carbon (TOC) content ranged from 1.4 % to 1.7 %, and total nitrogen (TN) ranged from 0.21 % to 0.27 % (Fig. 2). In the top sediment layer (0–5 mm) TOC was significantly lower in the cores with foraminifera than in the control cores (1.4 % ± 0.05 % standard deviation and 1.6 % ± 0.07 %, respectively), while no significant differences were observed in the 5–10 mm layer (interaction “treatment × sediment layer”, F_(1,8) = 35.6 and p < 0.05). Similarly, TN was significantly lower in the top layer of the cores with foraminifera than in the control cores (0.2 % ± 0.01 and 0.3 % ± 0.01 %, respectively), while no effect of foraminifera was observed in the deeper sediment layers (F_(1,8) = 21.1 and p < 0.05).

3.3 Oxygen distribution in the sediment

Replicated dissolved oxygen microprofiles were homogeneous within each sampling zone, and modelled oxygen profiles used for dissolved oxygen uptake (DOU) estimates showed a good fit with the measured data (R² > 0.97; Fig. 1 in the Supplement).

Figure 3 Changes of the oxygen penetration depth (a) and dissolved oxygen uptake (b) with sampling time in the control (black) and bioturbated cores (red). To visually differentiate the otherwise identical values, a small amount of noise was added to the data (with a jitter factor of 0.5 on both x and y axes). Lines were plotted based on the linear model estimates (see Table 1) and drawn as solid or dashed lines when the time × treatment variable was significant or insignificant (at a 0.05 threshold), respectively. Dashed vertical lines delimit the time intervals selected by the segmented analysis.

Table 1 Statistics of the effect of the experimental treatment on the sediment oxygen parameters. Results of the statistical analysis (linear mixed-effect models with “core” as random effect) for the oxygen penetration depth and diffusive oxygen uptake in the five time intervals selected by the segmented analysis. Explanatory variables showing a significant effect on the response variable (p < 0.05) are shown as bold characters.

During the first 36 d of the experiment, oxygen penetration depth (OPD) ranged from 2.1 to 3.6 mm in the control cores and from 2.4 to 4.2 mm in the cores with foraminifera (Fig. 3a). Linear mixed-effect models showed a significant effect of the treatment in the 9 to 22 d time interval (Table 1); i.e. OPD was on average 350 µm larger in cores with foraminifera than in control cores. In the 22 to 55 d interval, treatment and its interaction with time showed a significant effect, such as the average difference between control and cores with foraminifera being about 300 µm and tending to reduce with time to reach similar values at 55 d (Table 1, Fig. 3a).

After 55 d, OPD ranged from 3.6 to 4.5 mm (Fig. 3a) and did not show any significant differences between the cores with foraminifera and the control cores (Table 1).

DOU ranged from 2.0 to 11.7 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ (Fig. 3b) and was significantly influenced by treatment in the 9 to 22 d interval (such as average DOU being 7.7 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ in control cores and 4.9 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ in cores with foraminifera; Table 1) and by treatment and its interaction with time in the 22 to 55 d time interval, such as the maximal difference between the two treatments of 4.6 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ at 22 d being reduced to close to 0 at 55 d. After 55 d, DOU ranged from 2.0 to 5.2 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$ and did not significantly differ between treatments (Table 1).

3.4 Prokaryote community structures

Bacterial communities dominated prokaryotic communities, with more than 97 % of reads corresponding to bacterial OTUs and less than 3 % of reads related to archaeal OTUs. The relative abundance of bacterial OTUs in prokaryotic communities significantly increased with depth, with 97 % of bacteria in the 0–5 mm sediment layer and 99.5 % in the 5–10 mm sediment layer (ANOVA2, depth effect, F_(1,8)67.3, p < 0.001). Furthermore, bacterial richness was positively correlated to TOC (R² = 0.46, p < 0.01).

Figure 4 (a) Bacterial community structure (relative abundance) in control (C) and bioturbated (F) cores in two sediment layers. Community structure is represented at the phylum level. Only phyla representing at least 1 % of the reads in at least one sample are represented. (b) Non-metric multidimensional scaling of the bacterial communities in control (“C” and open black circles) and bioturbated (“F” and open red circles) cores in the 0–5 mm (“H1” labels) and 5–10 mm (“H2” labels) sediment layers. (c) Richness and diversity of bacterial communities in the different sediment layers in three control (open black circles) and three bioturbated (open red circles) cores.

Figure 5 (a) Proportion of sulfate-reducing bacteria (Desulfobacterales, Desulfovibrionales, and Syntrophobacterales) and (b) proportion of methanogenic archaea (Methanobacteriales, Methanosarcinales, and Methanomicrobiales). Values are shown in the different sediment layers for three control (open black circles) and three bioturbated (open red circles) cores.

The most abundant phyla in the sediment were Proteobacteria, Chloroflexi, Bacteroidetes, and Actinobacteria (Fig. 4a). The NMDS analysis and PERMANOVA tests showed significant differences in bacterial community structures between depths (Fig. 5b, sediment layer effect, PERMANOVA, F_(1,10) = 13.1, p < 0.005). Indeed, phylum-level analyses showed that the relative abundance of Bacteroidetes in the bacterial community was larger in the 5–10 mm than the 0–5 mm depth intervals, whereas the opposite pattern was observed for Proteobacteria (Fig. 4a). Although the presence of foraminifera did not significantly influence the bacterial community structures (PERMANOVA, foraminifera effect, F_(1,10) = 0.53, p > 0.6), the foraminiferal activity significantly reduced bacterial richness in the top sediment layer (Fig. 4c, ANOVA2, “sediment layer × foraminifera treatment” interaction, F_(1,8) = 6.3, p < 0.05). This effect of foraminifera on bacterial OTU numbers was not detected by the Shannon diversity index considering the relative abundance of each bacterial OTU (ANOVA2, F_(1,8) < 0.9 and p > 0.05 for both foraminifera treatment and “foraminifera treatment × sediment layer” interaction). It is also worth noting that bacterial diversity significantly decreased with depth for both control and bioturbated cores (Fig. 4c, ANOVA2, sediment layer effect, F_(1,8) = 106, and p < 0.0001).

Specific analyses performed on the main sulfate-reducing orders of bacteria (Desulfobacterales, Desulfovibrionales, and Syntrophobacterales) showed that the relative abundances (% of reads) of these three orders within bacterial communities increased with depth (Fig. 5a, ANOVA2, sediment layer effect, F_(1,8) = 54, and p < 0.0001). The relative abundance of sulfate-reducing orders in the 5–10 mm depth interval was significantly different in cores with foraminifera and in control cores (ANOVA2,“sediment layer × foraminifera effect” interaction, F_(1,8) = 6.5, and p < 0.05), such as there being a 20 % reduction of sulfate-reducing prokaryotes in cores with foraminifera.

Figure 6 (a) Archaeal community structure (relative abundance) in control (C) and bioturbated (F) cores in two sediment layers. Community structure is represented at the phylum level. Only phyla representing at least 1 % of the reads in at least one sample are represented. (b) Non-metric multidimensional scaling of the archaeal communities in control (“C” and open black circles) and bioturbated (“F” and open red circles) cores in the 0–5 mm (“H1” labels) and 5–10 mm (“H2” labels) sediment layers. (c) Richness and diversity of archaeal communities in the different sediment layers in three control (open black circles) and three bioturbated (open red circles) cores.

Archaeal communities were dominated by Thaumarchaeota in the 0–5 mm depth layer and by Woesearchaeota in the 5–10 mm depth layer (Fig. 6A). The pattern observed with depth for Thaumarchaeota was due to the genus Candidatus Nitrosopumilus which represented more than 80 % of reads of the archaeal community sampled in the 0–5 mm depth layer, whereas it corresponded to less than 15 % of reads from the 5–10 mm depth layer. Consequently, NMDS and PERMANOVA tests showed a clear influence of sediment depth on the structure of the archaeal community (Fig. 6b, PERMANOVA, F_(1,11) = 38.3, p < 0.005). This effect was likely due to a significant increase in archaeal richness and diversity between sampled sediment layers (ANOVA2, sediment layer effect, F_(1,8) > 100, and p < 0.0001 for archaeal richness and the Shannon diversity index). In comparison, no significant effect of the treatment was detected on archaeal community structure (PERMANOVA, F_(1,11) = 0.1815, p > 0.82), archaeal richness (ANOVA2, foraminifera effect, F_(1,8) = 1.1, p > 0.32), and archaeal diversity (ANOVA2, foraminifera effect, F_(1,8) = 1.6, p > 0.23). Taxa-specific analyses on relative abundances of methanogenic archaea in communities (Methanobacteriales, Methanosarcinales, and Methanomicrobiales) also revealed no significant influence of the presence of foraminifera (ANOVA2, foraminifera effect, F_(1,8) = 1.8, p > 0.21), whereas the proportion of methanogens in communities increased with depth (Fig. 5b, ANOVA2, layer sediment effect, F_(1,8) = 90.1, p < 0.0001).

4 Discussion

4.1 Oxygen and organic matter as main determinants of microbial communities in control sediments

The decreasing dissolved oxygen concentration measured in sediments usually determines the vertical distribution of microbial communities (Fenchel and Finlay, 2008). In the control cores of our experiment, non-metric dimensional scaling (NMDS) results clearly demonstrated that the bacterial and archaeal communities were structured by the sediment depth and the associated oxygen availability in pore water. For example, the archaeal genus Candidatus Nitrosopumilus, involved in the nitrification process, showed a preferential distribution in the 0–5 mm sediment layer because this genus needs oxygen to oxidize ${\mathrm{NH}}_{\mathrm{4}}^{+}$ into ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Walker et al., 2010). In addition, oxygen penetration depth ranged from 2 to 5 mm in control cores, and strictly anaerobic microorganisms like sulfate-reducing bacteria and methanogenic archaea were more represented in the communities found in the anoxic 5–10 mm sediment layer than those in the shallowest sediment layer (0–5 mm).

Without organic matter (OM) addition during the experiment, we also observed in the control cores that the total organic carbon (TOC) content was slightly lower in the upper sediment layer than in the deep layer, likely due to the positive influence of oxygen availability on the mineralization of OM in sediments. Indeed, the aerobic mineralization of sedimentary OM is known to be faster than anaerobic mineralization, irrespective of the degree of lability of OM (Kristensen et al., 1995). The vertical distribution of dissolved oxygen in sediments was thus governing OM dynamics and the structure of microbial communities. In turn, the vertical gradient of TOC and TN in sediments generated by OM mineralization could also shape the bacterial community. For example, the lower representation of phylum Bacteroidetes – which are abundant in nutrient-rich aquatic environments (Landa et al., 2013) – in the top sediment layer compared with the bottom layer could be due to the low OM measured in the control cores at the end the experiment.

Overall, in the control sediment, both oxygen and OM availability were the main parameters structuring microbial communities in the present experiment. In such conditions, we can expect that if foraminiferal activities modify these two determinants they would in turn modulate the microbial compartment.

Finally, we observed fluctuations in the oxygen penetration depth (OPD) and diffusive oxygen uptake (DOU) in control cores at the beginning of the experiment (from 0 to 22 d), which might be due to insufficient acclimation prior to the experiment that did not allow it to reach steady-state oxygen microdistribution in the sediment. In the second part of the experiment (after 22 d), we measured an increase of OPD and a decrease of DOU. Similar observations were made previously in sediment without meiofauna between 5 and 14 d of the experiment (Bonaglia et al., 2020). Although the kinetics are different (likely due to the OM-rich sediment used in their experiments), we may hypothesize that a decrease of available OM throughout the experiment led to non-linear changes in OPD and DOU in the control cores.

4.2 Foraminiferal motion activity

In our experiment, benthic foraminifera built up burrows down to 7 mm in the sediment. Although these burrows were not as deep as the centimetre-long burrows previously reported on miliolid and some deep-sea species (Severin et al., 1982; Groß, 2002; Heinz and Geslin, 2012), they were of the same order of magnitude as is known for the coastal species Ammonia beccarii (Green and Chandler, 1994) and the dominant species in our study Haynesina germanica (Deldicq et al., 2023). These shallow burrows confirm that the intertidal foraminiferal species used in the present experiment prefer oxygenated microhabitats (Bouchet et al., 2009; Cesbron et al., 2016). However, foraminifera could burrow 2 mm deeper than the maximal oxygen penetration depth measured in the experimental cores. Although foraminiferal mobility is known to be inhibited by low oxygen concentration (Maire et al., 2016), it seems that during our experiment the community dominated by H. germanica remained active even below the oxygen penetration depth, suggesting that their burrows might provide enough dissolved oxygen to sustain their activity.

Despite this tolerance to low oxygen concentration, observations showed that foraminifera mainly created their burrows during the first 3 weeks of the experiment, and no new burrow could be observed during the period lasting from 40 to 90 d (observations being made from the cores edges, it is possible that some burrows inside the cores were not visible on the core walls). This contrasts with previous reports suggesting that frequently fed deep-sea foraminifera can continuously generate new burrows over the course of several years (Hemleben and Kitazato, 1995). The difference could come from behavioural differences between deep-sea foraminifera and the coastal species used in our experiment or due to the lack of added food in our setup, which might have starved the foraminifera and hence limited their long-term activity.

Despite this potential limitation of foraminiferal activity by fresh OM, the TOC content measured at the end of the experiment in sediments (from 1.4 % to 1.7 %) was of the same order of magnitude as contents usually reported from sediments of the Authie Bay and Boulogne-Sur-Mer harbour (ranging from 1 % to 1.7 %; Francescangeli et al., 2020). Although their reduced activity at the end of the experiment may likely be due to the absence of fresh OM input, foraminiferal survival remained high with on average 75 % of the individuals found alive after 85 d of experiment, stressing that the experimental conditions were close to those observed in the field.

4.3 Foraminiferal bioturbation stimulates aerobic organic matter mineralization

Foraminiferal activity in the first month of the experiment resulted in a significant increase of OPD with a maximum average difference of about 0.7 mm between the bioturbated and control cores on day 36 (average difference was about 0.3 mm in the 9–55 d time interval). It therefore suggests that benthic foraminiferal burrowing activity increased the volume of oxygenated sediment by about 20 %, which is of the same order of magnitude as previously reported in other meiofaunal organisms (Bonaglia et al., 2020). In both foraminifera (this study, day 36) and meiofauna (Bonaglia et al., 2014, 2020), the OPD enhancement led to a decrease of DOU in bioturbated cores, suggesting that foraminifera affect dissolved fluxes in a similar way as meiofaunal ostracods, nematodes, copepods, and oligochaetes.

Nevertheless, macro-invertebrates and meiofaunal organisms seem to have different impacts on benthic oxygen fluxes. Bioturbating macro-invertebrates tend to increase both the DOU (Forster and Graf, 1995; Volkenborn et al., 2007; Lagauzère et al., 2009) and the total oxygen uptake (TOU; Kristensen, 1985; Pelegrí and Blackburn, 1994; Michaud et al., 2005; Politi et al., 2021). In contrast, meiofaunal bioturbation leads to a decrease in DOU (this study; Bonaglia et al., 2014, 2020) and an increase in TOU (Bonaglia et al., 2014). In the freshwater environment, bioirrigation by chironomid larvae increased DOU in organic-matter-poor sediment, whereas the same bioturbation activity decreased DOU in organic-matter-rich sediment (Stief and de Beer, 2002), suggesting that OM availability and benthic microbe respiration mitigate the effect of bioturbators on diffusive oxygen fluxes. In our experiment, the decrease of TOC in cores with foraminifera suggests an increase in OM mineralization. Hence, the decrease in DOU would likely be a consequence of the reduced OM availability in bioturbated cores.

In previous work, the reduced DOU was interpreted as an increase of meiofaunal predatory pressure on their bacterial preys, leading to a decrease in the population of aerobic prokaryotes (Bonaglia et al., 2014). In our study, bacterial richness was positively correlated to TOC, suggesting that the low bacterial richness in sediment layers bioturbated by foraminifera was due to low OM content rather than a top-down control by predation. A similar mechanism was described in freshwater sediments with tubificid worms which reduced the quantity and the quality of the sedimentary OM by stimulating OM mineralization, leading, in turn, to a decrease in bacterial richness and diversity (Cariou et al., 2021).

As the availability of fresh OM had significant control over bacterial community structures in marine sediments (Deng et al., 2020), foraminifera most likely reduced the quality (consuming the most labile fraction of OM) and the quantity of the OM in sediments by stimulating OM mineralization (i.e. total organic carbon loss) during the 3 months of the present experiment. Consequently, the availability and quality of OM were more limiting in bioturbated than in non-bioturbated sediments, hence reducing the ability of multiple bacterial taxa to coexist (increased competition with the reduction of trophic niches; Langenheder et al., 2010; Šimek et al., 2014). Such reduction of the number of trophic niches available in the sedimentary column would have then decreased the bacterial richness. Nevertheless, this effect was not observed by the Shannon bacterial diversity index because the reduction of OM associated with foraminifera activities probably affected low-abundance (rare) OTUs, which have a lower influence on the Shannon diversity index than on bacterial richness (e.g. Haegeman et al., 2013). It is also worth noting that the collection of samples for microbial communities was done after 85 d of the experiment when the effect of foraminifera on dissolved oxygen gradient was not significant. In these conditions, we can expect that microbial changes were less marked at this date than after 1 month of the experiment when foraminifera had the strongest effect on oxygen concentrations in sediments. As already mentioned for sulfate-reducing bacteria and methanogenic archaea, the availability of dissolved oxygen was recognized as a main structuring factor of microbial community structure and biogeochemical processes in marine sediments (Kristensen and Holmer, 2001; Bertics and Ziebis, 2009). Thus, future experiments should measure the dynamics of microbial communities during experiments to evaluate the potential time-dependent effects of foraminiferal bioturbation on the microbial compartment.

4.4 Foraminifera modulate anaerobic diagenetic processes

In our study, benthic foraminifera improved the pore-water oxygenation, and their burrows might also affect a series of diagenetic processes. Indeed, coastal foraminifera are known to accumulate large amounts of nitrate in their cells (Geslin et al., 2014; Langlet et al., 2014; LeKieffre et al., 2022), and deep-sea foraminifera can reduce nitrate and greatly contribute to benthic denitrification (Langlet et al., 2020; Choquel et al., 2021). Our results suggest that foraminiferal bioturbation also affected the benthic nitrogen cycle via enhancing microbial OM degradation, since lower total nitrogen (TN) content was measured in sediments bioturbated by foraminifera in comparison with control sediments. Similar decreases in TN have been reported in sediments bioturbated by macro-invertebrates (Shen et al., 2017; Cariou et al., 2021). Several bioturbating meiofaunal organisms (including rotifers, polychaetes, oligochaetes, crustaceans, ciliates, and nematodes) were also shown to affect the benthic nitrogen cycle by enhancing microbial denitrification (Rysgaard et al., 2000; Prast et al., 2007; Bonaglia et al., 2014). Although not quantified in this experiment, we can expect that foraminiferal bioturbation might affect microbial denitrification in a similar way as other meiofaunal organisms. Thus, further experiments using ¹⁵N-nitrate tracing methods (Bonaglia et al., 2019) will be necessary to determine whether foraminifera contribute to the benthic nitrogen cycle via enhancing the denitrifying activity of microorganisms by bioturbation.

Furthermore, the enhancement of oxygen penetration depth by meiofaunal bioturbation can accelerate sulfide removal (Bonaglia et al., 2020). Ventilation of ghost shrimp burrows was also reported to increase sulfate reduction in reduced micro-niches (Bertics and Ziebis, 2010). In addition, bioturbation can control the community composition of sulfate-reducing bacteria (as shown in meiofauna; Bonaglia et al., 2020) and the abundance of active sulfate-reducing bacteria (as shown in macro-invertebrates; Mermillod-Blondin et al., 2004). In our experiment, we observed a low relative abundance of Desulfobacterales, Desulfovibrionales, and Syntrophobacterales OTUs in the deepest layer (5–10 mm) of bioturbated cores. Although these orders have a flexible metabolism (Dörries et al., 2016), they are generally considered to be a good proxy for anaerobic sulfate reduction (Wasmund et al., 2017), suggesting that foraminiferal bioirrigation might inhibit sulfate reduction in the sediment. Foraminifera are known to be sensitive to free sulfide (Bouchet et al., 2007; Richirt et al., 2020), so the oxygenation of their burrows likely provides a sulfide-free microhabitat in deeper sediment layers.

Finally, our analysis on the proportion of methanogenic archaeal groups in the community did not support the hypothesis that foraminiferal bioturbation activity influenced methanogenic processes in sediments. This corroborates previous experiments showing no effect of bioturbating meiofauna on methane fluxes (Bonaglia et al., 2014). Methanogenesis usually occurs in deeper sediment layers in organic-matter-rich sediments (Froelich et al., 1979). Methane production is likely minimal in the top centimetre of the sediments used in the present experiment, as indicated by the low relative abundance of methanogenic archaea (< 2 % of all the archaea population). Further experiments using deep-dwelling foraminiferal species and organic-matter-rich sediment would be of great interest to evaluate the potential role of these organisms in the benthic methane cycle.

4.5 Foraminifera as ecosystem engineers

Our results clearly show that foraminifera, at densities commonly reported in coastal environments, affect oxygen distribution and fluxes in the sediment via their burrowing activity. Previous studies showed that foraminifera rework sediment (Groß, 2002; Deldicq et al., 2021), and the present study takes our knowledge a step further in showing that they can also perform bioirrigation; hence, foraminifera should now be considered to be bioturbators. We also report that foraminifera affect prokaryotic distribution and diversity, showing a broad impact on the benthic ecosystem functioning and suggesting that foraminifera might be single-celled ecosystem engineers (as defined by Jones et al. (1994): ecosystem engineers “directly or indirectly modulate the availability of resources to other species, by causing physical state changes in biotic or abiotic materials. In so doing they modify, maintain and/or create habitats”).

The foraminiferal vertical distribution pattern is known to be affected by macrofaunal bioturbation (Bouchet et al., 2009; Thibault de Chanvalon et al., 2015; Maire et al., 2016), and meiofaunal bioturbation processes are deeply interconnected with macrofaunal organisms (Nascimento et al., 2012; Bonaglia et al., 2014; Lacoste et al., 2018). To fully discuss the role of foraminiferal bioturbation in benthic ecosystem functioning, it now appears necessary to further study their interactions with other benthic compartments such as meiofauna and macrofauna.