An Integrative Model of Carbon and Nitrogen Metabolism in a Common Deep-Sea Sponge (Geodia barretti)

Deep-sea sponges and their microbial symbionts transform various forms of carbon (C) and nitrogen (N) via several metabolic pathways, which, for a large part, are poorly quantified. Previous flux studies on the common deep-sea sponge Geodia barretti consistently revealed net consumption of dissolved organic carbon (DOC) and oxygen (O2) and net release of nitrate (NO3-). Here we present a biogeochemical metabolic network model that, for the first time, quantifies C and N fluxes within the sponge holobiont in a consistent manner, including many poorly constrained metabolic conversions. Using two datasets covering a range of individual G. barretti sizes (10–3,500 ml), we found that the variability in metabolic rates partially resulted from body size as O2 uptake allometrically scales with sponge volume. Our model analysis confirmed that dissolved organic matter (DOM), with an estimated C:N ratio of 7.7 ± 1.4, is the main energy source of G. barretti. DOM is primarily used for aerobic respiration, then for dissimilatory NO3- reduction to ammonium (NH4+) (DNRA), and, lastly, for denitrification. Dissolved organic carbon (DOC) production efficiencies (production/assimilation) were estimated as 24 ± 8% (larger individuals) and 31 ± 9% (smaller individuals), so most DOC was respired to carbon dioxide (CO2), which was released in a net ratio of 0.77–0.81 to O2 consumption. Internally produced NH4+ from cellular excretion and DNRA fueled nitrification. Nitrification-associated chemoautotrophic production contributed 5.1–6.7 ± 3.0% to total sponge production. While overall metabolic patterns were rather independent of sponge size, (volume-)specific rates were lower in larger sponges compared to smaller individuals. Specific biomass production rates were 0.16% day–1 in smaller compared to 0.067% day–1 in larger G. barretti as expected for slow-growing deep-sea organisms. Collectively, our approach shows that metabolic modeling of hard-to-reach, deep-water sponges can be used to predict community-based biogeochemical fluxes and sponge production that will facilitate further investigations on the functional integration and the ecological significance of sponge aggregations in deep-sea ecosystems.


INTRODUCTION
Sponges are abundant and key ecosystem engineers of the deep sea that occur scattered on soft-and hard-bottomed surfaces and in multi-or mono-specific aggregations (Maldonado et al., 2017). These so-called sponge grounds create complex habitats and thereby support high local biodiversity (Klitgaard, 1995;Hogg et al., 2010;Beazley et al., 2013). As deep-sea sponges process large amounts of water for filter feeding, they are implicated to have an important role in the biogeochemical cycling and benthic-pelagic coupling of carbon (C), nitrogen (N), and silicon (Si) (Maldonado et al., 2012Kutti et al., 2013;Kahn et al., 2015).
Deep-sea sponges efficiently filter (preferably nano-and pico-) plankton (Pile and Young, 2006;Yahel et al., 2007;Kahn et al., 2015) and, similar to their shallow counter parts (de Goeij et al., 2013), consume dissolved organic carbon (DOC), constituting most of their diet (Bart et al., 2020b). The efficiency at which sponges process (e.g., assimilate, respire, release) organic and inorganic nutrients (e.g., C, N) is termed production efficiency, which is an important ecological and metabolic parameter to determine energy (re)cycling in organisms (also referred to as growth efficiency) (Thomassen and Riisgard, 1995;Sterner and Elser, 2002;Maldonado et al., 2012;de Goeij et al., 2017) and ecosystems (Keesing et al., 2013;Kahn et al., 2015;Porada et al., 2018). Despite its potential relevance, the C and N production efficiency of deep-sea sponges and most other deepsea benthos are largely unknown-especially in situ assessments of metabolic rates and conversions-because of the technical constraints of conducting (properly controlled) experiments in the deep sea. Recent advances in technology and accessibility of remotely operated vehicle have increased the ability to collect live deep-sea sponges for ex situ-controlled laboratory experiments (e.g., Kutti et al., 2015;Leys et al., 2018;Bart et al., 2020a,b) and even in situ metabolic rate measurements (Yahel et al., 2007;Kahn et al., 2015;Maldonado et al., 2020a).
Geodia barretti (Porifera, Demospongiae) is one of the best studied deep-sea sponge species. In the North-Atlantic Ocean, along the continental shelf and in fjords, the massive, ballshaped G. barretti can be present in large densities of up to 0.4-5 individuals per square meter, known as "Geodia grounds" (Murillo et al., 2012;Kutti et al., 2013;Beazley et al., 2015). It is considered to be a slow-growing sponge species that, if undisturbed, can reach a meter in diameter, although most specimens have a diameter of 20-30 cm (Klitgaard and Tendal, 2004;Kutti et al., 2013). Oxygen (O 2 ) consumption by G. barretti individuals, combined with biomass estimates from video imaging, indicated that Geodia sponge grounds have high metabolic activity and C demand (ex situ, Kutti et al., 2013;eddy correlation, Cathalot et al., 2015). To date, the metabolic C demands of G. barretti have been derived from O 2 consumption measurements and a supposed respiratory quotient (RQ) for organic matter (OM) mineralization (O 2 :CO 2 ). Despite the importance of RQ in estimating the C demand from O 2 measurements, its value has never been empirically determined for G. barretti. The RQ values for G. barretti might deviate from values for canonical organic matter mineralization because of the various metabolic pathways conducted by the endosymbionts.
G. barretti is considered as a high microbial abundance (HMA) species, containing a dense and diverse community of microbial symbionts (∼10 11 microbes per cubic centimeter) (Hoffmann et al., 2009;Leys et al., 2018). Microbial symbionts (i.e., the microbiome) and sponge host are a metabolically integrated functional unit, known as a "holobiont." This, generally considered mutualistic, symbiosis is known to benefit in a variety of ways, which include nutrition, development, defense, and immunity (Pita et al., 2018). In G. barretti, the microbiome is actively involved in the (re)cycling of C and N compounds and, therefore impacts the holobiont nutrition (Hoffmann et al., 2009;Leys et al., 2018;Bart et al., 2020a). As for many sponge species, anaerobic and aerobic N-transforming processes occur in G. barretti, thereby expanding its metabolic capacity (Hoffmann et al., 2009;Rooks et al., 2020). It has been hypothesized that G. barretti actively controls the oxygenation level of its tissue to maintain aerobic and anaerobic microenvironments that facilitate the co-existence of aerobic and anaerobic symbionts (Hoffmann et al., 2005). The nitrifying and denitrifying potential of the microbial populations within the tissue of G. barretti has been demonstrated by cutting sponges in small fragments (0.3-0.4 cm 3 ) and incubating those pieces under labeled substrates (Hoffmann et al., 2009;Rooks et al., 2020). However, N transformations within intact, pumping G. barretti individuals have yet to be resolved and quantified.
Two studies have examined the C, N, and O 2 exchange rates of intact G. barretti individuals. Leys et al. (2018) studied larger specimens (∼1-L volume) using the "in-ex" methodology, i.e., sponge filtration rates are combined with the composition of inhalant and exhalant water to obtain release/uptake fluxes. Bart et al. (2020b) used incubation flow chambers to quantify fluxes for smaller (∼100-ml volume) specimens. Both studies found that G. barretti effectively removed organic C and concluded that the observed organic matter consumption was sufficient to sustain the minimal metabolic needs. Bart et al. (2020b) found a relatively higher net organic C intake compared to Leys et al. (2018). This difference can be attributed to the smaller size of the individuals tested by Bart et al. (2020b), as volume-specific pumping rates may be related to sponge size (Morganti et al., 2019). Regarding N, intact G. barretti was found to release nitrate (NO − 3 ) rather than ammonium (NH + 4 ), indicative of net nitrification (Leys et al., 2018), yet the metabolic conversions of O 2 , C, and N within and exchanges among sponge and symbionts can be complex, and it is not straightforward to derive those from net fluxes obtained by the in-ex and chamber incubation techniques.
To integrate existing experimental metabolic flux data, we constructed a metabolic network model for the holobiont. Metabolic (or stoichiometric) network models are common tools to quantify and predict the intracellular metabolism of, for example, eukaryotic cells and microbes (Palsson and Varma, 1994;Edwards and Covert, 2002) based on net resource (e.g., C, N) fluxes. Such models are underexplored to quantify sponge-symbiont metabolism, despite their potential to elucidate the intricacies of the sponge physiology and their ecological significance. To date, only one study that employs a genomescale metabolic network model to reconstruct biochemical conversions in the sponge-symbiont system of the tropical sponge Amphimedon queenslandica exists (Watson, 2017). In our study, we applied the metabolic network model to quantify internal and unconstrained C, N, and O 2 fluxes of the common and abundant deep-sea G. barretti holobiont. To assess the size dependency of G. barretti metabolism, we compared and separately analyzed the metabolic datasets of Leys et al. (2018) and Bart et al. (2020b) because they assessed different size classes. The size dependency of metabolism is particularly relevant when estimating the metabolic fluxes and metabolic demands of sponge grounds from video images. The outcome of the model advances the understanding on G. barretti metabolism and its potential role in biogeochemical cycling within the vast, benthic deep-sea ecosystems where it abounds.

