Journal of Marine Sciences

Journal of Marine Sciences / 2014 / Article

Research Article | Open Access

Volume 2014 |Article ID 462529 |

Thijs Christiaan van Son, Rune Halvorsen, Karl Norling, Torgeir Bakke, Maria Kaurin, Fredrik Melsom, "Identification of Fine-Scale Marine Benthic Ecoclines by Multiple Parallel Ordination", Journal of Marine Sciences, vol. 2014, Article ID 462529, 23 pages, 2014.

Identification of Fine-Scale Marine Benthic Ecoclines by Multiple Parallel Ordination

Academic Editor: Horst Felbeck
Received09 Jun 2013
Revised12 Nov 2013
Accepted13 Nov 2013
Published16 Feb 2014


The species-environment relationship is a fundamental structural property of natural ecosystems. Marine sedimentary macrofauna is known to be structured by a range of environmental variables; however, the mechanisms by which environmental variables covary to form complex-gradients (i.e., groups of intercorrelated environmental variables), and how these are related to coenoclines (i.e., gradients in species composition), remain poorly understood. We classified our study area into geomorphological features that were used for stratified sampling of macrofaunal polychaetes, molluscs, and echinoderms. The resulting species-by-site matrix was subjected to indirect gradient analysis by a multiple parallel ordination strategy, using detrended correspondence analysis and nonmetric multidimensional scaling. One major and one minor coenocline were identified. Based on the correlation between complex-gradients and the main coenocline we hypothesise the existence of two ecoclines that we have termed Periodic hypoxia and Periodic physical forcing. We conclude that a combination of recurrent (periodical) and extreme events is likely to determine the variation found in the species composition of marine sedimentary ecosystems. Based on the results of our study, we conclude that indirect gradient analysis is a useful tool for enhancement of our basic mechanistic understanding of the processes governing the compositional structure of marine sediment communities.

1. Introduction

The species-environment relationship is a fundamental structural property of natural ecosystems. Several studies of this relationship in marine sedimentary systems demonstrate that marine benthic macrofauna is structured by a range of environmental variables such as sediment particle-size composition [1, 2], organic matter content [35], salinity [6, 7], temperature [810], and exposure to current or wave energy [11, 12].

Basic knowledge of species-environment relationships can be obtained from studies with a general purpose, that is, to summarise the main structure of a species-by-site matrix, to relate this structure to environmental variables, and to generate hypotheses about species-environment relationships [13, 14]. Among plant ecologists, but to a much lesser degree among marine ecologists, this is accomplished by application of indirect gradient analysis in a three-step procedure: (1) to summarise the main coenocline structure of a species-by-site matrix by use of unconstrained ordination methods such as detrended correspondence analysis (DCA) or nonmetric multidimensional scaling (NMDS). The extracted ordination axes represent the major gradients in species composition, that is, coenoclines [14]. A species-by-site matrix typically contains between one and four coenoclines [15]. (2) To relate recorded environmental variables to the extracted coenoclines. In this step, groups of correlated variables that covary with the coenoclines, that is, complex-gradients, are identified. (3) To use patterns of covariation identified in Step 2 to hypothesise about ecoclines [14], that is, gradients in species composition conditioned by environmental variation. It is important to note that both the ecoclines and the major complex-gradients [15] are abstract concepts or hypotheses resulting from inductive methods [16] which subsequently can be tested using hypothetico-deductive methods.

While plant ecologists have embraced indirect gradient analysis as a most valuable tool for exploring species-environment relationships, marine ecologists, especially marine benthic ecologists, have generally had a much more applied focus in their research. For this applied purpose, the PRIMER [17] software has been a valuable tool. However, PRIMER is not well suited for studies of indirect gradient analysis with the aim to identify coenoclines, complex-gradients, and ecoclines. One reason for this is that the axes of NMDS ordinations as implemented in PRIMER, are unscaled and thus not interpretable in terms of compositional turnover. Such interpretation is essential for identification of coenoclines, for relating these to complex-gradients, and for hypothesising about ecoclines [18].

The main aim of this study is to describe species-environment relationships of benthic invertebrates in a temperate fjord system and, implicitly, to address the application of indirect gradient analysis in marine benthic community studies. This is done by applying multiple parallel ordination (MPO [19]) to a data set consisting of records of polychaetes, molluscs, and echinoderms from stratified sampling in the Oslofjord, SE Norway. Multiple parallel ordination implies that two fundamentally different ordination methods are applied in parallel to the same set of (multiple) data and results checked for similarity of the identified compositional structure. The multiple data sets are derived from one original species-by-site matrix by differential weighting of species’ abundance. The MPO approach opens for examination of both qualitative and quantitative properties of the species-by-site matrix. To identify complex-gradients that are related to the coenoclines, a large set of potentially relevant environmental variables were recorded covering ocean current, terrain, and sediment profile properties, sediment particle composition, and organic content.

2. Material and Methods

2.1. Study Area

The study took place in the western part of the Inner Oslofjord (commonly referred to as the Vestfjord; see Figure 1) and covered an area of 1.05 km2 (UTM32 bounding coordinates; NW corner: 6631500, 584725; SE corner: 6630500, 585775). The study area was chosen because it covers a wide range of variation along several putatively important environmental gradients, because it is known to have high habitat diversity, and because it is situated at a sufficient distance from the centre of Oslo to avoid strong anthropogenic impact (freshwater runoff, supplies of contaminants, nutrients, and organic material). The wastewater plant VEAS is, however, located ca. 1200 m from the southern edge of the study area.

2.2. Bathymetry

In 2004 and 2005, the Geological Survey of Norway performed bathymetric surveys covering the whole inner part of the Oslofjord [20]. Onboard the research vessel (R/V) Seisma, an interferometric sonar system (GeoSwath, Geo-Acoustics) with 125 and 250 kHz transducers was used to provide a three-column xyz-file with georeferenced depth values. The xyz-file had a resolution of 1 m and was imported into GRASS 6.4.1 [21] as a gridded bathymetry raster using the command. GRASS was used for all GIS analyses. In order to reduce the noise near the swath edges of the interferometric survey, the gridded bathymetry raster was smoothed by application of the mode option in the method parameter of the r.neighbors command. A square smoothing window of 5 pixels × 5 pixels with the focal cell in the center was applied to each cell of the bathymetry raster [22].

2.3. Selection and Location of Stations

Sampling sites were selected and located based on a stratified sampling strategy, which, in turn, was based on a benthic terrain classification algorithm developed by Lundblad et al. [23]. This algorithm used fine- and broad-scaled rasters of slope and bathymetric position index (BPI) to classify the seabed into fine- and broad-scaled geomorphological features.

Given the size of the research vessel, the ability to keep the vessel in position, and the sampling depth, the expected radius of precision of sampling was estimated to be about 15 m. We therefore used moving windows of 33 pixels × 33 pixels (33 × 33 m) and 65 pixels × 65 pixels ( m), respectively, to obtain fine-scaled and broad-scaled input rasters of slope and BPI for use in the terrain classification algorithm.

Slope is a measure of terrain steepness, formally defined as the first-order derivative of the gridded bathymetry raster. Slope measures the magnitude of the strongest gradient from a raster cell to any of its neighbours. Fine- and broad-scaled slope were derived using the r.param.scale command of GRASS.

The BPI, the marine version of the topographic position index (TPI) [22, 23], is defined as the second-order derivative (the first derivative of slope) of the gridded bathymetry raster. BPI quantifies the rate of change of the slope and indicates the position of a cell relative to its neighbours in the marine landscape. Cells that form part of a crest feature, situated higher than their neighbouring cells, obtain positive values, while negative values are obtained for cells in depressions [22, 23]. Fine- and broad-scale BPI rasters, BPI33 and BPI65, respectively, were derived using a combination of the r.neighbors and r.mapcalc commands. By the former (with parameters: method = “average”, size = 33 or 65, and circular neighborhood), the focalmean33 and focalmean65 rasters were obtained, which were used by the latter command to obtain BPI33 = bathymetry − focalmean33 and BPI65 = bathymetry − focalmean65.

A slightly modified version of the algorithm by Lundblad et al. [23] was used to obtain a classification of the geomorphology of the study area. Our algorithm also took into account the mean depth for geomorphological features covering large proportions of the total study area. This was done by dividing the area covered by slopes and flats into shallow and deep sections (depth was converted into a categorical variable with two classes, separated at mean depth = 58.6 m). Just like Lundblad et al. [23], two classifications of geomorphological features were obtained; a fine-scale classification into a maximum of 14 classes that was used in stratified sampling to ensure that the full range of topographic variation in the study area was represented, and a broad-scale classification that is useful because it represents a level of generality that it is suitable for communicating relationships between faunal community structure and geomorphology. The two classifications were not fully nested (Figure 2; see Table 1 for an overview of relationships between fine-scaled and broad-scaled features, and for selection criteria for each feature).

GeomorphologyCover (%)SitesBPI65BPI33Slope65Slope33Depth

