Mechanistic niche modelling to identify favorable growth sites of temperate macroalgae

The European seaweed cultivation sector is in a transition phase with the rise of seaweed aquaculture due to an increased interest in seaweed resources. Identifying regions with optimal growth conditions for the cultivation of specific seaweed species contributes to the cultivation process. An understanding how these regions evolve under climate change is required to ensure favorable growth conditions on the long-term. In the present research, regions with favorable growth conditions for specific seaweed species were identified by combining physiological and environmental data in a mechanistic niche model. The outcome of the mechanistic model is a speciesspecific response, the habitat suitability, which quantifies growth as a function of the temperature, salinity, light and nutrient requirements of the seaweed species. Habitat suitability was quantified in European marine waters for brown seaweeds (i.e. Fucus serratus, F. vesiculosus, Ascophyllum nodosum, Saccharina latissima, Laminaria digitata, Laminaria hyperborea) and red seaweeds (i.e. Chondrus crispus, Gracilaria gracilis, Furcellaria lumbricalis). The model was validated using independent distribution data and the validation statistic was good (area under the curve; AUC > 0.8) for five out of nine species and fair (0.6 < AUC < 0.8) for the remaining four species. The warm-temperate region extending from the coast of Portugal to the south coast of Brittany is currently a suitable habitat for most of the studied species. Due to climate change, we predict that the most optimal environmental conditions will shift northwards, i.e. 110 to 635 km by 2100, depending on the climate scenario. The results of the present study can be used to: (1) select target species for seaweed aquaculture in a specific marine region; (2) select sites for long-term, optimal growth conditions of the specific seaweed species in this study. As such, our results contribute to the decision making process in marine spatial planning and blue growth prior-


Introduction
Aquaculture, involving the farming of fish, crustaceans, molluscs and seaweeds, is the fastest growing sector of global food production and one of the main themes of the European Union's Blue Growth Strategy [1,2]. Production of seaweeds for human food, pharmaceuticals, fertilizers, fodder and chemicals has doubled over the past ten years and, considering biomass, seaweeds are among the most important cultivated marine organisms [3,4]. Apart from the yield, the cultivation of seaweeds provides an environmental service by the seaweed's bioremediation capacity. This bioremediation capacity can be used to mitigate the main environmental challenge of animal aquaculture, its nutrient discharge in surface waters. By integrating seaweed growth and animal aquaculture, up to 90% of the nutrient discharge could be remediated from proximate surface waters, thereby lowering the risk of eutrophication [3,5]. In addition to the environmental services, the aquaculture of seaweeds is also profitable. The global, annual yield of seaweed cultivation is 28 million tonnes wet weight with an estimated value of USD 5 to 6 billion of which ca. 50% is destined for human consumption [4,6]. The European contribution to the global yield, however, is < 1% of the biomass [4].
European seaweed production is mainly based on the harvesting of wild stocks but the increasing market interest in seaweed resources and the need to assure the environmental sustainability of the harvesting process has stimulated the development of the seaweed aquaculture sector in Europe [7,8]. With seaweed cultivation being historically mainly an Asian tradition and relatively new to European waters, identifying regions with favorable growth conditions for specific seaweed species is needed for an optimal yield and could contribute to the further development of the European sector. Moreover, knowledge about the most optimal growth sites is indispensable to underpin marine spatial planning processes and marine policy. Cultivation of specific seaweed species at sub-or supraoptimal temperatures will result in a reduced yield [9]. For example cultivating specific populations of the red seaweed Porphyra umbilicalis at 5°C instead of its temperature optimum of 10°C reduces the daily growth rate by 50% [9]. For some seaweed species, favorable regions for growth have already been identified by inferring a species' abiotic niche from its distribution. However, optimal growth of a seaweed species may occur at a specific combination of environmental conditions that in the present distribution does not occur [10].
Predicting suitable habitat for specific seaweed species can be achieved by combining physiological and environmental data in a mechanistic niche model [11]. This habitat suitability is a quantification of the environmental restrictions affecting a species' physiology and, consequently, growth [12]. Temperature ranges for growth and survival have been determined experimentally for a large number of seaweed species during the last four decades [9,[13][14][15][16]. These temperature responses, combined with responses to other abiotic variables, can be integrated in a mechanistic model. Apart from predicting habitat suitability under present conditions, mechanistic models may also forecast range shifts and contractions under future climatic conditions [12]. At many temperate European coastlines, range shifts of seaweed species are already taking place [17,18]. Simulated future changes of seawater temperatures predict a poleward movement of biogeographic regions with a shrinking of the cold-temperate region at the southern border and a northward expansion at the expense of the polar region [19], the magnitude of the shift and the migration rate depending on the species [20,21].
In the present study, the species-specific habitat suitability of nine temperate macroalgae is quantified at a European scale by means of mechanistic species distribution modelling (SDM). This mechanistic SDM predicts the distribution of a species by using its physiological response to environmental conditions. The model might aid in the process of site selection for optimizing aquaculture activities, thereby stimulating blue growth. The tool may also support dynamic marine spatial planning by including climate change to allow the selection of sites that in the long-term will still have favorable conditions for seaweed growth. The shift in a species' distribution is quantified according to four climate change scenarios; the Representative Concentration Pathways, adopted by the Intergovernmental Panel on Climate Change [22]. The model's outcome is validated using independent distribution data to assess the predictive power in function of its sensitivity and specificity.

