Projected 21st‐century distribution of canopy‐forming seaweeds in the Northwest Atlantic with climate change

Climate change is predicted to alter the distribution and abundance of marine species, including canopy‐forming seaweeds which provide important ecosystem functions and services. We asked whether continued warming will affect the distribution of six common canopy‐forming species: mid‐intertidal fucoids (Ascophyllum nodosum, Fucus vesiculosus), low‐intertidal Irish moss (Chondrus crispus), subtidal laminarian kelps (Saccharina latissima, Laminaria digitata) and the invasive Codium fragile.

Marine range shifts will affect habitat-forming species in coastal habitats, with compounding impacts on the species who depend on their ecosystem structure, functions and services (Wernberg et al., 2016;Witman & Lamb, 2018). Significant losses of habitat-forming species have already occurred in coral reefs (Carpenter et al., 2008), mangrove forests (Polidoro et al., 2010), seagrass meadows (Short et al., 2011) and seaweed beds (Krumhansl et al., 2016) due to multiple factors, and climate change poses an additional threat to many of these foundation species.
One way to examine how species may respond to climate change is with species distribution models (SDM). SDM are designed around the niche concept, where a species persists and maintains a stable population size under a given set of abiotic and biotic conditions (Hutchinson, 1957), and link known records of species presence to environmental variables . SDM are powerful tools to understand how species may respond to climate change, particularly if model projections can be compared to physiological data which allow for more robust projections (Elith, Kearney, & Phillips, 2010;Talluto et al., 2016). Range shift predictions based on laboratory-experiments with temperature manipulation (Piñeiro-Corbeira, Barreiro, Cremades, & Arenas, 2018) have matched observed range shifts of seaweeds in the Northeast Atlantic (NE-Atlantic) due to increasing sea surface temperature (SST; Piñeiro-Corbeira, Barreiro, & Cremades, 2016). To date, despite the widespread use of SDM for habitat-forming species (e.g., Marcelino & Verbruggen, 2015;Record, Charney, Zakaria, & Ellison, 2013;Valle et al., 2014), only two studies have incorporated the use of physiological data (Franco et al., 2018;Martínez, Arenas, Trilla, Viejo, & Carreño, 2014). The incorporation of physiological data is particularly important for invasive species which are in violation of the assumption that a species is at equilibrium with the environment .

K E Y W O R D S
Chondrus crispus, Codium fragile, community composition, fucoid, kelp, maximum entropy, physiological threshold, range shift, rocky shore, species distribution model project future distributions under different climate-change scenarios; and (c) compare the SDM to physiological thresholds to assess the accuracy of projections and identify areas of stable growth versus reduced growth and/or mortality. We hypothesized species-specific northward range shifts, which would alter the composition of canopy-forming seaweeds along the NW-Atlantic coast. Any range shifts would have important consequences for ecosystem structure, functions and services, with implications for management and conservation.
We collected presence-only occurrence records from the literature and regional experts (Supporting Information Tables S1.1-S1.6 in Appendix S1) to determine each species present-day distribution.
Generally, records were collected from 1980 onwards, as this was the first year significant increases in SST were detected in the NW-Atlantic (Barnett et al., 2001). However, few occurrence records exist for the Arctic due to their inaccessibility. Therefore, additional records prior to 1980 were used north of Newfoundland (47°N), where seaweeds that were present prior to 1980 could still occur after 1980 as any Arctic warming would promote northward range shifts. As the environmental layers do not perfectly match the coastline, occurrence records that did not overlay the environmental data were moved to the closest pixel. All records were projected into the Behrmann cylindrical equal-area projection (see following section) using arcGIS ® 10.5 (ESRI, 2011), and duplicate records within a pixel were removed.

