Oceanographic variability drives the distribution but not the density of the aggregation forming deep-sea sponge

,


Introduction
Sponges (Porifera) are important benthic habitat-forming organisms, known to provide critical ecosystem services (Ottaviani, 2020).Aggregating deep-sea sponges increase the structural complexity of the environment resulting in a greater variety of fauna and enhanced local biodiversity compared to adjacent non-sponge habitats (Beazley et al., 2013, Maldonado et al., 2015).Sponges host sessile epifauna (Klitgaard, 1995) including echinoderms, cnidarians, molluscs, bryozoans and crustaceans that settle on their surface (Kazanidis et al., 2016), as well as residing within the channels and pits of the sponge's tissues (Wulff 2006;Buhl-Mortensen et al., 2010).
After hexactinellid sponges have died, spicule mats continue to provide an important substratum, with their own diverse faunal communities (Bett and Rice 1992;Gutt et al., 2013).Aggregations may also act as refugia, nurseries, and foraging sites for demersal fish species (Kenchington et al., 2013;Clark et al., 2016), including those of commercial importance such as Sebastes sp.(Freese and Wing 2003).Additionally, sponges host a diverse microbial community (Thomas et al., 2016).These microbes yield an equally diverse variety of bioactive metabolites which may have important biomedical applications in the future, with many compounds displaying antiviral, antimicrobial, and antiprotozoal qualities (Sharma et al., 2016;Koch et al., 2021).Sponge aggregations also play significant roles in biogeochemical cycling (carbon, nitrogen, and silicon) and benthic-pelagic coupling in particular sponge grounds found adjacent to coral reefs (Reiswig 1975;Yahel et al., 2003;Hadas et al., 2006;de Goeij et al., 2008;Hoffmann et al., 2009;Cathalot et al., 2015;De Clippele et al., 2021).The high carbon sequestration rates of some glass sponge reefs, achieved via their suspension-feeding and microbe associates, are comparable to that of "old-growth forests" and other marine carbon sequesters such as Macrocyctis kelp (Dunham et al., 2018).
Sponges are slow-growing and long-lived (Leys and Lauzon, 1998;Fallon et al., 2010) so they are sensitive to disturbance (Dinwoodie et al., 2021).Consequently, sponges are particularly vulnerable to bottom fisheries activity (e.g.beam trawls) which causes mechanical damage, dislodgement (Wassenberg et al., 2002), and sedimentation (Tjensvoll et al., 2013).The effects of trawling activity on deep-sea sponge aggregations have been observed on the Porcupine Seabight (Vieira et al., 2020), Flemish Cap (Pham et al., 2019), Schulz Bank (Morrison et al., 2020) and the Gulf of Alaska (Freese 2001).Vieira et al. (2020) re-surveyed Pheronema carpenteri aggregations first reported by Rice et al. (1990) in the Porcupine Seabight across four key transects in 1983/4.Since the original surveys, peak density has decreased from 1.09 to 0.03 indv/m 2 , whilst biomass density has decreased from 246 to 4 gwwt/m 2 .During this time, commercial bottom trawling occurred within the location of P. carpenteri aggregations, and multiple trawl scars were observed.Vieira et al. (2020) conclude that their observations suggest a substantial and long-lasting negative impact on P. carpenteri by bottom trawl fishing.
Sponges are also vulnerable to other anthropogenic activities including oil and gas exploration (Vad et al., 2018), and deep-sea mining (Simon-Lledó et al., 2019).The sensitivity of deep-sea sponge aggregations to human activity, as well as their functional significance and structural complexity, has resulted in their classification as vulnerable marine ecosystems (VMEs) by the United Nations General Assembly, Resolution 61/105 (UNGA, 2016;FAO 2009).Additionally, deep-sea sponge aggregations are also listed as an OSPAR (Convention for the Protection of the Marine Environment of the North-East Atlantic, OSPAR, 1998) threatened and/or declining habitat (https://www.ospar.org/documents?v=7234).The 'listed' status assigned to sponge aggregations has led to the implementation of a series of conservation and management measures within individual states' exclusive economic zones (EEZs), and in areas beyond national jurisdiction (ABNJ).Such measures are largely based upon restricting fishing activities, such as area-based closures to bottom fisheries, and the implementation of bycatch reporting frameworks by regional fisheries management organisations (RFMOs).Other fisheries-based measures included the "Move-on Rule" which is in place in RFMO regions.The Move-On Rule requires fishing vessels to cease activity and move a minimum distance away from a location where there is evidence a VME has been encountered, e.g. a specified threshold of a VME indicator has been caught in fishing gear, details of which are determined by individual RFMOs.In addition, area-based management tools are also used, such as marine protected areas (MPAs).Across the northeast Atlantic, the main policy drivers of deep-sea MPAs are the European Union's Habitats Directive (92/43/EEC; European Commission, 1992) and OSPAR Convention (OSPAR, 1998).
Central to area-based management is the need to understand the spatial distribution of a habitat or species.Critical to this is identifying and understanding the environmental drivers of a habitat or species' distribution as they characterise the fundamental niche (Hutchinson 1957;Soberón et al., 2017).This underlying need requires observational and environmental data, collection of which poses difficult challengesparticularly with the formerin the context of the deep sea.The deep sea has been historically under-sampled, so fundamental data describing species or habitat distribution is often lacking.Of that available data, there is a strong bias in the spatial distribution of data, with much deep-sea research being historically focussed in the North Atlantic and the EEZs of developed nations (Hughes et al., 2021).Habitat suitability models are tools that use statistical or machine learning means to create complete coverage maps of predicted distribution (presence-absence, density, or abundance-based) from relatively sparse observational data; these maps can then be used to aid marine spatial planning and conservation (Ross et al., 2015;Howell et al., 2016;Lauria et al., 2017;Ramiro-Sánchez et al., 2019;Burgos et al., 2020).However, critical to making accurate predictions is an understanding of what drives a species' or habitat's distribution or density, allowing modellers to select the most appropriate predictor variables.Despite the advances made in our understanding of sponge ecology and the ecosystem services they provide, drivers of deep-sea sponge distribution, in particular their density or abundance, are poorly understood.Important environmental drivers of deep-sea sponge distribution, identified so far from habitat suitability modelling studies, include temperature, depth, silicate concentration, sediment grain size, and minimum bottom salinity (Ross and Howell 2013;Knudby et al., 2013;Miller et al., 2015;Ross et al., 2015;Howell et al., 2016;Rooper et al., 2017;Ramiro-Sánchez et al., 2019;Beazley et al., 2021).
Studies have also observed relationships between sponge distribution and density, and hydrodynamic features and processes.Local hydrodynamics play an important role in supplying food particles (e.g.particulate organic matter) to sponges (Hanz et al., 2021), whilst also preventing the accumulation and subsequent smothering of sponges by sediment deposition (Whitney et al., 2005;Roberts et al., 2018).Rice et al. (1990) first proposed an association between the P. carpenteri density in the Porcupine Seabight and regions "downstream" of high near-bed currents, enhanced by internal tides (internal waves of tidal frequency).P. carpenteri forms the most extensive Hexactinellid aggregations in the North-East Atlantic.Communities composed of P. carpenteri have been described by Rice et al. (1990) from 1250 m on the Porcupine Seabight where peak densities of 1.5 indv/m 2 and biomass up to 10 g m − 2 ash-free dry weight reported; by Barthel et al. (1996) off of Morocco at depths of 740-1300 m (peak density of 6 indv/m 2 ); by Le Danois (1948) from Ireland to Spain in 1000-2000 m water depth; by Hughes and Gage (2004) from the Rockall-Hatton Basin at 1100 m (mean density of 1.53 indv/m 2 ); and from 1450 m on Goban Spur by Duineveld et al. (1997), Flach et al. (1998) and Lavaleye et al. (2002) (reviewed in Howell, 2010).There are indications that this species may also be common to the west of the Faroe Islands and south of Iceland at depths of between 800 and 1160m Copley et al. (1996).Rice et al. (1990) suggest that P. carpenteri are not able to withstand high exposure to strong currents, but from a peripheral location, could still utilise the increased organic particulates (i.e.marine snow) delivered via internal waves.Although few studies have been conducted to fully explore these relationships, White (2003) supported this hypothesis, finding that P. carpenteri were found in regions adjacent to enhanced near-bed tidal currents in the Porcupine Seabight.However, it would appear this hypothesis is not uniform across taxa nor locations as numerous sponge aggregations are known to occur within regions of tidally modulated near-bed currents (Klitgaard et al., 1997;Whitney et al., 2005;Roberts et al., 2018).For example, Ostur (demosponge) aggregations are found to withstand strong bottom currents caused by internal waves in the Faroe Shetland Channel (Davison et al., 2019).The current hypothesis is that this is in part linked to differences in external morphology between taxa.P. carpenteri are delicate and weakly anchored into the sediment by their spicule mats.Comparatively, ostur species (e.g.Geodia spp.) are characteristically attached to hard rock and gravel/cobble sedimentseven regularly incorporating gravel into their bodies (Klitgaard and Tendal 2001) which creates more secure anchorage, allowing better resistance to high current speeds (Davison et al., 2019).
Internal waves are detectable through fluctuations (variability) in temperature (Al Senafi and Anis 2020).Internal waves, driven by internal tides, propagate at the interface (isopycnal) between two water masses, before breaking onto sloping topography.As internal waves propagate, the stratified water masses (isopycnal) physically rise and fall along slopes by up to hundreds of metres.This motion results in fluctuations in local temperatures at the slope (see Fig. 1c in Reid et al., 2019).Additionally, breaking internal waves induce turbulent bottom currents that re-suspend sediments and particulate organic matter, as well as generate and maintain near-bottom nephaloid layers.As a result of the relationship between temperature fluctuations and internal wave activity, temperature variability has been previously used in benthic ecology studies as a proxy for internal wave activity (Davison et al., 2019), with success, to investigate its influence on driving species distribution and density.
This study investigates the relationship between temperature variability (as a proxy for internal wave activity), along with other environmental variables, and the distribution and density of P. carpenteri in the Porcupine Seabight and Goban Spur.Our wider objective was to also test whether gridded temperature variability successfully captures any relationships between internal wave activity and P. carpenteri at the same spatial scales used in habitat suitability models (e.g.200 × 200 m).If achieved, this would be an important step forward as these processes could be incorporated into future habitat suitability models.This would help improve model performance, and ultimately help map vulnerable marine ecosystems and inform marine spatial planning.

Study area
Central to the study area is the Porcupine Seabight, an amphitheatreshaped impression along the Irish continental slope, over 200 km to the southwest of Ireland (North-East Atlantic) (Fig. 1).The Porcupine Seabight is bound by the Porcupine Bank to the north and Goban Spur to the south, both of which are features of the Irish continental slope and fall within the area of study.The Porcupine Seabight slopes steadily down from 200 m at the Irish Shelf to 4000 m as it widens into the Porcupine Abyssal Plain (PAP) to the southwest.The oceanographic setting is dominated by a poleward flowing current along the continental slope, characterised by Eastern North Atlantic Water (ENAW) and Mediterranean Outflow Water (MOW).The current regime within the Porcupine Seabight is dominated by cyclonic flow with mean residual bottom flows between 2 and 5 cm s− 1 (White 2007).
This study area was chosen for several reasons; firstly, because of the availability of data.Several aggregations of P. carpenteri are known to occur within the Porcupine Seabight and Goban Spur region, for which data were available.High resolution (200 × 200 m), full coverage environmental data were also available for this region, including outputs from a regional oceanographic modelling system (Nagy et al., 2020) and bathymetric data.The complete spatial coverage of environmental data, including oceanographic data, means for this region it was possible to analyse all P. carpenteri aggregations where distribution data were available.This also builds upon the spatial limitations of previous work, e.g.White (2003), where only two in-situ moorings were able to be used across the study area.
Secondly, there is a conservation and management interest in the region.Historical trawling has occurred in the Porcupine Seabight and, in turn, had negative effects on P. carpenteri aggregations (Vieira et al., 2020).Other anthropogenic activity in the region also includes 12 oil and gas licences for 'exploration' (EMODnet) which span the northern and eastern flanks of the Porcupine Seabight.The combination of anthropogenic activity and VMEs in the region means that mapping VMEs and identifying underlying drivers of distribution is critical.There are currently no marine protected areas or any other area-based management designated within the study area for the protection of deep-sea sponge aggregations.However, the North-East Atlantic Fisheries Commission's (NEAFC) "Move-On Rule" does apply to this region; stipulating that if more than 400 kg of live sponge is present in a single trawl tow, an encounter with a VME has occurred.In this instance, vessels must cease fishing immediately and move at least two nautical miles away.Additionally, the European Union's Common Fisheries Policy regulations (2016/2336) state that "no fishing authorisation shall be issued for the purpose of fishing with bottom trawls at a depth below 800 m" throughout the EU EEZ, which applies to the study area, thereby Fig. 1.The red-bounding box in the map insert highlights the study region, 100 km to the southwest of Ireland.Central to the study area are the Porcupine Seabight (PSB) and Goban Spur (GS), as well as Porcupine Bank (PB) and the Porcupine Abyssal Plain (PAP).The main map illustrates the average annual temperature variability across the study area, overlaid with the presence-pseudo-absence distribution of P. carpenteri.Contours provided by GEBCO.A thick red contour denotes the 800 m contour.Map projection: Goode Homolosine Ocean.
protecting any sponge aggregations below 800 m from direct bottom trawling activity.

Environmental data
Publicly available bathymetric data were obtained from INFOMAR (The Integrated Mapping for the Sustainable Development of Ireland's Marine Resource Programme) at a 200 × 200 m resolution.Using the Benthic Terrain Modeller 3.0 (Walbridge et al., 2018) add-on application in ArcMap (Version 10.7), the following bathymetry-derived variables were calculated at 200 × 200 m resolution: broad-scale bathymetric position index (BBPI), fine-scale bathymetric position index (FBPI), terrain ruggedness (rugosity), slope, curvature (slope of slope), profile curvature and plan curvature.Inner and outer radiuses (cells) used to calculate BBPI and FBPI were 1/10 and 1/3 respectively.1/10 was selected to identify broader-scale (2 km) seabed features such as canyons, with 1/3 selected to identify finer-scale (<1 km) features such as gullies.These are key variables that have also been used in previous deep-sea sponge models, such as Ross and Howell (2013), Ross et al., (2015), Rooper et al., (2017), andRamiro-Sánchez et al., (2019).
Seabed temperature for the study area was extracted from the Marine Institute's (Galway, Ireland) Regional Ocean Modelling System (ROMS) which has an average horizontal resolution of 1.1 km (Nagy et al., 2020).Minimum and maximum temperature data were extracted and averaged for each month over one year, 2015.From these data, monthly, average seasonal, and average annual temperature variability (maximum minus minimum temperature) were calculated as indicators for internal wave activity (Nansen, 1902;Al Senafi and Anis 2020).Monthly temperature variability was calculated as the maximum minus minimum for each month.Average seasonal and annual (averaged across 2015) temperature variability was calculated as an average of the monthly temperature variabilities in the target timeframe.Seasons were classified as follows: Spring (March, April, May), Summer (June, July, August), Autumn (September, October, November), and Winter (December, January, February).The model data were then plotted in ArcMap.Using the Inverse Distance Weighted (IDW) tool in ArcMap, temperature variability was interpolated into a 200 × 200 m raster grid using a neighbourhood size of 25 and a power of two.
Previous studies (Davison et al., 2019) that have used temperature variability as a proxy for internal wave activity have done so using archived CTD casts and oceanographic moorings.Davison et al. (2019) collated 36 CTD casts from the British Oceanographic Data Centre from 1984 to 2012 for their study region (Faroe-Shetland Channel).CTD profiles were then overlapped to create a depth profile of temperature variation/range across the study site (see Fig. 3 within Davison et al., 2019).Three 'long-term' moorings were deployed over 132 days, as well as seasonal short-term moorings (four over April-May, five over September-October).Gridded temperature variability data outputted from a ROMS was chosen for this current study over the Davison et al. (2019) approach.Gridded data from ROMS are directly relevant for habitat suitability models because they require gridded environmental data across the entire study region which would not be compatible with the Davison et al. (2019) approach, and helps to achieve the wider aims of this study (improve habitat suitability model performance).

Pheronema carpenteri data
P. carpenteri datasets were created for both density (individuals per metre) and presence-pseudo-absence. P. carpenteri count data were obtained from six video ROV transects known to contain P. carpenteri aggregations from two research cruises: i) DeepMap cruise CE15011 (2015) with ROV Holland I funded by the Eurofleets2 programme, and ii) SeaRovers CE19015 (2019) with ROV Holland I, funded by the Irish Government and European Union.Each video transect was analysed in full by recording individual P. carpenteri, along with their latitudinal and longitudinal position obtained from ROV position data.A grid was overlaid throughout the entire video transects and paused in regions of high density to reduce the potential of double-counting individuals.This produced a dataset of P. carpenteri presence along each transect.
The density dataset was calculated at the same resolution as the environmental data -200 × 200m.To account for the variation in sampling effort in each cell, an average 'density per metre ROV transect' for each cell was calculated by dividing the number of individual P. carpenteri counted in one cell by the length (metres) of transect in the corresponding cell.As a large amount of ROV transects were zero density data (no P. carpenteri present), only 12 ROV transects of zero density were randomly selected and used in model construction, thus avoiding issues of low prevalence.The final density dataset was comprised of 116 data points; 100 were of zero density and 16 were of P. carpenteri densities ranging from 0.02 to 12.95 individuals per metre ROV transect.
The presence-pseudo-absence dataset was generated by collating archived datasets (including the density transects) to create a large meta-dataset.P. carpenteri meta-data comprised archived ROV video transects and trawl data (Supplementary Material 1).The dataset was reduced to one-point-per-cell, where each point's value was either 0 (pseudo-absence) or 1 (presence) for each 200 × 200 m cell.Absence data was considered 'pseudo-absences' because of the lack of certainty and therefore unable to be 'true-absences'.Absences are not 'true absences' because only a small proportion of each 200 × 200 m cell has been sampled (e.g.ROV or trawl transect) -not each entire cell.Therefore, it cannot be said with complete certainty that each 200 × 200 m cell is truly absent of any P. carpenteri individuals.Once the archived data had been combined, the final dataset contained 657 data points (80 P. carpenteri presence and 577 pseudo-absence) spanning the study area.This method uses the same approach as in previous studies, such as Ross and Howell (2013), and Ross et al. (2015).A summary of ROV transect metadata can be found in Supplementary Material 2 (Table S2.1).