Materials & methods
To quantify the habitat suitability of economically important macroalgae in the European marine environment, we made use of mechanistic models of which the outcome will be validated using current, species-specific distribution data. Hence, information on the abiotic conditions that influence seaweed growth, information on species-specific physiological responses as well as distributional data were required. Species included in this study are brown seaweeds (i.e. Fucus serratus, F. vesiculosus, Ascophyllum nodosum, Saccharina latissima, Laminaria digitata, Laminaria hyperborea) and red seaweeds (i.e. Chondrus crispus, Gracilaria gracilis, Furcellaria lumbricalis). Species were selected based on their potential economic importance and the availability of physiological data, using only seaweeds that are distributed throughout European shores.

Model & environmental data
For each species, the habitat suitability is quantified using a mechanistic model that integrates a set of abiotic variables that may limit macroalgal growth. Based on literature review, environmental layers on sea surface temperature, sea surface salinity, photosynthetically activate radiation (PAR), diffuse attenuation (DA), nitrate concentration, depth and temperature of the warmest month were included ( Table 1). All of these variables were obtained through Bio-ORACLE version 2 [23] using the sdmpredictors [24] package in R Studio [25], except for the mean temperature of the warmest month, which was obtained through MARSPEC [26]. The effect of these variables on macroalgal growth was deduced from physiological experiments available in literature. As such for each variable, a species-specific response that quantifies growth as a function of the variable was defined [27]. The growth rate is hereby defined as the ratio of mass change to the time elapsed during the experiment. This relative growth rate is measured without drying the specimens, however, surface moisture is removed. This growth-based, variable-specific response ranges from zero (not suitable) to one (highly suitable). The combined effect of all potentially limiting variables was calculated by multiplying the response of all variables. Responses were mapped in geographical space (−20°to 20°l ongitude, 40°to 70°latitude; Fig. 1) at a spatial resolution of 5 arcminutes (i.e. 9.2 km) using the R package raster [28]. Correlation between environmental variables was tested using the sdmpredictors package [24], setting the threshold for autocorrelation at 0.80 [27].
Seaweed cultivation is not necessarily restricted to coastal regions and can be, for example, integrated in coastal protection or offshore wind farms [67]. Therefore, and in addition to physiological limitations, spatial planning of the marine environment and depth of the water column was also taken into account as a practical and legislative limitation for deploying offshore aquaculture infrastructure. In this context, Lee et al. [30] reported that floating, offshore infrastructure can easily anchor between a depth of 0-70 m, but not deeper than 700 m. Hence, to take this into account, an anchoring depth of 0-70 m was modelled to be optimal and increasingly challenging towards the maximum depth of 700 m by means of a linear relation.
The anchoring depth of the cultivation infrastructure is independent from the cultivation depth, the latter common to be 2 m below the water surface [33]. To account for marine spatial planning, OSPAR's marine protected areas (MPAs) were added as an example of marine Table 1 Variables used in the mechanistic model with their units, response functions and an example of a physiological experiment. Linear response functions represent two or more curves that were fitted to physiological data to cover the sub-and supraoptimal conditions. Habitat suitability (hs) is defined as the product of hs sst , hs sss , hs par & hs nitrate . 'x' stands for the species-specific slope and 'c' for the species-specific constant of the linear response (x sst + c sst and x sss + c sss for the temperature and salinity response, respectively).  [31], therefore four out of six categories were included in this layer. To account for the use of artificial substrates, as part of the floating cultivation infrastructure, substrate preferences were not included in the model. The physiological derived limitations were translated in response functions ( Table 1). The growth responses to temperature and salinity were considered negatively skewed with a steeper slope for supraoptimal conditions compared to suboptimal conditions. Based on the physiological experiments, two or more curves were fitted to cover the sub-and supraoptimal conditions. To do so, the species-specific growth rate was normalized (growth rate within range 0-1) and expressed in function of temperature or salinity, as described in Eggert [32]. For example, Bird et al. [13] and Fortes & Lüning [9] quantified growth of C. crispus over a period of 14 days at a temperature gradient ranging from 0 to 25°C while keeping salinity, irradiance and period of illumination constant. For this species they found maximal growth at 15°C. Hence, a temperature of 15°C would result in a temperature response of 1, while a temperature of 17°C would result, based on the linear relation derived from physiological experiments, in a habitat suitability of 0.8. Annual averages were used for the temperature and salinity metrics ( Table 1). The response to DA was modelled as a reduction in PAR using Lambert-Beer's law: with PAR as irradiance (unit: mol photons m −2 s −1 ) at depth d (m) and x (m −1 ) as the DA coefficient. The response to PAR and nitrate was modelled as a functional response type II, also referred to as a Monod response [27]. This functional response for PAR was predicted by the formula: with μ as the predicted growth rate limited due to suboptimal light conditions and μ max as the maximum growth rate. The parameter ks PAR is the species-specific, half-saturation constant deduced from the physiological experiments. For the functional response to nitrate limitation, ks DIN was used as the rate-limiting parameter. Habitat suitability was assumed to be zero if the mean temperature of the warmest month was equal to, or exceeded the lethal temperature of a species. This factor was added to account for high summer temperatures, especially relevant for coastal areas, that could cause seaweed degradation. DA accounts for the turbidity of the water with a DA of 0 m −1 for completely clear waters and a DA of 1 m −1 for extremely turbid waters. A depth of 2 m in Lambert-Beer's law was maintained throughout the study as a common depth for floating cultivation infrastructure [33]. For example, an irradiance at the sea surface of 300 mol photons m −2 day −1 is at a depth of 2 m reduced to 245 mol m −2 day −1 with a common turbidity in the open ocean of 0.1 m −1 [34]. In more turbid waters (e.g. DA > 0.5 m −1 ), like sandy coasts or estuaries, the irradiance at a depth of 2 m decreases to 110 mol m −2 day −1 .
The relative contribution of a specific variable on the seaweed's growth was quantified as: with hs var as the response of a variable averaged over the study area and ∑(1-hs var ) as the response summed over all variables. For example, 1hs SST (1mean temperature response) for Laminaria hyperborea is predicted to be 0.36 and ∑(1-hs SST ) to be 0.67, resulting in a relative contribution of 0.53 ( 0.36 / 0.67 ). The latter can be interpreted as the contribution of surface temperature in affecting the habitat suitability of L. hyperborea throughout the study area.