| Environmental data
The global marine environmental dataset Bio-Oracle 2.0 was used to represent present-day conditions which provide long-term averages (LTA) from 2000 to 2014, in Behrmann cylindrical equal-area projection, at a 7 km pixel resolution (Assis, Tyberghein et al., 2018b;Tyberghein et al., 2012). All pixels >100 km from shore were excluded from analysis and environmental data were cropped from 32-84°N and 42-95°W within the statistical environment R (R Core Team, 2014) using the "sdmpredictors" package (Bosch, Tyberghein, & De Clerck, 2018). To incorporate the range of environmental conditions seaweeds encounter during the year, we incorporated LTA of the maximum (e.g., warmest) month (M-max) and the minimum (e.g., coolest) month (M-min) for SST, sea surface salinity (SSS), sea ice coverage (SIC), current velocity (CV) and sea surface nitrate (NO 3 − ) for the three intertidal species (fucoids, C. crispus) as well as yearly LTA minimum and maximum (Y-max, Y-min) diffuse attenuation (DA) for the subtidal species (kelps, C. fragile), all of which are known to impact seaweed growth and survival (Harley et al., 2012).
To project onto future climate, we choose two global climate models from CMIP5 (Taylor, Stouffer, & Meehl, 2012), one projecting lower (GFDL-ESM2M; GFDL hereafter; Dunne et al., 2012) and one higher (IPSL-CM5A-LR; IPSL hereafter; Dufresne et al., 2013) levels of SST warming (Bopp et al., 2013). In addition, we used two representative concentration pathways (RCP): a low carbon emission and strong mitigation scenario (RCP2.6) and a high emission business-as-usual scenario (RCP8.5; Moss et al., 2010). This combination of two climate models and two RCPs provides a range of best-to worst-case scenarios of warming over the next century (Bopp et al., 2013). LTAs of SST M-max and M-min were derived for present-day (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), mid-century (2040-2050) and end-century (2090-2100) for each climate model and RCP scenario. Present-day was chosen to begin in 2006 as this was the year CMIP5 climate models switch from historical runs to future projections forced by each RCP scenario. To account for the differences between present-day observed (from Bio-Oracle) and projected (from GFDL and IPSL) environmental data, we calculated relative changes in SST M-max and M-min for future periods (mid-century minus present-day; end-century minus present-day) for each climate model and RCP scenario.
This relative change (anomaly) in SST was then added to the presentday observed SST. All projections are based on ocean temperature and the climate model data were re-gridded to the 7 km grid used for present-day environmental layers.

| SDM building, evaluation and projection
The correlative SDM was built using Maxent 3.3.3k (Phillips, Anderson, & Schapire, 2006;Phillips, Dudík, & Schapire, 2004) within R using the "dismo" package (Hijmans, Phillips, Leathwick, & Elith, 2017). Maxent consistently performs well compared to other TA B L E 1 Climate variables included in the variable selection process with their percent contribution (PC; %), permutation importance (PI; %), average training and test gain and AUC, and the value for maximizing the sum of test sensitivity and specificity (MTSS) logistic threshold (LT; bolded value indicates binomial probability of p < 0.001) for the cross-validation runs (k = 5) and full run (k = 1) indicating the order of variable removal (OR; bolded variables included in best model). See methods and Supporting Information   correlative presence-only modelling algorithms (Elith et al., 2006).
Each species present-day distribution was built using k-fold crossvalidation with k = 5 (Kohavi, 1995), allowing for linear, quadratic and hinge feature classes Phillips et al., 2006), and set to 5,000 iterations. Continuous model output was turned into binary presence/absence data based on maximizing the sum of test sensitivity and specificity threshold (Liu, Newell, & White, 2016). Model performance was evaluated using a one-tailed binomial test, and with the area under the curve (AUC) of the receiver operator characteristics curve (Phillips et al., 2006). A jackknife test as implemented Once the best SDM was determined per species, each species' SDM was rerun using the full set of occurrence records to ensure the models were trained under the full set of environmental conditions they may experience (k = 1). This was repeated for present-day observed conditions, and relative projected mid-century and endcentury conditions for both climate models and RCP scenarios. The future projections used present-day environmental conditions for all variables but switching the SST terms for the climate data calculated above. We projected with SST terms for all six species as SST M-max was retained as an important climatic driver within the best model per species, contributing at least 5% to the models, and warming ocean temperatures are a major driver of climate change. Two species were also projected with SST M-min as this was also retained in their best models (see below). Models were clamped to avoid training with data outside the species range (Elith et al., 2010).