Datasets
The experimental data of Leys et al. (2018) and Bart et al. (2020b) were used to develop and constrain the metabolic network model for larger and smaller G. barretti, respectively, since they cover different size ranges. We refer to their publications (Leys et al., 2018;Bart et al., 2020b) for a full description of sponge collection, experimental procedures, data analyses, and data presentation and only briefly recap some relevant aspects here.
In short, Leys et al. (2018) maintained G. barretti individuals in ex situ (running seawater) conditions in the Institute of Marine Research aquarium facilities at Austevoll (Norway) and measured exchange fluxes using the "in-ex" method. The flow rates and O 2 uptake of G. barretti were collected simultaneously in 2012 (n = 17, sponge volume = 150-3,500 ml) ( Table 1). Dissolved inorganic N [DIN, consisting of NH + 4 , nitrite (NO − 2 ), and NO − 3 ] exchange rates and bacterial cell removal were measured in 2011 (n = 53, sponge volume = 125-3,400 ml), and DIN exchange rates and total OC (TOC = DOC + POC) removal were measured in 2014 (n = 29, sponge volumes unknown) ( Table 1). Leys et al. (2018) converted bacterial cell removal to bacterial C uptake using standard conversion factors and defined DOC uptake as TOC uptake minus bacterial C uptake (Table 1). Bart et al. (2020b) measured exchange fluxes in G. barretti individuals (n = 12, sponge volume = 9-210 ml) ( Table 1) using incubation flow chambers under ex situ (running seawater) conditions in aquarium facilities in Bergen (Norway) in 2017 (n = 9) and 2018 (n = 3). O 2 , DOC, and living POC (LPOC, consisting of bacteria and phytoplankton C) removal rates were calculated from changes in incubation water and used to calculate a C budget ( Table 1; Bart et al., 2020b). For both studies, the incubation water was collected from 160-m depth (Leys et al., 2018) and 200-m depth (Bart et al., 2020b) from nearby fjords, thus representing deep water. This corresponds with the depth range of 200-300 m where G. barretti grounds are typically found (Klitgaard and Tendal, 2004).

DIN Exchange Rates From Incubation Chambers
Changes in DIN (NH + 4 , NO − 2 , and NO − 3 ) were also measured by Bart et al. (2020b), but as these data were not presented before, the procedures and the results are presented in this paper. DIN fluxes were determined as follows: duplicate samples for DIN (NH + 4 , NO − 2 , and NO − 3 ) were collected at t = 0, 1, 2, 3, and 4 h in 2017 and t = 0, 2, 4, 6, and 8 h in 2018 from the incubation chambers with acid-washed 100-ml polycarbonate syringes. The samples were subsequently filtered over sterile 0.2-µm polyethersulfone syringe filters (Whatman Puradisc), collected in 10-ml high-density polyethylene vials, and stored at −20 • C until further analysis. Nutrients were analyzed in the laboratory with a QuAAtro Gas Segmented Continuous Flow Analyzer. Measurements were made simultaneously for NH + 4 (Helder and de Vries, 1979), NO − 2 , and NO − 3 combined with NO − 2 (Grasshoff et al., 2009). All measurements were calibrated with standards diluted in low-nutrient seawater. DIN exchange rates were calculated from a linear regression slope and incubation volume and corrected with changes in control (i.e., no sponge) incubations. Fluxes were normalized to sponge volume to obtain fluxes in µmol N cm −3 (sponge) day −1 (Table 1, Supplementary Figure 1).

