INTRODUCTION
The marine environment is of a transboundary nature and needs to be studied at relevant scales of space and time. For instance, along Europe a north to north-east shift in the distribution of several marine species has been observed (Hummel et al., Reference Hummel, Bogaards, Bachelet, Caron, Sola and Amiard-Triquet2000; Beaugrand et al., Reference Beaugrand, Reid, Ibanez, Lindley and Edwards2002; Mieszkowska et al., Reference Mieszkowska, Kendall, Hawkins, Leaper, Williamson, Hardman-Mountford and Southward2006; Jansen et al., Reference Jansen, Pronker, Kube, Sokolowski, Sola, Marquiegui, Schiedek, Wolowicz, Wendelaar Bonga and Hummel2007). Global climate change is said to be one of the causes. However, the degree and impact of such biodiversity changes and its causes and consequences remain largely unknown (Heip et al., Reference Heip, Hummel, van Avesaath, Appeltans, Arvanitidis, Aspden, Austen, Boero, Bouma, Boxshall, Buchholz, Crowe, Delaney, Deprez, Emblow, Feral, Gasol, Gooday, Harder, Ianora, Kraberg, Mackenzie, Ojaveer, Paterson, Rumohr, Schiedek, Sokolowski, Somerfield, Sousa Pinto, Vincx, Weslawski and Nash2009). The question thereby still is whether changes in the level of diversity in an ecosystem really matter, since despite a huge difference in level and type of diversity between different seas, e.g. in the Mediterranean (more than 17,000 plant and animal species), the North Sea (more than 1500 species) and the central Baltic (only 73 marine species), trophic relationships in these systems are assumed to be similar (Elmgren & Hill, Reference Elmgren, Hill, Ormond, Gage and Angel1995; Coll et al., Reference Coll, Piroddi, Steenbeek, Kaschner, Ben Rais Lasram, Aguzzi, Ballesteros, Bianchi, Corbera, Dailianis, Danovaro, Estrada, Froglia, Galil, Gasol, Gertwagen, Gil, Guilhaumon, Kesner-Reyes, Kitsos, Koukouras, Lampadariou, Laxamana, López-Fé de la Cuadra, Lotze, Martin, Mouillot, Oro, Raicevich, Rius-Barile, Saiz-Salinas, San Vicente, Somot, Templado, Turon, Vafidis, Villanueva and Voultsiadou2010; Magni et al., Reference Magni, Rajagopal, Como, Jansen, van der Velde and Hummel2013; Zettler et al., Reference Zettler, Karlsson, Kontula, Gruszka, Laine, Herkül, Schiele, Maximov and Haldin2014). Determining the patterns of biodiversity of benthic organisms and the factors that may explain them requires an integrated research strategy, which is beyond the tradition, capabilities and scales of classic research (Heip & Hummel, Reference Heip and Hummel2000). The first ideas for such integrated research, involving a large-scale pan-European network of marine observatories, have been formulated during the FP5 (EC 5th Framework Programme) project BIOMARE (Implementation and networking of large scale, long-term marine biodiversity research in Europe), which was initiated by the European Network of Marine Research Stations (MARS). The subsequent FP6 Network of Excellence MarBEF (Marine Biodiversity and Ecosystem Functioning) adopted, and focused on, the integration of datasets and joint research (Escaravage et al., Reference Escaravage, Herman, Merckx, Wlodarska-Kowalczuk, Amouroux, Degraer, Grémare, Heip, Hummel, Karakassis, Labrune and Willems2009; Heip et al., Reference Heip, Hummel, van Avesaath, Appeltans, Arvanitidis, Aspden, Austen, Boero, Bouma, Boxshall, Buchholz, Crowe, Delaney, Deprez, Emblow, Feral, Gasol, Gooday, Harder, Ianora, Kraberg, Mackenzie, Ojaveer, Paterson, Rumohr, Schiedek, Sokolowski, Somerfield, Sousa Pinto, Vincx, Weslawski and Nash2009).
These projects have led to recommendations for the selection of sites and indicators to monitor marine ecosystems and their biodiversity at a pan-European scale. The aforementioned activities have been taken under the COST Action ES1003 EMBOS (Development and implementation of a pan-European Marine Biodiversity Observatory System), which lasted from 2011 to 2015. EMBOS focused on the following goals: (1) to build a large-scale pan-European network of marine observatories for biodiversity to overcome fragmentation; (2) to facilitate the monitoring of changes in biodiversity at pan-European scales using harmonized methodologies; and (3) to assess the feasibility of the EMBOS system through pilot studies.
In a series of surveys during 2014 and 2015 the 40 members of EMBOS, representing 22 European countries, have measured at 28 marine sites along the European coastline (Baltic, Atlantic, Mediterranean) the degree and variation of the diversity and densities of hard and soft bottom communities using jointly agreed and harmonized protocols, tools and indicators. In this paper, we present an overview of the harmonized methods and tools used for sampling and analysis of soft-sediment macrozoobenthos along the European coastline, as well as the results of the first surveys on the geographic patterns of diversity. For the soft-sediment benthos the adopted hypothesis was that the diversity for all taxonomic groups would decrease with increasing latitude (see reviews on patterns and causes of the Latitudinal Diversity Gradient (LDG) in Stehli et al., Reference Stehli, McAlester and Helsley1967; Schopf et al., Reference Schopf, Fisher, Smith, Battaglia and Beardmore1978; Rohde, Reference Rohde1992; Roy et al., Reference Roy, Jablonski, Valentine and Rosenberg1998; Gaston, Reference Gaston2000; Willig et al., Reference Willig, Kaufman and Stevens2003; Hillebrand, Reference Hillebrand2004). The results of the surveys carried out in 2014 were used to examine whether this LDG could also be found from our results. Hillebrand (Reference Hillebrand2004) indicated that most studies on LDG comprise only few organism types and are often restricted to certain regions. Because of the harmonized set-up of the EMBOS Pilots we can address these problems and compare data on a wide range of organisms from a wide latitudinal range. Moreover, since we also included environmental data (temperature, salinity, wave height, organic and mud content of the sediment) in our analyses we can elucidate whether the trends in diversity might be explained by these factors regardless of latitude.
The results for the hard bottom benthic communities along the European coastline are reported in Kotta et al. (Reference Kotta, Orav-Kotta, Jänes, Hummel, Arvanitidis, van Avesaath, Bachelet, Benedetti-Cecchi, Bojanic, Como, Coppa, Coughlan, Crowe, Dal Bello, Degraer, Juanes de La Pena, Fernandes de Matos, Espinosa, Faulwetter, Frost, Guinda, Jankowska, Jourde, Kerckhof, Lavesque, Leclerc, Magni, Pavloudi, Pedrotti, Peleg, Pérez-Ruzafa, Puente, Ribeiro, Rilov, Rousou, Ruginis, Silva, Simon, Sousa-Pinto, Troncoso, Warzocha and Weslawski2016) focusing on relationships between cover, diversity and environmental factors, in Dal Bello et al. (Reference Dal Bello, Leclerc, Benedetti-Cecchi, Arvanitidis, van Avesaath, Bachelet, Bojanic, Como, Coppa, Coughlan, Crowe, Degraer, Espinosa, Faulwetter, Frost, Guinda, Jankowska, Jourde, Kerckhof, Kotta, Lavesque, de Lucia, Magni, Fernandes de Matos, Orav-Kotta, Pavloudi, Pedrotti, Peleg, Juanes de la Pena, Puente, Ribeiro, Rigaut-Jalabert, Rilov, Rousou, Rubal, Ruginis, Ruzafa, Silva, Simon, Sousa-Pinto, Troncoso, Warzocha, Weslawski and Hummel2016) with regard to scale-specific variability in community diversity and abundance, and in Puente et al. (Reference Puente, Guinda, Juanes de la Pena, Echavarri-Erasun, Ramos, de la Hoz, Degraer, Kerckhof, Bojanic, Rousou, Orav-Kotta, Kotta, Jourde, Pedrotti, Leclerc, Simon, Bachelet, Lavesque, Arvanitidis, Pavloudi, Faulwetter, Crowe, Coughlan, Benedetti-Cecchi, Dal Bello, Magni, Como, Coppa, de Lucia, Ruginis, Jankowska, Wesławski, Warzocha, Silva, Ribeiro, Fernandes de Matos, Sousa-Pinto, Troncoso, Peleg, Rilov., Espinosa, Pérez-Ruzafa, Frost, Hummel and van Avesaath2016) focusing on the role of physical variables in biodiversity patterns of intertidal macroalgae. The results for the functional diversity of the soft bottom benthic communities are reported in Pavloudi et al. (Reference Pavloudi, Faulwetter, Keklikoglou, Vasileiadou, Chatzinikolaou, Mavraki, Nikolopoulou, Bailly, Rousou, Kotta, Orav-Kotta, Bachelet, Lavesque, Benedetti-Cecchi, Dal Bello, Bojanic, Como, Coppa, Magni, Coughlan, Crowe, Degraer, Juanes de la Pena, Guinda, Puente, Fernandes de Matos, Ribeiro, Espinosa, Kerckhof, Jankowska, Weslawski, Peleg, Rilov, Pérez-Ruzafa, Ruginis, Jourde, Leclerc, Simon, Pedrotti, Silva, Sousa-Pinto, Rubal, Troncoso, Warzocha, van Avesaath, Frost, Hummel and Arvanitidis2016).
MATERIALS AND METHODS
Sampling procedure
To assess the trends in diversity in soft-sediment benthic communities along the European coastline, the individual densities of macrozoobenthos species were determined at 28 stations through field surveys from the south to the north of Europe (Figure 1; Table 1). The sampling stations covered a latitudinal cline extending from the Mediterranean, through the Eastern Atlantic and the North Sea to the Baltic (Figure 1). Sampling procedures and treatment and analyses of samples were harmonized for all participants as follows.
EMBOS soft substrate sampling locations.
Characteristics of the EMBOS soft substrate sample locations. Values of the Mud (% DW; <63 µm), Organic matter (%) content of sediment, and Salinity are averages of measurements at each site according to the EMBOS protocol. Values for wave height (m) and sea surface temperature (SST;°C) are the resultant of long-term Remote Sensing data interpretations for the coastal zone as specified in the Materials and methods.
Sampling was conducted in early spring, given that by being early into, or just before, the reproduction season, only low numbers of juveniles, if any, are present. Most samples were therefore taken in April, except for those from Crete and one location in Cyprus where samples were taken in May (Table 1). In intertidal areas, sampling was done at low tide, during the day, preferably at noon. In regions without tidal variation, samples were taken at the waterline during calm weather.
When multiple stations within a region were sampled, the stations had a distance of at least 2 km from each other (Figure 2). At each station, samples were taken at the lower intertidal (i.e. just above Low Water Level (LWL)) and/or upper subtidal (UST) level at a maximum depth of 2 m. At each level, three plots were chosen parallel to the LWL line, with a distance between plots of 100–200 m. At each plot, three replicate sediment samples with a depth of 30 cm, each at a distance of about 2 m, were taken with a hand-corer with a diameter of 13 cm (Figure 3) and sieved over a 0.5 mm mesh. The residues of the sieved samples were stained with Rose Bengal and preserved in 96% ethanol or 4% formaldehyde solution buffered with borax or hexamethylene tetramine. Thus, at each station a total of nine replicates per level were taken. Replicates at the three plots were taken in sediment of comparable granulometry, whereby the median grain size was 0.06–0.2 mm (sandy silt). The metadata describing the sampling campaign can be accessed at http://lifewww-00.her.hcmr.gr:8080/medobis/resource.do?r=embos_2014.
Harmonized EMBOS sampling scheme.
Standard EMBOS sampling corer for soft sediments.
Taxonomic analysis
The residues were sorted in the laboratory under a magnification lamp or stereomicroscope. The macrofauna taxa (>0.5 mm) were determined to the species level, when possible, according to up-to-date taxonomic literature, and their abundance recorded. The most optimal mesh-size for obtaining most macrofaunal species, and to have more realistic abundance values, while comparing different regions, is 0.5 mm (Rees, Reference Rees1984; Bishop & Hartley, Reference Bishop and Hartley1986; Bachelet, Reference Bachelet1990; Ferraro & Cole, Reference Ferraro and Cole2004; Couto et al., Reference Couto, Patricio, Neto, Ceia, Franco and Marques2010).
Species nomenclature followed the World Register of Marine Species (WoRMS) (Costello et al., Reference Costello, Bouchet, Boxshall, Fauchald, Gordon, Hoeksema, Poore, van Soest, Stöhr, Walter, Vanhoorne, Decock and Appeltans2013). Oligochaetes, turbellarians, sponges, nemerteans, insect larvae and meiofaunal taxa (mainly nematodes and copepods) were not identified to the species level but at the overarching group levels. Incomplete specimens were counted as single individuals if containing the head, whereas other body parts assigned to a given taxon were collectively counted as ‘1’ (one). Analysed samples were kept at the individual laboratories for (re)analyses, if needed, at a later stage.
Environmental factors
In order to analyse species diversity in relation to the main environmental factors, salinity was measured in the surrounding water with a CTD, and sediment temperature with the tip of a thermometer 1 cm below the surface. Moreover, additional sediment samples were taken for grain size and organic carbon content analyses. For this, two cores per plot were taken from the upper 5 cm top layer of the sediment with a 3 cm diameter corer. Samples were pooled per plot and mixed, resulting in three sediment samples per station. The samples were preserved at room temperature in aluminium foil or plastic bottles, after drying at 60°C for at least 24–48 h. Shell remains, which were abundant, were removed from the dried sediments with pincers before further sample processing. The sediment samples were centrally analysed by the Institute for Coastal Marine Environment in Torregrande, Oristano, Italy.
Sediment samples (~4 g) were pre-treated with H2O2 (20% volume) to eliminate organic material, washed with bi-distilled water to eliminate chlorides and then oven-dried at 40°C for 12 h. They were subsequently wet sieved in order to obtain the sand fraction (>63 µm). Samples were then pre-treated with Na-Hexametaphosphate 0.6% to avoid particle flocculation for 24 h and sonicated for 5 min before analysis. The analysis of mud content (<63 µm) was performed using a Galai CIS 1 laser instrument, with specific analytical size intervals of 0.5 µm (De Falco et al., Reference De Falco, Magni, Teräsvuori and Matteucci2004). The organic matter (OM) content in the sediments was determined from a subsample (~1 g) by loss on ignition (LOI) at 500°C for 3 h (Dean, Reference Dean1974).
An overview of the stations and the environmental factors used in the analysis is given in Table 1.
To estimate the variations of sea surface temperature (SST) and significant wave height (Hs), data from 1985 to 2013 were acquired by the Instituto de Hidráulica Ambiental de la Universidad de Cantabria (Fundación IH, Santander, Spain) from reanalysis sources at the nearest coastal point with information to the reference points. SST values were supplied, with daily temporal resolution, by the Operational Sea surface Temperature and sea-Ice concentration Analysis (OSTIA) dataset, which is under MyOcean2 project by UK-Met Office (NASA) (Stark et al., Reference Stark, Donlon, Martin and McCulloch2007). Specifically, the Group for High Resolution Sea Surface Temperature (GHRSST) L4 Gap-free gridded products have been used, with a spatial resolution of 0.05°. The wave data used in this work come from the Global Ocean Wave reanalysis database (GOW), reflecting wave height in the open coastal zone (Reguero et al., Reference Reguero, Menéndez, Méndez, Mínguez and Losada2012). Hourly significant wave height data were extracted with a spatial resolution of 0.125° for all sites. Salinity was obtained from in situ measurements provided by the World Ocean Database (WOD) of the National Oceanic and Atmospheric Administration (NOAA)-NESDIS National Oceanographic Data Centre (NODC) (Levitus et al., Reference Levitus, Antonov, Baranova, Boyer, Coleman, Garcia, Grodsky, Johnson, Locarnini, Mishonov, Reagan, Sazama, Seidov, Smolyar, Yarosh and Zweng2013). The salinity profiles used in this study were acquired between 1985 and 2014, and due to the sparse available data, for each station the salinity value was calculated as the average of all data points within a circle of 0.4° radius in the first 40 m around the reference points.
Statistical analysis
Differences in macrozoobenthic assemblages between regions and/or sub-regions were analysed using a non-metric multidimensional-scaling (nmMDS) ordination model based on the Bray–Curtis dissimilarity matrix (Clarke & Warwick, Reference Clarke and Warwick2001) calculated on the mean values among the three replicate samples within each plot. A one-way ANOSIM randomization/permutation test was used to check for the significance of differences among regions (i.e. Baltic, Atlantic and Mediterranean regions) in the ordination model (Clarke & Warwick, Reference Clarke and Warwick2001). A one-way ANOSIM was also used to check for the significance of differences among sub-regions (see Table 2A and 2B, respectively, for the identification of the regions and the subregions). All data were square-root transformed prior to the analyses, to minimize the effect of dominant species (Clarke & Warwick, Reference Clarke and Warwick2001). One of the plots from Cyprus (Softades-Larnaca1) was excluded from the analyses as no fauna was found in each of the three replicates. MDS 2D representations were considered acceptable when the stress factor did not transgress the value of 0.2 (Clarke & Warwick, Reference Clarke and Warwick2001). nmMDS and ANOSIM were conducted with the PRIMER v6.1.12 package (Clarke & Gorley, Reference Clarke and Gorley2006).
ANOSIM pairwise global tests on the differentiation between regions on basis of macrozoobenthic species and densities (2A – for the regions: global R = 0.439, significance level of sample statistic 0.1%; 2B – for the sub-regions: global R = 0.619, significance level of sample statistic 0.1%; n.s., not significant).
Hierarchical clustering and ordination methods were performed to identify, analyse and compare the multivariate pattern of soft sediment macrozoobenthic assemblages based on diversity and density descriptors (number of species (S), total density (n), Margalef species richness (d = (S−1)/ln(n)), Shannon diversity (H’ = −∑pi(ln(pi)) and Pielou evenness (J’ = (H’/ln(S)) per sample with pi being the proportion of species i in the sample) and environmental and geographic location characteristics. Multivariate analyses based on community descriptors were used because an initial analysis solely based on densities for each of the species resulted in poor relations with environmental descriptors, since the locations to compare at European scale each had a unique species composition with few species in common. A Principal Component Analysis (PCA), an indirect gradient analysis (i.e. ordination with optimal gradient distribution solely based on ‘species’ data, with plotting of environmental and location specific characteristics afterwards to identify potential relations), was performed on ‘ln(x + 1)’-transformed data, using the CANOCO for Windows version 4.5 software package (Ter Braak & Smilauer, Reference Ter Braak and Smilauer1998). A Detrended Correspondence Analysis (DCA) showed that the data had a short gradient length, which allowed the use of a linear ordination method. In the PCA analysis, environmental and location-specific characteristics with little explanatory power towards the observed patterns in species compositions, or highly co-varying with other characteristics, are not shown in the PCA results, as omitting (or adding) potential explanatory parameters does not have an influence on the ordination results in an indirect analysis.
The total macrozoobenthic densities (subdivided in abundances for larger taxonomic groups) were plotted in bar-graphs against the environmental factors that showed the best explanatory power in the PCA analysis (subdivided into classes from low to high values). Significant differences in total densities between identified classes were tested by using two-sided t-tests after testing for similarity of variances (F-tests) in Microsoft Excel 2010.
Additional graphs were plotted identifying relations of Shannon diversity (H’ log e) and total macrozoobenthos densities (the average of 3 replicates per sampling plot) with the environmental parameters. Because diversity indicators such as S, d, J’ and H’ are highly correlated, as can be expected since they are calculated on the basis of the same elements, the focus was on Shannon diversity (H’) as community descriptor. Best fitting regressions were calculated starting from linear regression through a quadratic and cubic, to maximally a quartic (polynomial) regression, assuming that r2 should increase at least 10% in order to adopt a higher polynomial level as being better fitting (r 2 will always increase at a higher-order term, but it is hard to imagine a biological meaning for exponents greater than 3) (McDonald, Reference McDonald2014).
RESULTS
Relations of major taxonomic groups with latitude and region
Multivariate analyses (MDS and ANOSIM) revealed differences in the densities-based species composition of macrozoobenthic assemblages among the Baltic, Atlantic and Mediterranean regions (Figure 4A; Table 2A).
MDS analysis of the on densities based macrozoobenthos species compositions for the major European regions (A) and their sub-regions (B).
We also found a significant differentiation between most sub-regions as indicated by the pairwise test (Figure 4B; Table 2B). In particular, all Mediterranean sub-regions (i.e. Eastern Mediterranean, Sea of Crete, Ionian Sea and Western Mediterranean) were significantly different from each other (Table 2), with the highest R value when comparing the Eastern with the Western sub-region (R = 0.879, P < 0.001; Table 2B).
In contrast, most Atlantic sub-regions (i.e. Atlantic Ocean, English Channel and North Sea) did not significantly differ (P > 0.10; Table 2B), except for the Bay of Biscay which is different from all other sub-regions.
The distribution of the mean densities of the major taxonomic groups (expressed in absolute and relative terms) varied with latitude (Figure 5), reflecting the major differentiation among regions and sub-regions as revealed by multivariate analyses. The following major trends were found: (1) at low latitudes (=Eastern Mediterranean at 32–36°N): very low densities of all groups; (2) at somewhat higher latitudes (=Western Mediterranean >36°N): high densities of mainly Malacostraca (crustaceans); (3) middle latitudes (=Atlantic and North Sea): intermediate densities, with a dominance of Polychaetes; (4) high latitudes (=Baltic): intermediate densities with mainly insects.
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with latitude. Average total densities ± standard error for clusters in steps of 4° latitude; significant differences in total densities are indicated with different characters.
Relations of major taxonomic groups with environmental factors and location characteristics
The densities of all major taxonomic groups together (i.e. total densities) did not show a consistent pattern in relation to salinity (Figure 6A). Yet, for each separate taxonomic group a specific relation with salinity was found (Figure 6B). Both in absolute as well as in relative terms, insects were mainly found at low salinities, and crustaceans (Malacostraca) at higher salinities. Polychaetes were found equally at all salinities, whereas bivalves and gastropods were not present at the lowest salinities. These results corroborate the above-described results on the distribution of taxa by region, as low salinities occur in the Baltic, where mainly insects were found, and high salinities are found in the Mediterranean favouring the crustaceans.
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with salinity. Average total densities ± standard error for clusters according to the Thalassic salinity series (0–5 = oligohaline, 5–18 = mesohaline, 18–30 = polyhaline, 30–40 = mixoeuhaline, >40 = metahaline); significant differences in total densities are indicated with different characters.
Densities of major taxonomic groups were higher at higher temperatures (Figure 7A). The insects, mainly occurring in the Baltic (Figure 7B), were found at low temperatures (as do occur in the Baltic in early spring). The crustaceans (Malacostraca), and especially the bivalves, were found mainly at higher temperatures, while polychaetes were found equally at all temperatures.
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with year average sea surface temperature (°C). Average total densities ± standard error for clusters in steps of 4°C, except for the lowest temperature class that covers 5°C; significant differences in total densities are indicated with different characters.
Densities of all taxonomic groups increased with increasing mud content (Figure 8A) as well as with increasing organic matter content (Figure 9A). Values levelled off at a mud content percentage higher than 12–24%.
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with mud (<63 µm) content (%). Average total densities ± standard error for clusters of times 2 increasing mud contents; significant differences in total densities are indicated with different characters.
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with organic matter content (% DW). Average total densities ± standard error for clusters of times 2 increasing organic matter contents; significant differences in total densities are indicated with different characters.
All taxonomic groups showed a relatively similar increase of densities at higher mud or organic matter content, so the relative distribution of taxonomic groups at different levels of mud or organic matter remained rather similar (Figures 8B & 9B).
The total density was highest at the lowest wave height (Figure 10A). That pattern was also observed for the bivalves, whereas the density of polychaetes showed an opposite trend with highest densities at the highest wave heights (Figure 10B).
Absolute (A) and relative (B) numbers of individuals per major taxonomic group with average wave height (m). Average total densities ± standard error for clusters of increasing wave height; significant differences in total densities are indicated with different characters.
The total density was on average two times higher in lagoons than in estuaries or the open coast – the last two having similar densities (Figure 11A). The densities were higher at subtidal than at intertidal stations (Figure 11A). At intertidal stations in estuaries the polychaetes are more dominant than in other stations, and in the open coast stations the crustaceans dominate (Figure 11B).
Absolute (A) and relative (B) numbers of individuals per major taxonomic group for each of the identified system types and strata where sampling has taken place. Average total densities ± standard error; significant differences in total densities are indicated with different characters where both classification systems (types and strata) have been tested separately as they consist of the same samples.
Variation in community diversity and densities, and environmental variables
The PCA, based on the diversity and abundances of benthic communities, showed that the (variance in) environmental and geographic location characteristics coincided strongly with the (variance in) macrozoobenthos diversity (Shannon H’) and densities (>90%; Table 3). Regarding geographic location most variance coincided with that of the main geographic regions Atlantic (positively correlated to that of diversity) and Baltic (negatively correlated) (Figure 12). Moreover, the (variance in) diversity showed a strong positive correlation with that of the (variance in) major environmental factors and location characteristics (salinity, temperature, mud content, organic matter content, wave height, the estuarine system type and intertidal stratum) (Figure 12, Table 3). At the other hand, the (variance in) diversity was negatively correlated with that of the factor ‘open coast’, which corresponds to the much lower diversity at open coasts than at the other locations (in estuaries or lagoons; Figure 13).
Results of a Principal Component Analysis (PCA) based on macrozoobenthic community descriptors (diversity and densities descriptors). Most important correlated environmental (continuous variables) and geographic (presence/absence data) characteristics are plotted as dashed vectors afterwards (as a PCA is an indirect gradient analysis) (Wave h = wave height).
Shannon diversity (average ± standard deviation) of macrozoobenthos for each of the identified system types and strata where sampling has taken place. Significant differences in Shannon diversity are indicated with different characters; both classification systems (types and strata) have been tested separately as they consist of the same samples.
Results of Principal Component Analysis (PCA) between the diversity and density based macrozoobenthic community descriptors (number of species (S), total density (n), Margalef species richness (d), Pielou evenness (J’) and Shannon diversity (H’) per sample) and environmental and geographic location characteristics (characteristics with little explanatory power towards the observed patterns in species compositions are not indicated; to these belong the environmental variables ‘CaCO3 content’ and ‘sampling depth’).
The variance in densities was only marginally explained by the measured environmental factors (Figure 12). Yet, the location characteristics ‘lagoon’ and ‘open coast’ were respectively strongly positively and negatively correlated with the (variance in) densities. Clearly the densities at open coast stations are much lower than in lagoons (Figure 13). The lower explanatory level of the measured environmental factors for the variance in densities, in comparison to that for diversity, also became clear from an analysis of the relationship between the individual factors, as indicated in the following section.
Relation of benthic diversity and total densities with environmental factors
Benthic diversity was highest at latitudes from 40 to 45°N (Figure 14A). The highest overall densities (including all species) were found at 40°N (southern Atlantic Ocean; Bay of Biscay) in comparison to lower as well as higher latitudes (Figure 14B). The highest soft-sediment benthic diversity and densities were found in areas with annual average SST of 15 to 20°C (Figure 15A, B). The curves of the diversity and density variation with increasing temperature (Figure 15) were the reverse (mirror images) of those found for latitude (Figure 14), as was to be expected since temperature is generally inversely related to latitude. The lowest diversity values were found at salinities below 10, whereas the higher diversities were mainly found above 15 ppt (Figure 16A). Although no clear trend was found for densities, it appears that high densities occurred primarily at salinities above 20 ppt (Figure 16B).
Diversity (A; Shannon H’; 3rd order polynomial, r 2 = 0.48, P < 0.001) and density (B; total number per m2; 4th order polynomial, r 2 = 0.34, P < 0.001) of macrozoobenthos with latitude.
Diversity (A; Shannon H’; 3rd order polynomial, r 2 = 0.39, P < 0.001) and density (B; total number per m2; 4th order polynomial, r 2 = 0.21, P < 0.001) of macrozoobenthos with sea surface temperature (SST, °C).
Diversity (A; Shannon H’; 4th order polynomial, r 2 = 0.44, P < 0.001) and density (B; total number per m2; 4th order polynomial, r 2 = 0.16, P < 0.01) of macrozoobenthos with salinity.
Diversity and densities sharply increased with mud content between 0–10% (Figure 17A, B) but slightly decreased at the highest levels (>40%). A similar relation was found between organic matter content and diversity (Figure 18A), while no decrease at higher levels was found for the densities (Figure 18B).
Diversity (A; Shannon H’; 2nd order polynomial, r 2 = 0.32, P < 0.001) and density (B; total number per m2, 2nd order polynomial. r 2 = 0.20, P < 0.001) of macrozoobenthos with mud content (%DW).
Diversity (A; Shannon H’; 2nd order polynomial, r 2 = 0.26, P < 0.001) and density (B; total number per m2, 3rd order polynomial. r 2 = 0.21, P < 0.001) of macrozoobenthos with organic matter content (%DW).
The diversity of macrozoobenthos increased at increasing average wave heights from 0.4–1.4 m in coastal areas (Figure 19A) and decreased again at higher wave heights. For the densities no clear relation with wave height could be found (Figure 19B).
Diversity (A; Shannon H’; 3rd order polynomial, r2 = 0.41, P < 0.001) and density (b; total number per m2, 4th order polynomial. r 2 = 0.37, P < 0.001) of macrozoobenthos with average wave height (%DW).
DISCUSSION
Patterns in diversity are often related to latitude, a phenomenon known as the latitudinal diversity gradient (LDG), whereby a decrease with increasing latitude is found (Roy et al., Reference Roy, Jablonski, Valentine and Rosenberg1998; Rex et al., Reference Rex, Stuart and Coyne2000; Attrill et al., Reference Attrill, Stafford and Rowden2001; Willig et al., Reference Willig, Kaufman and Stevens2003; Hillebrand, Reference Hillebrand2004). Several hypotheses for the underlying causes for such pattern are suggested, yet none of them is solely sufficiently convincing (Willig et al., Reference Willig, Kaufman and Stevens2003; Hillebrand, Reference Hillebrand2004), although solar energy input (and for the marine territory the sea surface temperature as its proxy) is most often mentioned as the main acting principle (Rohde, Reference Rohde1992; Roy et al., Reference Roy, Jablonski, Valentine and Rosenberg1998). In corroboration to this rule we indeed found a positive relation between temperature and diversity. However, the curve with temperature is bell-shaped, i.e. at higher temperatures the diversity decreases again. This decrease could be related to the physiological limit at higher temperatures in marine invertebrates (Newell, Reference Newell1979), a situation that occurred in our study in the South-Eastern Mediterranean.
This means also that in our study a low diversity at the lower (Mediterranean) latitudes was found. Therefore, in contrast to the general LDG pattern, in our study for benthic species of soft substrates the highest diversity (as well as densities) was found at latitudes around 40–50°N, while lower diversity (and densities) was found at lower and higher latitudes. Deviations from the general LDG pattern have also been reported before for European marine benthos by Renaud et al. (Reference Renaud, Webb, Bjørgesæter, Karakassis, Kędra, Kendall, Labrune, Lampadariou, Somerfield, Włodarska-Kowalczuk, Van den Berghe, Claus, Aleffi, Amouroux, Bryne, Cochrane, Dahle, Degraer, Denisenko, Deprez, Dounas, Fleischer, Gil, Grémare, Janas, Mackie, Palerud, Rumohr, Sardá, Speybroeck, Taboada, Van Hoey, Węsławski, Whomersley and Zettler2009) who found no or weakly positive relationships. The explanation for this kind of diverting trend is that the impact of (local variation in) environmental factors is stronger than that of latitude related factors (Gaston, Reference Gaston2000; Renaud et al., Reference Renaud, Webb, Bjørgesæter, Karakassis, Kędra, Kendall, Labrune, Lampadariou, Somerfield, Włodarska-Kowalczuk, Van den Berghe, Claus, Aleffi, Amouroux, Bryne, Cochrane, Dahle, Degraer, Denisenko, Deprez, Dounas, Fleischer, Gil, Grémare, Janas, Mackie, Palerud, Rumohr, Sardá, Speybroeck, Taboada, Van Hoey, Węsławski, Whomersley and Zettler2009).
One of the factors not related to latitude, yet influencing the diversity, is salinity. A higher diversity with higher average salinity, as observed in our dataset, is in line with the classical Remane pattern of larger species richness at marine conditions compared with brackish conditions (Remane, Reference Remane1934; Attrill, Reference Attrill2002). A slight increase of diversity near freshwater conditions was also found in our results. However, in contrast to the classical pattern, a depression in diversity was also found at salinities around 35–40 ppt, which is just above the regular seawater salinity. These higher salinities occurred in the South-Eastern Mediterranean and thus may help to explain the lower diversity at these latitudes.
Other environmental factors relating to diversity and density, that could also be causative for deviations from the LDG, were organic matter and mud content of the sediment and wave height. An increase of all these factors resulted in higher diversity, up to a certain maximum or threshold level above which the diversity decreased. The decrease at high levels, resulting in a bell-shaped curve, may be because sediments with e.g. high mud and/or organic content could be too swampy or anoxic and thus unsuitable for benthic animals (Pearson & Rosenberg, Reference Pearson and Rosenberg1978; Hyland et al., Reference Hyland, Balthis, Karakassis, Magni, Petrov, Shine, Vestergaard and Warwick2005; Magni et al., Reference Magni, De Falco, Como, Casu, Floris, Petrov, Castelli and Perilli2008).
In addition to latitudinal variation, the longitudinal differentiation also should be taken into account. This can be explained by the variation within regions, as the lowest salinities can be found towards the east in the Baltic (Bonsdorff, Reference Bonsdorff2006), and the lowest wave heights and highest temperatures in the east of the Mediterranean (Lejeusne et al., Reference Lejeusne, Chevaldonne, Pergent-Martini, Boudouresque and Perez2010). Each of these ‘extreme’ environmental conditions, coinciding with a lower diversity, exists in the eastern part for both regions (the Baltic and Mediterranean), thereby giving the result of the observed strong axis for longitude in the multivariate analysis. Moreover, the lower diversity in the east of the Mediterranean might also be related to lower primary production in the east compared with the west (Bricaud et al., Reference Bricaud, Bosc and Antoine2002), whereby the carrying capacity of the Eastern Mediterranean might be lower.
At an even smaller (sub-regional) scale, a remarkable result was the difference of the Bay of Biscay compared with the, among each other similar, NE Atlantic, North Sea and English Channel sub-regions. In the Bay of Biscay, SST is known to be higher than in the surrounding areas connected by the Gulf Stream (Jenkins et al., Reference Jenkins, Moore, Burrows, Garbary, Hawkins, Ingolfsson, Sebens, Snelgrove, Wethey and Woodin2008), which may explain the much higher community resemblance to the Moroccan coast and Mediterranean than the adjacent sub-regions (Jenkins et al., Reference Jenkins, Moore, Burrows, Garbary, Hawkins, Ingolfsson, Sebens, Snelgrove, Wethey and Woodin2008), as observed also for macroalgae (Sauvageau, Reference Sauvageau1897; Fischer-Piètte, Reference Fischer-Piètte1955).
In addition to the correlation with environmental factors and (sub)regions a relation was found with location characteristics. The higher diversity in estuaries and lagoons, than at open coast, corroborates other studies and is normally attributed to the higher diversity of habitats (ecotones) and the higher production in those systems (De Wit, Reference De Wit, Gillo and Venora2011; Miththapala, Reference Miththapala2013). The higher densities in lagoons are also corroborated by these authors, and attributed to the edge-effect, i.e. higher numbers at locations where different habitats/ecotones meet each other.
However, some of the above mentioned environmental as well as geographic influencing factors may coincide and cannot be fully teased apart with the current sampling scheme, as e.g. latitude and temperature, or salinity and longitudinal location in the Baltic. Further integrated studies, and thus an extension of the harmonized EMBOS approach and an even wider geographic coverage, may be needed to overcome this kind of multicollinearity and to reach a stronger conclusion.
In summary, the differentiation in diversity and densities of macrozoobenthos and in environmental factors coincides with strong interregional differentiation, especially between the Atlantic and the Baltic, as seen in the PCA analysis. In the Atlantic, much higher values of salinity, temperature, mud and organic matter content in the sediments, as well as wave heights are observed. Therefore, it can be suggested that the latitudinal and regional differences found in densities and diversity are strongly determined by the environmental differentiation between the Baltic (low salinities, no tide), the Atlantic (high mud and organic matter content in sediment, moderate temperature and salinity, high average wave heights) and the Mediterranean (high temperatures, high salinity, small tide).
Consequently, some soft substrata taxa are numerically more dominant in one region or the other, e.g. crustaceans in the Mediterranean, polychaetes in the Atlantic and insects in the northern Baltic. Thus, it may be argued that latitudinal trends and regional differences in diversity and densities are merely a result of including very different bodies of water in this survey, which is unlike testing a continuous gradient along a single coast line; i.e. the two extremes, the Baltic which is a semi-enclosed sea with very low salinities favouring a high proportion of insects (Bonsdorff, Reference Bonsdorff2006), vs. the Mediterranean, also a semi-enclosed sea yet with relatively high salinities and high temperatures favouring a high proportion of crustaceans. However, what should be taken into account is that these two enclosed bodies of water also have very different histories, the former being much younger than the latter, which can also affect diversity levels through evolutionary processes (Bonsdorff, Reference Bonsdorff2006; Nordström et al., Reference Nordström, Lindblad, Aarnio and Bonsdorff2010).
Taking all results together, it can be concluded that the hypothesis expecting a negative linear gradient in diversity with latitude was not supported, as the highest observed diversity was observed at latitudes between 40 and 45°N. This bell-shaped relationship is most probably driven by large differences in environmental factors that are not merely determined by latitude.
However, the latitudinal range of this study is restricted to 30 to 60°N, which should be preferably extended over the whole latitudinal gradient in order to obtain a better understanding of this hypothesis.
Moreover, the accuracy of the observed patterns between diversity or densities and environmental parameters can be improved by using in situ measurements instead of remote sensing data. In this study remote sensing data were probably not always sufficiently detailed to get accurate values for semi-enclosed areas, whereby the relationships with diversity and densities were not optimal. Further, to improve the resolution of relations between environmental factors and densities, seasonality should be taken into account, since in invertebrates the densities normally show strong annual fluctuations and thus much weaker relationships with environmental factors than will be (and thus were) found for (the less fluctuating) diversity.
CONCLUSION
Conclusions of this study are:
-
– By employing harmonized tools and methodologies, the EMBOS system delivered an extensive comparable dataset on the diversity and densities of the soft-sediment macrozoobenthic community over a large-scale gradient along the European coastline.
-
– Species diversity has no linear (negative) relationship with latitude, yet in deviation to general theory, a bell-shaped relation is found in the studied latitudinal range.
-
– In general, the diversity or densities of soft-sediment macrozoobenthos were positively related with environmental factors such as temperature, salinity, mud and organic matter content in sediment, and wave height. For some relationships an optimum curve (e.g. temperature from 15–20°C; mud content of sediment around 40%) or bimodal curve (e.g. salinity) was found. Densities were also higher at lagoon and subtidal stations (versus estuarine, open coastal and intertidal stations).
-
– Latitudinal trends and regional differences in diversity and densities are merely the result of including areas with specific sets and ranges of environmental factors, such as the Baltic with typical salinity clines (favouring insects) and the Mediterranean with higher temperatures (favouring crustaceans). Therefore, putative latitudinal trends are indirect.
FINANCIAL SUPPORT AND ACKNOWLEDGEMENTS
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 and Helen Orav-Kotta were supported by Institutional research funding IUT02-20 of the Estonian Research Council and the BONUS project BIOC3, 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 Estonian Research Council. Pedro Ribeiro and Puri Veiga were supported by the Portuguese Science Foundation (FCT) through a post-doctoral grant (ref. SFRH/BPD/69232/2010 and SFRH/BPD/81582/2011, respectively). Valentina de Matos was supported by the Portuguese Science Foundation (FCT) through a doctoral grant (ref. SFRH/BD/86390/2012). 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’). The authors want to thank Dimitra Mavraki, Matina Nikolopoulou, Nicolas Bailly (of HCMR), Alina Sousa, André Costa, Cristina Santo, João Castro, Marta Mamede, Nuno Mamede, Susana Celestino (of MARE-UE), Hugues Blanchet, Xavier de Montaudouin, Benoît Gouillieux, Michel Leconte (of Arcachon Marine Station), Laureen Beaugeard and Clément Chaigneau (of University of La Rochelle), Concepción Marcos, Jessica Titocci, Antonio Sala Mirete and Jose Gabriel Hernández (of University of Murcia) for support given during sampling and analyses. Emilia Jankowska and Jan Marcin Węsławski were supported by the Statutory Funds of the Institute of Oceanology, Polish Academy of Sciences. A. Satta, CNR-IAMC Oristano, is acknowledged for sediment analysis. Jan Warzocha, Sławomira Gromisz and Bartosz Witalis (NMFRI) were supported by Statutory Funds from the Polish Ministry of Science and Higher Education (ref. P1-3/2012-2014).
