What global biogeochemical consequences will marine animal–sediment interactions have during climate change?

The activities of macrobenthos (animals > 500 µm that live in or on sediments) significantly influence the storage and recycling of organic carbon (OC) and biogeochemically coupled elements in sedimentary deposits (Middelburg, 2018). Benthic ecology has a long history of studies on seafloor communities and animal–sediment relationships, from Petersen’s (1913) initial efforts to understand variation in fisheries to Sanders’s (1958) classic paper on feeding modes and sediments and many hundreds of related studies that have attempted to understand the complex relationships between biota and their sedimentary habitats (e.g., Gray, 1974; Rhoads, 1974; Snelgrove and Butman, 1994). More recently, researchers have begun to consider such relationships in the context of climate change (e.g., Smith et al., 2008; Snelgrove et al., 2018; Ehrnsten et al., 2020; Moore and Smale, 2020). Although several studies have documented benthic community response(s) to changing environmental conditions (e.g., warming, acidification, deoxygenation, food availability) in multiple seabed environments (Smith et al., 2008; McClain et al., 2012; Jessen et al., 2017; Sweetman et al., 2017), hereafter referred to as climate change, the implications of these responses for global-scale elemental cycling and OC storage remain underexplored in the study of seafloor sediments, Earth’s largest biome. Our article focuses heavily on sedimentary environments from the land–sea interface onto the continental shelf, where the majority of active carbon cycling with seafloor macrobiota occurs. However, we acknowledge the contributions of deep-sea environments that exchange less carbon on the short term but contribute significantly over long timescales through the vast seafloor area they encompass.

Expected functional responses to climate change by macrobenthic species include poleward distributional shifts (Morley et al., 2018), rearrangements of their communities (taxonomic and functional groups; Grilo et al., 2011), and alteration of their size distributions or behavior (Figure 1a), in addition to regional changes in the spatial (e.g., nearshore to offshore gradient changes) and temporal scales of these responses. The conceptual basis for how climate change will influence important physical drivers of benthic carbon cycling (e.g., substrate type, seabed geomorphology, and benthic boundary layer conditions; Mayer et al., 2018; Townsend et al., 2018) remains underdeveloped, and the predictive basis for doing so remains largely unresolved.

Informed marine management decisions, partly based on carbon credits associated with different carbon burial “hot spots” (locations in the coastal margin with the highest OC burial rates), climate mitigation strategies (e.g., preservation and enhancement of aquatic carbon sinks—“blue carbon”; Nellemann et al., 2009), and ocean biogeochemical models, all require some understanding of the role that ocean ecosystems play in sedimentary OC accumulation and long-term sequestration (Bianchi et al., 2018). Vegetated coastal habitats, such as seagrass (Kennedy et al., 2010), salt marshes, and mangroves (Atwood et al., 2017; Windham-Myers, 2018; Macreadie et al., 2019), have received considerable attention in recent years for their role in C sequestration and burial. However, C cycling in coastal systems is influenced by both above- and below-ground plant production and/or other land-margin accretionary soil characteristics, which differentiate them from the more extensive soft-sediment marine biomes reviewed here.

Past climate shifts have resulted in many large-scale redistributions and extinctions of marine organisms (Bond and Grasby, 2017), which we can use to predict adaptive responses or extinctions of organisms on long timescales (e.g., millennial to decamillennial) that integrate impacts of multiple stressors. The geological record provides a chronology of change that integrates many factors, which could help extrapolate experimental and mechanistic studies that typically address individual climate drivers over shorter timescales (e.g., decadal to centennial). For example, although causes of mass extinctions vary, most correlate strongly with dramatic increases in atmospheric CO₂ (Bond and Grasby, 2017), with corresponding ocean warming, acidification, and oxygen loss (Penn et al., 2018). Large-scale biodiversity loss, exemplified by increased dominance of simplified forms as demonstrated in shallow-tier trace fossils, decreased burrow size and bioturbation intensity, and nonrandom extinction (Twitchett and Barras, 2004; Buatois and Mángano, 2016; Caswell and Frid, 2017) profoundly influence bioturbation rates and patterns in shallow marine ecosystems (Twitchett and Barras, 2004). Global oceanic anoxic events produce similar responses (Caswell and Frid, 2017), suggesting that even localized, persistent stressors associated with climate change can significantly alter soft-sediment biological processes.

Here, we examine how the demographic and evolutionary responses of macrobenthic fauna to climate change could alter biogeochemical cycling in ocean sediments. Acknowledging the coarseness of estimates for size group contributions to total benthic biomass, we estimate that macrobenthic biomass amounts to about 1 Gt (Bar-On et al., 2018), recognizing the uncertainty and the incomplete inclusion of major groups such as echinoderms that dominate macrobenthic biomass in deep water.

Specifically, we (1) explore the potential importance of benthic macrofauna undergoing range expansion on biogeochemical cycling and carbon storage in margin sediments; (2) use trace-fossil diversity and milestones in animal–substrate interactions over geologic time to examine past large-scale macrofaunal community shifts correlated with environmental changes; (3) discuss recent projections of macrofaunal adaptations and/or distribution shifts in response to climate stressors, highlighting ongoing research needs; and (4) hypothesize how potential phenotypic, epigenetic, and genetic responses of macrobenthic populations to climate change could impact benthic biogeochemical cycles and suggest how large-scale biogeochemical models could incorporate this impact.

