In situ Pumping Rate of 20 Marine Demosponges Is a Function of Osculum Area

Sponges play a key role in the transfer of energy and nutrients into many benthic ecosystems, and the volume of water they process is an important regulator of these fluxes. Theoretical scaling relationships between sponge volume, osculum cross-sectional area, and pumping rates were recently proposed and confirmed for small sponge specimens in the lab. To examine how these relationships apply to field populations we measured, in situ, the pumping rate (PR) of 20 species representative of different morphologies and host types (high- and low-microbial-abundance, HMA and LMA) from temperate and tropical regions. The total oscula area (∑OSA) increased allometrically with sponge volume (V) exhibiting similar exponents (∑OSA=aVb, b ranging 0.6–0.7) for all species, except for tropical HMAs (b = 0.99). Osculum flow rate (OFR) also increased allometrically with OSA and oscula of the same size pumped at the same rate irrespective of sponge volume. As a result, and in contrast to former reports, the PR of most of the sponges increased allometrically (PR=a∑OSAb) with scaling exponent b≈0.75, whereas PR of tropical HMAs increased isometrically. Osculum jet speed declined with the increase in the OSA for most species. The number of oscula and their OSA were the best predictors of the PR in sponges, explaining 75–94% of the in situ variation in PR throughout the natural range of sponge size. The pumping rate of a sponge population can be estimated by measuring the osculum density and cross-sectional area distribution once the relationships between the OSA and OFR are established for each species.


