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

Authors: Thijs Christiaan van Son; Rune Halvorsen; Karl Norling; Torgeir Bakke; Maria Kaurin; Fredrik Melsom
Journal: Journal of Marine Biology
 (2014)

---

## Abstract

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 termedPeriodic 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.


---

## Body

## 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 [3–5], salinity [6, 7], temperature [8–10], 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 Figure1) 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.Figure 1
Map showing the location of the study area close to Oslo, Norway.


### 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 r.in.xyz 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 (65×65 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).Table 1
Sampling of sites based on stratified random sampling of geomorphological features. The classification into geomorphological features was based on a slightly modified version of an algorithm by Lundblad et al. [23], which takes into account fine- and broad-scaled geomorphological variation.
Geomorphology
Cover (%)
Sites
BPI65
BPI33
Slope65
Slope33
Depth
Fine-scaled features
Narrow depression
6.4
1-2, 26
≤−1 SD
≤–1 SD
Broad depression with open bottom
3.2
3-4
≤−1 SD
w/in 1 SD
Narrow crest
7.6
5-6, 25
≥1 SD
≥1 SD
Local depression on flat
0.1
7
w/in 1 SD
≤–1 SD
≤5
Lateral midslope depression
1.4
8-9
w/in 1 SD
≤–1 SD
>5
Local crest on flat
0.1
10
w/in 1 SD
≥1 SD
≤5
Lateral midslope crest
1.2
11-12
w/in 1 SD
≥1 SD
>5
Shallow flat
28.4
13–15, 27
>−1 SD
w/in 1 SD
≤5
≤58.6
Deep flat
15.7
16–18
>−1 SD
w/in 1 SD
≤5
>58.6
Shallow slope
21.9
19–21, 28
>−1 SD
w/in 1 SD
>5 and  ≤70
≤58.6
Deep slope
14.0
22–24
>−1 SD
w/in 1 SD
>5 and ≤70
>58.6
Total
100.0
Broad-scaled features
Depression
9.6
1–4, 26
≤−1 SD
Crest
11.1
5-6, 25
≥1 SD
Deep flat
16.1
7, 11, 16–18
w/in 1 SD
≤5
>58.6
Shallow slope
19.6
8-9, 19–21, 28
w/in 1 SD
>5
≤58.6
Shallow flat
29.4
10, 13–15, 27
w/in 1 SD
≤5
≤58.6
Deep slope
14.2
12, 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.Bathymetry raster (a) showing the depth in meters, map of fine-scaled geomorphological features (b), and broad-scaled geomorphological features (c) in the study area. The maps of the geomorphological features were obtained by applying a modified version of the classification algorithm in Lundblad et al. [23]. The numbers indicate the location of the 28 sampled sites. All maps are scaled in UTM 32N units.
(a)
(b)
(c)Of the 14 fine-scaled geomorphological features recognised by the classification scheme, 11 were present in the study area (Table1), 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 Table1 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 Table2). 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).Table 2
Recorded environmental variables: group affiliation, name, description, range, and mean (of raw, untransformed values).
Variable
Description
Range
Mean
DEM derived
Aspect_E
Eastness component of aspect
−1.00–1.00
0.10
Aspect_N
Northness component of aspect
−0.96–0.92
–0.35
BPI
Bathymetric position index
−1.52–1.34
–0.06
Depth [m]
Depth derived from DEM
35.8–93.2
58.6
Seabed_curv
Maximum curvature of the seabed (in any direction)
−0.011–0.042
0.005
Slope [deg]
The slope of the seabed
0.3–19.8
6.7
Current model
Current_bottom [m/s]
Maximum current speed at the bottom
0.026–0.077
0.049
Current_surface [m/s]
Maximum current speed at the ocean surface
0.061–0.091
0.068
Change_curr_bottom
The rate of instantaneous change of max bottom current
−1.71e-3–1.18e-3
1.96e-6
Change_curr_surface
The rate of instantaneous change of max surface current
−1.18e-3–3.61e-4
−1.92e-4
Direction_curr_E
Eastness component of the direction of the bottom current
−1.00–0.97
–0.47
Direction_curr_N
Northness component of the direction of the bottom current
−1.00–1.00
0.24
Particle-size
Mud content [%]
The silt-clay fraction of the sediment
8.4–77.2
55.6
Sed_kurt
Kurtosis of the sediment particle distribution
0.70–1.31
0.94
Sed_med [μm]
Median sediment particle-size
28.0–944.9
97.1
Sed_skew
Skewness of the sediment particle-size distribution
−0.78–0.36
0.09
Sed_sort
Sorting of the sediment particle-size distribution
4.2–8.6
5.8
Other variables
BMD [cm]
Biogenic mixing depth derived from sediment profile images
0.0–3.1
1.7
Penetration_depth [cm]
Penetration depth of SPI—a proxy for sediment shear strength
0.0–19.4
12.2
CN
Carbon-to-nitrogen-ratio (weight/weight)
7.9–19.3
11.5
TOC [%]
Percentage of total organic carbon
1.7–5.6
3.8
TN [%]
Percentage of organic nitrogen
0.12–0.63
0.34
VEAS_dist [m]
Distance from effluent of VEAS wastewater treatment plant
1143–2120
1640Three 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 R 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 R values from qualitative (presence/absence; R=1) 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 R values between R=1 and R=max⁡ (R=1, 2, 4, 8,…, 2n, …, 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 R 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 (R=1, presence-absence), intermediate weight given to abundance (intermediate R, see selection criteria below), and strong weight given to abundance (R=max⁡, 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 r 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.Table 3
Procrustes correlation coefficients (r) between parallel DCA and GNMDS ordinations, Kendall’s tau correlations (τ1 and τ2) between corresponding axes for each abundance range (R1–Rmax⁡), and percent stress for GNMDS ordinations. The abundance ranges indicated by bold-face numbers represent no (R1), intermediate (R32), and a strong (Rmax⁡ = 571) weight given to abundance (see Section 2.6 for explanation of selection criteria).
Range
R
τ
1
τ
2
Stress
R
1 (qual)
0.779
0.831
0.243
20.3
R
2
0.821
0.905
0.249
19.1
R
4
0.838
0.889
0.296
17.6
R
8
0.844
0.852
0.291
15.7
R
16
0.865
0.862
0.339
14.4
R
32
0.883
0.873
0.418
13.6
R
64
0.904
0.873
0.481
13.2
R
128
0.932
0.889
0.481
13.4
R
256
0.935
0.894
0.481
13.4
R
max
⁡ (raw)
0.923
0.915
0.423
12.9DCA 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 =1·e-7, 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 P value of approximately 0.10 for this data set with n=28) 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].


## 2.1. Study Area

The study took place in the western part of the Inner Oslofjord (commonly referred to as the Vestfjord; see Figure1) 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.Figure 1
Map showing the location of the study area close to Oslo, Norway.


