Water masses constrain the distribution of deep-sea sponges in the North Atlantic Ocean and Nordic Seas

Water masses are bodies of water with distinctive physical and biogeochemical properties. They impart vertical structure to the deep ocean, participate in circulation, and can be traced over great distances, potentially influencing the distributions of deep-sea fauna. The classic potential temperature−salinity (θ−s) diagram was used to investigate deep-sea sponge (demosponge genus Geodia) association with water masses over the North Atlantic Ocean and Nordic Seas. A novel analysis was conducted, based on sampling the curvature of climatological mean θ−s curves at sponge locations. Sponges were particularly associated with turning points in the θ−s curves, indicative of intermediate and deep water masses. Arctic geodiid species (G. hentscheli and G. parva) were associated with Arctic intermediate and deep waters in the Nordic Seas, and with dense overflows into the northern North Atlantic. Boreal species (G. atlantica, G. barretti, G. macandrewii, and G. phlegraei) were associated with upper and intermediate water masses in the Northeast Atlantic and with upper, Atlantic-derived waters in the Nordic Seas. Taken together with distributional patterns, a link with thermohaline currents was also inferred. We conclude that water masses and major current pathways structure the distribution of a key deep-sea benthic faunal group on an ocean basin scale. This is most likely because of a combination of the physical constraints they place on the dispersal of early life-history stages, ecophysiological adaptation (evolved tolerances) to specific water masses, and the benefits to filter-feeders of certain phenomena linked to water column structure (e.g. nepheloid layers, internal waves/ tides, density-driven currents).