Statistical analysis
Pearson's correlation coefficients were calculated between environmental variables using the "cor()" function in the core "stats" package in R, v4.1.0(R Core Team, 2021).Highly correlated variables (<− 0.8 or >0.8) were not used when constructing multivariable models to eliminate collinearity during the variable selection process.The eliminated variable in each correlated pair was chosen during the variable selection process.
Quasi-Poisson and Tweedie generalised additive models (GAMs; Wood, 2004) were fitted to P. carpenteri density data with individual environmental variables using spline smooths to determine statistically significant (p-value < 0.05) variables.GAMs were applied using the "gam()" function from the "mgcv" package (Wood, 2004) in R. To determine any more complex non-linear relationships, as well as a result of the poor statistical performance of the GAMs, a regression-based Random Forest was also applied to the density dataset.First, the Boruta algorithm (Kursa and Rudnicki 2010) was used to guide variable selection for Random Forest by identifying any variables deemed as "unimportant" in the model for removal.Next, the "least important" variables in correlated pairs were removed.Finally, the remaining variables were removed in turnbased upon "least importance" -until the most parsimonious model was found.Once the final variables were selected, the number of variables (two) to be randomly sampled as candidates at each split was determined using the "tuneRF" function from the "randomForest" package (Liaw and Wiener 2002).The final Random Forest algorithm was then run, with 1000 trees (ntree) grown.The model residuals were tested for spatial autocorrelation using Moran's I (Moran, 1948), with a focal window of 3 × 3. Moran's I produces a p-value and z-score, where the null hypothesis is that data is randomly distributed (no spatial autocorrelation), and the alternative hypothesis is that spatial structure exists within the dataset.In the latter scenario, a positive z-score suggests data are spatially clustered (similar values are spatially grouped), compared to a negative z-score which suggests that data are competitively clustered (similar values are spatially repelled).
Similarly to the density dataset, GAMs were also initially trialled on the presence-pseudo-absence dataset.Binomial GAMs were fitted to P. carpenteri presence-pseudo-absence data with individual environmental variables using spline smooths to determine statistically significant (p-value < 0.05) variables, and explore each variable's relationship with P. carpenteri presence-pseudo-absence.Through forward-and backwards-step selection, a final presence-pseudo-absence GAM was selected using the individually statistically significant variables.The final GAM residuals were tested for spatial autocorrelation using Moran's I (Moran, 1948), with a focal window of 3 × 3. To account for any detected spatial autocorrelation, a residual autocovariate (RAC) variable was calculated from the GAM residuals (Crase et al., 2012).The final model was re-run with the addition of the RAC but the resultant model was not statistically significant, nor were any model terms.For thoroughness, forward-and backwards-step selection was re-run to include the RAC, and again, no statistically significant model could be found.
As a result, a random forest algorithm was fitted to the presencepseudo-absence datasets using the same methodology as the density dataset.After variable selection and once a preliminary random forest model had been fitted, Moran's I was calculated on the model residuals and a new RAC value was calculated (Crase et al., 2012) from the model residuals.Subsequently, the preliminary random forest model was re-fitted with the addition of the RAC.A series of different focal sizes were trialled (3 × 3, 5 × 5, 7 × 7, 9 × 9) when calculating the RAC to strike a balance between effectively reducing spatial autocorrelation and not creating an overwhelmingly powerful predictor variable (Georgian et al., 2019).Moran's I was then recalculated on the final model's residuals to determine whether spatial autocorrelation had successfully been accounted for with the inclusion of the RAC variable.