Validation
The model's discriminatory power between presence and absence was assessed using the area under the curve (AUC) statistic of a receiver operating characteristic (ROC). Independent presence data of all nine species was obtained through the Ocean Biogeographic Information System (OBIS; www.iobis.org) and the Global Biodiversity Information Facility (GBIF; www.gbif.org) by using the R packages robis and rgbif. Questionable records were removed by setting a maximum depth of 50 m and excluding coordinates on land using the R package obistools [35]. Records that fell outside the credible range of a species, based on expert knowledge or literature review, were also filtered out. To reduce a bias towards sites with a high sampling density, the records were thinned using the R package spThin, setting the thinning distance to 50 km (as in Bosch et al. [36]). Finally, a total of 793 records were retained with on average 88 records ± 29 (mean ± sd) per species. Background points (pseudoabsences) were randomly selected from the entire study area, after applying a bathymetric mask of 50 m depth, using the R package dismo [37]. Ten thousand pseudoabsences were selected for each species, excluding grid cells that contain a species' presence. AUC values were calculated using the R package ROCR [38] and classified according to Swets [39]. However, because an AUC index does not account for the model's goodness-of-fit [40], testing for overfitting was done by comparing the classification of the occurrence records with the classification of the coastal grid cells (grid cell's distance to coast < 10 km). Therefore, the coastal cells were classified as optimal (hs ≥ 0.75), suboptimal (0.5 < hs < 0.75) and poor (hs ≤ 0.5). The distribution of habitat suitability over the coastal cells was subsequently compared with the distribution of the occurrence records. Additional performance metrics that were used are the true positive rate (TPR), which is the relative amount of correctly classified presences, the true negative rate (TNR), which is defined as the relative amount of correctly classified absences, and the true skill statistic (TSS). The true skill statistic is defined as TPR + TNR -1 and ranges from −1 to +1, where +1 indicates a perfect model and values of zero or less indicate a performance worse than random [41]. The TPR, TNR and TSS are obtained after applying a threshold on the continuous habitat suitability. The threshold was determined by selecting for a maximum TSS from all possible thresholds (range 0-1) [42]. Note that the AUC statistic does not rely on a threshold.