Linear Inverse Modeling Concepts
Metabolic network models are a type of linear inverse modeling (LIM), which are based on mass balances under steady-state conditions. This implies that the internal changes of a substance due to chemical reactions are balanced by the inflow and outflow into a system and that the system can be solved by a set of linear equations (Edwards and Covert, 2002). The steady-state metabolic situation is justified by the rationale that chemical reactions take place on shorter timescales than transport (exchange) fluxes (Edwards and Covert, 2002). As detailed through subsequent sections, LIM of stoichiometric (metabolic) networks requires the following: (1) selection of relevant substances and (chemical) processes and conceptualization of the metabolic network, (2) stoichiometrically balancing of all chemical reactions and construction of mass balances for each substance, (3) constraining model solutions with actual flux measurements or fluxes derived from literature (although LIM can solve for or constrain missing flux data), and (4) mathematically solving the model, including variability and uncertainty analysis (Edwards and Covert, 2002;Soetaert and van Oevelen, 2009;van Oevelen et al., 2010).

Sponge Metabolic Network Model Concept
The metabolic network model integrates the C, N, and O 2 metabolism of G. barretti and consists of both net holobiont metabolism and host versus symbiont-driven metabolism (Figure 1 and Table 2). In detail, sponges rely largely on heterotrophy as they filter the water column for LPOC, consisting predominantly of bacterioplankton in the deep sea (Bart et al., 2020b; Table 1) (r1), and consume dissolved organic matter (DOM) (r11). Sponges aerobically respire assimilated organic matter to obtain energy for basal metabolism (r1-r2).  (Leys et al., 2018;Bart et al., 2020b) which were analyzed with the metabolic network model.
Dissolved organic carbon, DOC (r11) 25.7 ± 12.3 b 5.2 ± 2.1 O 2 (r12) 11.2 ± 8.1 7.7 ± 5.8 Positive fluxes indicate uptake, and negative fluxes indicate release. a Phytoplankton-C assimilation was not included in the model because of its small contribution to organic C uptake. b The DOC range is larger than the range reported by Bart et al. (2020b) to match the variation in O 2 uptake rates. c Dissolved inorganic N flux data are analyzed as part of this study.
FIGURE 1 | Schematic representation of the metabolic network model for a sponge holobiont (large yellow box) containing nitrogen species (brown boxes), organic substances (green boxes), and O 2 /CO 2 (blue boxes). Production (filled green) is not a standing stock, but a closure term for production of sponge biomass. The arrows indicate metabolic transformations: solid arrows refer to heterotrophic processes, dashed arrows refer to chemoautotrophy, and dotted arrows indicate exchange with the environment. Reaction numbers correspond with those in Table 2: blue numbers indicate aerobic processes, and orange numbers indicate anaerobic processes.
Frontiers in Marine Science | www.frontiersin.org
The metabolic waste of N is mineralized to NH + 4 (r1-r2). The assimilated organic matter that is not required for basal metabolism is allocated to biomass production (r10). The C (or N) balance for sponges applied in the network model is: assimilation = production + respiration (mineralization in case of N) (Thomassen and Riisgard, 1995;Maldonado et al., 2012).
Microbial symbionts expand the metabolic capacity of the sponge host by exploiting metabolic pathways additional to mere aerobic respiration, such as heterotrophy coupled to denitrification (r3) (Hoffmann et al., 2005;Rooks et al., 2020) and dissimilatory NO − 3 reduction to NH + 4 (DNRA) (r4) in anoxic parts of the sponge. The latter process has not been measured, but the accompanying genes (napA, nrfA) were found in G. barretti (Gavriilidou, personal communication). The included chemoautotrophic processes are nitrification in oxic (r5, r6) and anaerobic NH + 4 oxidation (anammox) (r8) in anoxic parts of the sponge holobiont (Hoffmann et al., 2005;Rooks et al., 2020). N 2 fixation is not considered in the model since, at present, no genetic evidence exists for G. barretti (Rooks et al., 2020). All considered microbial processes yield production (growth) of symbionts, either via anaerobic heterotrophy (denitrification and DNRA) or chemoautotrophy (nitrification and anammox) (r7, r9; Table 2). The production of symbiont biomass is considered an integrative part of sponge holobiont production, as the sponge host can phagocyte and ingest (in an undetermined amount) microbial symbionts that may continuously grow within the sponge tissue, including chemoautotrophic bacteria. The process of NH + 4 assimilation by symbionts (for which genetic evidence exists; Gavriilidou, personal communication) is considered implicit in organic matter consumption (r1-r2), where it functions as an internal feedback loop (from NH + 4 to production). Because of the steady-state approach, production of sponge (holobiont) biomass is treated as a closure term (r10).