Ocean sediments represent Earth’s largest interactive reservoir of OC storage (Middelburg, 2018). Continental shelves and other ocean margin deposits cover approximately 16% of global seabed area but account for >90% of total ocean OC burial (Berner, 1982; Hedges and Keil, 1995; Burdige, 2007; Bianchi et al., 2018). Open ocean and deep-sea bottoms receive most of their organic matter from surface production, but coastal sediments, which also rely on surface production, still receive more OC from land (Hedges and Keil, 1995; Middelburg, 2019). The relative importance of terrestrially derived labile OC in coastal sediments varies widely between passive and active margins (Blair and Aller, 2012), with phytoplankton contributions becoming increasingly dominant seaward on wider margins as the influence of terrestrial sources decreases (Bauer et al., 2013). Balances between OC supply, remineralization, and burial can vary seasonally along land-sea gradients through changes in riverine inputs, oxygen availability, and physical forcing (Blair and Aller, 2012). In some locations, elevated nutrient inputs from land cause eutrophication and enhanced delivery of phytodetritus to the seabed, resulting in bottom water depleted of dissolved oxygen and depauperate benthic communities (Breitburg et al., 2018). Spatial distributions of OC in surface sediments of the global ocean reflect the preponderance of high OC inputs near continental margins, where past studies document significant ventilation rates and/or particle bioturbation rates (Figure 2). Intensified OC burial (Bianchi et al., 2018) occurs on shelf and adjacent regions, along with comparatively tight oceanic benthic–pelagic coupling. Observations from shelf systems are critical for understanding and predicting how changes in physical (temperature, circulation, oxygen, and sediment transport) and biological characteristics (production, metabolic rates, species abundance, composition, diversity, and functional groups) will affect ocean C cycling.

Multiple factors regulate OC stocks, including overall oxygen exposure time (Hedges and Keil, 1995) determined by physical and biological transport processes (e.g., erosion, resuspension, deposition; Boudreau and Jorgensen, 2001), biodeposition by suspension feeders (Herman et al., 1999), sediment texture and organic-mineral associations (permeable or muddy sediments; Mayer et al., 1996), and the OC source (e.g., phytoplankton vs. terrestrial OC; Burdige, 2007, 2006). Projections of how climate change will alter spatial–temporal patterns in benthic community structure, with corresponding influences on OC cycling across gradients in depth and latitude, must consider hot spots of burial and remineralization (Bianchi et al., 2018). Microbes and animals alter oxygen dynamics directly through respiration and indirectly via bioirrigation and particle movement (Aller, 1982, 1994; Kristensen, 2000). Burrows increase surface area available for solute exchange and oxic/anoxic interfaces (Kristensen, 2000) and alter porewater exchange with the overlying water (Figure 1a). Below the redox interface, between oxic and anoxic porewaters, a spectrum of anoxic microbial transformations drive the redox processes coupled to the turnover of OC, but often at lower rates than within the oxic layer because of the overall depletion of labile substrates. Electron acceptors (nitrate and sulfate) occur in high concentrations within the oxic zone, in contrast to high concentrations of electron donors (e.g., ammonium and sulfides) in the underlying reduced zone (Aller, 1994; Middelburg and Levin, 2009). The sediment redoxcline gets shallower with lowering bottom-water oxygen levels during periods of high respiration and low bioturbation, resulting in differences in the breakdown of organic compounds under oxic versus anoxic conditions (Middelburg and Levin, 2009; Birchenough et al., 2012).

Terrestrial landscape mapping using satellite remote sensing methods offers ways to measure global terrestrial primary production and carbon stocks (e.g., Running et al., 2004; Mulder et al., 2011; Dorigo et al., 2017; Fan et al., 2019). Ocean researchers have applied comparable approaches to assess net primary production in surface waters and estimate carbon flux to the seafloor (Lutz et al., 2007). These ocean measurements of carbon flux, however, do not address the spatial complexity of the seafloor substrata, benthic organisms, and habitats, which strongly influence the fate of the carbon reaching the seabed. Indeed, high spatial heterogeneity characterizes the continental margins where the highest net primary production occurs. Unfortunately, satellite remote sensing tools sample ocean environments to depths of tens of meters at best, leaving most of the seafloor beyond the reach of these methods. The currently hidden complex spatial patterns of OC cycling, macrofaunal dynamics, and the interplay of these factors in the benthos, therefore, limits current understanding of carbon cycling and sequestration at the seafloor. Many past studies relied on in situ measurements of animal–sediment characteristics to assess fine-scale biogeochemical processes (Figure 1a). However, scaling-up from fine-scale observations to regional and global scales (Figure 1a–c) is extremely challenging given the lack of seafloor environmental data (e.g., seabed substrata morphology, overlying oceanographic characteristics influencing C delivery) available at an appropriate resolution to capture spatial variability in the benthos (Lecours et al., 2015; Townsend et al., 2018).

Advances in seafloor mapping (Figure 1) provide reasonable proxies of seafloor characteristics (e.g., benthic macrofaunal species distributions, macrofaunal community characteristics, and sediment types; Reiss et al., 2014; Robert et al., 2014; Diesing, 2020), extending point sampling information to broader spatial scales based on remotely sensed environmental variables (Brown et al., 2011). Innovations in vessel-mounted acoustic remote sensing techniques over the past two decades, such as multibeam echosounders, offer promising approaches for benthic ecosystem mapping at comparable spatial resolutions to satellite remote sensing on land (tens of meters to submeter resolution depending on water depth; Figure 1b; Brown et al., 2011; Mayer et al., 2018). High-resolution bathymetric measurements can resolve seafloor topography, and extraction of information based on the strength of acoustic seafloor backscatter signals can support inferences on seafloor substrata and habitat characteristics (Figure 1b; Brown et al., 2011; Lamarche and Lurton, 2018). These remotely sensed acoustic data sets have demonstrated considerable potential to map seabed substrate (e.g., Diesing et al., 2014; Misiuk et al., 2018), single species distributions (e.g., Galparsoro et al., 2009; Brown et al., 2012), macrofaunal assemblage patterns (e.g., Lacharité and Brown, 2019), and seabed landscapes sometimes referred to as “seascape” (e.g., Boström et al., 2011; Shaw et al., 2014) or “benthoscape” maps (Zajac, 2008; Brown et al., 2012). The next step in the evolution of these methods is to develop models to extrapolate biogeochemical processes at the fine/millimeter-to-centimeter scale (Figure 1a) to broader spatial scales (Figure 1b), utilizing spatial data to improve estimates of OC burial. Very few studies have attempted this extrapolation (see Smeaton and Austin, 2019; Gogina et al., 2020), which is limited primarily by the paucity of remotely sensed acoustic data sets, with less than 9% of the ocean floor currently mapped by these modern, remote sensing systems (Mayer et al., 2018).