INTRODUCTION
Sponges are widespread across different ocean regions from intertidal to abyssal depths (Van Soest et al., 2020). In some areas they can constitute ∼90% of the benthic biomass, excluding benthic fishes (Klitgaard and Tendal, 2004;Murillo et al., 2012). Such dense sponge aggregations can filter hundreds of cubic meters of water per m 2 daily, with significant consequences to the benthic-pelagic coupling and local biogeochemical cycles (Kutti et al., 2013;Kahn et al., 2015). Their role in organic matter cycling has been equated to that of the microbial loop because some sponges can remove dissolved organic matter from the water column and make it available in the form of detritus to higher trophic levels (de Goeij et al., 2013;Rix et al., 2016aRix et al., ,b, 2017Rix et al., , 2020. There is growing evidence of the key role that sponges play in nutrient cycling in shallow and deep-water communities (Maldonado et al., 2016;Pawlik and McMurray, 2019;Zhang et al., 2019;Folkers and Rombouts, 2020). However, the magnitude of these processes is mainly driven by the amount of water they pump. Therefore, an efficient method to measure the volume of water processed by a sponge population is crucial to properly determine the relevance of the processes mediated by sponges in the ecosystems.
As filter feeding organisms, sponges rely on their ability to daily process copious amounts of water through their body [up to 35 mL min −1 (cm sponge) −3 , Weisz et al., 2008] to obtain the required oxygen and nutrients and excrete the metabolic waste products (Hadas et al., 2008;Morganti et al., 2017). The body of most demosponges, a class that includes >75% of the extant sponge species (>8,000 species according to the World Porifera Database, Van Soest et al., 2020) is entirely specialized for suspension-feeding through a complex and highly vascularized aquiferous system. This efficient fluid transport system is composed of choanocyte chambers (pumping units), ostia (incurrent apertures) and a complex network of incurrent and excurrent canals. The latter converge into an excurrent cavity (sometimes referred as atrium) from which the water is expelled through exhalant apertures called oscula (Simpson, 1984). The water current is generated by the synchronized beating action of choanocyte flagella gathered in the choanocyte chambers (Asadzadeh et al., 2019), and the volume of water pumped is positively correlated to their density (Massaro et al., 2012). After passing the pumping units, where the capture of particles occurs, the water is ejected via exhalant canals through a single exhalant aperture (osculum). The choanocyte chambers and inhalant and exhalant canals associated with one single osculum are considered an "aquiferous module" (Fry, 1970(Fry, , 1979Ereskovskii, 2003). As the sponge increases in size, the number and size of the aquiferous modules increase (Wiens et al., 2006). Consequently, their aquiferous system may become more complex (Gaino et al., 1995;Hammel et al., 2012) as reflected by the increase of number of oscula and their cross-sectional area in larger specimens (Gili et al., 1984;Dahihande and Thakur, 2019). Quantitative description of the aquiferous module and the relationships between the aquiferous module and sponge size are yet to be provided. Only few studies have explored the allometric scaling of the increase in the number of oscula and the osculum cross-sectional area and its consequences for sponge pumping (Dahihande and Thakur, 2019;Goldstein et al., 2019;Kealy et al., 2019).
Sponge pumping behavior has been documented to vary between individuals and species (Reiswig, 1971). Environmental factors such as temperature, food concentration, suspended sediment concentration, water flow, and viscosity have been observed to affect sponge pumping rates in laboratory and in situ experiments (Reiswig, 1971;Vogel, 1974;Frost, 1980;Riisgård et al., 1993;Riisgård and Larsen, 1995;Tompkins-MacDonald and Leys, 2008;Grant et al., 2018Grant et al., , 2019Dahihande and Thakur, 2019). Intrinsic factors, such as osculum cross-sectional area, microbial abundance, choanocyte density, reproductive stage, and sponge volume have also been shown to affect sponge pumping rates Massaro et al., 2012;McMurray et al., 2014;Strehlow et al., 2016;Dahihande and Thakur, 2019;Goldstein et al., 2019;Morganti et al., 2019). Pumping rate has been observed to increase with sponge volume in several studies (Thomassen and Riisgård, 1995;Kowalke, 2000;McMurray et al., 2014;Lewis and Finelli, 2015;Morganti et al., 2019), indeed the latter has been shown to be the major determinant of pumping rate. Volume-specific pumping rate has been documented to exhibit a pattern of decrease with sponge volume, but the rate of this decrease varied between sponges (Morganti et al., 2019). However, sponge volume might be a difficult parameter to estimate for some species and communities (e.g., excavating sponge species, Rosell and Uriz, 2002). Direct relationship between the osculum cross-sectional area and the amount of water pumped by a sponge has been recently demonstrated in laboratory experiments with small sponge explants (Strehlow et al., 2016;Kumala et al., 2017;Goldstein et al., 2019;Kealy et al., 2019) and in one in situ study (Gokalp et al., 2020), suggesting that osculum cross-sectional area and the number of oscula may be more practical predictors of sponge pumping. However, the relationships between sponge structure and function and, specifically between sponge pumping rates and the oscula number and their cross-sectional area are yet unresolved.
Our objective was to establish the relationships between the size structure (morphology, number of oscula, and osculum cross-sectional area) of marine sponges and their in situ pumping rates and to compare the observed relationships with the theoretical scaling recently suggested by Goldstein et al. (2019) and Kealy et al. (2019) in order to improve our understanding of the factors controlling this key process, and to allow accurate and efficient quantification of sponge pumping rate in nature. To that end, we examined the variation in the number of oscula and osculum cross-sectional area with sponge volume, and the relationship of these parameters with osculum flow rate and sponge pumping rate for a variety of dominant sponge species in three regions covering their natural size spectrum. In each region, we examined both high-and low-microbial-abundance sponges (HMA and LMA, respectively). Finally, we provide a quantitative comparison of different approaches to estimate the pumping rate of sponge populations in situ.

Study Sites
The study was conducted in situ in three different areas representing a wide range of depths and natural habitats:
All examined species belonged to the Demospongiae class and species from the Order Haplosclerida were present in the three regions. Species from the Orders Agelasida and Dyctyoceratida were present in the Mediterranean and Caribbean Sea and, species from the Order Poecilosclerida were present in the Mediterranean and Red Sea (Supplementary Table 1).
The study included three main different morphologies from the examined regions: 8 massive species from all regions, 4 encrusting species from the Mediterranean and Red Sea, 7 tubular species from the Red and Caribbean Sea, and a branching species from the Red Sea (Supplementary Figure 1,  Supplementary Table 2). HMA and LMA species were present in all regions: there were 3 HMA species from the Mediterranean, 3 from the Red Sea and 5 from the Caribbean, and 2 LMA species from the Mediterranean, 4 from the Red Sea and 3 from the Caribbean. HMA and LMA species occurred in the main morphologies, but there were no HMA encrusting species: there were 7 massive, 3 tubular and 1 branching HMA species, and 1 massive, 4 tubular and 4 encrusting LMA species (Supplementary Table 2).

Sponges Pumping Rate Measurements
Pumping rate measurements were conducted on visually healthy individuals over a broad size range of sponge species representative of the population in the study area. Pumping rate was measured using the dye front speed method (DFS) described by Yahel et al. (2005). To avoid deviations from ambient water density, the sodium fluorescein powder was loaded in the syringe prior to the dive. Underwater, the dye powder was mixed with ambient water drawn into the syringe next to the sponges and a disposable syringe filter (25 mm, 0.2 µm) was installed on the syringe to avoid the release of dye particles. A small amount of dyed seawater from the syringe was placed at the lower part of the tube (∼ the first centimeter) while covering the tube's upper end with the index finger (Supplementary Video 1). The transparent tube was positioned as close as possible above the sponge osculum, the index finger was removed from the tube's upper end and the movement of the dye inside the tube was recorded by a second diver using a video camera. This procedure was repeated 6-10 times for each measurement. A frame-byframe analysis was used to measure the speed of the dye front inside the tube in each repetition and an average dye front speed was calculated for all good repetitions (3-10). The tube internal diameter was selected to be as close as possible but slightly larger than the diameter of the osculum (no larger than 40% of the osculum diameter, Yahel et al., 2005) by using a set of 12.5 cm long glass tubes with internal diameter increments of 5 or 6 mm. The rate of water flow from the osculum was calculated following Yahel et al. (2005) as the product of tube cross-sectional area and the dye front speed, or, in the few cases where the tube was smaller than the osculum (∼16%), as the product of dye front speed and osculum cross-sectional area. Due to the large size of the oscula from the Caribbean Sea sponges, another method was used to estimate the pumping rate of the Caribbean species . The excurrent water velocities were determined by videotaping the movement of dye fronts in the excurrent plume. A ruler was positioned directly behind and parallel to the excurrent stream. The video also established the area of the excurrent plume, which was not always equal to the crosssectional area of the osculum. One diver recorded the movement of the dye with the video camera as a second diver released small puffs of a concentrated fluorescein dye solution into the excurrent stream in front of the ruler at the level of the oscular opening. Pumping rate of the Caribbean species was calculated as the product of the dye front speed in the excurrent plume and the area of the excurrent plume as detailed in Weisz et al. (2008). Analysis of three replicates of DFS measurements per osculum indicated an accuracy of 9 to 12% for the three methods.
For each species, we sampled at least 39 specimens from the Mediterranean study site, whereas fewer specimens were sampled at the Red Sea and the Caribbean's (2-6 and 2-16, respectively). While at the Mediterranean and Caribbean Sea each specimen was sampled only once due to logistic constraints, at the Red Sea the pumping rate of each specimen was measured at least twice during the day and at least once during the night. Such repetitive measurements allowed to assess the consistency of the pumping behavior and to detect differences along the daily cycle (see Moskovich, 2020 for more details). Nocturnal dye front speed measurements were made with violet-UV light as depicted in Supplementary Video 1 to minimize the disturbance to the studied sponges.
All sampled species from the Caribbean Sea had only one osculum. The sampled specimens from the Red Sea were selected to have a small number of oscula (1-3; except from a single specimen of Suberites clavatus that possessed 12 oscula). Due to the higher number of oscula and logistic constraints to oscula accessibility, for the multi-osculated species of the Mediterranean Sea the oscula were divided into three size classes based on the maximum osculum diameter of each species. Wherever possible, we sampled at least one osculum representative of each size class in each sampled specimen. The osculum diameter size classes were as follow: D. avara, C. crambe and A. oroides: small ≤ 2 mm, 2 > medium ≤ 4 mm, large > 4 mm. P. ficiformis: small ≤ 2 mm, 2 > medium ≤ 3 mm, large > 3 mm.
The volume of water processed by the whole sponge (hereafter pumping rate, PR, mL min −1 sponge −1 ) was calculated for each specimen as the sum of the average osculum flow rate (OFR, mL min −1 ) of each osculum class within that specimen, multiplied by the number of oscula in the respective size class within each specimen, using the following equation: Where i denotes the osculum size classes (large = 1, medium = 2, small = 3) within the sponge specimen, n i is the number of oscula from the ith osculum size class in the specimen, and OFR i is the mean osculum flow rate (mL min −1 ) of the ith osculum size class within the specimen. Due to logistic constraints, such as oscula accessibility, we could not sample all oscula from each Mediterranean sponge, but an effort was made to sample as many oscula as possible per each specimen in order to have a representative number of oscula per specimen. The sampled oscula were on average (±SD): 39 ± 24% of total osculum number in D. avara, 65 ± 28% in C. crambe, 54 ± 31% in P. ficiformis, 85 ± 20% in C. reniformis, and 84 ± 18% in A. oroides. Note that all oscula of each specimen were sampled in the Caribbean (number of oscula =1) and Red Sea species (number of oscula =1-3, except from a single specimen of Suberites clavatus where 4 oscula were sampled out of 12). Sponge and oscula dimensions were measured for each sponge specimen. Each osculum was photographed individually with an appropriate scale and measured with the software Image J (Ver. 64). Each sponge specimen was photographed vertically and horizontally close to a ruler. In addition, an outline sketch of the horizontal profile of the sponge was drawn in which the different heights were measured and noted. Subsequently, we used manual delineation with Image J to estimate the sponge area from the vertical images. Manual delineation was also used to estimate sponge volume from the horizontal images on the basis of the approximated single or combined geometrical shapes (such as cylinders and cones) that characterized each specimen, the heights of the profile and the area. Estimates of sponge volume generated by this method closely resembled estimates made by 3D video photogrammetry (Moskovich, 2020, Supplementary Figure 2).

Statistical Analyses
We explored the scaling laws between sponge volume (V), osculum cross-sectional area (OSA), number of oscula, osculum jet speed (U 0 ), osculum flow rate (OFR), and pumping rate (PR) using power and linear regressions. A t-test was used to compare whether the scaling exponent b of the allometric functions was significantly different from one. For the Mediterranean species, linear models with Tukey post-hoc analysis were applied to test for differences in OFR (mL min −1 ) across oscula size classes (levels: large, medium, small) by correcting for the sponge volume (V, cm 3 ).
To analyze the relationship between the sponge pumping and the allometric variables (sponge volume, total OSA, number of oscula, and average OSA), multiple and backward stepwise regression analyses were performed. Variables were removed from the model when F < 3.900 and p > 0.055. Variables retained in regression models were screened for collinearity via variance inflation factor (collinear measures with variance inflation factor > 4 were avoided). Variables were ln-transformed to satisfy the normality and/or heteroscedasticity assumptions and a complete residuals analysis was performed to validate the robustness of the resulting model.
Three specimens, two C. crambe (sponge volume > 30 cm 3 ) and one D. avara (sponge volume > 400 cm 3 ) were removed from these analyses because they were out of the specific species population volume range.
The exponents and coefficients of the allometric scaling, multiple and backward regressions were calculated using SigmaPlot 11. The linear model and figures were performed with R studio (version 3.2.1) using the lm function and ggplot2 package, respectively.

Determinants of Sponge Pumping Rate
Mean osculum jet speed (U 0 , cm s −1 ) decreased allometrically with increasing osculum cross-sectional area (OSA, mm 2 ) in all Mediterranean species ( Figure 1A) with negative scaling exponents ranging from −0.39 to −0.76. The U 0 of the smallest oscula we sampled from the different species (OSA ∼2 mm 2 ) ranged from 7.3 to 18.4 cm s −1 whereas, U 0 of the largest oscula sampled was 20 to 40% lower, ranging from 2.4 to 7.1 cm s −1 . For the tropical (Red and Caribbean Sea) species for which sample size was limited, we examined only the trend exhibited by U 0 plotted vs. OSA and tested whether it decreased (b<0) or increased (b>0). U 0 decreased with OSA in all five Mediterranean species, two out of the seven Red Sea species, and five out of the eight Caribbean species. Two tropical sponges showed no trend, potentially due to a low OSA sampled range, and an increase of U 0 with OSA was observed in six of the tropical species. All three main growth forms included species exhibiting both patterns in the relationship of U 0 ∼ OSA b (b < 0: 5 massive, 2 encrusting and 5 tubular; b > 0: 1 massive, 2 encrusting and 2 tubular). Similarly, both host types (HMA and LMA) included species exhibiting both patterns (b < 0: 5 LMA, 7 HMA; b > 0: 3 LMA, 3 HMA) ( Table 4). In all species, the osculum flow rate (OFR, mL min −1 ) increased as a function of OSA ( Figures 1B, 2A,B). This trend was highly significant for all species, but with different allometric parameters for each species. In the Mediterranean species, all scaling exponents (b) were < 1 (0.36-0.90, Figure 1B and Figure 2A with all Mediterranean species plotted together), indicating that the rate of the increase of the OFR with OSA decreases for larger oscula. The sponge species from the Red Sea and the Caribbean Sea were pooled together due to the small number of specimens studied within each species. In the tropical species, the relationships of the OFR with OSA were either isometric (a linear increase of the OFR with the OSA) or allometric, but with scaling exponent b > 1 ( Figure 2B). S. vesparium was the only tropical species for which we had a sufficient number of measurements to be analyzed separately. The OFR of S. vesparium increased as a function of OSA with an isometric, rather than allometric scaling (b was not significantly different from one, t-test p = 0.857, Supplementary Figure 3). Note that the a coefficients of the allometric equations of the different sponges could be markedly different (e.g., Figure 1B) and the smallest oscula we measured in the Mediterranean species (OSA ∼2 mm 2 ) showed distinct flow rates ranging between 6.7 mL min −1 for P. ficiformis to 22.4 mL min −1 for A. oroides.
Since OSA was found to be the best predictor of the OFR, we explored the relationship between the whole organismpumping rate (PR, mL min −1 sponge −1 ) and the sum of the cross-sectional area of all oscula in a sponge ( OSA, mm 2 ). PR of all temperate sponges and the tropical LMAs increased as an allometric function of the total oscula area (PR∼ OSA b ) with scaling exponents (b) significantly smaller than one (t-test p < 0.05 for all sponges but P. ficiformis, see Supplementary Table 3A) and ranging 0.67-0.76, whereas tropical HMAs increased isometrically (b was 1.07 and not significantly different from one, t-test p = 0.320) (Figures 3A,B; and see Supplementary Table 3A for allometric functions calculated for each Mediterranean sponge species separately). However, it should be noted that b = 0.78 for tropical HMAs if S. vesparium is not included in the analysis. FIGURE 2 | Allometric relationship (OFR = aOSA b ) between osculum flow rate (OFR, mL min −1 ) and osculum cross-sectional area (OSA, mm 2 ) for: (A) temperate species from the Mediterranean Sea (5 species), and (B) tropical species from Red Sea (7 species, redish color shades) and Caribbean Sea (8 species, blue color shades). Note that all Caribbean species have a single osculum. The allometric equations are shown in each panel and the shade areas represent the 95% confidence interval for the regression line. The b scaling exponents were statistically different from zero for all equations (p < 0.001). Note the log-log scale.
FIGURE 3 | The relationship between the pumping rate (PR, mL min −1 sponge −1 ) and the sum of the cross-sectional area of all oscula ( OSA, mm 2 ) are plotted using log scale in (A) the Mediterranean species and (B) tropical species from Red Sea (reddish color shades) and Caribbean Sea (blue color shades). For the regression analysis, sponges species were clustered accordingly to host type: HMA (closed symbols, solid regression lines) and LMA (empty symbols, dashed regression lines). The scaling exponents and coefficients calculated for HMA and LMA separately are shown in each panel and the shade areas represent the 95% confidence interval for the regression line. The b scaling exponents were statistically different from zero for all equations (p < 0.001). The scaling exponents and coefficients calculated separately for each Mediterranean species are shown in Supplementary Table 3A. Note different log-log scale between the two panels.
In comparison to PR∼ OSA b , the allometric relationships of pumping rate with sponge volume (PR∼V b ) showed higher variability among regions and host type (b range: 0.  Figures 3B, 4B), but not for tropical LMAs. These results were corroborated by a backward stepwise multiple FIGURE 4 | The relationship between the pumping rate (PR, mL min −1 sponge −1 ) and the sponge volume (V, cm 3 ) are plotted using log scale in (A) the Mediterranean species and (B) tropical species from Red Sea (reddish color shades) and Caribbean Sea (blue color shades). For the regression analysis, sponges species were clustered accordingly to host type: HMA (closed symbols, solid regression lines) and LMA (empty symbols, dashed regression lines). The scaling exponents and coefficients calculated for HMA and LMA separately are shown in each panel and the shade areas represent the 95% confidence interval for the regression line. The b scaling exponents were statistically different from zero for all equations (p < 0.001). The scaling exponents and coefficients calculated separately for each Mediterranean species are shown in Supplementary Table 3B. Note different log-log scale between the two panels.  . These values were calculated as 100 × R 2 of each regression model. Each column represents different predictor parameter: V, sponge volume (cm 3 ); A, sponge area (cm 2 ); OSA, total osculum cross-sectional area (mm 2 ), and the last column is the % variance explained by the multiple regression analysis using number of oscula (#oscula) and mean oscula size (average OSA) as predictor variables. Mean (± the 95% CI) is the average R 2 for all five species. % improvement was defined as the mean (± the 95% CI) increase of the percentage of variance (100 × R 2 ) explained using the different predictors and compared to the percentage of variance explained using V (sponge volume) as a predictor.
regression analysis which showed that removal of the sponge volume from the model did not reduce the model ability to explain the variance in PR in all Mediterranean species except in C. crambe (p < 0.001; Supplementary Table 4). A forward multiple regression analysis has yielded similar results.
The total OSA is a function of the number of oscula and their average cross-sectional area. Using these parameters instead of total OSA further refined the model. The number of oscula and the average OSA within each specimen explained between 75 and 94% of the in situ variations in pumping rates along with the natural sponge size range (multiple regression analysis, p < 0.01; Table 2), representing an improvement of ∼30% in the accuracy of the pumping rate estimates over the use of sponge volume as a predictor in Mediterranean species (Table 1).

Sponge Morphology
Sponge and oscula dimensions were measured for all sampled specimens and reflected the natural population size distribution at the study site (less so where <10 specimens were sampled). The volume of the sampled sponges varied considerably between the study sites (Supplementary Figure 5). Mediterranean species ranged from 0.8 to 440 cm 3 , Red Sea species ranged from 1 to 196 cm 3 , and Caribbean species were larger, with volume ranging from 50 to 32,552 cm 3 (Supplementary Table 2). The sponge volume varied with sponge species and host type: HMA sponges were on average 4 times larger than LMAs in the Mediterranean sponge species (87 ± 9 and 20 ± 3 cm 3 , respectively), 7 times larger in the Red Sea species (94 ± 52 and 13 ± 9 cm 3 , respectively), and 45 times larger in the Caribbean Sea species (7,186 ± 9,147 and 159 ± 109 cm 3 , respectively) (Supplementary Table 2; see also Supplementary Figure 5). In the Mediterranean species the average number of oscula was lower in C. reniformis and A. oroides (4 ± 0.4) and 3-folds higher in P. ficiformis and D. avara (12 ± 2 and 13 ± 2, respectively), and the OSA ranged from 0.1 to 64 mm 2 . In contrast, the specimens analyzed from the Red Sea typically had a low number of oscula (1-3) and the OSA ranged from 8 to Values are the regression coefficients ± the 95% CI for each variable (rows) and each species (columns). The intercept (a), R squared (R 2 ), and variance inflation factors (VIF) are also presented. ***p < 0.001; **p < 0.01; *p < 0.05. Table 2). All Caribbean species were single-osculated with a large OSA that ranged from 79 to 39,584 mm 2 (Supplementary Table 2). Over the range of sponge sizes analyzed in this study, the total OSA increased with sponge volume at a similar rate for both HMA and LMA species from the Mediterranean Sea (scaling exponent ± 95% CI, b (HMA) = 0.56 ± 0.10, b (LMA) =0.60 ± 0.12, Figure 5A). HMA and LMA Red Sea and Caribbean species were analyzed together and showed similar results (scaling exponent ± 95% CI, b (HMA) = 0.99 ± 0.38, b (LMA) =0.71 ± 0.29, Figure 5B). However, for both temperate and tropical species, LMAs showed higher total OSA than HMAs over the entire analyzed size range (Figures 5A,B). The osculum area that is associated with the sponge volume, expressed as the ratio between total OSA and sponge volume (mm 2 cm −3 ), decreased with the increase of sponge volume in all studied species (Figures 5C,D). Osculum area associated with sponge volume was, on average, about 3folds higher in LMAs compared to HMAs for the same size range specimens in temperate Mediterranean species ( Figure 5C) and about 12-folds higher in tropical species (Figure 5D).

mm 2 (Supplementary
The number of oscula and their average area (average OSA, mm 2 ) increased with the sponge volume (V, cm 3 ) in all Mediterranean species. However, the increase of the number of oscula with sponge volume (Table 3A) was more pronounced than the increase in average OSA that showed a marginal increase with sponge volume (Table 3B).

DISCUSSION
Traditionally sponge pumping rates are examined on the individual level (listed in Morganti et al., 2019). However, observations from in situ studies suggested that differences in pumping rate between individuals of the same species were not only related to sponge size, but also to the number and size of the oscula (Gili et al., 1984;Savarese et al., 1997). Recently, theoretical scaling relationships between sponge size, structure, and pumping rates were proposed and confirmed for small sponge specimens in the lab Kealy et al., 2019). To examine if, how, and to what extent these relationships apply to field populations we measured, in situ, the pumping rate of 20 species representative of different morphologies and host types (high-and low-microbial-abundance, HMA and LMA) from temperate and tropical regions. Using these data we examined how oscula number and osculum cross-sectional area varied with sponge volume and pumping rate in order to deepen our understanding of the scaling relationships between sponge size, structure, and functioning, and to test if they can serve as an efficient predictors of sponge pumping rate in the field.
Data acquisition in this study was based on punctual (snapshot) measurements, limited to calm seawater conditions, and only fully functional specimens with entirely opened oscula were sampled. Consequently, temporal effects such as contraction or expansion of oscula (Reiswig, 1971;Kumala et al., 2017;Goldstein et al., 2019) were not accounted for in this study. Although several previous studies documented a daynight cycle in pumping behavior (Reiswig, 1971;McMurray et al., 2014;Strehlow et al., 2016), our day-night pumping rate measurements of Red Sea species showed no clear daily pattern (Supplementary Figure 6). We also did not find such a trend in day-night comparison of the pumping rate of S. vesparium made with acoustic doppler velocimeters (Weisz, 2006).

Mechanisms Underlying the Relationship Between Pumping Rate and Sponge Size
In multi-osculated sponges, each is associated with a cluster of choanocyte chambers and the network of canals that feed and drain them. These functional units were defined as an "aquiferous module" (Fry, 1970(Fry, , 1979Ereskovskii, 2003). The volume of water processed by sponges has been observed to be positively correlated to the density of choanocyte chambers (Massaro et al., 2012). Due to the modular structure of sponges, pumping rate was expected to scale isometrically (i.e., linearly) with sponge volume, under the tacit assumption that the density of the choanocyte chambers is constant (Reiswig, 1975;McMurray et al., 2014;Goldstein et al., 2019;Kealy et al., 2019). Nevertheless, the rate of formation of new aquiferous modules and/or new choanocyte chambers might not increase at the same rate as that of sponge size. In fact, empirical data indicates that allometry rather than isometry dictates the common scaling in most sponges (Reiswig, 1981;Riisgård et al., 1993;Ribes et al., 1999;Kowalke, 2000;Morganti et al., 2019), suggesting that choanocyte chambers density may be reduced with the growth of sponge size, at least in some parts of the sponges body. This hypothesis is supported by observations in Cinachyrella cf. cavernosa where the number of choanocyte chambers increased only marginally with the increase in sponge volume (Dahihande and Thakur, 2019). The multi-osculated Mediterranean sponge species are plotted in the left panels (A,C), the tropical sponges from the Red Sea (reddish color shades) and the Caribbean Sea (blue color shades) are plotted in the right panels (B,D). For the regression analysis, sponge species were clustered accordingly to host type: HMA (closed symbols, solid regression lines) and LMA (empty symbols, dashed regression lines). The scaling exponents and coefficients calculated for HMA and LMA separately are shown in each panel and the shade areas represent the 95% confidence interval for the regression line. The b scaling exponents were statistically different from zero for all equations (p < 0.01) except for the tropical HMA species (n.s, p > 0.05).
Alternatively, structural changes associated with increased size such as an increased length of the canal system may also affect the system drag and hence pumping rate (Leys et al., 2011;Ludeman et al., 2017).
Our data indicate that pumping rate, as well as the osculum area associated with each unit of sponge volume, significantly decreased with sponge volume for many (but not all) of the sponges we studied (Figures 5C,D in this study, and Morganti et al., 2019). In other words, the volume-specific pumping rate (mL min −1 (cm sponge) −3 ) that is often treated as constant (e.g., Weisz et al., 2008;McMurray et al., 2018), is in fact volume dependent in most of the sponges we studied (Figure 5,  Supplementary Figure 7). The lower volume-specific pumping rate observed for larger sponges suggests that the density of the choanocyte chambers is reduced in their larger aquiferous systems. This decrease may be associated with an increase of the density of structural elements (e.g., spicules) that is required to stabilize and strengthen the larger and often higher structure and also to the need for larger canals. Structural differences were observed between different species (Turon et al., 1997) and empirical data showed that differences in the amount of spicules may vary with sponge volume (Barthel and Theede, 1986;Bavestrello et al., 2000;Dahihande and Thakur, 2019). Further studies are needed to explore the variation of structural and functional units over sponge size and age in relation to sponge pumping.

HMA and LMA Comparison
Bacterial community composition and density in sponge tissue vary over several order of magnitudes and the literature commonly differentiates between sponges with high-(HMA) and low-(LMA) microbial-abundance at the two end points of a spectrum of life strategies (Vacelet and Donadey, 1977;Erwin et al., 2011;Giles et al., 2012;Hentschel et al., 2012;Bjork et al., 2013;Poppell et al., 2013). These two guilds are also differentiated by their aquiferous systems  and   Data are expressed as the coefficient a and scaling exponent b ± the 95% CI. n, number of sampled specimens. p indicates the probability that the scaling coefficient b is significantly different from zero.
feeding strategies (reviewed by Maldonado et al., 2012;Morganti et al., 2017). It was previously suggested ) that the volume-specific pumping rate (pumping rate normalized to sponge volume, mL min −1 (cm sponge) −3 ) of HMA sponges is lower in comparison to that of LMA sponges.
In this study, we observed that the volume of HMA species in each habitat was considerably larger than that of the nearby LMA species (Supplementary Figure 5). Larger sponges (and hence HMAs) have a lower total oscula area to volume ratio (Figures 5C,D), suggesting each osculum is draining a larger aquiferous system. Allometric analysis (Figure 5) indicates that the relationships between the volume and total OSA of LMA and HMA sponges is described with similar scaling exponent (b) but very different coefficients (a). For example, the ratio between the a coefficients of LMA and HMA Mediterranean species (Figure 5A, 9.42/3.81 = 2.5) suggests that LMA sponges have 2.5 times the total OSA of similar size HMA sponges (White and Gould, 1965;Gould, 1971). The S similarity coefficient (S = [a LMA /a HMA ] 1/(1−b) , Gould, 1971) is 8.6, suggesting that to attain the same ratio of total OSA per sponge volume, an HMA sponge would need to be roughly nine times larger than LMA sponge (White and Gould, 1965;Gould, 1971).
In the Mediterranean, HMA species possessed on average an aquiferous system that is three times larger than that of the LMAs, and the OFR of HMA sponges was, on average, twice as much as the equivalent size oscula of the nearby LMA sponges. Similarly, the OFR of equivalent size oscula in the tropical sponges was higher in HMAs compared to LMAs species (compare the open and close symbols in Figure 2B), whereas the size of the aquiferous system associated with HMA oscula was >10-folds larger in comparison to that of the tropical LMA sponges (Figures 5B,D). These two strategies might reflect the different aquiferous system organization and development between sponge species and HMA-LMA host type , and differences in the functionality of the choanocyte cells.
Since the choanocyte chamber surface may vary between species (Turon et al., 1997), we cannot exclude the possibility that the different ratio observed between HMA and LMA species might be also due to differences in the function, size, and structure of choanocyte chambers. Interestingly, the morphological (number of oscula and dimensions) and physiological (osculum flow rate) measurements observed in P. ficiformis, which is defined as HMA species based on microbial abundances in its tissue, more closely resemble those from LMA species.
It should be noted that the relationship between total OSA and pumping rate did not differ between HMA and LMA sponges (Figures 3A,B), other than for the tropical HMA sponge S. vesparium (b = 0.78 for tropical HMAs without S. vesparium).

Oscula Number and Their Cross-Sectional Area
In both the temperate and tropical species analyzed in this study, the relationships between the total OSA (i.e., the total osculum cross-sectional area) and sponge volume conformed to a power function allometric scaling with a similar exponent (b = 0.6-0.7), except for tropical HMAs (b = 0.99). A theoretical allometric scaling of the increase in OSA with sponge volume has been recently suggested to be: OSA ∼V 2/3 and this suggestion was confirmed in small specimens in laboratory experiments (Halicondria panicea, b = 0.66, Goldstein et al., 2019). Our results agree with those and previous estimates (Aphrocallistes vastus, b = 0.84, Leys et al., 2011;Cinachyrella cf. cavernosa, b = 0.50, Dahihande and Thakur, 2019) and provide an in situ confirmation of the suggested theoretical exponent based on a very limited size range, on several species that encompass the natural range of both parameters (Table 4). An increase of total OSA with sponge volume can be achieved either by the increase of the number of oscula and/or by the increase of OSA with sponge volume. Due to the high morphological plasticity of sponges, changes in the number of oscula and aquiferous modules may be attributed to different environmental conditions (Plotkin et al., 1999;Ereskovskii, 2003). A recent study observed a depth induced change in osculum morphology in Chondrosia reniformis. This change was attributed to wave action and sediment loading and did not affect total OSA, bacterial clearance, respiration and growth (Gokalp et al., 2020; but see Lesser et al., 2020). As the sampled specimens inhabited the same area under similar environmental conditions, we postulate that the observed variability within each species should be attributed to inherent properties such as sponge size, metabolism, and internal structure, rather than to abiotic factors such as temperature and current (Riisgård et al., 1993).
The network of inhalant and exhalant canals changes considerably with the number of oscula (Gaino et al., 1995) and the number of oscula was observed to increase with sponge volume in this and previous studies (Gili et al., 1984;Dahihande and Thakur, 2019). However, larger sponges also require a more extensive network of canals system that contributes to the scaling b<1 exponent of the increase of pumping rate as a function of sponge volume observed here and in previous studies (Thomassen and Riisgård, 1995;Kowalke, 2000;McMurray et al., 2014;Lewis and Finelli, 2015;Dahihande and Thakur, 2019;Morganti et al., 2019). To what extent the dimensionality of the aquiferous system controls the scaling exponent, i.e., whether the number of choanocyte chambers (and potentially choanocyte number in each chamber) increases isometrically or allometrically with sponge volume requires further studies.

Exploring Pumping Rate at the Osculum Level
The theoretical allometric scaling parameters for osculum jet speed, osculum cross-sectional area, and pumping rate of U 0 ∼ OSA 1/2 and PR ∼ OSA 3/2 rely on the hypothesis suggested by Goldstein et al. (2019) that the pumping units (choanocyte chambers) in sponges are similar size, with similar individual pumping rate and of a similar uniform distribution over sponge volume. Our field data, encompassing a large number of species and a variety of morphologies, do not agree with the prediction of this hypothesis.
Several studies documented a positive correlation between the OSA and U 0 (Strehlow et al., 2016;Goldstein et al., 2019;Kealy et al., 2019), which were in agreement with the theoretical scaling (H. panicea, b=0.39-0.57, Goldstein et al., 2019;Kealy et al., 2019). Our results show an opposite trend of an allometric decrease in U 0 with OSA for two-thirds (12 out of 18) of the species for which we had sufficient available data, including all Mediterranean species (negative scaling exponent, b = −0.39 to −0.76). This discrepancy is not surprising since, as detailed by White and Gould (1965), and Gould (1966aGould ( ,b, 1971, allometric relationships apply only to the range of data from which they are derived. Therefore, such scaling relationship cannot be expected to provide predicting power beyond the size range for which it was established and extrapolation of equations derived from very small oscula as those employed by Goldstein, Kealy, and co-workers (OSA<4 mm 2 ) to sponges with much larger oscula like those in our field studies (OSA up to 39,000 mm 2 ) is unwarranted and in fact may be misleading. The different scaling relationships observed for the tiny explants and very small sponges compared to the fully grown sponges we studied suggest that the underlying mechanisms controlling the variation in pumping rate differ between the two cases. Moreover the discrepancy between our and previous studies can be attributed to a different approach since Goldstein, Kealy and co-workers focused on small explants that dynamically constrict and expand their oscula over periods of minutes to hours, while our measurements are attributed to fully open and static oscula.
It must be noted that one third (6) of the species we studied, exhibited an increase (b>0) in U 0 with OSA, but the number of replicates and OSA range tested was insufficient to draw statistically significant conclusions (Table 4). Overall, the pattern of decrease or increase in U 0 with OSA appeared unrelated to either growth form or host type (HMA-LMA, Table 4).
Similarly, the theoretical scaling of OFR with OSA of PR∼OSA 3/2 suggested by Goldstein et al. (2019) and confirmed for small explants in laboratory experiments (OFR∼OSA 1.45 , Goldstein et al., 2019;OFR∼OSA 1.76 , Kealy et al., 2019) deviated from the allometric exponents (b) we measured for a wide range of sponge size in the field. All but a single b exponent were significantly smaller than one, ranging from 0.36 up to 0.97 (mean b=0.69, Table 4). The tire-shaped Caribbean species S. vesparium was the only exception with essentially isometric scaling of OFR to OSA and an exponent not significantly different from one. The mean relationships we observed at the osculum level (OFR∼OSA 0.69 ) were similar to those observed at the sponge level (PR∼ OSA 0.75 ) for all species other than the tropical HMA species. These scaling relationship closely resembled the value expected for the allometric scaling in the metabolic theory of ecology (b = 0.75) derived from a biological branching network model (West et al., 1997;Brown et al., 2004), which could reflect that of the canals of the aquiferous system assuming its complexity increases with sponge size.
Again, the theoretical considerations that explained the relationships observed for small explants and tiny sponges with dynamic osculum are not expected to hold for much larger and more static sponges. The scaling b<1 exponents we measured, suggest that aquiferous system structure and choanocyte chambers density vary among species, with sponge volume, and potentially with sponge age. The distribution of the choanocyte chambers may also vary in different parts of the sponge body.
As the sponge grows, the aquiferous system may become more complex and each osculum might become associated with a larger sponge volume. Such a trend is indicated by the decrease of the ratio of total OSA to volume with the increase in sponge volume observed for all Mediterranean ( Figure 5C) and may suggest a constraint on the growth of these sponges. In the Mediterranean species, we have surveyed a large number of oscula over a wide range of sponge and oscula sizes. The only tropical sponge for which we have an equivalent sample size is the tire shape sponge S. vesparium that exhibited no variation of the relationship between the total OSA/volume ratio with sponge volume and, therefore, suggests no growth constraints. Different morphotypes (e.g., encrusting, massive, tube, vase-shape) may develop different aquiferous organizations (e.g., multi-osculated sponges and single osculated massive TABLE 4 | Analysis of published data showing the scaling exponents for the allometric function between total osculum cross-sectional area ( OSA), osculum jet speed (U 0 ), osculum flow rate (OFR), sponge volume (V) and pumping rate (PR). Scaling exponent (b) and R squared (R 2 ) are reported.

Predicting Pumping Rates by Osculum Cross-Sectional Area
For the five multi-osculated Mediterranean sponges for which a representative oscula size ranges were measured, we observed that oscula of the same size class pumped at the same rate irrespective of sponge volume. This indicates that the OFR is independent of sponge volume and that OSA should be a good descriptor for pumping rate estimates. Indeed, pumping rate estimates based on the average OFR for each osculum class within a single specimen, were in good agreement (R 2 ≥ 0.76) with pumping rate estimates that used the average OFR for each osculum class of each species (Supplementary Figure 8).
Using oscula number and their average size (average OSA) as predictor variables in a multiple regression analysis explained on average (±95% CI) 87 ± 9% of the variability in the pumping rate (Tables 1 and 2). This represents a significant improvement compared to the use of sponge volume as a predictor variable that explained 67 ± 9% of the variability of the pumping rate.
The variation of the oscula number and their OSA and its relationship to sponge volume differed between the species. The total OSA of a specimen integrates these two variables. Therefore, we explored the relationship between total OSA and sponge pumping rate and compared it with the relationship between sponge volume and pumping rate observed in our previous study on the same sponge species (Morganti et al., 2019). Pumping rates of the whole sponge greatly varied between specimens of the same volume from different species, particularly in large size specimens (see Figure 4 in Morganti et al., 2019). However, when the pumping rate was analyzed as a function of the total area of all its oscula (i.e., total OSA), these differences diminished and the allometric relationships between the pumping rate and total OSA were similar for all examined species (Figure 3).
The high predictive power of total OSA, or the combination of average OSA and oscula number, for sponge pumping rate also applies to sponges from the other examined regions. In the present study, the proportion of pumping rate variation explained by total OSA was higher than that explained by sponge volume for Red Sea and Caribbean Sea species (Table 4). Similarly, OSA appeared to be a better predictor than sponge volume in the hexactinellid sponge Aphrocallistes vastus (Leys et al., 2011) and a mix of demosponges (Southwell et al., 2008), but this is not always the case ( Table 4).
The scaling relationships between OSA and OFR of each species are an efficient and reliable approach to estimate the pumping rate of an entire sponge population or of a sponge community by simply measuring for each species the oscula density and their OSA. The pumping rate of a sponge community can be calculated as: PRpop i L pumped min −1 m −2 = j n j,i (a i OSA b i j,i ) Where, PRpop i is the pumping rate (L pumped min −1 m −2 ) of the population of the ith sponge species, n j,i is the number of oscula of size class j from the ith sponge species per square meter, OSA is the osculum cross-sectional area (mm 2 ) for size class j of the ith sponge species, a i is the allometric coefficient (mL min −1 ) of the ith sponge species and b i is the allometric exponent of the ith sponge species. Similarly, the average OSA and the oscula number per m 2 can be used. A similar approach has been used to estimate the whole reef fluxes of Red Sea sponges (Genin et al., 2009) and deep-sea glass sponge reef in British Columbia (Kahn et al., 2015;Dunham et al., 2018), but note that the allometric relationships between OSA and OFR were not established in those studies.

CONCLUSIONS
This study highlights the importance of the osculum as a pumping operational unit in sponges. Investigating the pumping rates of 20 sponge species from temperate and tropical regions at the osculum level demonstrates that the osculum crosssectional area is an excellent predictor of sponge pumping rate. Our results verified, in situ, the suggested theoretical allometric scaling parameter for OSA∼V 2/3 over a broad range of sponges sizes from different species and geographical regions. Our in situ measurements also indicate that the OFR scaled allometrically with OSA showing a scaling exponent b<1 (average b = 0.69). These data are at odds with previous suggestions for an allometric scaling parameter of OFR∼OSA 3/2 . Pooling together data from all species pointed to a general allometric scaling between total OSA and the pumping rate of PR∼ OSA 3/4 for a variety of morphologies and for both HMA and LMA species in healthy populations under optimal sea conditions. Tropical HMA sponges stood out as an exception and the relationships of pumping rates and OSA in these sponges seem to be isometric (b=1), but more data are required in order to establish this assertion.
Using the number of oscula and average osculum crosssectional area increased the predictive power for PR up to 94% of the explained variability. This represents an improvement of ∼30% over previous estimates based on sponge volume. Once the allometric relationships between OSA and OFR are established for a sponge species, and the effect of seasonality is resolved (see Morganti et al., 2019), simple surveys of the density and size distribution of the oscula provide an efficient and reliable means to estimate the volume of water filtered by a single sponge or a whole population. The error associated with such an estimate is likely smaller than an order of magnitude and thus such data will provide an important contribution to our understanding of the role of sponge in the ecology of benthic systems.
Further study is required to understand the intrinsic mechanisms that drive the relationships between sponge function and structure. In particular, these should elucidate how the aquiferous structure, the choanocyte density and function, and the proportion of the mesohyl material change over sponge size and differ between sponge species, host type, and morphology.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
TM, MR, RM, JW, GY, and RC designed the study, performed the experiments, and analyzed the data. All authors contributed critically to the draft and gave the final approval for publication.

FUNDING
Financial support was provided by the Spanish Government through the Grant (RTI2018-094187-B-100) and the 'Severo Ochoa Centre of Excellence' accreditation (CEX2019-000928-S), and by the Pure Oceans Foundation Grant 2019-Spoplastics to RC and MR; and by a FPU fellowship from Ministerio de Educación, Cultura y Deporte (MECD) and Max Planck Society to TM. Additional financial support was provided by ISF grant 1280/13 and BSF grants 2012089 and 2017622 to GY. This is a contribution from the Marine Biogeochemistry and Global Change research group from the Generalitat de Catalunya (2017SGR1011).

ACKNOWLEDGMENTS
We thank Manel Bolivar and Eduard Serrano for assistance in the fieldwork. We are grateful to the Parc Natural del