Fine-scaled features
 Narrow depression 6.41-2, 26≤−1 SD≤–1 SD
 Broad depression with open bottom 3.23-4≤−1 SDw/in 1 SD
 Narrow crest7.65-6, 25≥1 SD≥1 SD
 Local depression on flat0.17w/in 1 SD≤–1 SD≤5
 Lateral midslope depression1.48-9w/in 1 SD≤–1 SD>5
 Local crest on flat0.110w/in 1 SD≥1 SD≤5
 Lateral midslope crest 1.211-12w/in 1 SD≥1 SD>5
 Shallow flat 28.413–15, 27>−1 SDw/in 1 SD≤5≤58.6
 Deep flat 15.716–18>−1 SDw/in 1 SD≤5>58.6
 Shallow slope 21.919–21, 28>−1 SDw/in 1 SD>5 and ≤70≤58.6
 Deep slope 14.022–24>−1 SDw/in 1 SD>5 and ≤70>58.6
Broad-scaled features
 Depression 9.61–4, 26≤−1 SD
 Crest 11.15-6, 25≥1 SD
 Deep flat 16.17, 11, 16–18 w/in 1 SD≤5>58.6
 Shallow slope 19.68-9, 19–21, 28w/in 1 SD>5≤58.6
 Shallow flat 29.410, 13–15, 27w/in 1 SD≤5≤58.6
 Deep slope 14.212, 22–24 w/in 1 SD>5>58.6
 Total 100.0

Cover: estimated percentage of the study area covered by the feature; Sites: number of the sampled sites belonging to each fine- and broad-scaled feature class; BPI65 and BPI33: bathymetric position index rasters calculated from 65 and 33 m moving windows; Slope65 and Slope33: corresponding rasters for slope. The mean depth of the study, 58.6 m, was used to differentiate between shallow or deep slopes and flats, respectively.

Of the 14 fine-scaled geomorphological features recognised by the classification scheme, 11 were present in the study area (Table 1), while three, local crest in depression, depression on crest, and steep slope, were not found. Flats and slopes dominated the study area, while narrow depression and narrow crest each covered more than 5% of the total area. Some features such as local depression on flat, lateral midslope depression, local crest on flat, and lateral midslope crest each covered less than 1.5% and formed small patches only.

A total of 28 sites were selected for sampling. Each of the 11 fine-scaled and 6 broad-scaled geomorphological features encountered in the study area were represented by at least one site (see Table 1 and Figure 2). To accomplish this, 24 of the sites were distributed among the fine-scaled features according to the following criteria: features covering less than 1% of the total area were sampled once, features covering 1–10% of the area were sampled twice, while features covering more than 10% were sampled three times. Finally, one site was placed randomly in each of the broad-scaled features crest, depression, flats, and slopes.

2.4. Faunal Samples

Benthic macrofauna from the 28 sites was sampled from the R/V Trygve Braarud during March and April 2009. At each site, four samples were collected using a 0.1 m² van Veen grab. Each sample was washed through 5 mm and 1 mm sieves on board the research vessel. The residue (sediment and fauna) from these two sieves was fixed in ethanol and stained using Rose Bengal stain. In the laboratory the fauna was separated from the remaining sediment, sorted into main taxonomic groups (Polychaeta, Mollusca, Echinodermata, Crustacea, and Varia), and subsequently conserved in 70% ethanol. All polychaetes, molluscs, and echinoderms were identified to the species level or to the lowest taxon possible and the number of individuals of each taxon was counted. The total number of individuals in each sample was used as measure of abundance. The identified taxa represent the most species-rich taxonomic groups in the study area and are likely to represent more than 90% of the total species richness of the sampled fauna.

2.5. Environmental Variables

A total of 23 environmental variables were recorded for all 28 sites (see Table 2). From the bathymetry raster grid provided by the Geological Survey of Norway, we extracted aspect, BPI, depth, maximum curvature of the seabed, and slope using the fine-scaled moving window (33 × 33 pixels). Gridded ocean current rasters with 15 m resolution (A. Staalstrøm, unpublished observations) were used to derive variables for maximum speed and direction of surface and bottom currents. The rate of instantaneous change in the current speed was obtained by estimating the profile curvature [24] of a raster surface representation of the ocean current, applying the r.param.scale command (parameters: method = “profc” and size = 15) to a moving window of 15 pixels × 15 pixels (i.e., 225 m × 225 m).


DEM derived
 Aspect_EEastness component of aspect−1.00–1.000.10
 Aspect_NNorthness component of aspect−0.96–0.92–0.35
 BPIBathymetric position index−1.52–1.34–0.06
 Depth [m]Depth derived from DEM35.8–93.258.6
 Seabed_curvMaximum curvature of the seabed (in any direction)−0.011–0.0420.005
 Slope [deg]The slope of the seabed0.3–19.86.7
Current model
 Current_bottom [m/s]Maximum current speed at the bottom0.026–0.0770.049
 Current_surface [m/s]Maximum current speed at the ocean surface0.061–0.0910.068
 Change_curr_bottomThe rate of instantaneous change of max bottom current−1.71e-3–1.18e-31.96e-6
 Change_curr_surfaceThe rate of instantaneous change of max surface current−1.18e-3–3.61e-4−1.92e-4
 Direction_curr_EEastness component of the direction of the bottom current−1.00–0.97–0.47
 Direction_curr_NNorthness component of the direction of the bottom current−1.00–1.000.24
 Mud content [%]The silt-clay fraction of the sediment8.4–77.255.6
 Sed_kurtKurtosis of the sediment particle distribution0.70–1.310.94
 Sed_med [μm]Median sediment particle-size28.0–944.9 97.1
 Sed_skewSkewness of the sediment particle-size distribution−0.78–0.360.09
 Sed_sortSorting of the sediment particle-size distribution4.2–8.65.8
Other variables
 BMD [cm]Biogenic mixing depth derived from sediment profile images0.0–3.11.7
 Penetration_depth [cm]Penetration depth of SPI—a proxy for sediment shear strength0.0–19.412.2
 CNCarbon-to-nitrogen-ratio (weight/weight)7.9–19.311.5
 TOC [%]Percentage of total organic carbon1.7–5.63.8
 TN [%]Percentage of organic nitrogen0.12–0.630.34
 VEAS_dist [m]Distance from effluent of VEAS wastewater treatment plant1143–21201640

Three sediment samples from the upper 2 cm of the sediment were collected from each site using a Gemini corer. Sites with sediments too coarse for the Gemini corer were sampled using the van Veen grab. We determined the sediment particle-size fractions by wet-sieving through successively smaller sieves from 2048 μm to 63 μm. These fractions were subsequently subjected to sediment particle-size analysis using GRADISTAT [25]. Sediment profile images [26] were taken using a sediment profile imaging (SPI) device. These images were used to record penetration depth [27], which is a proxy for sediment compactness [28] or sediment shear strength (high penetration depth corresponds to low shear strength), and biogenic mixing depth [29]. The Euclidean distance from the effluent of a wastewater treatment plant (VEAS) was calculated using the GRASS command r.grow.distance. Total organic carbon (TOC) and total nitrogen (TN) were obtained from the sediment samples using a Thermo Finnigan EA1112 Series Flash Elemental Analyzer [30], and the C/N-ratio (weight/weight) was calculated from these two measures.

2.6. Data Analysis

The statistical computing software R, version 2.15 [31], was used for all statistical analyses. Except for the isoMDS function in the MASS package [32], the vegan package version 2.0-4 [33] was used for all multivariate analyses. Prior to analysis, all environmental variables were zero-skewness transformed and ranged onto a 0-1 scale in order to improve homogeneity of variances and equalise the weights attributed to them [34].

The species-by-site matrix was subjected to multivariate analysis according to the principles of the multiple parallel ordination (MPO) framework [19]. The term “multiple” refers to our application of the ordination methods to several species-by-site matrices, each with different weights given to abundance. Abundance weighting was performed by varying the range of the abundance scale, that is, the ratio of the maximum abundance value recorded for any species, and the minimum value, which was always set to 1. A sequence of values from qualitative (presence/absence; ) data to raw species-abundance data (R = the maximum number of individuals of one species recorded in any site) was used, including the full log2-series of values between and (, 2, 4, 8,…, , …, max) as obtained by weighting each abundance value using the power function [35]. The term “parallel” refers to the parallel use of two ordination methods, here detrended correspondence analysis (DCA) [36, 37] and global nonmetric multidimensional scaling (GNMDS) [38] (other combinations of methods can also be used when appropriate (e.g., [39])). The two methods were applied to data matrices for all values. Parallel ordination offers an opportunity for controlling if the two different methods reveal the same structure, which is regarded as an indication that the true underlying structure has been found [18]. Congruence between parallel DCA and GNMDS ordinations (i.e., ordinations based on a specific range of the abundance scale) was assessed by means of Procrustes analysis [38, 40], and by calculating pairwise Kendall’s rank correlation coefficients between corresponding DCA and GNMDS axes. Because the patterns expressed by ordinations obtained with different weighting functions often differed considerably and express different aspects of compositional variation in the (raw) data [19], we report ordination results for three abundance weights of the original species-by-site matrix: no weight given to abundance (, presence-absence), intermediate weight given to abundance (intermediate , see selection criteria below), and strong weight given to abundance (, the raw, or close-to-raw observed abundance). The choice of intermediate weighting function to be used in reporting of results was based on the following criteria: (1) the Procrustes correlation coefficient for the parallel DCA and NMDS ordinations was among the highest in the range of intermediate weightings (see results in Table 3); (2) Kendall’s between first axes of parallel DCA and GNMDS ordinations were among the highest (Table 3), and (3) also the second ordination axes confirmed each other (i.e., tau between second axes of parallel DCA and GNMDS ordinations ≥ 0.4; [41]). The strong weighting function was selected by the same criteria, but considering the strongest weight given to abundance (raw) first and then subsequently less strong weights until the criteria were met.