At global scales, ocean digital elevation models, produced using satellite-gravity models that calibrate the gravity-to-topography ratio using bathymetric soundings (Becker et al., 2009), must rely on considerably coarser data resolution than similar terrestrial models (e.g., Figure 1c). Nonetheless, this modeling offers some capability to assess global spatial patterns in marine sediments in the absence of higher resolution remotely sensed acoustic data. Broad-scale sediment mapping has utilized seafloor digital elevation models with in situ sediment samples to generate coarse-resolution regional and global substrate maps (e.g., Dutkiewicz et al., 2015; Mitchell et al., 2019; Diesing, 2020). Such maps clearly oversimplify spatial patterns and heterogeneity apparent at finer spatial scales but may be useful to determine the fate of particulate organic carbon (POC) as delivered globally (Lutz et al., 2007) once it reaches the seafloor.

Although the temporal and spatial scales of geological events differ from those of the Anthropocene, the fossil record can help to predict future bioturbation response to climate change. The fossil record demonstrates that profound changes in benthic macrofaunal composition repeatedly influenced geochemical cycles and fluctuations of the seafloor (Bond and Grasby, 2017). Evolutionary radiations in the marine environment reveal novel burrow architectures and feeding guilds (Mángano and Buatois, 2014) which, in turn, acted as major players on geobiological feedbacks including biogeochemical cycles or sudden episodes (Boyle et al., 2014). Arguably, the most important of these episodes took place around the Ediacaran-Cambrian transition, encompassing the establishment of novel forms of animal–sediment interactions, complex burrows, and intense and deeper bioturbation (Figure 3; Mángano and Buatois, 2014; Gougeon et al., 2018). These innovations may have resulted in the extinction of the Ediacara biota (Laflamme et al., 2013) and accelerated the demise of previously widespread stromatolites (Walter and Heys, 1985) and other microbially-stabilized substrates (Seilacher, 1999).

Subsequent to the Cambrian Explosion, the Ordovician Radiation saw major increases in the diversity of burrow and boring architectures, with innovations in styles of bioturbation and bioerosion, respectively (Levinton and Bambach, 1975; Buatois et al., 2016a, 2020). Novel and long-lasting strategies (Levinton and Bambach, 1975) of animal–substrate interactions first appeared in nearshore and shelf areas, but later expanded into brackish-water coastal environments and into deep marine environments; the latter represented by the establishment of diverse feeding strategies, including the appearance of trace fossils interpreted as taxa that enhance microbial farming and trapping of microorganisms (Seilacher, 1977; Orr, 2001; Buatois et al., 2009). Modern bioturbators and bioeroders rose to dominance during the Mesozoic Marine Revolution, including decapod crustaceans, irregular echinoderms, bivalves, and various worms (Vermeij, 1987; Buatois et al., 2016b). In particular, the polyphyletic adoption of mantle fusion in bivalves allowed greater control of internal water pressure and stronger burrowing (Stanley, 1968). A Mesozoic revolution in benthic predators may have enhanced taxonomic diversification and abundance of deeper-burrowing bivalves (Vermeij, 1987). The Mesozoic Marine Revolution saw remarkable increases in bioturbation depth and in the intensity and efficiency of biogenic reworking, resulting in increased sediment turnover, increased complexity of exploitation patterns in the infaunal ecospace, and the establishment of complex feedback mechanisms between predation and infaunalization (Thayer, 1979; Buatois et al., 2016b; Rojas et al., 2021).

Dramatic increases in oxygenation at about 200 and 400 Mya resulted in persistent overall ventilation of the shallow ocean (Lu et al., 2018). However, crises that strongly reduced biodiversity and bioturbation punctuated these long-term trends. The hallmark of major extinctions, commonly a sedimentary interval displaying little or no evidence of bioturbation, indicates a marked shutdown of biogenic sediment overturn and associated carbon burial. In the most severe case, the end-Permian mass extinction and reduced biodiversity led to the collapse of the biogenic reworking on the shelf and the reappearance of Ediacaran-style sediments (Figure 3; Buatois and Mangano, 2011; Hofmann et al., 2015). In modern oceans, homogenization of the uppermost centimeters of sediments (i.e., sediment mixed layer) can result in a mottled texture. However, bioturbation patterns in Lower Triassic postextinction deposits indicate a mixed layer strongly affected by biodiversity loss, with poorly mixed sediments at or immediately below the sediment–water interface (Buatois and Mangano, 2011; Hofmann et al., 2015; Luo et al., 2021). The collapse of the mixed layer presumably dramatically changed both the geochemistry of marine sediments and the water column, including OC burial, but effects likely altered shallow-shelf phytoplankton consumers and deposit feeders similarly (Levinton, 1996). Suppressed bioturbation may have created a large OC sink in continental shelves and margins, as well as adjacent nearshore areas (Hofmann et al., 2015). Estimates of biotic recovery times after the end-Permian mass extinctions vary, but recent estimations indicate approximately 4 Mya for full recovery of the sediment-mixed layer (Hofmann et al., 2015), underscoring the long-term effects of mass extinctions.

Analysis of the trace-fossil record (i.e., fossils of animal activity) demonstrates numerous interacting factors that influence the spatial and temporal distribution of bioturbators, including hydrodynamic energy, type and degree of substrate consolidation, oxygen content of bottom and interstitial water, salinity, sedimentation rate, food supply, water turbidity, acidification, temperature, and bathymetry, as well as modification of aspects of the physical environment (e.g., substrate consolidation) through bioturbation. Of these factors, bottom-water oxygen, acidification, warming, and carbon delivery are most relevant during biotic crises. In particular, extinctions involving ocean anoxia and acidification suggest habitat displacement, including latitudinal changes in the aftermath of the end-Permian mass extinction, with high-latitude settings providing refuges for benthic organisms from lower latitudes (Fraiser and Bottjer, 2007; Beatty et al., 2008). In the specific case of the end-Permian mass extinction, elevated atmospheric CO₂, increased temperatures, and latitudinal species shifts clearly interlinked (Beatty et al., 2008).

