Ocean acidification and warming modify stimulatory benthos effects on sediment functioning: An experimental study on two ecosystem engineers

Many macrofauna have a stimulatory effect on sediment functioning through their burrowing, feeding and irrigation activities. Here, we investigated the single and combined effect of ocean acidification and warming on the stimulatory effect of two key-species inhabiting sandy seabeds in the Southern Bight of the North Sea; the bivalve Abra alba and the polychaete Lanice conchilega. The species were separately incubated in natural sediment in the laboratory under ambient, low pH (pH: -0.3), warm (T: + 3°C) and mimicked climate change (pH: -0.3, T: +3°C) conditions. After six weeks of incubation, nutrient and oxygen exchange were measured at the sediment-water interface to estimate aerobic sediment metabolism and nitrogen cycling. Both species facilitate sediment community oxygen consumption, nitrification and denitrification under ambient conditions. The stimulatory effect of A. alba disappeared in a low pH environment and decreased over time in the warmer treatments along with increased mortality. In contrast, L. conchilega stimulated sediment biogeochemical cycling more when seawater becomes acidified (+ 8 to 41%, depending on the function) but warming had no effect. We explain these species-specific climate change effects by different behavioral and physiological coping strategies that cascade on to sediment biogeochemical cycling, especially through altered oxygenation the sediment matrix.