Reactions and Coefficients
All reactions are stoichiometrically coupled for N and C, with N as base and C derived from the stoichiometry of the substrates. For example, in reaction 2, assimilated DOM is converted to sponge (holobiont) biomass ("Spo") with a N production efficiency (f p_N DOM ), while the remaining fraction is excreted as NH + 4 (1−f p_N DOM ). The C production efficiency (f p_C DOM ) and the fraction that is respired to CO 2 (f r_C DOM = 1−f p_C DOM ) depends on the stoichiometry (molar C:N ratio) of DOM and sponge biomass: f p_C DOM = f p_N DOM · C:N spo / C:N DOM ( Table 2). Reactions 1, 3, and 4 have a similar structure as reaction 2 ( Table 2). The measured (molar) stoichiometry of G. barretti biomass "Spo" (C:N spo ) is 4.2 ± 0.29 [n = 20 specimens collected by Bart et al. (2020b) and Maier et al. (2020)]. The C:N ratios of assimilated food sources (bacteria and DOM) are unconstrained, but experimental studies indicated that G. barretti preferentially assimilates N-rich DOM, with C:N ratios close to or below Redfield (C:N < 6.6) (Maier et al., 2020;Bart et al., 2020a). The C:N ratios of marine bacteria (C:N bac ) are typically 4 to 6, thus below the Redfield ratio (Sterner and Elser, 2002;Zimmerman et al., 2014;White et al., 2019). The C:N ratios of seawater DOM (C:N DOM ) are generally higher, with average C:N ratios of labile DOM of 10.7 (Hopkinson and Vallino, 2005). The C production efficiencies of metazoans typically range from ∼20 to 30%, which would correspond to f p_N DOM of 32-48% on Redfield DOM (C:N DOM = 6.6). However, G. barretti might have a higher f p_N DOM as heterotrophic symbionts can assimilate NH + 4 into organic N (Hentschel et al., 2012;Gavriilidou pers. communication). The C production efficiency for denitrification (f p_C denit ) is ∼25% (Koike and Hattori, 1975;van den Berg et al., 2016), and DNRA has a comparable production efficiency (f p_C DNRA ) (Strohm et al., 2007;van den Berg et al., 2016). Depending on the stoichiometry of DOM, the N production efficiencies (f p_N denit ) and (f p_N DNRA ) are expected to range from 30 to 60%. The variability and uncertainty in C:N of the assimilated substrates and production efficiencies were considered in model solutions ( Table 2, "Section Model Implementation and Solutions").
Oxygen is used as electron acceptor in r1 and r2, and the respiratory quotient (RQ; O 2 consumed to CO 2 produced) depends on the elemental composition of the substrate and the fraction of substrate that is respired (Table 2 and Figure 1). RQ ranges from 1 for oxidation of glucose to 1.16 for complete mineralization of plankton Redfield OM (C:N = 6.6) to NH + 4 (Del Giorgio and Williams, 2005;Middelburg, 2019). In anoxic parts of the sponge tissue, NO − 3 may be used as an electron acceptor and is reduced to N 2 in denitrification (r3, Table 2) and to NH + 4 in DNRA (r4, Table 2). NO − 3 consumption per CO 2 produced for Redfield OM in denitrification (ϕ denit ) is 0.8 (r3; Soetaert et al., 1996;Strohm et al., 2007). In DNRA, this ratio (ϕ DNRA ) is determined as 0.56 (r4; Koeve and Kähler, 2010).
In the oxic compartments of the sponge holobiont, NH + 4 is oxidized to NO − 2 (r5) and NO − 2 is oxidized to NO − 3 (r6). NH + 4 oxidation in G. barretti is primarily conducted by archaea, while NO − 2 oxidizers are primarily Nitrospira bacteria (Radax et al., 2012). The energy and the corresponding chemoautotrophic yields are higher for NH + 4 oxidation than for NO − 2 oxidation, with an average reported molar C fix :N oxidized ratios of 0.073 for NH + 4 -oxidizing archaea (Zhang et al., 2020) and 0.022 for NO − 2oxidizing bacteria (Zhang et al., 2020). These chemoautotrophic yields were used to constrain the production of sponge holobiont "Spo" (r7) as a function of r5 and r6 (Table 2 and Figure 1). In the anoxic compartments of the sponge holobiont, anammox takes place (Hoffmann et al., 2009), in which N 2 is produced from NO − 2 and NH + 4 by anammox bacteria (r8). The growth and the energy yields of anammox bacteria are well studied; anammox bacteria fix inorganic C coupled to NO − 2 oxidation (r9), with ∼0.07 C fix :NH + 4 oxidized (Prosser, 2005;Thamdrup, 2012), which was used to constrain r9 as a function of r8 ( Table 2 and Figure 1).

Model Constraints
The metabolic network model is generic and can be applied to any sponge volume, depending on data availability. Since two almost complete metabolic flux datasets were available for G. barretti that covered distinct volumes [9-210 ml (Bart et al., 2020b) versus 150-3,500 ml (Leys et al., 2018)] and are based on different methods, we quantified metabolism of "smaller" (from incubation chambers) and "larger" (by in-ex) G. barretti individuals separately (Table 1). This not only allowed us to evaluate the potential size dependence of metabolic activities but also provided an independent assessment of the eventual metabolic patterns inferred.
In addition, the model was constrained with 13 inequalities (soft constraints), of which 10 stated that all metabolic reactions (r1-r10) must be positive. Furthermore, two inequalities limited DOM consumption in denitrification and DNRA (r3, r4) to be lower than aerobic DOM consumption (r2). Lastly, the anammox rates (r8) were constrained to be limited to a maximum of 20% of denitrification rates (r3) (Hoffmann et al., 2009;Rooks et al., 2020).

Model Implementation and Solutions
The metabolic network model of G. barretti was implemented and solved in the modeling environment R (R Core Team, 2018) with R-package LIM (van Oevelen et al., 2010). With 17 unknowns (reactions) and 16 knowns (eight mass balances and eight constraints), the model was almost even determined, resulting in small ranges of possible solutions for each reaction. The (small) ranges in reaction rates were assessed by Bayesian sampling of the solution space with 500 iterations using the function "xsample" (Van den Meersche et al., 2009). This resulted in reaction rate values with standard deviations < 4% of average. The variation in reaction rates is mostly depending on model input parameters, which comprises reaction coefficients (e.g., C:N DOM and f p_N DOM ) ( Table 2) and constraints (e.g., measured O 2 uptake and NO − 3 release; Table 1). However, the linear setup allows only single-parameter values rather than ranges. Inclusion of uncertainty and variability in parameters was achieved by creating one million model input files with values that were randomly sampled from parameter ranges. The ranges for reaction coefficients are as explained in "Section "Reactions and Coefficients"" ( Table 2), and the ranges for constraints are mean ± standard deviation (SD) for larger and smaller sponges ( Table 1). Only ∼0.2-1.0% of parameter combinations resulted in a feasible model solution, corresponding to 10,088 and 2,228 feasible model solutions for smaller and larger sponges, respectively. The feasibility of parameters to produce successful results gives an indication on the likelihood of these parameter values. A successful input file is available in the Supporting Information. The R code to run the model is available at Zenodo (doi: 10.5281/zenodo.4139792). The model results from all feasible model solutions, and the likelihood of parameter values in smaller and larger sponges was analyzed for mean, SD, and min and max based on 5 and 95% confidence intervals (CI).