Various trace-fossil models enable reconstruction of infaunal response to decreased oxygenation (Savrda, 2007). Overall, a reduction of bottom-water oxygen parallels decreases in intensity of bioturbation, diversity of trace fossils, depth of penetration, and burrow size. Hypoxia events today produce similar effects, although with clear exceptions (Levin, 2003). Furthermore, periods of oxygen depletion strongly affected nearshore bioturbators in comparison to deeper water benthos, as illustrated by the Late Devonian mass extinction (Buatois et al., 2013). Traces of deposit feeders, which usually characterize preextinction, slightly deeper offshore settings, dominated Upper Devonian postextinction nearshore sediments. In contrast, the notable absence of vertical burrows produced by a suspension-feeding infauna, typical of wave-dominated nearshore sandy settings, suggests differential selection of trophic types during mass extinctions (Buatois et al., 2013). Thus, the fossil record documents dramatic changes in burrow types, distributions, diversity, and functions; by extension, these changes must have altered rates of bioturbation, oxygenation, and OC cycling in seafloor sediments profoundly (Kristensen, 2000; Sweetman et al., 2017; Middelburg, 2018).

Anticipated impacts of climate change on seafloor communities and biogeochemical cycling include effects from rising temperature, oxygen loss (deoxygenation), and changes to the carbonate system (ocean acidification; Sweetman et al., 2017; Figure 4), in addition to influences of climate change associated with spatial–temporal responses in the overlying marine ecosystem (Figure 5 and Table 1). These changes can alter primary productivity and plankton community composition and thus the size, amount, timing, and transport mode of POC delivery to deeper waters, and ultimately to the benthos (Smith et al., 2008; Boyd et al., 2019). Reductions in food supply associated with reduced POC flux will influence macrobenthic communities directly, with expected declines in faunal abundance, biomass, and diversity. These declines will ultimately lead to reduced ecosystem function, including major changes in the contribution of the benthos to C cycling (Smith et al., 2008; Sweetman et al., 2017; Snelgrove et al. 2018).

Ocean acidification, associated with increasing atmospheric CO₂ concentrations, will also measurably alter benthic communities and biogeochemistry. Lowered aragonite saturation states will negatively impact the survival and recruitment of shell-bearing benthos (Green et al., 2013), whereas long-term exposure to elevated CO₂ can produce disparate effects on bioturbation and bioirrigation by coastal, soft-bodied infauna (Figure 4), which may experience natural exposure to porewaters with low pH (Godbold and Calosi, 2013; Widdicombe et al., 2013). These ocean acidification effects on animals can propagate to microbes that govern nutrient cycling via animal–microbe interactions (Laverock et al., 2013).

The scale of species responses to specific environmental variables and the strength of their association with specific benthic habitat characteristics remain critical gaps in our ability to predict benthic system response to climate change (Renaud et al., 2015). The lack of seafloor morphology, substrate, and benthic habitat maps at scales and resolutions required to resolve habitat suitability and migration spatially (e.g., barriers to dispersal such as benthic geomorphology and ocean circulation patterns) has fundamentally limited understanding of future shifts.

Changing oceanographic conditions that exceed physiological thresholds may reduce biodiversity (Figure 5c); however, functional redundancy in benthic communities may provide a buffer against total loss of key functional groups (Rosenfeld, 2002; Rossi et al., 2009; Belley and Snelgrove, 2016; Frid and Caswell, 2016). Moreover, adaptive and acclimation capacities of benthic species may mean less severe declines than expected (see next section; Rosenfeld, 2002; Frid and Caswell, 2016). Responses of macrobenthic species to environmental change can reflect changes in assemblage diversity; however, we generally lack accurate tools to predict such change. Declines in diversity (Figure 5e) can yield declines in ecosystem function (Figure 5d) (e.g., for bioturbation; Danovaro et al., 2008a, 2008b; Sperling et al., 2016), which may influence biogeochemical cycling. Improved biogeochemical model predictions require better understanding of the relationships between bioturbation, seafloor biogeochemical cycling, and biodiversity-ecosystem functioning relationships (Snelgrove et al., 2014).

Climate warming causes increased stratification and deoxygenation of deeper waters, with implications for animal communities and associated geochemical processes. Deoxygenation decreases oxygen penetration in sediments, shallowing the depth of the oxic/anoxic interface. Hypoxia shifts biogeochemical cycles toward anaerobic metabolic pathways driven by nitrate and sulfate influxes (Middelburg and Levin, 2009). Oxygen concentrations can also alter the composition, size, biomass, density, and diversity of benthic communities (Levin, 2003; Levin et al., 2009; Figure 4), reducing bioturbation and bioirrigation (Smith et al., 2000), shifting toward microbial C cycling under low oxygen conditions (Middelburg and Levin, 2009; Gallo, 2018), and, in most cases, enhancing C preservation (Emerson and Hedges, 1988; Hedges and Keil, 1995).

Most low-oxygen settings also experience elevated CO₂ because respiration draws down oxygen, both in coastal waters and bathyal oxygen minimum zones (Breitburg et al., 2015; Gobler and Baumann, 2016). Yet despite common co-occurrence, few studies have manipulated both parameters experimentally, although meta-analyses along natural gradients have addressed this problem (Sperling et al., 2016). Multiple shallow-water bivalve species exhibit reduced growth and/or survival under either low oxygen or low pH conditions; together, these stressors often induce additive or synergistic negative effects (Figure 4), although high tolerance may exist in some species or life stages (Gobler and Baumann, 2016). Elevated CO₂ typically raises energetic requirements of invertebrates, inducing greater OC consumption (Thomsen et al., 2013; Queirós et al., 2015); the addition of hypoxia causes metabolic depression, potentially impairing how sediment fauna cope with acidification (Ravaglioli et al., 2019). Studies of the combined effects of reduced oxygen and elevated CO₂ suggest slower microbial degradation and enhanced carbon burial, although additional studies of coupled abiotic stressors involving oxygen are needed (Sampaio et al., 2021). Given food limitation in many benthic communities, decreased POC production, timing (e.g., Figure 5d), and export from surface waters with increased ocean stratification will alter benthic species biomass, function, and geochemical processes (Figure 4; Smith et al., 2008; Jones et al., 2014).

