What and where? Predicting invasion hotspots in the Arctic marine realm.

The risk of aquatic invasions in the Arctic is expected to increase with climate warming, greater shipping activity and resource exploitation in the region. Planktonic and benthic marine aquatic invasive species (AIS) with the greatest potential for invasion and impact in the Canadian Arctic were identified and the 23 riskiest species were modelled to predict their potential spatial distributions at pan-Arctic and global scales. Modelling was conducted under present environmental conditions and two intermediate future (2050 and 2100) global warming scenarios. Invasion hotspots-regions of the Arctic where habitat is predicted to be suitable for a high number of potential AIS-were located in Hudson Bay, Northern Grand Banks/Labrador, Chukchi/Eastern Bering seas and Barents/White seas, suggesting that these regions could be more vulnerable to invasions. Globally, both benthic and planktonic organisms showed a future poleward shift in suitable habitat. At a pan-Arctic scale, all organisms showed suitable habitat gains under future conditions. However, at the global scale, habitat loss was predicted in more tropical regions for some taxa, particularly most planktonic species. Results from the present study can help prioritize management efforts in the face of climate change in the Arctic marine ecosystem. Moreover, this particular approach provides information to identify present and future high-risk areas for AIS in response to global warming.

INTRODUCTION

Invasive species are considered to be one of the main drivers of biodiversity loss and recent extinctions (Bellard, Cassey, & Blackburn, 2016; Blackburn, Bellard, & Ricciardi, 2019). In addition, they may cause great ecological and economic impacts and, in some cases, even impact human health (Cook, Payne, MacLeod, & Brown, 2016; Hewitt, Campbell, & Gollasch, 2006; Ricciardi, Hoopes, Marchetti, & Lockwood, 2013; Simberloff et al., 2013). Estimates suggest that their presence may entail costs of up to 12% of the gross domestic product of affected countries (Marbuah, Gren, & McKie, 2014).

The greatest numbers of aquatic invasions have been documented in temperate regions (Ruiz & Hewitt, 2009). However, high latitude areas are generally warming at a faster rate than other areas (Niemi et al., 2019; Overland et al., 2019) and various climate change scenarios predict that the ice‐free season will continue to lengthen (Barnhart, Miller, Overeem, & Kay, 2016; Jahn, 2018; Niederdrenk & Notz, 2018; Sigmond, Fyfe, & Swart, 2018). If warming is limited to 1.5°C, as proposed in the Paris Agreement, ice‐free summers in the Arctic will continue to rise even if global warming ceases (Screen, 2018). Moreover, temperature extremes at regional scales could be much greater than the global mean with some models predicting a rise greater than 7°C in the Arctic with a global average increase of only 1.5°C at peak warming (Seneviratne et al., 2018). Aquatic invasions in Arctic ecosystems have thus been identified as an emerging issue (Ricciardi et al., 2017) as the area is becoming increasingly susceptible due to regional climate change and associated increases in shipping activity and resource exploitation (Melia, Haines, & Hawkins, 2016; Miller & Ruiz, 2014; Niimi, 2004; Smith & Stephenson, 2013).

One of the leading causes of increasing invasion risk is the expansion of transportation networks (Early et al., 2016). Aquatic transportation, in particular, is highly relevant given that the shipping network is responsible for 90% of global trade (Kaluza, Kölzsch, Gastner, & Blasius, 2010; Xu et al., 2014), with ballast water and biofouling being the main vectors for transport of organisms (Bailey, 2015; Drake & Lodge, 2007; Ruiz, Fofonoff, Steves, & Carlton, 2015). Consequently, marine taxa are one of the dominant groups that are unintentionally introduced at a global scale (Turbelin, Malamud, & Francis, 2017), and are predicted to contribute a 3‐ to 20‐fold increase in global invasion risk (Sardain, Sardain, & Leung, 2019). In the Arctic, almost half of known marine invasions have been due to shipping activities as a single‐pathway means of transport (Chan et al., 2019). Greater opportunities for aquatic invasive species (AIS) transport are expected in the future as shipping traffic is projected to grow given predicted ice‐free conditions that will allow more direct shipping corridors between the Pacific and the Atlantic oceans (Melia et al., 2016; Miller & Ruiz, 2014; Smith & Stephenson, 2013), reducing transit distances by up to 24% relative to current shipping routes (Buixadé Farré et al., 2014).

A recent study tallied a total of 34 non‐indigenous species for the entire Arctic, with shipping activity being one of the main purported vectors (Chan et al., 2019). In addition, non‐indigenous species have been identified from both ships en route to the Arctic and in the environment to which they arrive, highlighting the risk of new invasions in this region (Ashton, Riedlecker, & Ruiz, 2008; Chan, MacIsaac, & Bailey, 2015; Dispas, 2019; Gíslason et al., 2014; Golder, 2018; Laget, 2017; Lambert, Shenkar, & Swalla, 2010; Svavarsson & Dungal, 2008; Tremblay, 2017; Zimina, 2014). Furthermore, the origins of many newly observed species in Arctic locations are unknown so they have often been classified as cryptogenic (Goldsmit, Howland, & Archambault, 2014; MacDonald, Bluhm, Iken, Gagaev, & Strong, 2010). These numbers are high compared to the other polar environment—the Antarctic—where there are no established AIS populations to date and only five non‐indigenous marine species that have potentially arrived through human activities (McCarthy, Peck, Hughes, & Aldridge, 2019). However, Hughes et al. (2020) suggest that marine taxa pose the greatest potential invasion risk to the Antarctic Peninsula.

Given the vast size and remoteness of the Arctic, few baseline studies have been conducted (Archambault et al., 2010; CAFF & PAME, 2017; Piepenburg et al., 2011) and the resident fauna are thus poorly described relative to more accessible locations elsewhere in the world (Costello et al., 2017). Consequently, simply identifying newly arriving species can be challenging. To improve this detection gap, new methods and identification techniques are being used in remote regions, such as environmental DNA and metabarcoding (Brown, Chain, Zhan, MacIsaac, & Cristescu, 2016; Chain, Brown, MacIsaac, & Cristescu, 2016; Lacoursière‐Roussel et al., 2018; Leduc et al., 2019). Species distribution modelling (SDM) may be used to predict the distribution of suitable habitat for AIS (Barbosa & Schneck, 2015; Barbosa, Schneck, & Melo, 2012), including for high‐risk AIS in areas of concern, such as high latitude regions (Byrne, Gall, Wolfe, & Agüera, 2016; de Rivera, Steves, Fofonoff, Hines, & Ruiz, 2011; Goldsmit et al., 2018; Ware et al., 2016). Modelled results can inform preinvasion management policies to avoid new arrivals and, potentially, subsequent species establishments and consequent needs for eradication actions (Floerl, 2014; Locke & Hanson, 2009).

