High climate velocity and population fragmentation may constrain climate‐driven range shift of the key habitat former Fucus vesiculosus

The Baltic Sea forms a unique regional sea with its salinity gradient ranging from marine to nearly freshwater conditions. It is one of the most environmentally impacted brackish seas worldwide, and the low biodiversity makes it particularly sensitive to anthropogenic pressures including climate change. We applied a novel combination of models to predict the fate of one of the dominant foundation species in the Baltic Sea, the bladder wrack Fucus vesiculosus.


| INTRODUC TI ON
Owing to the large and densely populated drainage area, the Baltic Sea is one of the most environmentally impacted seas in the world with habitat loss, eutrophication, pollution and overfishing (Diaz & Rosenberg, 2008;Halpern et al., 2008). In addition, recent invasions by non-native species form a strong pressure on the native biota (Ojaveer & Kotta, 2015). Moreover, the Baltic Sea is one of the fastest warming regions in the world, resulting in dramatic retreat of winter ice cover, rise of surface water temperature and reduction of salinity (IPCC 2013;Meier, Kjellström, & Graham, 2006;Meier et al., 2012a). Low functional diversity of the Baltic marine ecosystem exposes it to fast environmental deterioration and potential loss of essential ecosystem services due to low resilience capacity (Meier et al., 2012b;Österblom et al., 2007).
Climate-driven changes in temperature, winter ice cover, salinity and circulation pattern combined with changes in nutrient loading affect patterns of habitat availability. Species may track the displacement in habitat distribution, persist in the altered habitat due to plasticity or adaptation, or go locally extinct (Chevin, Lande, & Mace, 2010;Munday, Warner, Monro, Pandolfi, & Marshall, 2013).
Successful range shifts may depend on dispersal capacity of essential genotypes, but changes in distribution pattern, for example increased fragmentation, may also affect the overall metapopulation persistence (Ovaskainen & Hanski, 2003). Clearly, to track environmental change, populations must shift their range at least as fast as does the environment.
Most successful colonizers of the Baltic Sea consist of marine and freshwater species with a generally broad salinity tolerance, but there is growing evidence for many examples of local adaptations to the Baltic environment (Berg et al., 2015;Johansson, Pereyra, Rafajlovic, & Johannesson, 2017;Momigliano et al., 2017;Pereyra, Bergström, Kautsky, & Johannesson, 2009;Serrão, Kautsky, & Brawley, 1996;Sjöqvist, Godhe, Jonsson, Sundqvist, & Kremp, 2015;Väinölä & Johannesson, 2017). Genetic variation is considered an essential component for persistence in a changing environment where new genetic combinations or the dispersal of locally adapted genotypes may retard the range shift of a species.
Bladder wrack (Fucus vesiculosus) is the only canopy-forming seaweed on hard bottoms in the Baltic Sea, and it constitutes a key habitat with many associated algae, invertebrate and fish species (Kersen, Kotta, Bučas, Kolesova, & Dekere, 2011;Schagerström, Forslund, Kautsky, Pärnoia, & Kotta, 2014;Wikström & Kautsky, 2007). At present, F. vesiculosus extends into the southern part of the Gulf of Bothnia, being less common or partly replaced by its sister species F. radicans further north. Fucus vesiculosus also extends into the inner parts of the Gulf of Finland. Locally, reduced light availability mostly driven by elevated nutrient loading and phytoplankton production may limit the depth range of F. vesiculosus (Torn, Krause-Jensen, & Martin, 2006). Fucus vesiculosus thus has rather well-defined range margins in the Baltic Sea, and predicted climate-driven changes together with ongoing eutrophication and habitat fragmentation may all modify the present distribution of F. vesiculosus (Vuorinen et al., 2015). A loss of this habitat-forming species may have far-reaching consequences and lead to changes in biodiversity and functioning of major shallow ecosystems in the Baltic Sea.
In this study, we used a species distribution model to predict the current and future distribution of F. vesiculosus under the most plausible climate and nutrient loading scenarios. In addition, to explore whether locally adapted genotypes of F. vesiculosus (Ardehed et al., 2016) may track climate change, we applied a biophysical model to simulate dispersal, connectivity and ultimately habitat fragmentation at the receding margin where loss of connectivity may further endanger population persistence (Opdam & Wascher, 2004). The modelled rate of dispersal was compared with estimates of climate velocity, calculated as the migration of the isolines of temperature and salinity across the seascape in kilometre per decade (Burrows et al., 2011), to assess whether predicted spread of tolerant genotypes is sufficient to match the pace of climate change. Finally, we discuss potential ecosystem consequences related to shifting patterns in species distribution and connectivity and suggest some management actions to reduce possible ecosystem impacts.