Correlations
Pearson's correlation coefficient was calculated between all environmental variables, including depth (bathymetry), broad and fine bathymetric position index (B/FBPI), slope, rugosity, curvature, profile curvature, plan curvature and average annual, seasonal (spring, summer, autumn, winter) and monthly (Jan-Dec) temperature variability.Results suggest that within the density dataset, all temperature variability variables, as well as curvature-plan curvature, and curvatureprofile curvature are highly correlated (− 0.8 < or >0.8).Within the presence-pseudo-absence dataset, the correlated pairs were the same as in the density dataset, with the addition of BBPI-FBPI and profile curvature-plan curvature.A figure of correlation results can be seen in Supplementary Material 3 (Figure S3.1).

Pheronema carpenteri density and environmental variables
The overall results from fitted Quasi-Poisson and Tweedie GAMs suggested no relationships between P. carpenteri density and any of the individual environmental variables.The results of individually fitted GAMs (Supplementary Material 4) were either not statistically significant, or when plotted, the fitted models did not make ecological sense.As a result, an attempt was not made to fit a multi-variable GAM to the density dataset.In light of the GAM results, a regression-based Random Forest algorithm was applied to the density set to either (1) confirm the previous GAM results, or (2) fit more complex relationships.
For the selection of variables and model fitting, average annual temperature variability was used as the temperature variability component in the models because all temperature variability variables were correlated so could not be used in the same model.Mean annual temperature variability was selected because of its broader temporal scale, therefore better representing the prevailing oceanographic conditions and reducing any seasonal influences.The Boruta algorithm deemed curvature, plan curvature and profile curvature as unimportant terms in the model; therefore, these variables were removed.In addition, FBPI and average annual temperature variability were removed because of their low importance, which maximised model performance.The final variables selected were therefore depth, slope, BBPI and rugosity.The results of the final Random Forest model were poor; mean of square residual = 2.01, variance explained = 15.86%.Thus confirming, along with the GAM results, that there is no relationship between P. carpenteri density and the environmental variables used in this study at the scale of use.Moran's I was calculated on the residuals of the random forest model to test for spatial autocorrelation.Moran's I was not statistically significant (p-value = 0.95) and therefore the null hypothesis was accepteddata is randomly distributed and therefore no spatial autocorrelation has been detected.

