INTRODUCTION
As the demand on terrestrial resources (energy, food and space) has intensified (e.g. Hubbert, Reference Hubbert1949; Hoogwijk et al., Reference Hoogwijk, Faaij, van den Broek, Berndes, Gielen and Turkenburg2003), there has been an increasing interest in utilizing marine resources. This is because oceans and seas cover over 70% of the Earth's surface and host an extraordinarily rich biodiversity. Nearshore areas in particular are among the most important biomes with all major taxonomic groups being represented and responsible for 90% of the world's marine primary production (Kaiser et al., Reference Kaiser, Attrill, Jennings, Thomas, Barnes, Brierley, Hiddink, Kaartokallio, Polunin and Raffaelli2011). The consequent increase and diversification in use of coastal natural resources has resulted in a necessity for stock-taking and possible protective measures for a variety of marine species (Halpern et al., Reference Halpern, Walbridge, Selkoe, Kappel, Micheli, D'Agrosa, Bruno, Casey, Ebert, Fox, Fujita, Heinemann, Lenihan, Madin, Perry, Selig, Spalding, Steneck and Watson2008). Obtaining information on the distribution of marine species has become an important goal; nevertheless, large-scale patterns of biodiversity and factors causing change in biodiversity are still not well understood, mostly due to the inconsistency of data sources and differences in methodology (e.g. Goulletquer et al., Reference Goulletquer, Gros, Boeuf and Weber2014).
Marine biodiversity varies over large scales whereas most studies focus on local and partly on regional scales. As direct mapping of large-scale patterns is extremely costly, knowledge-based management across large spatial scales requires alternative, low-cost sampling models as well as research strategies beyond the scope and vision of current national and institutional monitoring activities. There is also a requirement for a standard way of measuring local community patterns over large-scale environmental gradients. The recently built network of research locations in Europe (EMBOS) has addressed this issue by establishing a cost-effective monitoring programme for the diversity of hard bottom intertidal communities at a European scale. The principal aim of EMBOS is to assess broad-scale patterns in marine biodiversity in relation to natural and anthropogenic gradients. The assumption is that non-random patterns of variability in richness and cover of species in communities probably reflect the scale of underlining processes. Measuring variability at several spatial scales therefore enables the identification of the range of processes likely to be behind the observed differences among communities.
The observed patterns of each eco-region are not just defined by large-scale biogeographic processes, however (e.g. Roy et al., Reference Roy, Jablonski, Valentine and Rosenberg1998), but are also (and often interactively) a product of regional (Ramos et al., Reference Ramos, Puente and Juanes2016b) and local-scale variability in biogeochemical conditions (Dal Bello et al., Reference Dal Bello, Leclerc, Benedetti-Cecchi, Arvanitidis, van Avesaath, Bachelet, Bojanić, Como, Coppa, Coughlan, Crowe, Degraer, Espinosa, Faulwetter, Frost, Guinda, Ikauniece, Jankowska, Jourde, Kerckhof, Kotta, Lavesque, de Lucia, Magni, Fernandes de Matos, Orav-Kotta, Pavloudi, Pedrotti, Peleg, de la Pena, Puente, Ribeiro, Rigaut-Jalabert, Rilov, Rousou, Rubal, Ruginis, Pérez-Ruzafa, Silva, Simon, Sousa-Pinto, Troncoso, Warzocha, Weslawski and Hummel2016; Ramos et al., Reference Ramos, Puente, Guinda and Juanes2016a). The EMBOS programme has allowed a thorough examination of scale (see Dal Bello et al., Reference Dal Bello, Leclerc, Benedetti-Cecchi, Arvanitidis, van Avesaath, Bachelet, Bojanić, Como, Coppa, Coughlan, Crowe, Degraer, Espinosa, Faulwetter, Frost, Guinda, Ikauniece, Jankowska, Jourde, Kerckhof, Kotta, Lavesque, de Lucia, Magni, Fernandes de Matos, Orav-Kotta, Pavloudi, Pedrotti, Peleg, de la Pena, Puente, Ribeiro, Rigaut-Jalabert, Rilov, Rousou, Rubal, Ruginis, Pérez-Ruzafa, Silva, Simon, Sousa-Pinto, Troncoso, Warzocha, Weslawski and Hummel2016) but is also designed to measure the key environmental gradients at multiple spatial scales. The current study design matches with theoretical expectations that (1) local biological assemblages are defined by several interacting environmental gradients (Whittaker, Reference Whittaker1974); (2) local benthic diversity represents the interplay of historical evolutionary processes (i.e. the number of species in the regional species pool) combined with separate and interactive effects of local environmental gradients (e.g. Ricklefs & Schluter, Reference Ricklefs, Schluter, Ricklefs and Schluter1993; Holyoak et al., Reference Holyoak, Leibold, Moquet, Holt, Hoopes, Holyoak, Leibold, Moquet and Holt2005; Ramos et al., Reference Ramos, Puente, Guinda and Juanes2016a); and (3) these relationships are scale dependent i.e. local communities assemble from the macro- to the microscale and within these constraints, environmental stress plays a central role in shaping community diversity (Menge & Sutherland, Reference Menge and Sutherland1987).
In general, three main local gradients/stresses have been recognized in intertidal hard bottom habitats: wetness, exposure to wave action and salinity (Kaiser et al., Reference Kaiser, Attrill, Jennings, Thomas, Barnes, Brierley, Hiddink, Kaartokallio, Polunin and Raffaelli2011). When moving away from the waterline toward the highest part of the shore the environment becomes progressively drier. However, this gradient is amplified by waves and tides. Specifically, gently sloping areas tend to dissipate the wave energy arriving at the shore while steeper sloping shores experience much greater physical impact. In addition the wetness gradient is expected to be more pronounced on macrotidal shores (North Atlantic) compared with microtidal shores (Mediterranean and Baltic Seas). Nevertheless, features of rocky shore zonation patterns are similar worldwide independent of biogeographic region suggesting that similar abiotic forces (e.g. a combination of wave and tide effect) can generate similar biotic patterns (Stephenson & Stephenson, Reference Stephenson and Stephenson1972; Bird et al., Reference Bird, Franklin, Smith and Toonen2013 but see also Underwood, Reference Underwood1978). The salinity gradient occurs in estuaries and/or semi-enclosed seas and is often characterized by a significant loss of marine taxa and/or lowered diversity along a gradient of salinity reduction (e.g. Bonsdorff & Pearson, Reference Bonsdorff and Pearson1999).
In addition to natural environmental gradients, anthropogenic impacts are becoming increasingly important. To date, human activities have globally been linked to a diverse range of ecological changes in marine ecosystems including profound changes in coastal biodiversity (Goulletquer et al., Reference Goulletquer, Gros, Boeuf and Weber2014). Among the most commonly reported effects, climate change is expected to induce shifts in air and seawater temperature and affect benthic communities globally but often these effects are very context-specific and depend on the local environmental setting (Lima et al., Reference Lima, Ribeiro, Queiroz, Hawkins and Santos2007a, Reference Lima, Ribeiro, Queiroz, Xavier, Tarroso, Hawkins and Santosb, Kotta et al., Reference Kotta, Möller, Orav-Kotta and Pärnoja2014; Thorner et al., Reference Thorner, Kumar and Smith2014). In addition eutrophication is considered to be a significant concern in some regional seas (e.g. the Baltic, North Sea) (Conley et al., Reference Conley, Carstensen, Aigars, Axe, Bonsdorff, Eremina, Haahti, Humborg, Jonsson, Kotta, Lännegren, Larsson, Maximov, Rodriguez Medina, Lysiak-Pastuszak, Remeikaité-Nikiené, Walve, Wilhelms and Zillén2011) but locally also in the Atlantic Ocean and Mediterranean Sea (Howarth et al., Reference Howarth, Billen, Swaney, Townsend, Jaworski, Lajtha, Downing, Elmgren, Caraco, Jordan, Berendse, Freney, Kudeyarov, Murdoch and Zhu1996; Karydis & Kitsiou, Reference Karydis and Kitsiou2012). In order to assess the relative ecological importance of human impacts, it is necessary to understand how marine communities respond to various natural and anthropogenic ecological disturbances either separately or as interacting drivers of change. Such knowledge is of the utmost importance for understanding biodiversity patterns and assessing the consequences of various human activities on ecosystem integrity which in turn informs the management of natural resources.
Intertidal community structure is shaped by natural as well as by anthropogenic variables (e.g. Zacharias & Roff, Reference Zacharias and Roff2001; Crowe et al., Reference Crowe, Cusson, Bulleri, Davoult, Arenas, Aspden, Benedetti-Cecchi, Bevilacqua, Davidson, Defew, Fraschetti, Golléty, Griffin, Herkül, Kotta, Migné, Molis, Nicol, Noël, Sousa Pinto, Valdivia, Vaselli and Jenkins2013), however, only a handful of studies have been conducted that examine and compare the importance of these variables across distant aquatic ecosystems (Rees et al., Reference Rees, Pendle, Waldock, Limpenny and Boyd1999; Danovaro et al., Reference Danovaro, Canals, Gambi, Heussner, Lampadariou and Vanreusel2009; Narayanaswamy et al., Reference Narayanaswamy, Coll, Danovaro, Davidson, Ojaveer and Renaud2013; Schiele et al., Reference Schiele, Darr, Zettler, Berg, Blomqvist, Daunys, Jermakovs, Korpinen, Kotta, Nygård, von Weber, Voss and Warzocha2016). In the current study, we identified the role of multiple environmental variables on the patterns of cover and richness of intertidal hard-bottom communities throughout Europe and established the functional form relationships between these drivers and biotic patterns. The following hypotheses were put forward: (1) interregional differences account for a large part of the variability in hard bottom intertidal communities; (2) higher shore levels are characterized by reduced richness of macroalgae and invertebrates compared with lower levels; (3) elevated wave exposure and increased tidal range will result in systematically higher benthic richness compared with sheltered habitats; (4) elevated regional-scale loads of nutrients are expected to increase both cover and richness of understorey algae compared with cover and richness of canopy-forming algae; (5) temperature anomalies (higher than regular temperatures) are expected to reduce the growth of canopy forming algae and facilitate the growth of alternative understorey algae.
MATERIALS AND METHODS
Study area
Europe covers an extensive range of coastal areas spanning from subtropical to subpolar climatic regions. Such a wide climatic gradient leads to a great variety of environmental conditions impacting the diversity and stability of intertidal benthic macrophyte and invertebrate communities (e.g. Bulleri et al., Reference Bulleri, Benedetti-Cecchi, Cusson, Maggi, Arenas, Aspden, Bertocci, Crowe, Davoult, Eriksson, Fraschetti, Golléty, Griffin, Jenkins, Kotta, Kraufvelin, Molis, Sousa Pinto, Terlizzi, Valdivia and Paterson2012). A remarkable interregional variability exists: for example the semi-enclosed Mediterranean Sea with its oligotrophic conditions, high biodiversity and rates of endemism is very different to the Baltic Sea, which is characterized by eutrophic conditions and low biodiversity. Furthermore, the coasts of the North-East Atlantic Ocean, Arctic Ocean, Mediterranean, North and Baltic Seas differ greatly in terms of temperature, light availability, wave exposure and tidal regime.
Sampling and supporting environmental variables
The large-scale network of locations established through the implementation of EMBOS (http://embos.info/) was utilized for the study. In each study location the representative field sites were chosen based on commonalities and integrity of intertidal habitat as well as the lack of mobile sediment impacting bedrock. The network of locations covered all major European basins (Figure 1).
Map of the sampling stations in the study area. The annual average SST isotherms (°C) are also indicated.
Sampling was carried out early in biological spring for each region in 2014 (i.e. in late March/early April in the Mediterranean, April in the temperate NE Atlantic, end of April to mid-May in the Baltic, and about one month after the loss of ice-cover in the Arctic). At each location two stations were sampled. Each station was 5–15 m wide and separated by a distance of 50–100 m. The stations were placed in representative areas of the shore where commonly occurring species were observed i.e. avoiding rockpools and other atypical features.
At each station (N = 39), the assessment was carried out at mid intertidal and lower intertidal levels. The mid intertidal was defined operationally as approximately 25% of the vertical extent of the shore centred on Mean Tidal Level. The lower intertidal was defined as 25% of the vertical extent of the shore upwards from Mean Spring Low Water. In most of the sites, however, high intertidal and upper subtidal were also sampled. In subtidal systems, stations were placed in areas dominated by macro-algae, rather than in urchin barren habitats.
During sampling the biological characterization involved the assessment of per cent coverage of different taxa of canopy-forming and understorey algae as well as sessile invertebrates within each of five quadrats placed haphazardly within a site at each shore height. The understorey algae mostly consisted of fleshy and filamentous species. In macrotidal shores (e.g. NE Atlantic) the quadrat size was 0.5 × 0.5 m (0.25 m2) and in microtidal systems (e.g. Mediterranean) this was 0.2 × 0.2 m (0.04 m2). The size of quadrats differed between the two regions to account for differences in tidal amplitude and size of organisms. The studied micro-tidal habitats are often very narrow and larger quadrats would not have allowed the differentiation between the low-shore and mid-shore habitats in this micro-tidal system. Smaller quadrats in the NE Atlantic would not have been appropriate to sample usually large-sized fucoid and kelp algae. The gathered taxonomic information was optimized for up-scaling, so that a minimum amount of information was lost when datasets were coupled and integrated. This improved the efficiency and performance of large-scale approaches and analyses.
Geographic coordinates (as latitude and longitude) were recorded during sampling. Tidal range was obtained from tide tables. Regional-scale proxies for weather, eutrophication and climate change variables were obtained from different online databases. A mean climatology for wind, solar energy and cloud cover was obtained from the Climatic Research Unit at the University of East Anglia (https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1073). The climatology collection was compiled from existing data sources and algorithms, and was designed to satisfy the needs of modellers and investigators of the global carbon, water and energy cycle. The site data were interpolated as a function of latitude, longitude and elevation using thin-plate splines (Stackhouse & Gupta, Reference Stackhouse, Gupta, Hall, Collatz, Meeson, Los, De Colstoun and Landis2012).
Sea ice concentration (expressed as a mean percentage of ocean area covered by sea ice) was obtained from the Goddard Space Flight Center, Sea Ice Concentrations from Nimbus-7 Scanning Multichannel Microwave Radiometer and the Defense Meteorological Satellites Program Special Sensor Microwave/Imager passive microwave data (https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=981). These original data were re-gridded by the National Snow and Ice Data Center from their original 25 km spatial resolution and EASE-Grid into equal angle Earth grids with quarter degree spatial resolutions in latitude/longitude (Armstrong & Knowles, Reference Armstrong, Knowles, Hall, Collatz, Meeson, Los, De Colstoun and Landis2010).
Nutrient pollution, sea surface temperatures and temperature anomalies data were obtained from Halpern et al. (Reference Halpern, Walbridge, Selkoe, Kappel, Micheli, D'Agrosa, Bruno, Casey, Ebert, Fox, Fujita, Heinemann, Lenihan, Madin, Perry, Selig, Spalding, Steneck and Watson2008). In their study nutrient pollution data came from Food and Agriculture Organization (FAO) national statistics (http://faostat.fao.org) on average annual use of fertilizers (nutrients) for the years 1993–2002 and were distributed across landscapes in agricultural lands with dasymetric techniques. Sea surface temperature (SST) data from 1985 to 2005 have been constructed from Advanced Very High Resolution Radiometer Pathfinder Version 5.0 by NOAA's National Oceanographic Data Center and the University of Miami's Rosenstiel School of Marine and Atmospheric Science (http://pathfinder.nodc.noaa.gov). Then it was assessed how many times the SST anomaly exceeded the standard deviation of SSTs for each 4 × 4 km grid cell and week of the year. By incorporating the standard deviation, this threshold-based approach accounts for natural variability at a given location, which can vary widely from place to place. Finally, a SST change metric (later referred to as temperature anomalies) was developed by subtracting the number of non-zero positive anomalies in the early period (1985–1990) from the number in the recent period (2000–2005). As the anthropogenic driver data (nutrient load and temperature anomalies) had extreme left-skewed distributions the data were log[X + 1]-transformed and then rescaled between 0–1, with the highest log-transformed value set = 1.
Boosted regression tree modelling
Environmental gradients are expected to affect biota interactively and when exploring potential influences on the diversity of intertidal hard bottom communities, the analytical framework should have capabilities to quantify the functional form of the relationships between separate environmental gradients and the biota, as well as their interactive effects. Modelling is the most appropriate tool, and several refined statistical approaches have already been applied in the field (Herkül et al., Reference Herkül, Kotta, Kutser and Vahtmäe2013; Sandman et al., Reference Sandman, Wikström, Blomqvist, Kautsky and Isaeus2013). The main focus of spatial modelling has often been directed towards predicting species distribution ranges (Peterson, Reference Peterson2003; Soberon & Peterson, Reference Soberon and Peterson2005) without investigating underlying mechanisms. Modelling can however improve our understanding of the relationship between the environment and biota and inform the theoretical understanding of biotic patterns and their driving forces (Kotta et al., Reference Kotta, Möller, Orav-Kotta and Pärnoja2014).
A novel predictive modelling technique, Boosted Regression Trees (BRT) does not use any predefined data model but instead uses an algorithm to analyse the relationship between the biota and environmental variables. BRT copes with different non-linear relationships which are common in ecological data but difficult to analyse using more traditional methods. The BRT also avoids overfitting the data, thereby providing robust estimates. Most importantly from an ecological perspective, BRT automatically detects and models interactive effects between predictors. Due to its strong predictive performance, BRT is increasingly used in ecology (Elith et al., Reference Elith, Leathwick and Hastie2008; Kotta et al., Reference Kotta, Kutser, Teeveer, Vahtmäe and Pärnoja2013).
Multicollinearity can be an issue with modelling when explaining relationships between environment and the biota. Thus, prior to modelling, a correlation analysis was conducted for environmental variables. Among the studied environmental variables only water temperature, solar radiation and cloudiness showed any correlation (r between 0.82 and 0.90, P < 0.01). Thus, in order to avoid multicollinearity issues, the final models included only water temperature.
The contribution of different environmental variables to the cover and richness of hard bottom intertidal macroalgae and invertebrates was explored using the BRT modelling. The BRT technique iteratively develops a large ensemble of small regression trees constructed from random subsets of the data. Each successive tree predicts the residuals from the previous tree to gradually boost the predictive performance of the overall model (Elith et al., Reference Elith, Leathwick and Hastie2008). In fitting a BRT the learning rate and the tree complexity must be specified. The learning rate determines the contribution of each successive tree to the final model, as it proceeds through the iterations. The selected value of a tree complexity determines whether only main effects (tree complexity = 1) or interactions are also included (tree complexity >1). Ultimately, the learning rate and tree complexity combined determine the total number of trees in the final model. Following the suggestions by Elith et al. (Reference Elith, Leathwick and Hastie2008) the model learning rate was kept at 0.1 and tree complexity at 5 for all models. It was also checked that the final models had more than 1000 trees. Nevertheless, a selection of model parameters had only marginal impact on model performance with optimal models improving predictions less than 1%. Model performance was evaluated using the cross validation statistics calculated during model fitting (Hastie et al., Reference Hastie, Tibshirani and Friedman2009). Thus, when running models, a random selection of 80% of the data was used for training the model and the rest of the data i.e. 20% was assigned for testing model accuracy. The BRT modelling was done in the statistical software R using the gbm package (RDC Team, 2013).
RESULTS
Outputs of the BRT models described a significant proportion of variability in the cover and richness of macroalgae and invertebrates inhabiting hard bottom intertidal habitats. For the majority of models the explained variance was higher than 80%. The model explained 71% of overall variability when only considering cover of canopy algae.
The results of the BRT models including region (i.e. basin) as an explanatory variable showed that interregional variability was very important for hard bottom intertidal macrophytes and invertebrates. There was a large difference however in the role of regional variability among types of organisms (understorey algae, canopy-forming algae or sessile invertebrates) and aggregated community variables (cover or richness). In general, regional variability was a better explanation for algal patterns rather than invertebrate patterns and total cover rather than richness. Regional variability explained 60% of the variability in canopy algae cover and 53% of the understorey algae cover but only 37% of the sessile invertebrate cover. In the richness models, the same values were about one-third lower than the levels in the cover models. Understorey algae had the highest cover in the English Channel, the Irish Sea and the Strait of Gibraltar and highest richness in the Bay of Biscay. Canopy algae had the highest cover in the Irish Sea, the Mediterranean Sea, the English Channel and the Baltic Sea and highest richness in the Bay of Biscay, the English Channel, the Mediterranean Sea and the Mid North Atlantic coasts. Both the cover and richness of sessile invertebrates were highly variable with no clear spatial patterns (Figure 2).
Average values of the cover and richness of understorey algae, canopy algae and sessile invertebrates in different regions. The code of tidal levels is as follows: H – high intertidal, M – mid intertidal, L – low intertidal, S – subtidal.
As our main goal was to quantify roles of different environmental gradients on intertidal algae and invertebrates and regional variability largely incorporates changes in various natural and anthropogenic forcing, the BRT models below did not include the variable ‘region’ among the predictors. Such BRT models identified tidal level, tidal range, nutrient loading and water temperature (anomalies) to be the most significant predictors of both cover and richness of macroalgae and invertebrates (Table 1).
Relative contribution of environmental predictors in the BRT models of total cover and richness of macroalgae and invertebrates.
The cover and richness of understorey algae was higher in subtidal and low intertidal areas compared with medium and high intertidal areas. Moreover, higher tidal range was associated with elevated values of cover and richness of understorey algae. Nutrient load was also among the most important variables predicting the understorey algal patterns. The cover of understorey algae increased with nutrient load whereas the richness of understorey algae was lowest at intermediate loads (Figure 3).
Standardized functional-form relationships showing effects of the four most important environmental variables on the cover and richness of understorey algae in the study area, whilst all other variables are held at their means. The variables are ordered by their relative contribution in the BRT model (shown in brackets). Upward tickmarks on x-axis show the frequency of distribution of data along this axis. The code of tidal levels is as follows: H – high intertidal, M – mid intertidal, L – low intertidal, S – subtidal. The values of nutrient loads and temperature anomalies are log-transformed and rescaled between 0–1. See the Methods section for further information on environmental variables.
As for the understorey algae, cover and richness of canopy algae was higher in subtidal and low intertidal areas compared with medium and high intertidal areas. Elevated tidal range reduced the cover of canopy algae but increased the associated taxonomic richness. Water temperature strongly contributed to the patterns of canopy algae with cover and richness being higher at reduced water temperatures. In contrast to understorey algae, the effect of nutrient load for patterns of canopy algae was moderate. Specifically, elevated nutrient loads were associated with higher cover of canopy algae whereas highest richness values were found at intermediate nutrient loads (Figure 4).
Standardized functional-form relationships showing effects of the four most important environmental variables on the cover and richness of canopy algae in the study area, whilst all other variables are held at their means. The variables are ordered by their relative contribution in the BRT model (shown in brackets). Upward tickmarks on x-axis show the frequency of distribution of data along this axis. The code of tidal levels is as follows: H – high intertidal, M – mid intertidal, L – low intertidal, S – subtidal. The values of nutrient loads and temperature anomalies are log-transformed and rescaled between 0–1. See the Methods section for further information on environmental variables.
The richness of canopy algae and sessile invertebrates showed similar patterns along the tidal levels and temperature gradient. However, no such similarities emerged for the cover of sessile invertebrates. Instead, the cover of sessile invertebrates was higher in high and medium intertidal areas compared with low intertidal and subtidal areas. Nutrient loads and tidal range highly contributed to the cover of sessile invertebrates with higher cover values expected at intermediate nutrient loads and tidal range (Figure 5).
Standardized functional-form relationships showing effects of the four most important environmental variables on the cover and richness of sessile invertebrates in the study area, whilst all other variables are held at their means. The variables are ordered by their relative contribution in the BRT model (shown in brackets). Upward tickmarks on x-axis show the frequency of distribution of data along this axis. The code of tidal levels is as follows: H – high intertidal, M – mid intertidal, L – low intertidal, S – subtidal. The values of nutrient loads and temperature anomalies are log-transformed and rescaled between 0–1. See the section of methods for further information on environmental variables.
The BRT analyses also identified a few important interactive effects. Specifically, the effect of nutrient loads on understorey algae and sessile invertebrates varied widely among tidal levels. Although the cover of understorey algae increased with nutrient load across tidal levels, the effect was stronger in subtidal and low intertidal areas compared with medium and high intertidal areas. Similarly, the cover and richness of sessile invertebrates increased with nutrient loads in contrast to understorey algae where the stronger effects were found in high and medium intertidal areas compared with low intertidal and subtidal areas. The richness of understorey algae decreased with nutrient loads with the effect being proportionally stronger in subtidal and low intertidal areas compared with medium and high intertidal areas (Figure 6).
Partial dependence plots of the BRT model for the cover and richness of understorey algae and sessile invertebrates at different tidal levels in the study area. The values of nutrient loads are log-transformed and rescaled between 0–1. See the Methods section for further information on environmental variables.
Water temperature interacted with other environmental variables in the models of canopy algae. The cover of canopy algae decreased at elevated temperatures; however, this effect was much stronger in subtidal and low intertidal areas compared with medium and high intertidal areas. The richness of canopy algae was reduced by elevated ice scour. However, the effect disappeared in areas characterized by higher temperature anomalies (Figure 7).
A partial dependence plot of the BRT model for the cover of canopy algae at different tidal levels and a three-dimensional partial dependence plot in the BRT model for the richness of canopy algae along the gradients of water temperature anomalies and ice cover in the study area. The values of temperature anomalies are log-transformed and rescaled between 0–1. See the Methods section for further information on environmental variables.
DISCUSSION
Based on the results of the analyses, initially proposed hypotheses were supported. The results confirmed that interregional variability accounted for about half of the variability in intertidal hard bottom macroalgal and invertebrate communities at the pan-European scale. Interregional variability displayed a systematically stronger effect on algal communities than on invertebrate communities and had a significantly stronger effect on cover patterns compared with richness of species. It was found that higher shore levels were characterized by reduced algal and invertebrate richness although elevated wave exposure and increased tidal range mitigated the effect. Elevated nutrient loads increased both cover and richness of understorey species but had a moderate effect on canopy algae and increases in seawater temperature and associated weather anomalies did not seem to have any effect on the cover of canopy algae whereas a very strong negative effect on canopy richness was observed. On the other hand, the patterns of understorey algae were essentially uncoupled with those of elevated temperature anomalies.
There exist broad differences in how biogeographic history affects local patterns of species cover and richness and these differences may explain why regional variability was more important for species cover than richness. Specifically, the number of species in the regional pool sets the maximum limit of species in a local community. At small spatial scales, however, local communities become quickly saturated by species and the number of species in the species pool is much higher than local communities can actually accommodate (Kotta & Witman, Reference Kotta, Witman and Wahl2009). Earlier studies have found clear evidence for saturation at the studied quadrat scale with the saturation effect being amplified towards the coastline characterized by high environmental stress (Russell et al., Reference Russell, Wood, Allison and Menge2006). Therefore, local species richness of the hard bottom intertidal habitats is defined either by limiting environmental factors or short-term biotic mechanisms that leads to the exclusion of some of the species from a community rather than regional speciation processes. As observed in this study, intertidal hard bottom habitats are often highly patchy, especially at the smallest spatial scales with different microhabitats being suitable for different sets of species. Such microhabitat heterogeneity also suggests the existence of a strong small-scale component of variability and a weak regional pool effects in intertidal species richness.
Total cover values, in contrast to species richness, are determined by environmental stressors such as desiccation, wave exposure, ice scour and nutrient loads that often play a central role in defining community structure (Menge & Sutherland, Reference Menge and Sutherland1987; Smale et al., Reference Smale, Burrows, Evans, King, Sayer, Yunnie and Moore2016). Such environmental stressors do not only operate at a regional scale, but can even introduce significant variability within the same region (Smale et al., Reference Smale, Burrows, Evans, King, Sayer, Yunnie and Moore2016) and thus synchronize cover values at multiple spatial scales, particularly where they are having a direct effect on recruitment or mortality (Dal Bello et al., Reference Dal Bello, Leclerc, Benedetti-Cecchi, Arvanitidis, van Avesaath, Bachelet, Bojanić, Como, Coppa, Coughlan, Crowe, Degraer, Espinosa, Faulwetter, Frost, Guinda, Ikauniece, Jankowska, Jourde, Kerckhof, Kotta, Lavesque, de Lucia, Magni, Fernandes de Matos, Orav-Kotta, Pavloudi, Pedrotti, Peleg, de la Pena, Puente, Ribeiro, Rigaut-Jalabert, Rilov, Rousou, Rubal, Ruginis, Pérez-Ruzafa, Silva, Simon, Sousa-Pinto, Troncoso, Warzocha, Weslawski and Hummel2016). For example, stressors such as wave exposure, ice scour and desiccation at upper intertidal zones have a notable impact on the abundance patterns of many sessile intertidal species (Dayton, Reference Dayton1971; Heaven & Scrosati, Reference Heaven and Scrosati2008). When the intensity or the severity of disturbances outweighs colonization, then the inhabited area remains smaller in size than potentially possible (Johnston et al., Reference Johnston, Keough and Qian2002). Regional variability in productivity and consumption may also have a fundamental influence on patterns of species cover. It has been shown that high productivity regions are more prone to overgrazing (functional types of grazers which cause barrens of diversity) than to nutrient loads, and low productivity regions are less affected by grazing (no barrens) and more affected by nutrients (Connell & Irving, Reference Connell and Irving2008).
We also expected that environmental gradients are major underlying components of the observed interregional variability. Apart from the potential effects of the regional species pool and species biogeography not resolved in this study, our BRT models showed that most variability in cover and richness of hard bottom macroalgae and invertebrates across Europe was explained by site-specific abiotic characteristics such as tidal regime and climate variables coupled with anthropogenic eutrophication.
More specifically, we expected higher shores to be characterized by reduced richness of macroalgae and invertebrates compared with lower shores. The BRT analyses confirmed this expectation as both cover and richness of understorey and canopy algae increased at lower intertidal and subtidal zones compared with mid- and high intertidal zones. Such findings confirm well-known principles of intertidal ecology whereby higher cover and richness of algal communities is found at lower intertidal and subtidal zones due to reduced risk of desiccation and thermal stress (Dayton, Reference Dayton1975) but in some cases probably due to reduced grazing or a combination of both (Underwood & Jernakoff, Reference Underwood and Jernakoff1981; Hawkins & Hartnoll, Reference Hawkins and Hartnoll1985). Species richness of sessile invertebrates also increased in line with the previously observed zonation patterns for algal communities which can be linked to physiological tolerances of species, foundation effects due to canopy-forming seaweeds and broader availability of suitable habitats as space, these being some of the most limiting factors for intertidal communities (Tomanek & Helmuth, Reference Tomanek and Helmuth2002). On the other hand, the cover of sessile invertebrates displayed the complete opposite i.e. cover increased at high and medium intertidal areas compared with low intertidal and subtidal areas. This most probably relates to the fact that competition and predation pressure decreases at higher intertidal zones allowing an increase in the cover of species that are able to effectively colonize high intertidal areas (Dayton, Reference Dayton1971).
Elevated wave exposure and increased tidal range were also expected to support higher benthic richness compared with sheltered and/or microtidal habitats. The BRT analyses showed that tidal range accounted for a large variability of the cover and richness patterns of benthic intertidal species. High-energy shores were characterized by high cover and richness of understorey algae; thus, an elevated tidal range probably reduced the desiccation stress of small non-corticated filamentous macroalgae, otherwise common in microtidal shores. The tide height also relates to the tidal stream strength and wave energy reaching the shore so lower cover of larger organisms (canopy algae and some invertebrates) is associated with higher tidal range possibly due to mechanical disturbance related to elevated disturbance by waves (Denny et al., Reference Denny, Miller, Stokes, Hunt and Helmuth2003; Helmuth & Denny, Reference Helmuth and Denny2003).
Elevated regional-scale loads of nutrients were expected to increase both cover and richness of understorey algae species, whereas relationship between nutrient loads and cover and richness of canopy-forming species is bell-shaped. Observed nutrient loads displayed a significant separate effect on the cover and richness of the intertidal communities. It was found that elevated nutrient loads increased both cover and richness of understorey species but had a moderate role on canopy algae. Such observations might relate to functional traits of the studied organisms. For example, understorey algae and sessile invertebrates are both characterized by a short lifespan, thus, elevated nutrient availability might be beneficial throughout their life cycle. Beneficial impacts of nutrients on filamentous understorey algae correlates well with previous findings suggesting shifts in floral dominance patterns from perennial to ephemeral species in case of prolonged events of nutrient influx (Gorostiaga & Diez, Reference Gorostiaga and Diez1996; Diez et al., Reference Diez, Secilla, Santolaria and Gorostiaga1999). However, the highest richness values of understorey algae were observed when nutrient loads were the lowest, minimum richness at intermediate loads and moderate richness at highest loads. This could potentially reflect shifts in taxonomic composition along the eutrophication gradient as communities under pristine conditions differ significantly from those under eutrophic conditions.
The increase in nutrient loads was positively correlated with the cover of canopy algae, but the slope of the function was not as steep when compared with understorey algae. Canopy algae have internal nutrient reserves and might therefore be less influenced by external nutrient availability (Pedersen & Borum, Reference Pedersen and Borum1996). Due to relative stability in internal nutrient reserves the canopy algae are favoured over smaller ephemeral species at low nutrient environments. On the other hand, the highest canopy richness was observed at intermediate nutrient levels that resulted in a bell-shaped relationship. In oligotrophic conditions an increment of nutrients promotes richness by directly improving algal growth conditions and indirectly increasing habitat heterogeneity (Field et al., Reference Field, Behrenfeld, Randerson and Falkowski1998; Kotta et al., Reference Kotta, Möller, Orav-Kotta and Pärnoja2014) whereas nutrient overload is detrimental for many canopy-forming algae and a steep decrease in cover and richness values can occur. The latter phenomenon might be a result of algal blooms or increased epiphyte overgrowth that reduces light availability for macroalgae and/or induces supplemental herbivory (Valiela et al., Reference Valiela, Mcclelland, Hauxwell, Behr, Hersh and Foreman1988; Orav-Kotta & Kotta, Reference Orav-Kotta and Kotta2004; Kraufvelin et al., Reference Kraufvelin, Ruuskanen, Nappu and Kiirikki2007). Ultimately, this leads to competitive exclusion where only the best competitors for light are favoured (Abrams, Reference Abrams1995).
Cover of sessile invertebrates was highest at intermediate nutrient load. Here, several mechanisms can generate the observed functional form relationship but the exact mechanism can only be validated experimentally. Most plausibly, filamentous understorey macroalgae outcompete invertebrates for space in areas of high nutrient loads.
Finally, it was expected that elevated temperatures and increased temperature anomalies would facilitate understorey algal species but hinder the growth of canopy-forming algae. The BRT models demonstrated that temperature described a large proportion of the variability of intertidal communities, especially those of canopy algae. It was found that climate-induced increases in seawater temperature and associated weather anomalies did not seem to have any effect on the cover of canopy algae whereas a very strong negative effect on canopy richness was observed. This suggests that the maintenance of macroalgae-related functions of understorey assemblages in coastal environments (Britton-Simmons, Reference Britton-Simmons2006; Wikström & Kautsky, Reference Wikström and Kautsky2007; Golléty et al., Reference Golléty, Migné and Davoult2008) may be threatened by climate change as species richness is often considered as a reliable measure to reflect the state of an ecosystem (Maestre et al., Reference Maestre, Quero, Gotelli, Escudero, Ochoa, Delgado-Baquerizo, García-Gómez, Bowker, Soliveres, Escolar, García-Palacios, Berdugo, Valencia, Gozalo, Gallardo, Aguilera, Arredondo, Blones, Boeken, Bran, Conceição, Cabrera, Chaieb, Derak, Eldridge, Espinosa, Florentino, Gaitán, Gatica, Ghiloufi, Gómez-González, Gutiérrez, Hernández, Huang, Huber-Sannwald, Jankju, Miriti, Monerris, Mau, Morici, Naseri, Ospina, Polo, Prina, Pucheta, Ramírez-Collantes, Romão, Tighe, Torres-Díaz, Val, Veiga, Wang and Zaady2012) and is regularly found as a basis of management advice. It is important to note that species richness incorporates taxonomic identity of species in an ecosystem that is also linked to functional traits of organisms. If the richness of canopy algae should decrease then it is highly questionable whether sustained rates of high canopy cover could provide similar functions required by the ecosystem as well as whether full recovery is likely. For example, in temperate Western Australia a canopy dominated temperate reef system switched to a turf dominated state due to an exceptional heatwave in 2011 and the recovery to the previous state is further hindered by high grazing rates of herbivores (Smale & Wernberg, Reference Smale and Wernberg2013; Wernberg et al., Reference Wernberg, Smale, Tuya, Thomsen, Langlois, de Bettignies, Bennett and Rousseaux2013). Increased frequency of disturbance events and rising temperatures are expected to cause species range shifts (Lima et al., Reference Lima, Ribeiro, Queiroz, Hawkins and Santos2007a) and pose a significant threat to the distribution of canopy-forming algae and jeopardize the stability of algal-dominated marine regions as well as the overall food-web stability (Wernberg et al., Reference Wernberg, Thomsen, Tuya, Kendrick, Staehr and Toohey2010; Smale et al., Reference Smale, Burrows, Moore, O'Connor and Hawkins2013). Unexpectedly, the patterns of understorey algae were essentially uncoupled with those of elevated temperature anomalies. One likely explanation for such a phenomenon could be that filamentous algal species dominate in understorey communities and rapidly replace one another when environmental conditions become unfavourable to some taxa. The major loss of canopy algae highlights the need to understand whether filamentous algae could compensate for such a loss and provide similar functions to an ecosystem as canopy algae. For example, it has been demonstrated that filamentous algae can be more productive than canopy-forming algae on a biomass basis whereas the annual net production per area is two to seven times lower than for canopy-forming algae (Copertino et al., Reference Copertino, Connell and Cheshire2005). Likewise, it has been shown that filamentous algae may host similar invertebrate species as perennial macroalgal species (Wikström & Kautsky, Reference Wikström and Kautsky2007). However, due to their seasonal occurrence, ephemeral species are not likely to replace the majority of functions provided by perennial species such as habitat formation and associated diversity (Råberg & Kautsky, Reference Råberg and Kautsky2007; Arnold et al., Reference Arnold, Teagle, Brown and Smale2016).
All components of the environment and biota are interlinked and depend on each other. Thus, any model that resolves only separate effects of the environment on biota is only a simplified representation of reality. In this study, we assessed the importance of different environmental interactions and found that, surprisingly, only a few interactive effects emerged as being important. For the ephemeral species, the effect of nutrient loads varied along a gradient of tidal level. This may be partly due to the fact that the sensitivity of habitats to eutrophication decreases with the intensity of water exchange (e.g. Tett et al., Reference Tett, Gilpin, Svendsen, Erlandsson, Larsson, Kratzer, Fouilland, Janzen, Lee, Grenz, Newton, Ferreira, Fernandes and Scory2003) and areas higher in the shore have lower exposure to water enriched with nutrients. Alternatively, elevated intensities of other stressors in the higher shore levels (including elevated rates of desiccation or ice scour, the effects of nutrient loads) are expected to be weaker (Jutterström et al., Reference Jutterström, Andersson, Omstedt and Malmaeus2014). For the canopy-forming species, tidal level moderated the response of temperature to species cover. The high inter-tidal level–temperature interaction had the weakest effect on the cover of canopy algae as these algae can be highly adapted to elevated temperatures/desiccation events (Guenther & Martone, Reference Guenther and Martone2014). Similarly, our study showed that elevated temperature anomalies may counterbalance the effects of ice scour on the richness of canopy algae, promoting an elevated richness of canopy algae at high latitudes (especially in mid and high-intertidal areas). This correlates with recent documentation of the current effects of seawater warming and ice retreat on littoral macrophytes at high latitudes (Weslawski et al., Reference Weslawski, Wiktor and Kotwicki2010). On the other hand, elevated temperature anomalies reduced algal cover at low latitudes (especially in low intertidal areas).
Correlative species distribution models are frequently applied as a sophisticated tool to improve our understanding of the relationship between environment and biota (Elith et al., Reference Elith, Graham, Anderson, Dudík, Ferrier, Guisan, Hijmans, Huettmann, Leathwick, Lehmann, Li, Lohmann, Loiselle, Manion, Moritz, Nakamura, Nakazawa, Mc Overton, Peterson Townsend, Phillips, Richardson, Scachetti-Pereira, Schapire, Soberón, Williams, Wisz and Zimmermann2006, Reference Elith, Kearney and Phillips2010). It is important that the analytical framework has the capacity to quantify the functional form of separate and interactive relationships in order to apply potential reasoning for causations. The machine learning methodology used in this study provides a robust framework to seek for the mechanisms that determine spatial distribution patterns of intertidal species under current environmental conditions. The established functions between environmental gradients and anthropogenic stressors can potentially be used to predict the cover and diversity of intertidal hard bottom communities at the pan-European scale and seek for potential biodiversity hotspots. Such knowledge provides better starting conditions to understand various aquatic ecosystems and contributes towards environmental policies such as marine ecosystem-based management that seeks to manage marine resources in accordance with ecosystem health while providing socio-economic services needed by people (Crowder & Norse, Reference Crowder and Norse2008; Espinosa-Romero et al., Reference Espinosa-Romero, Chan, McDaniels and Dalmer2011; Long et al., Reference Long, Charles and Stephenson2015). However, care should be taken when interpreting cause-effect relationships based on the results of this study as further experimental research is required to validate the obtained results. This is mainly because we were only able to include one temporal replicate of sampled sites which poses a notable limitation for interpretation of results, but is understandable given the scale of this study. Nevertheless, extreme weather conditions such as heavy winter storms in the NE Atlantic Ocean or explicitly warm summers in the Mediterranean Sea could result in atypically sparse algal canopies. Regional deviations from typical disturbance history may therefore affect the credibility of data obtained in one sampling season as such large-scale environmental stochasticity may override the effect of environmental variability operating at smaller spatial scales.
In this study we demonstrate how the distribution patterns of benthic hard bottom algal and invertebrate communities at pan-European scale are affected by various environmental gradients and human-induced stressors. For the majority of large-scale study initiatives the main difficulty is the amount of time required in order to overcome restrictions to access required data and/or to harmonize different databases. This amplifies the need for global science to become more unified and opened through data sharing as data compiled in this study is virtually based on the voluntary initiative of scientists and good will. Projects such as EMBOS, that go beyond the scope and vision of current national and institutional monitoring activities, promote a simple and standardized sampling framework required in order to address complications associated with large-scale studies that examine biological mechanisms of life on earth.
ACKNOWLEDGEMENT
We are thankful to Enrique Ostalé-Valriberas for his kind support during the sampling programme and Caroline Broudin, Fanny Noisette, Aline Migné, Céline Houbin, Laurent Lévêque, Dominique Davoult, François Vandenbosch, Pascal Riera, Elena Maggi, Ofrat Raveh and Elvira Ramos for their help with sampling.
FINANCIAL SUPPORT
This article is based upon work from COST Action ES1003 Development and implementation of a pan-European Marine Biodiversity Observatory System (EMBOS), supported by COST (European Cooperation in Science and Technology). Jonne Kotta, Holger Jänes and Helen Orav-Kotta were supported by Institutional research funding IUT02-20 of the Estonian Research Council and the BONUS project BAMBI, the joint Baltic Sea research and development programme (Art 185), funded jointly from the European Union's Seventh Programme for research, technological development and demonstration and from the Estonian Research Council. Pedro Ribeiro was supported by FCT through post-doctoral grant ref. SFRH/BPD/69232/2010 funded through QREN and COMPETE, and the strategic project UID/MAR/04292/2013 granted to MARE. Jérōme Jourde was financially supported by the Région Poitou-Charentes through CPER funding, La Rochelle University and CNRS. Martina Dal Bello and Lisandro Benedetti-Cecchi were supported by a grant from the Italian Ministry for Research and Education (PRIN project ‘Biocostruzioni costiere: struttura, funzione, e gestione’). Gil Rilov was partly supported by the Israeli Science Foundation (grant number 1217/10) and the Israel Ministry of Environmental Protection. Valentina Kirienko Fernandes de Matos was supported by the Portuguese Science Foundation (FCT) through a doctoral grant (ref. SFRH/BD/86390/2012). Emilia Jankowska and Jan Marcin Węsławski were supported by the Statutory Funds of Institute of Oceanology Polish Academy of Sciences.
CONFLICT OF INTEREST
None.
ETHICAL STANDARDS
The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional guides.
ETHICS POLICY
Authors complied with best practice and a high standard of ethics. Authors declare that the work has not been previously published or submitted for publication elsewhere, appropriate ethics conduct and approval were adhered to where field campaigns were conducted such as there were no destruction of habitat or killing of species, none of the text has been copied or paraphrased from other publications.
CONTRIBUTORS
The authors of this article are partners within the COST Action EMBOS on a European Marine Biodiversity Observatory System. First author of the paper has performed statistical analyses and led the writing. Co-authors have made a significant contribution to the conception, design, execution, or interpretation of the reported study. All authors have commented on drafts of the manuscript and have approved the final article to be true and included in the disclosure.