Although future invasions are expected to be enhanced by climate change (Bellard et al., 2013), predicted SDM‐modelled distributions of environmental suitability for invasive species have been understudied in this context (Bellard, Jeschke, Leroy, & Mace, 2018), particularly for marine ecosystems (Barbosa et al., 2012; Bellard et al., 2018). Nevertheless, SDM tools have been successfully applied to marine species (Lowen, McKindsey, Therriault, & DiBacco, 2016; Meißner et al., 2014; Reiss, Cunze, König, Neumann, & Kröncke, 2011; Robinson et al., 2011; Robinson, Nelson, Costello, Sutherland, & Lundquist, 2017; Weinert et al., 2016), including forecasting distributions under global change scenarios in high latitude regions (de Rivera et al., 2011; Goldsmit et al., 2018; Jueterbock, Smolina, Coyer, & Hoarau, 2016).

Regions where many species coincide, known as biodiversity hotspots, may provide information about species richness, endemism and/or threatened taxa, and may be of importance for conservation (Myers, Mittermeier, Mittermeier, Da Fonseca, & Kent, 2000). Based on this concept, potential invasion hotspots may be identified as areas that are hospitable to and thus at risk of invasion by exceptionally high numbers of AIS for a given region (Ruiz, Fofonoff, Steves, Foss, & Shiba, 2011). The use of a biological invasion hotspot approach has been limited to date (but see Bellard, Leroy, Thuiller, Rysman, & Courchamp, 2016; Catford, Vesk, White, & Wintle, 2011; Ibanez, Silander, Allen, Treanor, & Wilson, 2009; O'Donnell et al., 2012; Torres et al., 2018); with few studies in marine environments. Most analyses of AIS hotspots in marine ecosystems have focused on vectors (Davidson et al., 2010; Drake & Lodge, 2004; Pearce, Peeler, & Stebbing, 2012; Semmens, Buhle, Salomon, & Pattengill‐Semmens, 2004; Tidbury, Taylor, Copp, Garnacho, & Stebbing, 2016) and current AIS richness (Edelist, Rilov, Golani, Carlton, & Spanier, 2013; Katsanevakis et al., 2014; Kelly, Leach, Cameron, Maggs, & Reid, 2014; Rilov & Galil, 2009; Ruiz et al., 2011; Verlaque, 2001), although a few have also predicted hotspots of potential marine biological invasions (Cheung et al., 2009; Gallardo, Zieritz, & Aldridge, 2015; Jones & Cheung, 2015).

This study aims to identify potential high‐risk AIS and predict hotspots of invasion for those under current and projected future environmental conditions in Arctic environments. Cold‐tolerant planktonic and benthic AIS were scored to produce a list of the top species with the highest relative likelihood of invasion and impact, since these ecological groups include the dominant known highest risk marine invasive species (Molnar, Gamboa, Revenga, & Spalding, 2008). The distribution of suitable habitats for this set of species was then modelled to predict regions of overlap under current and future projected conditions, thus identifying hotspots of potential biological invasions. Although all analyses were done globally, there was an emphasis at the Canadian and pan‐Arctic scales.

MATERIALS AND METHODS

Study region

Distribution model outputs were produced with global coverage to allow evaluation of global and pan‐Arctic patterns although analyses were focussed on the potential highest‐risk species for the Canadian Arctic. The rationale for this approach is the vast expanse of the Canadian Arctic (Archambault et al., 2010), its vulnerability to invasion (Chan, Bailey, Wiley, & MacIsaac, 2013; Chan, MacIsaac, & Bailey, 2016; Goldsmit et al., 2018; Goldsmit, McKindsey, Archambault, & Howland, 2019), the ambiguous status of new species that have potentially arrived through human transport (Goldsmit et al., 2014) and the unprecedented warming that it has experienced over the last decade (Bush & Lemmen, 2019; Niemi et al., 2019). The region is immense, accounting for eight of the 19 Arctic marine ecoregions (Figure 1; Spalding et al., 2007). In fact, Canada has the longest coastline in the world, the majority of which is situated in the Arctic (accounting for almost 2,000,000 km2 of the territorial sea; Archambault et al., 2010). Even though few AIS have been identified in Canada relative to other Arctic regions, Chan et al. (2019) suggest that this accounts for ca. 20% of all AIS recorded from marine Arctic waters. Recently, additional species have been identified for the first time in the Canadian Arctic using molecular tools (Brown et al., 2016; Chain et al., 2016; Lacoursière‐Roussel et al., 2018) or identified as cryptogenic as their origin is unclear (Goldsmit et al., 2014). Moreover, Canadian Arctic ports (especially those in the Hudson Bay Complex; Figure 1) have been identified as being of moderate to high ecological risk of invasion since they could provide suitable habitat for various AIS (Goldsmit et al., 2018; Goldsmit, Nudds, et al., 2019) and there is evidence that a number of non‐indigenous species are being transported into the region by shipping (Chan et al., 2015; Dispas, 2019; Laget, 2017; Tremblay, 2017).

FIGURE 1

Arctic ecoregions as delineated by Spalding et al. (2007). Canadian Arctic ecoregions are coloured. 1, Hudson Bay Complex; 2, Northern Grand Banks‐Southern Labrador; 3, Northern Labrador; 4, Baffin Bay‐Davis Strait; 5, Lancaster Sound; 6, High Arctic Archipelago; 7, Beaufort‐Amundsen‐Viscount Melville‐Queen Maud; 8, Beaufort Sea‐continental coast and shelf; 9, Chukchi Sea; 10, Eastern Bering Strait; 11, East Siberian Sea; 12, Laptev Sea; 13, Kara Sea; 14, North and East Barents Sea; 15, White Sea; 16, North and East Iceland; 17, East Greenland Shelf; 18, North Greenland; 19, West Greenland Shelf

Species selection

A three‐step procedure was used to select potential marine AIS (zoobenthos, phytobenthos and zooplankton) for modelling (Figure S1). Most were species that are considered AIS in other regions of the world and may be transported via shipping to the Canadian Arctic; several others had been previously detected in other Arctic environments. These steps included: (a) prescreening analysis and selection; (b) ranking a subset of these species using an invasive species screening tool; and (c) selecting a final list of higher risk species based on tool results.

The first step of the prescreening analysis was to identify potential Arctic AIS. One hundred species were identified using a combination of data sources including published articles, grey literature and web‐based global invasive species databases (see details in Figure S1.1). Species' biological and ecological characteristics for survival in Arctic conditions (e.g. temperature and salinity tolerances) and their ability to arrive in the region via shipping were evaluated, identifying a total of 31 potential Arctic AIS (Figure S1.2a). These species were then ranked to evaluate invasion risk using the Canadian Marine Invasive Species Tool (CMIST; Drolet et al., 2016)—a rapid screening‐level risk assessment tool to quantify the risk of existing or potential marine invaders in a given area. The semiquantitative tool uses existing information and expert opinion to evaluate the potential for arrival, establishment, spread and impact by a given species and has been applied in a number of eco‐regions within Canada (DFO, 2017; Drolet et al., 2016; Moore, Lowen, & DiBacco, 2018; Therriault et al., 2018). CMIST scores are computed based on responses to 17 questions related to the likelihood and impact of invasion (see Figure S1.2a) according to the information found in the literature and other sources. Expert assessor knowledge on the risk assessment areas was used to score the potential impact of a given species based on its known impacts observed elsewhere and the availability of suitable habitats and environmental conditions. Uncertainty is scored for each risk question by assigning a qualitative score reflecting the quality of information available to answer each CMIST question. A Monte Carlo randomization procedure is then used to obtain adjusted risk scores that include uncertainty (Drolet et al., 2016).

CMIST‐ranked species with either medium or high mean risk scores (N = 18; Table 1) were retained for more detailed assessment using SDM. To this list, five potentially harmful phytoplankton species were added (Table 1; Figure S1.2b). These dinoflagellate species have been found in ballast water tanks and/or in ballast water exchange zones of Canadian domestic ships that discharge their ballast water into Canadian Arctic ports (Laget, 2017). All are known to have the capacity to produce toxins and have been implicated in harmful algal events throughout the world (Harmful Algal Information System, from the Intergovernmental Oceanographic Commission of UNESCO, http://haedat.iode.org/).

TABLE 1

List of modelled species and the methodology for species selection. Ecological groups were classified as zoobenthos, phytobenthos, zooplankton and phytoplankton. Selection methods were the Canadian Marine Invasive Screening Tool (CMIST; Drolet et al., 2016), and harmful dinoflagellate species found in ballast of vessels discharging in the Canadian Arctic (Laget, 2017)

Species	Common name	Taxa	Ecological group	Selection method	Predictors included in SDM
Amphibalanus eburneus	Ivory barnacle	Crustacea	Zoobenthos	CMIST	Bottom temperature Sea surface temperature Bottom salinity Sea surface salinity Ice thickness Depth Distance to land
Botrylloides violaceus	Violet tunicate	Tunicata	Zoobenthos	CMIST
Botryllus schlosseri	Golden star tunicate	Tunicata	Zoobenthos	CMIST
Carcinus maenas	Green crab	Crustacea	Zoobenthos	CMIST
Chionoecetes opilio	Snow crab	Crustacea	Zoobenthos	CMIST
Ciona intestinalis	Vase tunicate	Tunicata	Zoobenthos	CMIST
Littorina littorea	Common periwinkle	Mollusca	Zoobenthos	CMIST
Membranipora membranacea	Coffin box bryozoan	Bryozoa	Zoobenthos	CMIST
Molgula manhattensis	Sea grape	Tunicata	Zoobenthos	CMIST
Mya arenaria	Soft shell clam	Mollusca	Zoobenthos	CMIST
Paralithodes camtschaticus	Red king crab	Crustacea	Zoobenthos	CMIST
Codium fragile spp. fragile	Dead man's fingers	Chlorophyta	Phytobenthos	CMIST	Bottom temperature Sea surface temperature Bottom salinity Sea surface salinity Ice thickness Depth Distance to land Photosynthetically active radiation Dissolved oxygen pH Minerals and nutrients (calcite, iron, nitrate, phosphate and silicate)
Dumontia contorta	Dumont's tubular weed	Rhodophyta	Phytobenthos	CMIST
Sargassum muticum	Japanese wireweed	Phaeophycea	Phytobenthos	CMIST
Undaria pinnatifida	Wakame	Phaeophycea	Phytobenthos	CMIST
Acartia (Acanthacartia) tonsa	No common name found	Copepoda	Zooplankton	CMIST	Sea surface temperature Bottom salinity Sea surface salinity Ice thickness Depth Distance to land Chlorophyll concentration Dissolved oxygen pH
Aurelia limbata	Brown‐branded moon jelly	Cnidaria	Zooplankton	CMIST
Mnemiopsis leidyi	Warty comb jelly	Ctenophora	Zooplankton	CMIST
Alexandrium tamarense	No common name found	Dinoflagellata	Phytoplankton	Laget (2017)	Sea surface temperature Sea surface salinity Ice thickness Distance to land Photosynthetically active radiation Dissolved oxygen pH Minerals and nutrients (iron, nitrate, phosphate and silicate)
Dinophysis caudata	No common name found	Dinoflagellata	Phytoplankton	Laget (2017)
Dinophysis dens	No common name found	Dinoflagellata	Phytoplankton	Laget (2017)
Gonyaulax polygramma	No common name found	Dinoflagellata	Phytoplankton	Laget (2017)
Kryptoperidinium triquetrum	No common name found	Dinoflagellata	Phytoplankton	Laget (2017)

The final species list for modelling thus included a total of 23 known marine AIS or harmful algal species from four ecological groups (zoobenthos, phytobenthos, zooplankton and phytoplankton) that pose potential risks for invasion in the Canadian Arctic. To simplify terminology, hereafter when referring to this suite of species they will be termed as ‘AIS’. (See Table S1 for information on species' characteristics and impacts). Likewise, the term ‘invasion’ is used to make reference to the complete process of a species transitioning all invasion stages (transport, arrival, establishment and spread; Lockwood, Hoopes, & Marchetti, 2007).

Species data

Global scale occurrence data of marine invaders selected for modelling was compiled for both native and invaded ranges using global biodiversity databases such as the Global Biodiversity Information Facility (GBIF—www.gbif.org), Ocean Biogeographic Information System (OBIS—www.obis.org), and invasive species lists with coordinate location information and specific literature (Table S1). An effort was made to maximize the number and quality of occurrence records to best predict potential distributions by doing a vast and complete search of occurrence records and by checking information to find the original source of those records (García‐Roselló et al., 2015; Guisan, Graham, Elith, & Huettmann, 2007; Guisan, Zimmermann, et al., 2007; Lobo, 2008). Both native and invaded ranges were used for training and evaluating all models (Verbruggen et al., 2013), since invaded areas provide valuable information on species' tolerances to climatic conditions that may not be present in their native range (Marcelino & Verbruggen, 2015) and may improve predictions for the extent of suitable habitat (Broennimann & Guisan, 2008). A single presence record was counted per grid cell to decrease the possibility of overprediction, hence, occurrence points were considered at the same resolution as the corresponding environmental layers (García‐Roselló et al., 2015). All points were verified to ensure that they were in sea grids (and not over land). When necessary, points were moved to the closest sea grid using the Near Proximity tool in ArcMap v10.2.2.

Environmental data

Marine data layers, prepared specifically for ecological modelling, were used as environmental predictors and downloaded at a global scale from Bio‐ORACLE v2 (http://www.bio‐oracle.org/). These layers were produced with climate data describing monthly averages for the period from 2000 to 2014 representing recent (hereafter referred to as present) conditions (Assis et al., 2018). They were obtained from preprocessed global ocean re‐analyses combining satellite and in situ observations at regular two‐ and three‐dimensional spatial grids (Assis et al., 2018). A set of 37 environmental layers was used from this source comprising bottom and sea surface temperature and salinity, ice thickness, chlorophyll, dissolved oxygen, pH, photosynthetically active radiation and minerals and nutrients (iron, calcite, nitrate, phosphate and silicate; Table 1; Table S2). Long‐term maximum, minimum and mean values were used for most predictors when available (Table S2). Resolution of environmental layers was 5 arcmin (approximate 9.2 km at the equator). Bathymetry and land distance layers were obtained from Aquamaps (http://www.aquamaps.org/; Kaschner et al., 2016), but with a resolution of 30 arcmin (approximately 55 km at the equator; Table S2). Bottom substratum type was not included in the analysis due to the lack of availability of a high‐resolution global scale database; suitable benthic habitat was assumed to be present within each study region.

Habitat suitability model

Habitat suitability for selected AIS was modelled with MaxEnt v3.3.3k (Phillips, Anderson, & Schapire, 2006) and biomod2 v3.4.6 (Thuiller, Georges, Engler, & Breiner, 2020) within R v3.6.3 (R Core Team, 2019). Multiple modelling techniques were used so that the results from each approach could be combined into an ensemble model, which has been shown to reduce the biases that single models may have (Araújo & New, 2007).

MaxEnt is a machine learning method based on maximum entropy that predicts the potential geographic distribution of suitable habitat for species using species occurrence data and various combinations of environmental data layers as input. The method is one of the most widely used SDM algorithms due to its high predictive accuracy and efficiency in modelling range shifts under future climate change scenarios (Bucklin et al., 2015; Elith et al., 2006, 2011; Hijmans & Graham, 2006; Pearson, 2007). It has also recently been shown to outperform other modelling techniques in accurately predicting invasive species distributions (Battini, Farías, Giachetti, Schwindt, & Bortolus, 2019). This was complemented with SDM using a suite of other techniques within biomod2, another well‐known and widely used SDM tool (Hao, Elith, Guillera‐Arroita, & Lahoz‐Monfort, 2019) that was selected because it allows multiple models to be run on the same training and testing datasets to better ensure the consistency and accuracy of resultant ensemble models. The four models were run with biomod2: generalized linear model (GLM), random forest (RF), artificial neural network and generalized additive model (GAM) because each represents a fundamentally different modelling technique. Initially, BIOCLIM (a presence‐only model) was also included in runs but was later removed as its performance was very low, in agreement with other studies (Elith et al., 2006), and is thus the least frequently used SDM approach within biomod2 (Hao et al., 2019). More detailed information about each of these models included from biomod2 may be found in Thuiller, Lafourcade, Engler, and Araújo (2009).

MaxEnt generates background data to compare with known presence points. This study used the default option, which generates 10,000 random background points. In contrast, selected models from biomod2 require presence and absence data, hence it was necessary to create pseudo‐absence values before running the biomod2 models (Thuiller et al., 2020). For consistency with MaxEnt models, 10,000 pseudo‐absence points were selected randomly using the biomod2 default settings.

Model predictive power for MaxEnt and biomod2 was evaluated using cross‐validation with 70% of the occurrence points chosen randomly and used to train the model and the other 30% to test the complete set of models (Araújo, Pearson, Thuiller, & Erhard, 2005). For MaxEnt, data were partitioned by a random process of k = 500 training and validation iterations (Hijmans, 2012), and for biomod2 the data were partitioned randomly. Training points for MaxEnt were selected by random seeding with the convergence threshold set at 0.00001. The hinge feature in MaxEnt was used as it produces complex yet smooth and ecologically meaningful response curves and improves model performance (Merow, Smith, & Silander, 2013; Phillips & Dudík, 2008).

Environmental predictors for each species were selected based on initial models run with MaxEnt, which has been shown to have good discrimination accuracy regarding the identification of important variables compared to other models (Smith & Santos, 2019). Selected predictors included those that are most typically considered as important for related taxa and have been used in other modelling studies (Table 1; Barnes, 1999; Belanger et al., 2012; Chust et al., 2016; Cusson, Archambault, & Aitken, 2007; Gallardo et al., 2015; Jensen, Mousing, & Richardson, 2017; Leidenberger, Obst, et al., 2015; Wagner, 1977), or were identified by taxonomic experts as being important (G. Winkler, personal communication, 2017; A. Rochon, personal communication, 2018). After the first run, predictors with a relative contribution score of <4% were excluded (Jueterbock et al., 2016). Correlations may lead models to produce erroneous response curves to predictors that do not reflect species physiological tolerances (Marcelino & Verbruggen, 2015). Thus, highly correlated environmental predictors (correlation coefficient ≥ 0.7; Dormann et al., 2013), as determined using the SDMtoolbox (Brown, 2014; Table S2), were evaluated and the one with the highest contribution to the modelling exercise retained for model construction. As per Goldsmit et al. (2018), environmental predictors were identified by evaluating the combination of: (a) the response curves for each species—to evaluate if the predictor behaves in a biologically meaningful way in the model (Marcelino & Verbruggen, 2015); (b) a species‐specific Jackknife test—to evaluate the contributions of the various parameters and analyse the importance of predictors; and (c) the estimates of the contribution of each environmental layer to model prediction. A minimum of three environmental predictors were included in each model. Each species was then modelled again with the selected environmental predictors (Table S3). Model performance was evaluated as the area under the curve (AUC) of the receiver operating curve and true skill statistic (TSS). AUC is the probability that the model correctly ranks a random presence site versus a random site from the study area (Phillips et al., 2009). Values close to 1 indicate good prediction in site discrimination, a value of 0.5 indicates a prediction equal to random and values lower than 0.5 indicate a performance that is worse than random predictions. TSS assesses the accuracy of predictions using sensitivity (proportion of correctly predicted presences) and specificity (proportion of correctly predicted absences) in its calculation (TSS = sensitivity + specificity − 1) and is an appropriate evaluation alternative for model predictions converted to binary maps using a threshold (Allouche, Tsoar, & Kadmon, 2006). It ranges from −1 to +1, where values between 0 and −1 indicate performance no better than random, while a statistically reliable model performance is indicated by values >0.4, excellent models by a minimum of 0.7, and 1 indicates perfect agreement with the model (Allouche et al., 2006). Models with a TSS score of <0.7 and sensitivity = 0 were excluded from the final ensemble predictions.

MaxEnt and biomod2 models were each run five times. The continuous values produced by the models were transformed to binary values to identify predicted suitable and unsuitable habitat since continuous model projections may be difficult to interpret. In addition, binary presence/absence maps are more useful for risk assessment exercises and are required for some model evaluations (Liu, White, & Newell, 2013). For each species, the binary results for the five runs per model were averaged to create a single binary result per model. The resulting binary results for each model were then averaged to produce the final ensemble binary projection for each species. Binary transformation was done using the maximum training sensitivity plus specificity threshold, which maximizes TSS values to create binary maps and has been shown to produce the most accurate predictions (Andrade, Velazco, & Júnior, 2020; Jiménez‐Valverde & Lobo, 2007; Liu et al., 2013), and has been used in other studies with similar objectives (e.g. Bellard et al., 2013; Duffy et al., 2017; Wisz et al., 2015). Following transformation, all binary models were delimited using a threshold for the maximum depth each species could inhabit according to their ecological requirements (Goldsmit et al., 2018; Table S4). An exception was made for phytoplankton as they were all surface dinoflagellates; thus, there was no need to consider maximum depth for their distribution (A. Rochon, personal communication, 2018). Heat maps showing the total number of modelled AIS that may find suitable habitat in a region—hereafter ‘AIS richness’—were then created using combined maps representing the cumulative number of species (of the 23 modelled) predicted to find suitable habitat in a given grid cell at global and pan‐Arctic scales. It should be kept in mind that, at a global scale, richness includes potential native and invaded ranges of all species modelled, while at the pan‐Arctic scale, richness includes mainly predicted potentially suitable habitats.

Future projections

Once all models were run, tested and validated, predicted habitat suitability was evaluated under projected global change scenarios at global and pan‐Arctic scales. The same set of environmental layers that were used for contemporary models were used for projected future models although only temperature, salinity and ice layers were available for projected future scenarios (Table S3). This is common practice in similar modelling exercises (e.g. Goldsmit et al., 2018; Jueterbock et al., 2016; Leidenberger, De Giovanni, Kulawik, Williams, & Bourlat, 2015; Weinert et al., 2016). Future environmental layers were obtained from Bio‐ORACLE v2 for RCP4.5 emission scenario for the years 2050 and 2100. The average and range of climatic anomalies for temperature, salinity and sea ice thickness of this scenario at the end of the century are given in Assis et al. (2018). In short, it represents an intermediate greenhouse emission (temperature anomaly of 2.4°C by 2100; Moss et al., 2008) and stabilization scenario resulting from the implementation of appropriate mitigation strategies (Clarke et al., 2007). This anomaly also coincides approximately with the expected increase in temperature in Arctic regions under the proposed efforts of the 2015 Paris Agreement (Solecki et al., 2018). Future layers were produced for 2040–2050 and 2090–2100 by averaging data from distinct atmosphere–ocean general circulation models provided by the Coupled Model Intercomparison Project 5, which was developed by the World Climate Research Programme's Working Group on Coupled Modelling (see Assis et al., 2018 for more information).

Resulting distribution models of each individual species were compiled by ecological groups (zoobenthos, phytobenthos, zooplankton and phytoplankton) and these combined models were used to compare present with forecasted distributions for the years 2050 and 2100. This was achieved by spatially analysing species richness and distribution change between present and future climate scenarios using ArcMap v10.2.2 and SDMtoolbox (Brown, 2014). Overlap of regions was analysed at both global and pan‐Arctic scales, to evaluate predicted latitudinal shifts in suitable habitat. Changes in the distribution of suitable habitat over time were evaluated independent of natural or anthropogenic habitat overlap to identify regions of loss, gain and no change in time, for both time frames.

RESULTS

Analysis of environmental predictors and model evaluation

Varying combinations of environmental predictors contributed to explaining the species distribution models of the ecological groups analysed (Figure 2). Sea surface temperature and land distance contributed to models for all ecological groups, although in varying proportions. Bottom temperature and depth were important in explaining models for zoobenthos and phytobenthos, in addition to ice thickness for the latter. Other main predictors for zooplankton included sea surface salinity, ice and depth, whereas important phytoplankton predictors included nutrients and minerals (especially iron; Figure 2; Table S3). All environmental predictors used to construct the final model for each species were within training ranges. Analysis of unimodal environmental response curves indicated that environmental conditions were within suitable ranges for all modelled species.

FIGURE 2

Proportion of times individual predictors were used for model building per ecological group. Bot Sal, bottom salinity; Bot Temp, bottom temperature; DissO2, dissolved oxygen; LandDist, land distance; SSS, sea surface salinity; SST, sea surface temperature

After averaging the five replicates of each model, most AUC values were >0.9 and TSS > 0.8, indicating good model performance and good prediction in site discrimination with most of the models included in the ensemble results. There were only four species with somewhat lower AUC or TSS values for some of the averaged models (Chionoecetes opilio, Dinophysis caudata and Gonyaulax polygramma with AUC between 0.76 and 0.83 in MaxEnt models; and G. polygramma and Aurelia limbata with TSS between 0.75 and 0.79 in GAM, GLM and Maxent models), but still well above the thresholds for random site prediction (Table S5). A few individual replicated GAM, GLM and RF models (N = 19) for five species were excluded from model‐specific binary outputs (maximum 3/species) due to their having TSS < 0.7 and/or sensitivity = 0; which resulted in model averages based on fewer than five replicates for these species (Table S5).

Hotspots vulnerable to invasions

From the 23 modelled species, up to 20 were predicted to have overlapping suitable habitat in Arctic regions, mostly located near the coasts surrounding Hudson Bay, Northern Grand Banks/Labrador, Chukchi/Eastern Bering seas and Barents/White seas (Figure 3a). These hotspots are predicted to have potentially greater AIS richness compared to other Arctic regions in the present and both modelled future scenarios. The areal extent of hotspots vulnerable to invasion is predicted to increase over time; this increase is predicted not only in total area, but also in the number of species projected to encounter appropriate habitat there (Figure 3b).

FIGURE 3

Predicted total aquatic invasive species (AIS) richness at an Arctic scale: (a) predicted hotspots of AIS richness for present and future (2050 and 2100) conditions in the Arctic (1, Hudson Bay; 2, Northern Grand Banks/Labrador; 3, Chukchi/Eastern Bering seas; 4, Barents and White seas). Colours represent the number of overlapping species with predicted suitable habitat in a given area; (b) predicted future extension of suitable habitat by area for each category (natural breaks Jenks) of AIS richness. Values shown in the bars are the net percentages of suitable habitat at each level of AIS richness within each projected time period

These same regions were predicted to be hotspots for individual ecological groups, particularly zoobenthos and phytoplankton (Figure 4a,b,d). When evaluated independently, these groups showed the same pattern of future increased areal extent in the predicted suitable habitat, but at varying taxa‐dependent scales (absolute area extent: zoobenthos ~6 × 106 km2, phytobenthos, zooplankton and phytoplankton ~1–2 × 106 km2; Figure 4a–d). Despite these differences in magnitude for absolute predicted future suitable habitat, at a relative scale the percentage change in suitable habitat was predicted to be greater through time for some groups (e.g. zoobenthos). Furthermore, the relative change in predicted suitable habitat for various categories of AIS richness differed between ecological groups (Figure 4e–h). For example, zoobenthos were predicted to have greater future increases in the areal extent of habitats suitable for a high number of overlapping species (i.e. with high AIS richness).

FIGURE 4

(a–d) Predicted hotspots of aquatic invasive species (AIS) richness per ecological group under present and future (2050 and 2100) conditions at an Arctic scale. Colours represent the number of overlapping species with suitable habitat in a given area. (e–h) Predicted future change of suitable habitat for each category of AIS richness

Only one species of zooplankton (Aurelia limbata) was predicted to have extended suitable habitat in the Arctic; hence, no potential hotspots with overlapping species were observed (Figure 4c). It should be highlighted that a few zoobenthic species have native ranges included in the area of analysis: (a) C. opilio: Beaufort Sea, Bering Sea and Northern Grandbanks/Labrador; (b) Littorina littorea: White Sea; (c) Mya arenaria: Northern Grandbanks/Labrador; and (d) Paralithodes camtschaticus: Bering Sea (Figure 1; Table S1).

At a global scale, total species richness was higher and largely concentrated on the coasts of Northern Europe and the White Sea, Black Sea, Northwest Atlantic, some regions of the Northwest and Northeast Pacific, Iceland and southern regions of Temperate South America and Temperate Australasia (Figure S2).

Distribution change

Suitable habitat was projected to increase in the future for all AIS combined as well as for all individual ecological groups in the Arctic realm under both climate change scenarios (Figures 5 and 6), with zoobenthos showing the greatest potential habitat gains (more than double the gains of other ecological groups; Figures 5a and 7) under future climate scenarios. The overall predicted future suitable habitat changes at the pan‐Arctic scale for all AIS combined were net habitat gains of +68.8% and +121.4% by 2050 and 2100, respectively (Figures 6 and 7). However, the same analysis predicted an overall suitable habitat loss at a global scale in both future scenarios (−4.3% by 2050 and −3% by 2100; Figures 6 and 7).

FIGURE 5

Overall change in predicted suitable habitat for each ecological group at an Arctic scale between present and future scenarios (2050 and 2100): (a) zoobenthos, (b) phytobenthos, (c) zooplankton and (d) phytoplankton

FIGURE 6

Overall change in predicted suitable habitat of all aquatic invasive species combined at global and pan‐Arctic scales between present and future scenarios (2050 and 2100)

FIGURE 7

Percentage of predicted habitat suitability change (gain or loss) for the four ecological groups and all aquatic invasive species together at global and pan‐Arctic scales. PB, phytobenthos; PP, phytoplankton; ZB, zoobenthos; ZP, zooplankton

These results show that most modelled species are predicted to find suitable habitats in colder regions, with a trend towards poleward shifts in future distributions, particularly in northern latitudes (the poleward shift to the south is predicted to be lower; Figure 6). At a global scale, however, planktonic species are predicted to experience habitat loss (except for a small habitat gain predicted for phytoplankton towards the year 2100), while benthic species are predicted to have a positive gain, although much less than that at the pan‐Arctic scale (Figure 7; Figure S3).

DISCUSSION

This study predicted the distribution of suitable habitat for 23 high‐risk AIS at pan‐Arctic and global scales and identified hotspots of suitable habitats for multiple species. Overall, results indicate that suitability will increase over time in Arctic regions, particularly in Hudson Bay, Northern Grand Banks/Labrador, Chukchi/Eastern Bering seas and Barents/White seas regions. Gradually, these regions could become more vulnerable to invasions due to increases in the extent of suitable habitat and in the number of species that gain access to these habitats and establish reproducing populations (i.e. AIS richness). This trend was observed across all assessed ecological groups that showed gains in suitable habitat under future scenarios, with zoobenthos exhibiting the greatest future distributional changes (more than double). Nevertheless, all other ecological groups are predicted to increase their future suitable habitat by 30% up to more than 90%.

Interestingly, predicted shifts under future climate scenarios differed at pan‐Arctic and global scales. When all species were considered collectively, there was a great predicted increase in suitable habitat at the pan‐Arctic scale. In contrast, the net result at a global scale was a predicted loss of suitable habitat, suggesting a marked northward shift with global change. Projected future scenarios considered only the direct effect of climate change. However, projected increases of shipping to the area will likely exacerbate the risk as new transport routes open due to longer ice‐free conditions that increase opportunities for natural resource extraction (Melia et al., 2016; Miller & Ruiz, 2014; Smith & Stephenson, 2013). Likewise, the Antarctic could be affected by native and invasive species from northern regions, given that environmental conditions are similar. Arctic ports could act as possible sources of AIS to the Antarctic since some ships transit between poles and could act as a potential pathway of transportation, although survival and establishment in such scenarios are unknown (McCarthy et al., 2019).

The use of CMIST as a screening tool to identify high‐risk species for specific ecoregions is a recent approach to stay ahead of the invasion process (Drolet et al., 2016; Moore et al., 2018; Therriault et al., 2018). Although the use of CMIST combined with habitat modelling, as done in this study, places more emphasis on the likelihood of invasion than the likelihood of impact component of risk, such a tool may be used as a first step to develop a ranked watch list to guide early detection and monitoring efforts and to prioritize species for detailed risk assessments and potential regulations (Locke, Mandrak, & Therriault, 2011). CMIST provided the opportunity to identify potential AIS of greater risk to the Arctic by considering known species' characteristics and impacts in other regions of the world where they have already invaded. Eighty‐seven per cent of the species included in the present modelling assessment belong to groups with the greatest numbers of known invasions in the marine Arctic ecosystem (Chan et al., 2019). The uneven emphasis of organisms retained for modelling from the different ecological groups (maximum of 11 zoobenthic species versus a minimum of 3 zooplankton species) may, in part, be explained by the fact that most known introduced marine organisms are benthic species (Streftaris, Zenetos, & Papathanassiou, 2005). These proportions are comparable with previous assessment studies (Leidenberger, Obst, et al., 2015), although the uneven coverage does make general trends more difficult to interpret.

Sea surface temperature and land distance were retained in all models, consistent with other studies that found these predictors to contribute significantly to explaining species distributions (Bradie & Leung, 2017; Leidenberger, Obst, et al., 2015; Stelzer, Heyer, Bourlat, & Obst, 2013). Water temperature has been identified as the most relevant predictor of global marine species distribution and land distance of moderate importance for both benthic and planktonic species (Bosch, Tyberghein, Deneudt, Hernandez, & De Clerck, 2018; Bradie & Leung, 2017). Depth was an important predictor, as shown by Bosch et al. (2018) and Snickars et al. (2014) who suggested that bathymetry is of high relevance for modelling various taxa. Ice thickness was moderately important for phytobenthos and zooplankton, perhaps reflecting limits on distribution of the former due to ice abrasion, changes in light exposure and a preference of ice‐free waters for the latter (Clark et al., 2013; Kube, Postel, Honnef, & Augustin, 2007; Pascual et al., 2015; Richardson, 1979). Iron was an important predictor for all modelled phytoplankton species, likely as it is known to be important to phytoplankton growth, abundance, dominance and species distributions (Hecky & Kilham, 1988; Moore et al., 2013) and plays a role in the development of harmful algal bloom species (Doucette & Harrison, 1990; Wells, Mayer, & Guillard, 1991). In Arctic regions, meltwater can be a significant bioavailable source of iron to surrounding coastal oceans (Bhatia et al., 2013; Tovar‐Sánchez et al., 2010) and evidence suggests a link between ice and iron from glacial meltwater leading to blooms in some phytoplankton taxa (Aguilar‐Islas, Rember, Mordy, & Wu, 2008; Joli et al., 2018). This is important given that increased iron due to a shift from an ice‐covered to an open water Arctic Ocean (Screen, 2018; Seneviratne et al., 2018) could create favourable conditions for harmful species that may arrive in the region.

Four Arctic regions (i.e. Hudson Bay, Northern Grand Banks/Labrador, Chukchi/Eastern Bering seas and Barents/White seas) were identified as potential hotspots for invasions. Invasion hotspots could pose even greater risks if they coincide with major shipping routes, biodiversity hotspots or areas of special interest/importance, which is the case for most of the invasion hotspots predicted in the present study. The predicted invasion hotspots overlap with regions that have been highlighted for being special (with regard to their uniqueness, importance for species life histories, threatened species and/or habitats, biological productivity, diversity, etc.). For example, some Ecologically and Biologically Significant Areas (also known as EBSAs) and Marine Refuges coincide with the predicted invasion hotspots in Northern Grand Banks/Labrador, Chukchi/Eastern Bering seas and Barents/White seas (Kenchington et al., 2011; Speer & Laughlin, 2011; Templeman, 2007; www.dfo‐mpo.gc.ca/oceans/oeabcm‐amcepz/refuges/index‐eng.html), as well as in the southern Hudson Bay, which has been identified as being an area of high biological importance (Stephenson & Hartwig, 2010). Ecoregions with a high concentration of polynyas could also be at greater risk, given that they may act as important biological hotspots due to their increased productivity and biodiversity (Marchese, 2015). This may be the case for the Chukchi and Eastern Bering seas, which are known to have concentrations of polynyas (CAFF, 2013). This is consistent with previous studies that have identified biodiversity hotspots in marine and terrestrial ecosystems as being at particular risk to future invasions (Bellard et al., 2014; Li, Liu, Kraus, Tingley, & Li, 2016) and climate change (Ramírez, Afán, Davis, & Chiaradia, 2017). The presence of AIS in biological hotspots could endanger native and endemic species as AIS are considered to be one of the leading causes of biodiversity loss and contemporary extinctions (Bellard, Cassey, et al., 2016; Blackburn et al., 2019) with impacts on ecosystem structure and functioning, including changes in food webs, biomass, flux rates, etc. (Ehrenfeld, 2010). Relating where invasion hotspots overlap with biological hotspots could be crucial for prioritizing conservation efforts.

The Barents Sea has been highlighted as a region that already has a substantial number of invasions and environmental conditions that increase the probability of successful establishment (Chan et al., 2019). Two of the species modelled in the present study, C. opilio and P. camtschaticus, are already established there (Chan et al., 2019; Dvoretsky & Dvoretsky, 2009; Hansen, 2015). The former is thought to have been introduced by ballast water, while the latter was an intentional introduction (Alvsvåg, Agnalt, & Jørstad, 2009; Jørgensen & Nilssen, 2011). The Barents Sea appears to be in transition from a cold Arctic to warm Atlantic climate regime (Lind, Ingvaldsen, & Furevik, 2018), making it particularly vulnerable to invasion. Indeed, it is predicted to suffer one of the largest future habitat losses by endemic species in the Arctic (Renaud et al., 2019), leaving potential niches available for novel species to occupy. Additionally, Arctic ecoregions will be more exposed to potential future arrivals with further ice reduction and increased navigability of the Northern Sea Route and the Northwest Passage, although much greater investment in infrastructure, navigation and communications would be needed to this end (Buixadé Farré et al., 2014). Nonetheless, both human activities and AIS are likely to increase over time in the Arctic, as has been observed over the last few years (Chan et al., 2019; Dawson, Pizzolato, Howell, Copland, & Johnston, 2018).

At a global scale, planktonic organisms generally showed loss of suitable habitat that, in all but one case, exceeded predicted gains. This may be due to their being dispersive pelagic organisms which have the capability of expanding rapidly and show extensive distribution changes in response to temperature increases related to global warming (Poloczanska et al., 2013). Predicted future changes in sea surface temperatures, in particular, which are expected to outpace shifts in bottom temperature (Levitus et al., 2012), may be driving this pattern. Indeed, models in the present study showed that sea surface temperature was more important for planktonic organisms, whereas bottom temperature was important for benthic organisms. However, the pattern of predicted change by taxa is different at the pan‐Arctic scale where all ecological groups showed high positive suitable habitat gains coinciding with other modelling studies (de Rivera et al., 2011; Goldsmit et al., 2018; Jueterbock et al., 2013; Townhill et al., 2018; Villarino et al., 2015; Ware et al., 2016). This could be explained by the fact that temperature increases have been shown to be greater in the Arctic than in other areas (Bush & Lemmen, 2019; Pendleton et al., 2019) and thus predicted gains in suitable habitat there would be expected to be relatively greater than losses elsewhere.

A marked poleward shift was predicted, consistent with other forecasting studies for various types of marine organisms, including invertebrates, algae and fish (Byrne et al., 2016; Chust et al., 2013; Goldsmit et al., 2018; Jones & Cheung, 2015; Mackey et al., 2012; Townhill et al., 2018; Valle et al., 2014; Wisz et al., 2015). Poleward shifts may serve as an early warning signal of ecosystem change due to climate warming. Multiple studies across different taxa are showing consistency between observed responses of marine organisms to climate change, particularly in high‐latitude regions (Poloczanska et al., 2013). Cold‐adapted species typically have narrow thermal windows and low‐energy demand lifestyles, making them more sensitive to temperature changes (Poloczanska et al., 2016). Predicted poleward shifts in suitable habitat in the present study are disproportionately located in Arctic areas, as has been shown in Cheung et al. (2009), García Molinos et al. (2015) and Jones and Cheung (2015). Invasive species in the Antarctic region have been limited by physiology at cold temperatures rather than geographic limits, but global warming may remove those physiological barriers and alter distributions (Aronson et al., 2007). However, it should be noted that a substantial number of the species modelled here are already distributed in northern/cold regions, which may bias observed patterns. Although, three modelled species (C. intestinalis, B. schlosseri and U. pinnatifida) are among the most likely to become invasive in the Antarctic Peninsula according to a recent horizon scanning study (Hughes et al., 2020).

The locations with greater reductions of the sea ice season over the last 40 years (Stammerjohn, Massom, Rind, & Martinson, 2012) coincide with the regions of predicted invasion hotspots in the present study, suggesting that these areas are already experiencing changing environmental conditions and that such variations are projected to continue in time, including the probability of a complete ice‐free Arctic during summer (Screen, 2018; Sigmond et al., 2018). The combination of species invasions and predicted reductions in sea‐ice cover could alter reproduction/phenology timing, energy pathways and food‐web dynamics with subsequent impacts on production at higher trophic levels (Haug et al., 2017; Poloczanska et al., 2013 and references therein). In this context, AIS could take advantage of new habitats and resources and outcompete native species, which are generally expected to be more sensitive to temperature changes given that they live within a narrow low‐temperature range (Peck, 2005). Invasive, non‐native and native boreal species could thus expand their ranges as suitable habitats in polar regions become more common (de Rivera et al., 2011; Goldsmit et al., 2018; Renaud et al., 2019; Ware et al., 2016) with subsequent impacts on richness, community composition and Arctic ecosystems. However, it remains unknown whether biogeographic boundary locations will change (Costello et al., 2017).

The combination of methods used in this study permitted the inclusion of different modelling algorithms for more robust prediction results based on an ensemble model. Models had varying data requirements (e.g. presence‐background or presence–absence) and used varying statistical approaches (e.g. statistical regression models and machine‐learning decision trees). Modelling studies have restrictions and limits that should be recognized when interpreting results. A main limitation of the present study is that potential species occurrences and range changes are based solely on abiotic factors and do not reflect the realized distribution of a species, assuming that species are in equilibrium with their environment (Bellard et al., 2018; Guisan & Zimmermann, 2000). Biotic interactions, potential resources and life history traits were not included in the modelling exercise, even though they are known to be important in shaping species spatial distributions (Wisz et al., 2013). Moreover, including physiological data such as reproductive temperature to represent phenology, can significantly change SDM predictions. Including this type of data has been shown to strongly limit predicted northward shifts under future climate change scenarios (Chefaoui, Serebryakova, Engelen, Viard, & Serrão, 2019). However, this particular aspect was somewhat offset by the fact that environmental conditions for reproduction, together with conditions needed for different life stages, were considered in this study during the species selection using CMIST scores. Given that the present study modelled the potential distribution of suitable habitat in a new environment, biotic interactions may not be realistic if applied in a different context, such as different regions or time periods (Anderson, 2017 and references therein). Despite limitations and restrictions, SDM can provide valuable information to help manage resources in marine ecosystems that will face increasing anthropogenic pressures (Reiss et al., 2014) and has been shown to be a useful predictive tool to assess various taxa concurrently (Gallardo et al., 2015; Leidenberger, Obst, et al., 2015). SDM may also be particularly useful for regions such as the Arctic, where predicting biodiversity changes under global warming effects is challenging due to the paucity of baseline data for most organisms (Renaud et al., 2019; Wassmann, Duarte, Agusti, & Sejr, 2011).

Distributional change of AIS may increase risk in previously unimpacted areas, potentially creating new problems for wildlife and even human health (such as harmful algal blooms). These types of episodes are being facilitated by climate change (Hallegraeff, 2010; Kordas, Harley, & O'Connor, 2011; Poloczanska et al., 2013, 2016). The present study predicts changes over such an extensive area (~6 × 106 km2, equivalent to almost three times the size of Greenland) that altered community structure may be widespread. Results presented here provide information on AIS that pose some of the greatest threats to the Arctic, with areas at greatest risk identified as hotspots. This information is valuable given that aquatic invaders are understudied in the context of climate change (Bellard et al., 2018). The Arctic is predicted to be affected by increased habitat suitability for a number of potential AIS and, given its vast area, could be severely impacted by AIS accumulating in specific locations under both current and future environmental conditions. Identification of hotspots through SDM, predictions of habitat vulnerability and particular areas of concern could guide ballast water practices and other management actions, including prevention, early detection monitoring, rapid response and conservation planning. Information such as that provided by the present study should help guide how to prioritize management efforts in this unique and vast region.

Fig S1

Click here for additional data file.

Fig S2

Click here for additional data file.

Fig S3

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.

Table S3

Click here for additional data file.

Table S4

Click here for additional data file.

Table S5

Click here for additional data file.

Table S5

Click here for additional data file.