Pheronema carpenteri presence-pseudo-absence and environmental variables
Seabed depth sampled for the presence-pseudo-absence dataset range from 206 to 3110m.Within this, P. carpenteri distribution is not homogenous and was observed strictly between 974 and 1329m, a 355m range (Fig. 2).These regions were characterised by intermediate temperature variability (0.6-1.2 • C), with pseudo-absence data distributed across all temperature variability ranges (0.07-2.2 • C) (Fig. 2).
GAMs fitted to all individual environmental variables were significant (Supplementary Material 5) and therefore all considered in the model selection process.Individual fitted GAMs characterise the niche that P. carpenteri is occupying as regions of low rugosity (smooth surface roughness and low complexity) and gentle slopes.Forward-and backwards-step variable selection both produced the same final, bestperforming model containing depth (bathymetry), FBPI and average annual temperature variability (Table 1).All variables were fitted with spline smooths and a knot value of three.When Moran's I was calculated using the residuals of the selected GAM, spatial autocorrelation was detected (p-value = < 0.001, z-score = 0.64).Subsequently, new GAMs were fitted to the data including a RAC variable calculated from the preliminary GAM; however, no statistically significant model could be produced that included RAC.Subsequently, random forest models were fitted to the presence-pseudo-absence data.
The Boruta algorithm determined all variables as "important" terms in the model.Following the removal of correlated variables and stepwise removal, the most parsimonious model was found.The final variables selected were average annual temperature variability, depth, FBPI and rugosity.The results of the Random Forest model were good (Table 2).Spatial autocorrelation was detected in the residuals of the Random Forest model (p-value = < 0.001, z-value = 0.43).Subsequently, a new RAC variable was calculated from the Random Forest residuals.A final Random Forest (RAC-RF) was fitted with the preliminary and RAC variables; model performance markedly increased (Table 1).The most important variable was RAC, followed by rugosity, annual temperature variability, depth and FBPI.The addition of the RAC variable successfully accounted for spatial autocorrelation when Moran's I was calculated (p-value = 0.99).See Fig. 3 for partial dependence plots for each RAC-RF predictor.