Simulations
The Representative Concentration Pathways (RCP) were used for modelling future climate predictions by the year 2100 [43]. Climate scenarios RCP 2.6, RCP 4.5, RCP 6.0 & RCP 8.5 were implemented, with predicted CO 2 concentrations of 421 ppm, 538 ppm, 670 ppm & 936 ppm by 2100, respectively. The present atmospheric CO 2 concentration is 410 ppm (d.d. May 2018; Earth System Research Laboratorium; www.esrl.noaa.gov). The temperature rise for the study area by 2100 is 0.6°C, 0.9°C, 1.4°C and 2.7°C according to listed RCP scenarios, respectively. RCP 2.6 is a peak-and-decline scenario that assumes a maximum emission of greenhouse gasses before 2020. According to RCP 4.5, emissions peak around 2040 and steadily decline thereafter. RCP 6.0 is more pessimistic and the emission peak is predicted to occur around 2080. According to RCP 8.5, emissions will continue to increase and will not level before 2100 [43]. Besides simulating rising temperatures under climate change, the RCP scenario's also include simulations for salinity changes due to alterations in the hydrological cycle [44]. The scenarios for temperature and salinity were obtained through the Bio-ORACLE website (www.bio-oracle.org). However, for the variables nitrate, DA and PAR, future layers were not available and therefore layers with present conditions were used.
The change in habitat suitability under global warming was quantified by implementing the future environmental layers for each respective climate scenario in the mechanistic model and, subsequently, to identify how present favorable regions for seaweed aquaculture evolve in a changing climate. To quantify the potential northward range shift, habitat suitability was quantified under present and future conditions, along a latitudinal gradient from 30°N to 80°N and averaged for each degree of longitude (from −20.0°to 20.0°by 2.5°). A non-parametric Kruskal-Wallis test was used to assess the shift of the habitat suitability between present and future climate conditions. A one-sided Wilcoxon rank-sum served as post-hoc testing to detect pairwise differences.