Several meta-analyses have identified climate-related thresholds for sustaining populations of benthic macrofauna, which depend on species tolerance to environmental change, reliance on environmental cues, population abundance, and/or inherent ability to acclimate, disperse, or adapt (Foden et al., 2013; Levin and Le Bris, 2015). Because environmental conditions (gradients in oxygen and temperature, in particular) shift geographically, we anticipate that the distribution of many benthic species will move correspondingly in accordance with their physiological thresholds. Biogeographical pattern ultimately depends on a confluence of ecological and evolutionary constraints that define a species niche. In the simplest form, a species niche can be represented by a gradient of abundance, with an assumed peak in abundance corresponding to the geographic center of a species distribution (“abundant centre” hypothesis sensu; Brown et al., 1995). Genetic diversity, correspondingly, should track spatial patterns, with a decline in intrapopulation diversity toward range limits. However, we lack empirical evidence of the “abundant centre” in marine systems, particularly in coastal areas (Sagarin and Gaines, 2002), and patterns of genetic diversity likely depend strongly on dispersal phenotype, whereby species with high dispersal capabilities likely show no strong decline in intrapopulation diversity toward range limits (Cahill and Levinton, 2016; Ntuli et al., 2020). For meroplankton, dispersal of planktonic larvae (Scheltema, 1989) can contribute to high population connectivity and therefore to increased gene exchange and reduced genetic isolation by distance (Cowen and Sponaugle, 2009; Modica et al., 2017; Ntuli et al., 2020). Estimates of direct dispersal distances from ocean currents and planktonic larval life span and estimates derived from measures of genetic isolation by distance regressions converge on larval dispersal distances of shelf benthic invertebrates on the order of 100–1,000 km (Kinlan and Gaines, 2003; Shanks, 2009). Lecithotrophic larvae and near-bottom distributions of some deposit-feeding bivalve larvae may exempt them from strong coastal currents above the benthic boundary layer, potentially reducing dispersal distances to 10 km or less. But even modest dispersal distances combined with thermal limitations at low latitudes increase the likelihood of expansion of the ranges of shelf species to higher latitudes (Pinsky et al., 2019).

Larval dispersal distances of deep-sea species can range from 0.24 to 2,028 km (average 33 km; Baco et al., 2016), depending on taxonomy and life history (e.g., adult and larval mobility), suggesting that the dispersal capacity of some taxa will allow them to track horizontal shifts in thermal environments. However, with or without the capacity for dispersal, benthic species may need to adapt to new environmental conditions in order to persist (Boyd et al., 2016; Bay et al., 2017). Several lines of evidence suggest the need for such evolutionary change. Large standing genetic variation in large populations tends to facilitate genetic adaptation (Bay et al., 2017), although physiologically relevant genes may lack sufficient variation to respond to directional selection. Emerging evidence points to distinct evolutionary dynamics across space within species; for instance, rapid evolution within populations moving into new habitats (i.e., leading edges) may favor traits that enable range expansion, whereas rapid evolution in historic habitats (i.e., trailing edges) may favor traits that inhibit range contraction (e.g., expanding environmental tolerance; Beaugrand and Kirby, 2018; Pinsky et al., 2019). Warming at leading edges, however, naturally allows for expansion to higher latitudes. Finally, the unknown role of epigenetic and transgenerational plasticity in hindering or facilitating adaptive responses, among other factors, hinders prediction of the trajectory of benthic responses. Documenting these responses in all marine biomes presents an enormous challenge for biologists.

The dynamics of seafloor carbon cycling are characterized by multiple timescales in both the climate forcing (e.g., modified seasonality, frequency and duration of hypoxia, slow but steady decline of bottom-water oxygen over decades, steady warming and heat waves) and the benthic response (e.g., multiple organisms that differ in size and turnover, community changes, order-of-magnitude ranges in organic matter lability). Numerical models complement observational data by providing the required extrapolation means for quantitatively analyzing and predicting OC-macrofaunal dynamics over the entire spectrum of temporally changing conditions. The temporal scales addressed by models vary from short-lived events (Lessin et al., 2018) and recovery to the feedbacks between biogeochemical cycles and the evolution of bioturbation on geological timescales (Boyle et al., 2014; Section III). Models designed to address decadal to centennial timescale responses to contemporary climate change (Section IV) are particularly relevant to our article. The spatial span considered by marine sediment models also varies widely, from resolving vertical gradients within the few decimeters of deposits in diagenetic models, to applications focusing on regional seas in coastal ecosystem models, to global-scale assessments of benthic fluxes, recycling, and sedimentary C sequestration in global ocean biogeochemical models.

The propagation of transient, or time-varying, conditions within the sediment leads to temporally changing stocks and fluxes constrained, on the one hand, by the duration, magnitude, and frequency of perturbations on the seafloor (Dale et al., 2008a; Middelburg and Levin, 2009) and, on the other hand, by the nature and rates of transport and reaction (Soetaert et al., 1996a). Four broad categories of natural or climate perturbations occur on the seafloor: (1) changes in quantity and/or quality of OC delivery from waters above (e.g., lability of phytodetritus or zooplankton pellets; Smith and Baldwin, 1984; Ruhl et al., 2008; Smith et al., 2008, 2013; Figure 5d); (2) changes in depositional environments due to transient transport phenomena, such as turbidity currents and storms (Lessin et al., 2018); (3) changes in environmental conditions, notably temperature, oxygen, or pH; and (4) changes in benthic communities resulting from bottom-up effects of particulate organic matter supply (Dame, 1993; Griffiths et al., 2017) and food-web alteration, such as predation pressure (e.g., Figure 5d; Singh et al., 2013; Ollivier et al., 2018). Direct human disturbance can also interact with and exacerbate climate change effects on benthic ecosystems (Bindoff et al., 2019). In the coastal zone, these disturbances can take the form of pollution, trawling and dredging, coastal development leading to habitat squeeze and altered hydrography, and excess delivery of nutrients that cause eutrophication and lead to harmful algal blooms and hypoxia (Middelburg and Levin, 2009; Bindoff et al., 2019; Ravaglioli et al., 2019; Ehrnsten et al., 2020). Overfishing and species invasions can also change the composition of benthic species and their interactions with the seafloor (Bindoff et al., 2019). Overfishing will greatly exacerbate the consequences of climate change, with 60% of exploited species vulnerable to both climate change and overfishing by 2050 in climate models under representative concentration pathway (RCP) 8.5, especially in tropical areas (Bindoff et al., 2019). Understanding the direct and indirect consequences of these changes for benthos remains limited.