Temperature variability within the study area
Across the study area, there are distinct spatial differences in temperature variability throughout the Porcupine Seabight.As expected, steeper sloping regions within the study area coincide with regions of higher temperature variability, particularly in the North-East and North-West of the Porcupine Seabight, along with Porcupine Bank and the Southwest Approaches.Regions of low variability are the continental shelf, the north and south ends of the Porcupine Seabight and Goban Spur, with the Porcupine Abyssal Plain characterised by negligible temperature variability.Regions of higher temperature variability calculated from the ROMS coincide with regions described by Rice et al. (1990) where the "ray slope" is exceeded by the slope of the seafloor (see Fig. 2. In Rice et al., 1990), and therefore in the regions where near-bed enhanced currents are likely to occur.The assumptions of Rice et al. (1990) were confirmed by White (2003) using in-situ current measurements.Additionally, White (2003) evidenced that the enhanced near-bed current speeds are a result of baroclinic (internal) tidal motions that may be enhanced by reflecting internal (M 2 ) waves.Furthermore, areas of high-temperature variability in the northeast of the Porcupine Seabight (Fig. 1) coincide with documented internal wave activity.Within this region, internal waves occur at the interface between Mediterranean Outflow Water (MOW) and Eastern North Atlantic Water (ENAW).The interface forms a diffuse transition zone between 700 and 900 m that is marked by nephaloid layers (Dickson and McCave 1986), which include particulate organic matter.This region is characterised by strong residual near-seabed current flows (up to 15 cm s − 1 ; Mohn et al., 2014) but is also further enhanced by breaking internal waves and tides.As a result, strong along-slope transport of sediment and organic matter fluxes occurs (Wienberg et al., 2020 and references therein).The Belgica Mound Province, a collection of cold-water coral mounds, occurs within this region of high-temperature variability in the northeast of the Porcupine Seabight.Here, internal waves propagating along the MOW-ENAW interface enhance the energy supply and particle flux in the area (Wienberg et al., 2020), thus shaping the distribution of the cold-water coral mounds.