Results
The distribution of habitat suitability was visualized for all nine species, here depicted are the results of C. crispus as an example. Figures for all nine species are shown in Figs. S1-S9.

Spatial distribution
Overlaying the habitat suitability map solely based on abiotic restrictions ( Fig. 1A) with the non-environmental layers (depth & marine spatial planning), leaves only a fraction (6.2%) of the optimal habitat (Fig. 1B). Suitable environmental conditions for the cultivation of C. crispus are predicted at the northern part of the Bay of Biscay and the north-west coast of Spain. After taking depth and marine spatial planning into account, the south coast of Brittany is predicted as a favorable region for cultivation of C. crispus (Fig. 1B). The relatively large suitable area in Fig. 1A is caused by the broad temperature optimum of C. crispus (between 10 and 15°C and only a slight decrease in growth towards 20°C; [9]). The impact of climate change on habitat suitability (Fig. 1C) suggests a northward range shift, mostly driven by a temperature rise of ca. 2.7°C by 2100 (according to RCP 8.5; see further). The Bristol Channel and the Southern Bight of the North Sea are predicted to be potential regions with present and future favorable conditions for the cultivation of C. crispus. The warm-temperate region extending from the coast of Portugal to the center of the Bay of Biscay is predicted to evolve from a region with optimal habitat conditions (hs > 0.75) to a region with poor habitat conditions (hs < 0.50) with local lethal temperatures. Conditions along the coast of the North Sea are predicted to improve for cultivation of this species, according to the RCP 8.5 climate scenario (Fig. 1C).

Variable influence on habitat suitability
According to the mechanistic model, temperature is the most important variable in determining habitat suitability with a mean relative contribution of 0.45 ( Table 2). The combined effect of temperature and salinity was 0.76 ± 0.04 on average and never below 0.70. Temperature and salinity are almost equally influential in the models of F. vesiculosus and F. lumbricalis because the growth optima of these seaweeds is between 20 and 30 psu [9,45] while the average salinity of the North Sea is 34-35 psu [46]. These optima are lower compared to less euryhaline species such as C. crispus or S. latissima with an optimum around 30-35 psu [9,47]. Photosynthetically active radiation (PAR), diffuse attenuation (DA) and nitrate influenced habitat suitability one a more local scale, therefore their relative contribution was rather low compared to the impact of temperature and salinity. For example, habitat suitability in the Wadden Sea and the center of the North Sea was reduced due to turbid waters (DA > 0.5 m −1 ) and due to low nitrate concentrations (< 1.0 μmol ml −1 ), respectively.