Range Stress


DCA ordinations were obtained by use of the decorana function in the vegan package with default arguments: detrending by segments, nonlinear rescaling of axes, and no down-weighting of rare species. One hundred two-dimensional GNMDS ordinations with a Bray-Curtis dissimilarity matrix as input were obtained by first applying the isoMDS function with the following arguments: distance = initMDS (dissimilarity matrix), max number of iterations = 500 and convergence tolerance = , number of axes = 2. Subsequently, by application of the post-MDS function, all configurations were centered and rotated to principal components and all axes were individually rescaled linearly into half-change units.

Environmental variables were passively fit to the first two axes of ordinations using the envfit function. This function finds, by bivariate linear regression, vectors pointing in the direction of the most rapid increase in the ordination diagram for each environmental variable [33]. The envfit function uses a randomised permutation procedure to test the hypothesis that the goodness-of-fit of a given variable is not better than that of a random variable (each environmental variable was permuted 999 times). In order for graphs to remain readable, only environmental variables with Kendall’s with ordination axes above a threshold value of 0.225 (which corresponds to a value of approximately 0.10 for this data set with ) were displayed. The ordisurf function was used to display the variation in DCA ordination space of selected variables as well as the total abundance, species richness, and species evenness. The contours were fitted using isotropic smooth surfaces via generalized additive models (GAM), assuming Gaussian (for environmental variables and species evenness) or Poisson (for abundance and species richness) distribution of errors.

Because confirmed axes extract very similar gradient structures, only results by one of the methods, DCA, are shown. DCA was chosen because none of the DCA ordinations showed any sign of distortions such as the tongue effect [38, 42], because DCA provides an intuitively interpretable scaling in SD units [13], and because DCA species scores obtained by use of “Hill’s scaling” can be interpreted as estimates for species optima with respect to the same scaling of the axes as the site scores [36].

3. Results

3.1. Species Richness, Abundance Patterns, and Core Species

A total of 122 species with a total abundance of 12 672 individuals (polychaetes, molluscs, and echinoderms) were observed. Polychaetes were the most species-rich group with 76 taxa, followed by molluscs with 33 taxa and echinoderms with 13 taxa. Polychaetes were also numerically dominant (68.9% of individuals), followed by molluscs (29.2%) and echinoderms (1.9%). Spiophanes kroyeri (2389 individuals), Pseudopolydora paucibranchiata (746), Heteromastus filiformis (585), and Chaetozone setosa (534) were the most abundant Polychaete species; Ennucula tenuis (2216) and Thyasira equalis (465) were the most abundant species of molluscs, while Ophiocten affinis (105) was the most abundant echinoderm species. Sixteen species had a frequency in samples above 90% (presence at 26 or more sites), while a total of 39 species were rare, having a total abundance ≤3. Forty-two species (36.9%) were found in all six, while 13 species (10.7%) occurred in only one or two broad-scaled geomorphological features. Most species of the latter group were related to crest and shallow slope features. These two feature types were also particularly rich in rare species (for a complete overview of fauna and their relation to geomorphological features, see Table 4).

FrequencyAverage abundance in presence sitesfreqoccab.avrab.tot

 Pherusa falcata 1.001.670.1131.675
 Eumida sanguinea 0.670.172.501.000.1132.006
 Eunice pennata 0.670.172.502.000.1132.337
 Streblosoma bairdi 0.330.500.800.400.501.002.331.502.001.500.43121.7521
 Ophelina acuminata 0.330.330.200.400.
 Maldane sarsi 0.330.170.600.400.
 Anobothrus gracilis
 Glycera alba
 Goniada maculata
 Scoletoma fragilis
 Spiophanes kroyeri
 Typosyllis cornuta
 Pholoe inornata
 Polycirrus norvegicus
 Chaetozone setosa
 Pista cristata
 Polyphysia crassa
 Heteromastus filiformis
 Terebellides stroemi
 Pseudopolydora paucibranchiata
 Brada villosa
 Prionospio cirrifera 1.000.830.601.001.001.0060.677.0021.6714.6010.508.200.892517.52438
 Trichobranchus roseus 1.000.830.600.600.750.801.673.403.673.331.673.500.75212.9562
 Prionospio fallax 1.000.671.000.801.000.8025.3316.758.004.2511.005.000.862411.00264
 Diplocirrus glaucus 1.000.670.801.001.000.603.673.754.253.6013.504.000.82235.52127
 Polycirrus medusa 1.000.500.400.600.500.602.671.673.501.0012.001.330.57163.1951
 Melinna elisabethae 1.000.500.200.600.500.606.675.671.001.671.003.000.54153.6054
 Spiochaetopterus typicus 0.671.
 Aphrodita aculeata 0.670.670.400.200.500.403.001.751.
 Melinna cristata 0.670.500.400.801.001.002.502.673.004.252.753.800.71203.3066
 Syllidia armata 0.670.330.800.600.250.602.501.501.
 Scalibregma inflatum 0.670.330.400.600.500.401.501.
 Amaena trilobata 0.670.170.400.400.500.802.
 Trichobranchus glacialis 0.330.670.200.401.000.601.
 Pista lornensis 0.330.170.600.400.750.801.001.001.671.501.001.250.50141.2918
 Sige fusigera 0.330.170.200.400.250.401.
 Eupolymnia nebulosa 0.670.330.
 Glycera lapidum 1.000.330.400.250.2011.676.501.501.001.000.3295.8953
 Cossura longocirrata 0.670.
 Sphaerodorum gracilis 0.670.330.202.501.001.000.1851.608
 Pholoe assimilis 0.330.
 Chaetoparia nilssoni 0.330.600.501.
 Scoloplos armiger 0.670.2019.003.000.11313.6741
 Scoletoma cf. magnidentata 0.330.400.400.500.
 Amphicteis gunneri 0.170.400.200.250.402.
 Ophiodromus flexuosus 0.330.
 Phyllodoce rosea 0.330.400.500.401.501.002.001.500.2981.5012
 Phyllodoce groenlandica 0.330.
 Dipolydora cf. coeca 0.1711.000.04111.0011
 Petaloproctus borealis
 Sabellides octocirrata 0.200.400.
 Orbinia sertulata 0.200.400.
 Nephtys ciliata 0.200.500.401.
 Glycera unicornis 0.600.501.001.000.1851.005
 Rhodine loveni
 Hiatella arctica 0.670.334.002.000.1443.0012
 Lepeta caeca 0.670.1720.006.000.11315.3346
 Mya truncata 0.670.171.501.000.1131.334
 Arctica islandica 0.330.500.201.001.671.000.1851.407
 Thyasira equalis
 Astarte montagui
 Ennucula tenuis 1.000.671.
Pseudamussium peslutrae 1.000.500.400.
 Thyasira sarsi 0.670.830.801.
 Abra nitida 0.330.500.400.800.501.005.002.672.
 Thyasira flexuosa 0.330.330.801.
 Philine scabra 0.330.330.
 Corbula gibba 1.000.830.200.407.0020.605.001.500.391112.00132
 Parvicardium minimum 0.670.
 Macoma calcarea 0.670.400.400.500.607.
 Astarte elliptica 0.670.200.600.
 Thyasira gouldi 0.330.200.400.801.
 Nuculana pernula 0.330.
 Buccinum undatum 0.400.400.250.401.
 Strongylocentrotus droebachiensis 0.672.000.0722.004
 Amphiura chiajei 0.670.830.401.503.603.500.3293.1128
 Ophiocten affinis
 Ophiura albida 1.000.500.400.400.500.402.001.331.001.002.502.000.50141.6423
 Amphiura filiformis 0.670.500.600.400.250.401.505.673.
 Ophiura robusta 0.670.330.
 Amphipholis squamata 0.330.670.400.
 Brissopsis lyrifera 0.670.400.
 Echinocardium flavescens 0.330.400.400.250.601.001.501.001.001.670.36101.3013

 Euchone rubrocinta 0.331.000.0411.001
 Exogone verugera 0.331.000.0411.001
 Harmothoe cf. fragilis 0.332.000.0412.002
 Laonice cirrata 0.332.000.0412.002
 Nicolea zostericola 0.331.000.0411.001
 Pherusa plumosa 0.332.000.0412.002
 Pholoe pallida 0.331.000.0411.001
 Praxillura longissima 0.331.000.0411.001
 Rhoine gracilior 0.332.000.0412.002
 Gyptis rosea 0.670.
 Eteone cf. longa 0.670.
 Galathowenia oculata
 Laetmonice filicornis
 Orbinia norvegica
 Typosyllis armillaris
 Parougia eliasoni
 Flabelligera affinis
 Paraamphinome jeffreysii
 Amage auricola
 Mediomastus fragilis
 Neoamphitrite grayi
 Parvicardium pinnulatum 0.331.000.0411.001
 Tellimya ferruginosa 0.670.
 Aporrhais pespelecani 0.330.
 Cuspidaria cuspidata
 Kurtiella bidentata
 Abra alba
 Astarte sulcata
 Antalis entalis
 Diaphana minuta
 Cylichna alba
 Lunatia montagui
 Mytilus edulis
 Echinocardium cordatum
 Ophiura ophiura
 Ophiopholis aculeata
 Ophiothrix fragilis