Table 1
Evaluation metrics for GAM fitted to presence-pseudo-absence data, including overall adjusted R-squared (R-SQ), deviance explained (D), un-biased risk estimator (UBRE) score, and model term (variable) estimated degrees of freedom (EDF) and p-values.There are, however, potential limitations to using temperature variability as a proxy for internal wave activity in our study.Other oceanographic processes may attribute variability in temperature, for example, the breakdown and re-stratification of the seasonal thermocline.However, P. carpenteri is primarily distributed deeper than the permanent thermocline and winter mixed layer (White and Dorschel 2010).Therefore, at these depths, temperature variability due to seasonal stratification is likely to be minimal.Additionally, seasonal changes in water mass signatures, such as the Mediterranean Outflow Water (García Lafuente et al., 2007), may attribute to some temperature variability.However, the inclusion of tides in the ROMS, using 'average annual temperature variation' to reduce the seasonal influence, and overlap between critical slope angle, high-temperature variability and known internal wave activity in the study area means the authors have confidence that the temperature variability in this study area is largely attributed by internal wave activity.

Pheronema carpenteri presence-pseudo-absence
The depth distribution of P. carpenteri is confined to a small bathymetric range, 974-1,329m.This distribution pattern of being bathymetrically confined has been observed in other aggregating deep-sea sponges, including Ostur in the Faroe-Shetland Channel (Bett 2001;Davison et al., 2019, Kazanidis et al., 2019) and the wider North-East Atlantic (Klitgaard and Tendal 2004).The depth distribution of P. carpenteri in the Porcupine Seabight falls mainly below the base of the permanent thermocline (600 to 1000 m), which is associated with enhanced near-bed currents two to three times stronger than flows at other depths (White and Dorschel 2010).This depth-current relationship is reflected in P. carpenteri distribution within regions of intermediate temperature variability.These relationships suggest P. carpenteri is distributed deeper, but still proximal to regions of enhanced currents and therefore higher temperature variability (Fig. 2).Average annual temperature variability is the second most significant term within the full GAM after depth, suggesting that it is potentially an important driver of P. carpenteri distribution within the study area.GAM results support the conclusions of Rice et al. (1990) that P. carpenteri is not found in regions of enhanced currents.This is probably due to the morphological features of P. carpenteria weak spicule anchor system on soft, fine sediments means P. carpenteri is unlikely to be able to stay anchored in areas of high current speeds.
From proximal distances, P. carpenteri could still exploit an increase in suspended food particles associated with enhanced currents.Internal waves and tides are known to generate and maintain bottom and nearbottom nepheloid layers (Cacchione and Drake 1986;Puig et al., 2004), which are "zones of increased turbidity relative to surrounding waters" (McCave 2009).Nepheloid layers form both bottom and intermediate layers in the Porcupine Seabight.These contain dissolved and particulate matter of terrigenous and biological origin (Vermeulen 1996), have residence times in the order of weeks and months (Gardner et al., 1985) and are known to be associated with cold-water coral reefs (Mienis et al., 2007).It is plausible that nepheloid layers and other suspended material generated within enhanced current regions could be transported "downstream" to P. carpenteri aggregations in the Porcupine Seabight (White 2003).In the northern region of the Porcupine Seabight, the cyclonic (anti-clockwise) and down-slope current (Pingree and Le Cann 1990) could transport suspended material, including food particles, westward and down-slope ( Wienberg et al., 2020) where there are known P. carpenteri aggregations.Although current and suspended particulate data would be required to test this hypothesis.However, an association between internal waves/tides, nepheloid layers and Geodia spp.distribution was also drawn by Roberts et al. (2021).