Despite the potential importance of transient conditions for C cycling, most diagenetic models (benthic biogeochemical models addressing the dynamics at the typical decimeter scale of sediment deposits; Berner, 1980; Boudreau, 1996) assume a steady state (Arndt et al., 2013). Notable exceptions include models of the transient response of microbial communities to abrupt or periodic shifts in transport regimes (Wirtz, 2003; Dale et al., 2008b) or changes in biogeochemical fluxes induced by seasonal or decadal changes in bottom water temperature, oxygen, sediment, and OC delivery (Soetaert et al., 1996b; Katsev et al., 2007; van de Velde et al., 2018). However, only a few models (Katsev et al., 2007; Soetaert and Middelburg, 2009) include feedback of temporally changing oxygen levels on biologically induced transport using a highly parameterized approach. Many studies have investigated qualitatively the response of macrobenthos to changing environmental conditions at different spatial and temporal scales and the extent to which irrigation, particle-reworking scales, and intensity track community changes (Pearson and Rosenberg, 1978; Gerino et al., 1998; Diaz and Rosenberg, 2008; Middelburg and Levin, 2009; Van Colen et al., 2012). Conceptual models elucidating animal–sediment interactions therefore exist but, because suitable parameterizations are lacking, none have translated this knowledge into the quantitative and predictive framework of diagenetic models (Peña et al., 2010).

Assessment at the local scale of sediment deposits must consider several important challenges in order to address organism–sediment interactions and their transient responses. On the one hand, current ecological models can estimate meiofaunal and macrofaunal secondary production (Wei et al., 2010) but ignore the constraints imposed by biogeochemical dynamics and small-scale vertical sediment and sediment–water interface gradients. On the other hand, biogeochemical diagenetic models assume steady benthic biomass and functioning (Burdige, 2006; Middelburg, 2018), an unrealistic assumption given changing environmental conditions, biological diversity, and species distributions on the seafloor (Figure 5c and e). Fundamental differences in the characteristic spatial scales of ecological (typically centimeter and beyond) and biogeochemical (typically millimeter) processes (Meysman et al., 2002) and the multiple timescales inherent in sediment biogeochemistry (Katsev et al., 2007) complicate efforts to merge the two modeling approaches. Incorporation of microbial processes at the interface between sediment biogeochemistry and faunal dynamics adds an additional major challenge (Middelburg and Levin, 2009; Middelburg, 2018). This missing microbial link limits the incorporation and parameterization of key concepts in diagenetic models, such as the inverted microbial loop (e.g., microbes profiting from animal-induced sediment mixing and detrital shredding) and the associated positive priming interactions induced by biologically driven penetration of labile OC (Bianchi, 2011; Middelburg, 2018; Aller and Cochran, 2019). Theoretical models of priming exist (Canfield, 1994; Arndt et al., 2013) but lack real-case applications, especially given that biogenic transport, priming, and resulting OC decomposition vary in time. Inclusion of nonadditive effects of multiple macrofaunal species on net biogeochemical fluxes and microbial community composition and function constitutes another major challenge (Michaud et al., 2009).

Therefore, to fully couple descriptions of feedback mechanisms between biogeochemical, microbial, and animal dynamics (e.g., Figure 5) and their parameterization will require incorporation of several new theoretical concepts within traditional early diagenetic models. These modifications need to include mathematical representations on how moving organisms, treated as a separate phase, affect the sedimentary mass and momentum balance (Meysman et al., 2002) and on how faunal dynamics control the benthic microbial loop through the downward transport of labile organic matter and the subsequent degradation of recalcitrant organic matter by positive priming (Middelburg, 2018). Key processes to include in diagenetic models are representation of viral lysis as a dominant loss term for microbial communities (Danovaro et al., 2008a, 2008b; in addition to grazing) and the efficient reassimilation of dissolved OC produced by microbes in a closed loop at the bottom of the food chain; both processes may significantly limit C transfer to higher trophic levels (Middelburg, 2018).

Following the pioneering work by Soetaert et al. (2001), several benthic–pelagic models have been developed in parallel over the last two decades to help unravel the biogeochemical dynamics at the larger scale of coastal ecosystems. These ecosystem models have been applied to estuaries and bays (Arndt and Regnier, 2007) as well as to regional seas, most notably the North Sea and the Baltic Sea (Reed et al., 2011; Griffiths et al. 2017). However, few included dynamic benthic ecological elements, specifically coupling biogeochemical fluxes and microbes (Yakushev et al., 2017), or feedbacks between OC supply and macrofaunal biomass (Ehrnsten et al., 2019) and their responses to hypoxia and nutrient loading (Soetaert and Middelburg, 2009). Furthermore, ecosystem models only recently linked macrofaunal dynamics to changes in biological transport and their effects on biogeochemical fluxes (Daewel and Schrum, 2013; Butenschön et al., 2016; Zhang and Wirtz, 2017), inducing strong temporal dynamics resulting from high sensitivity to environmental drivers (Timmermann et al., 2012; Ehrnsten et al., 2019; Zhang et al., 2019). The major limitation of such coastal ecosystem models rests in their highly simplified representation of the vertical sedimentary structure, which cannot resolve strong biogeochemical and biological spatial gradients with sediment depth. Recent explicit representations of organism–sediment interactions as a function of depth in sediment cores (Zhang et al., 2019) nonetheless link dynamics to local macrobenthic biomass and food resources, both of which may vary temporally and spatially.