Introduction
The important role of macrofauna species as benthic ecosystem engineers has long been recognized. Infaunal macrobenthic species are known to alter both the physical structure of and the chemical fluxes within sediments (Laverock et al., 2011). The foraging movements of macrobenthic species rework the top layer of the sediment [i.e. bioturbation (Meysman et al., 2006;Kristensen et al., 2012)], facilitating transport of oxygen and organic matter into deeper layers and of excretion products in and out of the sediment towards the water column. This altered transport of solutes and solids that is driven by physiological and behavioral processes stimulate microbial activity and organic matter mineralization and thus increase turnover of nutrients (Mermillod-Blondin et al., 2004).
Over the last decades, human activities have strongly increased the atmospheric carbon dioxide (CO 2 ) levels. If no drastic CO 2 mitigation measures are taken, this increase is predicted to continue up to atmospheric CO 2 concentrations of ca. 940 ppm by the end of this century compared to a current value of ca. 416 ppm (NOAA; Collins et al., 2013). Oceans are counteracting this increase by absorbing CO 2 at a current rate of 10 6 tons per hour (Brewer, 2009). This absorbing capacity, however, has already led to a decrease in surface-ocean pH of 0.1 units since the start of the industrial revolution (Caldeira and Wickett, 2003). A further decline of about 0.29 pH units by the end of the 21 st century is predicted for the open ocean (Collins et al., 2013;Bindoff et al., 2019;Allan, 2021). In coastal shallow waters, acidification rates often follow different trends. These areas are subject to other stressors, such as freshwater inputs (Carstensen and Duarte, 2019), eutrophication (Cai et al., 2017) and seasonal variations (Hartman et al., 2019;Omar et al., 2019) which may lead to higher acidification rates . On top of their effect on ocean pH, higher concentrations of atmospheric CO 2 and methane enhance the greenhouse effect, thus impacting ocean surface temperature. This is expected to lead to an increase of global (0 -200 m depth) ocean temperature of 1.0 up to 3.2°C by the end of this century (Collins et al., 2013). These alterations in both seawater temperature and pH may challenge marine benthic species and the way in which they affect their environment (Fabry et al., 2008;Byrne, 2011;Dupont and Pörtner, 2013;Nagelkerken and Munday, 2016). Ocean acidification (OA) disrupts acid-base regulations of marine species through the intake of acidified water and calcifying organisms in particular are negatively affected as OA causes the dissolution of their CaCO 3 shells (Matoo et al., 2013;Thomsen et al., 2015;Lemasson et al., 2017;Hurd et al., 2019;Figuerola et al., 2021), whereas changing temperatures may move organisms away from their thermal optima (Pörtner et al., 2004;Sokolova et al., 2012;Reddin et al., 2022). The consecutive behavioral and physiological acclimation and adaptation to climate change is well documented for macrobenthic species (e.g. Tuomainen and Candolin, 2011;Ong et al., 2017;Vlaminck et al., 2022a;Vlaminck et al., 2022b). Alma et al. (2020) for example, found reduced shell strength when the purple-hinge rock scallop was exposed to a pCO 2 of 1050 μatm, while Steeves et al. (2018) predicted increased growth rates of two bivalve species; Crassostrea virginica and Mytilus edulis, in a warmer climate. Van Colen et al. (2020) identified a shift in feeding mode of Scrobicularia plana from predominantly suspension feeding to deposit feeding, in a warmer and acidified environment. Little research has investigated the cascading effect on ecosystem functioning. For example, previous research has demonstrated that OW increased and OA decreased microbial gene abundance and thus modify structure and functioning of nitrogen cycling microbial communities (Currie et al., 2017;Seidel et al., 2022). However, it remains largely unexplored to what extent climate change induced alteration in macrofauna behavior that shape the habitat conditions for the microbial community, i.e. ecosystem engineering, also affect functioning.
In this study, we aim to improve our understanding of the effects of climate change on single species contributions to ecosystem functioning, with a focus on sediment nitrogen cycling. We studied the effect of climate change on two widely spread key species that are abundant in the Southern Bight of the North Sea: the white furrow shell Abra alba and the sand mason, Lanice conchilega (Tebble, 1976;Hartmann-Schröder, 1996). Abra alba is a tellenid bivalve living in the upper layer of the sediment. In search for food, A. alba randomly reworks the upper layer of the sediment, impacting organic matter sequestration and resuspension (Bernard et al., 2012). Lanice conchilega is a tube-dwelling polychaete that can occur in high densities, forming reef-like structures that increase local biodiversity (Rabaut et al., 2007;De Smet et al., 2015). Lanice conchilega acts as a piston pumper that ventilates the sediment, oxygenating deeper layers and stimulating nitrogen cycling (Forster et al., 1999;Van Hoey et al., 2008). Recently the effect of ocean acidification on the behavior and ecophysiology of these species has been investigated. Vlaminck et al. (2022b) found a decreased energy intake and a higher metabolic loss through increased excretion rates for the bivalve. Another study shows that this decreased energy intake, amongst others resulted from a reduced feeding activity (Vlaminck et al., 2022a). Additionally, L. conchilega increased pumping frequency when pH was lowered (Vlaminck et al., 2022a). Whether this behavioral and physiological plasticity to OA, and possible warming-related acclimation, also translates into altered sediment functioning has hitherto not been investigated.
In the present work, two microcosm experiments were conducted, where sediments and Abra alba and L. conchilega individuals were incubated under four different scenarios; a control environment with ambient temperature and pH (CTRL), a scenario with elevated temperature (OW: ocean warming), a scenario with decreased pH (OA: ocean acidification) and a scenario with both an elevated temperature and a decreased pH (CC: climate change). After 6 weeks, nutrient and oxygen fluxes at the sediment-water interface were measured to investigate the effect of A. alba and L. conchilega on the benthic ecosystem metabolism and biogeochemical cycling under the four scenarios. When pH is lowered, we expect the organisms to cope with hypercapnia according to the previously observed behavioral and physiological changes (Vlaminck et al., 2022a;Vlaminck et al., 2022b). More specifically, as A. alba is expected to lower respiratory activity and feeding under low pH, we hypothesize that this species' contribution to sediment biogeochemistry will decrease. In contrast, as L. conchilega is expected to increase its pumping frequency when exposed to low pH seawater, we hypothesize that this species' contribution to sediment ecosystem functioning will increase. In addition, ocean warming may alter the contribution of both species, depending on their physiological responses to temperature variability and the warming effects on the metabolic rate of the sediment microbial community.