| Physiological thresholds
A previous laboratory experiment (Wilson et al., 2015) provided information on the growth and survival of all study species except S. latissima between 12 and 29°C from populations found in the middle of each species respective range in Atlantic Canada. Saccharina latissima was assumed to experience similar growth as L. digitata based on experimental studies completed on S. latissima populations from Halifax (Bolton & Lüning, 1982;Simonson, Scheibling, & Metaxas, 2015), Maine, and Long Island Sound (Redmond, 2013).
These data were used to create three physiological thresholds (PTs) for heat-related growth and survival to define the water temperatures at which each species experienced reduced growth, reduced growth and partial mortality, and complete mortality. Reduced growth was defined as the first temperature at which there was a statistically significant growth rate reduction, partial mortality was at least one replicate dying, and complete mortality was all replicates dying at a specific temperature. See Supporting Information Table   S2.1 and Figure S2.1 in Appendix S2 for the data supporting our PT.
Cold-related survival was not tested but all species exhibit measurable photosynthesis in water temperatures of at least 0°C (Lüning, 1990). Present-day and projected future SST M-max data were classified into these three categories, creating a layer for each species' PT (Martínez et al., 2014).
We compared our PTs to the results of the SDM which allowed us to confirm the accuracy of projected range shifts based on SST changes (Franco et al., 2018;Martínez et al., 2014). To do so, the predicted present-day and projected future SDM were overlaid onto the corresponding species' PT. Areas corresponding to PTs for stable and reduced growth were likely to contain suitable habitat, while areas corresponding to PTs for reduced growth and partial mortality, or complete mortality were less likely to contain suitable habitat and may indicate areas with greater susceptibility to multiple stressors in a warming climate. In the following written section, we classify all areas that do not indicate stable growth as "poor growth" for simplicity. The reader is directed to Table 2 and the associated figure per species for the exact PT category and associated temperature.
Seaweed thermal performance curves are useful tools to predict the response to projected warming scenarios (Harley et al., 2012), and our PTs were used to reflect growth optima and important inflection points. We displayed SSTs in 1°C increments to account for the possibility that an inflection point occurred within the 3°C temperature range for each PT, and for the possible acclimation to different SST across a species range.

| Variable importance
We considered 12 environmental variables to be included in the variable selection process for the subtidal species, and 10 for the intertidal species (Supporting Information Tables S2.2 and S2.3 in Appendix S2). Each variable was removed in a backward selection process based on the jackknife test implemented by Maxent (
Codium fragile is not known to extend southward to the predicted southern (warm edge) range limit of ~32°N (Figure 6a). The occurrence records used to train the SDM aligned well with the PTs, with no records occurring at SST M-max values indicating complete mortality ( Figure 7). Furthermore, the allocation of the predicted area to SST M-max for poor growth was <3% for fucoids and C. crispus, <7% for the kelps, and ~12% for C. fragile (Table 2).
F I G U R E 1 Distribution of Ascophyllum nodosum: (a) present-day occurrence records (black dots) and SDM, and projected (b, c) endcentury distribution based on SST changes under RCP8.5 from GFDL-ESM2M (left column) and IPSL-CM5A-LR (right column; Behrmann world equal-area projection). Physiological thresholds were overlaid to show areas of stable growth (green, 12-22°C) and reduced growth (yellow-orange, 23-25°C); no areas of reduced growth and partial mortality (pink-red, 26-28°C) or complete mortality (dark red, ≥29°C) were observed. Pixels that appear to be on land are in narrow fjords The northern distribution limit of A. nodosum extends along the Labrador coast to ~67°N (Figure 1a). Chondrus crispus ( Figure 2a) and C. fragile (Figure 6a) occur north to ~51°N in southern Labrador, and the predicted habitat areas aligned well with known occurrence records. Laminaria digitata's northern range was slightly underpredicted in western Hudson Bay, but based on occurrence records, its northern limit is ~62°N (Figure 3a). The northern range limits of F. vesiculosus (Figure 5a) and S. latissima

| Seaweed flora under a strong mitigation scenario (RCP2.6)
Sea surface temperature M-max (mean SST of the warmest month on average) was retained in the best models by all six species, was in the top four variables based on test gain and contributed >5% to permutation importance and percent contribution (  Figure S2.4) leading to an average increase of 30% (118,605 km 2 ) of total habitat area by 2100 (  Figure   S2.5), C. crispus (Supporting Information Figure S2.6), and C. fragile (Supporting Information Figure S2.7) leading to an average increase of total habitat area by 1%, 16% and 6% by 2100, respectively F I G U R E 3 Distribution of Laminaria digitata: (a) present-day occurrence records (black dots) and SDM, and projected (b, c) end-century distribution based on SST changes under RCP8.5 from GFDL-ESM2M (left column) and IPSL-CM5A-LR (right column; Behrmann world equal-area projection). Physiological thresholds were overlaid to show areas of stable growth (green, 12-19°C), reduced growth and partial mortality (pink-red, 20-22°C), and complete mortality (dark red, ≥23°C). Pixels that appear to be on land are in narrow fjords End-Century RCP8.5 End-Century RCP8.5 Table 2). Laminaria digitata (Supporting Information Figure S2.8) was projected to lose 4% of its average total habitat area while S. latissima (Supporting Information Figure S2.9) showed little change (Table 2).