DIN Exchange in Incubation Chambers
All G. barretti showed a net release of NO − 3 , corrected to control incubations in which NO − 3 concentrations remained constant (Supporting Figure 1), of, on average, 0.95 ± 0.79 µmol N cm −3 day −1 (mean ± SD, used throughout text, n = 12, Table 1). NO − 3 release was positively correlated to O 2 consumption (r = 0.82, p < 0.01, n = 12). The NO − 2 and NH + 4 concentrations in the chambers did not significantly change (i.e., the regression slope was not significant) during incubations with sponges and neither in control incubations without sponges (Supporting Figure 1). Although the concentration changes were insignificant, we calculated control-corrected NO − 2 and NH + 4 fluxes to obtain feasible ranges for the metabolic network model analysis. The average (control-corrected) NO − 2 release and NH + 4 release rates by G. barretti in incubation chambers were 0.023 ± 0.069 and 0.061 ± 0.47 µmol N cm −3 day −1 , respectively (n = 12, Table 1).

Size Relationships
The datasets of smaller G. barretti individuals (9-210 ml in incubation chambers) and larger G. barretti individuals (150-3,500 ml with in-ex) were analyzed together to explore for allometric relationships in data. A statistically significant allometric relationship was found between sponge volume and O 2 consumption rates, although only 21% of the variation could be explained with this regression (Figure 2). The allometric FIGURE 2 | Regression between sponge volume (cm −3 ) and measured volume-specific oxygen consumption (µmol O 2 cm −3 day −1 ) in Geodia barretti by Bart et al. (2020b) (blue circles) and Leys et al. (2018) (dark gray squares). Log-log linear models were used to fit all data (dashed black line). The gray zones indicate the 95% confidence interval. The green circle indicates that the two methods (incubation chambers and in-ex) can produce consistent results (despite the large variability). equation for volume-specific O 2 consumption rates (vOcr, µmol cm −3 day −1 ) as a function of volume (V, cm −3 ) for G. barretti is: vOcr = 25.5 · V −0.23 (p < 0.01, R 2 = 0.21, n = 29, Figure 2). No temperature correction was needed since all data were collected at a similar temperature of 8-9 • C. At overlapping sizes, O 2 consumption rates are comparable between the two datasets (green circle, Figure 2). No allometric relations were found for DIN net exchange data (data not shown).

Model Results
The equalities (hard constraints) and inequalities (soft constraints) were internally consistent, as the model was able to resolve an internal metabolic flux network for smaller and larger G. barretti sponges ( Table 3). As imposed with measured metabolic flux data (Table 1), almost all model-produced metabolic rates (µmol cm −3 day −1 ) were higher in smaller compared to larger specimens (Table 3), except the assimilation of bacteria (r1), which was higher in larger sponges (Tables 1, 3). However, the overall relative partitioning of internal C-and N-transforming processes was similar for smaller and larger sponge individuals (Figures 3, 4).

Sponge Host Metabolism
G. barretti metabolism was dominated by heterotrophy, with DOM as primary food source (Table 3 and Figure 3). DOM assimilation and aerobic respiration by the sponge (holobiont) (r2) were the main heterotrophic processes in G. barretti, independent of size (Table 3 and Figure 3). Because this process dominated sponge metabolism, it was also most sensitive to Release is indicated with a negative sign.
input parameters. Lower values of C:N DOM input range ( Table 2) resulted in more feasible model solutions ( Figure 5). Feasible C:N DOM values averaged 7.6 ± 1.4 in both smaller and larger sponges ( Figure 5). Higher values of f p_N DOM from the rather wide input range (20-80%, Table 2) resulted in more feasible model solutions in smaller sponges, while the opposite was found in larger sponges (Figure 5). The averages of feasible f p_N DOM values were 58 ± 18% in smaller sponges compared to 41 ± 17% in larger sponges (Figures 5, 6). The corresponding average C production efficiency (f p_C DOM ) was 33 ± 12% in smaller sponges and 24 ± 11% in larger sponges (Figures 5, 6). The associated C (biomass) production rate was 4.8 ± 2.8 and 1.0 ± 0.55 µmol C cm −3 day −1 in smaller and larger sponges, respectively (Figures 3, 6). The majority of assimilated DOC, i.e., 67 ± 12% in smaller sponges and 76 ± 11% in larger sponges, was aerobically respired to CO 2 (Figure 6).
The model-calculated exchange rates of O 2 and NO − 3 comprised a smaller range compared to the measured rates (compare Tables 1, 3) because not all combinations of flux measurements resulted in feasible model solutions. The average model-based NO − 3 release rates were lower compared to the averages of the measured rates because more feasible solutions were obtained at the low end of measured NO − 3 release rates (compare Tables 1, 3). Higher O 2 uptake values from the measured range resulted in more feasible solutions in smaller sponges, while the opposite occurred in large sponges (compare Tables 1, 3).

DISCUSSION
The aim of this study was to infer internal C and N metabolic conversions in a deep-sea sponge holobiont system. To this end, a metabolic network model integrating C and N metabolism of the sponge G. barretti and its microbial symbionts was developed. The presented model can serve as a valuable data analysis tool to quantify internal and intermediate routes in sponge metabolism given any metabolic dataset. In this study, we used two independent metabolic datasets to constrain internal G. barretti metabolism that encompassed a range of sponge volumes. The model results indicate that G. barretti has complex but flexible metabolism consisting of aerobic and anaerobic processes. The measured fluxes and model metabolic results suggest that specific metabolic rates decline with increasing G. barretti size, while the ratios between oxic and anoxic processes and between sponge host and microbial metabolism seem rather independent of size.

