Climate change winner in the deep sea? Predicting the impacts of climate change on the distribution of the glass sponge Vazella pourtalesii

Shallow-water sponges are often cited as being ‘climate change winners’ due to their re siliency against climate change effects compared to other benthic taxa. However, little is known of the impacts of climate change on deep-water sponges. The deep-water glass sponge Vazella pourtalesii is distributed off eastern North America, forming dense sponge grounds with enhanced biodiversity on the Scotian Shelf off Nova Scotia, Canada. While the strong natural environmental variability that characterizes these sponge grounds suggests this species is resilient to a changing environment, its physiological limitations remain unknown, and the impact of more persistent anthropogenic climate change on its distribution has never been assessed. We used Random Forest and generalized additive models to project the distribution of V. pourtalesii in the northwest Atlantic using environmental conditions simulated under moderate and worst-case CO2 emission scenarios. Under future (2046–2065) climate change, the suitable habitat of V. pourtalesii will increase up to 4 times its present-day size and shift into deeper waters and higher latitudes, particularly in its northern range where ocean warming will serve to improve the habitat surrounding this originally sub-tropical species. However, not all areas projected as suitable habitat in the future will realistically be populated, and the re duced likelihood of occurrence in its core habitat on the Scotian Shelf suggests that the existing Vazella sponge grounds may be negatively impacted. An ef fective monitoring programme will require tracking changes in the density and distribution of V. pourtalesii at the margins between core habitat and where losses and gains were projected. Species distribution models showed that the suitable habitat of the deep-water glass sponge Vazella pourtalesii will ex pand under climate change scenarios. Photo: Fisheries and Oceans Canada