| Seaweed flora under a business-as-usual scenario (RCP8.5)
By end-century, northward shifts of the southern (warm) range limits were projected for all species except for C. fragile (Table 2). Unlike   were projected to have northward range expansions leading to an overall increase of habitat by 10% and 74%, respectively (Table 2).
In comparison, C. fragile was projected to experience no southern range shift and minimal northward range expansion that corresponded to an overall habitat loss of 2% (Figure 6b,c).
These species-specific range shifts cumulated to a shift in sea-

| D ISCUSS I ON
Climate change will shift the distribution of many marine species, including foundation species or ecosystem engineers who play key roles in providing habitat structure and other ecosystem functions to associated communities as well as services for human well-being (Pecl et al., 2017). Using SDM and PTs, we projected species-specific range shifts of canopy-forming seaweeds which will alter rocky shore communities along NW-Atlantic coasts. This included a shift in dominance from native to invasive species and a transition to opportunistic life histories. The projected response was dependent on the magnitude of carbon emissions and associated warming, with larger range shifts under the business-as-usual (RCP8.5) compared to the strong mitigation scenario (RCP2.6), highlighting the benefits of climate mitigation. We discuss the projected restructuring of the habitat-forming seaweed community and the consequences this will have for the functions and services they provide.

| Current drivers of seaweed distribution
Our while others found the reverse (Franco et al., 2018;Jueterbock et al., 2013). Regardless of season, the dominance of SST as an important predictor was expected as seaweed distributions are primarily defined by SST, and closely follow SST isotherms (Lüning, 1990;van den Hoek, 1975). Increasing summer SST is the primary driver of shifts of the southern (warm) range limit, as SST surpasses growth and mortality thresholds. Whereas, increasing winter SST facilitates faster growth rates (Bolton & Lüning, 1982;Marbà et al., 2017) and greater recruitment (Filbee-Dexter et al., 2016). Diffuse attenuation and sea ice coverage were also important predictors for seaweed distribution. Diffuse attenuation is an indicator of light availability, which sets seaweeds maximum depth limit. Sea ice coverage influences both the upper vertical distribution limit through ice scour, which often results in bare intertidal zones in the Arctic, and the lower vertical limit by reducing light availability (Krause-Jensen et al., 2012;Küpper et al., 2016). Unlike SST, diffuse attenuation and sea ice coverage are not commonly used as predictor variables; only one study has incorporated diffuse attenuation for fucoids (Jueterbock et al., 2013) and sea ice coverage for kelps (Assis, Araújo et al., 2018a).
As mean SST of the warmest month on average (M-max) is an important driver of seaweed distribution and incorporated into each species' SDM, we compared SDM output to warm-temperature PTs to assess the accuracy of our projections. The PTs were based on a standardized laboratory experiment (Wilson et al., 2015) and agreed F I G U R E 5 Distribution of Fucus vesiculosus: (a) present-day occurrence records (black dots) and SDM, and projected (b, c) end-century distribution based on SST changes under RCP8.5 from GFDL-ESM2M (left column) and IPSL-CM5A-LR (right column; Behrmann world equalarea projection). Physiological thresholds were overlaid to show areas of stable growth (green, 12-25°C), and reduced growth and partial mortality (pink-red, 26-28°C); no areas of complete mortality (dark red, ≥29°C). Pixels that appear to be on land are in narrow fjords with model output for fucoids, C. crispus, and C. fragile but suggested some areas of projected future habitat for kelps corresponding to growth reductions and mortality. Therefore, we may be underestimating the southern edge shift for kelps. Unfortunately, our warmwater PTs cannot provide any information on the species northern range limit; thus, future research could focus on establishing coldwater and light requirement PTs to further refine northern range limits. Other studies have also observed strong agreement between PTs and correlative SDM output (Franco et al., 2018;Martínez et al., 2014). However, neither the SDM nor the PT incorporated the species' potential to acclimate or adapt to increasing SST. Seaweeds are very responsive to acclimatization which can result in greater tolerance to temperature stress (Harley et al., 2012). Less is known about their ability to adapt to long-term increases in SST, but the occurrence of regional ecotypes suggests it is possible across evolutionary time scales. Yet climate change may occur at rates too fast to allow for adaptive evolution via genetic changes and other adaptation mechanisms may be required . For reefbuilding corals, SDM projected minimal habitat loss by end-century if corals adapted to 1°C warmer SST (Cacciapaglia & van Woesik, 2015), and seaweeds may respond in a similar adaptive manner.
Our projections are based on changes in ocean temperature, which will impact the distribution of intertidal and subtidal species.
Yet, air temperature also impacts the distribution of intertidal species F I G U R E 6 Distribution of Codium fragile: (a) present-day occurrence records (black dots) and SDM, and projected (b, c) end-century distribution based on SST changes under RCP8.5 from GFDL-ESM2M (left column) and IPSL-CM5A-LR (right column; Behrmann world equalarea projection). Physiological thresholds were overlaid to show areas of stable growth (green, 12-25°C), and reduced growth (yellow-dark orange, ≥26°C). Pixels that appear to be on land are in narrow fjords
Other studies have examined the response of seaweeds to continued SST increases in the NW-Atlantic using different climate models and scenarios. Using older carbon emission scenarios (B1, A1B, A2), Jueterbock et al. (2013) (Jueterbock et al., 2013) and RCP8.5 (Assis et al., 2014). Similar to our study, projections for C. crispus based on changes in SST isotherms for the B1 scenario also projected a shift to the Gulf of Maine by 2100 (Müller et al., 2009). B1 projections for S. latissima (Müller et al., 2009) and RCP8.5 projections for S. latissima and L. digitata (Assis, Araújo et al., 2018a) suggested southern range shifts for kelp to Newfoundland, whereas our study projected their existence further south in the northern Gulf of Maine with possible extirpations in parts of Nova Scotia and Gulf of St. Lawrence.
No previous study has projected the response of C. fragile. Thus, most studies agree on a general northward shift of seaweed distributions, while there is some uncertainty about the magnitude of species-specific shifts. This uncertainty may be driven by the inherent differences and biases between different global climate models and emission scenarios (Bopp et al., 2013). It may also be driven by the intraspecific differences in the responses of NW-and NE-Atlantic populations, which are genetically distinct (Assis et al., 2014;Neiva et al., 2018;Olsen et al., 2010), yet our study trained and projected SDMs within the NW-Atlantic and not the entire North Atlantic as other studies. Lastly, this uncertainty may be due to our calculation of relative changes in SST prior to modelling future projections which had been done for only one of the above-mentioned studies (Müller et al., 2009).
Our SDM projected general expansions of suitable habitat in northern areas for most species; but we do not provide calculations of northern (leading) edge shifts as seaweed occurrence records are very limited in the Arctic. Systematic surveys are either >30 years old (Lee, 1980;Wilce, 1959) or lacking entirely (Filbee-Dexter et al., 2019), although recent surveys have documented seaweed compositions along Baffin Island (Küpper et al., 2016) and Greenland (Høgslund et al., 2014). By incorporating as many northern occurrence records as possible, including older and most up-to-date records, we predicted a much greater northern present-day and projected future distribution for fucoids than previous studies (Assis et al., 2014;Jueterbock et al., 2013). As more detailed Arctic occurrence records become available, our SDM and resulting projections could be further improved to provide higher certainty for species' northern edges.
The projected range expansions in the NW-Atlantic are possible as land masses connect temperate and Arctic seas thereby facilitating dispersal (Krause-Jensen & Duarte, 2014). Northward expansions of seaweeds will depend on each species' dispersal ability, as long-distance dispersal is possible via rafting (Fraser et al., 2018;Kalvas & Kautsky, 1998;Olsen et al., 2010;Trowbridge & Todd, 1999) or kelp zoospores (Reed, Laur, & Ebeling, 1988 (Assis, Araújo et al., 2018a;Franco et al., 2018;Jueterbock et al., 2013;Raybaud et al., 2013), the Arctic (Jueterbock, Smolina, Coyer, & Hoarau, 2016), and in Japan (Takao, Kumagai, Yamano, Fujii, & Yamanaka, 2015). Yet future studies are needed for the west coast of North America as well as Africa and South America. In the Southern Hemisphere, extinctions are more likely  as there are no land masses to promote poleward migration to Antarctica, and strong circumpolar currents limit dispersal abilities, but warming and increased storm events may promote long-distance rafting of kelps (Fraser et al., 2018).
Currently, the Gulf of Maine and Atlantic coast of Nova Scotia are already hot spots for climate change, with greater SST increases than other ocean regions (Pershing et al., 2015) leading to decreases in kelp abundance (Filbee-Dexter et al., 2016;Krumhansl et al., 2016;Witman & Lamb, 2018), shifts in fucoid composition (Ugarte et al., 2010) and commercial fish stock collapses (Pershing et al., 2015). The current commercial rockweed harvest in Gulf of Maine and Canadian Maritimes (Seeley & Schlesinger, 2012) may not be sustained with continued warming as our projections suggest a shift from the dominance of A. nodosum with lower tolerance to warmer SST to the more opportunistic F. vesiculosus (Wilson et al., 2015). In the subtidal zone of Maine and Nova Scotia, our projections suggested a shift in dominance from native kelp to invasive C. fragile, potentially allowing C. fragile to settle and form dense beds (Scheibling & Gagnon, 2006 In comparison, projected SST increases, decreases in ice cover, and changing salinity and turbidity patterns in the Arctic will generally have positive effects on the distributional range and production of fucoids and kelps (Filbee-Dexter et al., 2019;Marbà et al., 2017).
Yet the projected northward spread of temperate/boreal species may negatively affect established species such as F. evanescens (Küpper et al., 2016) and the endemic kelp L. solidungula (Filbee-Dexter et al., 2019). Negative interactions following the northward shift of temperate species to polar waters have already resulted in Arctic food web alterations (Kortsch, Primicerio, Fossheim, Dolgov, & Aschan, 2015) and structural changes in fish communities as Arctic species retreat poleward . Lack of baseline data makes it difficult to quantify if any range shifts, negative or positive interactions have already occurred between temperate and Arctic seaweeds and associated communities. Therefore, detailed studies and ongoing monitoring are needed to understand the consequences of climate change in Arctic ecosystems.
Overall, other habitat-forming species have also been projected to experience poleward range shifts, losses in total suitable habitat area and alteration of species richness patterns. SDM of mangroves projected species-specific poleward latitudinal range shifts of 2°, general losses of total suitable habitat, and shifts in species richness globally (Record et al., 2013). For seagrasses, SDM projected 81.4% of Zostera noltii's current range in the NE-Atlantic to remain as suitable habitat by 2100 with a 888 km northward shift of the range centre (Valle et al., 2014). Furthermore, analysis of carbon sources from coastal sediments in Greenland suggest large increases in seagrass meadows over the 21st century (Marbà, Krause-Jensen, Masqué, & Duarte, 2018). For reef-building corals, SDM projected a loss of 24%-50% of present-day habitat, mainly between 5 and 15° latitudes, poleward expansions, and decreases in coral richness by end-century (Cacciapaglia & van Woesik, 2015). This highlights that many coastal ecosystems around the world will experience shifts in habitat-forming species with ripple effects on dependent communities, including humans.