INTRODUCTION
A water mass is a body of water with distinctive ranges of temperature, salinity, and other physical and biogeochemical properties. These properties are acquired at the ocean surface in particular formation regions, where they are determined by the local climate (Pickard & Emery 1990, Tomczak 1999. Following formation, water masses sink and spread along appropriate density surfaces, and can be ABSTRACT: Water masses are bodies of water with distinctive physical and biogeochemical properties. They impart vertical structure to the deep ocean, participate in circulation, and can be traced over great distances, potentially influencing the distributions of deep-sea fauna. The classic potential temperature−salinity (θ−s) diagram was used to investigate deep-sea sponge (demosponge genus Geodia) association with water masses over the North Atlantic Ocean and Nordic Seas. A novel analysis was conducted, based on sampling the curvature of climatological mean θ−s curves at sponge locations. Sponges were particularly associated with turning points in the θ−s curves, indicative of intermediate and deep water masses. Arctic geodiid species (G. hentscheli and G. parva) were associated with Arctic intermediate and deep waters in the Nordic Seas, and with dense overflows into the northern North Atlantic. Boreal species (G. atlantica, G. barretti, G. macandrewii, and G. phlegraei) were associated with upper and intermediate water masses in the Northeast Atlantic and with upper, Atlantic-derived waters in the Nordic Seas. Taken together with distributional patterns, a link with thermohaline currents was also inferred. We conclude that water masses and major current pathways structure the distribution of a key deep-sea benthic faunal group on an ocean basin scale. This is most likely because of a combination of the physical constraints they place on the dispersal of early life-history stages, ecophysiological adaptation (evolved tolerances) to specific water masses, and the benefits to filter-feeders of certain phenomena linked to water column structure (e.g. nepheloid layers, internal waves/ tides, density-driven currents).
KEY WORDS: Water masses · Deep-sea · Sponges · Geodia · North Atlantic Ocean · Nordic Seas OPEN PEN ACCESS CCESS traced over great distances using their characteristic combinations of water properties (e.g. temperatures, salinities, dissolved oxygen concentrations, etc.) (Pinet 2003). Away from formation regions, different water masses interleave and can mix with each other (Tomczak 1999). In many parts of the ocean, the water mass structure is similar over large areas and can be relatively unchanging over time (Emery & Thomson 2004). The concepts of water masses and mixing are also key to understanding deep-ocean circulation (Pickard & Emery 1990, Stewart 2007. In a global ocean that consists of distinct bodies of water occupying finite volumes, mixing with each other, and participating in circulation, it follows that water masses are likely important to the biogeography of deep-sea fauna.
The influence of water masses on the spatial and bathymetric distributions of deep-sea biota has been proposed for various regions and taxonomic groups: pelagic microbial communities in the Nordic Seas and Arctic Ocean (Monier et al. 2013, Le Moine Bauer et al. 2018, Busch et al. 2020; benthic foraminifera in the South Atlantic Ocean off Africa (Schmiedl et al. 1997) and the Pacific Ocean off South America (Ingle et al. 1980); benthic ostracods in the Southwest Pacific (Corrège 1993); benthic amphipods near the Greenland−Iceland−Faroe Ridge of the Northeast Atlantic (Weisshappel & Svavarsson 1998); demersal fish in the North Atlantic (Koslow 1993) and larval fish assemblages off south-western Australia (Muhling et al. 2008); and benthic invertebrates generally, off the east coast of the USA (Rowe & Menzies 1969) and in the Rockall−Porcupine region of the Northeast Atlantic (Gage 1986, Tyler & Zibrowius 1992, Howell et al. 2002, Eerkes-Medrano et al. 2020. The link between water masses and biogeographic distribution is a recurring theme in the scientific literature on cold-water corals (CWCs) (e.g. Dullo et al. 2008, Rüggeberg et al. 2011, Radice et al. 2016, Kenchington et al. 2017. As with sponges (Phylum Porifera), the presently better-studied CWCs also form deep-sea benthic ecosystems (e.g. reefs and mounds) that are often hotspots of biodiversity and biomass, and that provide numerous important ecosystem goods and services (Freiwald & Roberts 2005 and references therein, cf. Hogg et al. 2010 andMaldonado et al. 2017). Notably, Dullo et al. (2008) investigated the relationship between hydrography (in this case, temperature, salinity, density, and dissolved oxygen) and the distribution of the scleractinian coral Lophelia pertusa. Lophelia reefs were found not to be randomly distributed with respect to water column structure, but rather confined to a temperature−salinity field bounded by a distinct water density envelope (Dullo et al. 2008). This finding highlighted the importance of water mass boundaries to the distribution of a deep-sea filter-feeder and ecosystem engineer, analogous to structure-forming sponges in key respects.
The influence of water masses on the distribution of deep-sea sponges has been suggested in a number of studies. Most often, these relate to members of the structure-forming demosponge genus Geodia and to the North Atlantic sponge grounds in which they play a key role, contributing to the provision of various ecosystem goods and services (e.g. habitat provision, benthic−pelagic coupling, biogeochemical cycling, and biotechnological potential) (Hogg et al. 2010, Maldonado et al. 2017). Klitgaard & Tendal (2004) considered the hydrographic conditions at sponge grounds in the Northeast Atlantic and identified several water masses as being potentially important. They also noted that sponge ground distributions appear to follow arcs defined by the main branches of the North Atlantic Drift. In that work, the terms 'warm water ('Atlantic')' (or, alternatively, 'boreal') and 'cold water ('Arctic')' were used to describe 2 types of sponge assemblage, with the link to water masses of different characteristics and origins implied.  acknowledged the idea of species−water mass associations in their taxonomic and biogeographic study of geodiid species in the North Atlantic; they also revealed sites where boreal and Arctic species are sympatric. Beazley et al. (2015) partially attributed the distribution of structure-forming sponges, including Geodia barretti, to water masses present in the Sackville Spur area off Newfoundland, Canada. Knudby et al. (2013) invoked differences in water mass structure to explain the poor predictive performance of certain species distribution models when extrapolated (i.e. models trained on data from one area and used to predict sponge presence in another, with different oceanographic conditions, performed relatively poorly).
In the deep-sea research community, a number of accounts have recently come to light relating sitespecific sponge distributions (e.g. along depth gradients) to local water mass structure (e.g. Roberts et al. 2018, Eerkes-Medrano et al. 2020, andothers as yet unpublished). Furthermore, the potentially critical role of water masses in the broader scale distribution of ecosystem-engineering deep-sea organisms was highlighted by Morato et al. (2020) and (with specific reference to sponges) in a review article by Puerta et al. (2020), which called for improved understanding of the drivers of deep-sea biogeography in the face of projected climate change. As yet, no rigorous oceanographic study demonstrating associations between deep-sea sponges and water masses at a broad (e.g. ocean basin) scale exists. This is a key knowledge gap, particularly as deep-sea sponges may be very slow-growing (Leys & Lauzon 1998), reproduce infrequently (Klitgaard & Tendal 2004), and take millennia to form grounds (Murillo et al. 2016). Reflecting their perceived vulnerability to disturbance and environmental change, deep-sea sponge grounds have been classified as a 'habitat under immediate threat and/or decline' by the OSPAR Commission (OSPAR 2008), and a 'vulnerable marine ecosystem (VME)' by the FAO (FAO 2009).
In the present study, the classic temperature−salinity (T−s) diagram was used to identify water masses throughout the water column at sites hosting deepsea sponges (members of the demosponge genus Geodia) across the North Atlantic Ocean and Nordic Seas. The potential association of sponges with water masses was investigated by plotting sponge occurrences as points in T−s space on a diagram containing T−s curves for sites across the study area, and conducting a novel analysis based on the curvature of the relevant curves sampled at those points. This new approach represents a methodological advance for deep-sea biogeographical studies: so far, observations that species distributions correlate with particular physical, chemical, or biological properties of the water column are confounded by strong collinearity between parameters and also with depth/hydrostatic pressure (see Puerta et al. 2020 and references therein). On this basis, suggested associations with water masses are largely qualitative, whereas here we have developed a means of quantitatively assessing associations with clear, accepted water mass signatures (i.e. based on the curvature of T−s curves).
The purpose of this study was threefold: (1) to test the hypothesis that deep-sea sponges are associated with water masses over broad oceanographic regions (e.g. ocean basins); (2) to test the traditionally held view that boreal and Arctic sponge ground species are associated with different water masses and major current pathways; and (3) to reconcile sites where boreal and Arctic species co-occur with the concept of water-mass-influenced distribution patterns.
The approach adopted here cannot resolve causal mechanisms (i.e. explain why water masses might influence the distribution of deep-sea sponges), although some possible mechanisms are considered in Section 4. It is also not a primary objective of this work to precisely map sponge occurrences onto specific water masses or current pathways, though pertinent observations of such linkages are presented where appropriate. Rather, this paper presents a first broad-scale, quantitative analysis to test the possible association between deep-sea sponges and water masses, and provides a basis upon which future biophysical and biogeographical studies can build.

Overview
First proposed by Helland-Hansen (1916), the T−s diagram is the archetypal 'property versus property' diagram in which subsurface temperatures are plotted against corresponding salinities for a given oceanographic station to produce a curve, the 'T−s curve', useful for interpreting water mass structure (Emery & Thomson 2004). In this work, in situ temperatures were converted to potential temperatures (at surface pressure), θ, to correct for the internal heating caused by the compressive effect of hydrostatic pressure (Emery & Thomson 2004). Therefore, we refer instead to θ−s (rather than T−s) diagrams throughout.
Source water types, sampled close to their sites of formation, are represented as points in θ−s space. Temporal and spatial variability in sampling, or mixtures of 2 adjacent water masses, leads to straight lines (or 'mixing lines') on θ−s diagrams, which are sometimes identified as water masses in their own right. When 3 or more water masses are present, the resulting θ−s curve has turning points, where curvature is relatively high, indicating the influence of intermediate water masses (and deep water masses, if so-called bottom waters are also present) (Fig. 1). The exact degree of curvature is a function of the differences in θ−s properties of the water masses present at a given station and the extent to which the water mass responsible for the turning point influences the water column structure at that station (e.g. the volume/dilution of a source water type present in a certain depth range). Note that a turning point represents the influence of the core of a water mass, whereas boundaries between adjacent water masses fall along the straight, mixing line segments of the θ−s curves, which connect the influence of distinct water masses in θ−s space.
A method was developed to infer mean water temperatures and salinities at sponge localities across the North Atlantic Ocean and Nordic Seas from conductivity-temperature-depth (CTD) data. This allowed sponge records to be plotted onto θ−s diagrams in order to assess possible water mass associations. Clustering of sponges about turning points in the θ−s curves (indicating intermediate and/or deep water masses) was tested for by comparing the frequency distribution of θ−s curvature sampled at sponge localities to that expected if the sponges were randomly distributed with respect to the curves (and thus water mass structure). Systematic differences between these 2 distributions would demonstrate an association between the sponges and water masses. Specifically, sponges should occur more frequently at higher curvatures, in such a case, than can be expected due to chance.
An additional θ−s diagram-based analysis was undertaken to explore near-bed water masses at sites where both boreal and Arctic sponge fauna co-exist. These sites experience high oceanographic variability, and so a near-bed (rather than full water column) focus was applied, with the aim of determining whether multiple water masses influence the bed over time and whether this could explain the occurrence of both assemblages. The key elements of these analyses are addressed in turn below.
The data set consisted of the occurrence records ( Fig. 2) of  http://dx.doi.org/ 10.5061/dryad.td8sb) and Cárdenas & Rapp (2015; data available in the supplementary material of that paper), supplemented with new records from Davis Strait (Bedford Institute of Oceanography unpubl. data), the Jan Mayen and Mohn Ridges (University of Bergen unpubl. data), and Rosemary Bank (Eerkes-Medrano et al. 2020). The records originate from a variety of sampling methods (e.g. trawls, dredges, box cores, remotely operated vehicle (ROV) observations and targeted sampling, etc.) and include at least positional coordinates and water depths. It is likely the most comprehensive and best curated data set of its kind; hence its use (in an earlier form) in the modelling study of Howell et al. (2016). Species identifications have been quality checked and agreed upon by the authors (P.C. and H.T.R.). Identifications were made based on morphology (external and, where possible, spicule morphology) and, in many cases, molecular data (e.g. sequences of the Folmer fragment of the mitochondrial cytochrome c oxidase subunit I [COI] gene). Where there was ambiguity in species identification (e.g. morphological distinction of G. phlegraei and G. parva can be particularly challenging), the record in question was excluded. This data set presents a befitting challenge for the approach adopted here, since it includes species characteristic of boreal (G. atlantica, G. barretti, G. macandrewii, and G. phlegraei) and Arctic sponge grounds (G. hentscheli and G. parva). These are species for which there exists a relatively high number of occurrence records (for deep-sea species) with broad distribution (i.e. 567 records globally, across all 6 species; 551 records in the study domain; Fig. 2). Potential uncertainty in the positions of these records varies from approximately ±10 m (for ROV sampling) to several km (for scientific trawls). Any associated depth uncertainty is assumed to be small relative to the vertical extent of major features in water column structure. Where in situ measurements of environmental parameters, such as temperature and salinity, were available they were mostly single time-point measurements, precluding their use in this study on the basis that temporal variability had not been accounted for. The new data set has been made available at https://doi.org/10.1594/PANGAEA.924531.

Study site delimitation
Study site boundaries were defined according to the grid cells of a 2° × 2° grid which encompassed sponge occurrence records or clusters thereof (Fig. 2). Major bathymetric features (e.g. straits, sea mounts, ridges, fracture zones, gullies, etc.) were used to help discriminate sites. The 2° × 2° grid cell size was selected as a compromise between CTD data availability and the averaging of conditions over unacceptably large areas (e.g. areas larger than oceanographic sub-regions characterised by meso scale horizontal extents on the order of 100s of km). Smaller grid cell sizes yielded few to no CTD profiles in less frequently sampled regions; larger ones entailed integrating conditions over diverse hydrography (and thus failing to characterise the site of interest with sufficient precision). Away from the ocean surface, temperature and salinity behave as conservative properties: there are no appreciable sources or sinks of heat and salt in the interior of the ocean, and thus they can be considered to be changed only by mixing between water masses. Near the surface, this is not the case and the properties are non-conservative, being strongly affected by meteorological conditions and associated surface processes. The θ−s diagram approach adopted here would not typically be applied in waters shallower than ~200 m, and it is also conventional to omit data relating to shallow depths (e.g. < 200 m) from constructed θ−s diagrams (Pickard & Emery 1990). For this reason, several areas with reported sponge occurrences were not defined as study sites and were excluded from further analysis (see Fig. 2). These were shallow continental shelf areas (e.g. the southern part of the Western Barents Sea) and the Mediterranean Sea (where all 4 G. barretti records were shallower than 200 m).
The water column structure of inshore areas that are nevertheless deep-water in nature (e.g. the Norwegian fjords and Norwegian Trench) may not be adequately characterised given the spatial resolution of the current approach (i.e. relevant grid cells incorporate portions of the adjacent shelf). These areas were also excluded from the analyses.

Extracting and preparing CTD data
High-resolution CTD profile (downcast) data were extracted for each site from NOAA's World Ocean Database (WOD18; Boyer et al. 2018). Data were extracted for the period 1 January 1990 to 10 October 2018, which is contemporaneous with the frequent use of instrumentation with high depth resolution. CTD cast data had been interpolated (during the extraction process) at the standard depth levels recommended by the International Association of Physical Oceanography and at additional depth levels as specified by WOD18 (i.e. 137 levels in total, spanning 0−9000 m water depth range, with higher depth resolution closer to the surface) (Boyer et al. 2018).
It should be noted that a small number of sponge records in our data set pre-date the CTD data extrac-tion period (having been collected in the late 1800s/early 1900s). A necessary assumption implicit in our approach is that there has been no major restructuring of key water masses in the North Atlantic and Nordic Seas over this period. An anomalously weak Atlantic Meridional Overturning Circulation (AMOC) has prevailed since approximately AD 1850 (Thornalley et al. 2018) and, superimposed on this state, there is natural decadal and multidecadal variability, and evidence of continued AMOC weakening more recently (Caesar et al. 2018). Whilst our assumption is likely fair for many sites, it may be less so for those undergoing measurable change associated with trends and variability in the AMOC. Nevertheless, we are currently limited to this assumption by a focus in the scientific literature on temporal changes to the properties of specific water masses and to the strength of the AMOC, and an apparent dearth of attention to any associated consequences for water column (vertical) structure in particular locations.
Extracted CTD data were loaded in Ocean Data View (ODV v.5.1.5; Schlitzer 2018). ODV was used to convert measured, in situ temperatures to potential temperatures at sea surface water pressure (reference pressure = 0 dbar), and also for preliminary visualisation and cleaning of the data. Duplicate casts were removed (casts on grid cell boundaries were occasionally duplicated when amal gamating data across several grid cells; duplicates represented 0.8% of the originally extracted casts [average value over all sites]). Casts with erroneous positional coordinates (e.g. placing them on land) were removed (this was extremely rare; 2 casts were removed across the entire data set on this basis). Casts containing obviously spurious temperature and salinity data were removed (this was rare; < 5.0% of the originally extracted casts were identified as spurious and removed [average value over all sites]).
Each site-specific CTD data set was filtered to remove casts conducted in water shallower than 200 m (see Section 2.3). This filtering excluded more casts for sites at the continental shelf edge and slope (where site boundaries incorporated some shallow shelf areas) than for sites that were truly deep water (e.g. seamounts, ridges, and fracture zones). For example, 78.9% of originally extracted CTD casts were filtered out of the Shetland Islands (SHET) data set, on the basis that they were conducted in water shallower than 200 m, whereas 0 casts were removed from the Chaucer and Minia Seamount sites (CHAU and MIN, respectively) on the same basis (NB: aver-age percentage filtered out across all sites was 27.1%). An assumption implicit in this filtering is that the typical water column structure determined just off-shelf applies at adjacent sponge localities on the shelf edge or upper slope. This is a reasonable firstorder approximation for most sites, particularly as those with the most complex hydrography are excluded from the principal (θ−s curvature-based) analysis (as explained in Section 2.5 below).

Depth envelopes of potential temperature and salinity
Envelopes of θ and s versus depth, z, were constructed for all sites shown in Fig. 2, and are presented in Fig. S1 in the Supplement at www.intres.com/articles/suppl/m659p075_supp.pdf. For each site, climatological means (i.e. mean averages of all data available for the extraction period), θ ‾ and s ‾ , were calculated at each depth level present. Standard deviations (SD) of the means were also determined. The depth profiles of θ ‾ and s ‾ were linearly interpolated to 1 m depth resolution and smoothed using a 100 m moving average. Climatological mean potential temperatures and salinities at sponge localities (i.e. θ ‾ Sp and s ‾Sp , respectively) could then be inferred based on their reported depths (Fig. 3, upper panels). The 1 m depth resolution achieved by interpolation simplified this process, since the depths of sponges are often reported to the nearest metre.
The moving average window size applied above (100 m) was selected to remove noise from the climatological θ−s curves constructed later (Section 2.6) whilst preserving persistent features of the water column structure, ahead of curvature analysis (Section 2.7). Near to the ocean surface, this window size is coarse relative to the vertical gradients in θ and s, but these data are not used in further analysis (Section 2.6) and are shown in Fig. S1 and Fig. 3 for completeness only. The number of CTD casts used for each site is shown in Table 1.
The choice to work with climatological (and spatial) mean data necessarily represents some loss of information on temporal (and spatial) variations that could be relevant to sponge distribution. However, away from the ocean surface, the observed variability in the temperature and salinity profiles was very low for most sites (Fig. S1), with little deviation from 'typical' profile shapes (justifying the use of climatological means). The envelopes for Denmark Strait (DEN), Davis Strait (DAV), and the Faroe Islands (FAR) study sites (coloured grey in Fig. 2) indicated higher temperature and salinity variability over much of their depth ranges. Given that these sites require the averaging of data over large areas with diverse hydrography, the climatological mean curves generated were deemed uninstructive and the sites were not in cluded in the primary θ−s curvature analysis which follows.
A separate analysis was designed for and applied to the DEN and DAV sites on the basis that species characteristic of both boreal and Arctic sponge grounds occur together or in close proximity over relatively narrow depth ranges. We outline the rationale and methods for this additional analysis in Section 2.8. The occurrence data presently available for FAR does not indicate extensive areas where boreal and Arctic species may be sympatric, and thus the site was not investigated further.

θ−s diagrams
θ−s scatterplots were constructed for each site (Fig. S2) using the CTD data sets described above (Sections 2.4 and 2.5). Typical, or climatological, θ−s curves were estimated and superimposed (see also the θ−s diagram in Fig. 3) using corresponding θ ‾ and s ‾ value pairs determined from the depth envelopes (Section 2.5). All data relating to the top 200 m of the water column were omitted for reasons explained in Section 2.3.
A 'master' θ−s diagram (see Fig. 4a) was populated with the climatological curves for every site, except DEN, DAV, and FAR (as explained above). This was overlaid with sponge occurrences plotted in θ−s space using θ ‾ Sp and s ‾Sp values inferred as described in Section 2.5 and illustrated in Fig. 3. For consistency with the removal of temperature and salinity data from near the ocean surface, any sponge records shallower than 200 m were also excluded from the diagram and subsequent analyses. They were similarly excluded from site-specific θ−s diagrams presented (e.g. see

θ−s curvature analysis
Curvature values were calculated at frequent, equally spaced (with respect to arc length) intervals along each site-specific, climatological θ−s curve. θ−s diagrams are often presented similarly to that shown in Fig. 4a ( -and y-axes, respectively).
The equation of curvature, κ, at any point (x, y) on a plane curve parameterised by arc length, t, is given by: (1) where x = x(t) and y = y(t) (representing transformed salinity and potential temperature values, respectively) and primes denote first and second derivatives with respect to t.
Parameterising the curves by arc length has the benefit that, since t is monotonically increasing, x(t) and y(t) are single-valued functions, the first and second derivatives of which can be evaluated with relative ease. Arc length, t, from one end of each curve to every x, y data point along that curve, was approximated by cumulative chordal distance (i.e. summation of the lengths of linear segments, determined using Pythagoras' Theorem). Cubic spline interpolation of x versus t and y versus t was then used to determine x and y values at frequent, equal intervals of t (i.e. at 100 positions along the longest curve and proportionately fewer for shorter curves). Using these values, derivatives were evaluated numerically, and curvatures calculated along each curve according to Eq. (1). Standardising the arc length interval for these calculations, in the way described, ensured all curvature values were directly comparable both within and between curves. Normalisation of the data dur-ing its transformation to x−y coordinates had the effect of making the variables in Eq. (1) dimensionless, and all curvatures calculated were therefore likewise dimensionless.
The curvature value assigned to each sponge record was that falling closest (along the relevant site's curve) to the sponge in x−y space. Multiple records of the same species occurring at similar depths (± 200 m) at any given site were removed (leaving a single record) to mitigate the possibility of spatial sampling bias inflating the frequency of certain curvature values. Unusually high curvature values had been assigned to a small number of sponges, reflecting the occasional presence of some high-curvature noise in the curves despite smoothing (Section 2.5). These values were spurious, in that they did not represent persistent water column features, and were removed from the analysis. The final histogram char- acterising curvature at sponge localities was constructed from 151 remaining curvature values/ sponge records (n = 151). The frequency distribution of curvature expected if sponges were distributed randomly with respect to the curves -equivalent to the 'background' curvature distribution -was determined by repeated random sampling; 151 sample curvatures (corresponding to the number of values in the 'sponge localities' case) were drawn at random (and with replacement) from the entire set of calculated curvature values across all curves. This procedure was repeated 10 000 times. The histogram characterising the 'background' case was constructed from the mean averages, over all repetitions, of frequencies for each curvature bin (hence, error bars could be determined for this case; see Fig. 5). The few incidences of very high curvature relating to noise, as identified above, were also removed from the full curvature data set before this repeated sampling was conducted.

Special-case sites: Denmark Strait and Davis Strait
DEN and DAV are special-case sites in that species characteristic of both boreal and Arctic sponge grounds occur together or in close proximity over relatively narrow depth ranges. In a regime whereby sponges are associated with water masses, these groups of species would more usually be separated geographically and/or by depth, and therefore this situation must be investigated with respect to the main hypothesis of the study. Both sites have complex hydrography, hosting water exchanges between major ocean basins and being places where ocean currents interact. Depth envelopes of θ and s (Section 2.5; Fig. S1) confirmed that these sites exhibit high oceanographic variability over much of their water columns. As a result, DEN and DAV were necessarily excluded from the main analyses described above (Sections 2.2−2.7) and (in order to reconcile their boreo-Arctic species mix with the dichotomous situation elsewhere) they were accorded their own distinct physical analysis of near-seabed conditions, as follows.
Near-seabed conditions in specific sub-areas and depth ranges where boreal and Arctic Geodia species co-exist were characterised using subsets of the CTD data described above (see Table 2 for details). Near-bed measurements from CTD casts in these areas were used to populate θ−s diagrams (see Fig. 6). These characterise bottom conditions directly, rather than inferring them based on the adjacent water column structure (as was necessary and appropriate for other sites in this study based on the limited available data at many, and their lower apparent variability at depth). Key water types were indicated on the diagrams to visualise the influence of different water masses at the bed. So-called 'mixing triangles' were outlined to guide interpretation of data points representing mixtures of 2 or more water masses. The corners of each mixing triangle were defined by the 3 major water masses prevalent at that site; any mixture of these water masses will have θ and s values lying within the triangle. The relative proportion of each water mass in the mixture can be estimated (not shown) based on its position within the triangle (see Thomas & Bowers 2012).

Sponge associations with water masses
Some separation of boreal and Arctic Geodia species is apparent in Fig. 4a. G. atlantica, G. barretti, G. macandrewii, and G. phlegraei are largely associated with water less dense than approximately 27.9 kg m −3 (potential density anomaly, σ θ ), whereas G. hentscheli and G. parva are found principally in waters more dense than this. There are few exceptions to this thresholding; notably, a number of G. parva and G. hentscheli records from the Northwest Atlantic (at REYK, SGRN, NLAB, and FLEM) occurring in waters less dense than is typical and a single, deep G. barretti record off the Shetland Islands in water denser than appears typical for this species.
Closer inspection of the data shown in Fig. 4a (e.g. Figs. 4b,c,d & A1) reveals an apparent clustering of sponge records about turning points in their respective θ−s curves. A conservative estimate suggests 18 (out of 28) sites have sponge records occurring at or very near turning points, and 13 of these have clusters of records (i.e. 2 or more) at the turning points.
Apparent associations between sponges and water masses, consistent across several sites (in Fig. A1), reveal likely broad-scale patterns. Water mass definitions are given in Fig. A1 with details in Table S1, and literature sources consulted when interpreting site-specific water mass structures in Fig. A1 are given in Table S2. Strikingly, G. hentscheli and G. parva are associated largely with Arctic Intermediate Water (ArIW) and Arctic Ocean Deep Water (ArODW) in the Nordic Seas (see the JMFZ, JMR, MOHN, SHLZ, SSVAL, and ICEP sites in Fig. A1). In the Northwest Atlantic, these species are also found associated with water masses derived from dense overflows from the Nordic Seas over the Green-  Fig. A1). Full site names are given in Fig. 2  For both Arctic and boreal species, the observed sponge associations with water masses north and south of the Greenland−Scotland Ridge are consistent with the major current pathways of the Atlantic/ Nordic Seas' thermohaline circulation, particularly when considered together with sponge geographical distributions along shelf-edges, slopes, and other major deep bathymetric features hosting boundary currents (see schematics presented in Mauritzen 1996, Curry & Mauritzen 2005, Rudels et al. 2005, Våge et al. 2013. We return to this point in Section 4.1. Fig. 5 shows the results of the θ−s curvature analysis (see Section 2.7), which quantitatively demonstrates the clustering of sponges about high-curvature turning points (i.e. intermediate and deep water masses) on an ocean basin scale. The frequency distribution of the θ−s curvature sampled at sponge localities has less counts at low curvature compared to the power law distribution representing the randomly sampled 'background' case. It also shows a secondary peak at higher curvature values, which is not present in the 'background' case, indicating that more sponges are found at high-curvature turning points than would be expected if they were distributed randomly with respect to the curves (and thus the water mass structure).

Special-case sites: Denmark Strait and Davis Strait
Fig . 6 shows the results of the distinct θ−s analysis designed to investigate the occurrence of both boreal and Arctic sponge ground species in close proximity/together at these 2 high-variability, dynamic sites. Despite the CTD data for the sub-areas of interest (i.e. those hosting boreal and Arctic sponge ground species; Table 2) within these sites being irregularly distributed spatially and temporally, they demonstrate the influence of both Atlantic-and Arcticderived water masses at the seabed over time and relatively small spatial scales.
For the sub-area within DEN (Fig. 6a) (Harden et al. 2016), cannot be excluded. The majority of data points fall more or less centrally within a proposed mixing triangle for these 86 Fig. 5. The frequency distribution of θ−s curvature sampled at sponge localities (n = 151) compared with that from repeated random sampling (with replacement) of the set of climatological mean θ−s curves (n = 151, 10 000 repetitions). Sponges occur less frequently at lower curvature and more frequently at higher curvature segments (i.e. the turning points associated with distinct water masses) than would be expected if they were randomly distributed with respect to the curves. Grey 'background' bars re present mean frequencies determined over 10 000 repetitions of random sampling; error bars are ±1 SD, calculated according to the Poisson distribution of the underlying data for each bin (i.e. ± ) mean freq.
water masses, suggesting that mixtures might be commonplace near the bed at this dynamic and hydrographically complex site. The relative importance of Re-circulating Atlantic Water (RAW) (sensu Mauritzen 1996) cannot be discerned, but it may contribute to the influence of Atlantic-derived water masses in the region.
The equivalent data set for the sub-area of interest in DAV is more sparse (Fig. 6b). The majority of data points indicate the primary influence on the bed is Irminger Atlantic Water, though in this case supplied instead by the West Greenland Current (i.e. [WG]IrAW) (Tang et al. 2004, Azetsu-Scott et al. 2012. Mixtures of (WG)IrAW with other water masses are also indicated: namely, mixtures with colder, fresher Arctic Water, derived from Arctic Ocean outflow through the Canadian Arctic Archipelago (i.e. [C]ArW), and with Baffin Bay Deep Water (BBDW), the origin of which is yet to be settled unequivocally (Tang et al. 2004).   Table 2. Details of the analysis to characterise water masses influencing the bed at sites where boreal and Arctic Geodia species co-occur. Sub-areas of interest are specified using the opposite corners (southwestern-most to northeastern-most) of rectangular bounding boxes. The number of CTD casts includes only those that terminated in the depth range of interest, from which bottom values of temperature and salinity were extracted for the analysis 4. DISCUSSION

Deep-sea geodiid sponges constrained by water masses and currents
In this study, a novel analysis was applied to a data set incorporating hydrographic data for 28 sites across the North Atlantic Ocean and Nordic Seas (capturing many different water masses) and occurrence records for 6 species of deep-sea geodiid sponge. The records for some species spanned temperature and salinity envelopes almost as broad as those characterising the whole study region. Over the entire data set, there was a statistical relationship between deep-sea sponges and high-curvature turning points in θ−s curves, which indicate the influence of intermediate and deep water masses on water column structure. Sponges are not uniformly distributed over their temperature and salinity ranges, but are instead clustered at points where there are clear water mass intrusions affecting the water column. The results suggest that, within the broader temperature and salinity envelopes/tolerances of these species, there is a role for water masses in structuring sponge biogeography. Whilst it remains difficult to disentangle which environmental factors are most important to the sponges, our findings highlight the need to consider deep-sea sponge distributions as being constrained by water masses (and current pathways), albeit potentially through a plethora of associated mechanisms (of which ecophysiological adaptation to factors such as temperature, which has received much attention historically, is but one -see Section 4.2).
So-called Arctic geodiid species (Geodia hentscheli and G. parva) are indeed associated with Arctic intermediate and deep waters in the Nordic Seas, and with dense water overflows into the North Atlantic Ocean over the Greenland−Scotland Ridge. In contrast, boreal species (G. atlantica, G. barretti, G. macandrewii, and G. phlegraei) are generally associated with upper and intermediate water masses in the Northeast Atlantic and with upper, Atlanticderived waters in the Nordic Seas. In the Northwest and Central North Atlantic Ocean, they are also found in possible association with so-called deep water masses.
The observed sponge associations with water masses north and south of the Greenland−Scotland Ridge (taken together with sponge geographical distributions along shelf-edges, slopes and other major deep bathymetric features) reflect the major current pathways of the Atlantic/Nordic Seas' thermohaline circulation (Mauritzen 1996, Curry & Mauritzen 2005, Rudels et al. 2005, Våge et al. 2013. Boreal geodiid types in the Northeast Atlantic and Nordic Seas appear to correspond to the various branches of the North Atlantic and Norwegian Atlantic Currents. Arctic types in the Nordic Seas correspond approximately to the East Greenland, Jan Mayen, and East Icelandic Currents, and in the Northwest Atlantic both boreal and Arctic species appear associated with deep slope/boundary currents (in particular, the Deep Western Boundary Current off the eastern Canadian margin). Questions remain, however, over the existence and nature (e.g. direction) of mechanisms connecting boreal geodiid species at depth in the eastern, central, and western regions of the North Atlantic. Rigorous genetic and physical connectivity studies will be required to resolve such questions.
Where Arctic and boreal species co-occur, the regions are invariably ones of complex hydrography, hosting water masses of both Atlantic and Arctic origin and being characterised by turbulent mixing/ entrainment (e.g. DEN and DAV). On this basis we predict that, other conditions being favourable (e.g. substrate, sedimentation, etc.), a number of other areas may support currently unknown or additional aggregations of sponges incorporating both Arctic and boreal geodiid species: the Faroe−Shetland Channel/Faroe Bank Channel/Wyville Thomson Ridge area; the North and Northeast Iceland Shelf; the Fram Strait (particularly where Re-circulating Atlantic Water or the West Spitsbergen Current interact with Arctic Ocean outflow); the Barents Sea (in the vicinity of the Polar Front and between Svalbard, Franz Josef Land, and Novaya Zemlya); and east Baffin Bay (along the pathway of Atlanticderived water penetrating into the bay).

Causal mechanisms and directions for future work
Why then are sponges associated with specific water masses and currents? North Atlantic Geodia species are known to be gonochoristic and oviparous (Spetland et al. 2007, Koutsouveli et al. 2020. Unfortunately, Geodia larvae have never been observed and their characteristics are unknown -the larvae in the whole Astrophorina suborder are unknownmaking it difficult to estimate dispersal capabilities. G. parva and G. hentscheli are also known to reproduce asexually by budding , with buds (fragments of the parent sponge) of diam-eter ~1 cm or less being observed. The authors have observed Geodia juveniles, with diameters on the order of 1 mm, free-living in sediment samples from within sponge grounds in the Nordic Seas and Arctic Ocean (authors' unpubl. obs.), although it is unknown whether these juveniles result predominantly from local sexual or asexual reproduction. It seems likely that the larvae (or resuspended juveniles) are constrained within water masses and/or major ocean currents by regions of density stratification or current shear. The implication is that the main agents of Geodia dispersal are close to neutrally buoyant and passive, or have weak swimming capability (in the case of larvae), and thus their distribution is dominated by physical oceanographic processes and settings. Alternatively, larvae may possess the ability to move away from such gradients (in order to remain within optimal conditions), and are thus likewise bound to water masses and/or currents. McManus & Woodson (2012) and Gary et al. (2020) highlighted the potential of larval characteristics (behavioural and buoyancy) to determine patterns of distribution across large areas of ocean. Based on the results presented here, we concur, and note that studies of deep-sea sponge larvae are required if we are to further refine our understanding of the link between water masses, currents, and deep-sea sponge biogeography.
The physiological tolerances of early and adult lifehistory stages with respect to temperature, salinity, pressure, pH, oxygen, nutrient, and food availabilities have also been implicated in structuring the geographic and bathymetric ranges of certain deep-sea taxa (Young et al. 1996, Arellano & Young 2011, Brooke et al. 2013. Deep-sea sponges may become adapted to environmental conditions prevailing within given water masses, their products and derivatives, and given current systems, experiencing these conditions as optimal, with larval stages realising enhanced recruitment success and/or adult stages thriving. Conditions in adjacent water masses may be experienced as sub-optimal, with potentially lethal or sub-lethal adverse consequences (e.g. Tyler et al. 2000, Strand et al. 2017. The long-term persistence of deep-sea sponge grounds in the face of environmental change will likely depend on how that change compares to the levels of 'typical' or natural variability the sponges may have adapted to tolerate (e.g. the enhanced temperature variability tolerated by a dense aggregation of sponges in the Faroe− Shetland Channel, as reported by Davison et al. 2019). We return to this point in Section 4.3 below.
A further consideration is that aggregations of benthic filter-feeders benefit from a supply of food and recruits, and (generally) from enhanced currents (Genin et al. 1986, Rice et al. 1990, White 2003. Nepheloid layers, turbid layers containing elevated concentrations of suspended matter, provide a feasible source of food and larvae/very small juveniles. They are primarily the result of resuspension of bottom sediments by strong currents, for example those associated with 'critically incident' internal waves and tides at the continental slope/shelf-edge (Dickson & McCave 1986, McCave 1986, Thorpe & White 1988. Intermediate nepheloid layers (INLs) (i.e. those detached from the seabed and extending from the slope into the adjacent water column), in particular, are known to act as vehicles for episodic larval transport (Ryan et al. 2010). Such layers tend to spread laterally along isopycnal surfaces and to be carried by large-scale current systems (sometimes achieving considerable spatial extents). Thus, their distribution patterns are intimately linked with the density structure of the deep ocean and with its circulation (Dickson & McCave 1986, McCave 1986, Thorpe & White 1988). Our study cannot resolve whether sponges are associated with the core of a water mass or with a water mass boundary supporting internal waves and nepheloid layers. However, it seems likely that mechanisms (e.g. internal waves or density-driven currents, plus nepheloid layers) that can transport larvae (or very small juveniles), supply food, and enhance currents locally will structure the distribution of certain benthic invertebrates (including sponges; see Davison et al. 2019) in a way that reflects water column structure and other oceanographic processes.
In reality, the observed association between deepsea sponges and water masses is likely the result of a combination of the factors mentioned above. Directions for future research should aim to populate analyses, such as those presented here, with more sponge records, and to establish whether sponge presence/absence observations, known sponge ground locations, or sponge density data correlate with the following: levels of stratification in the water column (buoyancy frequency measurements); persistent or recurring INLs (optical turbidity measurements); patterns of internal wave/tidal energy dissipation (turbulent kinetic energy dissipation measurements and models); and major thermohaline current pathways (velocity measurements and hydrodynamic models). Targeted sampling of nepheloid layers is required to assess food quality and to check for early life-history stages amongst any associated zooplankton, perhaps using modern, high-technology sampling methods (see Ryan et al. 2010). Efforts to investigate the physi-ological tolerances of deep-sea sponges should continue (e.g. after Maldonado & Young 1998), as should those to improve understanding of interactions between sponges and flow regimes over various spatial and temporal scales (e.g. after the numerical modelling work of Culwick et al. 2020). In particular, we echo the sentiments of Maldonado (2006) in stating that our understanding of sponge larval ecology is still lacking generally. Whilst it may be possible to draw inferences based on shallow-water studies, dedicated research on deep-sea sponge larval characteristics is preferable and should be prioritised.

Biogeographic implications under climate change
Observations and model predictions of the variability of water masses and thermohaline circulation reflect both natural, multidecadal oscillations (Bryden et al. 2003, Latif et al. 2006) and trends potentially attributable to anthropogenic climate change (Böning et al. 2008). Subsurface temperature and salinity are likely sensitive to anthropogenic change (Banks & Wood 2002), though Banks et al. (2000) concluded that, in Northern Hemisphere oceans where climate variability is inherently large, the signals of anthropogenic climate change were not readily detectable above natural variations. With a further 2 decades passing since the publication of that work, it seems likely that the conclusion is no longer valid (see Thornalley et al. 2018, Bilbao et al. 2019, Spooner et al. 2020. Although model predictions diverge on the response of the thermohaline circulation to global warming (Schweckendiek & Willebrand 2005), we note that the distributions of our target sponge species overlap with the major sites of water mass formation in the Nordic Seas and subpolar North Atlantic. Convection and water mass formation/transformation in these areas is undergoing change (Dickson et al. 1996, Schweckendiek & Willebrand 2005, Brakstad et al. 2019) and could be responding to anthropogenic warming (Yang et al. 2016 and references therein). In agreement with Kazanidis et al. (2019) and Puerta et al. (2020), our results suggest that changes to water mass properties, distributions, and transport under climate change will likely affect present-day sponge populations and influence their future distributions. The capacity of deep-sea sponges to respond to such environmental change is largely unknown. However, recent work from the Northwest Atlantic has shown that geodiid sponge grounds on the slopes of the Flemish Cap and Grand Bank (Murillo et al. 2016) and aggregations of the glass sponge Vazella pourtalesi on the Scotian Shelf, off Canada (Beazley et al. 2018), have persisted in the face of considerable environmental variability (since at least 130 ka BP, in the case of the geodiid grounds). Dullo et al. (2008) proposed using an envelope of potential density, associated with the occurrence of CWC reefs in the Northeast Atlantic and Norwegian Sea, as a tool to help locate other reefs at less explored sites. Our findings suggest that turning points in θ−s curves (and the presence of major current pathways) could be used in a similar way for aggregations of deep-sea sponges. Furthermore, an Arctic or boreal sponge assemblage may be predicted based on the nature and origin of the intermediate or deep water mass responsible for a given turning point.

Ecological significance and utility in palaeoceanography and conservation
This concept could be extended to use in species distribution (or habitat suitability) models (SDMs), whereby basin-scale predictive mapping of potential niches could be achieved using θ−s curvatures calculated from widely available temperature and salinity profile data. Ideally, however, this would be preceded by work to establish the environmental factors and mechanisms associated with water masses that are of direct ecological relevance to the sponges (as discussed in Section 4.2), and the choice of predictor variables (environmental layers) for modelling would be informed by this work. For example, if internal waves/tides were found to be important generally, buoyancy frequency could be included as a predictor in SDMs, or a 'bed slope criticalness factor' (a measure of likely bottom current intensification) could be used, relating the bed slope to internal wave ray slopes (after Frederiksen et al. 1992). If nepheloid layers were determined to be key, environmental layers characterising near-bed turbidity or suspended sediment concentration could be employed as predictor variables, and so on for the other potentially important factors identified in Section 4.2.
We also propose an inversion of the above concept, whereby the presence of particular deep-sea sponges has utility in providing information on water masses or boundary current systems influencing the seabed. By extension, the sponge spicule record from seabed cores could be examined to reconstruct palaeoenvironmental conditions within, and spatial adjustments of, specific water masses and current systems. A ten-dency of deep-sea sponges to have distributions which track particular water masses or currents increases their value as palaeoceanographic archives (e.g. Ellwood et al. 2006, Hendry et al. 2010, Jochum et al. 2012, and others) because they allow the study of specific components of the global ocean's density structure and thermohaline circulation over time. The same tendency, and its potential for validation into the past, makes deep-sea sponges good candidates for the prediction of future range shifts under climate change scenarios.
Finally, deep-sea sponge association with water masses, their products and derivatives over large geographical ranges and relatively narrow depth ranges implies that barriers to dispersal in the deep sea may be as significant (if not more so) in the vertical (e.g. because of density or current stratification) as they are in the horizontal (e.g. because of topographic obstruction and prevailing current direction). This is in agreement with the findings of several other authors (Zardus et al. 2006, Miller et al. 2011, Decker et al. 2012, Jennings et al. 2013, van Soest & de Voogd 2015, employing different methodologies. Areas where several different water masses are brought into close proximity by topographic constrictions and oceanographic processes (e.g. at sea straits, ridges, and sills), and which support turbulent mixing and entrainment, likely have special ecological significance in terms of diversity and genetic exchange, and may warrant international conservation effort.

CONCLUSIONS
In the North Atlantic Ocean and Nordic Seas, deep-sea sponges are particularly associated with turning points in θ−s curves indicating the influence of intermediate and deep water masses. The observed water mass associations (taken together with sponge distribution patterns) indicate a probable link with the major current pathways of the Atlantic Ocean's thermohaline circulation. Where Arctic and boreal species co-occur, the regions are invariably ones of complex hydrography, hosting water masses of both Atlantic and Arctic origin (e.g. DEN and DAV). The causal mechanism is proposed to be a combination of the following: (1) sponge larvae (or small juveniles) are constrained within water masses and/or major ocean currents by density stratification or current shear; (2) sponges have become adapted to tolerate/thrive in the prevailing environmental conditions within given water masses and current systems; and (3) aggregations of benthic filter-feeders benefit from a supply of food and recruits, and from enhanced currents, and thus persist at depths hosting internal waves/tides, density-driven currents, and nepheloid layers (phenomena intimately linked to the density [water mass] structure of the deep ocean).
Our results imply that present-day sponge populations and their future distributions will be affected by changes in water mass properties, distributions, and transport under climate change. We note that turning points in θ−s curves could be used to identify currently unknown sites of deep-sea sponge aggregations and, conversely, that deep-sea sponges and spicule deposits have utility in studying water masses or boundary current systems influencing the seabed now and in the past. The broadest implication of this work is that barriers to dispersal in the deep sea may be as significant (if not more so) in the vertical (e.g. because of density or current stratification) as they are in the horizontal (e.g. because of topographic obstruction and prevailing current direction). Highly dynamic areas, where several water masses are brought into close proximity (e.g. at sea straits, ridges, and sills) and the constraints on dispersal are somewhat mitigated by mixing, likely represent important hotspots of genetic exchange and diversity in the deep sea. Science, Aberdeen, UK) are thanked for providing additional Geodia records and material; samples from the Northwest Atlantic were collected and shipped to E.K. by Tim Siferd and Wojciech Walkusz, Fisheries and Oceans Canada (Winnipeg, Manitoba), as part of the research project 'Identification and mapping through predictive modelling of coldwater coral and sponge species in the Sub-Arctic/Eastern Arctic ' (2014−2017). That project was funded through Fisheries and Oceans Canada's International Governance Strategy (IGS) Research Fund to E.K. Samples were identified by Gabrielle Tompkins-MacDonald, NSERC Visiting Fellow in the Kenchington lab, and technician Emily Baker, also of the Kenchington lab. Jack Middelburg and 3 anonymous reviewers are thanked for their enthusiastic support and comments, which led to improvements to the manuscript. E.M.R. is grateful to the Johnson family and to Malen Roberts for personal support whilst writing the manuscript during the 2020 Coronavirus pandemic. This research has been performed in the scope of the SponGES project, which received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 679849. 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. Fig. A1. Potential temperature-salinity (θ−s) diagrams showing climatological mean curves for sites across the North Atlantic Ocean and Nordic Seas, with sponge occurrences overlaid. Site name abbreviations are defined in Fig. 2 and Table 1. Key water masses have been labelled where possible, and their full names can be found in the key, which defines water mass abbreviations used throughout this article. Water mass descriptions are provided in Table S1. Water mass names could not be assigned to the upper water column in ICEP and NEIS because of variability at these depths. Note the different θ and s axes scales applied from panel to panel