Field sampling
Sediment for all experiments was collected in June 2019 at station 780 (N:51°28.292', E:3°3.523') in the Belgian part of the North Sea; this station is characterized by fine sandy sediments with an organic carbon content of 0.40% and nitrogen content of 0.04-0.07% (Braeckman et al., 2014). Multiple Van Veen grabs were deployed from the Research Vessel Simon Stevin. On board of the ship, sediment was first sieved over a 1-mm mesh to eliminate all macrofauna. Subsequently, the sediment was thoroughly mixed before placing it in the experimental set-up. Abra alba individuals were collected at the same station and with the same gear on August 27 th 2019. Lanice conchilega individuals were collected at an intertidal sandy beach in Boulogne-sur-Mer (Nord Pas-de-Calais, France: N: 50.7345, E:1.5881) on July 4 th of the same year; we decided to collect the latter species from intertidal sediments since Van Veen grabs damage these organisms and we wanted to ensure collection of complete, undamaged and similarly sized tubes. Lanice conchilega specimens were transported to the laboratory in buckets with sediment and aerated seawater within 4 h after collection, and were incubated in circulating aerated seawater at in situ conditions in a climate-controlled room (Table 1). Water flow (approximately 500 L.h -1 ) was established between two tanks; one incubation tank (120 L) holding the cores with sediment and macrofauna (see 2.2), and one storage tank (80 L) with the pump. On a weekly basis, 30 L of seawater was refreshed to maintain ambient salinity and avoid nutrient build-up. All systems were subjected to a 12:12 h light:dark regime in a temperature-controlled room.

Experimental set-up and seawater manipulation
Separate experiments were conducted with A. alba (from 01/07/ 2019 to 12/08/2019) and with L. conchilega (from 02/09/2019 to 18/ 10/2019). 32 Plexiglas cores (10 cm internal diameter) were used in each experiment in order to establish a total of 4 replicates for the 4 seawater condition treatments for each of the two sediment types (with and without macrofauna). All cores were filled with sieved and mixed sediment up to approximately 10 cm and distributed over 16 tank systems, 4 replicate tanks for each of the 4 treatments. Sediment (volume = 0.772 ± 0.006 L) was allowed to settle and War 3.9 ± 0.2 2.5 ± 0.4 4.5 ± 0.3 2.9 ± 0.4 Wca 5.7 ± 0.3 3.6 ± 0.5 6.6 ± 0.4 4.1 ± 0.6 Standard deviation is given for each value. All tanks were kept at control values before manipulation. stabilize for 7 days, after which A. alba or L. conchilega individuals were added at natural densities (A. alba: n = 9, L. conchilega: n = 6; Van Hoey et al., 2014) to one core in each tank ). A fully crossed experimental design with the factors temperature (ambient or elevated), pH (ambient or lowered) and exposure time (3 or 6 weeks), was followed to investigate the effects of climate change on macrofauna and their effect on sediment biogeochemistry. After one (L. conchilega) and two (A. alba) weeks of acclimatization to ambient conditions, seawater pH was lowered by 0.1 unit per day and temperature was elevated by 1°C per day over a period of three days following a commonly applied gradual change for continental shelf benthos from temperate regions (e.g. Needham 2007, Matoo et al., 2013;Voet et al., 2022). Seawater pH manipulation (-0.3 units) was achieved using an IKS Aquastar Industrial control system with pH glass electrodes through the bubbling of pure CO 2 . PH electrodes were 2-point calibrated on a weekly basis with Hanna instrument buffer solutions (pH 4.01 and 7.01). The climate room was set at ambient temperature and an external heating system (Teco Refrigeration technologies; Model: TK2000) was used to manipulate seawater temperature (+3°C). Temperature and pH were logged every 15 minutes with the IKS Aquastar system. Twice a week, organisms were fed 5 mL of commercial Shellfish Diet 1800, corresponding roughly to 125 mgC m -3 (Reed Mariculture Inc., composed of 40% Isochrysis, 15% Pavlova, 25% Tetraselmis and 20% Thalassiosira weissflogii) diluted in 16 L of seawater and distributed equally in the 16 tanks. On a weekly basis, salinity was measured with a portable conductivity meter (YSI 30M salinity system). Temperature and pH were measured every week with a portable pH/conductivity meter (Metrohm; Model: 914). This pH electrode was calibrated using a TRIS buffer to allow conversion of the continuously measured pH (NBS) data to total (pH) scale (Dickson, 2010). Weekly water samples were taken from each tank and filtered through GF/C filter papers for total alkalinity (TA) measurements. The samples were kept at 4°C until analysis via titration using a HydroFIA TA (CONTROS Systems & Solutions GmbH) at the Flanders Marine Institute (VLIZ). Carbonate chemistry parameters (pH T , pCO 2 , HCO − 3 , CO 2− 3 , Wc and Wa) were derived from the measured values for temperature, pH, salinity and TA via CO2SYS software (Pierrot et al., 2011) (Table 1). Mortality was checked daily and dead organisms (L. conchilega: n = 2, A. alba: n = 61) were removed to avoid nutrient build-up.

