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 termed Periodic hypoxia and Periodic physical forcing. We conclude that a combination of recurrent (periodical) and extreme events is likely to determine the variation found in the species composition of marine sedimentary ecosystems. Based on the results of our study, we conclude that indirect gradient analysis is a useful tool for enhancement of our basic mechanistic understanding of the processes governing the compositional structure of marine sediment communities.

1. Introduction

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

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

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

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

2. Material and Methods

2.1. Study Area

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

2.2. Bathymetry

In 2004 and 2005, the Geological Survey of Norway performed bathymetric surveys covering the whole inner part of the Oslofjord [20]. Onboard the research vessel (R/V) Seisma, an interferometric sonar system (GeoSwath, Geo-Acoustics) with 125 and 250 kHz transducers was used to provide a three-column xyz-file with georeferenced depth values. The xyz-file had a resolution of 1 m and was imported into GRASS 6.4.1 [21] as a gridded bathymetry raster using the 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 ( 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).

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

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

2.4. Faunal Samples

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

2.5. Environmental Variables

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

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

2.6. Data Analysis

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

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

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

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

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

3. Results

3.1. Species Richness, Abundance Patterns, and Core Species

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

3.2. Faunal Analysis: Correspondence between DCA and GNMDS Ordinations

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

3.3. Description of DCA Ordination Patterns

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

3.4. Relationships between DCA Axes and Recorded Environmental Variables

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

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

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

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

4. Discussion

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

4.1. Ecoclines Related to the Main Coenocline (DCA1)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.3. Mechanisms and Processes That Account for Variation along Ecoclines

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

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

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

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

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

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

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

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

5. Conclusion

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

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

Conflict of Interests

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

Acknowledgments

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