| The study area-the Baltic Sea
The Baltic Sea was formed from a post-glacial freshwater lake and is only about 8,000 years old. More than 600 river and stream basins drain into the Baltic Sea that opens to the marine North Sea through the shallow Danish straits (Hannerz & Destouni, 2006).
The non-tidal circulation creates a relatively stable salinity gradient through the Baltic Sea from its entrance to the inner parts ( Figure   S1). The latitudinal extent of the Baltic Sea also produces a strong temperature gradient with a significant development of sea ice in the winter-spring season in its northern parts, whereas its southern areas are typically ice-free. A short geological history, low salinity and low temperature compared to nearby coastal waters in the east Atlantic make the Baltic Sea a species-poor ecosystem with key ecosystem functions often provided by single species (Bonsdorff, 2006;Bonsdorff & Pearson, 1999).

| Fucus vesiculosus occurrence
The distribution data of F. vesiculosus were compiled from several sources ( Table 1). The HELCOM dataset covered the whole Baltic Sea with a grid of 10 km × 10 km in which the presence or absence of F. vesiculosus was assessed. Only records from 1995 to the present (2015) were used. The resulting 10-km grid served as the input of current species distribution in the modelling tasks. Recently, a sister species to F. vesiculosus was described as F. radicans (Bergström, Tatarenkov, Johannesson, Jonsson, & Kautsky, 2005), which is endemic to the Baltic Sea. These two species may occur in the same habitat in the northern Baltic Sea (Schagerström et al., 2014) and they can be difficult to separate based on only morphological traits, and we acknowledge that occurrence data on F. vesiculosus may include also F. radicans.