At a global scale, few diagenetic model applications have investigated how biological transport influences OC processing and associated benthic nutrient fluxes. These applications typically rely on mathematical formulations for bioturbation and bioirrigation intensities (Meile and Cappellen, 2003) that are similar to those of local diagenetic model applications but at the global scale, are based on a series of model realizations. In particular, global simulations are constrained by empirical relationships for the (assumed time-invariant) parameters controlling biologically induced mixing, such as OC deposition flux and benthic oxygen levels (Soetaert et al., 1996a; Thullner et al., 2009; Dale et al., 2015). Because they typically capture spatial variability along a global hypsometry, these approaches only provide first-order quantitative estimates of benthic processes. Results, nonetheless, clearly point to the importance of biologically induced transport as a key process in quantifying sedimentary biogeochemical cycles in sediments at the global scale and offer predictive potential for new boundary conditions.

In principle, these global-scale applications could be extended to simulate the biogeochemical response of sedimentary systems to anthropogenic perturbations or natural variations at seasonal to decadal timescales, for example, by using variable temperature, OC input, and oxygen, nitrate or carbonate ion in the bottom waters as forcing fields, extracted from global data sets such as the National Oceanic and Atmospheric Administration World Ocean Atlas or from outputs of existing regional and global ocean models (e.g., Krumins et al., 2013; Yakushev et al., 2017). Similarly, several studies partly addressed the response of the benthos to these same global environmental drivers (Jones et al., 2014; Yool et al., 2017). For instance, due to decreased export production, models combined with statistical approaches (Wei et al., 2010) predict a 5.2% overall decrease in global ocean benthic biomass by the end of the century under RCP 8.5, although biomass may increase in regions such as the polar oceans and areas of increased temperature (i.e., Arctic temperature anomalies in Figure 5) and productivity (e.g., upwelling zones; Jones et al., 2014). In addition to declining seafloor biomass, overall body size is predicted to generally decrease (e.g., Figure 5f), with greater declines anticipated for macrofaunal biomass than in smaller groups (e.g., meiofauna), particularly in deeper waters (Yool et al., 2017). However, because no global models currently link these biogeochemical and ecological responses dynamically, we cannot project the coupled evolution of sedimentary OC cycling and macrobenthos activity at the ocean basin and global scales. This conclusion has even greater validity for biogeochemical and ecosystem components of global ocean biogeochemical models embedded in Earth system models, which either neglect or oversimplify sedimentary biogeochemical cycles (Soetaert et al., 2000; Hülse et al., 2018). Their representation typically assumes reflective boundaries, or conservative semi-reflective boundaries, or simulates benthic processes in a thin, homogeneous bioturbated layer (Soetaert et al., 2000). The Norwegian Earth System Model (Tjiputra et al., 2013) and Hamburg Ocean Carbon Cycle model (Ilyina et al., 2013) are exceptions, their configuration relying on a vertically resolved representation of biogeochemical fluxes and (time-invariant) bioturbation intensity. Moreover, the influence of climate change on macrobenthos remains largely conceptual (Figure 4), and few studies explicitly link environmental conditions to the distribution and adaptation of macrobenthos, thus limiting the ability of ocean biogeochemical models to resolve expected species responses to climate change spatially beyond indirect links to known ecosystem responses (Figure 5). Consequently, no global ocean biogeochemical model currently explicitly represents macrofaunal activity and its spatial and temporal variability across the global seafloor.

The vertically resolved coupling of OC degradation, macrobenthos dynamics, and biologically induced transport recently developed at the scale of coastal ecosystem models could provide an important building block for a unified modeling framework (Zhang et al., 2019). These ecosystem models do not resolve the sharp biogeochemical gradients necessary to predict the impact of changing macrobenthos and associated transport regimes on biogeochemical cycles. However, the newly developed, numerically efficient diagenesis modules, designed for coupling to large-scale ocean models, resolve the full redox sequence while including controls by bioturbation and irrigation using a highly parameterized approach (Soetaert et al., 2000; Hülse et al., 2018), thus offering a converging, complementary basis for further model development. We must develop approaches that apply such models to systems characterized by fluctuating boundary conditions and changes in transport regimes; such techniques exist with successful application to simulate, for example, benthic nutrient regeneration in regional seas (Ruardij and Van Raaphorst, 1995).

The ocean floor and its inhabitants provide multiple key ecosystem services, including C degradation, storage and burial, and nutrient transformations. Global climate change could potentially disrupt many of these services through species, population, and community-level responses to climate and oceanographic conditions. Physiology and adaptive variation will drive many such responses, where species or populations adapted to specific environmental conditions (ecological niche) track environmental change. Some benthic organisms may change spatially through (generally) poleward redistributions (Stanley et al., 2018) or altered bathymetric distributions (Figure 5a and c). Others will change temporally through the timing of biological processes (Edwards and Richardson, 2004; Figure 5b), compositionally through reorganization of benthic communities and biodiversity (Danovaro et al., 2004; Figure 5e), or phenotypically, where warmer environments favor smaller-bodied individuals (Cheung et al., 2012; Figure 5f), which will affect animal–sediment biochemical kinetics. Changes in the pelagic environment, including timing of phytoplankton blooms (Figure 5d), can profoundly influence primary and secondary productivity (Lotze et al., 2019) and function within these ecosystems, in turn altering benthic–pelagic exchange (coupling) and carbon subsidies (POC flux) to benthic ecosystems. Collectively, these responses influence marine food webs through changes to predation, trophic structure (cascades) and, ultimately, energy transfer. Despite generally limited examples of spatially explicit responses for macrobenthos, these emblematic examples serve as functional surrogates to spatial–temporal changes expected of macrobenthic species, in addition to indirect pathways of effect that will augment the functional role of benthos for global biogeochemical cycling.