Flux measurements and sampling
After three and six weeks of incubation, the sediment community oxygen consumption (SCOC) and nutrient exchange at the sediment water interface were quantified, using closed dark incubations of the cores. Sediment cores were closed airtight, cutting them off from the circulation system and allowing us to measure changes in oxygen and nutrient concentrations in the water column. The overlaying sea water was mixed to avoid stratification at a rate low enough to avoid resuspension of the sediment. Oxygen concentration in the overlaying sea water was continuously monitored with Robust Oxygen Optodes (OXROB 10, Pyroscience sensor technology). An optical oxygen meter (FirestingO 2 , Pyroscience) was used to control the sensors. Incubations were run for 3 h, sufficient to measure changes in oxygen and nutrient concentrations while the oxygen concentration of the seawater remained above 60% saturation. Water samples (10 mL) were taken from the total volume of water (0.728 ± 0.006 L) in the cores, and filtered through Whatman GF/F filters, at the start and end of the closed incubation, to determine nutrient concentrations.
Oxygen fluxes were calculated from the regression slopes of the oxygen concentration over time, according to Eq. (1).
where dC dt is the change in oxygen concentration in the overlying water (in mmol L −1 d −1 ), V is the volume of the overlying water (in L), and A is the sediment surface area (in m 2 ).
Nutrient fluxes were calculated as the difference between the nutrient concentrations at the end and start of the incubation. After the closed-core incubations after 3 weeks of incubation, the cores were reconnected to the normal water circulation, in order to be sampled again after 6 weeks of incubation.
Oxygen Penetration Depth (OPD) were measured in the sediment after 5 weeks of incubation. Oxygen microsensors (Unisense OX100, 100 μm tip size) were vertically inserted into the sediment and moved with a Unisense microprofiler at increments of 250 μm (two replicate profiles in each core). SensorTracePro software (Unisense) was used to visualize the output. The O 2 microsensor was 2-point calibrated with a 0% oxygen solution (Na-ascorbate solution) and a 100% oxygen solution (obtained through oxygenation of seawater through bubbling).