## 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 r.in.xyz 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 (65×65 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).Table 1
Sampling of sites based on stratified random sampling of geomorphological features. The classification into geomorphological features was based on a slightly modified version of an algorithm by Lundblad et al. [23], which takes into account fine- and broad-scaled geomorphological variation.
Geomorphology
Cover (%)
Sites
BPI65
BPI33
Slope65
Slope33
Depth
Fine-scaled features
Narrow depression
6.4
1-2, 26
≤−1 SD
≤–1 SD
Broad depression with open bottom
3.2
3-4
≤−1 SD
w/in 1 SD
Narrow crest
7.6
5-6, 25
≥1 SD
≥1 SD
Local depression on flat
0.1
7
w/in 1 SD
≤–1 SD
≤5
Lateral midslope depression
1.4
8-9
w/in 1 SD
≤–1 SD
>5
Local crest on flat
0.1
10
w/in 1 SD
≥1 SD
≤5
Lateral midslope crest
1.2
11-12
w/in 1 SD
≥1 SD
>5
Shallow flat
28.4
13–15, 27
>−1 SD
w/in 1 SD
≤5
≤58.6
Deep flat
15.7
16–18
>−1 SD
w/in 1 SD
≤5
>58.6
Shallow slope
21.9
19–21, 28
>−1 SD
w/in 1 SD
>5 and  ≤70
≤58.6
Deep slope
14.0
22–24
>−1 SD
w/in 1 SD
>5 and ≤70
>58.6
Total
100.0
Broad-scaled features
Depression
9.6
1–4, 26
≤−1 SD
Crest
11.1
5-6, 25
≥1 SD
Deep flat
16.1
7, 11, 16–18
w/in 1 SD
≤5
>58.6
Shallow slope
19.6
8-9, 19–21, 28
w/in 1 SD
>5
≤58.6
Shallow flat
29.4
10, 13–15, 27
w/in 1 SD
≤5
≤58.6
Deep slope
14.2
12, 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.Bathymetry raster (a) showing the depth in meters, map of fine-scaled geomorphological features (b), and broad-scaled geomorphological features (c) in the study area. The maps of the geomorphological features were obtained by applying a modified version of the classification algorithm in Lundblad et al. [23]. The numbers indicate the location of the 28 sampled sites. All maps are scaled in UTM 32N units.
(a)
(b)
(c)Of the 14 fine-scaled geomorphological features recognised by the classification scheme, 11 were present in the study area (Table1), 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 Table1 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 Table2). 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).Table 2
Recorded environmental variables: group affiliation, name, description, range, and mean (of raw, untransformed values).
Variable
Description
Range
Mean
DEM derived
Aspect_E
Eastness component of aspect
−1.00–1.00
0.10
Aspect_N
Northness component of aspect
−0.96–0.92
–0.35
BPI
Bathymetric position index
−1.52–1.34
–0.06
Depth [m]
Depth derived from DEM
35.8–93.2
58.6
Seabed_curv
Maximum curvature of the seabed (in any direction)
−0.011–0.042
0.005
Slope [deg]
The slope of the seabed
0.3–19.8
6.7
Current model
Current_bottom [m/s]
Maximum current speed at the bottom
0.026–0.077
0.049
Current_surface [m/s]
Maximum current speed at the ocean surface
0.061–0.091
0.068
Change_curr_bottom
The rate of instantaneous change of max bottom current
−1.71e-3–1.18e-3
1.96e-6
Change_curr_surface
The rate of instantaneous change of max surface current
−1.18e-3–3.61e-4
−1.92e-4
Direction_curr_E
Eastness component of the direction of the bottom current
−1.00–0.97
–0.47
Direction_curr_N
Northness component of the direction of the bottom current
−1.00–1.00
0.24
Particle-size
Mud content [%]
The silt-clay fraction of the sediment
8.4–77.2
55.6
Sed_kurt
Kurtosis of the sediment particle distribution
0.70–1.31
0.94
Sed_med [μm]
Median sediment particle-size
28.0–944.9
97.1
Sed_skew
Skewness of the sediment particle-size distribution
−0.78–0.36
0.09
Sed_sort
Sorting of the sediment particle-size distribution
4.2–8.6
5.8
Other variables
BMD [cm]
Biogenic mixing depth derived from sediment profile images
0.0–3.1
1.7
Penetration_depth [cm]
Penetration depth of SPI—a proxy for sediment shear strength
0.0–19.4
12.2
CN
Carbon-to-nitrogen-ratio (weight/weight)
7.9–19.3
11.5
TOC [%]
Percentage of total organic carbon
1.7–5.6
3.8
TN [%]
Percentage of organic nitrogen
0.12–0.63
0.34
VEAS_dist [m]
Distance from effluent of VEAS wastewater treatment plant
1143–2120
1640Three 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 R 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 R values from qualitative (presence/absence; R=1) 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 R values between R=1 and R=max⁡ (R=1, 2, 4, 8,…, 2n, …, 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 R 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 (R=1, presence-absence), intermediate weight given to abundance (intermediate R, see selection criteria below), and strong weight given to abundance (R=max⁡, 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 r 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.Table 3
Procrustes correlation coefficients (r) between parallel DCA and GNMDS ordinations, Kendall’s tau correlations (τ1 and τ2) between corresponding axes for each abundance range (R1–Rmax⁡), and percent stress for GNMDS ordinations. The abundance ranges indicated by bold-face numbers represent no (R1), intermediate (R32), and a strong (Rmax⁡ = 571) weight given to abundance (see Section 2.6 for explanation of selection criteria).
Range
R
τ
1
τ
2
Stress
R
1 (qual)
0.779
0.831
0.243
20.3
R
2
0.821
0.905
0.249
19.1
R
4
0.838
0.889
0.296
17.6
R
8
0.844
0.852
0.291
15.7
R
16
0.865
0.862
0.339
14.4
R
32
0.883
0.873
0.418
13.6
R
64
0.904
0.873
0.481
13.2
R
128
0.932
0.889
0.481
13.4
R
256
0.935
0.894
0.481
13.4
R
max
⁡ (raw)
0.923
0.915
0.423
12.9DCA 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 =1·e-7, 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 P value of approximately 0.10 for this data set with n=28) 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).Table 4
Frequency and average abundance in presence sites per species as a function of affiliation to broad-scaled features.
Frequency
Average abundance in presence sites
freq
occ
ab.avr
ab.tot
C
SS
SF
D
DS
DF
C
SS
SF
D
DS
DF
Polychaeta
Pherusa falcata
1.00
1.67
0.11
3
1.67
5
Eumida sanguinea
0.67
0.17
2.50
1.00
0.11
3
2.00
6
Eunice pennata
0.67
0.17
2.50
2.00
0.11
3
2.33
7
Streblosoma bairdi
0.33
0.50
0.80
0.40
0.50
1.00
2.33
1.50
2.00
1.50
0.43
12
1.75
21
Ophelina acuminata
0.33
0.33
0.20
0.40
0.25
2.00
2.00
1.00
1.00
1.00
0.25
7
1.43
10
Maldane sarsi
0.33
0.17
0.60
0.40
0.25
9.00
5.00
3.67
8.50
3.00
0.29
8
5.62
45
Anobothrus gracilis
1.00
1.00
1.00
1.00
1.00
1.00
51.33
18.67
7.60
6.40
6.50
3.80
1.00
28
13.61
381
Glycera alba
1.00
1.00
1.00
1.00
1.00
1.00
10.00
8.83
10.60
11.20
8.25
6.20
1.00
28
9.14
256
Goniada maculata
1.00
1.00
1.00
1.00
1.00
1.00
17.67
11.00
12.00
14.00
12.75
10.60
1.00
28
12.61
353
Scoletoma fragilis
1.00
1.00
1.00
1.00
1.00
1.00
11.67
11.50
14.60
17.40
10.25
12.20
1.00
28
13.07
366
Spiophanes kroyeri
1.00
1.00
1.00
1.00
1.00
1.00
35.33
47.00
33.60
141.40
94.25
149.80
1.00
28
85.32
2389
Typosyllis cornuta
1.00
1.00
1.00
1.00
1.00
1.00
13.33
7.83
11.40
10.00
7.00
8.20
1.00
28
9.39
263
Pholoe inornata
1.00
1.00
1.00
1.00
1.00
0.80
5.67
2.50
4.00
3.00
2.25
3.75
0.96
27
3.37
91
Polycirrus norvegicus
1.00
1.00
1.00
1.00
1.00
0.80
6.00
7.67
7.80
9.20
2.00
5.50
0.96
27
6.63
179
Chaetozone setosa
1.00
1.00
1.00
1.00
1.00
0.60
28.00
8.33
10.60
12.80
53.00
23.67
0.93
26
20.54
534
Pista cristata
1.00
1.00
1.00
1.00
0.75
1.00
14.00
12.83
13.60
12.40
9.67
9.40
0.96
27
12.04
325
Polyphysia crassa
1.00
1.00
1.00
1.00
0.75
0.80
23.33
16.83
7.80
6.40
4.00
8.75
0.93
26
11.12
289
Heteromastus filiformis
1.00
1.00
0.80
1.00
1.00
1.00
19.33
5.67
7.50
46.40
44.00
11.00
0.96
27
21.67
585
Terebellides stroemi
1.00
1.00
0.80
1.00
1.00
1.00
9.33
8.83
9.00
6.20
10.00
4.60
0.96
27
7.81
211
Pseudopolydora paucibranchiata
1.00
1.00
0.40
0.40
0.75
0.60
194.33
9.50
30.50
6.00
6.33
4.67
0.68
19
39.26
746
Brada villosa
1.00
1.00
0.40
0.40
0.25
0.40
1.33
4.50
5.00
2.00
3.00
4.00
0.57
16
3.50
56
Prionospio cirrifera
1.00
0.83
0.60
1.00
1.00
1.00
60.67
7.00
21.67
14.60
10.50
8.20
0.89
25
17.52
438
Trichobranchus roseus
1.00
0.83
0.60
0.60
0.75
0.80
1.67
3.40
3.67
3.33
1.67
3.50
0.75
21
2.95
62
Prionospio fallax
1.00
0.67
1.00
0.80
1.00
0.80
25.33
16.75
8.00
4.25
11.00
5.00
0.86
24
11.00
264
Diplocirrus glaucus
1.00
0.67
0.80
1.00
1.00
0.60
3.67
3.75
4.25
3.60
13.50
4.00
0.82
23
5.52
127
Polycirrus medusa
1.00
0.50
0.40
0.60
0.50
0.60
2.67
1.67
3.50
1.00
12.00
1.33
0.57
16
3.19
51
Melinna elisabethae
1.00
0.50
0.20
0.60
0.50
0.60
6.67
5.67
1.00
1.67
1.00
3.00
0.54
15
3.60
54
Spiochaetopterus typicus
0.67
1.00
1.00
1.00
1.00
1.00
1.00
4.00
6.60
3.60
7.25
7.20
0.96
27
5.26
142
Aphrodita aculeata
0.67
0.67
0.40
0.20
0.50
0.40
3.00
1.75
1.00
4.00
1.00
1.50
0.46
13
1.85
24
Melinna cristata
0.67
0.50
0.40
0.80
1.00
1.00
2.50
2.67
3.00
4.25
2.75
3.80
0.71
20
3.30
66
Syllidia armata
0.67
0.33
0.80
0.60
0.25
0.60
2.50
1.50
1.25
2.00
1.00
1.33
0.54
15
1.60
24
Scalibregma inflatum
0.67
0.33
0.40
0.60
0.50
0.40
1.50
1.00
1.00
1.00
1.00
1.00
0.46
13
1.08
14
Amaena trilobata
0.67
0.17
0.40
0.40
0.50
0.80
2.00
2.00
1.00
5.00
1.00
1.50
0.46
13
2.00
26
Trichobranchus glacialis
0.33
0.67
0.20
0.40
1.00
0.60
1.00
2.25
5.00
3.00
1.25
1.33
0.54
15
2.00
30
Pista lornensis
0.33
0.17
0.60
0.40
0.75
0.80
1.00
1.00
1.67
1.50
1.00
1.25
0.50
14
1.29
18
Sige fusigera
0.33
0.17
0.20
0.40
0.25
0.40
1.00
1.00
3.00
2.50
1.00
1.00
0.29
8
1.62
13
Eupolymnia nebulosa
0.67
0.33
0.20
0.20
0.20
3.50
2.00
2.00
1.00
1.00
0.25
7
2.14
15
Glycera lapidum
1.00
0.33
0.40
0.25
0.20
11.67
6.50
1.50
1.00
1.00
0.32
9
5.89
53
Cossura longocirrata
0.67
0.17
0.20
0.25
1.50
1.00
1.00
1.00
0.18
5
1.20
6
Sphaerodorum gracilis
0.67
0.33
0.20
2.50
1.00
1.00
0.18
5
1.60
8
Pholoe assimilis
0.33
0.17
0.20
4.00
1.00
1.00
0.11
3
2.00
6
Chaetoparia nilssoni
0.33
0.60
0.50
1.00
1.00
1.00
0.21
6
1.00
6
Scoloplos armiger
0.67
0.20
19.00
3.00
0.11
3
13.67
41
Scoletoma cf. magnidentata
0.33
0.40
0.40
0.50
0.20
1.00
1.00
1.00
1.00
1.00
0.32
9
1.00
9
Amphicteis gunneri
0.17
0.40
0.20
0.25
0.40
2.00
1.00
1.00
2.00
2.00
0.25
7
1.57
11
Ophiodromus flexuosus
0.33
0.20
0.20
0.25
1.50
1.00
2.00
1.00
0.18
5
1.40
7
Phyllodoce rosea
0.33
0.40
0.50
0.40
1.50
1.00
2.00
1.50
0.29
8
1.50
12
Phyllodoce groenlandica
0.33
0.20
0.20
1.00
1.00
1.00
0.14
4
1.00
4
Dipolydora cf. coeca
0.17
11.00
0.04
1
11.00
11
Petaloproctus borealis
0.17
5.00
0.04
1
5.00
5
Sabellides octocirrata
0.20
0.40
0.25
0.20
1.00
1.00
1.00
1.00
0.18
5
1.00
5
Orbinia sertulata
0.20
0.40
0.25
2.00
1.00
1.00
0.14
4
1.25
5
Nephtys ciliata
0.20
0.50
0.40
1.00
2.00
1.00
0.18
5
1.40
7
Glycera unicornis
0.60
0.50
1.00
1.00
0.18
5
1.00
5
Rhodine loveni
0.20
5.00
0.04
1
5.00
5
Mollusca
Hiatella arctica
0.67
0.33
4.00
2.00
0.14
4
3.00
12
Lepeta caeca
0.67
0.17
20.00
6.00
0.11
3
15.33
46
Mya truncata
0.67
0.17
1.50
1.00
0.11
3
1.33
4
Arctica islandica
0.33
0.50
0.20
1.00
1.67
1.00
0.18
5
1.40
7
Thyasira equalis
1.00
1.00
1.00
1.00
0.75
1.00
9.00
14.67
12.20
17.20
37.00
18.40
0.96
27
17.22
465
Astarte montagui
1.00
1.00
0.80
0.80
0.75
0.60
13.00
5.00
6.25
5.25
6.33
3.00
0.82
23
6.22
143
Caudofoveata
1.00
0.83
0.60
0.80
0.50
0.40
1.67
1.80
2.00
1.75
1.50
2.00
0.68
19
1.79
34
Ennucula tenuis
1.00
0.67
1.00
1.00
1.00
1.00
39.00
68.00
72.40
106.80
130.00
82.20
0.93
26
85.23
2216
Pseudamussium peslutrae
1.00
0.50
0.40
0.20
0.25
0.20
2.67
2.00
1.00
1.00
1.00
1.00
0.39
11
1.73
19
Thyasira sarsi
0.67
0.83
0.80
1.00
1.00
1.00
3.50
3.20
6.75
10.00
11.25
12.40
0.89
25
8.28
207
Abra nitida
0.33
0.50
0.40
0.80
0.50
1.00
5.00
2.67
2.00
3.00
5.00
2.60
0.61
17
3.06
52
Thyasira flexuosa
0.33
0.33
0.80
1.00
1.00
1.00
3.00
7.00
5.50
11.40
11.50
12.80
0.75
21
9.81
206
Philine scabra
0.33
0.33
0.20
0.20
0.25
0.20
2.00
1.00
1.00
2.00
2.00
1.00
0.25
7
1.43
10
Corbula gibba
1.00
0.83
0.20
0.40
7.00
20.60
5.00
1.50
0.39
11
12.00
132
Parvicardium minimum
0.67
0.17
0.20
1.00
1.00
1.00
0.14
4
1.00
4
Macoma calcarea
0.67
0.40
0.40
0.50
0.60
7.00
6.00
2.00
7.00
8.67
0.39
11
6.36
70
Astarte elliptica
0.67
0.20
0.60
0.20
3.00
2.00
3.33
1.00
0.25
7
2.71
19
Thyasira gouldi
0.33
0.20
0.40
0.80
1.00
1.00
3.00
2.00
0.32
9
1.89
17
Nuculana pernula
0.33
0.25
1.00
2.00
0.11
3
1.33
4
Buccinum undatum
0.40
0.40
0.25
0.40
1.00
1.00
1.00
1.00
0.25
7
1.00
7
Echinodermata
Strongylocentrotus droebachiensis
0.67
2.00
0.07
2
2.00
4
Amphiura chiajei
0.67
0.83
0.40
1.50
3.60
3.50
0.32
9
3.11
28
Ophiocten affinis
1.00
1.00
0.60
0.60
0.50
0.20
13.33
6.17
4.33
3.67
1.00
2.00
0.64
18
5.83
105
Ophiura albida
1.00
0.50
0.40
0.40
0.50
0.40
2.00
1.33
1.00
1.00
2.50
2.00
0.50
14
1.64
23
Amphiura filiformis
0.67
0.50
0.60
0.40
0.25
0.40
1.50
5.67
3.00
1.00
7.00
2.00
0.46
13
3.23
42
Ophiura robusta
0.67
0.33
0.20
2.00
1.00
1.00
0.18
5
1.40
7
Amphipholis squamata
0.33
0.67
0.40
0.20
0.20
2.00
1.25
1.00
4.00
2.00
0.32
9
1.67
15
Brissopsis lyrifera
0.67
0.40
0.20
1.00
1.00
1.00
0.18
5
1.00
5
Echinocardium flavescens
0.33
0.40
0.40
0.25
0.60
1.00
1.50
1.00
1.00
1.67
0.36
10
1.30
13
Polychaeta
Euchone rubrocinta
0.33
1.00
0.04
1
1.00
1
Exogone verugera
0.33
1.00
0.04
1
1.00
1
Harmothoe cf. fragilis
0.33
2.00
0.04
1
2.00
2
Laonice cirrata
0.33
2.00
0.04
1
2.00
2
Nicolea zostericola
0.33
1.00
0.04
1
1.00
1
Pherusa plumosa
0.33
2.00
0.04
1
2.00
2
Pholoe pallida
0.33
1.00
0.04
1
1.00
1
Praxillura longissima
0.33
1.00
0.04
1
1.00
1
Rhoine gracilior
0.33
2.00
0.04
1
2.00
2
Gyptis rosea
0.67
0.20
1.00
1.00
0.11
3
1.00
3
Eteone cf. longa
0.67
0.25
1.00
1.00
0.11
3
1.00
3
Galathowenia oculata
0.17
1.00
0.04
1
1.00
1
Laetmonice filicornis
0.17
1.00
0.04
1
1.00
1
Orbinia norvegica
0.17
1.00
0.04
1
1.00
1
Typosyllis armillaris
0.17
1.00
0.04
1
1.00
1
Parougia eliasoni
0.20
0.25
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Flabelligera affinis
0.20
1.00
0.04
1
1.00
1
Paraamphinome jeffreysii
0.25
0.20
2.00
1.00
0.07
2
1.50
3
Amage auricola
0.25
2.00
0.04
1
2.00
2
Mediomastus fragilis
0.25
1.00
0.04
1
1.00
1
Neoamphitrite grayi
0.20
1.00
0.04
1
1.00
1
Mollusca
Parvicardium pinnulatum
0.33
1.00
0.04
1
1.00
1
Tellimya ferruginosa
0.67
0.17
1.00
1.00
0.11
3
1.00
3
Aporrhais pespelecani
0.33
0.17
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Cuspidaria cuspidata
0.17
1.00
0.04
1
1.00
1
Kurtiella bidentata
0.17
1.00
0.04
1
1.00
1
Abra alba
0.17
0.20
1.00
1.00
0.07
2
1.00
2
Astarte sulcata
0.20
0.20
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Antalis entalis
0.20
1.00
0.04
1
1.00
1
Diaphana minuta
0.20
1.00
0.04
1
1.00
1
Cylichna alba
0.25
1.00
0.04
1
1.00
1
Lunatia montagui
0.20
1.00
0.04
1
1.00
1
Mytilus edulis
0.20
2.00
0.04
1
2.00
2
Rissoidae
0.20
2.00
0.04
1
2.00
2
Echinodermata
Echinocardium cordatum
0.17
0.25
2.00
1.00
0.07
2
1.50
3
Ophiura ophiura
0.17
0.25
1.00
1.00
0.07
2
1.00
2
Ophiopholis aculeata
0.20
2.00
0.04
1
2.00
2
Ophiothrix fragilis
0.20
1.00
0.04
1
1.00
1
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 (R1) and 12.9 (Rmax⁡; see Table 3). The Procrustes analysis revealed high configurational concordance between parallel DCA and GNMDS ordinations: Procrustes correlation coefficients r 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 (R1) representation of the species-by-site matrix, the second axis was not confirmed according to the “rule-of-thumb” of τ≥0.4 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 (Figure3) exhibited a dense clustering of sites at the intermediate abundance range (R32), while, for both the qualitative data set (R1) and for data with strong weight given to abundance (Rmax⁡), 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 R was observed. Gradient lengths of DCA1 increased with increasing weight given to abundance (1.97 SD units for R=1 to 3.00 SD units for R=max⁡). 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.Figure 3
DCA site diagram for axes 1 and 2 for three abundance ranges with no (R1), intermediate (R32), and strong (Rmax⁡) weight given to abundance. Positions of sampled sites 1–28, with affiliation to broad-scaled geomorphological feature indicated by the color of the dot, are shown. DCA axes are scaled in SD units of compositional turnover.