We strongly advocate that large-scale ocean models of C cycling must better represent benthic biological activity and associated transport modes for several reasons, despite limited current understanding of animal responses. First, recently published (observation-based) spatial mapping of bioturbation and ventilation, including some benthic faunal types (Solan et al., 2019), may significantly improve modeling of sedimentary processes across scales from local to global. Second, recent developments of system-scale ecosystem models and simulation of OC-macrobenthos interactions provide a valuable theoretical framework for global modelers interested in the dynamic coupling of secondary benthic production and biogeochemical fluxes. Third, ocean biogeochemical models applied at regional and/or global scales can now run at spatial resolutions that improve representation of biogeochemical dynamics on continental shelves and slopes, for example, 0.2°–0.5° globally (Bourgeois et al., 2016; Lacroix et al., 2021) and 0.1° for the Tropics (Busecke et al., 2019). In addition, higher resolutions can be obtained for specific shelf systems or ocean basins with regional models or via downscaled versions of global models (Busecke et al., 2019). Given the importance of large areas of the ocean margins as hot spots of intense biological–physical and benthic–pelagic coupling, and where changing environmental and climatic conditions will likely impact benthic fauna most, we should now expand the capabilities of models toward dynamic representation of organism–sediment couplings on the shelves and upper slopes. Management applications for such models extend beyond improved understanding of climate feedbacks to use in scenario planning for climate adaptation, design of marine conservation networks, fisheries management, development of coastal water quality criteria, and regulation of nutrients and sewage (Roberts et al., 2017).

To enable fully coupled descriptions of biogeochemical, microbial, and animal dynamics as input to models, we must evaluate how environmental drivers regulate the distribution of seabed macrofauna, how adaptive variation regulates their responses to environmental change (Stanley et al., 2018), how species functions regulate sedimentary OC cycling, and, ultimately, the cumulative effect of these drivers on microbial and biogeochemical processes. Experiments that consider cumulative impacts are limited (Belley and Snelgrove, 2016; Hale et al., 2017) and detailed data on functional group composition and activity for most seabed environments remain far from reach, but we could address the larger need through change proxies and indicator species. Addressing major unanswered questions about how climate modification of plankton communities will influence the export of organic matter to the seabed (e.g., Figure 5d), how infauna will adapt or respond in terms of biomass and function, and what the consequences will be for OC remineralization and sequestration in the coastal ocean will all require further mechanistic research. For macrobenthos, we particularly lack a quantitative evaluation of the relationships among environmental conditions, species distribution, and function, as well as how genetic adaptation regulates responses to climate change. Understanding these causal relationships will underlie our ability to resolve spatially how macrobenthos and their associated functional roles regulate C cycling.

The recent advances in seafloor habit mapping offer much potential for upscaling local observations on sediment characteristics and animal community composition to the basin or global scale, but they can only be used for estimation of biogeochemical processes at regional or global scale through better understanding of the two-way interactions between benthic animals and microbes. Environmental conditions regulate the performance of both animals and microbes via direct (e.g., predation) and indirect (e.g., inverted microbial loop) interactions as organisms compete for common resources. Moreover, quantifying the holistic response of benthic ecosystems to climate perturbations and resulting large-scale biogeochemical feedbacks requires further elucidation of interaction mechanisms. Achieving this goal will require coordinated efforts among ecologists, biogeochemists, and climate modelers (Bianchi et al., 2021). The grand challenges identified here could be addressed through integration into ongoing international programs including the Global Ocean Observing System and the Decade for Ocean Science for Sustainable Development.

All data generated for this work are stored on a University of Florida server cluster that has off-site replication and is available upon request. Data derived from previously published work are available online:

• Animal–substrate interactions through geologic time available in Buatois and Mángano (2018).

• Habitat centroid projections for benthic species on the North American Shelf available as online supplement for Morley et al. (2018).

• Projected changes in diapause duration of Calanus finnmarchicus available in Wilson et al. (2016).

• Projected rates of change in the depth distribution for fish and invertebrate specdies in the North Sea available in Dulvey et al. (2008).

• Projected changes in biomass aggregated by species functional groups available in Lotze et al. (2019).

• Projected changes in species richness available as online supplement for García Molinos et al. (2015).

• Projected median changes in maximum fish body size available in Cheung et al. (2012).

• Data were derived from Coupled Model Intercomparison Project 5 (CMIP5) and downloaded from the National Oceanic and Atmospheric Administration (2019) ESRL Climate Change Web Portal.

We thank Craig Smith and an anonymous reviewer whose detailed comments substantially improved this article.

This work was supported in part by the National Science Foundation’s (NSF) Chemical Oceanography, Low Temperature Geochemistry and Geobiology, Arctic Natural Sciences, Biological Oceanography, and Ecosystem Ecology programs; NASA Carbon Cycle Science program; NOAA; and DOE. LAL is supported by NOAA CHRP award NA18NOS4780172 and NSF awards OCE 1634172 and 1829623; PVRS and LAB, by the Natural Sciences and Engineering Research Council of Canada; and JJM, by the Netherlands Earth System Science Center. PR received funding from the VERIFY Project from the European Union’s Horizon 2020 research and innovation program under grant agreement No. 776810 and from the F.R.S.-FNRS (sabbatical leave). TBA was funded by an Early Career Research Fellowship from the Gulf Research Program of the National Academies of Sciences, Engineering, and Medicine. The content is solely the responsibility of the authors and does not necessarily represent the official views of the Gulf Research Program of the National Academies of Sciences, Engineering, and Medicine.

The authors declare no competing financial interests.

All coauthors contributed equally to the conceptual development and writing of this article. In particular, figures and tables were developed by the following coauthors: Figure 1 (CJB), Figure 2 (MJS and TBA), Figure 3 (LAB), Figure 4 (ESM and LAL), and Figure 5 and Table 1 (RRES).

How to cite this article: Bianchi, TS, Aller, RC, Atwood, TB, Brown, CJ, Buatois, LA, Levin, LA, Levinton, JS, Middelburg, JJ, Morrison, ES, Regnier, P, Shields, MR, Snelgrove, PVR, Sotka, EE, Stanley, RRE. 2021. What global biogeochemical consequences will marine animal-sediment interactions have during climate change? Elementa: Science of the Anthropocene 9(1). DOI: https://doi.org/10.1525/elementa.2020.00180

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Knowledge Domain: Ocean Science