Validation with distribution records
Validation of the mechanistic model shows a good fit of the quality- Table 2 Variable's impact on habitat suitability, quantified as relative contribution (Eq. (3)). For each species, temperature has the highest relative contribution to the habitat suitability. checked occurrence records on the predicted habitat suitability (Fig. 2). Of the 793 distribution records, pooled for all species, 392 were predicted to be in an optimal habitat (hs ≥ 0.75), 323 records in a suboptimal habitat (0.5 < hs < 0.75) and only 78 in a poor habitat (hs ≤ 0.50). The habitat quality of these occurrence records located in euhaline waters, an annual sea surface temperature in the range of 10-15°C and no limitation due to unfavorable nitrate or light conditions, was generally predicted to be optimal. Model performance varied between species with a minimum average predicted suitability of 0.60 ± 0.13 for the distribution records of F. lumbricalis and a maximum prediction for the records of S. latissima (hs 0.80 ± 0.10). The habitat quality of occurrences in the Baltic Sea was generally predicted to be poor due to low salinity levels in this marine region (mean salinity 7.26 psu ± 1.7) while the growth optimum for all macroalgae included in this study is between 20 and 30 psu [13,48,49]. The suitability of localities above the Arctic Circle was predicted to be poor due to suboptimal temperatures (mean temperature 5.8°C ± 1.3°C) compared to the growth optimum of 15°C of most temperate macroalgae [9]. For example, of the 100 distribution records of F. serratus ( Fig. 2A), 89 records were predicted to be in an optimal (68 records) or suboptimal habitat (21 records). The remaining 11 records were predicted to be in a poor habitat (hs < 0.5) mainly due to the suboptimal temperatures in the Arctic region (5 records) or due to the low salinity in the Baltic (6 records). Optimal habitat of this brown seaweed was mainly predicted to be in the North Sea or on the shores of the United Kingdom and Brittany.
Validation of the mechanistic model with the ROC shows that the prediction of the mechanistic model is accurate (Fig. 3). The AUC value for five out of nine species was good (i.e. > 0.8) and fair (i.e. 0.6 < 0.8) for the remaining four, (according to Swets [39]; Table 3). The steeper the ROC curve, the better the relation between true positive rate and false positive rate and the higher the AUC. The AUC statistic correlates with the average habitat suitability of a species' distribution records (Pearson's r = 0.79, p = 0.01). All these performance metrics indicate a reliable model with a slight variability (i.e. AUC min 0.70, AUC max 0.86) in predictive power between species. A 'fair' instead of a 'good' AUC statistic is primarily caused by a relatively low true negative rate (Table 3). This low true negative rate is caused by background samples (pseudoabsences) located in a suitable habitat and are therefore classified as false positive. For example, 3200 background samples of F. vesiculosus were located in an suitable habitat for this species and were therefore classified as false positives. TSS metrics correlated with the AUC values (Pearson's r = 0.93, p = 0.0003) and were all well above zero. The threshold on habitat suitability for presence/absence was set at 0.55.

Northward shift
The southern boundary of C. crispus is predicted at 42°latitude under present conditions but due to global warming, a northward shift of the species' distribution is expected (Fig. 4). Depending on the climate scenario, the pole-ward shift ranged from 157 km (RCP 2.6) to 635 km (RCP 8.5). The intermediate climate scenarios 4.5 & 6.0 predicted conditions causing a pole-ward shift of 230 km and 350 km, respectively. Based on pairwise testing, with present conditions as a reference, a northward shift of C. crispus' distribution was found for all climate scenarios (p < 0.001). The shift was quantified to be 300 km ± 27 km on average for brown seaweeds and 319 km ± 29 km for red seaweeds, under the RCP 6.0 scenario (Table 4).

Northward shift
The potential northward shift of nine temperate seaweeds was for the first time quantified by means of mechanistic modelling. The warmtemperate region extending from the coast of Portugal to the south coast of Brittany is currently a suitable habitat for most studied species but is predicted to evolve to a poor habitat suitability or to an environment with lethal conditions, depending on the climate scenario. This prediction agrees with the model of Müller et al. [19] who related seaweed distribution and temperature requirements to changing sea surface isotherms and found a northwards migration of cold-temperate, macroalgal populations from the coast of the Iberian Peninsula towards northern Brittany. The distribution of the studied seaweeds is predicted to shift 110-635 km north by 2100 according to the respective RCP scenarios (Table 4). This shift is in the same order of magnitude as described by Southward, Hawkins, & Burrows [50], who predicted a shift of 300-600 km north of intertidal, benthic species based on a temperature rise of 2°C and a northward migration of the SST isotherm. Lima & Wethey [51] studied past range shifts between 1955 and 2004, driven by a temperature rise of 0.74°C along the Portuguese coast, and found for seaweeds an average migration of 4.5 km per year (metaanalysis by Sorte, Williams, & Carlton [20]). Although the studies mentioned in this paragraph used correlative modelling to assess the impact of climate change on a species' distribution, the results of this study are similar regarding the magnitude of the pole-ward shift of temperate seaweeds, induced by global warming.
An assumption of the quantification of the northward shift is that the adaptive power of macroalgal species to the changing environment is limited or absent. However, mechanisms for adaptation to changing temperatures are physiological regulation (timescale seconds to Fig. 2. Predicted habitat suitability combined with distribution records. The model performs well on presences in cold-temperate, euhaline waters. The habitat suitability of occurrences in the Baltic Sea and more polar waters is predicted to be poor due to salinity and temperature conditions below the growth optimum of the species, respectively. The distribution of habitat suitability of F. serratus restricted to coastal grid cells was primarily poor (50%), followed by suboptimal (40%) and optimal (10%). Although just 10% of the coastal area was predicted to be an optimal habitat, 68% of the distribution records of this brown seaweed species were located in these grid cells. For C. crispus (panel B), only 8 out of 98 presences were predicted to be in a poor habitat. minutes), phenotypic acclimation (hours to days) and genetic adaptation (evolutionary timescale) [32]. Studies specifically on temperate seaweeds confirm the low adaptive potential to the warmer surface waters and describe a northward migration [18,52,53]. The variance in the magnitude of the northward shift among the species in this study is small (Table 4) because all nine species have similar temperature responses [13,14]. Indeed, all nine species are temperate macroalgae with a growth optimum between 10 and 15°C and eurythermal characteristics due to the large temperature seasonality in their distribution [32,54].

Mechanistic niche modelling
Previous studies on seaweeds' range-shifts predict the climatechange induced migration mainly by means of correlative niche modelling (e.g. Takao et al. [55]; Jueterbock et al. [56]) while the prediction of this study is based on the seaweed's physiological response. Regression-based, correlative models infer a species' ecological niche from its distribution. However, a species may be adapted to a combination of environmental conditions that in the present distribution does not occur. If new combinations of these conditions arise in the future, a statistical model may erroneously classify such environments as unsuitable [10]. Mechanistic models are independent of present climatic conditions because their parameters are not trained on the (present) distribution of a species. With a mechanistic approach, the distribution is defined by the environmental constraints acting on a species' physiology and is therefore considered useful in modelling future distributions under climate change [12]. The approach of mechanistic niche modelling also allowed to develop a model assigning an equal importance to coastal and offshore conditions, and to avoid training the model on distributions restricted to intertidal conditions. For all species, temperature has the highest relative contribution to the mechanistic model's prediction (i.e. > 0.38; Table 2). This is conform the general perception that temperature is the main driver in a seaweed's geographic distribution [57,58]. Temperature and salinity combined, determined over 75% of the habitat suitability. Habitat suitability for the Baltic Sea proved difficult to predict due to its low salinity. For example, Fucus vesiculosus occurs in the Baltic Sea while the habitat suitability for this area was classified as poor. A probable reason is that Atlantic populations of F. vesiculosus are not adapted to grow in the brackish water (salinities: 0-10 psu) of the Baltic Sea [48]. Physiological responses integrated in this model were obtained from Atlantic populations. Habitat quality of the Baltic Sea is classified as poor because intraspecific variation, in this case for the salinity response, is not included in the ecological model. Thus, if a species' salinity response is parametrized by a physiological experiment on Atlantic populations, the habitat quality of e.g. the Baltic Sea may be predicted to be poor. Overall, the reason why Baltic populations are able to grow in strongly reduced salinities is due to local adaptation of these populations. If intraspecific variation would be incorporated in the model, for example by including the physiological responses of multiple strains, the habitat suitability of regions such as the Baltic Sea is likely to increase and consequently the predicted distribution is also likely to expand. This underlines the importance of strain selection in harmonizing the seaweed's physiological requirements and environmental conditions at the cultivation site in order to optimize the cultivation activities. In the process of optimizing seaweed cultivation, it should be taken into account that biotic interactions were not included in the ecological model. Therefore, it is possible that the abiotic conditions are predicted to be favorable while the biotic conditions (e.g. competition, grazing pressure, parasitism, pests) potentially restrict seaweed growth. Accounting for biotic interactions in ecological models is challenging due to the complexity of the interaction and requires a species-specific parametrization of this interaction and the full distribution of the associated biota [12,59]. Although biotic interactions  Table 3 Validation statistics including area under the curve (AUC), true positive rate (TPR), true negative rate (TNR) and true skill statistic (TSS). AUC values were classified as good (a) , fair (b) , poor (c) and fail (d) , as described in Swets [39]. influence seaweed growth on a more local scale [12], it has the potential to severely affect the success of the cultivation activities [3,7]. Another restriction for mechanistic modelling is a potential lack of physiological data [12]. For example, under the RCP 6.0 & RCP 8.5 scenario's ocean surface pH is expected to drop from 8.1 to 7.8 by the end of the century due to anthropogenic CO 2 emission [22], consequently affecting intertidal communities [60]. Studies on the physiology of temperate macroalgae (Ecklonia radiata & Macrocystis pyrifera) suggest an unaffected growth under this projected acidification of 0.3 pH units [61,62]. To date it remains uncertain which species are most vulnerable to such pH drop [60]. Physiological data on the growth rate reduction of specific seaweeds in this study due to ocean acidification would have allowed to parametrize this change. Next to the intraspecific variation that is not accounted for, another restriction of the species distribution model is the spatial resolution of the environmental data. The spatial resolution of marine data layers (Bio-ORACLE: 9.2 km 2 ; [23]) is much smaller than the resolution of terrestrial layers (WorldClim: 1 km 2 ; [63]). Therefore, it is possible that a region that appears to have favorable physical conditions for seaweed cultivation contains microhabitats that are less suitable for seaweed farming due to unfavorable local abiotic or biotic conditions [64]. With higher resolution environmental layers it is not only more likely to detect these microhabitats but it would also allow to work on a more regional scale, for example the Southern Bight of the North Sea, to predict where aquaculture activities will be most beneficial.

Conclusion
We assessed the habitat suitability of nine temperate seaweeds by means of mechanistic niche modelling on a European scale. The outcome of the mechanistic model is a species-specific response, the habitat suitability, which quantifies growth as a function of the temperature, salinity, light and nutrient requirements of the seaweed species. The mechanistic model was validated using independent distribution data and the validation statistic was good (area under the curve; AUC > 0.8) for five out of nine species and fair (0.6 < AUC < 0.8) for the remaining four species. The warm-temperate region extending from the coast of Portugal to the south coast of Brittany is currently a suitable habitat for most of the studied species. Fig. 4. Predicted habitat suitability for C. crispus under four climate change scenarios along a latitudinal gradient from 30°N to 80°N. The higher the expected emission, RCP2.6 (·········) < RCP4.5 (--) < RCP6.0 (−·−) < RCP8.5 (−··), the greater the northward shift. The distribution of C. crispus is predicted to shift under climate change, not to contract. The surface of optimal habitat (hs > 0.75) in de study area remains relatively constant under the future climate scenarios: 6.7 million km 2 under present conditions and 6.9, 6.8, 6.7 and 6.5 million km 2 for the future scenarios RCP 2.6, RCP 4.5, RCP 6.0 & RCP 8.5, respectively. The poleward shift is based on longitudinal averages and it is due to the high variation in habitat suitability on higher latitudes, caused by the topography of Europe, that the magnitude of the northward extension does not show in this figure (see Fig. 1C).

Table 4
Northward shift of the species' distribution (in kilometers) under four climate change scenarios (RCP), compared to present latitudinal distributions. The higher the expected emission (RCP 2.6 < RCP 4.5 < RCP 6.0 < RCP 8.5), the greater the northward shift. The mean temperate rise for the study area under the respective climate scenarios is also given, compared to the present sea surface temperatures. RCP  Accounting for future climate changes using the RCP scenarios allowed to predict the evolution of a species' distribution over time and to quantify the potential northward shift. Depending on the climate scenario and the species, the magnitude of this northward shift ranged from 110 km to 635 km. Our results can be used to select target species for seaweed aquaculture in a specific marine region. Similarly, the model may aid in the process of site selection for long-term, optimal growth conditions of the specific seaweed species included in this study. As such, our results contribute to the decision-making process in marine spatial planning of aquaculture activities.