Freq: total site frequency; occ: total site occurrences; ab.avr: overall average abundance in presence sites; ab.tot: overall total abundance in presence sites.
3.2. Faunal Analysis: Correspondence between DCA and GNMDS Ordinations

The eigenvalues of DCA1 and DCA2 were 0.184 and 0.0981, respectively. The percent stress for the GNMDS ordinations ranged between 20.3 () and 12.9 (; see Table 3). The Procrustes analysis revealed high configurational concordance between parallel DCA and GNMDS ordinations: Procrustes correlation coefficients ranged between 0.78 and 0.93 (see Table 3). All first axes of parallel DCA and GNMDS ordinations were highly correlated, with Kendall’s values ranging between 0.83 and 0.92 (Table 3). In contrast, the corresponding second axes exhibited substantial variation in Kendall’s between abundance ranges (Table 3). For the qualitative () representation of the species-by-site matrix, the second axis was not confirmed according to the “rule-of-thumb” of between parallel ordinations. The second axes of ordinations with intermediate and strong weights given to abundance were confirmed. Only confirmed axes were subjected to environmental interpretation.

3.3. Description of DCA Ordination Patterns

The DCA ordinations of the species-by-site matrix (Figure 3) exhibited a dense clustering of sites at the intermediate abundance range (), while, for both the qualitative data set () and for data with strong weight given to abundance (), this cluster was less tight. Crest sites always obtained high scores on DCA1, although a shift of crest site 6 towards the main cluster with increasing was observed. Gradient lengths of DCA1 increased with increasing weight given to abundance (1.97 SD units for to 3.00 SD units for ). This difference in gradient length was mainly due to variation in the degree of separation of crest and shallow slope sites from the main cluster at the high-score end along DCA1. The shallow slope sites exhibited relatively high scores on both axes, irrespective of abundance range. The remaining broad-scaled features remained clustered at low scores although there was a tendency for shallow flatsites to move towards higher scores with increasing abundance range.

3.4. Relationships between DCA Axes and Recorded Environmental Variables

There was a general trend for environmental variables to be more strongly correlated with DCA axes obtained for strong weighting of abundance (Figure 4; see Table 5 for a complete overview of the relationship between DCA axes and environmental variables). Some environmental variables were, however, most strongly correlated with the DCA axes obtained for qualitative data (), for example, BPI (DCA1) and sediment sorting (DCA2), while other variables were most strongly correlated with DCA axes obtained with intermediate weighting of abundance (), for example, depth and sediment kurtosis (DCA1), and the northing component of the direction of the bottom current and TOC (DCA2). The overall highest correlation coefficients were obtained for depth, which was highly, negatively, correlated with DCA1 regardless of abundance range (0.66 ≤ ≤ 0.82). Biogenic mixing depth and maximum surface current speed were also highly, negatively, correlated with DCA1 axes ( ≥ 0.3). Biogenic mixing depth, sediment skewness, and maximum curvature of the seabed all had correlation coefficients above 0.3 also with DCA2.

(qual) (raw) (qual) (raw)

DEM derived
 Aspect_E0.011−0.026−0.0900.148 0.0210.196
 Aspect_N0.1010.032−0.074−0.069 −0.143−0.085
 BPI0.2280.1690.148−0.016 −0.164−0.201
 Depth−0.656−0.820−0.757−0.053 −0.0630.079
 Seabed_curv0.1690.1320.079−0.127 −0.169−0.312
 Slope0.127−0.005−0.153−0.148 −0.021−0.132
Current model
 Current_bottom−0.026−0.095−0.085−0.037 −0.1430.085
 Current_surface−0.307−0.407−0.471−0.212 −0.106−0.090
 Change_curr_bottom0.0740.1010.132−0.032 −0.053−0.005
 Change_curr_surface−0.058−0.138−0.1160.090 0.1430.254
 Direction_curr_E0.024−0.024−0.083−0.008 −0.110−0.143
 Direction_curr_N0.1790.0940.083−0.035 −0.265−0.099
Sediment features
 Mud content0.0530.1220.069−0.116 0.0950.164
 Sed_kurt0.2400.2720.245−0.191 −0.046−0.073
 Sed_med−0.027−0.074−0.0950.090 0.000−0.032
 Sed_skew0.0990.1840.152−0.056 0.0880.318
Other variables
 Penetration_depth−0.174−0.227−0.2860.184 0.0670.206
 CN0.2650.2860.318−0.170 −0.202−0.260
 TN−0.257−0.349−0.393−0.149 −0.1490.008
 VEAS_dist0.3490.4070.4070.138 −0.021−0.122

Bold values represent ≥ 0.225, which corresponds to a P value of approximately 0.10 in test of the hypothesis that when . The unconfirmed second axis (DCA2. ) according to criteria given in Section 2.6 is shown in italic.

The DCA ordination with intermediate weighting of abundance () was selected for particularly more careful interpretation, which showed that sediment particle-size variables were in general weakly correlated with the main coenocline. However, contour-plots for individual environmental variables (Figure 5) revealed that particle-size variables such as median sediment particle-size and mud content were both unimodally related to DCA1, with extreme values (minimum and maximum, resp.) just to the right of DCA1 = 0. This motivated separate correlation analyses for the positive, Subset+, and negative, Subset−, ends of DCA1. We allowed for a slight overlap between the subsets, and included in Subset+ all sites with DCA1 > 0 and in Subset− all sites with DCA1 ≤ 0.1. In Subset+, depth, median sediment particle-size, mud content, and the CN-ratio were highly correlated with DCA1 (Table 6) and made up a group of intercorrelated variables (Table 7), while in Subset− depth, which was not strongly correlated with any other variable, was the only variable highly correlated with DCA1 (Tables 6 and 7; see Table 8 for pair-wise correlations between all environmental variables). Sediment profile imagery (SPI) revealed that there was substantial difference in sediment particle-composition between extreme ends of DCA1 (Figure 6).

τ τ

DEM derived
Current model
Sediment features
 Mud content–0.640.0310.280.071
Other variables

DepthBMDPenetration depthCNSed_medMud contentTNTOC

Depth *****0.430.36−0.57−0.710.710.470.36
BMD0.28 *****0.93−0.29−
Penetration_depth−0.160.38 *****−0.21−0.360.210.040.00
CN−0.17−0.55−0.12 *****0.43−0.57−0.33−0.21
Sed_med0.140.380.19−0.26 *****−0.86−0.26−0.21
Mud content−0.26−0.30−0.010.11−0.68 *****0.330.21
TN0.310.28−0.12−0.380.35−0.44 *****0.91
TOC0.230.04−−0.410.56 *****

Mud contentMed_

BPI ***** 0.6460.079 0.185 0.2010.164 0.1900.163 0.077−0.223 0.106
Depth0.101 ***** 0.0740.051 0.481 0.032 0.4070.1310.3400.387−0.271−0.439
<0.0010.201 *****0.0530.0780.319 0.0900.223 0.1430.019 0.014 0.101
0.5700.5970.710 ***** 0.233 0.286 0.159 0.032 0.1600.0890.0050.254
0.2270.7070.5660.593 *****0.3190.089 0.0770.003 −0.234−0.239 0.0270.3490.1870.2560.024
0.1720.6210.0180.0850.020 ***** 0.2970.316−0.2440.1740.095−0.252 0.0220.0890.0600.0780.142
0.356<0.0010.2630.1390.5130.373 ***** 0.011 0.2060.0450.3830.338 −0.524
0.1390.5180.5180.0340.3940.0280.518 *****0.202−0.2650.1960.158−0.232 0.1110.196
Sed_kurt0.2270.1170.1000.4160.5780.0210.4400.137 *****−0.5670.5680.457−0.256−0.342 0.0670.073
Sed_med0.3230.8130.2520.2360.9840.0720.9370.048<0.001 *****−0.796−0.5220.4500.114 0.1760.209
Mudcontent0.1620.4210.2980.3170.3730.1980.5700.150<0.001<0.001 *****0.724−0.249 0.136−0.277−0.257 0.005
Sed_skew0.2270.1720.8900.4400.0880.4880.3130.2430.001<0.001<0.001 *****0.041 0.248−0.352−0.268
Sed_sort0.8120.2670.1220.8120.0840.0670.6060.0880.0640.0010.0680.766 *****0.1220.136
BMD0.3020.0030.1040.5000.7650.3500.1320.2040.0140.4040.5780.9210.381 *****0.5340.0900.258−0.514−0.385
0.5660.3320.7070.8280.8430.8740.7370.7370.2660.7670.3130.0680.321<0.001 ***** −0.228
TOC0.0970.0110.6070.2350.0110.5130.0040.4520.7510.1920.0400.0090.3610.5120.921 *****0.6670.000−0.277
TON0.4630.0050.9210.5130.1760.6620.0130.2750.6620.1260.0590.0520.7960.0640.937<0.001 *****−0.361−0.306
CN0.3040.0440.7220.9680.0600.5660.5530.4060.6200.9680.4770.2430.405<0.0010.0921.0000.008 *****0.207
VEAS_dist0.4450.0010.4680.0600.8580.294<0.0010.1500.5930.7220.9840.7670.7810.0050.2590.0400.0250.123 *****