### 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 (Figure4; 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 (R1), for example, BPI (DCA1) and sediment sorting (DCA2), while other variables were most strongly correlated with DCA axes obtained with intermediate weighting of abundance (R32), 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.Table 5
Kendall'sτ correlations between environmental variables and DCA axes 1 and 2 with no (R1), intermediate (R32), and strong (Rmax⁡) weight given to abundance.
Variable/range
DCA1
DCA2
R
1 (qual)
R
32
R
max
⁡ (raw)
R
1 (qual)
R
32
R
max
⁡ (raw)
DEM derived
Aspect_E
0.011
−0.026
−0.090
0.148
0.021
0.196
Aspect_N
0.101
0.032
−0.074
−0.069
−0.143
−0.085
BPI
0.228
0.169
0.148
−0.016
−0.164
−0.201
Depth
−0.656
−0.820
−0.757
−0.053
−0.063
0.079
Seabed_curv
0.169
0.132
0.079
−0.127
−0.169
−0.312
Slope
0.127
−0.005
−0.153
−0.148
−0.021
−0.132
Current model
Current_bottom
−0.026
−0.095
−0.085
−0.037
−0.143
0.085
Current_surface
−0.307
−0.407
−0.471
−0.212
−0.106
−0.090
Change_curr_bottom
0.074
0.101
0.132
−0.032
−0.053
−0.005
Change_curr_surface
−0.058
−0.138
−0.116
0.090
0.143
0.254
Direction_curr_E
0.024
−0.024
−0.083
−0.008
−0.110
−0.143
Direction_curr_N
0.179
0.094
0.083
−0.035
−0.265
−0.099
Sediment features
Mud content
0.053
0.122
0.069
−0.116
0.095
0.164
Sed_kurt
0.240
0.272
0.245
−0.191
−0.046
−0.073
Sed_med
−0.027
−0.074
−0.095
0.090
0.000
−0.032
Sed_skew
0.099
0.184
0.152
−0.056
0.088
0.318
Sed_sort
0.130
0.135
0.114
0.227
0.146
0.205
Other variables
BMD
−0.412
−0.434
−0.445
0.233
0.228
0.385
Penetration_depth
−0.174
−0.227
−0.286
0.184
0.067
0.206
CN
0.265
0.286
0.318
−0.170
−0.202
−0.260
TOC
−0.202
−0.303
−0.330
−0.245
−0.250
−0.106
TN
−0.257
−0.349
−0.393
−0.149
−0.149
0.008
VEAS_dist
0.349
0.407
0.407
0.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 τ=0 when n=28. The unconfirmed second axis (DCA2.R1) according to criteria given in Section 2.6 is shown in italic.Figure 4
DCA diagrams for axes 1 and 2 with vector fitted environmental variables related to the three selected abundance ranges (R1, R32, and Rmax⁡). Variables were fitted to the DCA axes by application of the envfit function with selection criteria |τ| ≥ 0.225. The length of the vectors indicates their degree of correlation. Sites are shown as black dots. Environmental variables are described in Table 2. DCA axes are scaled in SD units of compositional turnover.The DCA ordination with intermediate weighting of abundance (R32) 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).Table 6
Kendall’s rank correlation coefficientsτ between environmental variables and site scores along DCA axis 1, calculated separately for subsets of sites near the positive (Subset+; n=8) and negative ends (Subset−; n=22; sites 4 and 6 are included in both subsets) of the axis. Correlation coefficients τ > 0.6 are shown in bold. See Table 2 for description of environmental variables.
Var
Subset+
Subset−
τ
P
τ
P
DEM derived
Aspect_E
–0.21
0.548
–0.01
0.956
Aspect_N
0.14
0.720
–0.04
0.824
BPI
0.07
0.905
0.19
0.239
Depth
–0.93
0.000
–0.71
0.000
Seabed_curv
0.07
0.905
0.14
0.371
Slope
0.07
0.905
–0.02
0.911
Current model
Current_bottom
0.14
0.720
–0.06
0.696
Current_surface
–0.50
0.109
–0.19
0.239
Change_curr_bottom
–0.07
0.905
0.13
0.434
Change_curr_surface
–0.43
0.179
–0.22
0.160
Direction_curr_E
0.11
0.708
–0.11
0.479
Direction_curr_N
–0.14
0.720
0.21
0.166
Sediment  features
Mud content
–0.64
0.031
0.28
0.071
Sed_kurt
–0.43
0.179
0.30
0.054
Sed_med
0.64
0.031
–0.20
0.185
Sed.skew
–0.55
0.061
0.38
0.014
Sed.sort
0.11
0.708
0.14
0.381
Other variables
BMD
–0.50
0.109
–0.33
0.036
Penetration_depth
–0.43
0.179
0.00
1.000
CN
0.64
0.031
0.19
0.215
TOC
–0.43
0.179
–0.17
0.283
TN
–0.55
0.061
–0.25
0.113
VEAS_dist
0.50
0.109
0.26
0.091Table 7
Kendall’s rank correlations between selected important environmental variables, calculated separately for subsets of sites near the positive (Subset+;n=8; upper triangle) and negative ends (Subset−; n=22; lower triangle) of the main coenocline (DCA1). Sites 4 and 6 are included in both subsets. Correlation coefficients between variables strongly correlated with depth in Subset+ are shown in bold. See Table 2 for description of environmental variables.
Depth
BMD
Penetration depth
CN
Sed_med
Mud content
TN
TOC
Depth
*****
0.43
0.36
−0.57
−0.71
0.71
0.47
0.36
BMD
0.28
*****
0.93
−0.29
−0.29
0.14
0.11
0.07
Penetration_depth
−0.16
0.38
*****
−0.21
−0.36
0.21
0.04
0.00
CN
−0.17
−0.55
−0.12
*****
0.43
−0.57
−0.33
−0.21
Sed_med
0.14
0.38
0.19
−0.26
*****
−0.86
−0.26
−0.21
Mud content
−0.26
−0.30
−0.01
0.11
−0.68
*****
0.33
0.21
TN
0.31
0.28
−0.12
−0.38
0.35
−0.44
*****
0.91
TOC
0.23
0.04
−0.12
0.09
0.25
−0.41
0.56
*****Table 8
Matrix of Kendall’s correlation coefficientsτ between pairs of environmental variables (upper right triangle) and corresponding  P values for tests of the hypotheses that τ = 0 (lower left triangle). P values ≤ 0.10 and the corresponding τ values are shown in bold.
BPI
Depth
Seabed_curv
Current_bottom
Direction_curr_E
Direction_curr_N
Current_surface
Change_curr_bottom
Sed_kurt
Sed_med
Mud content
Med_skew
Med_sort
BMD
Penetration_depth
TOC
TON
CN
VEAS_dist
BPI
*****
-
0.222
0.646
0.079
-
0.164
0.185
-
0.127
0.201
0.164
-
0.133
0.190
0.163
-
0.032
-
0.141
0.077
−0.223
-
0.100
-
0.138
0.106
Depth
0.101
*****
-
0.175
0.074
0.051
-
0.067
0.481
-
0.090
-
0.213
0.032
-
0.111
-
0.184
-
0.151
0.407
0.131
0.340
0.387
−0.271
−0.439
Seabed_curv
<0.001
0.201
*****
0.053
0.078
0.319
-
0.153
0.090
0.223
-
0.154
0.143
0.019
-
0.211
-
0.222
-
0.051
-
0.069
0.014
-
0.048
0.101
Current_bottom
0.570
0.597
0.710
*****
-
0.073
0.233
-
0.201
0.286
-
0.110
0.159
-
0.138
-
0.104
0.032
-
0.092
-
0.029
0.160
0.089
0.005
0.254
Direction_curr_E
0.227
0.707
0.566
0.593
*****
0.319
0.089
-
0.116
0.077
0.003
-
0.121
−0.234
−0.239
-
0.041
0.027
0.349
0.187
0.256
0.024
Direction_curr_N
0.172
0.621
0.018
0.085
0.020
*****
-
0.120
0.297
0.316
−0.244
0.174
0.095
−0.252
-
0.129
0.022
0.089
0.060
0.078
0.142
Current_surface
0.356
<0.001
0.263
0.139
0.513
0.373
*****
-
0.090
-
0.105
0.011
-
0.079
-
0.136
-
0.070
0.206
0.045
0.383
0.338
-
0.080
−0.524
Change_curr_bottom
0.139
0.518
0.518
0.034
0.394
0.028
0.518
*****
0.202
−0.265
0.196
0.158
−0.232
-
0.173
-
0.045
-
0.101
-
0.149
0.111
0.196
Sed_kurt
0.227
0.117
0.100
0.416
0.578
0.021
0.440
0.137
*****
−0.567
0.568
0.457
−0.256
−0.342
-
0.152
-
0.043
-
0.061
0.067
0.073
Sed_med
0.323
0.813
0.252
0.236
0.984
0.072
0.937
0.048
<0.001
*****
−0.796
−0.522
0.450
0.114
-
0.040
0.176
0.209
-
0.005
-
0.048
Mudcontent
0.162
0.421
0.298
0.317
0.373
0.198
0.570
0.150
<0.001
<0.001
*****
0.724
−0.249
-
0.076
0.136
−0.277
−0.257
-
0.095
0.005
Sed_skew
0.227
0.172
0.890
0.440
0.088
0.488
0.313
0.243
0.001
<0.001
<0.001
*****
0.041
-
0.014
0.248
−0.352
−0.268
-
0.158
-
0.040
Sed_sort
0.812
0.267
0.122
0.812
0.084
0.067
0.606
0.088
0.064
0.001
0.068
0.766
*****
0.122
0.136
-
0.125
-
0.036
-
0.114
-
0.038
BMD
0.302
0.003
0.104
0.500
0.765
0.350
0.132
0.204
0.014
0.404
0.578
0.921
0.381
*****
0.534
0.090
0.258
−0.514
−0.385
Penetration_depth
0.566
0.332
0.707
0.828
0.843
0.874
0.737
0.737
0.266
0.767
0.313
0.068
0.321
<0.001
*****
-
0.013
-
0.011
−0.228
-
0.152
TOC
0.097
0.011
0.607
0.235
0.011
0.513
0.004
0.452
0.751
0.192
0.040
0.009
0.361
0.512
0.921
*****
0.667
0.000
−0.277
TON
0.463
0.005
0.921
0.513
0.176
0.662
0.013
0.275
0.662
0.126
0.059
0.052
0.796
0.064
0.937
<0.001
*****
−0.361
−0.306
CN
0.304
0.044
0.722
0.968
0.060
0.566
0.553
0.406
0.620
0.968
0.477
0.243
0.405
<0.001
0.092
1.000
0.008
*****
0.207
VEAS_dist
0.445
0.001
0.468
0.060
0.858
0.294
<0.001
0.150
0.593
0.722
0.984
0.767
0.781
0.005
0.259
0.040
0.025
0.123
*****Figure 5
Contour-plots of selected environmental variables on DCA.R32, axes 1 and 2, showing patterns of variation in ordination space. The contour-plots are fitted by use of the ordisurf command fitting GAM with Gaussian errors. The size of the bubbles is proportional to the magnitude of the target value and are comparable between plots because environmental variables are ranged between 0 and 1. Axes are scaled in SD units of compositional turnover. The equidistance between contours differs between plots. See Table 2 for description of environmental variables.Sediment profiles images (SPIs) of sites 7 (a) and 25 (b) showing the difference in sediment particle composition and environment at the two extremes of DCA1.
(a)
(b)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 theThyasira 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.Table 9
Species associated with the negative (Neg) and positive (Pos) extremes of DCA axes 1 and 2, with indication of feeding strategy, species loading on axes (“Load” column), and raw total abundance for each species (“Abund” column). Only species with DCA.R32 species loadings with absolute values >1.25 and with a raw total abundance above 10 were chosen.
Axis
End
Load
Species
Su
SD
SSD
G
Sc
P
Abund
DCA1
Neg
−2.00
Thyasira gouldi(M)
x
x
17
−1.96
Sige fusigera(P)
x
x
13
−1.81
Pista lornensis(P)
x
18
−1.47
Thyasira flexuosa(M)
x
x
206
−1.46
Melinna cristata(P)
x
66
−1.36
Thyasira sarsi (M)
x
x
207
−1.34
Ennucula tenuis(M)
x
2216
−1.26
Spiophanes kroyeri(P)
x
x
2389
Pos
2.78
Lepeta caeca(M)
x
x
46
2.59
Hiatella arctica(M)
x
x
12
2.46
Scoloplos armiger(P)
x
41
1.94
Pseudopolydora pauchibranchiata(P)
x
x
746
1.76
Corbula gibba(M)
x
132
1.75
Glycera lapidum(P)
x
x
53
1.66
Eupolymnia nebulosa(P)
x
15
1.61
Amphiura chiajei(E)
x
x
28
1.58
Pseudamussium peslutrae(M)
x
19
1.44
Ophiocten affinis(E)
x
105
1.41
Melinna elisabethae(P)
x
x
54
1.31
Anobothrus gracilis(P)
x
381
DCA2
Neg
−3.38
Scoloplos armiger (P)
x
41
−1.79
Astarte elliptica(M)
x
19
−1.71
Macoma calcarea(M)
x
x
x
70
−1.65
Aphrodita aculeata(P)
x
24
−1.56
Ophiura albida(E)
x
x
x
23
−1.35
Lepeta caeca(M)
x
x
46
Pos
4.81
Dipolydora coeca(P)
x
x
11
1.77
Amphicteis gunneri(P)
x
x
11
1.65
Amphiura chiajei(E)
x
x
28
1.59
Phyllodoce rosea(P)
x
12
1.54
Streblosoma bairdi(P)
x
x
21
1.38
Corbula gibba(M)
x
132
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, P<0.01). 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)).Contour-plots showing patterns of variation in ordination space for (a) total abundance (total number of individuals of all species), (b) species richness, and (c) species evenness for axes 1 and 2 of DCA.R32. The contour-plots are fitted by use of the ordisurf command fitting GAM with Poisson errors. The sizes of the bubbles are proportional to the magnitude of the target value, but not comparable between plots. Axes are scaled in SD units of compositional turnover. The equidistance between contours differs between plots.
(a)
(b)
(c)