INTRODUCTION
Over the past century, upper-ocean temperatures have increased across much of the northwest Atlantic, consistent with the global trend of increasing sea surface temperature due to anthropogenic climate change (Bush & Lemmen 2019). The effects of ocean warming have already been observed in marine ecosystems ABSTRACT: Shallow-water sponges are often cited as being 'climate change winners' due to their resiliency against climate change effects compared to other benthic taxa. However, little is known of the impacts of climate change on deep-water sponges. The deep-water glass sponge Vazella pourtalesii is distributed off eastern North America, forming dense sponge grounds with enhanced biodiversity on the Scotian Shelf off Nova Scotia, Canada. While the strong natural environmental variability that characterizes these sponge grounds suggests this species is resilient to a changing environment, its physiological limitations remain unknown, and the impact of more persistent anthropogenic climate change on its distribution has never been assessed. We used Random Forest and generalized additive models to project the distribution of V. pourtalesii in the northwest Atlantic using environmental conditions simulated under moderate and worst-case CO 2 emission scenarios. Under future (2046-2065) climate change, the suitable habitat of V. pourtalesii will increase up to 4 times its present-day size and shift into deeper waters and higher latitudes, particularly in its northern range where ocean warming will serve to improve the habitat surrounding this originally sub-tropical species. However, not all areas projected as suitable habitat in the future will realistically be populated, and the re duced likelihood of occurrence in its core habitat on the Scotian Shelf suggests that the existing Vazella sponge grounds may be negatively impacted. An ef fective monitoring programme will require tracking changes in the density and distribution of V. pourtalesii at the margins between core habitat and where losses and gains were projected. across the region, through northward shifts in the range of commercially harvested fish and their catch distribution (Nye et al. 2009, Pinsky & Fogarty 2012. Given the inevitable socio-economic impacts of a changing climate on the valuable fisheries that operate in the northeast USA and on the Scotian Shelf in Canadian waters, several studies have aimed to identify particular geographic areas of commercially harvested species considered vulnerable to present and future climate change (e.g. Hare et al. 2016, Rheuban et al. 2018. In the absence of information on physiological thresholds, species distribution models (SDMs) are often used to estimate potential climate change effects by predicting the degree of habitat loss or gain under various emission scenarios and time periods, based on the statistical association between occurrence records and spatial environmental data. Such correlative tools have been recently utilized to predict changes in the future distribution of commercially harvested groundfish (e.g. Kleisner et al. 2017, Stanley et al. 2018, McHenry et al. 2019) and motile invertebrate species (e.g. Tanaka et al. 2017, Stanley et al. 2018 in the northwest Atlantic. However, much less is known about the impacts of climate change on the distribution of long-lived, sessile benthic invertebrates of conservation, but not commercial, importance in the region (although see Morato et al. 2020). Such communities may be more vulnerable than motile species, as latitudinal shifts in distribution are unlikely to occur over short time scales.
Sponges are key components of benthic marine ecosystems worldwide, from intertidal to abyssal depths (Maldonado et al. 2017). In the deep sea, sponge species can be widely distributed across vast areas as isolated individuals, but in certain locations where environmental conditions are favourable, they may form dense aggregations commonly known as sponge grounds. The importance of such aggregations in deep-sea environments has only recently emerged, through demonstrations of their role in biodiversity enhancement and habitat provision (Beazley et al. 2013, Hawkes et al. 2019, Murillo et al. 2020, and their importance in benthic−pelagic coupling and the cycling of nutrients (Kutti et al. 2013, Pham et al. 2019. The Scotian Shelf off Nova Scotia, eastern Canada, is home to a globally unique aggregation of the glass sponge Vazella pourtalesii. This rosellid species has a wide but fragmented distribution along the continental margin of eastern North America betweeñ 100 and 750 m depth (Beazley et al. 2018, NOAA 2019, occurring in low densities off Florida in the southeastern USA, and forming the densest known monospecific aggregations of their kind in Emerald Basin on the Scotian Shelf. Through an analysis of in situ benthic imagery inside the 'Vazella sponge grounds' (sensu Beazley et al. 2018), Hawkes et al. (2019) found that on average, mean species density (i.e. richness) and abundance of other epibenthic megafauna was ~3−4 times higher in the presence of V. pourtalesii, demonstrating the importance of these sponge grounds in enhancing local biodiversity. Beazley et al. (2018) used classification Random Forest modelling to predict the present-day distribution of V. pourtalesii across the Scotian Shelf. A suite of 35 predictor variables relating to ocean terrain, fishing effort, and biological and physical oceanographic characteristics considered to represent the current oceanographic climate were utilized. Additionally, historical trends in bottom temperature and salinity, 2 variables that were identified as important determinants of its present-day distribution, were examined through the hindcast Simple Ocean Data Assimilation model, revealing that the Vazella sponge grounds have historically been subjected to strong inter-annual and multi-decadal variability in water mass characteristics, with annually averaged bottom temperatures varying by up to 8°C (range 4−12°C). This variability was related to periods of larger than average transport of cold, fresh Labrador Slope Water by the Labrador Current, consistent with the influence of the Atlantic Multi-Decadal Oscillation (AMO), a mode of natural variability in the North Atlantic responsible for warming and cooling during positive and negative phases of the Atlantic Meridional Overturning Circulation, respectively. The persistence of V. pourtalesii in Emerald Basin since its discovery in the late 1800s (Honeyman 1889) suggests that this species may be inherently adapted to a highly dynamic physical environment, and able to withstand the particularly pronounced cold conditions that prevailed in Emerald Basin and across the Scotian Shelf in the 1960s (Beazley et al. 2018). While the impact of persistently warmer conditions that exceed the upper thermal range of the AMO (~12°C; Beazley et al. 2018) on V. pourtalesii has not been quantitatively assessed, the presence of this species in the subtropical waters off the mid-and southeastern USA suggests it is adapted to withstand a warmer climate. However, its ability to adapt to such changes over relatively short time scales and the physiological tolerance of individuals and populations in different parts of its range, remain unknown.
The response of sponges to ocean warming and acidification is well studied in shallow-water ecosys-tems (see overviews in Carballo &Bell et al. 2018). In these environments, sponges have been referred to as potential 'climate change winners' (Bell et al. 2013(Bell et al. , 2018, due mainly to the physiological tolerance of some species to high thermal anomalies that have caused bleaching events in corals. Observations from several experimental and fieldbased studies suggest that some shallow-water sponges have an upper thermal tolerance of ~2−3°C above mean monthly values (Bennett et al. 2017, Ramsby et al. 2018. However, virtually nothing is known of the physiological effects of climate change on glass sponges, which are thought to be relatively sensitive to changes in temperature (Leys & Meech 2006), and how their physiological thresholds may drive community or ecological responses (e.g. distribution shifts) under a changing climate. In an atmospheric CO 2 doubling experiment, upper ocean (0−300 m depth) temperatures on the northwest Atlantic shelf between 35 and 45°N were predicted to warm by ~3°C (Saba et al. 2016), a rate 2−3 times faster than the global average. This change may have important implications for the distribution of the Vazella sponge grounds and the biodiversity that they support. Furthermore, understanding the effects of climate change on this vulnerable marine ecosystem indicator species (Fuller et al. 2008) will be necessary to determine whether the 2 bottom-fishery closures currently in place on the Scotian Shelf for the conservation of this species re quire more adaptive approaches to ensure their effectiveness into the future.
In the absence of information on the physiological tolerance of this long-lived, deep-water glass sponge, we here aim to deduce the potential impacts of future climate change on the distribution of V. pourtalesii using correlative SDM techniques. We employed the machine-learning method Random Forest, and generalized additive modelling (GAM) to project the distribution of this species in the northwest Atlantic using environmental variables simulated for the bi-decadal period 2046−2065 under moderate (representative concentration pathway [RCP] 4.5) and worst-case (RCP 8.5) CO 2 emission scenarios. This period was chosen instead of the commonly modelled end-of-century (2100) to allow for consideration of the results in current conservation management strategies. We assessed whether this long-lived, deep-water sponge species is a potential 'winner', 'loser', or is neutral to the effects of climate change based on whether its suitable habitat expanded, contracted, or showed no change under future climatic conditions. We assume that V. pourtalesii is locally adapted to the unique water mass characteristics that influence the northern and southern portions of its range, and have applied models separately to the occurrences located on the Scotian Shelf in Canadian waters and those off the mid -southeast USA. Collectively, our ap proach will provide a better understanding of the potential upper thermal tolerance limit of this species and how it may affect its distribution in the northwest Atlantic into the future. To our knowledge, our study represents the first application of SDMs to evaluate climate-induced changes in the distribution of a deepwater, ground-forming sponge, and has important implications for the conservation of these relatively understudied ecosystems.

Distribution of Vazella pourtalesii and study area
Occurrences of V. pourtalesii in the northwest Atlantic ( Fig. 1) come from a variety of sources (Table 1). On the Scotian Shelf (see Breeze & Horsman 2005), significant effort has been made to map the distribution of the Vazella sponge grounds for conservation management purposes, and as a result, the occurrences ( Fig. 1) are spatially biased towards this part of its range. There, records are collated from the Fisheries and Oceans Canada (DFO) multispecies research vessel trawl survey, DFO and Natural Resources of Canada (NRCan) optical (camera/video) benthic surveys (comprising the majority of observations; Table 1), and commercial bycatch records from the Fisheries Observer Program. The DFO multispecies research vessel trawl survey is a depthstratified random survey conducted primarily using Western IIA trawl gear. Tows are 30 min in duration, and the average linear distance covered is 3.17 km (Beazley et al. 2018). Commercial trawls in the region are often much longer in duration and length, collating catch data over distances of 10+ km. For both the re search vessel and commercial trawl data, the start position of each tow was used to represent the presence of V. pourtalesii in the catch.
Records located off the mid-southeast USA (Table 1) were extracted from a combination of sources including the primary literature and online repositories, and through direct video annotations collected for the purpose of this study. Three specimens collected in 1868 off Florida that were re-examined by Reiswig (1996) are included here, as well as remotely operated vehicle (ROV) and submersible records extracted from the NOAA Deep-Sea Coral Data Portal (https://deepseacoraldata.noaa.gov/). In both 2018 and 2019, NOAA conducted 2 ROV surveys off the mid-and southeastern USA ('Windows to the Deep' Okeanos Explorer EX1806 and EX1903 missions in 2018 and 2019, respectively) and reported the presence of V. pourtalesii. In order to increase the number of presences in this relatively under-sampled (compared to the Scotian Shelf) portion of its range, all dives within its expected depth range (<1000 m) from these 2 missions were examined and presences recorded using the Ocean Networks Canada SeaTube V2 annotation software (https://data.oceannetworks. ca/ SeaTubeV2). A total of 1051 V. pourtalesii individuals were observed from 9 different stations (see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m657 p001 _ supp .pdf, for details). Most observations were from the Blake Plateau where V. pourtalesii was found in localized, dense aggregations in association with the reef-building coral Lophelia pertusa (and its rubble). These observations ex tend the maximum depth of V. pourtalesii in the northwest Atlantic to 935 m (from 754 m previously reported from observations on Cape Canaveral North off Florida; NOAA 2019). Finally, a search for V. pourtalesii records not encompassed by the sources described above was made in the Ocean Biogeographic Information System (https://obis.org/) to ensure that all records were considered. No new records resulted from this search. Given the wide latitudinal gradient and presumed environmental heterogeneity over which V. pourtalesii is distributed, key environmental variables (depth, bottom temperature, and salinity) shown to be important for the distribution of this species (Beazley et al. 2018) were evaluated at the V. pourtalesii presence localities in order to assess the potential generality and future transferability of a model trained on its entire range (see Wenger & Olden 2012). Temperature−salinity plots and bivariate plots of the relationships between depth− temperature and depth− salinity at the location of the occurrences were examined (see Fig. S1). At comparable depths, records of V. pourtalesii on the Scotian Shelf were located in a colder, less saline environment than those located off the mid -southeast USA. We therefore chose to model the V. pourtalesii occurrences in the northern portion its range separately from its southern distribution. Two modelling subareas were created (see Fig. 1) and split at 40°N, which roughly corresponds to the boundary of the 2 regimes of the Ekman portion of the Atlantic Meridional Overturning Circulation (see Fig. 9 of Wang et al. 2019), reflecting the anticyclonic and cyclonic atmosphere circulations over the sub-tropical and subpolar gyres. The subarea north of 40°N is referred to as the 'Northeast US/ Atlantic Canada' subarea, which in cludes the northeast USA and Atlantic Canada regions, but also the Flemish Cap situated in international waters off New foundland, Canada (Fig. 1). The spatial extent south of 40°N is referred to as the 'Mid -Southeast US' subarea. Both subareas extend seaward to 2000 m depth (based on the resampled General Bathymetric Chart of the Oceans [GEBCO] 2019 bathymetry layer), which was the upper depth limit of model extrapolation based on previous modelling exercises of this species (see Beazley et al. 2016Beazley et al. , 2018. A 5 km land buffer was applied to all coastline within each subarea, and any environmental grid cell (see Section 2.2, below) falling within its boundaries was removed to avoid the inclusion of land points.

Environmental predictor layers
A total of 10 environmental predictor layers (Table 2), including 2 static (stable in time, i.e. depth, slope) and 8 dynamic (varying in time, e.g. temperature, salinity) layers were extracted for consideration in the models. Although previous modelling works (Beazley et al. 2018) highlighted the importance of both hydrographic conditions (e.g. bottom and surface temperature, bottom salinity) and ocean productivity (e.g. sea surface primary production) in predicting the distribution of this species, only a limited  . Monthly data for the time period 1990−2015 were extracted from BNAM and averaged across months and years to represent present-day mean climatic conditions. A version of BNAM was developed to simulate future climatologies for the bi-decadal period 2046− 2065 (referred to hereafter as the '2055 climatology') under 2 different IPCC (IPCC 2013) emission scenarios: RCP 4.5 and RCP 8.5 (Brickman et al. 2016). RCP 4.5 is considered an emission-stabilizing scenario in which radiative forcing increases to 4.5 watts m −2 and stabilizes around 2100, when trajectories of atmospheric CO 2 concentration and the median global temperature anomaly reach ~650 ppm and 2.4°C, respectively (Moss et al. 2010). RCP 8.5 is a highemission scenario in which radiative forcing in creases to 8.5 watts m −2 by 2100 (~1370 ppm CO 2 and a median temperature anomaly of 4.9°C; Moss et al. 2010), but does not stabilize until 2300 when it reaches 12 watts m −2 (van Vuuren et al. 2011). Averaged annual anomalies representing the difference between present day and future (2055) conditions were extracted from BNAM simulations (see Brickman et al. 2016) and applied (added) to the averaged present-day climatology layers modelled in this study to create predictor layers representing future forecasted conditions. Depth was extracted from the 2019 GEBCO grid (https://www.gebco.net/data_and_products/gridded_ bathymetry_data/gebco_2019/gebco_2019_info.html). Slope was derived from this bathymetry using the 'Slope' tool in ArcGIS Desktop 10.7 (ESRI 2019). Prior to the generation of slope, depth was projected in an Albers equal area conic projection where units are represented in metres.
Continuous raster surfaces of the BNAM point data were created in ArcGIS using the 'Point to Raster' tool, with point data displayed using geographic coordinates and a WGS 1984 datum. The final raster cell size was 0.088°, which is approximately equivalent to 7.5 km horizontal resolution at the centre of the study area (~40°N). Bilinear interpolation was used to resample both depth and slope to match the resolution of the climatic variables (0.088°).

Pseudo-absence generation
The most widely used SDM techniques can be broadly categorized into those that utilize only presences, and those that require information on both presences and absences (Barbet-Massin et al. 2012). In cases where real absence data are not available, absence points ('background data' or 'pseudo-absences') can be simulated and extracted from the model domain or 'background'. Despite popular belief, MaxEnt, which is often falsely referred to as a 'presence-only' method, requires information on where a species is considered absent (i.e. pseudoabsences) in order to contrast the conditions at those locations against where it is present (VanDerWal et al. 2009, Barbet-Massin et al. 2012, Guillera-Arroita et al. 2015. Such models are more accurately considered 'presence−background' models (Guillera-Arroita et al. 2015), and should closely resemble the results of a presence−absence model in cases where a species is well-sampled but rare, and background points closely emulate the distribution of true absences (Ward et al. 2009).
As the majority of the V. pourtalesii presence data were extracted from video surveys and discrete sampling locations, real absence data were not available over the majority of the subarea extents and at a scale relevant to the environmental variables modelled in our study. However, null catches from the DFO multispecies trawl survey operating on the Scotian Shelf in the Northeast US/Atlantic Canada subarea could be used to approximate absence locations, an approach which was recently used by Beazley et al. (2018) in a presence−absence Random Forest classification modelling exercise. We evaluated the results of Random Forest and GAM models generated using absences from null catches collected between 2007 and 2017 from the DFO multispecies trawl survey on the Scotian Shelf, augmented with randomly generated pseudo-absences across the remainder of the Northeast US/Atlantic Canada subarea (see Fig. S2). Absences were augmented with pseudo-absences to avoid potential extrapolation of models built on data collected only from the Scotian Shelf to other, environmentally dissimilar locations in the Northeast US/Atlantic Canada subarea.
The Random Forest and GAM models built on absences augmented with pseudo-absences had lower accuracy measures (see Table S2) compared to Beazley et al. (2018) and to the final models presented in this study (see Table 3), performing particularly poorly in their ability to correctly predict the absence data (lower 'specificity'). This was possibly due to the lower number of environmental variables used here, and close proximity of trawl survey absences to the presence data in its core habitat in Emerald Basin and the inability of the models to consistently classify or fit a model to the presences and absences that co-occur in the same water mass. The resulting predicted/projected distribution surfaces (see Fig. S3) appeared to over-extend the suitable habitat of this species into more coastal and shelf depths, where absences from the trawl survey are distributed. Therefore, we chose to run our models using randomly generated pseudo-absences across the entire Northeast US/Atlantic Canada subarea. A description of the methods used for pseudo-absence generation can be found in Text S1.

Species distribution modelling, evaluation, and prediction
Prior to model fitting, correlation between environmental variable pairs was examined using R statistical software version 3.6.1 (R Core Team 2019), and those variables (Depth and Mean Bottom Shear) considered highly correlated (Spearman's rho > 0.7) with others were not considered further. Two different SDM approaches were employed to model the distribution of V. pourtalesii in each subarea: Random Forest classification (Breiman 2001a) and GAMs (Hastie & Tibshirani 1986). Both techniques have been recently used in applications dealing with presences and pseudo-absences (González-Irusta et al. 2015, Morato et al. 2020, and were employed here to maximize both model prediction and inference, and to increase reliability of model predictions in areas outside the environmental domain of the training data. Black box approaches such as Random Forest tend to optimize model prediction (i.e. the ability to predict a response variable based on the set of independent variables; Breiman 2001b) without particular focus on identifying the relationship between the predictors and response, so long as accurate predictions of the response are yielded (James et al. 2013). In contrast, regression methods such as GAMs are better at identifying the type (linear or non-linear), strength, and direction of the relationship between the re sponse and predictors (James et al. 2013), and are therefore more useful in applications requiring biological inference (i.e. understanding how changes in the predicts variables affect the response). Furthermore, when extrapolating outside the domain of the training data, where different physical conditions from those used to train the model may exist, Random Forest predicts the same value as it would for the closest value in the tree for which it had training data (Breiman et al. 1984). Regression techniques are better at predicting trends in the data (Hengl et al. 2018), which is important when extrapolating species−environment relationships into the future.
Random Forest models were fitted using R package 'randomForest' (Liaw & Wiener 2002) with default parameters and 500 trees. GAMs were fitted with a binomial error distribution and cloglog link (see Zuur et al. 2009) using the 'gam' function in R package 'mgcv' (Wood 2011), after a series of generalized linear models indicated non-linear patterns between the response and some predictor variables. The final GAMs in each subarea were selected using a forward stepwise approach, which involved fitting single-variable models, identifying the model with the lowest Akaike's information criterion (AIC) and highest significance, and iteratively adding the other variables and selecting the best combination that resulted in the lowest AIC, and repeating the process until there was no further improvement in AIC. The final GAM built in the Northeast US/Atlantic Canada subarea included 7 variables (bottom and surface temperature and salinity, MLD, bottom current velocity, and slope), while 4 variables (bottom and surface temperature, surface current velocity, and slope) were included in the GAM built in the Mid -Southeast US subarea. Smoothers were initially applied to all terms, and the effective degrees of freedom (edf) were evaluated. Following Wood (2001), smooth terms were replaced with parametric linear terms for variables with edf close to 1. Smoothers were applied to bottom and surface temperature and salinity, and MLD in the Northeast US/Atlantic Canada subarea GAM, and to bottom temperature in the Mid -Southeast US subarea GAM to account for their nonlinear relationship with the response, with knots constrained to a maximum of 4 to avoid overfitting.
Model accuracy metrics and threshold probabilities were derived using 5-fold spatial block crossvalidation (CV), which reduces the effects of spatial autocorrelation on model prediction error , and has been applied in recent deepsea modelling efforts (Guinotte & Davies 2014). Using the 'spatialBlock' and 'spatialAutoRange' functions in the R package 'blockCV' (Valavi et al. 2019), each subarea was partitioned into regularly-shaped 'blocks', and the presence/pseudo-absence data falling within each block were randomly assigned to either the training or testing dataset of each of the 5 folds. To ensure similar prevalence rates between folds, latitudinal and longitudinal offsets were ap plied to shift the blocks on their longitudinal and latitudinal axes in the Northeast US/ Atlantic Canada subarea, and the 'numLimit' parameter was set to 0 to ensure that the most evenly dispersed number of records was chosen (with 500 iterations). The number of presences and pseudo-absences assigned to the training and testing data sets of each fold is shown in Table S3.
For each subarea, Random Forest models and GAMs were fitted to the training datasets from 4 folds and validated on the testing dataset from the fifth fold. The process was repeated 5 times, resulting in 5 CV runs from which accuracy metrics were derived. The threshold-independent area under the receiver operating characteristic curve (AUC) was calculated for each CV run, from which the mean and standard deviation were derived. The 'optimal. thresholds' function in the R package 'PresenceAbsence' (Freeman & Moisen 2008) was used to calculate several common threshold values (e.g. maximum Kappa, sensitivity = specificity, ob served prevalence; Liu et al. 2005) above which a given relative likelihood of occurrence is considered a presence. We selected the threshold which maximizes the sum of sensitivity and specificity (MSS), where sensitivity and specificity represent the pro portion of accurately pre dicted presences and absences, respectively. This threshold was used to convert the likelihood of oc currence outcomes from each CV run into predicted outcomes (0 or 1) that are then summarized into a 2 × 2 confusion matrix. Sensitivity, specificity, and the true skill statistic (TSS; Allouche et al. 2006) were derived and used to assess model performance along with AUC.
From presence−background data, it is not possible to determine whether a species is well surveyed but rare, or common but under-surveyed, a function which requires presence−absence or occupancy detection data (Guillera-Arroita et al. 2015). Consequently, predictions from models that use pseudoabsences over real absence data are more indicative of a species' relative likelihood of occurrence than probability of presence (Guillera-Arroita et al. 2015).
Therefore, relative likelihood of occurrence was adopted herein to refer to the Random Forest and GAM outputs. The final predictions/projections of the relative likelihood of V. pourtalesii occurrence under present-day and future climatic conditions were generated from models trained on the full presence/ pseudo-absence datasets in each subarea: 118 presences and 3373 pseudo-absence in the Northeast US/Atlantic Canada subarea for both Random Forest and GAM (see Text S1) and 18 presences and 1164 pseudo-absences for GAM, and 18 presences and 237 pseudo-absences for Random Forest in the Mid -Southeast US subarea. As all the data are used (important for small datasets), this approach was chosen over averaging the predictions generated from models built in each CV run. Its disadvantage is that the error estimates from CV do not apply perfectly to the final models, as slightly different datasets were used to train them. However, the error rates from CV are likely conservative, with better performance expected from models built on the full dataset .
Model uncertainty was evaluated by calculating the standard deviation of the predictions generated from the 5 CV models. This estimate of standard de viation highlights areas of higher or lower variability between models built on different, spatially partitioned subsets of data, and areas of higher uncertainty may be reflective of local non-stationarity (Nephin et al. 2020). Similar to the accuracy metrics described above, this measure of uncertainty is considered conservative, as standard deviation was derived from the 5 CV models and not the final models used for prediction. To aid in identifying where model outputs may be less certain, especially those of Random Forest, areas of model extrapolation were identified by highlighting each grid cell across both subareas where at least 1 environmental variable had values above or below the present-day variables used to train the final models.
Predictions of relative likelihood of occurrence were thresholded into a binary depiction of suitable vs. unsuitable habitat using the MSS values used to threshold the confusion matrices of each model built in each subarea. Suitable habitat defined here is based solely on the relationship of the V. pourtalesii presences and pseudo-absences with the 8 presentday physical oceanographic and ocean terrain environmental predictor variables used in the models of this study. The relative likelihood of occurrence projected to the 2 future climatic scenarios was also thresholded using MSS, and areas that experienced a gain, loss, or no change in suitable habitat from present-day con ditions were evaluated.

Variable importance and functional response curves
The importance of the predictor variables in Random Forest models was assessed using the mean decrease in Gini index extracted using the 'importance' function in the package 'randomForest'. For GAMs, the 'varImp' function in the package 'caret' (Kuhn 2020) was used to extract the importance of the predictor variables based on their p-value. Following Lopes et al. (2019), 'functional response curves' were generated to evaluate changes in the relative likelihood of occurrence across environmental variable gradients. This method allowed for a direct comparison of the relative likelihood of occurrence− environment relationships for the present-day, RCP 4.5, and RCP 8.5 scenarios. Relative likelihood of occurrence and the value of environmental predictors for each raster cell across each study area were extracted and smoothed using the 'loess' method with span = 0.5 applied to each curve. Curves were displayed with 95% standard error confidence intervals.

Changes in environmental conditions between present-day and future climatologies
The overall change in environmental conditions between present-day and the 2 future climatic scenarios was more prominent in the northern portion of Vazella pourtalesii's range than across its southern distribution (see Table S4, Figs. S4 & S5). Bottom and surface temperature showed the largest change, with average Mean Bottom Temperature increasing by 0.80 and 0.38°C between present-day and RCP 8.5 in the Northeast US/Atlantic Canada and Mid -Southeast US subareas, respectively, and Mean Surface Temperature by 0.90 and 0.58°C (Table  S4). This was followed by maximum MLD, which became shallower across both subareas in the future (change in average MLD: −0.59 and −0.35 m in the Northeast US/Atlantic Canada and Mid -Southeast US subareas, respectively). Changes in bottom salinity, bottom and surface current, and bottom shear between the present-day and future forecasted conditions were negligible in both sub areas, while sur-face salinity decreased in both sub areas but more drastically in the northern subarea (change in Mean Surface Salinity −0.25 and −0.09 in the Northeast US/Atlantic Canada and Mid -Southeast US subareas, respectively; Table S4).

Model performance
Accuracy metrics (AUC, sensitivity, specificity, and TSS) of the Random Forest models and GAMs built in each subarea are shown in Table 3, while additional summaries of the parametric and smoothed terms of the GAMs are shown in Table S5. In the Northeast US/Atlantic Canada subarea, both Random Forest and GAM had an excellent performance rating based on AUC (> 0.90), with GAM marginally outperforming Random Forest based on all 4 metrics. Models generated in the Mid -Southeast US subarea had a comparatively poor performance, with AUC values in the low 0.80 range, and low TSS values (0.52 and 0.49 for Random Forest and GAM, respectively). Models in both subareas performed similarly in terms of their ability to accurately predict the presences (i.e. sensitivity) and pseudo-absences (i.e. specificity) correctly, with specificity lower than sensitivity for all models ( Table 3). The environmental variables considered most important for the present-day distribution of V. pourtalesii (see Fig. S6) were Mean Bottom Temperature and Mean Surface Salinity for Random Forest and GAM, respectively in the Northeast US/ Atlantic Canada sub area. Mean Surface Temperature and Slope were the most im portant variables for these models in the Mid -Southeast US subarea (Fig. S6 Table 3. Accuracy measures for the Random Forest models and generalized additive models (GAMs) trained and tested on the Vazella pourtalesii presence/ pseudo-absence data in both the Northeast US/Atlantic Canada and Mid -Southeast US subareas. Cross-validation was done via 5-fold spatial blocking with random assignment of blocks into folds. Sensitivity, specificity, and the true skill statistic (TSS) were generated from a confusion matrix of tabulated outcomes that was thresholded using the maximum of sensitivity + specificity (MSS) identified for each model and subarea. AUC: area under the receiver operating characteristic curve

Present-day and future distribution -Northeast US/Atlantic Canada
The areas predicted with a higher relative likelihood of occurrence of V. pourtalesii under present-day conditions were relatively consistent between Random Forest and GAM in the Northeast US/ Atlantic Canada subarea (Fig. 2). In accordance with the spatial distribution of the occurrence data ( Fig. 1), the Scotian Gulf, an inlet which opens up into the LaHave and Emerald Basins on the inner shelf, and the Northeast Channel, a deep channel at the mouth of the Gulf of Maine (see Beazley et al. 2018 for further details), were predicted with high relative likelihood of occurrence of V. pourta lesii by both models. However, under future environmental conditions, the spatial distribution of the projections diverged between models. Random Forest projected a higher relative likelihood of occurrence of V. pourtalesii in the Gulf of Maine than GAM, while GAM projected a higher relative likelihood of occurrence in the Laurentian Channel/ Gulf of St. Law rence than Random Forest. The spread into the Laurentian Channel/Gulf of St. Lawrence intensified under the RCP 8.5 scenario for both models, although GAM projected a higher relative likelihood of occurrence further north, into the mouth of the St. Lawrence estuary.
The binary depiction of suitable versus unsuitable habitat and comparison of the percent change in areas projected as suitable habitat between the different climatic scenarios (Figs. 3 & 4) also exemplify the overall range expansion of V. pourtalesii under future climate change. However, certain areas of suitable habitat in the Laurentian Channel and Gulf of St. Lawrence projected by GAM were also asso- ciated with high model uncertainty (Fig. 5). Evaluation of the prediction/ projection outputs from each of the 5 CV models (not shown) showed that projections of relative likelihood of occurrence from the model based on the data partitioned into one particular fold (Fold 3) were erroneously high in the lower Laurentian Channel, mouth of the Gulf of St. Lawrence estuary, and northeast Scotian Gulf, congruent with the areas of higher uncertainty in Fig. 5. The training data of this fold had the lowest number of presences (70 compared to 89−114) and the highest number of pseudoabsences (2804 compared to 2645− 2716) of all folds (see Table S3), and in contrast to the other folds, the presence data selected by spatial blocking for Fold 3 were located in the lower Scotian Gulf and Northeast Channel, south of the main sponge grounds in Emerald Basin. While the CV models for Random Forest were built on the same subsets of data as the GAMs, their model projections were not as greatly affected by the varying numbers of presences and pseudoabsences in each fold, nor by their spatial distribution. This was reflected by the relatively low standard deviation across model predictions (Fig. 5). Random Forest is more prone to poor model results when predicting/projecting outside the environmental envelope of the training data. However, those areas where model uncertainty was high for both Random Forest and GAM were not considered extrapolated (see Fig. S7) with the exception of the area where GAM projected suitable habitat in the mouth of the Gulf of St. Lawrence estuary under the RCP 4.5 and 8.5 scenarios.
Another notable difference between model projections in this subarea was the decrease in rela-11 Fig. 3. Vazella pourtalesii suitable habitat (red) predicted/projected from Random Forest models (top row) and GAMs (bottom row) under present day, and representative concentration pathway (RCP) 4.5 and RCP 8.5 future climatic conditions. Relative likelihood of occurrence surfaces from Fig. 2 were thresholded using the maximum of sensitivity + specificity in Table 3. 40°N is indicated by the grey dashed line tive likelihood of occurrence of V. pourtalesii in its core habitat in the Scotian Gulf and Northeast Channel projected by Random Forest for both RCP scenarios (Fig. 2). In contrast, GAM projections of relative likelihood of occurrence intensified in the future. Although these areas were generally associated with higher model uncertainty (Fig. 5), this high uncertainty was due to the projections resulting from the GAM fit to the training data in Fold 3. With Fold 3 excluded (not shown), the standard deviation across the GAM predicted/ projected surfaces was relatively low, and highest in the Gulf of St. Lawrence estuary. The slight northward shift and contraction in the areas of highest relative likelihood of occurrence in the Scotian Gulf and Northeast Channel, respectively, as projected by GAM under RCP 8.5 suggests that these model results collectively indicate that some impacts to the future distribution of V. pourtalesii may be incurred with further climate change. Fig. 6 shows where Random Forest and GAM projected a gain, loss, and no change (i.e. the core suitable habitat) in suitable habitat between present-day and future forecasted conditions. For both models, the core of this species' habitat in the Scotian Gulf and Northeast Channel in the Northeast US/Atlantic Canada subarea remained relatively stable between the present-day and future forecasted conditions, suggesting these areas may serve as refugia against climate change in the future. Areas projected to gain suitable habitat under the RCP 4.5 and 8.5 scenarios were, on average, deeper than the present-day habitat of V. pourtalesii (see Table S6), indicating an extension of this species into deeper waters as a result of climate change. For Random Forest, these areas were associated with warmer mean bottom temperatures (7.39 ± 1.63°C SD, up to 12.48°C; RCP 8.5) compared to the present-day habitat in the Scotian Gulf and Northeast Channel/Gulf of Maine (Fig. 6), where temperatures were on average 6.97 ± 0.78°C and reached a maximum of 9.20°C. In contrast, areas where GAM projected a gain in suitable habitat were, on average, colder than the areas projected by Random Forest (6.04 ± 0.72°C for RCP 8.5, maximum: 9.33°C), and also slightly colder than average bottom temperature in its present-day habitat (6.77 ± 0.92°C in present day; although note the aforementioned higher uncertainty in GAM projections into the Gulf of St. Lawrence estuary). This difference is likely due to Random Forest projecting a larger gain in habitat in the Gulf of Maine, where simulated bottom temperatures are warmer than in the Laurentian Channel where GAM projected a larger gain than Random Forest (Fig. S4, Table S7). Table S6 shows that while the average depth of present-day suitable habitat predicted by both Random Forest and GAM is similar (203 ± 125 and 219 ± 177 m, respectively), the average depth of suitable habitat projected by GAM under RCP 8.5 is over 50 m deeper (335 ± 164 m) than that projected by Random Forest (285 ± 178 m), likely explaining why average bottom temperatures are slightly cooler in its future suitable habitat.
The slight peak in relative likelihood of occurrence at 6°C and larger mode peaking at ~9°C in the Random Forest functional curve for Mean Bottom Temperature (Fig. 7a) also reflects this pattern, with the smaller mode corresponding to temperatures in the Laurentian Channel/Gulf of St. Lawrence and the larger mode corresponding to the Scotian Gulf and Northeast Channel/Gulf of Maine where this model projected a larger gain in suitable habitat. For GAM, relative likelihood of occurrence under RCP 4.5 and 8.5 peaked at cooler temperatures (~6°C; Fig. 7a), corresponding to the Laurentian Channel and Gulf of St. Lawrence where it projected a larger gain than Random Forest. Bottom temperatures in both the Gulf of Maine and Laurentian Channel/Gulf of St. Lawrence warmed from their presentday con ditions (Fig. S4, Table S7), suggesting that the habitat of V. pourtalesii may respond positively to warming bottom waters. Patterns in the functional response curves presented here were similar to the presence−absence Random Forest partial dependence plots of Beazley et al. (2018) for similar environmental variables. For instance, both Random Forest and GAM models showed an increase in relative likelihood of occurrence at ~6°C and ~34 along the gradients in Mean Bottom Temperature and Salinity, respectively, similar to the peaks in Minimum Bottom Temperature (5−6°C) and Salinity (33.5) as shown by Beazley et al. (2018), suggesting that our use of pseudo-absences closely emulates the predicted distribution of V. pourtalesii based on real absence data.

Present-day and future distribution -Mid -Southeast US subarea
Patterns in the present-day and future predicted/projected distribution of V. pourtalesii in the Mid -Southeast US subarea were not as clearly de fined as in the northern portion of its range. Under presentday environmental conditions, both Random Forest and GAM predicted low to moderate V. pourtalesii relative likelihood of occurrence in a band following the Florida peninsula (Fig. 2). This band is consistent with the location of V. pourtalesii occurrences around the peninsula and on Blake Plateau (Fig. 1). While uncertainty associated with both models was relatively low (Fig. 5), predictions/ projections were low at the location of presence localities from both models (5% relative like- Forest models (top row) and GAMs (bottom row) generated in 5-fold spatial block cross-validation lihood of occurrence by GAM and up to 50% by Random Forest), reiterating their poor performance in this subarea (Table 3; Table S5). Changes in the intensity and distribution of the relative likelihood of occurrence under the 2 future emission scenarios were relatively small for Random Forest and virtually un changing for GAM (Figs. 2 & 3). Random Forest projected a slight gain in suitable habitat in the deeper waters of the Blake Plateau under RCP 4.5 conditions (Fig. 3). Similar to the pattern observed in the northern portion of its range, these areas corresponded to slightly colder temperatures than the present-day habitat of this species (see Table S8). However, the total area projected as suitable habitat by Random Forest showed only a negligible increase from RCP 4.5 to 8.5 (Fig. 4), suggesting that the warmer portion of this species' range will not benefit from additional ocean warming. GAM also projected a slight increase in suitable habitat under RCP 4.5, followed by a decrease back to its presentday size under RCP 8.5 (Fig. 4). Fig. 6 shows a loss of suitable habitat in the deeper waters of the Blake Plateau as well as along the flanks of the refugia areas, which was slightly more pronounced for RCP 8.5 than RCP 4.5.
Relative likelihood of occurrence predicted/ projected along the present-day and future gradients in Mean Bottom Temperature were highest at ~10− 11°C for both Random Forest and GAM in the Mid -Southeast US subarea (Fig. 7b). However, Mean Bottom Temperature was not the top environmental predictor variable in either model (Fig. S6), suggesting that any inferences based on temperature should be made with caution. The selection of Mean Surface Temperature by Random Forest as its top predictor is likely due to the spatial congruence between the location of V. pourtalesii occurrences (Fig. 1)  salinity is unitless of the Gulf Stream whose influence extends to the seafloor both off Florida in water depths of 700 m and further north on the Blake Plateau where depths reach 1000 m (Richardson 2019). While temperatures associated with this current are projected to increase in the future (Fig. S5), the associated projected occurrences remained relatively stable (Fig. 2), suggesting that other variables or interactions between variables may have a greater influence in this subarea. The functional response curves for both models (Fig. 7b) showed little change in the relative likelihood of occurrence predicted/projected for the present-day and future forecasted conditions, especially for GAM. Of the 4 environmental predictors included in GAM (Mean Bottom Temperature, Mean Surface Current, Slope, and Mean Surface Temperature), the static variable Slope was considered the most important ( Fig. S6), possibly explaining why spatial patterns in predicted/projected relative likelihood of occurrence and suitable habitat (Figs. 2 & 3) remained relatively unchanged between present-day and future forecasted conditions. Overall, the performance of both Random Forest and GAM models in this subarea was relatively poor (Table 3; Table S5), with poor congruence be tween the spatial predictions of relative likelihood of occurrence and the location of presence points under present-day environmental conditions (especially for GAM; see Figs. 1 & 2), suggesting that any inferences derived from these models should be taken with caution.

DISCUSSION
The results of our study show that under future climatic conditions, the potential suitable habitat of the glass sponge Vazella pourtalesii will in crease in the northwest Atlantic. Using Random Forest and GAM techniques, we projected a gain in habitat of up to~4 times its present-day size using environmental conditions simulated under moderate (RCP 4.5) and worst-case (RCP 8.5) CO 2 emission scenarios for 2046-2065. In the northern portion of this species' range, we projected a shift into deeper waters and higher latitudes, a pattern consistent with previous studies of the effects of climate change on marine species in the region (Nye et al. 2009, Morato et al. 2020). This shift in distribution was likely due to the warming of the colder waters surrounding its current habitat, and the availability of new, unoccupied niche for this deep-water, originally subtropical species (Schmidt 1870).

Effects of climate change on sponges: a deep-water example
Fitness and adaptability of individuals within the same population often differ between individuals located at the core and leading/trailing edges of its range, resulting in observed contractions and expansions at the trailing (southern) and leading (northern/ colder) edges of its distribution, respectively, in response to climate change (Rilov et al. 2019). Our projected gain in V. pourtalesii suitable habitat in the future, particularly in the northern portion of its range where a colder climate prevails, is not entirely unexpected. Its current distribution on the Scotian Shelf is highly associated with areas that experience regular ingression of Warm Slope Water (Beazley et al. 2018), a warmer, saltier water mass originating from the Gulf Stream. Under a CO 2 doubling scenario, Saba et al. (2016) predicted that the enhanced warming on the northwest Atlantic shelf was the result of increased incursion of Warm Slope Water into the region due to a northerly shift in the Gulf Stream and the retreat of the cold and fresh Labrador Current. The Laurentian Channel/Gulf of St. Lawrence, where significant gains in suitable habitat are projected, is projected to warm at least 1°C (see Fig. S4 and Table S7), shifting the climate closer to the temperature preference of V. pourtalesii as defined by its present-day distribution (~6°C; Fig. 7a). The southern range of V. pourtalesii did not contract under further warming as expected, and was instead projected (by Random Forest) to nearly double in size under RCP 4.5 conditions. However, no further expansion was projected by this model under RCP 8.5. The area predicted/projected as suitable habitat by GAM was similar between the present-day and future climatologies, but showed a slight decrease in size from RCP 4.5 to 8.5, suggesting that this species does not benefit from further warming in the warmer part of its range.
In contrast to the Northeast US/Atlantic Canada subarea, the mechanism for the projected shift of V. pourtalesii into deeper waters in the southern portion of its range may be due to the warming of waters beyond its upper thermal tolerance limit, causing it to retreat into deeper, cooler waters. Although both Random Forest and GAM models showed a loss of suitable habitat along the flanks of its core distribution, a clear extension of the leading (cooler) edge and decline of the trailing (warmer) edge of its core habitat was not apparent in this subarea (Fig. 6). Both models predicted poorly in that area, possibly due to the low number of presence points (18 presences). Wisz et al. (2008) evaluated the predictions of 12 different modelling algorithms applied to 46 species occurring at 3 different frequencies (10, 30, and 100 presence records), and found that no model predicted consistently well when sample size was less than 30. This suggests that model predictions/projections in this subarea should be considered for exploratory purposes only until additional information on the distribution of V. pourtalesii is collected.
Species affiliated with narrow environmental envelopes are often easier to model than those distributed over wide environmental gradients or mosaics (McPherson & Jetz 2007, Tsoar et al. 2007). High environmental variability associated with only a few V. pourtalesii presences (18) and/or local heterogeneity that is difficult to capture in SDMs could possibly explain the poor model results in the Mid -Southeast US subarea. For instance, the densest aggregations of V. pourtalesii off Florida co-occurred with vast regions of dead Lophelia pertusa reef. Niche overlap between L. pertusa reefs and sponges has been previously noted (van Oevelen et al. 2018), and the use of L. pertusa as settlement substrate by V. pourtalesii may be to take advantage of the current-induced flow and enhanced food supply due to interactions between reef topography and local hydrodynamic conditions (Davies et al. 2009, Mohn et al. 2014). Such small-scale associations would not be accurately captured by our coarser-scale data (Austin & Van Niel 2011).
Studies of shallow-water sponges have shown that while many appear to be relatively resilient to the combined effects of ocean warming and acidification compared to other benthic taxa, their tolerance to various climate stressors is highly species-specific (Bell et al. 2018). Bell et al. (2018) conducted a review of studies on the singular and combined effects of ocean warming and acidification on sponges and found that of the 44 studies focused solely on ocean warming, 30 demonstrated negative effects to host physiology, including pumping rate, filtration efficiency, choanocyte chamber size, and density; to gene expression, feeding ecology, and reproductive output; and to microbial community composition and function. Ramsby et al. (2018) examined the effects of incrementally increasing ambient seawater temperature on the bioeroding sponge Cliona orientalis and its symbiotic dinoflagellates (Symbiodinium), and found little effect of raising temperatures 2°C above monthly mean values. However, at 3°C above average (consistent with temperatures predicted to occur in year 2100 under RCP 8.5), C. orientalis bleached and reduced its energy reserves, consistent with the response of sympatric corals. Bennett et al. (2017) demonstrated via ex situ experimen-tation a tolerance of several Great Barrier Reef sponges to ocean temperatures forecasted under moderate (RCP 6.0; 3°C above ambient), but not high (RCP 8.5; 4.5°C above ambient) CO 2 emission scenarios.
We modelled changes in the distribution of the deep-water sponge V. pourtalesii under moderately strong climate policy (RCP 4.5) and worst-case (RCP 8.5) CO 2 emission scenarios, where forecasted anomalies in Mean Bottom Temperature in the Northeast US/Atlantic Canada subarea in 2046-2065 were, on average, + 0.57 ± 0.37°C (SD) and + 0.80 ± 0.56°C under RCP 4.5 and 8.5, respectively. In the core of this species' habitat in the Scotian Gulf and Northeast Channel, such anomalies would raise average bottom temperatures by a maximum of ~2°C (tõ 10°C under RCP 8.5) from their present-day average maximum (~8°C), which is below the thermal change often shown to cause deleterious effects in shallow-water sponges (~3°C above ambient). Emerald Basin, where V. pourtalesii forms the densest known aggregations of its kind, is subjected to strong inter-annual and multi-decadal variability in water mass characteristics, where empirical maximum bottom temperatures of 12°C were recorded from CTD and Argo float data collected between 1950 and 2015 (Beazley et al. 2018). Exposure to short-term fluctuations in environmental conditions may have important consequences for species' ranges (Wethey et al. 2011), and have shown to cause mass mortality events in deep-water sponge grounds elsewhere (Guihen et al. 2012). Our use of average conditions over monthly or seasonal data likely do not capture the overall temperature range ex perienced by these sponge grounds historically, nor the extremes that will likely be incurred in the future. Nonetheless, our correlative SDMs using mean climatic conditions show that while the population may persist into the future, further warming of V. pourtalesii's core habitat in Emerald Basin may result in a reduction in its relative likelihood of occurrence there compared to the present day. While it has long been assumed that the deep sea is a highly stable environment relatively buffered from the effects of short-term changes in atmospheric or surface conditions, recent studies suggest that deep-sea ecosystems may be relatively sensitive to even seasonal shifts in upper ocean conditions (Glover et al. 2010). As the bottom temperatures forecasted under RCP 4.5 and 8.5 in the core habitat of V. pourtalesii are not outside the range experienced historically by these sponge grounds, or warmer than the conditions associated with V. pourtalesii presence in the southern portion of its range, the reduction in relative likelihood of occurrence suggests that even marginal in creases in temperature in the deep sea may have deleterious impacts on animal physiology and habitat quality, and/or that variables other than bottom temperature, or possibly synergies between variables, may be responsible for this decline.
The occurrence of V. pourtalesii in certain localities off the mid -southeast USA, where empirical maximum bottom temperatures of 11.2°C were recorded (Table S1), points to an ability of this species to acclimate to temperatures higher than the average conditions experienced on the Scotian Shelf, suggesting some hope for the persistence of the dense sponge grounds that reside there. However, the densest aggregations of V. pourtalesii in the Mid -Southeast US subarea (up to 663 individuals on a single transect) were typically associated with bottom temperatures no higher than 8.6°C, suggesting that individuals found at ~12°C may be more representative of marginal than central populations, on the brink of their thermal tolerance limit (Bennett et al. 2019). The formation of sponge grounds at these higher temperatures is likely less probable. Evaluation of the Random Forest and GAM response curves for the Mid -Southeast US subarea (Fig. 7b) also support a preference for bottom temperatures <12°C (maximum relative likelihood of occurrence at ~10°C). Interpretation of our model results would greatly benefit from laboratory-based experiments where the physiological effects of incrementally increasing temperatures beyond 10°C are measured, and the upper thermal tolerance of V. pourtalesii definitively identified.

Limitations and importance of environmental interpretation
While we attempted to ensure a similar prevalence (i.e. proportion of observed presences) between data partitions (folds) used in 5-fold spatial-block crossvalidation, the lower number of presences and higher number of pseudo-absences, along with their spatial distribution, had an impact on the model outputs for GAM (but not for Random Forest). For irregularly sampled data, spatially clustered data, or those with a highly unbalanced prevalence such as ours, Roberts et al. (2017) suggested non-gridded and/or irregularly-shaped blocks (i.e. pie slices) to ensure even sampling of presence and absence data, but also cautioned that non-regular blocks may not address the issue of autocorrelation consistently. Caution must also be taken when using spatial blocking, as blocking structures that follow environmental gradients can lead to entire portions of envi-ronmental predictor space (i.e. ranges and/or combinations of predictor variables) being held out in the testing dataset, resulting in extrapolation between folds . The CV GAM model trained on the data partitioned in Fold 3 led to erroneously high projections of relative likelihood of occurrence in the lower Laurentian Channel, mouth of the Gulf of St. Lawrence estuary, and northeast Scotian Gulf, consistent with the locations of high model uncertainty. Similar applications of Random Forest modelling in the Gulf of St. Lawrence using groundfish trawl survey presences and absences (see Murillo et al. 2016) indicated an absence of V. pourtalesii in those areas. While we did not evaluate whether the CV model was extrapolating beyond the data used to train the model, unlike the data partitioned in the other folds, the dataset of Fold 3 excluded presences from the core sponge grounds in Emerald Basin. We recommend that alternate approaches to rectangular spatial blocking be investigated in the future for highly clustered species datasets such as those modelled here, to ensure uniform data selection between folds across the environmental domain of the species.
While our model results are an important first step in understanding the impacts of climate change on deep-water sponges, the caveats and uncertainties inherent to this particular type of habitat suitability modelling should be considered prior to the development of strategies and policies for the management of this species in the future. Climate trajectories inherently include a broad range of assumptions on future greenhouse gas emissions, population and economic growth, and technological change (van Vuuren et al. 2011), the uncertainties of which are carried forward into ecological models (Payne et al. 2016). RCP 8.5 trajectories are associated with the highest level of uncertainty, and changes in projected distribution under this RCP should be considered the worst-case scenario. Many factors can influence the response of a species to changes in the environment, including recruitment dynamics, physiological tolerances, food availability, competition, and community composition of existing or receiving communities (Poloczanska et al. 2016), which are difficult to capture in single-species models such as those presented here.
Correlative SDMs such as those applied here are based on the statistical association between occurrence records and spatial environmental data, and inherently capture complex interactions between a species and its current environment. When projecting future distributions, these interactions are assumed to be preserved in the new environmental space, an assumption that may not be valid under future climate change where novel environmental scenarios and non-equilibrium species distributions are likely (Kearney et al. 2010). Furthermore, we defined the habitat of V. pourtalesii in relation to only a few dynamic environmental variables (bottom and surface temperature, salinity, current velocity, and MLD). Other abiotic factors not captured here may be equally as important for the distribution of this species, such as nutrient availability and primary production (Beazley et al. 2018, Kazanidis et al. 2019, both of which are expected to decrease on the Scotian Shelf in the future (Pepin et al. 2013, Lavoie et al. 2019. Further interpretation of the results of these models with other information important for species' distributions is necessary in order to evaluate their uncertainty. For instance, while both Random Forest and GAM were consistent in projecting gains in suitable habitat in the deeper waters of the Laurentian Channel and Gulf of St. Lawrence, the direction of flow of oceanographic currents in the future, which was not captured by our models, does not support the transport of larvae beyond the edge of the Scotian Shelf. Although the east-to-west flow of the Labrador Current along the edge of the Scotian Shelf is predicted to weaken (Saba et al. 2016), no significant transport from the Scotian Gulf to the Laurentian Channel is anticipated under future climate change (see Fig. 6 of Saba et al. 2016). In contrast, the proportion of warm and salty Slope Water entering the Gulf of Maine through the Northeast Channel is predicted to increase (Saba et al. 2016), suggesting that the gain in suitable habitat in the Gulf of Maine projected by Random Forest (Fig. 6) is a much more likely event. Furthermore, the presence of hard substrate is crucial for the settlement of V. pourtalesii, which attaches to substrate ranging in size from pebbles to boulders (Hawkes et al. 2019). The Laurentian Channel is dominated by mud and clay substrate (DFO 2011), further reiterating that the settlement and proliferation of this species there is unlikely. This highlights the need for further interpretation with other factors not included in these models that may be important for the spatial distribution of a species.

Deep-water sponges as 'climate change winners'?
The predicted impacts of climate change on the home ranges of other marine species residing on the continental shelves of the northeast USA and Atlantic Canada are, to date, highly species-specific, with emergence of clear 'winners' and 'losers' , McHenry et al. 2019, Morato et al. 2020. For the distribution of 125 benthic invertebrates, and pelagic and demersal fish on the northeast USA shelf, McHenry et al. (2019) predicted more losers than winners, with those species highly associated with the seabed (benthic invertebrates and demersal fish) experiencing the most severe habitat loss and/ or fragmentation. In contrast, highly mobile species such as pelagic fish, cephalopods, and elasmobranchs showed a greater ability to shift their home ranges north to help mitigate climate-induced effects. For long-lived, habitat-forming cold-water corals in the North Atlantic (including our study area), Morato et al. (2020) predicted severe habitat losses ranging between 28 and 100%. These results suggest that the impacts on benthic or sedentary species may be more severe than on highly mobile species, varying by taxa likely as a consequence of longevity, reproductive periodicity, and dispersal capacity. Nonetheless, the projections of McHenry et al. (2019) and Morato et al. (2020) were made using a different suite of environmental variables and for a later time period than modelled here, when the severity of climate-induced effects on marine habitats will likely be higher.
If based solely on the projected expansion in suitable habitat in the future, our results may lend further weight to the concept that sponges, even in deep-sea environments, are potential 'winners' against climate change compared to other benthic groups. Whether the distribution of boreal or arctic sponge grounds in the North Atlantic will benefit from ocean warming in the same way remains unknown. Boreal assemblages may have some ability to shift their distribution north to help mitigate the effects of ocean warming. However, a laboratory-based study of the ef fects of increasing temperatures on the transcriptomic response of the extrem o phile Antarctic demosponge Isodictya sp. (González-Aravena et al. 2019) suggests that cold-adapted sponges may have an even more limited ability to tolerate increased temperatures compared to warm-adapted species. This suggests that Arctic assemblages located at their maximum latitudinal limit may therefore experience an overall loss of habitat due to contraction of their southern range as temperatures warm. Nonetheless, the local formation of sponge grounds has been strongly linked to factors other than temperature, such as the occurrence of internal tidal waves at the seabed (Klitgaard et al. 1997, Klitgaard & Tendal 2004 and the presence of water mass fronts (Klitgaard & Tendal 2004, Roberts et al. 2018, Meyer et al. 2019, which should be accounted for when making inferences of climate-induced distribution changes. Such characteristics are often poorly resolved in climate forecasting models.

Implications for conservation management
The size of the V. pourtalesii core habitat in the Scotian Gulf and Northeast Channel remained relatively unchanged between our present-day predictions and future projections (Fig. 6), suggesting that these areas may serve as important refugia into the future. However, the projected relative likelihood of occurrence of V. pourtalesii under both RCP scenarios was lower in these areas compared to presentday conditions (see Random Forest outputs; Fig. 2), including at the location of the 2 bottom-fishery closures implemented by DFO in 2013 for the protection of this species from bottom-fishing activities. While relative likelihood of occurrence or presence probability is not a direct correlate of abundance or biomass, it could be considered a proxy for habitat quality, where a higher-quality en vironment would presumably support higher species' abundances. In 2013, DFO closed 2 areas equating to nearly 260 km 2 in Emerald Basin to all bottom-tending fishing gears to protect 2 of the 5 densest concentrations of V. pourtalesii in the Scotian Gulf. As the purpose of the 2 sponge conservation areas is to protect the most significant concentrations (based on biomass/density) of this species, range shifts and potential abundance reductions under future climate change present a significant challenge for maintaining the effectiveness of these fishery closures. Effective conservation management of V. pourtalesii will require an iterative monitoring programme designed to track changes in its density and distribution within and outside these closed areas, with emphasis placed on the margins between core areas and those where a loss or gain of habitat are projected.

CONCLUSIONS
The results of our correlative SDMs show that the suitable habitat of the deep-water glass sponge V. pourtalesii will expand in the northwest Atlantic under future climate change, particularly in the northern portion of its range (Atlantic Canada), where ocean warming will serve to improve the conditions surrounding the current habitat of this species and increase the availability of currently unoccupied niche. While not all of these areas are likely to be populated in the future, any expansion of V. pourtalesii's niche will likely have positive implications for other benthic species, which were shown to be more abundant and diverse in the presence of this habitat-forming species (Hawkes et al. 2019), a function they have in common with other sponge grounds across the North Atlantic (Klitgaard 1995, Bo et al. 2012, Beazley et al. 2013. Although our projected expansion of suitable V. pourtalesii habitat suggests that this species may be a 'winner' against future climate change, the reduction in relative likelihood of occurrence in its core habitat on the Scotian Shelf suggests that the Vazella sponge grounds that re side there may experience losses as the northwest Atlantic continues to warm. We recommend that future studies aimed at elucidating the effects of climate change on the distribution of V. pourtalesii combine mechanistic approaches based on functional traits and physiological constraints with correlative SDM techniques (see Kearney et al. 2010) such as those presented here, in order to strengthen their results. Nonetheless, our results are an important first step in evaluating the impacts of climate change on a ground-forming deep-water sponge, serve as a basis for hypothesis testing in future laboratory-based physiological studies, and provide a starting point for its effective conservation management in light of climate change.
While the impacts of climate change on the distribution, abundance, and composition of other sponge grounds in the North Atlantic remain unknown, given their importance to various ecological functions, any positive or negative impacts of anthropogenically induced change are likely to have cascading effects throughout these benthic ecosystems.