The negative end of the main coenocline (DCA1) was made up by deep flats sites (7 and 16) and deep slope sites (22 and 23). Subsurface deposit feeders such as bivalves belonging to the Thyasira genus, the bivalve Ennucula tenuis, and the polychaete Spiophanes kroyeri were associated with this end (Table 9). Surface deposit feeders such as the polychaetes Melinna cristata and Pista lornensis also had abundance peaks near this end. The positive end of the main coenocline (crest sites 5 and 25 and the shallow-slope sites 9, 20, and 28) was dominated by suspension and surface deposit feeders. Suspension feeders are exemplified by the molluscs Corbula gibba, Pseudamussium peslutrae, and Hiatella arctica and surface deposit feeders by the polychaetes Pseudopolydora paucibranchiata, Anobothrus gracilis, and Melinna elisabethae. Crest sites 5 and 6 made up the negative end of the secondary coenocline (DCA2). No feeding habit was dominating, and a mixture of suspension feeders (e.g., the mollusc Astarte elliptica), predators (e.g., the polychaete Aphrodita aculeata), and subsurface deposit feeders (e.g., Macoma calcarea) was recorded in these sites. At the opposite, positive end of the secondary coenocline, surface deposit feeders such as the polychaetes Dipolydora coeca, Amphicteis gunneri, and Streblosoma bairdi (all of which are also facultative suspension feeders) were recorded as important species.


DCA1Neg−2.00Thyasira gouldi (M)xx17
−1.96Sige fusigera (P)xx13
−1.81Pista lornensis (P)x18
−1.47Thyasira flexuosa (M)xx206
−1.46Melinna cristata (P)x66
−1.36Thyasira sarsi (M)xx207
−1.34Ennucula tenuis (M)x2216
−1.26Spiophanes kroyeri (P)xx2389
Pos2.78Lepeta caeca (M)xx46
2.59Hiatella arctica (M)xx12
2.46Scoloplos armiger (P)x41
1.94Pseudopolydora pauchibranchiata (P)xx746
1.76Corbula gibba (M)x132
1.75Glycera lapidum (P)xx53
1.66Eupolymnia nebulosa (P)x15
1.61Amphiura chiajei (E)xx28
1.58Pseudamussium peslutrae (M)x19
1.44Ophiocten affinis (E)x105
1.41Melinna elisabethae (P)xx54
1.31Anobothrus gracilis (P)x381

DCA2Neg−3.38Scoloplos armiger (P)x41
−1.79Astarte elliptica (M)x19
−1.71Macoma calcarea (M)xxx70
−1.65Aphrodita aculeata (P) x24
−1.56Ophiura albida (E)xxx23
−1.35Lepeta caeca (M)xx46
Pos4.81Dipolydora coeca (P)xx11
1.77Amphicteis gunneri (P)xx11
1.65Amphiura chiajei (E)xx28
1.59Phyllodoce rosea (P)x12
1.54Streblosoma bairdi (P)xx21
1.38Corbula gibba (M)x132

Feeding strategies: Su: suspension; SD: surface deposit; SSD: subsurface deposit; G: grazer; Sc: scavenger; P: predator. Taxon: E: echinoderm; M: mollusca; P: polychaeta. Bold symbols indicate the feeding strategy that species with multiple feeding strategies have the highest affinity for. Data on feeding strategy collected from the Norwegian Institute of Water Research's Biological Traits Database [75] and World Register of Marine Species [76].

The total number of individuals (total abundance) was correlated with species richness ( = 0.37, ). Furthermore, the total number of individuals increased from the midpoint of the main coenocline towards its positive end and also, more weakly, towards its negative end. Total abundance also decreased from the negative to positive end of the secondary coenocline (Figure 7(a)). Species richness increased from the upper left (negative DCA1, positive DCA2) to the lower right (positive DCA1, negative DCA2) in the ordination diagram (Figure 7(b)). Species evenness was unimodally related to the main coenocline (Figure 7(c)).

4. Discussion

We found high similarity between the first axes obtained for the species-by-site matrix by the two ordination methods (DCA and GNMDS). This clearly indicates the existence of a main gradient in the composition of polychaetes, molluscs, and echinoderms, that is, a main coenocline, in the study area. The similarity between second ordination axes of our parallel ordinations was not as strong and consistent, indicating a considerably weaker secondary coenocline.

4.1. Ecoclines Related to the Main Coenocline (DCA1)

Based upon the contrasting patterns of correlations between ordination axes and environmental variables, and between the environmental variables themselves, observed for the two ends of the main coenocline (Subset+ and Subset−), we hypothesise that the main coenocline is related to two ecoclines. This hypothesis is also underpinned by the change in the nature of species-environment relationships observed with change in the weight attributed to abundance.

Near the positive end of the main coenocline (Subset+), depth forms a major complex-gradient together with mud content, median sediment particle-size, and the CN-ratio. This complex-gradient is the most important for explaining variation in species composition in Subset+. We hypothesise that periodically occurring, extreme cases of strong currents is the main structuring process that brings about this pattern. Accordingly, we have named the term periodic physical forcing for the ecocline formed by the covariation of this major complex-gradient and the corresponding coenocline.

Support for this hypothesis comes from the extreme sites in Subset+ being crest and slope sites found at shallow depth with low mud content, high median particle-size, and low organic quality (i.e., high CN-ratio). Both sediment profile images and surface images confirm that these sites are dominated by coarse sediments (including cobbles and stones in the most extreme sites). In the images (Figure 6) a thin layer of flocculated sediment was observed, indicating that some sedimentation or bed load transport from areas with softer sediments may take place between intermittent, extreme cases of physical forcing.

Such extreme cases of physical forcing is likely to take place when the inner basin of the Oslofjord is completely renewed by more saline, dense, and colder water from the deep water outside the fjord’s sills. Such deep-water renewals are likely to occur every winter [43]. The degree of deep-water renewal is dependent on the physical forcing and extent of northerly winds that push surface water out of the fjord [43] that allows deeper water to pass the sill and increase bottom current velocity in the inner fjord.

The presence of the surface deposit-feeding polychaetes A. gracilis, M. elisabethae, and P. paucibranchiata (see Table 9) also gives support to the physical forcing hypothesis: their presence suggests good supply of organic material to these sites, while the particle-size distribution suggests that this material is removed during periods of strong physical forcing.

Sediment particle-size variables have traditionally been considered the most important factors explaining the distribution of macrofauna in marine sediments [1, 2]. In our study we do, however, observe a nonlinear relationship between particle-size variables and the main coenocline. Although particle-size variables, such as mud content and median particle-size, are important for explaining variation near the positive end of the main coenocline, our results comply with the hypothesis that sediment particle-size composition is not the crucial factor in explaining gradients in species composition per se [4] but rather a result of other processes that affect the sediment properties at a particular site. It should be noted that the method used for estimating particle-size distributions by wet or dry sieving disrupts the structure of coagulated and flocculated sediment particles and may therefore fail to reflect the conditions experienced by organisms in the wild [44]. Sediment particle-size properties, recorded to describe conditions experienced by the macrofauna, should therefore preferably be recorded in situ [45, 46].

Depth, which was not strongly correlated with any other variable in Subset–, was the only variable strongly correlated with the main coenocline in that subset. Given the restricted variation in depth encountered in the study area, depth as such is hardly directly important for the species composition. Instead, we interpret depth as a proxy for many variables that are difficult to measure [47, 48]. Although not directly correlated with depth, other variables such as concentrations of total carbon and nitrogen as well as biogenic mixing depth also peaked near the negative end of this coenocline (see Table 5). These auxiliary variables do not form a complex-gradient but may separately explain some of the variation in species composition along this coenocline.

We hypothesise that duration of periodic hypoxia, for which depth acts as a proxy, is the main structuring process that brings about the observed pattern. Accordingly, the term periodic hypoxia is used for this ecocline. The periodic hypoxia hypothesis is supported by total carbon concentrations above 4% being recorded in the deepest-situated sites near the low-score end of the main coenocline. The total carbon concentrations recorded at these sites are thus above threshold concentrations (approximately 3.5%) for possible negative impacts of organic input on species richness as reported by Hyland et al. [5]. Accordingly, the observed reduction in species richness and abundance at these sites may indicate organic matter enrichment, which in turn brings about hypoxia in periods of water stagnation. The presence of T. flexuosa and T. sarsi, regarded as typical of organically enriched environments [49], further supports this hypothesis.

Additional support for this hypothesis comes from measurements of oxygen concentrations at sites close to the study area, showing that the deepest sites are likely to experience periods of hypoxia during periods of stagnant water masses [5053]. It is worth noticing that sampling of species and recording of environmental variables for the present was performed in early spring (March and April), after the spring circulation of water masses, in a period of normoxic conditions [52, 53]. This may explain the relatively low CN-ratios, between 7.9 and 10.5, at these sites, which indicates a relatively recent influx of fresh planktonic matter.

The observed reduction in species richness and abundance near the negative end of the main coenocline also accords with the hypothesis that the fauna is negatively impacted by recurrent periods with hypoxic conditions [54]. This opens for the possibility that the time between recurrent periods of hypoxia is too short to allow for complete faunal recovery before the next hypoxic event takes place [48, 55]. However, observations of high rates of biogenic mixing combined with high organic concentrations at these sites indicate that at the time of sampling, these sites were not exposed to hypoxia. Therefore, both biogenic mixing and concentrations of organic content appear to have important structuring effects on associated fauna between episodes of hypoxia.