## 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).Table 4
Frequency and average abundance in presence sites per species as a function of affiliation to broad-scaled features.
Frequency
Average abundance in presence sites
freq
occ
ab.avr
ab.tot
C
SS
SF
D
DS
DF
C
SS
SF
D
DS
DF
Polychaeta
Pherusa falcata
1.00
1.67
0.11
3
1.67
5
Eumida sanguinea
0.67
0.17
2.50
1.00
0.11
3
2.00
6
Eunice pennata
0.67
0.17
2.50
2.00
0.11
3
2.33
7
Streblosoma bairdi
0.33
0.50
0.80
0.40
0.50
1.00
2.33
1.50
2.00
1.50
0.43
12
1.75
21
Ophelina acuminata
0.33
0.33
0.20
0.40
0.25
2.00
2.00
1.00
1.00
1.00
0.25
7
1.43
10
Maldane sarsi
0.33
0.17
0.60
0.40
0.25
9.00
5.00
3.67
8.50
3.00
0.29
8
5.62
45
Anobothrus gracilis
1.00
1.00
1.00
1.00
1.00
1.00
51.33
18.67
7.60
6.40
6.50
3.80
1.00
28
13.61
381
Glycera alba
1.00
1.00
1.00
1.00
1.00
1.00
10.00
8.83
10.60
11.20
8.25
6.20
1.00
28
9.14
256
Goniada maculata
1.00
1.00
1.00
1.00
1.00
1.00
17.67
11.00
12.00
14.00
12.75
10.60
1.00
28
12.61
353
Scoletoma fragilis
1.00
1.00
1.00
1.00
1.00
1.00
11.67
11.50
14.60
17.40
10.25
12.20
1.00
28
13.07
366
Spiophanes kroyeri
1.00
1.00
1.00
1.00
1.00
1.00
35.33
47.00
33.60
141.40
94.25
149.80
1.00
28
85.32
2389
Typosyllis cornuta
1.00
1.00
1.00
1.00
1.00
1.00
13.33
7.83
11.40
10.00
7.00
8.20
1.00
28
9.39
263
Pholoe inornata
1.00
1.00
1.00
1.00
1.00
0.80
5.67
2.50
4.00
3.00
2.25
3.75
0.96
27
3.37
91
Polycirrus norvegicus
1.00
1.00
1.00
1.00
1.00
0.80
6.00
7.67
7.80
9.20
2.00
5.50
0.96
27
6.63
179
Chaetozone setosa
1.00
1.00
1.00
1.00
1.00
0.60
28.00
8.33
10.60
12.80
53.00
23.67
0.93
26
20.54
534
Pista cristata
1.00
1.00
1.00
1.00
0.75
1.00
14.00
12.83
13.60
12.40
9.67
9.40
0.96
27
12.04
325
Polyphysia crassa
1.00
1.00
1.00
1.00
0.75
0.80
23.33
16.83
7.80
6.40
4.00
8.75
0.93
26
11.12
289
Heteromastus filiformis
1.00
1.00
0.80
1.00
1.00
1.00
19.33
5.67
7.50
46.40
44.00
11.00
0.96
27
21.67
585
Terebellides stroemi
1.00
1.00
0.80
1.00
1.00
1.00
9.33
8.83
9.00
6.20
10.00
4.60
0.96
27
7.81
211
Pseudopolydora paucibranchiata
1.00
1.00
0.40
0.40
0.75
0.60
194.33
9.50
30.50
6.00
6.33
4.67
0.68
19
39.26
746
Brada villosa
1.00
1.00
0.40
0.40
0.25
0.40
1.33
4.50
5.00
2.00
3.00
4.00
0.57
16
3.50
56
Prionospio cirrifera
1.00
0.83
0.60
1.00
1.00
1.00
60.67
7.00
21.67
14.60
10.50
8.20
0.89
25
17.52
438
Trichobranchus roseus
1.00
0.83
0.60
0.60
0.75
0.80
1.67
3.40
3.67
3.33
1.67
3.50
0.75
21
2.95
62
Prionospio fallax
1.00
0.67
1.00
0.80
1.00
0.80
25.33
16.75
8.00
4.25
11.00
5.00
0.86
24
11.00
264
Diplocirrus glaucus
1.00
0.67
0.80
1.00
1.00
0.60
3.67
3.75
4.25
3.60
13.50
4.00
0.82
23
5.52
127
Polycirrus medusa
1.00
0.50
0.40
0.60
0.50
0.60
2.67
1.67
3.50
1.00
12.00
1.33
0.57
16
3.19
51
Melinna elisabethae
1.00
0.50
0.20
0.60
0.50
0.60
6.67
5.67
1.00
1.67
1.00
3.00
0.54
15
3.60
54
Spiochaetopterus typicus
0.67
1.00
1.00
1.00
1.00
1.00
1.00
4.00
6.60
3.60
7.25
7.20
0.96
27
5.26
142
Aphrodita aculeata
0.67
0.67
0.40
0.20
0.50
0.40
3.00
1.75
1.00
4.00
1.00
1.50
0.46
13
1.85
24
Melinna cristata
0.67
0.50
0.40
0.80
1.00
1.00
2.50
2.67
3.00
4.25
2.75
3.80
0.71
20
3.30
66
Syllidia armata
0.67
0.33
0.80
0.60
0.25
0.60
2.50
1.50
1.25
2.00
1.00
1.33
0.54
15
1.60
24
Scalibregma inflatum
0.67
0.33
0.40
0.60
0.50
0.40
1.50
1.00
1.00
1.00
1.00
1.00
0.46
13
1.08
14
Amaena trilobata
0.67
0.17
0.40
0.40
0.50
0.80
2.00
2.00
1.00
5.00
1.00
1.50
0.46
13
2.00
26
Trichobranchus glacialis
0.33
0.67
0.20
0.40
1.00
0.60
1.00
2.25
5.00
3.00
1.25
1.33
0.54
15
2.00
30
Pista lornensis
0.33
0.17
0.60
0.40
0.75
0.80
1.00
1.00
1.67
1.50
1.00
1.25
0.50
14
1.29
18
Sige fusigera
0.33
0.17
0.20
0.40
0.25
0.40
1.00
1.00
3.00
2.50
1.00
1.00
0.29
8
1.62
13
Eupolymnia nebulosa
0.67
0.33
0.20
0.20
0.20
3.50
2.00
2.00
1.00
1.00
0.25
7
2.14
15
Glycera lapidum
1.00
0.33
0.40
0.25
0.20
11.67
6.50
1.50
1.00
1.00
0.32
9
5.89
53
Cossura longocirrata
0.67
0.17
0.20
0.25
1.50
1.00
1.00
1.00
0.18
5
1.20
6
Sphaerodorum gracilis
0.67
0.33
0.20
2.50
1.00
1.00
0.18
5
1.60
8
Pholoe assimilis
0.33
0.17
0.20
4.00
1.00
1.00
0.11
3
2.00
6
Chaetoparia nilssoni
0.33
0.60
0.50
1.00
1.00
1.00
0.21
6
1.00
6
Scoloplos armiger
0.67
0.20
19.00
3.00
0.11
3
13.67
41
Scoletoma cf. magnidentata
0.33
0.40
0.40
0.50
0.20
1.00
1.00
1.00
1.00
1.00
0.32
9
1.00
9
Amphicteis gunneri
0.17
0.40
0.20
0.25
0.40
2.00
1.00
1.00
2.00
2.00
0.25
7
1.57
11
Ophiodromus flexuosus
0.33
0.20
0.20
0.25
1.50
1.00
2.00
1.00
0.18
5
1.40
7
Phyllodoce rosea
0.33
0.40
0.50
0.40
1.50
1.00
2.00
1.50
0.29
8
1.50
12
Phyllodoce groenlandica
0.33
0.20
0.20
1.00
1.00
1.00
0.14
4
1.00
4
Dipolydora cf. coeca
0.17
11.00
0.04
1
11.00
11
Petaloproctus borealis
0.17
5.00
0.04
1
5.00
5
Sabellides octocirrata
0.20
0.40
0.25
0.20
1.00
1.00
1.00
1.00
0.18
5
1.00
5
Orbinia sertulata
0.20
0.40
0.25
2.00
1.00
1.00
0.14
4
1.25
5
Nephtys ciliata
0.20
0.50
0.40
1.00
2.00
1.00
0.18
5
1.40
7
Glycera unicornis
0.60
0.50
1.00
1.00
0.18
5
1.00
5
Rhodine loveni
0.20
5.00
0.04
1
5.00
5
Mollusca
Hiatella arctica
0.67
0.33
4.00
2.00
0.14
4
3.00
12
Lepeta caeca
0.67
0.17
20.00
6.00
0.11
3
15.33
46
Mya truncata
0.67
0.17
1.50
1.00
0.11
3
1.33
4
Arctica islandica
0.33
0.50
0.20
1.00
1.67
1.00
0.18
5
1.40
7
Thyasira equalis
1.00
1.00
1.00
1.00
0.75
1.00
9.00
14.67
12.20
17.20
37.00
18.40
0.96
27
17.22
465
Astarte montagui
1.00
1.00
0.80
0.80
0.75
0.60
13.00
5.00
6.25
5.25
6.33
3.00
0.82
23
6.22
143
Caudofoveata
1.00
0.83
0.60
0.80
0.50
0.40
1.67
1.80
2.00
1.75
1.50
2.00
0.68
19
1.79
34
Ennucula tenuis
1.00
0.67
1.00
1.00
1.00
1.00
39.00
68.00
72.40
106.80
130.00
82.20
0.93
26
85.23
2216
Pseudamussium peslutrae
1.00
0.50
0.40
0.20
0.25
0.20
2.67
2.00
1.00
1.00
1.00
1.00
0.39
11
1.73
19
Thyasira sarsi
0.67
0.83
0.80
1.00
1.00
1.00
3.50
3.20
6.75
10.00
11.25
12.40
0.89
25
8.28
207
Abra nitida
0.33
0.50
0.40
0.80
0.50
1.00
5.00
2.67
2.00
3.00
5.00
2.60
0.61
17
3.06
52
Thyasira flexuosa
0.33
0.33
0.80
1.00
1.00
1.00
3.00
7.00
5.50
11.40
11.50
12.80
0.75
21
9.81
206
Philine scabra
0.33
0.33
0.20
0.20
0.25
0.20
2.00
1.00
1.00
2.00
2.00
1.00
0.25
7
1.43
10
Corbula gibba
1.00
0.83
0.20
0.40
7.00
20.60
5.00
1.50
0.39
11
12.00
132
Parvicardium minimum
0.67
0.17
0.20
1.00
1.00
1.00
0.14
4
1.00
4
Macoma calcarea
0.67
0.40
0.40
0.50
0.60
7.00
6.00
2.00
7.00
8.67
0.39
11
6.36
70
Astarte elliptica
0.67
0.20
0.60
0.20
3.00
2.00
3.33
1.00
0.25
7
2.71
19
Thyasira gouldi
0.33
0.20
0.40
0.80
1.00
1.00
3.00
2.00
0.32
9
1.89
17
Nuculana pernula
0.33
0.25
1.00
2.00
0.11
3
1.33
4
Buccinum undatum
0.40
0.40
0.25
0.40
1.00
1.00
1.00
1.00
0.25
7
1.00
7
Echinodermata
Strongylocentrotus droebachiensis
0.67
2.00
0.07
2
2.00
4
Amphiura chiajei
0.67
0.83
0.40
1.50
3.60
3.50
0.32
9
3.11
28
Ophiocten affinis
1.00
1.00
0.60
0.60
0.50
0.20
13.33
6.17
4.33
3.67
1.00
2.00
0.64
18
5.83
105
Ophiura albida
1.00
0.50
0.40
0.40
0.50
0.40
2.00
1.33
1.00
1.00
2.50
2.00
0.50
14
1.64
23
Amphiura filiformis
0.67
0.50
0.60
0.40
0.25
0.40
1.50
5.67
3.00
1.00
7.00
2.00
0.46
13
3.23
42
Ophiura robusta
0.67
0.33
0.20
2.00
1.00
1.00
0.18
5
1.40
7
Amphipholis squamata
0.33
0.67
0.40
0.20
0.20
2.00
1.25
1.00
4.00
2.00
0.32
9
1.67
15
Brissopsis lyrifera
0.67
0.40
0.20
1.00
1.00
1.00
0.18
5
1.00
5
Echinocardium flavescens
0.33
0.40
0.40
0.25
0.60
1.00
1.50
1.00
1.00
1.67
0.36
10
1.30
13
Polychaeta
Euchone rubrocinta
0.33
1.00
0.04
1
1.00
1
Exogone verugera
0.33
1.00
0.04
1
1.00
1
Harmothoe cf. fragilis
0.33
2.00
0.04
1
2.00
2
Laonice cirrata
0.33
2.00
0.04
1
2.00
2
Nicolea zostericola
0.33
1.00
0.04
1
1.00
1
Pherusa plumosa
0.33
2.00
0.04
1
2.00
2
Pholoe pallida
0.33
1.00
0.04
1
1.00
1
Praxillura longissima
0.33
1.00
0.04
1
1.00
1
Rhoine gracilior
0.33
2.00
0.04
1
2.00
2
Gyptis rosea
0.67
0.20
1.00
1.00
0.11
3
1.00
3
Eteone cf. longa
0.67
0.25
1.00
1.00
0.11
3
1.00
3
Galathowenia oculata
0.17
1.00
0.04
1
1.00
1
Laetmonice filicornis
0.17
1.00
0.04
1
1.00
1
Orbinia norvegica
0.17
1.00
0.04
1
1.00
1
Typosyllis armillaris
0.17
1.00
0.04
1
1.00
1
Parougia eliasoni
0.20
0.25
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Flabelligera affinis
0.20
1.00
0.04
1
1.00
1
Paraamphinome jeffreysii
0.25
0.20
2.00
1.00
0.07
2
1.50
3
Amage auricola
0.25
2.00
0.04
1
2.00
2
Mediomastus fragilis
0.25
1.00
0.04
1
1.00
1
Neoamphitrite grayi
0.20
1.00
0.04
1
1.00
1
Mollusca
Parvicardium pinnulatum
0.33
1.00
0.04
1
1.00
1
Tellimya ferruginosa
0.67
0.17
1.00
1.00
0.11
3
1.00
3
Aporrhais pespelecani
0.33
0.17
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Cuspidaria cuspidata
0.17
1.00
0.04
1
1.00
1
Kurtiella bidentata
0.17
1.00
0.04
1
1.00
1
Abra alba
0.17
0.20
1.00
1.00
0.07
2
1.00
2
Astarte sulcata
0.20
0.20
0.20
1.00
1.00
1.00
0.11
3
1.00
3
Antalis entalis
0.20
1.00
0.04
1
1.00
1
Diaphana minuta
0.20
1.00
0.04
1
1.00
1
Cylichna alba
0.25
1.00
0.04
1
1.00
1
Lunatia montagui
0.20
1.00
0.04
1
1.00
1
Mytilus edulis
0.20
2.00
0.04
1
2.00
2
Rissoidae
0.20
2.00
0.04
1
2.00
2
Echinodermata
Echinocardium cordatum
0.17
0.25
2.00
1.00
0.07
2
1.50
3
Ophiura ophiura
0.17
0.25
1.00
1.00
0.07
2
1.00
2
Ophiopholis aculeata
0.20
2.00
0.04
1
2.00
2
Ophiothrix fragilis
0.20
1.00
0.04
1
1.00
1
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 (R1) and 12.9 (Rmax⁡; see Table 3). The Procrustes analysis revealed high configurational concordance between parallel DCA and GNMDS ordinations: Procrustes correlation coefficients r 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 (R1) representation of the species-by-site matrix, the second axis was not confirmed according to the “rule-of-thumb” of τ≥0.4 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 (Figure3) exhibited a dense clustering of sites at the intermediate abundance range (R32), while, for both the qualitative data set (R1) and for data with strong weight given to abundance (Rmax⁡), 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 R was observed. Gradient lengths of DCA1 increased with increasing weight given to abundance (1.97 SD units for R=1 to 3.00 SD units for R=max⁡). 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.Figure 3
DCA site diagram for axes 1 and 2 for three abundance ranges with no (R1), intermediate (R32), and strong (Rmax⁡) weight given to abundance. Positions of sampled sites 1–28, with affiliation to broad-scaled geomorphological feature indicated by the color of the dot, are shown. DCA axes are scaled in SD units of compositional turnover.