Pheronema carpenteri density
Within the Porcupine Seabight, the density of P. carpenteri has no relationship with any of the environmental variables at 200 m resolution.In comparison, substrata have been found elsewhere to drive sponge densitya pattern also observed across the macrobenthos more generally in the Faroe-Shetland Channel (Bett 2001).Maps of substrata at resolutions relevant to habitat suitability (~200 m) are not available for our study areaas well as more generallyand therefore could not be considered.Klitgaard et al. (1997) and Bett (2003) also first noted a relationship between the distribution and density of sponge grounds in the Faroe-Shetland region with internal wave processes.More recently, this theory has been followed up by studies from Davison et al. (2019) and Kazanidis et al. (2019).They found that in the Faroe-Shetland Channel sponge density peaked within a narrow bathymetric range of 80 m.This distribution is likely due to the depth specific hydrographic setting of the channel; more specifically, the strong stratification that occurs between the North Atlantic and Arctic water masses.The physical interaction of these water masses along the slope generates internal waves that would then concentrate and supply particulate organic matter (food) to the aggregating sponges.Davison et al. (2019) demonstrated a positive linear relationship between the density of Geodia sp.(Demospongiae) and temperature variability (used as a proxy for internal waves) in the Faroe Shetland Channel (FSC).They concluded that high densities of Geodia sp.occupy regions of high internal wave activity in the FSC.As discussed by Davison et al. (2019), this is likely to be a result of Geodia sp.morphology; attaching directly onto hard substrata, and often found growing over and incorporating gravel/cobbles into their base.Davison et al. (2019) results suggest that Geodia sp. can withstand the high current and turbulence associated with internal wave activity, and therefore benefit from being in regions of increased availability of food.The results of Davison et al. (2019) differ significantly from our study and is therefore apparent that the relationship between sponge density and regions of temperature variability is not uniform across different sponge taxa.These relationships need further examination, including from new study areas outside of the Porcupine Seabight and Faroe-Shetland Channel.
It is of notable interest that the environmental variables used in this study do not display relationships with P. carpenteri density, but do with distribution.The results, therefore, suggest that these environmental variables are important characteristics of the niche P. carpenteri occupy, but are not drivers of dense aggregations, or at least at a 200 × 200m resolutionthere may be a mismatch between the resolutions of environmental and biological data (Rengstorf et al., 2014;Ross et al., 2015;Dolan et al., 2021).It is also important to note that other variables not used in this study could be drivers of density.For example, biological, physiological and ecological processes and interactions may influence density but are not captured by environmental variables.
The reproductive strategy of P. carpenteri is currently unknown but Hexactinellids are assumed to reproduce both asexually and sexually with lecithotrophic larvae (Leys and Ereskovsky 2006;Teixidó et al., 2006).Asexual reproduction occurs via budding or by the fragmentation of individuals (Teixidó et al., 2006).Both of these asexual strategies would have limited dispersal capability so may result in the retention of new buds or fragments close to the parent, thereby increasing local density.This process may also be highly localised/occur on small scales.As a result, there is likely to be a 'mismatch' between the spatial scales at which these ecological processes occur and environmental data is available.The observations of Barthel et al. (1996) describe small juveniles at the periphery of P. carpenteri patches off Morocco.This suggests that there may be some retention mechanisms of either buds, fragments or larvae within the aggregation or persistent recruitment of larvae into the aggregation patch.Larval dispersal models from within the Porcupine Seabight study area suggest that known P. carpenteri aggregations are connected via single generational dispersal events (Ross et al., 2019).Therefore, these aggregations are likely to be part of the same meta-population, and persistent recruitment between patches is plausible.However, until genetic studies have better established P. carpenteri population dynamics and reproductive strategy, this is a hypothesis.Additionally, there is a wider need for more research in understanding deep-sea sponge larvae behaviour (Gary et al., 2020) and the role of oceanography/hydrography (Busch et al., 2021;Ross et al., 2020) in shaping the distribution and settlement of larvae in the deep sea.Without this fundamental knowledge, models of larval distribution and connectivity will remain limited in their accuracy and reliability.
As well as reproductive strategy, feeding strategy is also likely to play a role in shaping the distribution and density of P. carpenteri, despite having a low microbial abundance (Colaço et al., 2020).Bart et al. (2020) assessed the assimilation of dissolved organic matter (DOM) by sponges and their associated bacterial community.In this study, Bart et al. (2020) found that the low microbial abundance sponge species (Hymedesmia paupertas and Vazella pourtalesii) had a four-five times higher uptake rate of DOM compared to high microbial abundance sponges (Geodia barretti) and that hexactinellids assimilated DOM more efficiently and at a higher rate compared to demosponges.Kazanidis et al. (2018) also demonstrated how Spongosorites coralliophaga has a versatile and flexible feeding strategy, assimilating carbon and nitrogen by all four food sources ( 15 N-ammonium chloride, 13 C-glucose, 13 C/ 15 N-labelled microalgae, 13 C/ 15 N-labelled bacteria) tested in the study.Kazanidis et al. (2018) conclude that the feeding strategy of S. coralliophaga plays an important role not just in the sponge's survival in food-limited conditions, but also its formation of dense aggregations with high biomass.Given the results of Bart et al. (2020) and Kazanidis et al. (2018), the high density of P. carpenteri observed within the Porcupine Seabight and Goban Spur may also be a result of its feeding strategy.