Organic Matter Assimilation by G. barretti
The measured DOC assimilation rates (Bart et al., 2020b) as well as the estimated DOC assimilation rates (i.e., not directly measured but derived from TOC; Leys et al., 2018) show that DOC is the dominant C source ( > 90%) that fuels G. barretti metabolism. Our metabolic network analysis showed that assimilated DOM serves multiple purposes, since it is used as energy source for aerobic and anaerobic respiration and as a substrate to produce sponge (holobiont) biomass (Table 3 and Figures 3, 6). While the DOC assimilation rates were available to constrain the model, the stoichiometry (molar C:N ratio) of assimilated DOM was not known a priori and was inferred by uncertainty analysis (Figure 5). DOM assimilation was feasible when DOM had a C:N < 9 (Figure 5). This indicates a preference for N-rich DOM because the reported bulk C:N ratios for DOM range from 9 for fresh (labile) to 17 for more refractory material (Hopkinson and Vallino, 2005). Sponges are recently shown to be more selective feeders than traditionally assumed on both particulate (e.g., plankton) (Maldonado et al., 2010;McMurray et al., 2016) and dissolved organic food (Fiore et al., 2017). Selective assimilation (and retention) of N-rich organic matter by G. barretti and other deep-sea sponges has been demonstrated in tracer assimilation studies using dualisotope-labeled substrates (Kazanidis et al., 2018;Bart et al., 2020a). Most marine heterotrophs preferentially assimilate N, as they have C:N ratios which are lower than the C:N ratios of available food sources (Steinberg and Landry, 2017). The average C:N ratio of 4.2 for G. barretti tissue matches well with the C:N ratios of sponges and other deep-sea North Atlantic benthos that had an average C:N of 3.9 (Parzanini et al., 2018) and with those of marine bacteria (C:N = 4-6; Table 2). Whether the preferential assimilation of N-rich DOM is opportunity driven or an actual metabolic requirement of G. barretti metabolism remains to be determined. Because sponges filter water for DOM, they will increase the C:N ratio of surrounding seawater DOM as they selectively retain N, which might influence its bioavailability for surrounding seawater bacteria.
As the C and N production efficiencies for G. barretti and other deep-sea sponges were largely unknown, they were inferred by sensitivity analysis (Figure 5). Our estimated C production efficiency values of 24 ± 8% in larger individuals and 31 ± 14% in smaller individuals (Table 3 and Figures 3,  5, 6) are within the reported range of 20-30% for metazoans (Sterner and Elser, 2002) and ∼30% in pelagic marine FIGURE 6 | Schematic representation of internal C and N metabolism in Geodia barretti. Metabolic processes are indicated with solid arrows with reaction (r) numbers corresponding to those of Figure 2 and Table 2. Aerobic and anaerobic processes are separated in the blue and yellow part, respectively. Exchange fluxes are shown with dotted arrows. Model-calculated average metabolic rates in µmol C cm −3 day −1 or µmol N cm −3 day −1 , shown adjacent to arrows, differ between smaller and larger G. barretti (right). The partition of dissolved organic matter to oxic and anoxic processes and to (C and N) sponge biomass production and respiration are shown with percentages. Dashed arrows indicate chemoautotrophic processes. Internal NH 4 + excretion and re-assimilation, shown by the green arrow, is expected to take place but is not explicitly modeled. Bacteria uptake (r1) and NO 2 − (r15) and NH 4 + (r13) exchange are not shown for simplicity. microbes in the same temperature range (Rivkin and Legendre, 2001). High production efficiencies for G. barretti (∼80%) have been reported in isotope tracer studies with labeled diatom DOM and marine bacteria (Maier et al., 2020;Bart et al., 2020b). In contrast, the encrusting deep-sea sponge Hymedesmia coricea had a maximum production efficiency of 10% when fed with stable-isotope-labeled algae and bacteria for 10 days . Although isotope tracer experiments are advantageous to identify and quantify certain metabolic routes, the obtained values might depend on the incubation time and are representative for the added isotopically labeled substrates only, making it hard to compare the uptake efficiencies to natural diet C and N processing rates. Compared to C production efficiencies, very little information is available on N production efficiencies, especially in sponges. Our estimated N production efficiencies (57 ± 14% in smaller G. barretti and 45 ± 11% in larger G. barretti) ( Table 3 and Figures 5, 6) are comparable to (larger sponges) or above (smaller sponges) those of 40-50% reported for marine zooplankton and bacteria (Touratier et al., 1999). The presence of genes for NH + 4 assimilation through glutamine metabolism in G. barretti (Gavriilidou, personal communication) and several other sponges (Hentschel et al., 2012) is a possible explanation for the slightly higher N production efficiencies in G. barretti compared to those in other organisms. Overall, our estimated values for the C:N ratio of assimilated DOM and the C and N production efficiencies agree with previous estimates for marine organisms, indicating that our model produces reasonable results and might serve as a valuable tool to estimate the efficiencies by which assimilated food is turned into production.

Microbial Nitrogen-Transforming Processes
In contrast to most marine animals, several sponge species (in particular, HMA demosponges) release NO − 3 rather than NH + 4 (Jiménez and Ribes, 2007;Keesing et al., 2013). NO − 3 release also seems a consistent feature of G. barretti as observed in incubation studies with intact G. barretti (Bart et al., 2020b; also presented in this study), explants (Strand et al., 2017;Fang et al., 2018), and tissue fragments (Hoffmann et al., 2009) and by direct inex measurements in intact G. barretti (Leys et al., 2018). The average measured NO − 3 release by intact G. barretti in incubation chambers (0.95 ± 0.79 µmol N cm −3 day −1 ; Table 1) was at the high end compared to average values of 0.20-0.65 µmol N cm −3 day −1 in explants and fragments (Radax et al., 2012;Strand et al., 2017;Fang et al., 2018). A plausible explanation for the higher NO − 3 release in intact sponges is that the functional integration between the host sponge and the microbiome has not been disrupted and that active pumping (not possible within tissue fragments) results in better tissue oxygenation. Indeed a positive relation between O 2 consumption and NO − 3 release in incubation chambers was found to be in agreement with Fang et al. (2018). The NH + 4 and NO − 2 exchange rates were much smaller compared to the NO − 3 release rates as observed by direct in-ex measurements and incubation experiments with G. barretti (this study; Fang et al., 2018) (Table 1). This implies that most substrate for nitrification is produced internally either via sponge cell excretion of NH + 4 as metabolic waste or via microbial DNRA (Figures 2, 6).
Genes (napA, nrfA) involved in all steps of DNRA have been found in metagenomes of G. barretti (Gavriilidou, personal These rates are derived from model results (µmol C cm −3 day −1 ) per carbon biomass (mmol cm −3 day −1 ).
communication) and other deep-sea sponges (Li et al., 2014), but the DNRA process rates in Geodia (or other sponges) have not been quantified yet. The model-based DNRA rates were higher than the model-based denitrification rates (Table 3 and Figures 4, 6). The high organic C concentrations compared to NO − 3 concentrations would indeed favor DNRA (van den Berg et al., 2016), but whether these conditions are found in G. barretti remains to be validated. Studies on DNRA and denitrification gene expression might help to understand the relative role of each process, and that information can be used to constrain the rates in the model. It is however likely that DNRA is an important nitrogen-transforming process in G. barretti. In an earlier version of the model without DNRA, more NH + 4 had to be produced by excretion, which was balanced by the assimilation of DOM with an even lower C:N, lower nitrification rates, and NO − 3 release. Dark carbon fixation rates in G. barretti (or other Geodiidae spp.) associated with nitrification and anammox have not yet been experimentally quantified. The first model-based estimates of dark carbon fixation rates for G. barretti are presented here ( Table 3 and Figure 3), acknowledging the different energy yields from each nitrogen-transforming process (e.g., NH + 4 oxidation, NO − 2 oxidation, anammox). CO 2 fixation rates contributed only a small fraction of 5.1-6.5% of G. barretti production and ∼1.5% of total C assimilation. These contributions are very similar to the range of 0.2-2.1% fixation relative to assimilation for the deep-sea encrusting sponge Hymedesmia coriacea (van Duyl et al., 2020). The benefit of hosting nitrifiers by the sponge holobiont is hypothesized to decrease or rapidly inactivate the internal concentration of NH + 4 produced as metabolic waste, which is possibly toxic to the host cells. In addition, it provides the substrate (NO − 3 ) for heterotrophic symbiont production in anoxic parts of the sponges. Combined with the production of DNRA and denitrification, N-transforming processes enhance the metabolic capacity of G. barretti and contribute significantly to G. barretti production [28 ± 16% (smaller individuals) and 32 ± 16% (larger individuals) of total production] ( Table 3 and Figure 3). The substantial production by endosymbionts supports the idea that G. barretti, which is a HMA sponge, uses "microbial farming" and phagocytes and digests part of the microbiome as part of heterotrophic feeding.
The model covers N-based metabolism but not sulfurbased metabolism because sulfur data were not available. Both genes for sulfate reduction (heterotrophy) and sulfur oxidation (chemoautotrophy) are found in G. barretti (Jensen et al., 2017). The contribution of microbial production to G. barretti might become even higher if these processes would be included in the metabolic network model.