4.2. The Ecocline Related to the Secondary Coenocline (DCA2)

The secondary coenocline was less consistently related to environmental variables and the results partly failed to provide direct indications of an important complex-gradient likely to be responsible for variation along this coenocline. It should, however, be noted that some of the variation along this ecocline is also accounted for by the main ecoclines. This is evident from the vectors for single environmental variables that contribute to the main ecocline, notably biogenic mixing depth and total organic carbon (Figure 4), which are represented by vectors that are correlated with both the main and the secondary coenoclines. However, the results also indicate existence of one additional complex-gradient, which accounts for residual compositional variation mainly expressed along DCA2. With reference to the group of variables most strongly correlated with the secondary coenocline, which consists of acceleration/deceleration of the bottom current, maximum bottom current, and the northing component of the bottom current, we hypothesise that modelled bottom-current strength is the factor contributing the most to this ecocline. Running from low maximum bottom-current speed at high DCA2 scores to high maximum bottom-current speed at low DCA2 scores, we term this ecocline the bottom current gradient.

A weak “ecocline of residual variation” related to current strength was not unexpected because the hydrodynamic regime of the study area is determined by tidal currents and variation in wave exposure, both of which have major impacts on sediment stability [56] and influx of organic material [4]. This makes the hydrodynamic regime an important determinant of the sedimentary characteristics of an area [57]. What is surprising though is that the species composition in the study area indicates variation related to a bottom current ecocline that is more or less unrelated to variation related to the periodic physical forcing ecocline.

The shallow-depth crest sites 5 and 25 and the shallow-slope site 9, all situated on coarse sediments, are evidently exposed to physical forcing beyond the erosion threshold. The physical forcing taking place at these sites is most probably caused by tidal currents, given that the shallowest sites are situated at depths below 30 m where wave exposure is likely to be minimal [57, 58]. Our current model predicts bottom currents of 5-6 cm/s at these sites which, according to the Hjulström diagram [59], are above the erosion threshold for coarse sand in noncohesive sediments [56]. The observed general lack of correlation between current variables and the main coenocline in our data may, however, result from a mismatch between the spatial resolution of the current model and the spatial resolution at which species respond to the structuring processes, in this case erosion caused by physical forcing.

4.3. Mechanisms and Processes That Account for Variation along Ecoclines

Our results indicate differentiation along the main coenocline between communities dominated by suspension and surface deposit-feeders in shallow waters where physical forcing is likely to be the main structuring factor, and communities dominated by subsurface deposit-feeders in deeper waters where periodic hypoxia and high organic input are likely to be important. Predators are evenly distributed in the area, seemingly not strongly responsive to any of the two ecoclines.

Near the physical forcing end of the main coenocline, a suspension-feeding community is observed, as expected, but also other feeding strategies are observed there. Low compositional turnover (relatively short coenoclines in terms of SD units used for scaling of DCA axes) may explain the lack of a clear-cut difference in trophic community structure. Gray [60] suggests that mixed communities of deposit-feeders and suspension-feeders is the general rule in marine benthic communities. Another factor that may contribute to making coenoclines short is the high degree of patchiness in marine benthic communities [61, 62]. Such patchiness results from the fine spatial and temporal scales on which the disturbances that affect the benthos operate. For example, Kelaher and Levinton [63] demonstrates that spatiotemporal variation in organic enrichment generates complex patchiness both in species assemblages and trophic structure, and Wieking and Kroncke [64] find the highest trophic diversity (number of feeding types present) at sites characterized by relatively low but highly variable input of intermediate-quality organic matter.

Our results show coincident patterns of total abundance and species richness, both of which reach the lowest values near the periodic hypoxia end and peak near the physical forcing end of the main coenocline. Low total abundance and low species richness at sites with periodic hypoxia and elevated concentrations of organic matter accord with the hypothesis of Pearson and Rosenberg [3] that influx, and increase, of organic content will lead to an increase in abundance until the point is reached where oxygen comes in short supply, beyond which abundance rapidly declines. This explanation is also supported by the results of Ramey and Snelgrove [65] who reported a range of TOC concentrations (0.9–8.1%) similar to ours (1.7–5.6%) and found a negative correlation between total abundance and organic content, and by Grebmeier et al. [9] and  Ambrose Jr. and Renaud [66] who found positive correlations between total abundance and organic matter content in data sets with substantially lower TOC concentrations than ours (0.1–1.9% and 0.5–1.5%, resp.). Furthermore, Buhl-Mortensen et al. [54], in a comprehensive study of 11 fjords, found a clear negative impact of hypoxia on species richness.

The high total abundance and high species richness near the physical forcing end of the coenocline may be due to higher niche diversity, brought about by increasing structural complexity with increasing sediment coarseness, as predicted by Gray [1]. Increasing species richness and more diverse communities with increasing sediment particle-size is also reported by, for example, Grebmeier et al. [9], Buhl-Mortensen and Høisæter [67], Newell et al. [68], and Trannum et al. [69]. Higher species richness at structurally more complex sites because of higher niche diversity also accords with general ecological theory [70]. The positive correlation between total abundance and species richness indicates that competition, that is, mutually negative interspecific interactions [71], is less important than environmental stress in the studied ecosystem.

4.4. Indirect versus Direct Gradient Analysis and Weight Given to Abundance

The methodological approach used in this study is to apply indirect gradient analysis (i.e., unconstrained ordination methods such as DCA and NMDS) in a general-purpose, inductive study [18] with the main aim to generate hypotheses about species-environment relationships. Our results demonstrate that the chosen method gives us new insight into the structuring of species composition by environmental complex-gradients. An alternative analytic approach could have been to use direct gradient analysis, based on constrained ordination methods such as Canonical Correspondence Analysis (CCA) and Redundancy Analysis (RDA) or based upon the BIO-ENV procedure in PRIMER, which also belongs to direct gradient analytic methods because the response of the community as a whole (represented by the dissimilarity matrix, not an ordination) is correlated with subsets of recorded environmental variables.

In unconstrained ordination, variation in species composition extracted on the axes is not constrained by the recorded environmental variables. In contrast, such constraining is an essential feature of constrained methods. While in unconstrained ordination the story about gradient relationships is told exclusively from the perspective of the species, constrained ordination summarises the responses of species to the recorded environmental variables. Like other models for relating response variables to independent variables, constrained ordination thus only summarises variation in species composition that is explained by the supplied variables. This has several important implications [18]: (1) compositional variation due to unmeasured, but important, variables is lost in a black box of “unexplained variation”; (2) compositional variation not linearly related to the supplied variables is lost as residual variation; and (3) being (linear) combinations of environmental gradients, axes of direct gradient analysis are complex-gradients, not gradients in species composition. Økland [18] therefore argues that indirect gradient analysis should be used to generate hypotheses about species-environment relationships, while direct gradient analysis should be used for testing hypotheses or for partitioning the variation in species abundances in a set of sites or on sets of explanatory variables [18].

Relationships between the main coenoclines, complex-gradients, and auxiliary variables were consistent across weighting functions, however; a clear trend is observed for species-environment relationships to become stronger with increasing weight given to abundance. This result accords with Cushman and McGarigal [72] who compared species-environment relationships of birds and Rueda and Salas [73] who studied molluscs. They found that, overall, the gradient structure of their data set was more strongly related to the environment when raw abundance data were used. However, examples also exist of studies (e.g., [74]) in which species-environment relationships are better explained by presence-absence data, emphasizing presence of rare and low-abundant species. In general, the probability that a data set contains distinctive quantitative abundance information is likely to decrease with increasing compositional and environmental variation in the data [13].

5. Conclusion

We provide strong indications that the main gradient in species composition (i.e., the main coenocline) in our temperate fjord system is related to periodic events of hypoxia and strong physical forcing (Figure 8). Furthermore, we find complex relationships between the main coenocline and several other environmental variables, including the influx of organic matter, the mud content, and biogenic mixing depth. The results are consistent with a hypothesis that different ecological processes contribute to explaining variation in species composition along the main coenocline and that these processes vary in their importance along the coenocline. We conclude that the ecological processes that have the strongest effect on the species composition in our study area are likely to include both periodically recurrent and rare, extreme events. A minor coenocline related to modelled current strength is also identified.

Based on our results we argue that identification of gradients in species composition by indirect gradient analysis, followed by interpretation of these coenoclines in terms of important complex-gradients and inference of ecoclinal relationships, is a useful strategy for basic, general-purpose ecological studies in which focus is on generating rather than testing of hypotheses. We also argue that application of indirect gradient analysis methods on species-by-site matrices from a wide variety of ecosystems and spatial scales (grains and extents) will enhance our basic mechanistic understanding of the ecology of marine sediment communities.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


