Predicted shifts in suitable habitat of interacting benthic species in a warmer and invaded Canadian Arctic

Species interactions are an integral component of every ecosystem. Together with environmental variables, these interactions help to shape community structure with respect to species composition and abundance. However, interactions between species may also be modulated by changes in the environment (i.e., Van der Putten et al., 2010; Durant et al., 2019) and similar species (e.g., of the same trophic level or functional group) may react differently to environmental change (Durant et al., 2007; Schweiger et al., 2008). These changes not only may affect species and modulate species interactions but also may impact higher trophic levels, with potential consequences to the wider ecosystem (Winder and Schindler, 2004; Durant et al., 2007).

Benthic habitats and organisms are essential components of aquatic ecosystems. In addition to often creating habitat, as in the case of kelp forests and eelgrass meadows that provide shelter for other species due to their 3-dimensional structures, benthic organisms may also play an important role in the food web and energy transfer to higher trophic levels (Dunton and Schell, 1986; Debenham, 2005; Endrédi et al., 2021). Indeed, benthic organisms play a fundamental role in the transfer of energy from the benthic environment to the water column and in energy cycling (McMeans et al., 2013; Endrédi et al., 2021).

This transfer and cycling of energy could be of particular importance in the Arctic, where benthic environments and organisms are thought to use a wider range of carbon sources and are characterized by energy cycling processes that are influenced by the high seasonality of the ecosystem (McMeans et al., 2013). This flow and functionality may be greatly impacted if environmental conditions change, as is already happening in the Arctic due to global change. This region is warming faster than the global average, with ice-free summers predicted in less than 3 decades (Ballinger et al., 2020; Hwang et al., 2020; Previdi et al., 2021). Analyses may need to be projected over much longer time scales in the Arctic because surface air temperature is predicted to warm much more than in other areas of the globe, to greater than 4°C by 2100 and 8°C by 2500 (Lyon et al., 2022), and winter temperatures are expected to increase by around 13°C by 2100 under the most extreme global change scenario (Overland et al., 2019). In particular, the Canadian Arctic is experiencing greater impacts than other Arctic regions, where more rapid increases in temperature and sea ice loss (duration and concentration) have been recorded (Mudryk et al., 2018; Derksen et al., 2019; Flato et al., 2019).

Environmental change affects native species directly and may also impact the rate of arrival of new species (e.g., new shipping routes), thereby enhancing their impact (Bennett et al., 2021). Invasive species are known to be one of the main causes of biodiversity loss and to play an important role in the extinction of native species (Bellard et al., 2016; Blackburn et al., 2019). Invasive species are also known to engage in strong competitive interactions with native species, competing for habitat and resources, and affecting the density and biomass of species in the receiving environment (e.g., Forrest and Taylor, 2002; Wallentinus and Nyberg, 2007; Kotta et al., 2010; Sorte et al., 2010; White et al., 2010; Drouin et al., 2012). Moreover, the impact caused by invasive species has been shown to be more detrimental relative to the effect of warmer temperatures in the ecosystem (Lopez et al., 2022). The rate of introduction of invasive species is lower in the Arctic than in temperate regions (Casas-Monroy et al., 2014; Chan et al., 2019). However, they have been identified as an emerging issue for the Arctic (Ricciardi et al., 2017), where new introductions are being reported or detected with increased survey effort together with the use of technologies such as environmental DNA and metabarcoding (MacDonald et al., 2010; Mathieson et al., 2010; Goldsmit et al., 2014; Brown et al., 2016; Chain et al., 2016; Golder, 2018; Grey et al., 2018; Lacoursière-Roussel et al., 2018; Dispas, 2019; Dhifallah et al., 2022). The observed trend of growing shipping traffic in the region over the past few years (Dawson et al., 2018; Protection of the Arctic Marine Environment, 2020) is predicted to continue due to ice retreat and the opening of new shipping routes (Mudryk et al., 2021) creating new possibilities for further introductions.

The sum of these changes will inevitably modify energy flow in Arctic marine systems with impacts to the food web (Winder and Schindler, 2004; Durant et al., 2019; Florko et al., 2021). Regime shifts in Arctic benthic ecosystems due to increased sea surface temperature (and consequent decreases in ice cover) can reorganize benthic communities, including increasing macroalgae cover (Kortsch et al., 2012). Such changes generally happen faster and more frequently in marine environments than in terrestrial ones since the dispersal of offspring is easier in water (Pinsky et al., 2013; Pinsky et al., 2020; and references therein).

Such a changing Arctic marine ecosystem may impact the spatial distribution of species. Tools such as species distribution models (SDMs) can provide predictions related to the potential changes in habitat suitability for different species. These types of models are correlative approaches that predominantly consider abiotic factors in their predictions and can provide estimates of ecological niches or availability of suitable habitat based on expected environmental conditions (Phillips et al., 2006). However, one of the main drawbacks of this methodology is that it does not account for biotic interactions, which precludes use of SDMs to evaluate the role that these interactions may play even at macroecological scales (Araújo and Luoto, 2007). The integration of different variables that are related to biotic components can generate more advanced predictions than do simple SDM models. Some previous work has included biotic variables either explicitly in the modeling process or a posteriori by overlaying predictions of species that are known to interact (e.g., Araújo and Luoto, 2007; Hof et al., 2012; Giannini et al., 2013; Friedland et al., 2021; Davis et al., 2022). Evaluation of changes in overlap of interacting species can improve SDM predictions and make projections more realistic (Guisan and Thuiller, 2005; Schweiger et al., 2012), an approach that has seen increasing applicability to marine environments (Jones et al., 2013; Rockwood et al., 2020).

In this context, understanding the implications of predicted habitat suitability shifts of organisms that may interact in the natural environment is important. We suggest that a more thorough assessment needs to be made at the ecosystem scale (i.e., not simply at the individual species level) to anticipate future ecological changes. To this end, the objective of this study is to predict the distribution of an interconnected set of benthic native and non-native primary producers and primary and secondary consumers in the context of a warmer and invaded Arctic. The general hypothesis behind this study is that predicted shifts in native species habitat suitability will lead to a mismatch of regions where producers and consumers overlay, and that the free habitat left by native species will become available for non-natives. Non-native species will also be able to find suitable areas where natives are distributed, potentially increasing the chances of competition or other interactions. Projected suitable habitats are overlapped to predict spatial changes (e.g., match-mismatch) of interacting trophic groups to better understand the potential for future effects of climate change and biological invasions on Arctic benthic species.