Size Dependency of Metabolism
Based on the metabolic datasets of Leys et al. (2018) and Bart et al. (2020b), all measured and modeled volume-specific rates, except the assimilation of external bacteria, were lower in larger than smaller sponges (Tables 1, 3). The measured O 2 consumption by G. barretti showed large variability within and across the two metabolic datasets, which can be partly attributed (21% based on R 2 ) to sponge volume (Figure 2; Leys et al., 2018;Bart et al., 2020b). The large variability within each metabolic dataset (Figure 2) might be linked to the intrinsic variability between individuals, e.g., the physiological and metabolic state of the holobiont (and analytical uncertainty), rather than external conditions as individuals were measured under the same experimental conditions. We cannot exclude that differences across the two metabolic datasets are partially caused by methodology [non-invasive in-ex, but snapshot (minutes) sampling vs. incubations, longer time integration (hours), but more prone to incubation effects]. However, comparable estimates of O 2 consumption rates were obtained in sponges of similar size (green circle, Figure 2), indicating that the two complementary methods give consistent results.
O 2 consumption in G. barretti is directly linked to pumping since pumping supplies sponge tissue with O 2 and G. barretti that reduce their pumping activity rapidly become anoxic (Hoffmann et al., 2005;Leys et al., 2018). Lower specific pumping and clearance rates in larger individuals have been shown for some massive (i.e., non-encrusting) sponge species (Reiswig, 1974;Morganti et al., 2019), including G. barretti (Leys et al., 2018). Reduced pumping and O 2 uptake in larger-sized ball-shaped sponges might be an effect of longer flow paths in larger individuals, resulting in increasing constraints to oxygenate the whole tissue. Lower O 2 consumption in larger individuals did not result in a shift from aerobic to anaerobic processes. This is surprising because we had presumed that larger sponges would run easily into anoxic conditions at the deepest and less water-irrigated portions of their mesohyl. In contrast, the results suggest that all mesohyl regions are similarly, well irrigated by the network of aquiferous canals. However, it could be that any metabolic differences were masked by the large uncertainty associated with aerobic and anaerobic N-transforming processes (Figure 4). It will be interesting to compare microbial community structure between smaller-and larger-sized G. barretti to see how that changes with size.
The lower production efficiencies (24 ± 8% for C, 45 ± 11% for N; Table 3 and Figure 5) and production rates (0.067% day −1 ; Table 4) in larger G. barretti compared to smaller G. barretti (31 ± 14% for C, 57 ± 12% for N; Table 3 and Figure 5; 0.19% day −1 , Table 4) fit with theoretical predictions that production efficiencies decrease with size (Banse, 1979;Sterner and Elser, 2002) because organisms cannot grow infinitely. Larger individuals will relatively use more energy for maintenance, while smaller individuals will allocate more energy to growth. Reiswig (1973) measured oscula growth in tubular sponge populations and observed that smaller individuals grow faster than larger individuals of the same sponge species and that a sharp decline in growth rate occurred when sponges reached maturity and started reproducing, thus allocating energy to gametes.

Sponge (Biomass) Production
Production by the sponge holobiont is allocated to various processes, which are mainly growth, reproduction, and renewal. The magnitudes and portions of each flow are largely unquantified for Geodia (and other deep-sea sponges). Therefore, we did not go any further with our model beyond estimating production by the sponge holobiont. Our production rates for G. barretti of 0.16% day −1 in smaller sponges and 0.067% day −1 in larger sponges (Table 4) seem low compared to modeled estimates for deep-sea sponges (0.23% day −1 ; Gontikaki et al., 2011), but not unrealistic as Geodiidae are generally considered to be slow growing. The yearly production would be 0.15 year −1 for larger sponges and 0.57 year −1 for smaller sponges, although temporal upscaling comes with large uncertainty as sponge metabolism, related to environmental conditions, might vary seasonally (Koopmans and Wijffels, 2008;Morley et al., 2016). Several sponge species (e.g., de Goeij et al., 2013;Alexander et al., 2014;Maldonado, 2016), including deep-sea species (Witte et al., 1997;Rix et al., 2016), release a substantial (up to 10% day −1 ) part of biomass as sponge detritus, but other studies show no or only limited production by (mostly non-encrusting) species (e.g., Yahel et al., 2007;Maldonado, 2016;McMurray et al., 2018;Pawlik and McMurray, 2020). Sponge detritus production by G. barretti was found to be low (Leys et al., 2018;Maier et al., 2020). Maier et al. (2020) estimated detritus release by G. barretti from tracer experiments as 0.03% day −1 , which would still be 19-44% of our estimated sponge production. However, because of their efficiency in turning metabolic waste into new products, it could also be that any organic waste gets recycled within the holobiont rather than excreted. The low production rates, from which only an unconstrained fraction will be used for net growth, can explain why this sponge species is typically considered slow growing.