## 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 (Figure4; 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 (R1), for example, BPI (DCA1) and sediment sorting (DCA2), while other variables were most strongly correlated with DCA axes obtained with intermediate weighting of abundance (R32), 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.Table 5
Kendall'sτ correlations between environmental variables and DCA axes 1 and 2 with no (R1), intermediate (R32), and strong (Rmax⁡) weight given to abundance.
Variable/range
DCA1
DCA2
R
1 (qual)
R
32
R
max
⁡ (raw)
R
1 (qual)
R
32
R
max
⁡ (raw)
DEM derived
Aspect_E
0.011
−0.026
−0.090
0.148
0.021
0.196
Aspect_N
0.101
0.032
−0.074
−0.069
−0.143
−0.085
BPI
0.228
0.169
0.148
−0.016
−0.164
−0.201
Depth
−0.656
−0.820
−0.757
−0.053
−0.063
0.079
Seabed_curv
0.169
0.132
0.079
−0.127
−0.169
−0.312
Slope
0.127
−0.005
−0.153
−0.148
−0.021
−0.132
Current model
Current_bottom
−0.026
−0.095
−0.085
−0.037
−0.143
0.085
Current_surface
−0.307
−0.407
−0.471
−0.212
−0.106
−0.090
Change_curr_bottom
0.074
0.101
0.132
−0.032
−0.053
−0.005
Change_curr_surface
−0.058
−0.138
−0.116
0.090
0.143
0.254
Direction_curr_E
0.024
−0.024
−0.083
−0.008
−0.110
−0.143
Direction_curr_N
0.179
0.094
0.083
−0.035
−0.265
−0.099
Sediment features
Mud content
0.053
0.122
0.069
−0.116
0.095
0.164
Sed_kurt
0.240
0.272
0.245
−0.191
−0.046
−0.073
Sed_med
−0.027
−0.074
−0.095
0.090
0.000
−0.032
Sed_skew
0.099
0.184
0.152
−0.056
0.088
0.318
Sed_sort
0.130
0.135
0.114
0.227
0.146
0.205
Other variables
BMD
−0.412
−0.434
−0.445
0.233
0.228
0.385
Penetration_depth
−0.174
−0.227
−0.286
0.184
0.067
0.206
CN
0.265
0.286
0.318
−0.170
−0.202
−0.260
TOC
−0.202
−0.303
−0.330
−0.245
−0.250
−0.106
TN
−0.257
−0.349
−0.393
−0.149
−0.149
0.008
VEAS_dist
0.349
0.407
0.407
0.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 τ=0 when n=28. The unconfirmed second axis (DCA2.R1) according to criteria given in Section 2.6 is shown in italic.Figure 4
DCA diagrams for axes 1 and 2 with vector fitted environmental variables related to the three selected abundance ranges (R1, R32, and Rmax⁡). Variables were fitted to the DCA axes by application of the envfit function with selection criteria |τ| ≥ 0.225. The length of the vectors indicates their degree of correlation. Sites are shown as black dots. Environmental variables are described in Table 2. DCA axes are scaled in SD units of compositional turnover.The DCA ordination with intermediate weighting of abundance (R32) 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).Table 6
Kendall’s rank correlation coefficientsτ between environmental variables and site scores along DCA axis 1, calculated separately for subsets of sites near the positive (Subset+; n=8) and negative ends (Subset−; n=22; sites 4 and 6 are included in both subsets) of the axis. Correlation coefficients τ > 0.6 are shown in bold. See Table 2 for description of environmental variables.
Var
Subset+
Subset−
τ
P
τ
P
DEM derived
Aspect_E
–0.21
0.548
–0.01
0.956
Aspect_N
0.14
0.720
–0.04
0.824
BPI
0.07
0.905
0.19
0.239
Depth
–0.93
0.000
–0.71
0.000
Seabed_curv
0.07
0.905
0.14
0.371
Slope
0.07
0.905
–0.02
0.911
Current model
Current_bottom
0.14
0.720
–0.06
0.696
Current_surface
–0.50
0.109
–0.19
0.239
Change_curr_bottom
–0.07
0.905
0.13
0.434
Change_curr_surface
–0.43
0.179
–0.22
0.160
Direction_curr_E
0.11
0.708
–0.11
0.479
Direction_curr_N
–0.14
0.720
0.21
0.166
Sediment  features
Mud content
–0.64
0.031
0.28
0.071
Sed_kurt
–0.43
0.179
0.30
0.054
Sed_med
0.64
0.031
–0.20
0.185
Sed.skew
–0.55
0.061
0.38
0.014
Sed.sort
0.11
0.708
0.14
0.381
Other variables
BMD
–0.50
0.109
–0.33
0.036
Penetration_depth
–0.43
0.179
0.00
1.000
CN
0.64
0.031
0.19
0.215
TOC
–0.43
0.179
–0.17
0.283
TN
–0.55
0.061
–0.25
0.113
VEAS_dist
0.50
0.109
0.26
0.091Table 7
Kendall’s rank correlations between selected important environmental variables, calculated separately for subsets of sites near the positive (Subset+;n=8; upper triangle) and negative ends (Subset−; n=22; lower triangle) of the main coenocline (DCA1). Sites 4 and 6 are included in both subsets. Correlation coefficients between variables strongly correlated with depth in Subset+ are shown in bold. See Table 2 for description of environmental variables.
Depth
BMD
Penetration depth
CN
Sed_med
Mud content
TN
TOC
Depth
*****
0.43
0.36
−0.57
−0.71
0.71
0.47
0.36
BMD
0.28
*****
0.93
−0.29
−0.29
0.14
0.11
0.07
Penetration_depth
−0.16
0.38
*****
−0.21
−0.36
0.21
0.04
0.00
CN
−0.17
−0.55
−0.12
*****
0.43
−0.57
−0.33
−0.21
Sed_med
0.14
0.38
0.19
−0.26
*****
−0.86
−0.26
−0.21
Mud content
−0.26
−0.30
−0.01
0.11
−0.68
*****
0.33
0.21
TN
0.31
0.28
−0.12
−0.38
0.35
−0.44
*****
0.91
TOC
0.23
0.04
−0.12
0.09
0.25
−0.41
0.56
*****Table 8
Matrix of Kendall’s correlation coefficientsτ between pairs of environmental variables (upper right triangle) and corresponding  P values for tests of the hypotheses that τ = 0 (lower left triangle). P values ≤ 0.10 and the corresponding τ values are shown in bold.
BPI
Depth
Seabed_curv
Current_bottom
Direction_curr_E
Direction_curr_N
Current_surface
Change_curr_bottom
Sed_kurt
Sed_med
Mud content
Med_skew
Med_sort
BMD
Penetration_depth
TOC
TON
CN
VEAS_dist
BPI
*****
-
0.222
0.646
0.079
-
0.164
0.185
-
0.127
0.201
0.164
-
0.133
0.190
0.163
-
0.032
-
0.141
0.077
−0.223
-
0.100
-
0.138
0.106
Depth
0.101
*****
-
0.175
0.074
0.051
-
0.067
0.481
-
0.090
-
0.213
0.032
-
0.111
-
0.184
-
0.151
0.407
0.131
0.340
0.387
−0.271
−0.439
Seabed_curv
<0.001
0.201
*****
0.053
0.078
0.319
-
0.153
0.090
0.223
-
0.154
0.143
0.019
-
0.211
-
0.222
-
0.051
-
0.069
0.014
-
0.048
0.101
Current_bottom
0.570
0.597
0.710
*****
-
0.073
0.233
-
0.201
0.286
-
0.110
0.159
-
0.138
-
0.104
0.032
-
0.092
-
0.029
0.160
0.089
0.005
0.254
Direction_curr_E
0.227
0.707
0.566
0.593
*****
0.319
0.089
-
0.116
0.077
0.003
-
0.121
−0.234
−0.239
-
0.041
0.027
0.349
0.187
0.256
0.024
Direction_curr_N
0.172
0.621
0.018
0.085
0.020
*****
-
0.120
0.297
0.316
−0.244
0.174
0.095
−0.252
-
0.129
0.022
0.089
0.060
0.078
0.142
Current_surface
0.356
<0.001
0.263
0.139
0.513
0.373
*****
-
0.090
-
0.105
0.011
-
0.079
-
0.136
-
0.070
0.206
0.045
0.383
0.338
-
0.080
−0.524
Change_curr_bottom
0.139
0.518
0.518
0.034
0.394
0.028
0.518
*****
0.202
−0.265
0.196
0.158
−0.232
-
0.173
-
0.045
-
0.101
-
0.149
0.111
0.196
Sed_kurt
0.227
0.117
0.100
0.416
0.578
0.021
0.440
0.137
*****
−0.567
0.568
0.457
−0.256
−0.342
-
0.152
-
0.043
-
0.061
0.067
0.073
Sed_med
0.323
0.813
0.252
0.236
0.984
0.072
0.937
0.048
<0.001
*****
−0.796
−0.522
0.450
0.114
-
0.040
0.176
0.209
-
0.005
-
0.048
Mudcontent
0.162
0.421
0.298
0.317
0.373
0.198
0.570
0.150
<0.001
<0.001
*****
0.724
−0.249
-
0.076
0.136
−0.277
−0.257
-
0.095
0.005
Sed_skew
0.227
0.172
0.890
0.440
0.088
0.488
0.313
0.243
0.001
<0.001
<0.001
*****
0.041
-
0.014
0.248
−0.352
−0.268
-
0.158
-
0.040
Sed_sort
0.812
0.267
0.122
0.812
0.084
0.067
0.606
0.088
0.064
0.001
0.068
0.766
*****
0.122
0.136
-
0.125
-
0.036
-
0.114
-
0.038
BMD
0.302
0.003
0.104
0.500
0.765
0.350
0.132
0.204
0.014
0.404
0.578
0.921
0.381
*****
0.534
0.090
0.258
−0.514
−0.385
Penetration_depth
0.566
0.332
0.707
0.828
0.843
0.874
0.737
0.737
0.266
0.767
0.313
0.068
0.321
<0.001
*****
-
0.013
-
0.011
−0.228
-
0.152
TOC
0.097
0.011
0.607
0.235
0.011
0.513
0.004
0.452
0.751
0.192
0.040
0.009
0.361
0.512
0.921
*****
0.667
0.000
−0.277
TON
0.463
0.005
0.921
0.513
0.176
0.662
0.013
0.275
0.662
0.126
0.059
0.052
0.796
0.064
0.937
<0.001
*****
−0.361
−0.306
CN
0.304
0.044
0.722
0.968
0.060
0.566
0.553
0.406
0.620
0.968
0.477
0.243
0.405
<0.001
0.092
1.000
0.008
*****
0.207
VEAS_dist
0.445
0.001
0.468
0.060
0.858
0.294
<0.001
0.150
0.593
0.722
0.984
0.767
0.781
0.005
0.259
0.040
0.025
0.123
*****Figure 5
Contour-plots of selected environmental variables on DCA.R32, axes 1 and 2, showing patterns of variation in ordination space. The contour-plots are fitted by use of the ordisurf command fitting GAM with Gaussian errors. The size of the bubbles is proportional to the magnitude of the target value and are comparable between plots because environmental variables are ranged between 0 and 1. Axes are scaled in SD units of compositional turnover. The equidistance between contours differs between plots. See Table 2 for description of environmental variables.Sediment profiles images (SPIs) of sites 7 (a) and 25 (b) showing the difference in sediment particle composition and environment at the two extremes of DCA1.
(a)
(b)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 theThyasira 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.Table 9
Species associated with the negative (Neg) and positive (Pos) extremes of DCA axes 1 and 2, with indication of feeding strategy, species loading on axes (“Load” column), and raw total abundance for each species (“Abund” column). Only species with DCA.R32 species loadings with absolute values >1.25 and with a raw total abundance above 10 were chosen.
Axis
End
Load
Species
Su
SD
SSD
G
Sc
P
Abund
DCA1
Neg
−2.00
Thyasira gouldi(M)
x
x
17
−1.96
Sige fusigera(P)
x
x
13
−1.81
Pista lornensis(P)
x
18
−1.47
Thyasira flexuosa(M)
x
x
206
−1.46
Melinna cristata(P)
x
66
−1.36
Thyasira sarsi (M)
x
x
207
−1.34
Ennucula tenuis(M)
x
2216
−1.26
Spiophanes kroyeri(P)
x
x
2389
Pos
2.78
Lepeta caeca(M)
x
x
46
2.59
Hiatella arctica(M)
x
x
12
2.46
Scoloplos armiger(P)
x
41
1.94
Pseudopolydora pauchibranchiata(P)
x
x
746
1.76
Corbula gibba(M)
x
132
1.75
Glycera lapidum(P)
x
x
53
1.66
Eupolymnia nebulosa(P)
x
15
1.61
Amphiura chiajei(E)
x
x
28
1.58
Pseudamussium peslutrae(M)
x
19
1.44
Ophiocten affinis(E)
x
105
1.41
Melinna elisabethae(P)
x
x
54
1.31
Anobothrus gracilis(P)
x
381
DCA2
Neg
−3.38
Scoloplos armiger (P)
x
41
−1.79
Astarte elliptica(M)
x
19
−1.71
Macoma calcarea(M)
x
x
x
70
−1.65
Aphrodita aculeata(P)
x
24
−1.56
Ophiura albida(E)
x
x
x
23
−1.35
Lepeta caeca(M)
x
x
46
Pos
4.81
Dipolydora coeca(P)
x
x
11
1.77
Amphicteis gunneri(P)
x
x
11
1.65
Amphiura chiajei(E)
x
x
28
1.59
Phyllodoce rosea(P)
x
12
1.54
Streblosoma bairdi(P)
x
x
21
1.38
Corbula gibba(M)
x
132
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, P<0.01). 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)).Contour-plots showing patterns of variation in ordination space for (a) total abundance (total number of individuals of all species), (b) species richness, and (c) species evenness for axes 1 and 2 of DCA.R32. The contour-plots are fitted by use of the ordisurf command fitting GAM with Poisson errors. The sizes of the bubbles are proportional to the magnitude of the target value, but not comparable between plots. Axes are scaled in SD units of compositional turnover. The equidistance between contours differs between plots.
(a)
(b)
(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 termperiodic 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 (Figure6) 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 polychaetesA. gracilis,M. elisabethae, andP. 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 termperiodic 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 [50–53]. 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 (Figure4), 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 currentgradient.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 wherephysical forcing is likely to be the main structuring factor, and communities dominated by subsurface deposit-feeders in deeper waters where periodic hypoxiaand 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 thephysical 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 theperiodic 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 thephysical 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].


## 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 termperiodic 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 (Figure6) 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 polychaetesA. gracilis,M. elisabethae, andP. 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 termperiodic 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 [50–53]. 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 (Figure4), 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 currentgradient.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 wherephysical forcing is likely to be the main structuring factor, and communities dominated by subsurface deposit-feeders in deeper waters where periodic hypoxiaand 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 thephysical 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 theperiodic 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 thephysical 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 ofhypoxia 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.Figure 8
Conceptual model summarising major complex-gradients and auxiliary variables related to the main and secondary coenoclines. Hypothesised ecoclines related to the main coenocline are given in quotes.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.



---

*Source: 462529-2014-02-16.xml*