Canada has the longest coastline in the world, primarily in the Arctic, covering approximately 2 million km² (Archambault et al., 2010). Most of this coast is rocky, with shorelines dominated and eroded by waves and affected by ice scour (Stewart and Lockhart, 2005; Frederick et al., 2016; Nyberg and Howell, 2016). This study focuses on the Eastern Canadian Arctic, delineated by 4 marine ecoregions, according to Spalding et al. (2007; Figure 1): Hudson Bay Complex (HBC), Baffin Bay-Davis Strait, Lancaster Sound, and Northern Labrador. Of these, HBC is the largest and is composed of the following subregions: James Bay, Hudson Bay, Hudson Strait, and Foxe Basin. HBC is relatively shallow (150 m mean depth in Hudson Bay) and characterized by a significant input of freshwater from tributary rivers and an important outflow of Arctic marine waters that pass through the system (Prinsenberg, 1986; Stewart and Lockhart, 2005; and references therein). The region north of the Labrador Sea through central Davis Strait and Baffin Bay (west into Lancaster Sound) is 650 m deep, on average, the deepest in the study area (Tang et al., 2004; Jørgensen et al., 2005).

The Eastern Canadian Arctic is known to have distinct and irreplaceable habitats due to the unique combination of characteristics, where sea ice plays an integral role in shaping the general ecology of the area, affecting the biota associated with the air-ice-water interface (Stewart and Lockhart, 2005). The region is characterized by large eelgrass beds and kelp forests that stabilize sediments and mitigate wave action, while also serving as shelter and feeding areas for a diverse range of marine biota (Stewart and Lockhart, 2005; Filbee-Dexter et al., 2022).

Native species were selected according to several characteristics. These included their known importance in coastal ecosystems as a source of habitat or as a food item for other organisms and/or their importance as a dietary item to northern indigenous peoples in the HBC region (Egeland et al., 2009; Landry et al., 2018; Filbee-Dexter et al., 2019; Landry et al., 2019; Murphy et al., 2021; Filbee-Dexter et al., 2022). The focus was on ecologically interacting species from linked trophic categories: primary producers (PP), primary consumers (PC), and secondary consumers (SC). Care was taken in the selection of native species to include at least 2 species in each category and that a sufficient amount and quality of information on these species (taxonomy, known distribution, known biological interactions) were available in the literature. These selection criteria led to the choices of 4 kelp species and 1 eelgrass as PP, 2 intertidal amphipods as PC, and 2 sculpins as SC (Tables 1 and S1).

The ecological interactions between these species are outlined graphically in Figure 2: kelp and eelgrass are known to be one of the food items consumed by intertidal amphipods (mainly when they degrade into particulate fractions; Harrison, 1977; Hop et al., 2002; Crawley and Hyndes, 2007; Legeżyńska et al., 2012; Murphy et al., 2021) and also provide habitat and shelter for these species (Dean et al., 2000; Kotta et al., 2010; Murphy et al., 2021). Likewise, these PP act as refuge and as a source of different prey items for small benthic fish such as sculpins (Dean et al., 2000; Murphy et al., 2021). Finally, sculpins feed on different prey items, and although their diet is diverse, intertidal amphipods are among their top known food items in Arctic ecosystems (Dick et al., 2009; Cui et al., 2012; Gray et al., 2017; Landry et al., 2018; Pedro et al., 2020; Figure 2).

Once native species were selected, the second step involved selecting non-native species that can be considered analogues of native ones. The selection criteria were (i) taxonomic proximity (species that are members of the same family for amphipods and sculpin, and same class or kingdom for kelp and eelgrass), (ii) invasion history and risk (species known to be invaders or to be expanding their range northward, or species identified as posing a relatively high risk for the region with respect to likelihoods of invasion and impact), and (iii) species with a high potential for introduction and survival/establishment in the Canadian Arctic (tolerance to environmental conditions and presence of introduction vectors; Mathieson et al., 2010; Wisz et al., 2015; Goldsmit et al., 2020; Goldsmit et al., 2021a; Goldsmit et al., 2023). Based on these criteria, 4 PP, 1 PC, and 2 SC non-native analogue species were identified (Tables 1 and S1).

Information on the known distribution (latitude and longitude) of each species of interest was compiled from various sources, mainly biodiversity databases such as the Global Biodiversity Information Facility (GBIF, www.gbif.org) and Ocean Biogeographic Information System (OBIS, www.obis.org). A thorough literature search was also performed to include more detailed information and fill distribution gaps for each species. The search was performed using Google Scholar during 2020, with the following main key words: “scientific name of the species” AND “distribution” AND “range.” As this review was not systematic, the search was completed, as needed, by adapting the search to include other terms (e.g., name of the species and a particular country or region of interest). Data points from the literature were used to fill gaps (Lee, 1980; Carlton and Scanlon, 1985; Josselyn and West, 1985; Borum et al., 2002; Bulleri and Airoldi, 2005; Provan et al., 2005; Chavanich et al., 2006; Kelly et al., 2006; Scheibling and Gagnon, 2006; Mathieson et al., 2008; Piscart et al., 2008; Semushin and Novoselov, 2009; Cheang et al., 2010; Drouin et al., 2012; Sabour et al., 2013; Sfriso and Facca, 2013; Armitage et al., 2014; Olesen et al., 2015; El Atouani et al., 2016; Hop et al., 2016; Alfonso et al., 2018; Schoenrock et al., 2018; Tempestini et al., 2018; Filbee-Dexter et al., 2019; Wilson and Lotze, 2019; Krause-Jensen et al., 2020; Ronowicz et al., 2020; Goldsmit et al., 2021b; Filbee-Dexter et al., 2022). Points that were uncertain (e.g., only one point in an isolated location) were verified and removed when confirming the original source of information was not possible. Details on the complete set of occurrence points prior to data cleaning and the list of sources used are provided in Table S2.