Ecosystem Context Using Model Estimations
The metabolic activity of G. barretti influences local C and N cycling at the benthic boundary layer, given the enormous abundance of this species over the extensive deep-sea area of the boreal region, forming the so-called Geodia grounds (Pham et al., 2019). G. barretti and other Geodiidae can make up > 90% of total benthic biomass (Klitgaard and Tendal, 2004;Murillo et al., 2012) in these regions at up to 6 kg WW m −2 (Klitgaard and Tendal, 2004;Kutti et al., 2013), although their distribution is patchy (Pham et al., 2019). Our sponge metabolic model results may aid in the first estimates of biogeochemical fluxes of C and N in these areas, specifically when only a subset of fluxes or environmental settings is known (e.g., only O 2 fluxes, no C or N). For example, Geodia grounds have been identified as hotspots of metabolic activity compared to bare sediments (Kutti et al., 2013;Cathalot et al., 2015), showing one order of magnitude higher O 2 consumption at high-abundance sponge grounds than for bare sediments (54 versus 6 mmol O 2 m −2 day −1 ; based on 3 kg G. barretti WW m −2 ; Cathalot et al., 2015). Using our modelestimated RQ and C:N conversions, future studies can now approximate C and N fluxes based on in situ-assessed O 2 fluxes alone. For example, G. barretti enhances nutrient recycling via DIN (NO − 3 ) release (0.4-1.1 mmol N m −2 day −1 for a biomass of 1-3 kg WW m −2 ), which is comparable to or in excess of values typical for marine sediments [0.06 mmol N m −2 day −1 in the north Norwegian coast (Glud et al., 1998) to 0.64 mmol N m −2 day −1 in Logachev mound province, North-Atlantic (de Froe et al., 2019)]. However, sponge grounds are also a sink for N, as they release N 2 at rates of 0.1-0.3 mmol N m −2 day −1 , comparable to fluxes found in shelf sediments (0.1-1.0 mmol N m −2 day −1 ; Middelburg et al., 1996;Seitzinger and Giblin, 1996).
Care must be taken, however, since, on a larger ecosystem scale, the average G. barretti biomass in the south-west Barents Sea region (i.e., Trømso Plateau, North Cap bank, and areas in between; approximately 125,000 km 2 ) is ∼50 g WW m −2 (Anisimova et al., 2010;Kêdra et al., 2017) and is even lower within the entire Barents Sea (∼15 g WW sponge m −2 , 1.5 million km 2 ; Wassmann et al., 2006). This points to ordersof-magnitude differences in local (see above) versus ecosystemscale biogeochemical fluxes driven by sponges, including many species other than G. barretti from different phyla. In addition to (largely unknown) variations in metabolic rates between sponge species, the effect of metabolic rates under varying environmental conditions, including seasonality or hydrodynamics, is not yet well understood, considering the two available datasets that we used for our model.
Another important aspect of our model that can improve the ecological context of deep-sea sponges and sponge grounds is the size dependence on metabolic rates. Currently, benthic biomass estimations in deep-sea ecosystems, such as the Barents Sea, are largely based on trawl catch data (e.g., Wassmann et al., 2006;Jørgensen et al., 2014). Trawling is an invasive and only semiquantitative (but cost-effective) method to assess benthic biomass and is generally skewed to larger individuals. The abundance and size distributions of deep-sea sponges can be better quantified on a larger scale by the non-invasive method of video transects (e.g., Kutti et al., 2013;Kêdra et al., 2017;Kazanidis et al., 2019;Maldonado et al., 2020b). Those data, combined with our modeled biogeochemical rates and efficiencies for different size classes of sponges, will further allow (first estimates of) quantifications of biogeochemical fluxes in sponge grounds of the NorthAtlantic Ocean and assess the functional integration of these sponge-dominated communities into deep-sea ecosystems.

DATA AVAILABILITY STATEMENT
The dataset of Leys et al. (2018) is available at University of Alberta Education and Research Archive: doi: 10.7939/ R3057D48V. The dataset of Bart et al. (2020b) will be made available at Pangaea. The model input files and R scripts are available on Zenodo: doi: 10.5281/zenodo.4139792.

AUTHOR CONTRIBUTIONS
AK developed the model, performed the data analyses and model analyses, prepared the figures, and wrote the manuscript. MB and JG designed and conducted the incubation experiments, performed nutrient analyses, and contributed to the writing of the manuscript. SV contributed to model conceptualization and development. JM acquired funding and contributed to the conceptualization of the model, interpretation, and writing of the manuscript. MM contributed to the conceptualization of the model, interpretation, and writing. SL and SM contributed data for model development. KS contributed technically with model implementation and analyses. DO designed the model, contributed to model development and interpretation, and writing of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This research has been performed in the scope of the EU SponGES project, which received funding from the European Union's Horizon 2020 Research and Innovation Program under grant agreement no. 679849. Further support included ERC starting grant agreement no. 715513 to JG and the Netherlands Earth System Science Center to JM. This document reflects only the authors' views, and the Executive Agency for Small and Medium-sized Enterprises (EASME) is not responsible for any use that may be made of the information it contains.

ACKNOWLEDGMENTS
Titus Rombouts and Pieter Slot (UvA) and Sharyn Ossebaar (NIOZ) are acknowledged for their analytical assistance with the nutrient measurements of the incubation experiments. We thank Asimenia Gavriilidou and Detmer Sipkema (WUR, SponGES) for their input on the microbial genome in G. barretti. We thank last Hans Tore Rapp (UiB) for the excellent project coordination. We would like to thank CF and CR for their valuable and constructive reviews.