Mass budget calculations
Nitrification, denitrification and total mineralization were estimated using the sediment-water exchange fluxes of O 2 , NO x and NH x . A mass budget of oxygen, nitrate and ammonium in the sediment was constructed as a function of source and sink processes (Soetaert et al., 2001;Braeckman et al., 2010). Oxygen, nitrate and ammonium are exchanged across the sediment-water interface. Oxygen is consumed through the oxidation of organic carbon (oxic mineralization) and of ammonium to nitrate (nitrification), or through the re-oxidation of reduced substances formed through anoxic mineralization. The O:C ratio is assumed to equal 1 (Redfield), while 2 moles of oxygen are required for the oxidation of ammonium to nitrate (Soetaert et al., 2001). Ammonium is the end product of organic nitrogen mineralization and is consumed during denitrification. Nitrification produces nitrate but 0.8 moles of NO 3 are being consumed for each mole of carbon that is denitrified. Part of the reduced substances remain buried in the sediment and are therefore not available for re-oxidation (Soetaert et al., 2001).
The rates of change are set to zero, based on the assumption that a geochemical steady state is reached approximately 1 week after the addition of organisms. We assume that the burial of anoxic substances is negligible, allowing us to combine oxic and anoxic mineralization (Braeckman et al., 2010). To solve the model, an extra equation is added, using the N:C ratio (mean value = 0.1351 mol N mol C -1 ) to appoint the relationship between nitrogen and carbon mineralization, resulting in a set of equations (that were solved using the limSolve package (Soetaert and Van den Meersche, 2017) in the open source software R version 4.0.4 (R Core Team, 2014).

Statistical analyses
Differences in survival of study organisms between treatments at the end of incubation were analyzed using two-way analysis of variances (ANOVA) with temperature and pH as fixed factors, with two levels each (ambient and increased temperature; ambient and lowered pH). Normal distribution and equal variances were assumed and tested with Shapiro and Levene's test respectively. If significant effects were found, Tukey HSD tests were performed to identify pairwise differences. To test the effect of A. alba or L. conchilega (present or absent), temperature, pH, and time (three or six weeks of incubation) and their interaction on biogeochemical fluxes, linear mixed effects models were constructed for each factor. We allowed random intercepts for tank identity to take the nonindependence of the two cores in each tank into account.
Step-bystep deletion of single model terms that were not involved in any significant interaction led to the minimal adequate model. The Kenward-Roger method was used to estimate degrees of freedom and calculate F-statistics and associated p-values. When no warming x pH interaction effects are present, the model compares the combined acidified treatments (OA and CC) with the nonacidified treatments (OW and CTRL) combined and/or the warmed (OW and CC) with the non-warmed (CTRL and OA) treatments. Model assumptions were checked graphically; the normality of residuals was numerically confirmed with a Shapiro-Wilk test. Statistical analysis was conducted in R-software version 4.0.4 (R Core Team, 2014) with packages lme4 (Bates et al., 2014), lmerTest (Kuznetsova et al., 2017) and MuMIn packages (Barton, 2020).

Nutrient and oxygen fluxes and cycling 3.2.1 Lanice conchilega
The presence of L. conchilega significantly increased OPD by 31 (± 1) %, compared to cores with no L. conchilega individuals present ( Figure 2). Furthermore, L. conchilega increased SCOC by 28% under ambient conditions, whereas stimulation in an acidified environment was even greater with a 37 (± 8) % increase ( Table 3,  Effects of temperature and pH on the survival rate of A. alba and L. conchilega during the 6-week incubation period. Error bars represent standard errors (± SE, n = 36). Vlaminck et al. 10.3389/fmars.2023.1101972 Frontiers in Marine Science frontiersin.org Temperature and pH effect on the oxygen penetration depth in the presence (orange circles) and absence (grey triangles) of L. conchilega and A. alba. Error bars represent the standard errors (± SE, n = 4).  Effect of Lanice conchilega presence (present: orange bars, absent: grey bars) on sediment community oxygen consumption (SCOC, mmol O 2 m -2 d -1 ), nitrification, and denitrification rates (mmol N m -2 d -1 ), in 4 different treatments, with different times pooled because experiment duration did not affect the fluxes and processes. Error bars represent standard errors (± SE). %, when comparing all non-acidified treatments with all acidified treatments ( Table 3, interaction effect of OA and organism; Figure 3; pairwise p-values< 0.001). On average the highest stimulation was found in the CC treatment, although this effect was not significant (p-value = 0.842). Experimental duration did not affect fluxes, therefore time was pooled in the graphic representation of the results (Figure 3).

Abra alba
Experimental duration significantly affected nitrification and denitrification rates, with respective decreases of 23 (± 6) % and 22 (± 1) % after six weeks of incubation compared to three weeks of incubation (Table 2, Time effect; Figure 1).
Abra alba presence increased OPD by 27 (± 6) % in all treatments, when compared to all cores without A. alba presence ( Table 2, organism effect; Figure 2). Irrespective of the presence of A. alba, the interaction of seawater pH and temperature significantly affected OPD ( Table 2, interaction effect of OA and OW; Figure 2). On average, OPD in sediments with Abra was deeper in CTRL when compared to OW, OA or CC (17%, 19%, and 6%, respectively).
In non-acidified conditions, A. alba significantly stimulated SCOC resulting in an increase of 35 (± 15) % averaged over time ( Table 2, interaction effect of OA and organism; Figure 4; pairwise p-value< 0.001). Although the stimulation was on average larger in the CC treatment with A. alba, the statistical comparison revealed that SCOC rates did not differ significantly between the sediment with and without A.alba in the CC treatment (p = 0.086). This stimulatory effect was absent under acidified conditions. Nitrification and denitrification were both stimulated by the presence of A. alba in a control environment (pairwise p-value = 0.001 and 0.010 respectively). This stimulatory effect was absent when either seawater temperature was raised (OW) or pH was lowered (OA) (pairwise p-values all< 0.05). However, nitrification and denitrification rates were stimulated by the presence of A. alba under combined warming and acidification, i.e. the CC scenario, with an increase of 61 (± 5) % and 110 (± 1) % averaged over time respectively (pairwise p-values: 0.008 and 0.010 respectively; Table 2: OW, OA and organism interaction effect; Figure 4).

Discussion
This experimental study investigated climate change effects on sediment ecosystem functioning, in particular how macrofauna bioturbation affect nitrogen cycling under simulated climate change conditions.". Interestingly, the effect of climate change differed between the two investigated key species. Both species facilitate biogeochemical cycling under ambient conditions. Lanice conchilega stimulates sediment aerobic metabolism and Effect of Abra alba presence (present: orange bars, absent: grey bars) on sediment community oxygen consumption (SCOC, mmol O 2 m -2 d -1 ), nitrification and denitrification rates (mmol N m -2 d -1 ), in 4 different treatments at 2 different times. Error bars represent standard errors (± SE).
biogeochemical cycling more when seawater becomes acidified, while warming has no significant effect on the species contribution. Contrary, Abra alba contributes less in an acidified environment and the contribution in the warmer treatments decreases over time along with lower densities. The measured N cycling process rates in this study, using strongly manipulated and reduced communities, are generally slightly lower than the natural rates described by Toussaint et al. (2021) for the same location of sediment collection. Nevertheless, the chosen approach does allow for empirical proof and understanding on how different aspects of climate change, i.e. ocean warming and acidification, modify the contribution of marine benthic fauna to nitrogen cycling.
Laboratory and climate change conditions had no effect on Lanice conchilega survival while Abra alba survival was reduced both by experiment duration and ocean warming. Although both species are widely spread in temperate continental shelf sea beds, the more abundant occurrence of L. conchilega in warmer regions such as the coastlines of Mauritius, Madagascar and Mozambique (Day, 1968;Lewinsohn and Fishelson, 1967) as well as this species' presence in the lower intertidal zone of estuaries and beaches (Alves et al., 2017) that are prone to strong diurnal fluctuations in temperature and dryness confirms that this species is more tolerant to suboptimal laboratory conditions, as well as to increased temperature in general.
Abra alba and Lanice conchilega individuals increased oxygen penetration depth in an environment with ambient pH and temperature which corroborates the found stimulatory effect on biogeochemical cycling of both species under those conditions. Pumping activity of Lanice conchilega increased oxygen penetration by 1.32 (± 0.16) mm which is similar to the 1.2 mm increase found for similar organism densities by Braeckman et al. (2010). Interestingly, even though L. conchilega is known to increase pumping frequency in a similarly acidified environment by 30% (Vlaminck et al., 2022a), oxygen penetration depth did not differ across treatments suggesting that in this experiment seawater conditions did not affect tube structure or ventilation magnitude which could lead to change in oxygenation depth (Kristensen et al., 2012). Such increased pumping frequency allows the worm to renew the O 2 pool faster in order to cope with the higher metabolic cost associated with inhabiting low pH environments (Briffa et al., 2012;Turner et al., 2015). This was for example also observed for the polychaete Alitta virens (Godbold and Solan, 2013) while a similar bio-irrigation response to other environmental stressors, like hypoxia and increased temperature, has been found for different polychaete species (Kristensen, 1983;Forster and Graf, 1995;Ouellette et al., 2004). The pumping of L. conchilega transports nitrate into the sediment and increases the availability of oxygen in deeper layers of the sediment (Forster and Graf, 1995), increasing the surface for coupled nitrification-denitrifcation, stimulating these processes. Foshtomi et al. (2018) furthermore proofs that L. conchilega influences the composition and density of denitrifying communities by changing the availability of nitrate. Higher diversity at depth was observed in L. conchilega patches, resulting in high denitrification rates. Increased irrigation activity under low pH conditions can thus explain the observed increased nitrification and denitrifcation rates under simulated ocean acidification and climate change in our study. Warming did however not significantly affect this species contribution to biogeochemical cycling, most likely due to the abovementioned tolerance to temperature fluctuations.
A. alba did not stimulate nitrogen cycling in ocean acidification or ocean warming conditions as was the case under ambient and climate change conditions. After 6 weeks of incubation, the stimulatory effect was greatest under ambient conditions where sediment oxygenation was, on average, the deepest in comparison to the other treatments. In a similarly acidified environment, A. alba reduced respiration and feeding, especially suspension feeding which is known to create the greatest subsurface advective flows (and hence stimulatory effects via sediment oxygenation) by tellinid bivalves (Volkenborn et al., 2012;Vlaminck et al., 2022a). These behavioral changes reduce the uptake of low-pH seawater, thereby limiting possible physiological costs related to hypercapnia but at the same time also reduce particle mixing and bio-irrigation, which can therefore explain the observed absence of A. alba contribution on all measured functions (SCOC, nitrification and denitrification) under ocean acidification in this study. Interestingly, when combined with warming (i.e. the climate change treatment), A.alba did stimulate nitrogen cycling under reduced pH. This difference can possibly be attributed to the lower costs for calcification with higher temperature (Clarke, 1993;Mancuso et al., 2019) as demonstrated by the higher saturation state of calcite and aragonite in both warmer treatments (Table 1). This may leave more energy for feeding and maneuvering in comparison to the single pH treatment where warming was not present, thereby stimulating biogeochemical cycling. This effect however seems to decrease with time when densities dropped at warmer temperatures, i.e. the ocean warming and climate change treatments. Braeckman et al. (2010) studied the effect of A.alba density on bioturbation and its contribution to sediment biogeochemistry. They found decreasing bioturbation rates and ammonium effluxes with decreasing density of A. alba. After six weeks of incubation only four A. alba individuals remained in the warmer treatments in our study, which is a similar density as the lowest that was used by Braeckman et al. (2010), where bioturbation was 78% lower compared to natural densities. Our study shows that at such reduced densities the stimulatory effect can disappear especially in the warming treatment where biogeochemical cycling is boosted in the absence of A. alba.
In conclusion, we show that ocean warming and acidification change the contribution of macrofauna to benthic ecosystem functioning, i.e. biogeochemical cycling. Declines in macrofauna activity or density, as observed for A. alba, will ultimately lead to major changes in the contribution of benthos to global nitrogen and carbon cycling (Sweetman et al., 2017;Snelgrove et al., 2018). However other effects, such as observed for L. conchilega, where ocean acidification increased the contribution of this species to overall sediment metabolism and biogeochemical cycling, could mitigate such reduced sediment functioning related to climate change. The variable results for both key-species depending on physiological tolerance and behavioral plasticity highlight the difficulty to upscale population studies to the community and ecosystem level (Ehrnsten et al., 2020).

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions
EV and CC contributed to conception and design of the study. EV acquired the data and performed the statistical analysis with the input of UB and CC. EV, CC, UB and TM performed the interpretation of the data. EV wrote the first draft of the manuscript. All authors contributed to the article and approved the submitted version.

Funding
Funding for the research leading to this publication was provided through the BRAIN-be project BR/175/A1/PERSUADE of the Belgian Science Policy Office. Additional funding was obtained from the Special Research Fund (BOF) of Ghent University through GOA project 01G02617. Experiments used infrastructure funded by the Flemish Science Fund FWO through project GOH3817N in the framework of EMBRC-Belgium.