Duplicate records were removed to prevent overprediction, retaining only one record per grid cell (5 arcmin, the same resolution as environmental layers; García-Roselló et al., 2015). Following these steps, the final number of occurrence points per species varied from 43 to 2,394 (Table S3). To avoid including spatially autocorrelated points, the Spatial Filtering by Climatic Heterogeneity approach was used. This method uses a process whereby the scale of environmental variability in an area determines the scale at which the distribution data are considered, such that areas of low environmental heterogeneity are filtered at large spatial scales and areas characterized by high environmental heterogeneity are filtered at smaller spatial scales (Brown et al., 2017) in ArcGIS (ArcMap v10.7.1, SDMtoolbox; details on this procedure are given in Table S3). Because many of the species included in the study are coastal, there were a few cases where occurrence points were located within inland grids. Hence, the complete set of points was verified to include only points distributed in sea grids, except for points that were within 2 grid cells inland; these points were relocated to the closest sea grid using the Near Proximity tool in ArcMap. Points that were further inland were removed. The cleaned occurrence data points can be found in the Github repository (https://github.com/robwschlegel/HBCproject).

Environmental information (global marine data layers) was downloaded from Bio-ORACLE v2.1 (http://www.bio-oracle.org/; Assis et al., 2018). The resolution of these data layers was 5 arcmin (approximately equivalent to 9.2 km at the equator). They represent climatologies of monthly averages from a combination of satellite and in situ measurements over a 15-year period (2000–2014) that were used to characterize “present” environmental conditions in the context of this study (Assis et al., 2018). Mean surface values were used for all variables selected (except for temperature, which was species-specific; see “sensitivity approach” explained below). Temperature, salinity, ice thickness, and chlorophyll a were used for all species. To this list, nitrate was added for fish and primary producer models, while photosynthetically available radiation (PAR) and iron were added for primary producer models. All predictors were selected according to organism type based on knowledge that they are biologically and ecologically important variables and have been used in previous similar modeling studies (Barnes, 1999; Cusson et al., 2007; Zacher et al., 2009; Belanger et al., 2012; Gallardo et al., 2015; Wisz et al., 2015; Bosch et al., 2018; Goldsmit et al., 2018; Goldsmit et al., 2020; Goldsmit et al., 2021b). Correlations between variables were evaluated, and those that were highly correlated (≥0.7) were removed (Dormann et al., 2013) as their inclusion could lead to erroneous predictions and reflect misleading physiological tolerances (Marcelino and Verbruggen, 2015).

A sensitivity analysis was performed to determine the best temperature layer to include (i.e., bottom versus surface) for each species (Goldsmit et al., 2021b). This approach allowed for the analyses to be more precise because temperature is a primary variable that limits species distributions (Spence and Tingley, 2020), is usually the best explanatory variable for SDMs (Bosch et al., 2018), and is the main variable driving measures in climate change scenarios (Intergovernmental Panel on Climate Change, 2022). To determine the best temperature layer for a given species, mean and long-term maximum and minimum values of bottom and surface temperature layers were considered. This sensitivity analysis was performed as per Goldsmit et al. (2021b), where individual generalized linear models (GLM) were run for each species with each of the above-described temperature layers. Model prediction accuracy was assessed via the true skill statistic (TSS; further details below; Allouche et al., 2006). For a given species, the GLM with the highest TSS was considered to perform best, and the associated temperature layer was selected for inclusion in SDMs. Further details on this approach are provided in Tables S3 and S4.

Bathymetric information was obtained from Aquamaps (Kaschner et al., 2016) at 30 arcmin resolution (approximately 56 km at the equator). As bathymetry is a distal variable (it will have an indirect effect on species distribution and will covary with other variables such as PAR or temperature), it was used only as a mask once models were run to delimit the threshold for the maximum depth known for each species (Guisan et al., 2017; Goldsmit et al., 2018; Goldsmit et al., 2020; Goldsmit et al., 2021b; Table S3).

The collated information was used to build predictions of habitat suitability for each species using ensemble models with the biomod2 v3.4.13 package in R v4.0.3 (Thuiller et al., 2021). Ensemble models are generated by combining multiple modeling techniques that are then aggregated in a final mean prediction. This method has been shown to have lower bias than using single models (Araújo and New, 2007). A total of 5 modeling algorithms were used in the ensemble calculations for this study: GLM, random forest, artificial neural network, generalized additive model (GAM), and multivariate adaptive regression splines (more information on these models can be found in Thuiller et al., 2021). In addition to environmental data, these models require presence/absence data related to species distributions. Species occurrence data points represent the required presence data; pseudo-absence data were generated randomly to represent the absence data points needed to run the models. A total of 5,000 pseudo-absence points were generated randomly for each species, based on adapting default features (normally 10,000 random pseudo-absence points are generated but this study only concentrated on the Northern Hemisphere, and so only 5,000 random pseudo-absence points were generated). This modification of default features was used to adapt to the present case study, where existing recommendations are given based on model type and the number of presence data points (Barbet-Massin et al., 2012). To prevent sampling bias, 3 replicate arrays of random pseudo-absence points were generated (Guisan et al., 2017). Further information on this process can be found in Table S3. Model evaluation was performed using cross-validation, where models were trained with 70% of the occurrence points and the remaining 30% used for validation. Each model type was run 5 times and all replicates from all models were averaged to generate a single prediction for the final outcome. These predictions were then converted to binary (i.e., suitable/not suitable) predictions. This conversion was done based on the maximum training sensitivity plus specificity threshold considering that priority was given to the proportions of correctly predicted presences (sensitivity) and absences (specificity) during model training. This method has been shown to produce accurate predictions (Jiménez-Valverde and Lobo, 2007; Liu et al., 2013). Thresholds used for each species are provided in Table S5. To assess prediction accuracy, several evaluation metrics were used: TSS, receiver operating curve, kappa, accuracy, and sensitivity (more information on these metrics can be found in Table S5). Only statistically well-performing models (i.e., those with TSS ≥0.7) were retained for final ensemble predictions (Allouche et al., 2006; Guisan et al., 2017).

A detailed description of the modeling steps and input data are provided in Table S3 following the ODMAP protocol (Overview, Data, Model, Assessment and Prediction; Zurell et al., 2020). Visualization of results and figures were made with ArcGIS (ArcMap v10.7.1), Excel Microsoft 365, and Canva (Free version, 2023).

Once models were run, tested, and validated under current environmental conditions, predictions were made based on projected future environmental conditions simulating climate change scenarios. Layers representing the most extreme climate change scenarios (RCP 8.5) were obtained from Bio-ORACLE v2.1 (Assis et al., 2018) for years 2050 and 2100. This scenario predicts that, under greenhouse gas emission without climate policies, projected temperature anomalies of 4.9°C (relative to preindustrial levels) are expected by the end of the century (Riahi et al., 2007; Moss et al., 2010). This scenario also foresees a complete open water period in the Arctic during the summer by 2050 (Hwang et al., 2020), and was chosen to visualize predictions that would provide non-native species a greater chance of finding suitable habitat in the area. Available layers included in the analysis of future projections were temperature, salinity, and ice thickness. The other variables used to build models in the previous step remained unchanged.

Binary habitat suitability based on current and future projections of all species were used to analyze predicted species area overlaps between trophic groups (PP, PC, and SC) and species with differing status (native versus non-native) to understand potential spatial overlap. First, spatial predictions were overlaid to create a weighted sum, with the area covered calculated to identify all regions where one or more species coincided in the same pixel. Area overlap was determined using tools from the Spatial Analyst toolbox in ArcMap. Subsequently, other quantifications of spatial overlap were applied following Carroll et al. (2019), where spatial affinity and spatial niche similarity were calculated for all species type combinations and times (present, 2050, 2100). “Spatial affinity” measures the spatial independence of occurrence distributions (affinity between two distributions) and is calculated using the Bhattacharyya coefficient, a statistical method to quantify density functions for spatial use and their spatial independence as a null assumption (Bhattacharyya, 1943). This metric was applied to the predicted suitable habitats among native trophic groups. In contrast, “spatial niche similarity” measures spatial independence by calculating how similar are the spatial niches occupied by different species. It is calculated by applying Schoener’s D metric, a statistical method to determine the equivalence of spatial niches occupied by species pairs (Schoener, 1970). This metric was applied to the predicted suitable habitats between native and analogue non-native species. Both metrics range from 0 to 1, with higher values interpreted as species having higher spatial affinity and spatial niche similarity, respectively (for more information, refer to Carrol et al., 2019). These metrics were calculated in R v4.0.3 following the script provided in Carroll et al. (2019), with the assumption that sharing the same space indicates that species potentially interact, given that the analyzed species are known to interact in trophic relationships and non-native species are analogue to native species and known to use similar resources and habitats (Pinkster, 1975; Schebling and Gagnon, 2006; Landry et al., 2018; Figure 2 and Table S1).

Temperature and chlorophyll a had the greatest influence on habitat suitability for the 3 types of trophic groups for both native and non-native species, ranging from 13% to 74% of variable importance. These variables were followed by ice thickness and salinity, which ranged in importance from 4% to 28% (Figure 3A). The other variables (nitrate, iron, and PAR) ranged between 6% to 20% of average importance for their respective groups (Figure 3A).

Averaged response curves were included to analyze general patterns among and between sets of species. A clear observation was that, as expected for cold-adapted organisms, the temperature range for native Arctic species was much narrower than for non-native ones. Another pattern is related to ice, such that native species had higher probabilities of occurrence related to higher levels of ice thickness when compared to non-native species (Figure 3B). However, considering that these patterns are generalities, there were also species within a few trophic groups that showed contrasting response curves. These results show that some organisms may have a narrower niche with respect to certain environmental factors than others, despite having similar ecological roles (e.g., Laminaria solidungula and Zostera marina, related to temperature response curves in Figure S1.1A). Further detail on specific response curves for each species and the complete set of variables used are provided in Figure S1.1–S.1.6).

In general, the performance of ensemble models was good, with only a few models (mainly related to GAM for some primary producers) with TSS <0.7. The complete set of values obtained for each evaluation metric across species is shown in Table S5 and Figure S2, which provides a graphic demonstration of the range of scores obtained for each model and species, highlighting the >0.7 TSS threshold.

A total of 120 maps of binary predictions resulted from the analysis of area overlap of the complete set of ecologically interacting species and time scenarios. Figure 4 shows a summary of different overlap schemes for native species with habitat suitability predictions in the present versus the end of the century (for the complete set of predictions through time, refer to Figure S3.1–S3.3). The proportion of suitable habitat for different groups and species may be expected to vary over time. Increases or decreases in the overlap of the 3 trophic groups through time (PP + PC + SC) were predicted to trigger different types of habitat suitability overlap responses. When the 3 trophic groups overlapped and were predicted to change through time (increase or decrease by the year 2100), the other 2 trophic groups were predicted to do the opposite (see arrow directions in Figure 4). Most predictions of overlapping groups maintained a similar spatial extent through time (Figure 4A–D), although the combination between certain pairs of species resulted in the prediction of a drastic reduction of the spatial extent in the HB region (e.g., as depicted by the size of the arrows in Figure 4E). Resulting overlaps between trophic groups at multiple time scenarios (present, 2050 and 2100) are presented in Figure S3.1–S3.3.

Analyses were pushed further to identify variation in patterns of area overlap between interacting species. Three categories of habitat suitability predictions were identified through time: (i) suitability increased: Zostera marina, kelp species and sculpin group 1 (Gymnacanthus tricuspis and Myoxocephalus scorpius); (ii) suitability decreased: Laminaria solidungula and sculpin group 2 (M. scorpioides and M. quadricornis), and (iii) suitability remained stable (amphipods). Area of suitable habitat overlap of the 3 trophic groups (PP + PC + SC) was projected to increase through time for predictions related to Z. marina (Figures 4C and S3.1) and kelp with sculpin group 1 (Figures 4B and S3.2) throughout the HBC ecoregion. This increase varied among species: Z. marina was predicted to have a much higher proportion of habitat suitability change through time (Figure 4C) relative to kelp species (Figure 4B). However, in both cases, the increase in overlap of the 3 trophic groups was detrimental to overlapping trophic groups (e.g., decrease of PP + PC in Figure 4B, decrease of PC + SC in Figure 4C). In contrast, overlap of the 3 trophic groups decreased for predictions related to L. solidungula (Figures 4D, 4E, and S3.3) and the 3 kelp species combined with sculpin group 2 (Figures 4A and S3.2B), with HBC being most affected by this mismatch. Given that the pattern of shifts in suitability of habitat through time for sculpin groups are in opposition (group 1 predicted to increase and group 2 to decrease), all overlapping combinations including these groups were affected by these differences (Figure 4A and E). Similarly, the predicted decrease of L. solidungula affects not only the proportion of area overlap but also trophic combinations involving this species (Figures 4D, 4E, and S3.3). This species is predicted to have an almost complete disappearance of suitable habitat in Hudson Bay and Hudson Strait, leaving few suitable coastal regions around Foxe Basin and Baffin Bay-Davis Strait, with even fewer areas in Northern Labrador. Due to this mismatch, the most common predicted overlap with this combination includes only the 2 groups of consumers (PC + SC) that cover a similar spatial extent, largely centered in HBC (Figure 4D), or a considerable reduction in spatial extent in the study region (Figure 4E).

A total of 90 maps of binary predictions resulted from the analysis of area overlap of the complete set of native and analogue non-native species and time scenarios. The analysis was pushed further to identify patterns of area overlap between native and non-native species. Three different categories were found related to habitat suitability predictions for native and non-native species through time (i) decreased predicted suitable habitat for the native species with a concurrent increase for non-native species, (ii) increase in predicted suitable habitat for native and non-native species, or (iii) no (or very little) predicted non-native habitat suitability in the assessment region (Figures 5 and S4). A few clear examples of these types of overlaps are presented here, but given that many were repetitive, the remainder are shown in Figure S4. In addition, overlapping maps of non-native primary producers (except for Dumontia contorta) were not generated, as suitability was only predicted in very limited regions and only in a few future scenarios.

Opposing predicted shifts in suitable habitat of native versus non-native species (i.e., a predicted decrease in the former and an increase in the latter) were generally restricted to combinations including native species that are expected to decrease their range under climate change scenarios (Laminaria solidungula and sculpin group 2; Figures 5A and S4.3B). In these situations, non-natives could fill open niches left by the native analogue species. In contrast, when both native and non-native species habitat suitability was predicted to increase (as is the case for Zostera marina, kelp, amphipods, and sculpin fish group 1), the area overlap increased through time such that native species could experience competition with their non-native analogue species (Figures 5B, S4.1–4.3A). For PP and PC, the most affected area of overlap with non-native species or loss of suitable habitat for analogue native species was concentrated in the HBC (Figures 5, S4.1, and S4.2), whereas the area of suitable habitat loss for native SC extended through the complete Eastern Canadian Arctic region (Figure S4.3).

Spatial affinity and niche similarity (Schoener’s and Bhattacharyya’s coefficients, respectively) provide a measure of the differences in space use by different populations (but not the strength of interaction in that space). The spatial affinity of prey and predators (PP and PC as prey with PC and SC as predators, respectively) increased substantially in relation to Zostera marina and amphipods through time, whereas it decreased noticeably between Laminaria solidungula and amphipods (Figure 6A). These predicted changes are clearly related to the respective predicted increase and decrease of these primary producers. Spatial niche similarity of native and analogue non-native species was higher for SC than for PP and PC (Figure 6B). In addition, this relationship decreased or slightly increased for SC, whereas it increased in higher proportions through time for the other 2 groups (Figure 6B).

This study highlights that the direction of predicted shifts in native Arctic species distributions varies primarily as a function of species identity (for the specific species under consideration) rather than trophic level. Suitable habitats for individual native species are predicted to shift independently (and in different ways) with expected variations in the environment, driven by climate change. In contrast, the overall suitable habitat of each trophic group is predicted to remain more stable through time, even though some species within those groups may shift (some may cope better with the new environmental conditions than others, as reflected by habitat suitability predictions). This trophic group stability may result in positive outcomes for overall ecosystem functioning, in general, although the match or mismatch between individual species may vary as a function of identity leading to changes in food web dynamics. Likewise, findings from this study demonstrate the possibility for future spatial overlaps (or lack thereof) between native and analogue non-native species that have the potential of being introduced to the Arctic. Such shifts through time could enable novel types of interactions (e.g., competition) due to the coincidence of habitat suitability in the same space and similarity in spatial niche, or replacement of native species by non-native analogues in cases where a reduction in spatial overlap of suitable habitats is predicted.

Temperature is the most important environmental variable for predicting the distribution of species in the marine environment (Bosch et al., 2018). In the context of Arctic climate change, where temperature increases are nearly certain, species identity and trophic group, together with their known thermal range/tolerances, can help to identify those that are most vulnerable. For example, contrasting 2 species in the same trophic group (PP: Laminaria solidungula and Zostera marina) shows that the range of the response curve of L. solidungula is much narrower than that for Z. marina, which translates into predictions of habitat loss and gain, respectively. L. solidungula has been identified as a cryophilic species (Bringloe et al., 2022) predicted to be sensitive to temperature changes, which will result in habitat loss, particularly in the south (Goldsmit et al., 2021b). However, species identity may influence expressed patterns and combinations of habitat change/overlap between trophic groups (e.g., predicted patterns within the trophic groups PP and SC varied with species composition). In addition to being a direct stressor, temperature change in Arctic environments can also trigger indirect changes. For example, decreased sea ice cover (which influences primary production by ice algae alters light and nutrient availability, etc.) can alter trophic dynamics due to changes in food availability (Bartsch et al., 2016; Bonsell and Dunton, 2018). Changes in the distribution of predators or prey may also directly affect species distributions of trophically linked species (Louthan et al., 2015; Ellingsen et al., 2020).

There was a marked difference in the response curves of native and non-native species related to sea ice thickness. Not surprisingly, native species from polar regions have a greater probability of occurrence in the presence of sea ice. This effect is supported by field observations where amphipods from Arctic coastal locations are often associated with sea ice because it provides them with protection from predators such as sculpin (Hop et al., 2002; Landry et al., 2018). Sea ice may also act as a physical barrier, restraining the survival and establishment of new species (e.g., introduced non-natives). However, as temperature increases, sea ice thickness, cover, and duration will inevitably be impacted (Derksen et al., 2019; Hwang et al., 2020).

Chlorophyll a was identified as one of the 3 most important variables for predicting the distribution of appropriate habitat for all species evaluated, contributing between 10% and 30% of average importance for the complete set of species evaluated. Other studies have shown chlorophyll a to contribute similarly when there are several predictors included in models and to have a greater influence on benthic polar and temperate species (Bosch et al., 2018; Guillaumot et al., 2018). The relative importance of this environmental variable may be related to the fact that chlorophyll a is a proxy for primary productivity and often associated with areas of upwelling where water is rich in a wide spectrum of nutrients and where organisms may congregate. This variable can thus provide important information on productivity and its influence on the distribution of native and non-native species in addition to improving model predictions (Gallardo et al., 2015; Goodyear et al., 2017).

There were 3 different categories of change in habitat suitability predictions through time: increase, decrease, or remain stable. Predicted mismatches between species in the various trophic groups can have varying impacts, depending on species life history traits and abilities to cope with a given change. For instance, organisms that are already found at their temperature extremes (i.e., poleward edges of their distributions) may shift with unpredictable patterns (no shift at all, greater shifts, or shifts in the opposite direction than expected) when facing temperature increases, as shown by Fredston et al. (2021) using a 5-decade dataset for fish and invertebrates. Analyzing potential changes of the selected species is important not only to the ecosystem (as a source of habitat or as a food item; Landry et al., 2018; Filbee-Dexter et al., 2019; Landry et al., 2019; Murphy et al., 2021; Filbee-Dexter et al., 2022) but also as dietary items to indigenous peoples in the HBC region (Egeland et al., 2009).

Organisms with more limited distributions and endemic to extreme polar environmental conditions may have limited ability to cope with future thermal increases as they have adapted physiologically to these extreme conditions (e.g., Gordillo et al., 2022). While species with greater mobility at the adult stage (e.g., fish) may be able to adapt behaviorally by moving to more suitable habitats under conditions of environmental change, sculpin are largely sedentary (Landry et al., 2019; Barton et al., 2020; but see Hermann, 2021). Even though most species included in this study are not likely to migrate great distances, the sculpin Myoxocephalus quadricornis is known to have a flexible diet and may be able to adapt behaviorally by exploiting ephemeral resources (Hermann, 2021). However, a decline in the distribution of this species has been observed in higher Arctic regions where a borealization of fish species is already happening (von Biela et al., 2022). Hence, the impact of mismatches between certain species may be reduced if there are potential matches with other species.

In contrast to more polar-adapted organisms, widely distributed species may be able to cope better with climate change. Such is the case of Myoxocephalus scorpius, a species with a wide thermal tolerance that is distributed from temperate to Arctic waters (Robins and Ray, 1986; Grans et al., 2013) and was predicted to increase its potential range of suitable habitat in this study under a climate change scenario. While amphipods are also widely distributed, their habitat suitability was predicted to remain nearly stable through time. They are versatile feeders and an important link between primary and secondary producers (Bradstreet and Cross, 1982; Grebmeier and Harrison, 1992), therefore providing stability in the system.

Like Myoxocephalus scorpius, Zostera marina was predicted to have considerable increases in suitable habitat under climate change scenarios. Growth of marine vegetation in the Arctic is generally stimulated by warming and reduced sea ice cover (Krause-Jensen and Duarte, 2014). Z. marina has a wide distribution in temperate to Arctic regions, and exposure to higher temperatures is likely to stimulate growth because its optimum growth temperature is higher than that of endemic Arctic species (Müller et al., 2009; Olesen et al., 2015; Beca-Carretero et al., 2018; Wilson and Lotze, 2019). However, particularly in the Hudson Bay region, eelgrass cover and health has declined, perhaps linked to blooms of microorganisms related to global warming, and the health of this species is impacted by reductions in salinity with summer freshwater discharge, early ice breakup, reduced water clarity, and overgrowth by other types of marine vegetation (Short, 2019; Leblanc et al., 2022). Given these observations, together with the fact that sediment type was not included as an environmental predictor in the model (following normal practice for this type of study, given the lack of high-resolution data; see, for example, Goldsmit et al., 2018; Goldsmit et al., 2020; Goldsmit et al., 2021b; Assis et al., 2022), predicted habitat suitability of eelgrass in the present study could be overestimated. According to Nyberg and Howell (2016), the regional focus of this study is largely erosional (rocky) and not the type of substrate generally associated with eelgrass. Even though eelgrass is typically found on sandy and muddy bottoms in regions that are less exposed, it has been observed sporadically on a range of substrate types, including cobble (Dean et al., 2000).

Pinsky et al. (2020) suggested that Arctic ecosystems may be particularly susceptible to changes in species composition, leading to community reorganization and altered network properties (Kortsch et al., 2012; Kortsch et al., 2015). Thus, stable Arctic food webs will depend on retaining a heterogeneity of resources (McMeans et al., 2013). Species identity may change but the ecological roles of different groups (i.e., food web properties) will likely be maintained (Pinsky et al., 2020). However, considering how predicted changes will impact complete Arctic food webs is important. Climate change may affect seasonal predator-prey mismatches at higher latitudes (Durant et al., 2019) with upward and downward multi-trophic cascading effects. Modifications in the spatial distribution of energy-rich prey may have large implications for keystone species in ecosystems and for higher trophic Arctic predators, such as marine mammals (Florko et al., 2021). Even changes in abundance may have large effects on multi-trophic interactions, as outlined by Lorentsen et al. (2010), who found changes in kelp abundance affected seabird foraging efficiency. Lower trophic level consumers, such as sculpin, can couple numerous energy pathways; however, knowledge of the role of these generalists is scarce in many ecosystems, including the Arctic (McMeans et al., 2013). Sculpin were the secondary consumers considered in this trophic analysis and are known to contribute to both pelagic and benthic energy pathways, serving as important regulators of trophic levels via predation (McCann et al., 2005). They are also known to be very abundant and provide an important, seasonally stable prey source for larger predators (Lowry et al., 1980; Dyck and Romberg, 2007; Kortsch et al., 2015; Breton-Honeyman et al., 2016). In addition to these modifications, effects of species change may extend beyond where species are physically located. For example, blue carbon (produced by kelp) can be exported in high proportions to deeper waters (Krause-Jensen and Duarte, 2016; Krause-Jensen et al., 2018; Pessarrodona et al., 2023). Hence, a change in habitat suitability for kelp may greatly impact marine carbon sequestration in both coastal and deeper food webs.

Spatial affinity between trophic groups (species-shared preferences for specific areas) remained relatively stable, except for the case of Zostera marina and Laminaria solidungula, given that these 2 species were predicted to considerably increase and decrease, respectively, their future ranges. This index measures how equally predators and prey share available resources; hence, spatial overlap between these 2 species and their primary consumers may be compromised. This compromised overlap is not surprising given that changes in dominant PP can have an important effect on species composition and abundance of related benthic biota (Barnes and Hendy, 2015).

Analysis of spatial niche similarity (Fieberg and Kochanny, 2005) showed that most analogue non-native species are expected to share an increasing proportion of the suitable areas where native species are also projected to be distributed under future scenarios of climate change. This measure of area where spatial resources are predicted to be shared could translate into ecological interactions (e.g., competition) between pairs of analogue species (e.g., Diederich, 2006; Vye et al., 2015; Britton et al., 2018; Zwerschke et al., 2018). Interactions such as competition can be negative, although that may not necessarily always be the case. Facilitation has been shown in certain analogue benthic invertebrates, where invasion success is not only determined by environmental conditions but also by other factors, including habitat type (Zwerschke et al., 2018).

An exception to the increased spatial niche similarity is seen for sculpin group 2 (Myoxocephalus scorpioides and M. quadricornis), where the predicted decrease in overlap through time is likely due to predicted loss of suitable habitat for these 2 species. Resources are unlikely to be shared spatially, as these species are not predicted to find appropriate suitable habitat (as per SDM results) in areas where their non-native analogues could expand.

The match/mismatch or overlap of predicted habitat suitability between native and analogue non-native species is of special interest given that the performance of cold-adapted native species often declines when exposed to higher temperatures or more stressful environmental changes, as opposed to the positive response that is often observed for non-native species (Kuprijanov et al., 2017; Torn et al., 2020; McKnight et al., 2021). Warmer temperatures will likely facilitate range expansion, arrival, survival, and establishment of non-indigenous species (Hellmann et al., 2008; Sardain et al., 2019; Essl et al., 2020), and the Arctic region is no exception (Jones and Cheung, 2015; Goldsmit et al., 2018; Goldsmit et al., 2020; Assis et al., 2022; Bringloe et al., 2022). Indeed, new shipping routes are expected to open in the region, increasing probabilities of species arrival to new areas (Mudryk et al., 2021) with trans-Arctic routes already in use under current conditions as they become safer to navigate at a much faster rate than model predictions (Cao et al., 2022). However, increased temperature and decreased ice cover does not mean that non-native species will not encounter further constraints. Decreased sea ice may increase light availability at the seabed, but this increase could be counteracted by increased run-off contributing to greater suspended sediments attenuating light (Bonsell and Dunton, 2018). Further, daylight is not present during winter at high latitudes, affecting photoperiods and potentially restricting non-native species (Gordillo et al., 2022). Another consequence of sea ice melting is changes in salinity, which could also impact biological processes of native and non-native species (McKnight et al., 2021).

The ecological impacts of invasive species in a simple trophic web, as the ones studied here, can be direct (habitat) or indirect (competition, predation, and grazing) and contribute to increased extinction rates and other significant ecosystem changes (Blackburn et al., 2014; Gallardo et al., 2016; Blackburn et al., 2019). The present study predicted many cases of spatial overlap, including amphipods, most PP (with the exception of Laminaria solidungula), and sculpin group 1 with their analogue non-native species. Many of the non-native species included in this study are known to have out-competed native species in their invaded ranges elsewhere (e.g., Curiel et al., 1998; Curiel et al., 2001; Orav-Kotta et al., 2009; Engelen et al., 2015; Jänes et al., 2015; Montgomery et al., 2022), reinforcing the notion that native Arctic species are not exempt from climate change-related incursions of non-native species. However, extinctions caused exclusively by competition are uncommon. Rather, extinction threats are explained most commonly by the combination of habitat loss with shifts in inter-trophic relationships (Davis, 2003). Competition-driven extinctions may take a long time to occur and may end up being affected by other factors such as different types of disturbances, diseases, or environmental fluctuations (Davis, 2003). The implication is that the native and analogue non-native species that were predicted to share the same space would likely continue competing for resources, with possible reductions in fitness of native congeners as opposed to extinction. Moreover, competition is not exclusively between native and non-native species, since native species may also compete among themselves for the same space (e.g., macroalgae compete with seagrass in coastal regions; Wang et al., 2021).

In certain cases, such as for Laminaria solidungula or sculpin group 2, where native species were predicted to have reduced suitable habitat and their non-native analogues to have increased suitable habitat, the latter theoretically could make use of the open niche space that was left by the native species. The theory of fluctuating resource availability in the context of invasibility (Davis et al., 2000) could play an important role such that unused resources may be used by non-native species. If native species decrease their range, then analogue non-native species could take advantage of the available resources to facilitate their establishment. From an ecosystem perspective, if a given native species is predicted to decrease its range due to its sensitivity to temperature change, analogue non-native species could serve equivalent or similar roles in the ecosystem. There may also be positive implications for trophic cascades where invasive prey ends up being abundant and available to native predators (Liversage et al., 2021). An example is the case of Gammarus tigrinus in the Baltic Sea which, although invasive to the area, has become an important prey item for native fish (Daunys and Zettler, 2006; and references therein). However, this situation may be oversimplified as many other types of interactions could occur. For example, different kinds of macroalgae vary with respect to the food or shelter they provide to other organisms (e.g., blade width or type of cover may affect associated organisms) with potential implications for associated fish distributions (Dean et al., 2000). Trophodynamics and energy budgets may also be affected, as was observed following the introduction of the invasive seaweeds to the diets of herbivorous amphipod species (Vázquez-Luis et al., 2013), considering that amphipods prefer certain types of phytobenthos (Goecker and Kåll, 2003; Crawley and Hyndes, 2007) and vary in their grazing rates (Orav-Kotta et al., 2009). In addition, the fecundity of non-native species may differ from that of native ones (e.g., Jänes et al., 2015), which could affect interactions with, and the abundance of, associated biota.

Several modeling algorithms were used for the present analysis to produce more robust predictions based on the averages of multiple modeling techniques. All models have their own limitations, and these are important to keep in mind when interpreting results. However, ensemble models are known to be powerful due to the reductions in the biases that each individual algorithm can create when making future projections (Araújo and New, 2007). Some of the limitations are data-related, specifically with respect to abiotic climate variables (e.g., temperature, salinity, ice) in remote regions. The variables used were obtained at a global scale and are considered to have a good general coverage. However, some places, such as Arctic coastal regions, are characterized by limited direct measurements as they are logistically difficult to obtain. Hence, there is a relatively high uncertainty, as is the case for some nutrients and biophysical variables (e.g., nitrate and chlorophyll a), for which there is no complete high-resolution cover around the globe (Assis et al., 2018). Another source of uncertainty was the lack of future projections for the complete set of environmental variables. While this lack is common in these type of modeling exercises (e.g., Jueterbock et al., 2016; Goldsmit et al., 2018; Goldsmit et al., 2020; Goldsmit et al., 2021b), certain variables (e.g., nutrients and minerals, chlorophyll a, and PAR) are assumed to remain stable through projections in time, even though some will likely vary considerably in the future (e.g., Steinacher et al., 2010; Frölicher et al., 2016).

Among the main critiques of SDM are that biotic interactions are not included in the analyses and that the models do not represent the realized distribution of organisms because they assume an equilibrium with their environment (Guisan and Zimmermann, 2000; Bellard et al., 2018). Biotic interactions, together with potential resources, are important factors that help shape the spatial distribution of species (Wisz et al., 2013; Louthan et al., 2015). The approach developed in this study combines independent predicted habitat suitability with known biotic interactions to evaluate the overlap between species. This approach facilitates analyzing species interchangeably and identifying factors that may affect individual organisms. It also enables more extensive analyses at a larger geographic scale. Similar approaches have been applied previously where species in a trophic relationship or where native and non-native species were modeled using other methods (e.g., Leidenberger et al., 2015; Herkül et al., 2016; Durant et al., 2019; Torn et al., 2020; Florko et al., 2021), but their corresponding habitat suitability predictions were not overlayed to analyze areas of match/mismatch, as was done here.

Climate change and invasive species have been identified as emerging issues and key drivers of biodiversity loss (Bellard and Blackburn, 2016; Ricciardi et al., 2017; Blackburn et al., 2019). Arctic marine fauna has already been affected by one or both of these factors (e.g., Durant et al., 2019; Mérillet et al., 2022). While one of the aims of the present study was to quantify the spatial overlap of predicted habitat suitability of native with analogue non-native species based on known ecological interactions, the results can also contribute to the development of hypotheses of expected ecological changes, amplifying modifications to trophic cascades, resource competition, and increased habitat suitability associated with the establishment of non-native species. The greatest risk to Arctic food web stability has been suggested to be the removal of resource heterogeneity (McMeans et al., 2013). Results from this study suggest that while species composition and identity may vary with future climate shifts and species introductions, trophic relationships and roles will likely be maintained over time, in agreement with the findings of other studies (McMeans et al., 2013; Pinsky et al., 2020). The approach developed here could be applied to a much greater number of species that are ecologically interconnected, including top predators, to analyze predictions at a larger scale. Studies on the predicted spatial match and mismatch of ecologically interacting species can help assess risks of climate change and the introduction of new species, while also informing management of potential modifications to biodiversity, species interactions, and ecosystem services at relevant spatial scales.

Sources of all occurrence points used to feed the model are included in supporting information, together with R scripts. The R scripts are also available via a GitHub repository: https://github.com/robwschlegel/HBCproject.

The supplemental files for this article can be found as follows:

See Tables S1–S5 and Figures S1–S4 attached as Supplementary Material.

The authors would like to thank Dr. Gemma Carroll for discussions and guiding them on the application of the overlap indexes that were used. They are grateful to Colin Jubinville for technical support and his incredible willingness to help recover information related to this project after facing major technical issues. The authors would like to acknowledge the ArcticKelp project from ArcticNet (special thanks to Karen Filbee-Dexter, Kathleen McGregor, Christopher J. Mundy, Ladd Johnson, and Philippe Archambault) for the invaluable kelp data that they collected in the field and to the Canadian Museum of Nature (Amanda Savoie) for the museum data that was also used to feed the models. They thank Kathryn Shoenrock, Iris Hendriks, Marta Ronowicz, and Gary Saunders for sharing species occurrence information. They would like to thank Anthony Asselin for his help on collating and verifying some of the occurrence records that were used in this study. They appreciate comments and exchanges performed during the last couple of years with all the members of the Hudson Bay Change project. This work benefited from feedback received during presentations of draft results to the International Council for Exploration of the Seas (ICES) Expert Working Groups on Ballast Water and Other Shipping Vectors and Introductions and Transfers of Marine Organisms. They appreciate the time and comments of the 2 reviewers; thanks to them the final version of the manuscript has been considerably improved. The author sequence follows the “first-last author-emphasis” norm.

Fisheries and Oceans Canada Results Fund and Aquatic Invasive Species Program Funding.

The authors declare that there are no competing interests.

Contributed to conception and design: JG, CWM, DD, KLH.

Acquired funding: KLH, DD.

Contributed to acquisition of data: JG.

Contributed to analysis and interpretation of data: JG, RWS.

Drafted the article: JG.

Revised and approved the submitted version for publication: CWM, RWS, DD, KLH.