The authors would like to thank Rita Amundsen and Berit Kaasa for assistance in the lab; Øystein Stokland for assistance in identifying mollusc species; Andre Staalstrøm for providing the ocean current model; the late Frode Olsgard, Bjørnar Beylich, Hasse Nilsson, Ketil Hylland, Gunhild Borgersen, Master students, and the crew of R/V Trygve Braarud for assisting in the field work; Eivind Oug for giving very insightful comments to the paper; and Norwegian Institute of Water Research (NIVA) for providing access to the SPI device and the Gemini corer. An earlier version of this paper benefitted substantially from comments by three anonymous reviewers.


  1. J. S. Gray, “Animal-sediment relationships,” Oceanography and Marine Biology, vol. 12, pp. 223–261, 1974. View at: Google Scholar
  2. D. P. Weston, “Macrobenthos-sediment relationships on the continental shelf off cape hatteras, North Carolina,” Continental Shelf Research, vol. 8, no. 3, pp. 267–286, 1988. View at: Google Scholar
  3. T. H. Pearson and R. Rosenberg, “Macrobenthic succession in relation to organic enrichment and pollution of the marine environment,” Oceanography and Marine Biology, vol. 16, pp. 229–311, 1978. View at: Google Scholar
  4. R. Rosenberg, “Benthic marine fauna structured by hydrodynamic processes and food availability,” Netherlands Journal of Sea Research, vol. 34, no. 4, pp. 303–317, 1995. View at: Publisher Site | Google Scholar
  5. J. Hyland, L. Balthis, I. Karakassis et al., “Organic carbon content of sediments as an indicator of stress in the marine benthos,” Marine Ecology Progress Series, vol. 295, pp. 91–103, 2005. View at: Google Scholar
  6. R. N. Hughes and M. L. H. Thomas, “Classification and ordination of benthic samples from Bedeque Bay, an estuary in Prince Edward Island, Canada,” Marine Biology, vol. 10, pp. 227–235, 1971. View at: Google Scholar
  7. A. J. Hirst, “Broad-scale environmental gradients among estuarine benthic macrofaunal assemblages of south-eastern Australia implications for monitoring estuaries,” Marine and Freshwater Research, vol. 55, no. 1, pp. 79–92, 2004. View at: Publisher Site | Google Scholar
  8. J. B. Buchanan, M. Sheader, and P. F. Kingston, “Sources of Variability in the Benthic Macrofauna off the South Northumberland Coast, 1971–1976,” Journal of the Marine Biological Association of the United Kingdom, vol. 58, pp. 191–209, 1978. View at: Google Scholar
  9. J. M. Grebmeier, H. M. Feder, and C. P. McRoy, “Pelagic-benthic coupling on the shelf of the northern Bering and Chukchi Seas. 2. Benthic community structure,” Marine Ecology Progress Series, vol. 51, pp. 253–268, 1989. View at: Google Scholar
  10. F. Olsgard, “Do toxic algal blooms affect subtidal soft-bottom communities?” Marine Ecology Progress Series, vol. 102, no. 3, pp. 269–286, 1993. View at: Google Scholar
  11. A. Mclachlan, A. C. Cockcroft, and D. E. Malan, “Benthic faunal response to a high energy gradient,” Marine Ecology Progress Series, vol. 16, pp. 51–63, 1984. View at: Google Scholar
  12. P. C. Fleischack and A. J. de Freitas, “Physical parameters influencing the zonation of surf zone benthos,” Estuarine, Coastal and Shelf Science, vol. 28, no. 5, pp. 517–530, 1989. View at: Google Scholar
  13. R. Økland, “Vegetation ecology: theory, methods and applications with reference to Fennoscandia,” Sommerfeltia, vol. 1, pp. 1–233, 1990. View at: Google Scholar
  14. R. H. Whittaker, “Gradient analysis of vegetation,” Biological Reviews of the Cambridge Philosophical Society, vol. 42, no. 2, pp. 207–264, 1967. View at: Google Scholar
  15. R. Halvorsen, “A gradient analytic perspective on distribution modelling,” Sommerfeltia, vol. 35, pp. 1–165, 2012. View at: Google Scholar
  16. R. H. Økland, “Wise use of statistical tools in ecological field studies,” Folia Geobotanica, vol. 42, no. 2, pp. 123–140, 2007. View at: Publisher Site | Google Scholar
  17. K. R. Clarke and R. N. Gorley, PRIMER v6: User Manual/Tutorial. PRIMER-E, Plymouth, 2006.
  18. R. H. Økland, “Are ordination and constrained ordination alternative or complementary strategies in general ecological studies?” Journal of Vegetation Science, vol. 7, no. 2, pp. 289–292, 1996. View at: Google Scholar
  19. T. C. van Son and R. Halvorsen, “Multiple parallel ordination: the importance of choice of ordination method and weighting of species abundance data”,” Sommerfeltia. In press. View at: Google Scholar
  20. A. Lepland, R. Bøe, A. Lepland, and O. Totland, “Monitoring the volume and lateral spread of disposed sediments by acoustic methods, Oslo Harbor, Norway,” Journal of Environmental Management, vol. 90, no. 11, pp. 3589–3598, 2009. View at: Publisher Site | Google Scholar
  21. GRASS Development Team, Geographic Resources Analysis Support System (GRASS) Software, Open Source Geospatial Foundation, 2012.
  22. M. F. J. Wilson, B. O'Connell, C. Brown, J. C. Guinan, and A. J. Grehan, “Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope,” Marine Geodesy, vol. 30, no. 1-2, pp. 3–36, 2007. View at: Publisher Site | Google Scholar
  23. E. R. Lundblad, D. J. Wright, J. Miller et al., “A benthic terrain classification scheme for American Samoa,” Marine Geodesy, vol. 29, no. 2, pp. 89–111, 2006. View at: Publisher Site | Google Scholar
  24. V. Olaya, “Basic land-surface parameters,” Developments in Soil Science, vol. 33, pp. 141–169, 2009. View at: Publisher Site | Google Scholar
  25. S. J. Blott and K. Pye, “Gradistat: a grain size distribution and statistics package for the analysis of unconsolidated sediments,” Earth Surface Processes and Landforms, vol. 26, no. 11, pp. 1237–1248, 2001. View at: Publisher Site | Google Scholar
  26. D. C. Rhoads and S. Cande, “Sediment profile camera for in situ study of organism-sediment relations,” Limnology and Oceanography, vol. 16, pp. 110–114, 1971. View at: Google Scholar
  27. D. C. Rhoads and J. D. Germano, “Characterization of organism-sediment relations using sediment profile imaging: an efficient method of remote ecological monitoring of the seafloor (Remots system),” Marine Ecology Progress Series, vol. 8, pp. 115–128, 1982. View at: Google Scholar
  28. D. H. Wilber, G. L. Ray, D. G. Clarke, and R. J. Diaz, “Responses of benthic infauna to large-scale sediment disturbance in corpus christi bay, Texas,” Journal of Experimental Marine Biology and Ecology, vol. 365, no. 1, pp. 13–22, 2008. View at: Publisher Site | Google Scholar
  29. M. Solan, B. J. Cardinale, A. L. Downing, K. A. M. Engelhardt, J. L. Ruesink, and D. S. Srivastava, “Extinction and ecosystem function in the marine benthos,” Science, vol. 306, no. 5699, pp. 1177–1180, 2004. View at: Publisher Site | Google Scholar
  30. J. I. Hedges and J. H. Stern, “Carbon and nitrogen determinations of carbonate-containing solids,” Limnology & Oceanography, vol. 29, no. 3, pp. 657–663, 1984. View at: Google Scholar
  31. R Development Core Team, R: A Language and Environment for Statistical Computing, Vienna, Austria, R Foundation for Statistical Computing, 2012.
  32. W. N. Venables and R. D. Ripley, Modern Applied Statistics with S, Stanford University, 2002.
  33. J. Oksanen, F. G. Blanchet, R. Kindt et al., “Community Ecology Package, ‘Vegan’,” Version 2. 0-4, 2012. View at: Google Scholar
  34. R. Økland, T. Økland, and K. Rydgren, “Vegetation-environment relationships of boreal spruce swamp forests in Østmarka Nature Reserve, SE Norway,” Sommerfeltia, vol. 29, pp. 1–190, 2001. View at: Google Scholar
  35. E. van der Maabel, “Transformation of cover-abundance values in phytosociology and its effects on community similarity,” Vegetatio, vol. 39, no. 2, pp. 97–114, 1979. View at: Publisher Site | Google Scholar
  36. M. O. Hill, DECORANA—A FORTRAN Program for Detrended Correspondence Analysis and Reciprocal Averaging, Ithaca, NY, USA, 1979.
  37. M. O. Hill and H. G. Gauch Jr., “Detrended correspondence analysis: an improved ordination technique,” Vegetatio, vol. 42, no. 1–3, pp. 47–58, 1980. View at: Publisher Site | Google Scholar
  38. P. R. Minchin, “An evaluation of the relative robustness of techniques for ecological ordination,” Vegetatio, vol. 69, no. 1–3, pp. 89–107, 1987. View at: Publisher Site | Google Scholar
  39. T. C. van Son, E. Oug, R. Halvorsen, and F. Melsom, “Gradients in traits composition and their relation to environmental complex-gradients and structuring processes: a study of marine sediment species communities,” The Open Marine Biology Journal, vol. 7, pp. 14–27, 2013. View at: Google Scholar
  40. P. H. Schönemann and R. M. Carroll, “Fitting one matrix to another under choice of a central dilation and a rigid motion,” Psychometrika, vol. 35, no. 2, pp. 245–255, 1970. View at: Publisher Site | Google Scholar
  41. H. Liu, T. Økland, R. Halvorsen et al., “Gradient analyses of forests ground vegetation and its relationships to environmental variables in five subtropical forest areas, S and SW China,” Sommerfeltia, vol. 32, pp. 1–196, 2008. View at: Google Scholar
  42. R. H. Økland, “A phytoecological study of the mire Northern Kisselbergmosen, SE Norway. II. Identification of gradients by detrended (canonical) correspondence analysis,” Nordic Journal of Botany, vol. 10, pp. 79–108, 1990. View at: Google Scholar
  43. H. G. Gade, “Horizontal and vertical exchanges and diffusion in the water masses of the oslo fjord,” Helgoländer Wissenschaftliche Meeresuntersuchungen, vol. 17, no. 1–4, pp. 462–475, 1968. View at: Publisher Site | Google Scholar
  44. P. V. R. Snelgrove and C. A. Butman, “Animal-sediment relationships revisited: cause versus effect,” Oceanography and Marine Biology, vol. 32, pp. 111–178, 1994. View at: Google Scholar
  45. J. M. Phillips and D. E. Walling, “The particle size characteristics of fine-grained channel deposits in the River Exe Basin, Devon, UK,” Hydrological Processes, vol. 13, pp. 1–19, 1999. View at: Google Scholar
  46. G. Wharton, J. A. Cotton, R. S. Wotton et al., “Macrophytes and suspension-feeding invertebrates modify flows and fine sediments in the Frome and Piddle catchments, Dorset (UK),” Journal of Hydrology, vol. 330, no. 1-2, pp. 171–184, 2006. View at: Publisher Site | Google Scholar
  47. R. Rosenberg, H. C. Nilsson, B. Hellman, and S. Agrenius, “Depth correlated benthic faunal quantity and infaunal burrow structures on the slopes of a marine depression,” Estuarine, Coastal and Shelf Science, vol. 50, no. 6, pp. 843–853, 2000. View at: Publisher Site | Google Scholar
  48. R. Rosenberg, S. Agrenius, B. Hellman, H. C. Nilsson, and K. Norling, “Recovery of marine benthic habitats and fauna in a Swedish fjord following improved oxygen conditions,” Marine Ecology Progress Series, vol. 234, pp. 43–53, 2002. View at: Google Scholar
  49. P. R. Dando and A. J. Southward, “Chemoautotrophy in bivalve molluscs of the genus Thyasira,” Journal of the Marine Biological Association of the United Kingdom, vol. 66, pp. 915–929, 1986. View at: Google Scholar
  50. F. B. Mirza and J. S. Gray, “The fauna of benthic sediments from the organically enriched Oslofjord, Norway,” Journal of Experimental Marine Biology and Ecology, vol. 54, no. 2, pp. 181–207, 1981. View at: Google Scholar
  51. E. Bagøien, S. Kaartvedt, and S. Øverås, “Seasonal vertical migrations of Calanus spp. in Oslofjorden,” Sarsia, vol. 85, no. 4, pp. 299–311, 2000. View at: Google Scholar
  52. J. A. Berge, T. Andersen, R. Amundsen et al., OvervåknIng Av ForurensnIngssItuasjonen I Indre Oslofjord 2008, Norsk Insititutt for Vannforskning (NIVA), Oslo, Norway, 2009.
  53. J. A. Berge, R. Amundsen, K. Bergland et al., OvervåknIng Av Indre Oslofjord I 2011-Vedleggsrapport, Norsk Insititutt for Vannforskning (NIVA), Oslo, Norway, 2012.
  54. L. Buhl-Mortensen, E. Oug, and J. Aure, “The response of hyperbenthos and infauna to hypoxia in fjords along the Skagerrak: estimating loss of biodiversity due to eutrophication,” in Integrated Coastal Zone Management, E. Moksnes, E. Dahl, and J. Støttrup, Eds., pp. 79–96, Blackwell, 2009. View at: Google Scholar
  55. A. B. Josefson and B. Widbom, “Differential response of benthic macrofauna and meiofauna to hypoxia in the Gullmar Fjord basin,” Marine Biology, vol. 100, no. 1, pp. 31–40, 1988. View at: Publisher Site | Google Scholar
  56. R. C. Grabowski, I. G. Droppo, and G. Wharton, “Erodibility of cohesive sediment: the importance of sediment properties,” Earth-Science Reviews, vol. 105, no. 3-4, pp. 101–120, 2011. View at: Publisher Site | Google Scholar
  57. S. J. Hall, “Physical disturbance and marine benthic communities: life in unconsolidated sediments,” Oceanography and Marine Biology, vol. 32, pp. 179–239, 1994. View at: Google Scholar
  58. T. Bekkby, P. E. Isachsen, M. Isæus, and V. Bakkestuen, “GIS modeling of wave exposure at the seabed: a depth-attenuated wave exposure model,” Marine Geodesy, vol. 31, no. 2, pp. 117–127, 2008. View at: Publisher Site | Google Scholar
  59. F. Hjulström, Studies of the Morphological Activity of Rivers as Illustrated by the River Fyris: Inaugural Dissertation, Almqvist & Wiksells Boktryckeri, 1935.
  60. J. S. Gray, “The ecology of marine sediments: an introduction to the structure and function of benthic communities,” Cambridge Studies in Modern Biology, vol. 2, pp. 1–185, 1981. View at: Google Scholar
  61. S. F. Thrush, “Spatial patterns in soft-bottom communities,” Trends in Ecology and Evolution, vol. 6, no. 3, pp. 75–79, 1991. View at: Google Scholar
  62. C. Kraan, J. van der Meer, A. Dekinga, and T. Piersma, “Patchiness of macrobenthic invertebrates in homogenized intertidal habitats: hidden spatial structure at a landscape scale,” Marine Ecology Progress Series, vol. 383, pp. 211–224, 2009. View at: Publisher Site | Google Scholar
  63. B. P. Kelaher and J. S. Levinton, “Variation in detrital enrichment causes spatio-temporal variation in soft-sediment assemblages,” Marine Ecology Progress Series, vol. 261, pp. 85–97, 2003. View at: Google Scholar
  64. G. Wieking and I. Kröncke, “Is benthic trophic structure affected by food quality? The Dogger Bank example,” Marine Biology, vol. 146, no. 2, pp. 387–400, 2005. View at: Publisher Site | Google Scholar
  65. P. A. Ramey and P. V. R. Snelgrove, “Spatial patterns in sedimentary macrofaunal communities on the south coast of Newfoundland in relation to surface oceanography and sediment characteristics,” Marine Ecology Progress Series, vol. 262, pp. 215–227, 2003. View at: Google Scholar
  66. W. G. Ambrose Jr. and P. E. Renaud, “Benthic response to water column productivity patterns: evidence for benthic-pelagic coupling in the Northeast Water Polynya,” Journal of Geophysical Research, vol. 100, no. 3, pp. 4411–4421, 1995. View at: Publisher Site | Google Scholar
  67. L. Buhl-Mortensen and T. Høisæter, “Mollusc fauna along an offshore-fjord gradient,” Marine Ecology Progress Series, vol. 97, no. 3, pp. 209–224, 1993. View at: Google Scholar
  68. R. C. Newell, L. J. Seiderer, and J. E. Robinson, “Animal: sediment relationships in coastal deposits of the eastern english channel,” Journal of the Marine Biological Association of the United Kingdom, vol. 81, no. 1, pp. 1–9, 2001. View at: Publisher Site | Google Scholar
  69. H. C. Trannum, H. C. Nilsson, M. T. Schaanning, and K. Norling, “Biological and biogeochemical effects of organic matter and drilling discharges in two sediment communities,” Marine Ecology Progress Series, vol. 442, pp. 23–36, 2011. View at: Publisher Site | Google Scholar
  70. P. H. Klopfer and R. H. MacArthur, “Niche size and faunal diversity,” The American Naturalist, vol. 94, pp. 293–300, 1960. View at: Google Scholar
  71. D. E. Goldberg, “Components of resource competition in plant communities,” in Perspectives on Plant Competition, J. B. Grace and D. Tilman, Eds., pp. 27–49, Academic Press, San Diego, Calif, USA, 1990. View at: Google Scholar
  72. S. A. Cushman and K. McGarigal, “Patterns in the species-environment relationship depend on both scale and choice of response variables,” Oikos, vol. 105, no. 1, pp. 117–124, 2004. View at: Publisher Site | Google Scholar
  73. J. L. Rueda and C. Salas, “Molluscs associated with a subtidal Zostera marina L. bed in southern Spain: linking seasonal changes of fauna and environmental variables,” Estuarine, Coastal and Shelf Science, vol. 79, no. 1, pp. 157–167, 2008. View at: Publisher Site | Google Scholar
  74. J. B. Wilson, “Species presence/absence sometimes represents a plant community as well as species abundances do, or better,” Journal of Vegetation Science, vol. 23, pp. 1013–1023, 2012. View at: Google Scholar
  75. E. Oug, A. Fleddum, B. Rygg, and F. Olsgard, “Biological traits analyses in the study of pollution gradients and ecological functioning of marine soft bottom species assemblages in a fjord ecosystem,” Journal of Experimental Marine Biology and Ecology, vol. 432-433, pp. 94–105, 2012. View at: Google Scholar
  76. W. Appeltans, P. Bouchet, G. A. Boxshall et al., Eds., “World register of marine species,” 2013, View at: Google Scholar

Copyright © 2014 Thijs Christiaan van Son et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Related articles

No related content is available yet for this article.
 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

No related content is available yet for this article.

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.