Consequences for the development of density-based habitat suitability
modelling to support spatial management P. carpenteri aggregations are considered both a Vulnerable Marine Ecosystem (VME) under UNGA Resolution 61/105 and a threatened and declining habitat under the OSPAR convention.The ability to predict and locate areas of suitable habitat for this species is critical for effective management.However, as it is aggregations of the species that are considered VME, not individual occurrences, we must move towards predicting the areas of high densities of individuals that would constitute an aggregation.
There are few published presence-absence HSMs of P. carpenteri (Ross and Howell 2013;García-Alegre et al., 2014;Ross et al., 2015;Howell et al., 2016;Howell et al., 2022), all of which have been built across varying spatial extents and resolutionfrom low-resolution basin-wide models, e.g.Howell et al. (2016), to high resolution local or regional models, e.g.García-Alegre et al. (2014).The explanatory variables used in these models are limited to topographical variables, temperature and depth, with the addition of some chemical variables, such as dissolved oxygen and silicate concentration, used by Howell et al. (2016) at a low resolution (0.016 • ).However, these published models do not incorporate any physical oceanographic/hydrographic variables or use a proxy such as temperature variability.The results of this study show that there is a strong relationship between P. carpenteri distribution and temperature variability.This is likely due to temperature variability acting as a proxy for internal wave activity, which subsequently plays a role in the food supply.Therefore, where such data are obtainable, temperature variability, or a direct metric for internal wave activity, should be used as an explanatory variable in HSMs to obtain accurate predictions of P. carpenteri distribution (presence-absence).However, the accessibility of this type of data may be limited.Global and basin-wide models, such as the HYbrid Coordinate Ocean Model (HYCOM), now incorporate internal tides (Arbic et al., 2012) but are restricted to coarse resolutions.Conversely, high-resolution models incorporating internal tides will be geographically limited to regional/local study areas.However, through the UN's Decade of Ocean Science and its affiliated programmes, e.g. the Deep Ocean Observing Strategy (https://www.deepoceanobserving.org/), efforts are being made to better characterise, model, and predict the physical oceanography of the deep ocean.Publicly available biological/distribution data is also lackinga critical component needed for not only building/training HSMs but also for independently validating HSMs (Anderson et al., 2016;Rooper et al., 2016Rooper et al., , 2018;;Howell et al., 2022).Where possible, distribution and density data on sponge aggregations should be made publicly available, such as through the ICES VME Database (https ://www.ices.dk/data/data-portals/Pages/vulnerable-marine-ecosystems.aspx) for consideration by the Working Group on Deep-Water Ecology (ICES, 2021).Furthermore, a more coordinated approach is required across the international deep-sea community to better target biological observing efforts.The UN Ocean Decade affiliated programme, Challenger 150 (https://challenger150.world/), sets out a blueprint for how a global-scale field programme can achieve this through standardising sampling methods, with targeted sampling across biogeographic regions, depth and latitude (Howell et al. 2020a(Howell et al. , 2020b)).As a result of these existing programmes (DOOS, Challenger 150), higher resolution, more powerful HSMs may be possible through the increased availability of high-resolution environmental variables and biological/distribution data.
Our findings suggest there is currently a limited capability to predict the density of P. carpenteri, and thus identify areas of dense aggregations that can be considered VMEparticularly at the fine-scale resolution required for area-based management, e.g.≤ 200 m.It is likely that environmental drivers that operate on much finer spatial and temporal scales than traditionally measured/modelled, as well as environmental drivers not commonly recorded, play a role.Additionally, there are mismatches in coarse resolution of oceanographic model outputs, e.g.0.25 • , at global or basin scale, compared to the biology which operates on much finer scales in the order of metres.Although coarse oceanographic data is useful for broad-scale habitat mapping, it is not useful in helping to understand and predict the distribution/density of sponges at the fine scale required for area-based management.However, targeted investigation and comparison of areas of modelled suitable habitat, ranging from zero to high density of sponges, would help determine key drivers and their scales of operation and inform the development of new, and measurement of existing, essential ocean variables over the next decade.
Additionally, a general lack of knowledge of P. carpenteri biology and physiology makes it challenging to understand the possible role these traits may have in driving the density of P. carpenteri.In particular, a better understanding of the reproduction, recruitment, and population dynamics of both local, and wider P. carpenteri meta-population is needed.This requires new data on known aggregations but also from new, undiscovered aggregations.The discovery of undocumented aggregations can be facilitated by the development of more accurate HSMs that incorporate temperature variability.With further understanding, existing modelling methods, such as the combination of HSM and dispersal modelling, trait-based modelling, or even novel methods, may enable the successful prediction of P. carpenteri density at a fine-scale resolution.

Conclusion
This study has shown that as individual variables and within multivariable models, temperature variability, depth and FBPI are significant drivers of Pheronema carpenteri distribution (presence-pseudoabsence) within the Porcupine Seabight.It would therefore be appropriate for future studies to incorporate temperature variability, as a proxy for internal wave activity, into habitat suitability modelling to build better performing models and make more accurate predictions of P. carpenteri distribution.On the contrary, this study has also shown that none of the environmental variables used has any relationship with P. carpenteri density, or at least at 200 × 200 m resolution.This may reflect a mismatch in the spatial scales in which environmental and biological processes are occurring.Biological processes likely occur at scales much finer than the resolution that environmental variables are traditionally observed/modelled.This study also establishes that the previously reported relationship between internal wave activity and sponge density is not universal across sponge taxa and that the role of internal waves in shaping sponge density requires further study.Furthermore, this study identifies a knowledge gap in the underlying drivers of P. carpenteri density.A greater understanding of this species is needed to better model species density and ultimately provide the data to support the effective management of this habitat-forming species.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig. 2. Presence (light blue) and pseudo-absence (dark blue) distribution of P. carpenteri with respect to depth (m) and average annual temperature variability ( • C).The red-bounding box represents the permanent thermocline.

Table 2
Comparison of Random Forest model performance, with and without the RAC variable, broken down by variance explained (%) and mean of square residuals.Additionally, results of Moran's I by model.