fmars-06-00538 August 30, 2019 Time: 17:21 # 1ORIGINAL RESEARCHpublished: 03 September 2019doi: 10.3389/fmars.2019.00538Edited by:Elizabeth Fulton,Commonwealth Scientificand Industrial Research Organisation(CSIRO), AustraliaReviewed by:Yong Jiang,Ocean University of China, ChinaKristina Øie Kvile,University of Oslo, Norway*Correspondence:Paul E. Renaudpr@akvaplan.niva.noSpecialty section:This article was submitted toGlobal Change and the Future Ocean,a section of the journalFrontiers in Marine ScienceReceived: 30 March 2019Accepted: 15 August 2019Published: 03 September 2019Citation:Renaud PE, Wallhead P, Kotta J,Włodarska-Kowalczuk M,Bellerby RGJ, Rätsep M, Slagstad Dand Kuklin´ski P (2019) ArcticSensitivity? Suitable Habitatfor Benthic Taxa Is SurprisinglyRobust to Climate Change.Front. Mar. Sci. 6:538.doi: 10.3389/fmars.2019.00538Arctic Sensitivity? Suitable Habitatfor Benthic Taxa Is SurprisinglyRobust to Climate ChangePaul E. Renaud1,2* , Phil Wallhead3, Jonne Kotta4, Maria Włodarska-Kowalczuk5,Richard G. J. Bellerby3,6, Merli Rätsep4, Dag Slagstad7 and Piotr Kuklin´ski51 Akvaplan-niva, Tromsø, Norway, 2 Department of Arctic Biology, The University Centre in Svalbard, Longyearbyen, Norway,3 Norwegian Institute for Water Research, Bergen, Norway, 4 Estonian Marine Institute, University of Tartu, Tartu, Estonia,5 Institute of Oceanology, Polish Academy of Sciences, Sopot, Poland, 6 SKLEC-NIVA Centre for Marine and CoastalResearch, State Key Laboratory of Estuarine and Coastal Research, East China Normal University, Shanghai, China,7 SINTEF, Trondheim, NorwayArctic marine ecosystems are often assumed to be highly vulnerable to ongoing climatechange, and are expected to undergo significant shifts in structure and function.Community shifts in benthic fauna are likely to result from changes in key physico-chemical drivers, such as ocean warming, but there is little ecological data on mostArctic species to support any specific predictions as to how vulnerable they are, orhow future communities may be structured. We used a species distribution modelingapproach (MaxEnt) to project changes over the 21st century in suitable habitat areafor different species of benthic fauna by combining presence observations from theOBIS database with environmental data from a coupled climate-ocean model (SINMOD).Projected mean % habitat losses over taxonomic groups were small (0–11%), andno significant differences were found between Arctic, boreal, or Arcto-boreal groups,or between calcifying and non-calcifying groups. However, suitable habitat areas for14 of 78 taxa were projected a change by over 20%, and several of these taxaare characteristic and/or habitat-forming fauna on some Arctic shelves, suggestinga potential for significant ecosystem impacts. These results highlight the weaknessof general statements regarding vulnerability of taxa on biogeographic or presumedphysiological grounds, and suggest that more basic biological data on Arctic taxa areneeded for improved projections of ecosystem responses to climate change.Keywords: benthic invertebrates, climate warming, multiple stressors, ocean acidification, species-distributionmodelingINTRODUCTIONOngoing climate change affects multiple environmental drivers, with direct and indirectimplications for marine ecosystems. For example, reductions in sea-ice volume and extent,increasing air and sea-surface temperatures, and freshwater and biogeochemical discharge in theArctic have already been noted to be occurring at rates higher than the global average (IPCC,2014; AMAP, 2017) with consequences for regional ocean acidification rates (Bates et al., 2012;Bellerby, 2017). Anthropogenic impacts from fossil fuel CO2 emissions are linked with changes inthese drivers and have led to changes in ocean chemistry. These impacts include reductions in bothFrontiers in Marine Science | www.frontiersin.org 1 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 2Renaud et al. Species Distributions in a Changing Arctic Oceanocean pH and carbonate mineral saturation states, which havebeen implicated in decreased performance of many marinetaxa (Pörtner, 2008; Kroeker et al., 2011). Models projectthat ongoing ocean acidification (OA) will continue over thecoming decades, and the degree of OA will be dependentmainly on social and policy decisions (e.g., Bellerby et al., 2014,2018). The amplification of both warming and acidificationobserved in the Arctic, in combination with other expectedphysical and biological impacts, has led to the commonviewpoint that Arctic ecosystems are particularly vulnerable toon-going climate changes (Smetacek and Nicol, 2005; Renaudet al., 2007; Brierley and Kingsford, 2009; Doney et al., 2012;Hoppe et al., 2018). Changes in community composition thatmay result are not merely of academic interest: they mayreflect underlying functional responses with broad-reachingimplications for ecosystem processes and services.Laboratory experiments and field sampling have generatednumerous but often contradictory predictions of how changesin climate-change-related drivers will affect marine life, fromindividuals to ecosystem processes. Physiology, phenology andlife-history, community structure, and system productivity(Pörtner et al., 2005; Ardyna et al., 2014; Renaud et al., 2015,2018) are just several examples of impacts that may be felt: fromthe very basic processes of life to the provisioning of ecosystemservices. Until now, however, most studies have focused oneffects of single factors, and usually in laboratory or isolated fieldsettings (Pedersen et al., 2016; AMAP, 2018). Recent syntheseshave called for investigations of broader thematic and spatialscope, including studying multiple factors in combination forarriving at more realistic predictions (e.g., AMAP, 2018), andadopting a pan-Arctic perspective when identifying trends inecosystem response (Wassmann, 2015). These perspectives areneeded for developing mitigation, adaptation, and managementstrategies for marine ecosystems across the Arctic (Halpernet al., 2012). Whereas some Arctic taxa are clearly threatened bychanging climatic conditions (Kovacs et al., 2011; Doney et al.,2012; Grebmeier, 2012; Beaugrand and Kirby, 2018), a multi-factor investigation of vulnerability of Arctic taxa in general hasnot been performed.One of the most intuitive changes that could be expecteddue to climate warming is the poleward expansion of borealorganisms. This has already been observed for both benthicand pelagic organisms, and is suggested to be largely relatedto thermal tolerance (e.g., Wethey and Woodin, 2008; Sundayet al., 2012; Neukermans et al., 2018). Indeed, dominance ofboreal and Arctic taxa in seafloor (benthic) communities of theBarents Sea have been shown to fluctuate with temperature inthe region over the past 100+ years (Blacker, 1965; Matishovet al., 2012). Niche theory is useful for making predictions asto historical or future distributions in the context of climatechange. Single-variable models (thermal-envelope models) havebeen used to predict, for example, changes in fish distributionwith climatic warming (Cheung et al., 2009). In Arctic systems,this technique is problematic due to high uncertainty ineither fundamental or realized niches of nearly all taxa inthe ecosystem. On-line databases can provide one way aroundthis problem by identifying the characteristics of an organism’srealized niche. This technique has been used to identify benthictaxa that may expand or contract under warming scenarios(Renaud et al., 2015). Although biased by the geographicalrange of data included in databases, and non-random samplingin general, the method can also provide information formulti-dimensional realized niches. Once current distributionsare mapped, a suite of species distribution models may beapplied to project future distributions in response to multipleenvironmental drivers. Such modeling activities have beenidentified as critical next steps in understanding climate changeeffects on benthic community structure in Arctic environments(Renaud et al., 2015).A key requirement for implementing these models is a wellresolved and validated hydrographic/ocean chemistry modelwith realistic climate forcings. This allows for accurate physicalcharacteristics to be assigned to each sampling location in bothhindcasting and forecasting modes. One such model has recentlybeen expanded to include the needed spatial and temporalperspectives for the carbonate system, and has been used toupscale results of a mesocosm experiment to other Arctic areas(Bellerby et al., 2012).This study investigates the vulnerability of benthic faunafound in the Arctic and its marginal seas to end-of-centuryclimate change. We use outputs from a state-of-the-art ecosystemmodel to drive species-distribution modeling for a broad rangeof characteristic macro- and mega- benthic taxa characteristicof the Arctic’s shelf seas. We analyse results to determine ifArctic taxa are more susceptible to changes in a combinationof environmental drivers than boreal or more widespreadtaxa. We also investigate whether calcification status orhigher taxonomic classification (phylum, class) could indicatevulnerability. Implications of the presence or absence of suchpatterns are discussed.MATERIALS AND METHODSEnvironmental ProjectionsWe ran the ocean biogeochemical model SINMOD (Slagstadet al., 2011, 2015), which includes a CO2 system module(Bellerby et al., 2012), to estimate bottom water temperatureand aragonite saturation state ((ar)). The SINMOD domain ispan-Arctic (see Figure 1) with 20 km horizontal grid resolutionand 25 fixed depth levels, and calculations were based onbiweekly saved output. Saturation state describes the seawaterconcentrations of carbonate and calcium ions relative to theequilibrium concentrations for that mineral; conditions near orbelow saturation ( = 1) are generally considered unfavorablefor production or maintenance of normal calcification. SINMODwas first run to produce a hindcast simulation (1979–2008),and then to produce a projection (2001–2099) under the SRES(Special Report on Emissions Scenarios) scenario A1B (seeSlagstad et al., 2015; Wallhead et al., 2017, for more details).SRES A1B is a mid-range business-as-usual scenario assuming a“balanced” use of fossil vs. non-fossil energy sources (Nakicenovicet al., 2000) that results in end-of-century global responsesthat are comparable to the more recent RCP6.0 scenario andFrontiers in Marine Science | www.frontiersin.org 2 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 3Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 1 | Bias-corrected SINMOD hindcasts and projections for bottom water temperature (A–C) and bottom water aragonite saturation state (D–F) under SRESA1B scenario. Left and middle figures show averages over the hindcast period (1978–2008) and the projection period (2090–2099), respectively; figures on the farright show the differences. Black lines show 50, 200, and 500 m depth contours.substantially more conservative than the SRES A2 or RCP8.5scenarios (Collins et al., 2013). SINMOD bottom-water outputwas corrected using bias estimates, calculated as a function ofspatial position (temperature) or of model bathymetry ((ar))using a compilation of in situ observations and matched modeloutput (see Wallhead et al., 2017). Water depth over theSINMOD grid was estimated by interpolating high-resolutionbathymetry products (the 500 m resolution IBCAOv3, Jakobssonet al., 2012, where available, otherwise the 2 min ETOPOv2,National Geophysical Data Center, 2006).Niche Description and Climatic ChangeOccurrence records for 95 benthic taxa characteristic of Arcticshelf infaunal and epifaunal communities were extracted fromthe Ocean Biogeographic Information System database (OBIS1;extracted February 14, 2017; Supplementary Table S1). We chosetaxa that were common, ecologically relevant (e.g., ecosystemengineers) and/or characteristic taxa of Arctic shelves, andwhere biogeographic affinities are well-known (see below). Inaddition, taxa were chosen to cover the range of factors studiedhere (biogeographic affiliation, calcification status, and highertaxonomic level). Further, we required that each taxon have atleast 30 records within the geographic domain of SINMOD inorder to maintain MaxEnt model performance and confidencelevel. This cutoff appeared to be sufficient since the modelperformed quite well for taxa with close to 30 observations(e.g., Parastichopus tremulus and Serripes groenlandicus, seeSupplementary Tables S2, S3). Bottom water temperature andaragonite saturation state from the bias-corrected SINMOD1www.iobis.comhindcast output were matched to the georeferenced benthic faunadistributional records by interpolating from the model grid to thefaunal data positions and averaging over the 5 years precedingeach faunal observation date. Water depth at the samplingpoints was mostly (61%) from measurements when plausiblevalues were provided; the remaining 39% were estimatedusing the high-resolution bathymetry products. Each specieswas assigned status as calcifying or not, biogeographic affinity(Matishov et al., 2012; World Register of Marine Species)2, andhigher taxonomic grouping (mollusk, echinoderm, polychaete,crustacean, bryozoan).In locations where species data have been collectedsystematically, for example through biological monitoring,both presence and absence of species at each site have beenrecorded. However, most observations of the OBIS databasehave been collected non-systematically and are available aspresence-only records, and the different gear types deployedand study foci make it difficult to infer absence data from theserecords. We therefore employed a species distribution modelingapproach that requires only presence data, in order to maximizethe utility of the database.From the water depth, temperature, and aragonite saturationstate associated with each species record, we estimatedthree-dimensional realized niches and related them to theprojected changes in these parameters through the end of thecentury. The contribution of each environmental variable tospecies occurrence probabilities in the Arctic was calculatedusing the MaxEnt method. MaxEnt is a machine-learningalgorithm for modeling species distributions from presence-only2www.marinespecies.orgFrontiers in Marine Science | www.frontiersin.org 3 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 4Renaud et al. Species Distributions in a Changing Arctic Oceanspecies records. In brief, MaxEnt identifies the environmentalcharacteristics of locations with species occurrences that differfrom the environmental setting across the whole geographicalregion of interest. Based on the observed mismatch, a species’suitable habitat is defined. More specifically, the MaxEntmodel minimizes the relative entropy between two probabilitydensities (one estimated from the presence data, and onefrom the landscape) defined in covariate space compared toa uniform distribution null model (see Merow et al., 2013for details). MaxEnt’s predictive performance is consistentlycompetitive with the highest performing methods. Sincebecoming available in 2004, it has been utilized extensivelyfor finding correlates of species occurrences, mapping currentdistributions, and predicting to new times and places acrossmany ecological, evolutionary, conservation and biosecurityapplications (Elith et al., 2006).Sampling BiasesBy default, MaxEnt models assume that: (1) all locations on thelandscape (or “background locations”) were equally likely to besampled, and (2) any focal species would have been recorded atthe sampled landscape points if it had actually occurred there.Assumption (1) is likely violated in our occurrence datasets,because of logistical and practical constraints on sampling inthe Arctic (e.g., due to ice cover). Given such spatiotemporalbiases, one cannot differentiate whether species are observedin a particular environmental niche because those locationsare preferable or because they receive the largest search effort.Therefore, in order to relax assumption (1), landscape pointswere only drawn from the “target group sampling” (TGS)locations, which we defined as the total set of sampled locationsfor all species. However, this does not relax assumption (2),which may be violated where sampling gear and study focushave excluded certain species (Phillips et al., 2009). In theory thegear/study focus biases could be accounted for by partitioningthe TGS, but this was not considered feasible for our particulardataset due to insufficient samples.Model Fitting and ValidationMulticollinearity can be an issue with MaxEnt when answeringif and when environmental variables are of ecological interest.Thus, prior to modeling, a correlation analysis was conducted forenvironmental variables and the final MaxEnt models includedvariables that were not significantly correlated with each other atp < 0.05 (aragonite saturation state, temperature, water depth).Calcite saturation state was excluded as it correlates highly witharagonite saturation state.In this study MaxEnt models were fitted as combinationsof basic functions and features. MaxEnt had six feature classes:linear, product, quadratic, hinge, threshold, and categorical.Products were all possible pairwise combinations of covariates,allowing simple interactions to be fitted. Threshold featuresallowed a “step” in the fitted function, hinge features were similarexcept they allowed a change in gradient of the response. Manythreshold or hinge features were fitted for one covariate, giving apotentially complex function. Nevertheless, the MaxEnt programwas allowed to simplify the associations between species and theenvironment and in many models only one or a few features wereused, e.g., hinge, linear, and quadratic.Segment-based (non-gridded) data were modeled using SWD(samples-with-data) format in MaxEnt for both presence andbackground sites. A separate model was run for each species.A 10-fold cross-validation was used to obtain out-of-sampleestimates of predictive performance and estimates of uncertaintyaround fitted functions. In order to reduce model overfitting, abalance between accurate prediction (model fit) and generality(model complexity) was sought by maximizing the penalizedmaximum likelihood function, i.e., the gain function. Whendoing so, regularization or the LASSO penalty was appliedby exploring a range of regularization parameter values andchoosing a value that maximizes measures of fit on a cross-validation data set. The LASSO penalty is based on the rationalethat features with larger variance should incur a larger penaltyand, thus be less likely to be included in the model (Hastie et al.,2009). For model validation a random selection of 25% of theoverall localities of species occurrences were used. The percentcontributions of individual variables to the final model wereidentified with jackknife tests. The jackknife test evaluates howeach variable contributes to the “gain” of the MaxEnt‘s model (i.e.,improvement in penalized average log likelihood compared tonull model) (Elith et al., 2011). We used the raw output (default)as this output is the closest estimate of the probability that thespecies is present (Elith et al., 2011).Model Performance and PredictionsWe used the area under the receiver-operating-characteristiccurve (AUC) as a single measure of overall model accuracy that isnot dependent upon a particular threshold. The value of the AUCis between 0.5 and 1.0 with AUC = 1.0 indicating that the modelhas a perfect match and AUC = 0.5 indicating that model is nobetter than random (Fielding and Bell, 1997).The MaxEnt model was used to predict suitable habitat foreach species under current and future environmental conditions,and differences between these predictions were used to assessthe response to climate change of potential habitat area for eachspecies. We are aware that projecting future species distributionsusually involves extrapolating models to novel combinations ofenvironmental variables, and such projections should be treatedwith extreme caution. However, a comparison of current andfuture environmental niche space indicated that only 20% of thewhole study area would enter a completely novel environment, asdefined by moving beyond the 3D convex hull of the present-dayniche space distribution. We therefore consider any artifacts ofextrapolation will have only a minor effect on our results.Statistical Analyses of Model ResultsChange in suitable habitat area from the present (1978–2008) tothe future (2090–2099) was calculated from the MaxEnt modelfor each taxon. One-way analyses of variance (ANOVA) wereperformed for each of the three factors [calcification (as a t-test),biogeography, higher taxonomic level] on the change in habitatarea (both percentage and absolute) after testing for homogeneityof variances and normality among factor levels (in each casethese results were non-significant). Since single-factor analysesFrontiers in Marine Science | www.frontiersin.org 4 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 5Renaud et al. Species Distributions in a Changing Arctic Oceanmay obscure interaction effects or even effects of one factor wheneffects of a second factor dominate the variability, we also ran a 3-factor ANOVA on the data. Since we lacked degrees of freedom torun a fully factorial test, we used a main-effects model. Analyseswere performed in Statistica v. 13.RESULTSEnvironmental Settings ProjectionsProjections for bottom water temperature (T) and aragonitesaturation state ((ar)) indicate substantial changes over the21st century (Figure 1). The warmer Atlantic Water-influencedzone expands northward and eastward from its current core inthe southern and southwestern Barents Sea (Figures 1A,B), andbenthic habitat in the Barents Sea warms by up to 6◦C between1979–2008 and 2090–2099 (Figure 1C). Over the same period,bottom (ar) on the shelves (<500 m depth) decreases, mostlyby 0.6–1.1 units from around 1–2 (weakly saturated) to 0–1(undersaturated) (Figures 1D–F). Exceptions include the shallowRussian shelf regions (<50 m depth), which are already stronglyundersaturated ((ar) < 0.4; Figure 1D) and cannot becomemuch more undersaturated in the future (Figures 1E,F). Thelargest decreases in (ar) (∼1 unit) are in the North, Bering,and Chukchi Seas, and on the East Siberian Shelf at ∼50–500 m(Figure 1F). Arctic water-dominated areas of the northern andeastern Barents/Kara Seas remain relatively stable in temperature(Figure 1C), but experience (ar) reductions as large as in thecentral Barents Sea (Figure 1F). Deeper areas (>500 m) in theArctic and North Atlantic basins, and in Baffin Bay, see relativelysmall changes in benthic environment (Figures 1C,F).Depicting Realized and Predicted Nichesof Benthic TaxaOf the 78 taxa with sufficient (>30) records for evaluation, themean change in suitable habitat was small: approximately a 5%loss in suitable habitat based on MaxEnt. The range, however,was high, from a 42% increase in habitat (the amphipod Byblisgaimardii) to a 53% decrease (the sea cucumber P. tremulus)according to the 3-variable model (Table 1, SupplementaryTable S2, and Figures 2, 3). MaxEnt performed very well formost (∼70%) taxa, with AUC > 0.9 (Supplementary Table S3).Temperature contributed most to defining the niche space for72% of the taxa, while depth and aragonite saturation state werethe largest contributors for 19 and 8% of the taxa, respectively(Supplementary Table S3).Surprisingly, both calcifiers (in general and by phylum)and Arctic (compared to Arcto-boreal or Boreal) taxa werefound to be resilient to the projected environmental changesin terms of group mean response. Calcifiers were slightlymore vulnerable on average than non-calcifiers (6.1 vs. 2.8%mean habitat loss) but the difference was not significant dueto high intra-group variability (t-test, P = 0.18, Figure 4).Neither biogeographic affinity nor higher taxonomic level hada significant influence on the group-mean change in habitat(one-way ANOVA tests, P = 0.39 and 0.80, Figure 4). WhenArcto-Boreal and Arctic taxa were combined, their mean habitatTABLE 1 | Winners (>20% increase in suitable habitat, upper panel) and losers(>20% decrease in suitable habitat, lower panel) as predicted bythe MaxEnt model.Change in HigherTaxon habitat (%) N taxon Calcify? BiogeographyByblis gaimardii +42 356 Crustacean No ArctoBorealYoldiella frigida +40 41 Bivalve Yes ArcticAstarte crenata +22 115 Bivalve Yes ArctoBorealParastichopustremulus−53 30 Echinoderm Yes BorealNotoproctus sp. −36 210 Polychaeta No ArctoBorealHaploops tubicola −28 294 Crustacean No ArctoBorealMya truncata −27 655 Bivalve Yes ArctoBorealBalanus balanus −26 764 Crustacean Yes BorealPagurus pubescens −25 222 Crustacean Yes ArctoBorealYoldiella solidula −25 64 Bivalve Yes ArcticScoletoma impatiens −24 292 Polychaeta No BorealAstarte montagui −23 409 Bivalve Yes ArctoBorealHaploops sp. −22 579 Crustacean No ArctoBorealYoldiella nana −21 342 Bivalve Yes ArctoBorealThe number of records extracted from OBIS for each taxon (N), higher taxon,calcification status, and biogeographical zones are indicated for each taxon. SeeFigures 2, 3 for maps of present (1978–2008) and future (2090–2099) distributionsfor the two biggest winners and losers.loss was less than for Boreal taxa (2.9 vs. 10.9% habitat loss,P = 0.043). Typically, heavily calcified taxa from the Molluscaand Echinodermata phyla have been suggested to suffer greatestfrom acidification (Hale et al., 2011). Our analysis suggested thatonly 24 and 25% of the taxa in these phyla, respectively, wouldexperience range reductions, compared to 45–100% for otherphyla (Supplementary Table S2). Similarly, only 12% of Arctictaxa lost habitat, compared with nearly 50% of Arcto-boreal andBoreal taxa (Supplementary Table S2). Results were the same forall analyses whether the absolute or percent change in suitablehabitat was analyzed. The 3-factor main effects ANOVA alsoindicated no significant factor effects (biogeography: F = 2.1,df = 2, P = 0.14; calcification: F = 1.9, df = 1, P = 0.17; highertaxonomic level: F = 0.70, df = 4, P = 0.59).The data we extracted from the OBIS database were at boththe species and genus levels. In most cases, genus-level taxawere assigned to broader biogeographic zones (e.g., ArctoBoreal)because they include several species. We tested for this potentialbias by analyzing only the data for taxa identified to the species-level or where it was reasonably certain that the taxa found couldbe ascribed to either Boreal or Arctic zones. Results of theseanalyses with 23 fewer taxa were nearly identical to those forthe entire data set (3-factor main-effects ANOVA: biogeography:F = 2.4, df = 2, P = 0.10; calcification: F = 0.89, df = 1, P = 0.35;higher taxonomic level: F = 0.23, df = 4, P = 0.92).Frontiers in Marine Science | www.frontiersin.org 5 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 6Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 2 | Habitat suitability maps for the two biggest winners [the Arcto-boreal amphipod Byblis gaimardii (+42%): (A,B); and the Arctic bivalve Yoldiella frigida(+40%): (C,D)] in the study area based on the results from the MaxEnt model (see Table 1). Maps show present (A,C) and predicted future (B,D) distributions for thetwo species. All models predict the expected probability of occurrence from 0 (0% probability; blue) to 1 (100% probability; red).DISCUSSIONWe found little evidence to support the belief that Arctic taxa aremore sensitive than taxa with more southerly or cosmopolitandistributions, nor did we find that calcifiers as a group are moresensitive to projected climate-related changes than non-calcifiers.Models did indicate, however, that individual taxa can experienceconsiderable changes in suitable habitat under future scenarios,with potentially significant ecological consequences.Roles of Niche-Defining ParametersWhen fitting MaxEnt models, temperature was the mostimportant factor determining current suitable habitat for morethan two-thirds of the taxa studied. This should come as littlesurprise as temperature is one of the main factors used todelineate biogeographical provinces (Longhurst, 1998), so studiesover a large spatial scale with strong temperature gradientssuch as this are likely to show an important temperatureeffect. Temperature was also identified (by MaxEnt modeling)as the most important factor for mollusk and crustacean speciesdistributions at local scales in an Arctic fjord (Drewnik et al.,2017). Temperature can act in many ways, and may be viewedas a “master parameter” that works directly through thermaltolerance for survival and/or reproduction (e.g., Wethey andWoodin, 2008), or indirectly through its impacts on, e.g., primaryproduction, ice cover, and zooplankton and microbial grazingpressure, all of which could affect benthic food supply (Maarand Hansen, 2011). Shifts in community structure of BarentsSea benthos and fish have been noted, with the southern limitof Arctic communities moving northward during warming andretreating again southward during cooling periods (Blacker, 1965;Matishov et al., 2012; Fossheim et al., 2015). These changeswere attributed to temperature, but it is unclear whether theywere direct tolerance effects or something more indirect. It isimportant to note, however, that some high-latitude taxa havevery narrow temperature ranges over which they occur, makingthem particularly vulnerable to climate-driven changes (Pörtneret al., 2014; Morley et al., 2019).We used aragonite saturation state to represent the degreeof potential species response to ocean acidification, but weacknowledge that this may not be the only relevant carbonatesystem driver. Whereas most of the calcifying taxa we studiedproduce aragonitic skeletons (Mollusca: Kuklin´ski, unpublisheddata), there is also a large number of organisms included in ourstudy which use calcite with varying degrees of Mg incorporation(Echinodermata, Bryozoa, Cirripedia: Kuklinski and Taylor,2009; Iglikowska et al., 2017, 2018). In addition, some species maybe more sensitive to pH or pCO2 or combinations of OA stress.Frontiers in Marine Science | www.frontiersin.org 6 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 7Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 3 | Habitat suitability maps for the two biggest losers [the Boreal sea cucumber Parastichopus tremulus (-53%): (A,B); and the Arcto-boreal polycheteNotoproctus spp. (-36%): (C,D)] in the study area based on the results from the MaxEnt model (see Table 1). Maps show present (A,C) and predicted future (B,D)distributions for the two species. All models predict the expected probability of occurrence from 0 (0% probability; blue) to 1 (100% probability; red).However, we suspect that using only one parameter (aragonitesaturation state) provides a reasonable first approximation, sincemost changes are expected to be quite strongly correlatedbetween different carbonate system parameters (e.g., Figure 9 inWallhead et al., 2017).We included depth in our distribution models, even thoughit will not change significantly in the near future. Depthcan often reflect a number of linked parameters importantfor an organism’s habitat state, such as sediment grain size,bottom currents, and vertical flux of particulate organiccarbon (i.e., a proxy for food for many benthic organisms).Maintaining homeostasis under conditions of ocean acidification,for example, is thought to require more energy, and, assuch, may only be possible when sufficient resources areavailable for affected organisms (Pansch et al., 2014). Increasedmetabolic rates due to higher temperatures may enhancethis challenge (Pörtner, 2008). Primary production in theBarents Sea is predicted to change, both quantitatively andspatially, due to climate change (Slagstad et al., 2015). Clearly,inclusion of depth in the model will not account for all ofthis, but given the strong relationship between depth andvertical flux, we expect use of depth as one of the nicheparameters to some extent accounted for spatial differencesin food supply, in addition to partially controlling for otherhabitat parameters.Accuracy of Environmental DataThe pan-Arctic SINMOD model, used here to provide bottomwater temperature and aragonite saturation state, has beenextensively tested and bias-corrected using field measurements(Slagstad et al., 2015; Wallhead et al., 2017). We thereforeconsider that the hindcast data used to define faunal nichesare unlikely to have been a major source of bias. Projectedchanges for the 2090s are, however, subject to numerous modeluncertainties and an ensemble approach would likely be requiredto quantify these. Nevertheless, we know that SINMOD givesbottom water projections that are comparable with global andother downscaling models, and it is sufficient for this studythat the model gives plausible projections for the chosen climatechange scenario (SRES A1B). For end-of-century simulations,the dominant source of uncertainty is likely to be the climatechange scenario itself (Hawkins and Sutton, 2009). Waterdepth estimates for the pan-Arctic region are also not withoutuncertainty (e.g., the absolute differences between plausiblemeasured values and the interpolated bathymetry productshad median = 6 m, 95th percentile = 103 m) but it seemsFrontiers in Marine Science | www.frontiersin.org 7 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 8Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 4 | Mean (±1 SE) percent habitat loss for benthic taxa between the present (1978–2008) and the end of the 21st century (2090–2099). P-values of analysesof variance are included for comparisons of species grouped into: (A) biogeographic zone, (B) calcification status, and (C) higher taxonomic group. Number of taxain each group is indicated by the number in parentheses under the group name.unlikely that such errors could significantly bias our estimates ofmean habitat loss.Surprising Robustness of HabitatSuitabilityGiven the reported higher vulnerability of polar organismsto ocean warming (IPCC medium confidence, Pörtneret al., 2014) and of highly calcified organisms to oceanacidification (IPCC high confidence, Ibid), the fact thatour projections for benthic taxa at high northern latitudesshow no clear increase in habitat loss for Arctic orcalcifying species may seem surprising. It leads us toquestion whether something could have been lost in thecomplexities of the MaxEnt analysis, and if the results mightbe understood or corroborated by simpler analyses basedon the statistics of the environment where the differentspecies were found.A first reason to expect higher Arctic vulnerability is thatenvironmental changes for Arctic species may be larger underclimate change, due to rapid loss of ice cover/thickness andstronger warming in polar regions (Smetacek and Nicol, 2005;Brierley and Kingsford, 2009; Doney et al., 2012). However, ifwe consider the median projected bottom water temperaturechange at the sampled presence points for each species (seeFigure 5A), the species-mean warming is slightly higher fornon-Arctic species (2.3 vs. 1.8◦C, P = 0.04, t-test, unequalvariances). Figure 5A also shows that the warming is muchmore consistent for Boreal species than for non-Boreal species(standard deviation = 0.15 vs. 0.87◦C; P < 10−11, F-test).Regarding the median decreases in aragonite saturate state atsampled points (Figure 5B), the mean decreases are again slightlylarger for non-Arctic species (0.87 vs. 0.63, P< 10−3, t-test). So atleast for bottom water temperature and saturation state, it seemsthe environmental changes are not larger for Arctic species, asFrontiers in Marine Science | www.frontiersin.org 8 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 9Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 5 | Environmental statistics over sampled presence locations vs. percent habitat loss from MaxEnt for Boreal (red), Arcto-boreal (green), and Arctic species(blue), and for both calcifiers (triangles) and non-calcifiers (circles). X-axes show the median projected changes in temperature (A) and aragonite saturation state (B),the 95% ranges in hindcast temperature (C) and aragonite saturation state (D), the 95th percentiles of hindcast temperature (E), and the 5th percentiles of hindcastaragonite saturation state (F). Projected changes are from the SINMOD projections (2090s vs. 1978–2008 under SRES A1B) interpolated to the presence locations;hindcast variables are from the bias-corrected SINMOD hindcast, interpolated to the presence locations and averaged over the 5 years preceding each sample date.one might have anticipated from Figure 1. Furthermore, there areno significant correlations between the median environmentalchanges and the % habitat loss (Figures 5A,B).A second basis for high Arctic sensitivity is the idea thatpolar organisms have narrow environmental tolerance ranges,e.g., a “stenothermal” characteristic due to low energy-demandlifestyles and specific adaptations to cold water (Pörtner et al.,2014; Morley et al., 2019). We can test this by considering the2.5−97.5% ranges in bottom water temperature and aragonitesaturation state at the sampled presence points as tolerancemetrics for each species (Figures 5C,D). Here we find thatArctic species do have somewhat narrower ranges that non-Arcticspecies on average (6.0 vs. 9.6◦C, 1.2 vs. 1.8; P < 10−3, t-tests).However, variability in tolerance ranges is large within each group(Figures 5C,D), and there are no significant correlations between% habitat loss and tolerance ranges for temperature or saturationstate (or indeed water depth, not shown). The rather moderatedifferences in mean tolerance range may in part reflect a higherlevel of historical variability in Arctic vs. Antarctic systems, atleast for temperature (Pörtner et al., 2014). Arctic shelf systemshave experienced, even in the past 150 years, several cycles ofwarming and cooling within the range we are observing today,and Arctic benthic fauna have expanded or retreated on decadalscales accordingly (e.g., Blacker, 1965).A third basis for high Arctic sensitivity is the idea that polarorganisms will have their habitat “squeezed out” under globalchange, either because they cannot move to higher latitudes inresponse to warming (Pörtner et al., 2014; Renaud et al., 2015), orbecause critical thresholds of saturation state will be approachedfirst in the Arctic under ocean acidification (Feely et al., 2009;Steinacher et al., 2009; Ciais et al., 2013). If we consider the 95thpercentile of bottom water temperature at the presence locationsas a simple metric of upper thermal limit (Figure 5E), we findthat Arctic species do indeed have a significantly lower upperthermal limit on average (3.3 vs. 9.2◦C, P < 10−12, two-samplet-test). However, there is only a marginally significant correlationbetween % habitat loss and upper thermal limit (Pearson’sr = 0.23, P = 0.04). Considering the 5th percentile of aragonitesaturation state at presence locations as a lower tolerance limitand potential predictor of vulnerability to acidification, wefind no significant difference between Arctic vs. non-Arcticspecies, and no significant overall correlation with % habitatloss (Figure 5F).With regard to the paradigm of higher sensitivity of calcifyingspecies, we find no significant difference in species-mean 5thpercentiles of aragonite saturation state at presence locations forcalcifiers vs. non-calcifiers (0.71 vs. 0.67, P = 0.76, t-test), andboth categories contain species with 5th percentiles much lessthan 1, and in several cases < 0.2 (Figure 5F). Such toleranceof corrosive water may be partly explained by physiologicalmechanisms associated with the calcification process, particularlyin these often heavily calcified taxa, which are already adaptedfor internal control of carbonate saturation state at the site ofprecipitation (Hendriks et al., 2015). Considerable variability inresponse to ocean acidification has been shown in experimentalmanipulations (Kroeker et al., 2011; Vihtakari et al., 2016), andFrontiers in Marine Science | www.frontiersin.org 9 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 10Renaud et al. Species Distributions in a Changing Arctic Oceannaturally high environmental variability in some habitats(shallow waters, vegetated substrate) may also help to promotetolerance via phenotypic plasticity and local adaptations(Hofmann et al., 2010). Such mechanisms may make the speciesless susceptible to modest decreases in saturation state, even atunder-saturation levels, although food supply may mediate theability of an organism to regulate conditions for calcification(Ramajo et al., 2016).Analysis of Habitat Losses in NicheSpaceIt is worth noting that neither large environmental changes nornarrow tolerance ranges will necessarily result in high specieshabitat loss; if climate change produces enough new habitat,the net loss of habitat may be zero (or indeed negative) nomatter how large the changes or narrow the tolerances. Rather,large changes and narrow tolerances will likely result in large(negative or positive) changes in suitable habitat area. To betterunderstand the projections, it can be helpful to consider thecalculation of % habitat loss in niche space – here a 3D space of(temperature, saturation state, water depth). The niche space canbe divided into a grid of subvolumes, where the ith subvolume“contains” a physical area of habitat A(1)i in the present andA(2)i in the future. The MaxEnt algorithm essentially determinesan occupancy probability Pji for the jth species and ith niche-space subvolume by fitting a complex preference function to thepresent-day niche occupancy data (Merow et al., 2013). Assumingno change over time in this niche-space preference function(ruling out e.g., genetic adaptation), the net change in habitat areaHj for species j is given by:1Hj =∑iPji(A(2)i − A(1)i)It follows that the projected habitat losses can in large part beexplained by the changes in niche subvolume area content1Ai =A(2)i − A(1)i . These latter are entirely independent of species andthe MaxEnt procedure, and depend only on the present andfuture environment (here based on the SINMOD biogeochemicalmodel and the high-resolution bathymetry products).The habitat area changes over niche space can be visualizedon a 2D grid over temperature (T) and aragonite saturationstate ((ar)) by integrating over the water depth dimension(Figure 6A). The largest losses are in the red cells around (-2 ≤ T < 1◦C, 1 ≤ (ar) < 1.5), while the largest gains are inthe blue cells around (−2 ≤ T < 1◦C, 0 ≤ (ar) < 0.2) and(0 ≤ T < 5◦C, 0.4 ≤ (ar) < 1). In order to be a “winner”under the projected environmental change, a species needs tohave an occupancy probability or niche preference function Pjithat favors the niche subvolumes that gain habitat area. Thisrequires that its present-day sampling distribution over nichespace gives evidence that the species can tolerate the winning bluecells. For example, the biggest winner, the amphipod B. gaimardii(+42% habitat area), has a widespread niche space distribution(magenta crosses). Comparing this to the sampling distributionover all species (the TGS, black dots), the MaxEnt algorithminfers a strong tolerance of low (ar) and a thermal window thatcomfortably contains the winning range (0 ≤ T < 5◦C). Thisallows B. gaimardii to make large area gains in the future BarentsSea, Bering Strait, Chukchi Sea, and on the East Siberian shelfbetween 50 and 500 m (Figures 1B,E, 2A,B). By contrast, thebiggest loser, the sea cucumber P. tremulus (-53% habitat area),has a much more confined sampling distribution in niche space(magenta circles), suggesting limited tolerance of undersaturatedwater ((ar) < 1) and a thermal window of roughly 4–10◦C.Consequently, P. tremulus loses habitat in the North Atlantic andbecomes increasingly confined to the Norwegian shelf shallowerthan 500 m (Figures 3A,B).Figure 6A also sounds a note of caution regarding ourprojections: the sampling distribution of presences (black dots)is highly non-uniform over niche space, and in particular thewinning niche cells (blue colors) are strongly undersampled. Forseveral winning cells, the present-day niche area is zero (indicatedby a gray “x”); these are new niches that do not exist yet inthe pan-Arctic region. Although the MaxEnt procedure is rarelystrictly extrapolating with respect to the 3D convex hull (seesection “Materials and Methods”) it is often being used to makebold interpolations based on limited data. This issue, plus otherpotential sources of bias (e.g., gear biases), reduce our confidencein the projections for individual species habitat loss, although wesuspect a lesser impact of sampling bias on the significance ofdifferences between group means. Our first attempt at runningthe analysis used the default MaxEnt assumption of uniformsampling effort over the pan-Arctic grid (a clearly untenableassumption, see Figure 6B); this led to much larger group-meanhabitat losses (59–82%), but mean losses of Arctic species (75%)were still not significantly different to those of boreal (59%) orArcto-boreal species (82%). The results presented in this paperaccount for this first assumption.Further insights into vulnerability and sampling needs areprovided by interpolating the 3D niche space area changesonto the present-day (1978–2008 average) values over theSINMOD grid (Figure 6B). Here, any species whose present-day geographical distribution is confined to areas of niche arealoss (yellow/red colors) is set to be a loser under the projectedchanges. Hence, we can project, independently of any MaxEntanalysis, that species with environmental tolerances that confinethem to e.g., the northern Barents Sea, or the East GreenlandShelf, are likely to lose habitat under climate change. Conversely,any species that tolerates present-day conditions in the KaraSea near the Gulf of Ob, or the shallow areas (<50 m) of theeastern Laptev Sea and East Siberian Sea, is set to win. It isthus unfortunate that relatively few species presence records fromthese latter areas have made it into the OBIS database, and futureattempts to project future species distributions would likely bewell served by more sampling effort in and/or data recovery fromthese particular regions.Ecological Consequences ofDistributional ShiftsApproximately half the taxa were projected to experiencesubstantial changes in suitable habitat area (>10% increaseor decrease), and some of these species are characteristicFrontiers in Marine Science | www.frontiersin.org 10 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 11Renaud et al. Species Distributions in a Changing Arctic OceanFIGURE 6 | Niche habitat area loss in niche space (A) and mapped to the present-day environmental conditions over the SINMOD grid (B). Area losses (color scale,same in both subplots) are calculated from bias-corrected SINMOD output (2090s averages under SRES A1B vs. 1978–2008 averages from hindcast). Area lossesin niche space (A) are integrated over the water depth dimension. Black dots show the sampling distributions of the OBIS presence data in niche space (A),calculated from bias-corrected SINMOD averages over 5 years preceding each sample date, and the presence locations in geographic space (B). Magenta crosses(+) show the sampling distributions of the biggest winner, the Arcto-boreal amphipod Byblis gaimardii (+42% habitat area), while magenta circles (o) show thesampling distributions of the biggest loser, the Boreal sea cucumber Parastichopus tremulus (-53% habitat area). Gray crosses (x) in (A) show cells in niche spacethat contain no physical area in the present-day environmental conditions (from bias-corrected SINMOD averages over 1978–2008, summing over all water depths).or dominant on Arctic continental shelves. The polychetesChaetozone setosa and Spiochaetopterus typicus, the bivalveMedicula ferruginosa, the barnacle Balanus balanus, and thedecapods Pagurus pubescens, Hyas araneus, and Pandalus borealisare common and characteristic taxa in the open Barents Sea andadjacent fjords. S. typicus (12% increase in habitat) and B. balanus(26% loss) are key ecosystem engineers in that the physicalstructure they provide is important for ecosystem function.Pandalus is an important commercial shrimp in Arctic and borealwaters, and a 17% increase in habitat was predicted. If these taxaare as strongly affected as our analyses suggest, the impact onecosystem structure and services in the Barents region could besignificant over the 21st century.Surprisingly, one of the biggest “winners” in our analyses,the Arctic bivalve Yoldiella frigida (40% habitat increase), wasdescribed as one likely to experience substantial range retractionin an analysis based solely on temperature ranges (Renaudet al., 2015). However, the Renaud et al. (2015) projectionconsidered only the Barents Sea. The MaxEnt algorithm predictsa much more modest increase for the Barents Sea (9%), withmost of the pan-Arctic habitat increase driven by gains onthe Chukchi/Russian shelves and coastal waters off Canada,Greenland, and Iceland (Figure 2D). Also, Renaud et al.(2015) assumed a hard upper temperature limit of 4◦C; theMaxEnt algorithm infers some tolerance of warmer temperatures(thus allowing for the slight increase in the future BarentsSea) based on the limited presence records (41) within thehindcast timeframe (1978–2008). In fact the histogram shownin Figure 5 of Renaud et al. (2015) also suggests some toleranceto warmer waters, although this latter is based on August-only SINMOD output for 2000–2009 without bias correction ormatching to sample date.One caveat of our study is the exclusive focus on bottomwater environment, while many benthic organisms reproducevia planktonic larvae that live for hours to months in thepelagic zone. A recent study of an Antarctic barnacle indicatedthe need to consider both stages of multi-phasic species sinceconditions are likely to change differently in the two oceanrealms (Gallego et al., 2017). Registered presence in the OBISdatabase, as used in our study, integrates environmental influenceon multiple life-history stages/processes from fertilization, todevelopment and dispersal, to recruitment and juvenile stages,any of which can be acted upon by ecological (includingacidification) interactions in both the water column and sea floor.Another possible caveat is our use of only two physical/chemicalparameters. Many other such parameters (e.g., carbon flux,bottom-water oxygen content) may define where an organismcan exist, and ecological interactions cannot be expected tobe unchanged as species ranges change. Kroeker et al. (2012)found ocean acidification to alter competitive hierarchies ofmacroalgae, and this could have consequences for associatedfauna. In contrast to what has been observed for Barents Seafish communities (Fossheim et al., 2015), benthic communitychange is not likely to progress via wholesale replacement ofArctic fauna with boreal fauna. Highly mobile fish are expectedto move northward and southward in response to changingextent of Atlantic Waters, whereas benthic taxa have limitedmobility and range changes are more likely to vary by taxa asa consequence of life-span, dispersal ability, etc. New speciescombinations will likely lead to novel and unpredictable speciesinteractions, further affecting species distributions (Pecl et al.,2017). Finally, our models are based on fixed niches so thatpopulations follow the environmental change by shifting theirgeographical distributions. However, populations may adapt toFrontiers in Marine Science | www.frontiersin.org 11 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 12Renaud et al. Species Distributions in a Changing Arctic Oceantolerate novel conditions through selection by evolving traitsthat provide tolerance in the new environment. Although it istempting to project future distributions that take into accountevolutionary responses, scientists seldom know if, or how quickly,the climate-sensitive traits of populations can evolve (Merilä andHendry, 2014; Jezkova and Wiens, 2016).CONCLUSIONDespite acknowledged issues and uncertainties with using speciesdistribution models, especially in poorly sampled areas like theArctic, such models are accepted as useful in providing a firstapproximation of how climate change will affect biodiversity(Pearson and Dawson, 2003; Wiens et al., 2009). This study isnot meant to emphasize fine-scale species-specific changes, butrather to delineate expected shifts in broad patterns. The findingsreported here are surprising in the context of current paradigmsregarding susceptibility to climate change. In particular, we foundthat for benthic species under end-of-century ocean warmingand acidification: (i) mean habitat losses over taxonomic groupswere generally small (0–11%), (ii) Arctic benthic species were notsignificantly more vulnerable than boreal or Arcto-boreal species,(iii) calcifying species were not significantly more vulnerablethan non-calcifiers. The lack of sensitivity of our results tosuch taxonomic groupings was consistent with simpler analysesbased on magnitude of environmental change and observedtolerance windows, although such criteria were individually poorpredictors of the habitat loss results from MaxEnt. Samplingbias in the presence data was a major issue, and in generalwe expect projections from benthic species distribution modelsfor high northern latitudes to be strongly sensitive to how thisis dealt with. Analysis of niche space area changes, dependentonly on the environmental data, suggested that largest futurehabitat losses will be for taxa now inhabiting cold/oversaturatedniches (−2 < T < 3◦C, 1 < (ar) < 1.5, e.g., the present-daynorthern Barents Sea and East Greenland Shelf) while the largestgains will be for those currently found in cold/undersaturatedniches (−2 < T < 5◦C, (ar) < 1, e.g., the present-day easternLaptev and East Siberian Seas). The latter “winning” niches wereundersampled in our dataset, and we expect that further fieldsampling and data recovery for these niches will be especiallyhelpful for projecting future species distributions.There are few similar studies from high latitudes (but seeMorley et al., 2019 for Antarctic waters), mainly due to limitedinformation about species ranges and physiological tolerances,and regional models powerful enough to project environmentalparameters into the future. As more data become availableand models are developed, the tools we used here can beimproved, and predictions refined. As in other predictivemodeling approaches, MaxEnt is a process that uses data miningand probability to forecast outcomes. When doing so, theestablished relationships are only based on statistical dependencebetween environmental and distributional data, and thereforemay fail to account for physiological limits and biologicalinteractions. Development of hybrid modeling platforms thatintegrate species distribution models and experimental resultscan facilitate higher resolution distribution models for individualtaxa, which can be particularly useful for economically orecologically influential species (Kotta et al., 2019). This work isimportant as changes in ranges of some of these species may havesubstantial impacts on ecosystem services provided by benthiccommunities (Pecl et al., 2017).DATA AVAILABILITYThe datasets generated for this study are available on request tothe corresponding author. Model output data are available onZenodo (https://zenodo.org/record/3245446).AUTHOR CONTRIBUTIONSPR, MW-K, and PK developed the study idea. JK and MRperformed MaxEnt modeling. RB and DS developed theSINMOD simulations. PW analyzed and bias-corrected theSINMOD output. PK was the project leader. PR led work onthe manuscript. All authors contributed to the interpretation ofresults and editing of the manuscript.FUNDINGThis work was supported by the Polish-Norwegian ResearchProgramme operated by the National Centre for Research andDevelopment under the Norwegian Financial Mechanism2009–2014 in the frame of Project Contract No. Pol-Nor/196260/81/2013. Additional funding was provided by theNorwegian Foreign Ministry and its Arctic 2030 program (CoArcproject), the Norwegian Research Council (“Adapting CoastalZone Management to Ocean Acidification (ACIDCOAST)”;Grant No. 255748), the FRAM – High North Research Centrefor Climate and the Environment under the Ocean AcidificationFlagship, and the BONUS MARES project funded by theEuropean Union’s Seventh Framework Programme for research,technological development, and demonstration through BONUS,the Joint Baltic Sea Research and Development Programme (Art185), and Akvaplan-niva.ACKNOWLEDGMENTSThe authors are grateful for technical assistance provided byHector Andrade. Species occurrence and bias-corrected modeloutput data used in this manuscript are openly accessible atZenodo (https://zenodo.org/record/3245446). Further SINMODoutput from the runs used in this manuscript can be found athttps://doi.org/10.5281/zenodo.886077.SUPPLEMENTARY MATERIALThe Supplementary Material for this article can be found onlineat: https://www.frontiersin.org/articles/10.3389/fmars.2019.00538/full#supplementary-materialFrontiers in Marine Science | www.frontiersin.org 12 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 13Renaud et al. Species Distributions in a Changing Arctic OceanREFERENCESAMAP, (2017). SNOW, Water, Ice and Permafrost in the Arctic (SWIPA) 2017. Oslo:Arctic Monitoring and Assessment Programme (AMAP).AMAP, (2018). AMAP Assessment 2018: Arctic Ocean Acidification. Tromsø: ArcticMonitoring and Assessment Programme (AMAP).Ardyna, M., Babin, M., Gosselin, M., Devred, E., Rainville, L., and Tremblay,J. -É (2014). Recent Arctic Ocean sea ice loss triggers novel fallphytoplankton blooms. Geophys. Res. Lett. 41, 6207–6212. doi: 10.1002/2014GL061047Bates, N. R., Orchowska, M. I., Garley, R., and Mathis, J. T. (2012). Seasonalcalcium carbonate undersaturation in shelf waters of the Western Arctic Ocean;how biological processes exacerbate the impact of ocean acidification. Biogeosci.Discuss. 9, 14255–14290. doi: 10.5194/bgd-9-14255-2012Beaugrand, G., and Kirby, R. R. (2018). How do marine pelagic species respondto climate change? Theories and observations. Ann. Rev. Mar. Sci. 10, 169–197.doi: 10.1146/annurev-marine-121916-063304Bellerby, R., Anderson, L. G., Osborne, E., Steiner, N., Pipko, I., Chierici, M., et al.(2018). Arctic Ocean Acidification an Update. AMAP Assessment 2018: ArcticOcean Acidification. Tromsø: Arctic Monitoring and Assessment Programme(AMAP).Bellerby, R. G. J. (2017). Ocean Acidification without borders. Nat. Clim. Chang. 7,241–242. doi: 10.1038/nclimate3247Bellerby, R. G. J., Miller, L., Croot, P., Anderson, L., Azetsu-Scott, K.,McDonald, R., et al. (2014). Acidification of the Arctic Ocean in Arctic OceanAcidification Assessment. Oslo: Arctic Monitoring and Assessment Programme(AMAP).Bellerby, R. G. J., Silyakova, A., Nondal, G., Slagstad, D., Czerny, J., De Lange,T., et al. (2012). Marine carbonate system evolution during the EPOCAArctic pelagic ecosystem experiment in the context of simulated Arctic Oceanacidification. Biogeosci. Discuss. 9, 15541–15565. doi: 10.5194/bgd-9-15541-2012Blacker, R. W. (1965). Recent changes in the benthos of the west spitsbergen fishinggrounds. Special Pub. Int. Comm. Northw. Atl. Fish 6, 791–794.Brierley, A. S., and Kingsford, M. J. (2009). Impacts of climate change on marineorganisms and ecosystems. Curr. Biol. 19, R602–R614. doi: 10.1016/j.cub.2009.05.046Cheung, W. W., Lam, V. W., Sarmiento, J. L., Kearney, K., Watson, R., andPauly, D. (2009). Projecting global marine biodiversity impacts under climatechange scenarios. Fish Fish. 10, 235–251. doi: 10.1111/j.1467-2979.2008.00315.xCiais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., et al. (2013).“Carbon and other biogeochemical cycles,” in Climate Change 2013: Thephysical Science Basis. Contribution of Working Group I to the Fifth AssessmentReport of the Intergovernmental Panel on Climate Change, eds T. F. Stocker, D.Qin, G. K. Plattner, M. B. Tignor, S. K. Allen, J. Boschung, et al. (Cambridge,UK: Cambridge University Press), 465–570.Collins, M., Knutti, R., Arblaster, J., Dufresne, J.-L., Fichefet, T., Friedlingstein, P.,et al. (2013). “Long-term climate change: projections, commitments andirreversibility,” in Climate Change 2013: The Physical Science Basis. Contributionof Working Group I to the Fifth Assessment Report of the Intergovernmental Panelon Climate Change, eds T. F. Stocker, D. Qin, G.-K. Plattner, M. Tignor, S. K.Allen, J. Boschung, et al. (New York, NY: Cambridge University Press).Doney, S. C., Ruckelshaus, M., Duffy, J. E., Barry, J. P., Chan, F., English, C. A.,et al. (2012). Climate change impacts on marine ecosystems. Ann. Rev. Mar. Sci.4, 11–37.Drewnik, A., We¸sławski, J. M., and Włodarska-Kowalczuk, M. (2017). Benthiccrustacea and mollusca distribution in arctic fjord–case study of patterns inhornsund, svalbard. Oceanologia 59, 565–575. doi: 10.1016/j.oceano.2017.01.005Elith, J., Graham, C. H., Anderson, R. P., Dudík, M., Ferrier, S., Guisan, A.,et al. (2006). Novel methods improve prediction of species’ distributionsfrom occurrence data. Ecography 29, 129–151. doi: 10.1111/j.2006.0906-7590.04596.xElith, J., Phillips, S. J., Hastie, T., Dudik, M., Chee, Y. E., and Yates, C. J. (2011).A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57.doi: 10.1111/j.1472-4642.2010.00725.xFeely, R. A., Doney, S. C., and Cooley, S. R. (2009). Ocean acidification: presentconditions and future changes in a high-CO2 world. Oceanography 22, 36–47.doi: 10.1098/rstb.2012.0442Fielding, A. H., and Bell, J. F. (1997). A review of methods for the assessment ofprediction errors in conservation presence/absence models. Environ. Conserv.24:38e49.Fossheim, M., Primicerio, R., Johannesen, E., Ingvaldsen, R. B., Aschan, M. M.,and Dolgov, A. V. (2015). Recent warming leads to a rapid borealization offish communities in the Arctic. Nat. Clim. Chang. 5, 673–677. doi: 10.1038/nclimate2647Gallego, R., Dennis, T. E., Basher, Z., Lavery, S., and Sewell, M. A. (2017). Onthe need to consider multiphasic sensitivity of marine organisms to climatechange: a case study of the Antarctic acorn barnacle. J. Biogeogr. 44, 2165–2175.doi: 10.1111/jbi.13023Grebmeier, J. M. (2012). Shifting patterns of life in the Pacific Arctic and sub-Arctic seas. Ann. Rev. Mar. Sci. 4, 63–78. doi: 10.1146/annurev-marine-120710-100926Hale, R., Calosi, P., McNeill, L., Mieszkowska, N., and Widdicombe, S. (2011).Predicted levels of future ocean acidification and temperature rise could altercommunity structure and biodiversity in marine benthic communities. Oikos120, 661–674. doi: 10.1111/j.1600-0706.2010.19469.xHalpern, B. S., Longo, C., Hardy, D., McLeod, K. L., Samhouri, J. F., Katona, S. K.,et al. (2012). An index to assess the health and benefits of the global ocean.Nature 488, 615–620. doi: 10.1038/nature11397Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of StatisticalLearning. New York, NY: Springer.Hawkins, E., and Sutton, R. (2009). The potential to narrow uncertainty in regionalclimate predictions. Bull. Amer. Meteorol. Soc. 90, 1095–1107.Hendriks, I. E., Duarte, C. M., Olsen, Y. S., Steckbauer, A., Ramajo, L., Moore,T. S., et al. (2015). Biological mechanisms supporting adaptation to oceanacidification in coastal ecosystems. Est. Coast. Shelf Sci. 152, A1–A8.Hofmann, G. E., Barry, J. P., Edmunds, P. J., Gates, R. D., Hutchins, D. A., Klinger,T., et al. (2010). The effect of ocean acidification on calcifying organisms inmarine ecosystems: an organism-to-ecosystem perspective. Ann. Rev. Ecol. Evol.Syst. 41, 127–147. doi: 10.1146/annurev.ecolsys.110308.120227Hoppe, C. J. M., Wolf, K. K. E., Schuback, N., Tortell, P. D., and Rost, B.(2018). Compensation of ocean acidification effects in Arctic phytoplanktonassemblages. Nat. Clim. Chang. 8, 529–533. doi: 10.1038/s41558-018-0142-9Iglikowska, A., Najorka, J., Voronkov, A., Chełchowski, M., and Kuklin´ski, P.(2017). Variability in magnesium content in arctic echinoderm skeletons. Mar.Environ. Res. 129, 207–218. doi: 10.1016/j.marenvres.2017.06.002Iglikowska, A., Ronowicz, M., Humphreys-Williams, E., and Kuklin´ski, P. (2018).Trace element accumulation in the shell of the arctic cirriped Balanus balanus.Hydrobiology 818, 43–56. doi: 10.1007/s10750-018-3564-5IPCC, (2014). Climate Change 2014: Synthesis Report. Contribution of WorkingGroups I, II and III to the Fifth Assessment Report of the Intergovernmental Panelon Climate Change, eds Core Writing Team, R. K. Pachauri, and L. A. Meyer,(Geneva: IPCC).Jakobsson, M., Mayer, L., Coakley, B., Dowdeswell, J. A., Forbes, S., Fridman, B.,et al. (2012). The international bathymetric chart of the Arctic Ocean (IBCAO)version 3.0. Geophys. Res. Lett. 39:L12609. doi: 10.1029/2012GL052219Jezkova, T., and Wiens, J. J. (2016). Rates of change in climatic niches in plant andanimal populations are much slower than projected climate change. Proc. R.Soc. London B 283:20162104. doi: 10.1098/rspb.2016.2104Kotta, J., Vanhatalo, J., Jänes, H., Orav-Kotta, H., Rugiu, L., Jormalainen, V.,et al. (2019). Modelling spatial distribution patterns of macroalgae and theirgrazers under current and future environmental conditions by integration ofdistribution and experimental data. Sci. Rep. 9:1821.Kovacs, K. M., Lydersen, C., Overland, J. E., and Moore, S. E. (2011). Impactsof changing sea-ice conditions on arctic marine mammals. Mar. Biodiv. 41,181–194. doi: 10.1007/s12526-010-0061-0Kroeker, K. J., Micheli, F., and Gambi, M. C. (2012). Ocean acidification causesecosystem shifts via altered competitive interactions. Nat. Clim. Chang. 3:156.doi: 10.1038/nclimate1680Kroeker, K. J., Micheli, F., Gambi, M. C., and Martz, T. R. (2011). Divergentecosystem responses within a benthic marine community to ocean acidification.Frontiers in Marine Science | www.frontiersin.org 13 September 2019 | Volume 6 | Article 538fmars-06-00538 August 30, 2019 Time: 17:21 # 14Renaud et al. Species Distributions in a Changing Arctic OceanProc. Natl. Acad. Sci. U.S.A. 108, 14515–14520. doi: 10.1073/pnas.1107789108Kuklinski, P., and Taylor, P. D. (2009). Mineralogy of Arctic bryozoan skeletons ina global context. Facies 55, 489–500. doi: 10.1007/s10347-009-0179-3Longhurst, A. (1998). Ecological Geography of the Sea. San Diego: Academic Press.Maar, M., and Hansen, J. L. (2011). Increasing temperatures change pelagictrophodynamics and the balance between pelagic and benthic secondaryproduction in a water column model of the Kattegat. J. Mar. Syst. 85, 57–70.doi: 10.1016/j.jmarsys.2010.11.006Matishov, G., Moiseev, D., Lyubina, O., Zhichkin, A., Dzhenyuk, S., Karamushko,O., et al. (2012). Climate and cyclic hydrobiological changes of the BarentsSea from the twentieth to twenty-first centuries. Pol. Biol. 35, 1773–1790. doi:10.1007/s00300-012-1237-9Merilä, J., and Hendry, A. P. (2014). Climate change, adaptation, and phenotypicplasticity: the problem and the evidence. Evol. Appl. 7, 1–14. doi: 10.1111/eva.12137Merow, C., Smith, M. J., and Silander, J. A. Jr. (2013). A practical guide toMaxEnt for modeling species’ distributions: what it does, and why inputsand settings matter. Ecography 36, 1058–1069. doi: 10.1111/j.1600-0587.2013.07872.xMorley, S. A., Barnes, D. K., and Dunn, M. J. (2019). Predicting which speciessucceed in climate-forced polar seas. Front. Mar. Sci. 5:507.Nakicenovic, N., Alcamo, J., and Davis, G. (2000). Special Report on EmissionsScenarios, Intergovernmental Panel on Climate Change. (New York, NY:Cambridge University Press), 599.Neukermans, G., Oziel, L., and Babin, M. (2018). Increased intrusion of warmingatlantic water leads to rapid expansion of temperate phytoplankton in theArctic. Global Chang. Biol. 24, 2545–2553. doi: 10.1111/gcb.14075Pansch, C., Schaub, I., Havenhand, J., and Wahl, M. (2014). Habitat traits andfood availability determine the response of marine invertebrates to oceanacidification. Glob. Chang. Biol. 20, 765–777. doi: 10.1111/gcb.12478Pearson, R. G., and Dawson, T. P. (2003). Predicting the impacts of climate changeon the distribution of species: are bioclimate envelope models useful? GlobalEcol. Biogeogr. 12, 361–371. doi: 10.1046/j.1466-822x.2003.00042.xPecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen,I. C., et al. (2017). Biodiversity redistribution under climate change: impacts onecosystems and human well-being. Science 355:eaai9214. doi: 10.1126/science.aai9214Pedersen, M. W., Kokkalis, A., Bardarson, H., Bonanomi, S., Boonstra, W. J.,Butler, W. E., et al. (2016). Trends in marine climate change research in theNordic region since the first IPCC report. Clim Chang. 134, 147–161. doi:10.1007/s10584-015-1536-6Phillips, S., Dudík, M., Elith, J., Graham, C. H., Lehmann, A., Leathwick, J.,et al. (2009). Sample selection bias and presence-only distribution models:implications for background and pseudoabsence data. Ecol. Appl. 19, 181–197.doi: 10.1890/07-2153.1Pörtner, H. O. (2008). Ecosystem effects of ocean acidification in times of oceanwarming: a physiologist’s view. Mar. Ecol. Prog. Ser. 373, 203–217. doi: 10.3354/meps07768Pörtner, H.-O., Karl, D. M., Boyd, P. W., Cheung, W. W. L., Lluch-Cota, S. E.,Nojiri, Y., et al. (2014). “Part A: global and sectoral aspects. contributionof working group II to the fifth assessment report of the intergovernmentalpanel on climate change ocean systems,” in Climate Change 2014: Impacts,Adaptation, and Vulnerability, eds C. B. Field, V. R. Barros, D. J. Dokken, K. J.Mach, M. D. Mastrandrea, T. E. Bilir, et al. (Cambridge: Cambridge UniversityPress).Pörtner, H. O., Langenbuch, M., and Michaelidis, B. (2005). Synergistic effects oftemperature extremes, hypoxia, and increases in CO2 on marine animals: fromearth history to global change. J. Geophys. Res. Oceans 110:C09S10.Ramajo, L., Pérez-León, E., Hendriks, I. E., Marbà, N., Krause-Jensen, D.,Sejr, M. K., et al. (2016). Food supply confers calcifiers resistance to oceanacidification. Sci. Rep. 6:19374. doi: 10.1038/srep19374Renaud, P. E., Carroll, M. L., and Ambrose, W. G. Jr. (2007). “Effects of globalwarming on arctic sea-floor communities and its consequences for highertrophic levels,” in Impactos Del Calentamiento Global Sobre Los EcosistemasPolares, ed. C. Duarte, (Bilbao: FBBVA Press).Renaud, P. E., Daase, M., Banas, N. S., Gabrielsen, T. M., Søreide, J. E., Varpe, Ø,et al. (2018). Pelagic food-webs in a changing arctic: a trait-based perspectivesuggests a mode of resilience. ICES J. Mar. Sci. 75, 1871–1881. doi: 10.1093/icesjms/fsy063Renaud, P. E., Sejr, M. K., Bluhm, B. A., Sirenko, B. I., and Ellingsen, I. H.(2015). The future of Arctic benthos: expansion, invasion, and biodiversity.Prog. Oceanogr. 139, 244–257. doi: 10.1016/j.pocean.2015.07.007Slagstad, D., Ellingsen, I. H., and Wassmann, P. (2011). Evaluating primaryand secondary production in an Arctic Ocean void of summer sea ice: anexperimental simulation approach. Prog. Oceanogr. 90, 117–131. doi: 10.1016/j.pocean.2011.02.009Slagstad, D., Wassmann, P. F., and Ellingsen, I. (2015). Physical constrains andproductivity in the future Arctic Ocean. Front. Mar. Sci. 2:85.Smetacek, V., and Nicol, S. (2005). Polar ocean ecosystems in a changing world.Nature 437:362. doi: 10.1038/nature04161Steinacher, M., Joos, F., Frolicher, T. L., Plattner, G., and Doney, S. (2009).Imminent ocean acidification in the Arctic projected with the NCAR globalcoupled carbon cycle-climate model. Biogeosci. 6, 515–533. doi: 10.5194/bg-6-515-2009Sunday, J. M., Bates, A. E., and Dulvy, N. K. (2012). Thermal tolerance andthe global redistribution of animals. Nat. Clim. Chang. 2:686. doi: 10.1038/nclimate1539Vihtakari, M., Havenhand, J., Renaud, P. E., and Hendriks, I. E. (2016). Variableindividual-and population-level responses to ocean acidification. Front. Mar.Sci. 3:51.Wallhead, P. J., Bellerby, R. G. J., Silyakova, A., Slagstad, D., and Polukhin, A. A.(2017). Bottom water acidification and warming on the western Eurasian Arcticshelves: dynamical downscaling projections. J. Geophys. Res. 122, 8126–8144.doi: 10.1002/2017JC013231Wassmann, P. (2015). Overarching perspectives of contemporary and futureecosystems in the Arctic Ocean. Prog. Oceanogr. 139, 1–12. doi: 10.1016/j.pocean.2015.08.004Wethey, D. S., and Woodin, S. A. (2008). “Ecological hindcasting of biogeographicresponses to climate change in the european intertidal zone,” in Challenges toMarine Ecosystems, eds J. Davenport, G. Burnell, T. Cross, M. Emmerson, R.McAllen, R. Ramsay, et al. (Dordrecht: Springer).Wiens, J. A., Stralberg, D., Jongsomjit, D., Howell, C. A., and Snyder, M. A.(2009). Niches, models, and climate change: assessing the assumptions anduncertainties. Proc. Natl. Acad. Sci. U.S.A. 106, 19729–19736. doi: 10.1073/pnas.0901639106Conflict of Interest Statement: The authors declare that the research wasconducted in the absence of any commercial or financial relationships that couldbe construed as a potential conflict of interest.Copyright © 2019 Renaud, Wallhead, Kotta, Włodarska-Kowalczuk, Bellerby,Rätsep, Slagstad and Kuklin´ski. This is an open-access article distributed under theterms of the Creative Commons Attribution License (CC BY). The use, distributionor reproduction in other forums is permitted, provided the original author(s) andthe copyright owner(s) are credited and that the original publication in this journalis cited, in accordance with accepted academic practice. No use, distribution orreproduction is permitted which does not comply with these terms.Frontiers in Marine Science | www.frontiersin.org 14 September 2019 | Volume 6 | Article 538

The thermophilization of Arctic marine ecosystems from climate change

Abstract

Extracted data

The thermophilization of Arctic marine ecosystems from climate change

Abstract

Extracted data

Related items

Related items