| Environmental layers from regionalized general circulation model
The environmental layers used in the species distribution model (SDM, see below) were produced from simulations of a coupled physical-biogeochemical model of the Baltic Sea (Meier et al., 2012a) based on one climate change scenario (A1B) and two nutrient loading scenarios. The A1B climate scenario is a scenario proposed in the Special Report on Emissions Scenarios by the Intergovernmental Panel on Climate Change (Nakićenović et al., 2000). The A1 scenario group describes a future world of very rapid economic growth, a global population that peaks in mid-century and declines thereafter, and the rapid introduction of new and more efficient technologies.
The two nutrient emission scenarios included "business-as-usual (BAU)" with increased river nutrient concentrations (HELCOM, 2007) and current levels of atmospheric deposition, and a "reference scenario (REF)" where nutrient concentrations in rivers and atmospheric deposition remain at current levels (Eilola, Meier, & Almroth, 2009 Meier, Döscher, & Faxen, 2003) was forced by regionalized atmospheric data from the global climate model ECHAM5 (Roeckner et al., 2006) (Webb, Coward, de Cuevas, & Gwilliam, 1997). For the present simulations, RCO-SCOBI used a horizontal resolution of 3.7 km (2 NM) and with 83 vertical levels with a layer thickness of 3 m. Run-off changes were calculated with a hydrological model for the Baltic catchment area (Lind & Kjellström, 2009). The RCO-SCOBI model produced the following physical and chemical data layers: Secchi depth, salinity, temperature, nitrate (NO 3 ) and phosphate (PO 4

| Calculations of climate velocities
Climate velocity (Burrows et al., 2011) was calculated from the climate scenario layers of projected temperature and salinity as a measure of the speed and direction of change in the Baltic seascape. Climate velocity has been suggested as a useful metric for comparing shifts in species distribution with the rate of climate change (Burrows et al., 2011;Pinsky, Worm, Fogarty, Sarmiento, & Levin, 2013). The grids (3.7 km) of temperature and salinity were first aggregated by calculating the average of 3 × 3 original grid cells. The gradient of scalars E (°C/km or psu/km) between neighbouring grid cells was then calculated using MATLAB (2016a, MathWorks Inc) for the whole domain as follows: where x and y are orthogonal coordinates in space. The climate velocity (km/decade) was then calculated for each grid cell as follows: where change is the local change in temperature or salinity per decade as calculated from the control period  and future layers (2070-2099).

| Species distribution model
The and wind data together with empirically derived algorithms to mimic diffraction (Isaeus, 2004). The average values of all RCO-SCOBI variables (season × layer × scenario combinations) were added to each cell of the species distribution data of 10 km. This dataset served as the input for model calibration and prediction.
The general concept of species distribution modelling (Elith & Leathwick, 2009) was applied to predict the spatial distribution of F. vesiculosus under current and future climate conditions in the Baltic Sea. To investigate which season (summer and winter) and vertical layer (surface and bottom) of RCO-SCOBI data yields the best prediction accuracy, candidate models of all season × layer combinations were calibrated. When building up models, 85% of the input data (4,094 data points) were randomly selected and used as model calibration data, while 15% of the data (722 data points) were reserved for external validation. In the model calibration, the binary distribution data of F. vesiculosus (presence = 1, absence = 0) were used as the dependent variable and the RCO-SCOBI layers together with the depth and wave exposure served as independent predictors. The probability of occurrence of F. vesiculosus was predicted for each of the data points in the validation dataset, and the predictive performance of the candidate models was validated by calculating the area under the receiver operating curve (AUC) (Fielding & Bell, 1997 Random forest (RF) was used as the modelling algorithm. RF is a machine learning method that generates a large number of regression trees, each calibrated on a bootstrap sample of the original data (Breiman, 2001). Each node is split using a subset of randomly selected predictors, and the tree is grown to the largest possible extent without pruning. For predicting the value of a new data point, the data are run through each of the trees in the forest and each tree provides a value. The model prediction is then calculated as the average value over the predictions of all the trees in the forest (Breiman, 2001). The package "randomforest" (Breiman, Cutler, Liaw, & Wiener, 2015) was used to run RF models in the statistical software r 3.2.2 (The R Foundation for Statistical Computing 2015). The parameter mtry (number of variables randomly sampled as candidates at each split) was set to the square root of the number of variables in the model calibration dataset, which is the recommended default value of mtry in the randomforest package (Liaw & Wiener, 2002), and the number of trees was set to 1,000 (Liaw & Wiener, 2002). The importance of predictor variables in the PRE model was assessed using the internal method of the package "randomforest" (mean decrease in accuracy) using 10 permutations (Breiman et al., 2015). Partial dependence plots (sensu Friedman, 2001) were produced using par-tialPlot function in "randomforest" to illustrate the dependence of model predictions on individual covariates. All geographical maps were produced in QGIS or ArcGIS (Esri).

| Biophysical model of dispersal and connectivity
Fucus vesiculosus may disperse in several ways including shortdistance dispersal of gametes and zygotes to occasional longdistance dispersal of dislodged and rafting buoyant sexually mature plants (Rothäusler, Corell, & Jormalainen, 2015;Thiel & Gutow, 2005). Also, dispersal on intermediate scales may occur through bed load transport of negatively buoyant plants and in some cases vegetative, adventitious branches (Johansson et al., 2017;Pereyra et al., 2013). Genetic studies indicate population structure on scales ranging from, in the extreme case 10 m, but more usually 1 km, which suggest relatively short dispersal distances (Muhlin, Engel, Stessel, Weatherbee, & Brawley, 2008;Tatarenkov, Jonsson, Kautsky, & Johannesson, 2007). We estimated the dispersal and connectivity of F. vesiculosus with a biophysical model based on a  (Levitus & Boyer, 1994). The model had a free surface, and the atmospheric forcing was the re-analysis dataset ERA40 (Uppala et al., 2005). Run-off was based on climatological data based on a number of different databases for the Baltic Sea and the North Sea. Validation of the NEMO-Nordic showed that the model is able to correctly represent the SSH, both tidally induced and wind-driven (Hordoir et al., 2015). For each grid cell, the dispersal probability of all other grid cells was calculated from particle positions at the release and end of the simulation. Probabilities were entered into a connectivity matrix where each element represents the probability to disperse from column i to row j (Jonsson, Nilsson Jacobi, & Moksnes, 2016). Elements along the diagonal represent local retention.

| Population model to assess expansion rate and fragmentation
A simplistic metapopulation model for F. vesiculosus was used to explore (1) the rate of spread of locally adapted genotypes and (2) the effect of future fragmentation caused by the predicted changes in habitat distribution and connectivity. The population model con- For the analysis of fragmentation, the population was initiated with an equal number of individuals in each grid cell predicted as habitat from the SDM for the present distribution and the predicted future for the BAU nutrient scenario. The population was projected for 25 generations, spanning 100 years for F. vesiculosus with an assumed generation time of 4 years (Lüning, 1985). We arbitrarily considered a grid cell populated by F. vesiculosus over a range of population sizes covering six orders of magnitudes, for example from 1 individual per km 2 to 1 individual per 1 m 2 , and as absent below this range. As the model of the spread of locally adapted genotypes did not include mortality or Allee effects, the predicted range expansion is likely an upper bound.

| Predicted environment from the RCO-SCOBI scenario model
The coupled physical-biogeochemical model of the Baltic Sea forced with the climate change scenario A1B predicts an increase of summer surface sea temperature of about 3-5°C at the end of this century (Figure 1a). Surface salinity is predicted to decrease with 1.5-2 psu in most of the Baltic Sea (Figure 1b). The change in nutrient loading of nitrogen and phosphorous was evaluated for two scenarios where the business-as-usual scenario produced higher nutrient concentrations in many areas although phosphorous is expected to decrease in the Gulf of Finland ( Figure 2). With higher nutrient loading resulting in increased primary production and plankton biomass, the model predicted a general decrease in water transparency with a reduction in Secchi depth of 1-2 m (Figure 3). However, little change in water transparency is expected at the margins of the distribution of F. vesiculosus.

| Dispersal of F. vesiculosus
The biophysical model showed that a 5-day drift in surface waters is expected to result in dispersal distances typically between 5 and 15 km along the Baltic Sea coast (Figure 7a). The dispersal direction reflects the general circulation pattern with southward dispersal along the Swedish east coast and northward along the Finnish west coast (Figure 7b).

| Metapopulation model of fragmentation and genotype range shift
In a simplistic metapopulation model, we combined the distribution  (Table 2).

| D ISCUSS I ON
The species distribution model (SDM) based on a regional climate model predicts a dramatic 30% reduction of the inhabitable area for the habitat-forming seaweed F. vesiculosus in the Baltic Sea. The mechanism behind this range shrinkage is a dramatic southward and westward expansion of low-salinity water, most notably in the Gulf of Bothnia and the Gulf of Finland. Although the projected range of salinity decline is uncertain, most models in an ensemble simulation of 16 scenarios showed a decrease in salinity (Meier et al., 2006). In addition to the dilution effect, an increasing nutrient load may enhance plankton primary production leading to reduced water transparency as well as trigger blooms of epiphytic filamentous algae (Valiela et al., 1997). In either way, F. vesiculosus populations could be weakened dramatically. We further show that the loss of suitable habitat may result in reduced connectivity leading to lower metapopulation persistence and to further fragmentation, especially at the distribution margins where significant losses of local populations F I G U R E 3 Projected future change in water transparency expressed as the Secchi depth, assuming a business-as-usual scenario for the nutrient loading Recent studies have pointed out that predictions of SDMs may both under-and overestimate future habitat loss if the population consists of locally adapted subpopulations (e.g., Hällfors et al., 2016;Oney, Reineking, O'Neill, & Kreyling, 2013). Predictions may be overly optimistic if locally adapted genotypes are dispersal-limited. This may apply to F. vesiculosus where poor dispersal ability with zygotes or clonal adventitious branches is generally expected (Johansson et al., 2017;Pereyra et al., 2013;Väinölä & Johannesson, 2017). Our biophysical model with drift duration of 5 days predicted relatively short dispersal distances, generally less than 10 km. Dislodged and drifting adult thalli may travel for longer distances and release gametes far from their native area (Rothäusler et al., 2015), which is also indicated for other fucoids (Buonomo et al., 2017). However, dislodged plants may only rarely lead to realized dispersal as spawning is often constrained by calm weather (Serrão et al., 1996), the moon phase (Andersson, Kautsky, & Kalvas, 1994) and the risk of beaching (Muhlin et al., 2008). At present, it is not possible to estimate the significance of drifting thalli to overall gene flow. Moreover, the oceanographical circulation models (RCO-SCOBI and NEMO-Nordic) cannot accurately simulate currents close to the coast due to their relatively coarse spatial resolution. The values in the coastal grid cells used in this study are hence likely an overestimate of actual coastal current velocities. In the current paper's context, the modelled dispersal should then be viewed as an upper bound and F. vesiculosus may be even more dispersal-limited than the model predicted.
Studies of genetic differentiation and morphology suggest that the population structure of F. vesiculosus in the Baltic Sea may be very complex with a mix of sexual and asexual recruitment (Ardehed et al., 2016). In addition, a closely related species F. radicans has evolved as an endemic species from F. vesiculosus inside the Baltic Sea (Ardehed et al., 2016;Bergström et al., 2005;Pereyra et al., 2013). Although genomic information about geographical differences in adaptive loci F I G U R E 5 Species distribution model (SDM) of the change in the distribution of Fucus vesiculosus from a control period  to the projected future climate period . (a) Modelled present distribution, (b) modelled distribution assuming a future continuation of present nutrient loading and (c) modelled distribution assuming a business-as-usual scenario for the nutrient loading F I G U R E 6 Importance of predictor variables in the final random forest model. The x-axis shows the decrease in accuracy as the increase of mean-squared error of the model when a variable is permuted is lacking, there are indications of trait differences that suggest the presence of locally adapted populations along the Baltic Sea gradient in both these species (Johansson et al., 2017;Pearson, Kautsky, & Serrão, 2000;Serrão et al., 1996).
If locally adapted genotypes will successfully track an advancing or receding optimal environment during a climate change scenario, the dispersal to and colonization of new areas must occur at a rate similar to the spatial change of critical environmental factors.
Failing to track a receding critical environment boundary, for example salinity, may lead to local or even global extinction of unique genetic variants and adaptations. The calculated climate velocities for salinity in this study indicate that the dispersal rate is lower or of the same order as the rate of range shift for F. vesiculosus (Table 2; Figure 4). Only if long-distance dispersal of drifting thalli is an important contribution to gene flow, are locally adapted genotypes expected to outpace the receding boundary of critical salinity. This also assumes relatively rapid selective sweeps of locally beneficial alleles (Morjan & Rieseberg, 2004). In addition, García Molinos, Burrows, and Poloczanska (2017) recently showed that many ongoing range shifts are correlated not only with the climate velocity but also with the directional agreement between climate velocity and dispersal, especially for species with propagules strongly influenced by physical water transport. In a qualitative comparison between the direction of decreasing salinity ( Figure S3) and mean dispersal direction (Figure 7b), there is a match between directions of climate velocity of salinity and dispersal along the Swedish coast and within the Gulf of Finland, but generally a mismatch along the Finnish and Estonian coasts. A mismatch for the sessile F. vesiculosus should further obstruct a range shift of locally adapted genotypes. The increasing fragmentation of the marginal populations may also reduce successful dispersal and increase local extinction as predicted by the metapopulation model ( Figure 8). Consequently, there is relatively high risk that locally adapted populations, evolved to tolerate low salinity in the Bothnian Sea and the Gulf of Finland, may go extinct under current IPCC climate scenarios. Such inability of adapted genotypes to track the spatial change in the environment may lead to an even more severe decline of F. vesiculosus in the future Baltic Sea than predicted by the SDM as there will be no locally tolerant genotypes when the critical salinity boundary sweeps south-west. This may lead to a reduced tolerance to low salinity at the population level with F. vesiculosus requiring higher salinity for persistence than today. This effect may be long-lasting unless evolution of new locally adapted populations is rapid. But the latter is not likely considering the long generation time and the potentially small size of fragmented populations of F. vesiculosus. There is also a need for more detailed information about the relative contribution of plasticity (Johansson A dramatic decline of F. vesiculosus due to reduced salinity will likely be accompanied by a severe structural and functional regime shift of the coastal ecosystem including a dramatic loss of biodiversity in the Baltic Sea (Kotta, Möller, Orav-Kotta, & Pärnoia, 2014;Vuorinen et al., 2015). Much of brackish and marine vegetation will be replaced by freshwater higher plants . However, none of the freshwater habitat-building species can replace F. vesiculosus on rocky bottoms. Consequently, such habitats will likely lose much of their extant vegetation, become less complex and speciespoor and be deprived of critical ecosystem services such as nutrient uptake, food source, refuge for juvenile fish and loss of aesthetic values for recreational activities. It is possible that if dispersal and colonization are sufficiently fast, the sister species F. radicans may replace F. vesiculosus in areas where salinity decreases over time, which could maintain canopies on rocky substrates. Further studies are here needed.
Conservation actions may mitigate the effect of low dispersal of locally adapted populations. Assisted dispersal and colonization through translocation of individuals of Fucus or zygotes is one option (Hiddink, Ben Rais Lasram, Cantrill, & Davies, 2012;Thomas, 2011).
It may be even possible to produce zygotes or seedlings through assisted spawning and fertilization (Serrão, Brawley, Hedman, Kautsky, & Samuelson, 1999). For some populations, as in the northern Baltic Sea with high rates of asexual reproduction (Ardehed et al., 2016), it may be possible to translocate clonal copies using adventitious branches of tolerant clones (Johannesson, Smolarz, Grahn, & André, 2011). Ultimately, even evolution could be assisted where genotypes are selected that are more tolerant to a rapidly changing environment (van Oppen, Oliver, Putnam, & Gates, 2015). Artificial selection may be a challenge for F. vesiculosus considering its long generation time of over 4 years.
F I G U R E 8 Distribution of Fucus vesiculosus produced by the metapopulation model (eq. 3) where habitat prediction from the species distribution model is combined with the connectivity estimated with a biophysical model. (a) Population density when metapopulation model is based on the present habitat prediction in Figure 6a and (b) population density when metapopulation model is based on the future habitat prediction in Figure 6c. Population density is given in relative units but was modelled as number of thalli per km 2 TA B L E 2 Rate of population expansion from three selected locations in the Gulf of Bothnia and the Gulf of Finland Rates are shown for two dispersal strategies with drift duration of either 5 or 30 days in surface waters (0-2 m) as a mean for the period May-August.