Long-term changes in North Adriatic soft-bottom polychaete assemblages following a dystrophic crisis

Abstract

In the Northern Adriatic Sea, the combination of natural and anthropogenic processes causes periodic hypoxia and anoxia, which can cause mass mortality of benthic organisms. Here, we report the results of a 20-year monitoring programme carried out after an anoxic event in 1989 that caused benthic mass mortality over an area of about 1200 km². The recovery dynamics of the polychaete fauna was followed at three stations (including one sampled a year before the anoxic event). Samples (0.1 m² Van Veen grab, 2-mm sieve) were collected from 1989 to 1994 and from 2003 to 2008. A total of 6519 specimens belonging to 186 species were identified. The analyses of univariate biodiversity indexes highlighted higher temporal variation of assemblage diversity in the period following the dystrophic crisis, followed by higher stability in the next decade. PERMANOVA suggested that the highest component of assemblage variation belonged to the factor year. However, differences across years varied in magnitude, depending on stations and were portrayed by non-metric multidimensional scaling plots of each year's centroids for each of the three stations. SIMPER analysis identified the species most characterizing the assemblages in each year for each of the three stations. Our results highlight the importance of analysing long-term data sets in order to understand assemblage dynamics following strong disturbance events. Assuming the 1989 anoxia as the event determining the biodiversity change, and considering our outcomes, which suggest that the assemblages still show a pattern of non-random changes through years, our data indicate that the polychaete fauna is still recovering without having reached a pattern of among-years stability. Focusing on the biology of the species characterizing the different periods might help in understanding the ecological processes that have determined the observed pattern.

Keywords:

• Long-term biomonitoring
• Northern Adriatic
• soft bottoms
• Polychaeta
• anoxia
• taxonomy

Introduction

Despite increasing concerns among ecologists, policy-makers and the general public about anthropogenic environmental changes in the marine ecosystem, their possible long-term effects on marine biodiversity are still poorly addressed quantitatively. The need to widen the spatial and temporal scales of monitoring programmes is widely recognized, as is the importance of taxonomy in quantifying the actual impact of human activities on marine biodiversity (Terlizzi et al. 2005; Bevilacqua et al. 2009; Musco et al. 2009, 2011; Mikac & Musco 2010). Hypoxia, defined by dissolved oxygen concentrations less than 2 ml/l O₂, and anoxia, defined as the complete absence of dissolved oxygen (0 ml/l O₂) (Diaz et al. 2004), are among the most widespread disastrous anthropogenic impacts on estuarine and marine environments. There is no other variable of such ecological importance to coastal marine ecosystems that has changed so drastically over such a short period of time as dissolved oxygen (Diaz & Rosenberg 2008). Starting in the 1960s, the number of systems reporting hypoxia-related problems has increased, with the main factor leading to these events being the input of excess nutrients leading to eutrophication. Hypoxia can strongly affect compositional and structural properties of assemblages by discriminating sensitive species in favour of tolerant ones and interfering with recruitment and growth. Ultimately, hypoxia can lead to the mass mortality of benthic organisms.

The Northern Adriatic Sea (hereafter NA) (Mediterranean Sea) is a sensitive ecosystem combining many features that promote the development of hypoxia: it is a semi-enclosed basin, with relatively shallow depth, characterized by soft substrates, high riverine input, high primary production, strong water stratification and long water residence time. Events of low dissolved oxygen (DO) and benthic mortalities have been noted here periodically for centuries (Crema et al. 1991), but following increasing anthropogenic input of nutrients, their frequency and intensity have increased markedly from the 1960s (Justić 1987). Since the 1980s, oxygen crises have been detected frequently in various areas of the NA (Stachowitsch 1984, 1991; Degobbis et al. 1995; Malej & Malačić 1995; Kollmann & Stachowitsch 2001). Such repeated events have caused destabilization and long-term changes in benthic communities. The NA is therefore a case study for long-term decrease in DO concentration and associated benthic community changes and mortalities (Kollman & Stachowitsch 2001; Stachowitsch et al. 2007).

While mass mortality caused by oxygen depletion proceeds very rapidly (Stachowitsch 1991), assemblage recovery can be a very slow process that varies considerably, ranging from years to decades (Rosenberg 1976; Stachowitsch & Fuchs 1995; Kollman & Stachowitsch 2001; Wu 2002). Long-term studies are thus necessary to detect the impact of such disturbances on benthic assemblages and to determine how various acute perturbations can influence the overall system (Stachowitsch & Fuchs 1995; Dauvin, 1998, 2000). A limit for interpreting long-term changes in the marine benthic assemblages is the lack of long-term monitoring studies based on sound design and fine taxonomy (i.e. identification at species level).

Considering their abundance, species richness and functional diversity, polychaetes are among the most important macrobenthic groups in soft bottoms (Knox 1977). They show remarkable levels of adaptation to a wide range of habitats and environmental conditions and play key roles in the functioning of benthic communities. Polychaetes might be considered as surrogates for the estimation of diversity and dynamics of the whole benthic assemblages (Olsgard & Somerfield 2000; Olsgard et al. 2003; Giangrande et al. 2005).

In autumn 1989, severe bottom oxygen depletion caused benthic mass mortality over an area of about 1200 km² in the offshore NA (Zavodnik et al. 1994). Here we report the results of the 20-year monitoring programme following that event, with the aim of analysing patterns of change in polychaete assemblages over the long term.

Materials and methods

Research area

The study was carried out in the Northern Adriatic Sea on three offshore monitoring stations of the Centre for Marine Research (Ruđer Bošković Institute, Rovinj). Station SJ 005 (hereafter 005) (45°18.4' N; 13°18.0' E; 31 m depth) and station SJ 007 (007) (45°17.0' N; 13°16.0' E; 31 m depth) are situated on the profile Poreč (Croatia)–Venice-Lido (Italy). Station SJ 107 (107) (45°02.8' N; 13°19.0' E; 37 m depth) is situated on the profile Rovinj (Croatia)–river Po delta (Italy) (Figure 1). Long-term abiotic data in terms of high-frequency estimations of oceanographic parameters are also available for this latter station. Granulometric analyses showed high similarity in sediment composition among the three stations. The sediment is weakly sorted silty sand composed mainly of sand (ranging 66.3–70.6%, depending on the station), followed by silt and clay fractions (27.8–32.1%) and gravel (1.4–1.6%).

Figure 1. Map of the research area and studied stations. The lined area is the one affected by anoxia in 1989.

Field and laboratory work

Water samples for the estimation of the oceanographic parameters were collected approximately monthly during the period 1986–2008 at station 107 at the layer 2 m above the bottom with 5-litre Niskin bottles. The DO concentration was determined by the Winkler titration method. Standard methods (Parsons et al. 1984) were used for the measurement of nutrients (total inorganic nitrogen, orthophosphate and orthosilicate), temperature, density, salinity, pH and chlorophyll a.

Macrofaunal samples were undertaken during two main periods, 1988–1994 and 2003–2008. Sampling at station 107 also included the pre-anoxic period (June and October 1988 and August 1989), while sampling at the other two stations started after the anoxia event, from December 1989. After the anoxic event (from 1990 to 1994 and from 2003 to 2008), samples were taken from late October to the beginning of January, once a year, at all stations.

At each station, 4–5 replicates were collected with 0.1 m² Van Veen grab. In the first research period (1988–1994) samples were sieved through a 2-mm mesh, while in the second they were sieved through 1- and 2-mm meshes. Samples were fixed in 4% buffered formaldehyde–seawater solution. The use of 1-mm mesh size in the second period was performed according the up-to-date recommendation in the macrobenthic fauna monitoring studies (i.e. Castelli et al. 2004). In this study, only 2-mm mesh size samples were used to analyse the whole research period, in order to allow long-term comparisons. However, 1-mm mesh size data from the recent most period were used in order to test the representativeness of the 2-mm mesh size samples.

In the laboratory the material was rinsed and macrobenthic organisms sorted and preserved in 70% ethanol. Polychaetes were determined to the species level using stereo and light microscopes. Nomenclature followed Castelli et al. (2008).

Data analyses

A descriptive, univariate analysis was performed on mean number of polychaete species, mean abundance, diversity indices and mean abundances of dominant species by including data about the pre-anoxic period from station 107 (June and December 1988 and August 1989) and, for all three stations, data from the anoxic event in December 1989 and post-anoxic period from 1990 to 2008.

A distance-based permutational multivariate analysis of variance (PERMANOVA; Anderson 2001; McArdle & Anderson 2001) was performed to test the differences in the polychaetes assemblages' structure among years and stations. PERMANOVA was based on Bray–Curtis similarity matrix of untransformed data. The experimental design included two factors: Year (Y, fixed, 10 levels) and Station (S, random, with 3 levels and crossed to Y). The quantitative multivariate description of post anoxic long-term changes in polychaete assemblages included sampling times from 1990 to 2008.

In order to assess the potential influence of mesh size in the characterization of multivariate structure of the polychaete assemblages, a Spearman's correlation test was done between the Bray–Curtis similarity matrices for the polychaete assemblages' data in the period 2003–2005 from 2-mm mesh size and data including both 1- and 2-mm mesh sizes.

As PERMANOVA highlighted the significance of the interaction term Y×S, thus indicating that patterns of temporal changes varied across stations (see results), three separate non-metric multidimensional scaling (nMDS) ordination plots of the year's centroids were produced to visualize differences across years at each of the three stations. Centroids were obtained using principal coordinates (see Terlizzi et al. 2005).

The SIMPER procedure (Clarke 1993) was employed to select, for each of the three stations, the species most characterizing the assemblages in each of the considered years. For these species feeding guilds analysis was done following Fauchald and Jumars (1979) and Gambi and Giangrande (1985).

A distance-based multivariate analysis for a linear model using forward selection, (DISTLM-forward; Anderson 2003), was performed to investigate possible correlations between the spatial–temporal distribution of the polychaetes and the environmental predictors variables at station 107. A separate analysis was performed in order to investigate individually the correlation between the assemblage distribution and the lowest oxygen values recorded every year at station 107 before the macrofaunal sampling.

All multivariate analyses were performed using the computer program PRIMER v6 (Clarke & Gorley 2006), including the add-on package PERMANOVA + (Anderson et al. 2008).

Results

Descriptive analyses

Anoxia was recorded in 1989 with a minimal oxygen value of 0.23 ml/l. Before the anoxic event, hypoxia occurred in September and November 1988 and in October 1989 (Figure 2). After the anoxic event, hypoxia was recorded only in September and October 1991 and never evidenced afterwards until the end of the research period. Considering the whole study period, a total of 6519 specimens belonging to 186 species were collected. Overall, the most abundant species were Owenia fusiformis (17.1%), followed by Pectinaria auricoma (11.7%), Notomastus latericeus (7.8%), Eunice vittata (6.6%) and Nothria conchylega (6.2%).

Figure 2. Oxygen concentrations measured in the bottom layer water on station 107 in the period 1986–2008.

The anoxic event corresponded to an abrupt decrease in the average number of species and individuals (Figure 3a,b). In December 1989, at station 005 only 4 species, namely E. vittata (3 ind.), Malmgreniella andreapolis (3 ind.), Marphysa bellii (1 ind.) and O. fusiformis (41 ind.) survived the anoxia. At station 007, only 4 specimens of E. vittata were found, while at station 107, 7 species, namely Drilonereis filum (1 ind.), E. vittata (1 ind.), Galathowenia oculata (12 ind.), Lumbrineris sp. (1 ind.), Neanthes succinea (2 ind.), Nephthys hystricis (4 ind.) and Nephthys incisa (1 ind.), were found. The average number of species grew rapidly one year after the anoxia and after one to four years (depending on the station) reached the levels comparable to those characterizing the 2003–2008 period (Figure 3a). In both research periods (1989–1994 and 2003–2008), the species richness appeared quite variable from year to year. At all stations there was an abrupt increase in the average polychaete abundance during 1990 and 1991 followed by a rapid decline in 1992 (Figure 3b). At station 107, the maximum abundance was reached in 1990, while at stations 005 and 007 it was reached in 1991. Four years after the anoxia (in January 1994) the average number of individuals reached the levels characterizing the 2003–2008 period. The Shannon–Wiener diversity index (H') decreased rapidly at station 107 during the anoxia and was very low in December 1989 at all stations (Figure 3c). The average value of H' increased in 1990 and three to four years after the anoxia it reached the levels characterizing the 2003–2008 period.

Figure 3. Temporal variation of the univariate indexes of polychaete assemblages at three research stations. (a) Number of species; (b) Number of individuals; (c) Shannon–Wiener diversity index (H').

Overall patterns of mean abundances of the 12 dominant polychaete species are shown in Figure 4. Among them only Ampharete acutifrons and Sternaspis scutata were not present in the polychaete assemblage on station SJ 107 before the anoxic event, while, taking into consideration all three stations, only O. fusiformis, G. oculata and E. vittata survived the anoxia. Species such as N. latericeus, A. acutifrons and S. scutata became dominant only in the most recent period. Other species such as Amphicteis gunneri showed the opposite trend of higher abundance in the pre-anoxic period and the first period following the anoxia, and a decrease in abundance in the most recent period. Most of the species (particularly O. fusiformis, P. auricoma, E. vittata and N. conchylega, and to a lesser extent also G. oculata and Aponuphis brementi) showed rapid population growth one or two years after the anoxia, followed by a rapid decline in the third year. Those species mostly influenced the observed rapid growth followed by the rapid decline of the whole polychaete fauna.

Figure 4. Trends of mean abundances of the 12 most abundant polychaete species overall. For each species, the total number of individuals is given in brackets.

Multivariate analyses

The Spearman correlation test between polychaete assemblages from both 1- and 2-mm and only 2-mm mesh size showed high level of significant correlation (Rho = 0.79; significance level of sample statistic: 0.1%).

PERMANOVA showed the significance of the interaction term Y×S, suggesting that the multivariate structure of the polychaete assemblages changed across years with a pattern varying across stations (Table I). The analysis also suggested that the highest component of variation belonged to the Year factor. How the patterns of assemblage variation through years changes inconsistently with stations is well represented in the nMDS plots in Figure 5, where interannual trajectories were over-imposed to show the overall pattern of change in assemblages through time. Different temporal patterns among the stations were also confirmed by the results of the SIMPER analyses (Table II). The first two years following anoxia (1990, 1991) were characterized by a high level of similarity among replicates at all three stations, contributed by very few species. Those early recovery assemblages were characterized by the tubiculous species O. fusiformis (deposit/filter-feeder) and P. auricoma (deposit-feeder). In 1991 the importance of motile species such as the carnivorous E. vittata (at stations 007 and 107) and the omnivorous N. conchylega (at station 107) became apparent. Beside those four species, the end of the first research period was also characterized by the tubiculous deposit-feeder A. gunneri at station 107. During the whole research period O. fusiformis was always among the species mostly characterizing station 007, while for the other two stations no species constantly characterizing the assemblages were observed. The most recent period (2003–2008) was marked by a different set of species characterizing the polychaete assemblages. In 2003 the importance of the motile burrowing deposit-feeder N. latericeus was evident at all stations. Although declining with time, this species remained among those characterizing station 107 until the end of the research period, while its importance rapidly declined at station 005, and, particularly, at station 007. In 2004 the contribution of tubiculous deposit-feeders, namely G. oculata (at stations 007 and 107) and A. grubii (at station 007), was evident. By 2005 to the end of the research, A. acutifrons was a further tubiculous deposit-feeder, particularly characterizing the assemblage at the station 005. This latter species also characterized station 107 in 2006 and station 007 in 2007. In the last three years, new deposit-feeder species, such as the motile Sternaspis scutata, Glycera unicornis and P. fauchaldi as well as the tubiculous S. kroyeri and T. stroemi, became important in characterizing the assemblages. It is worth noting that during the first research period (1990–1994) only 5 species were contributing to 50% similarity among the replicates at all 3 stations, while the recent most period counted 12 such species.

Figure 5. nMDS plots of Year's centroids for each of the three stations.

Table I. Results of PERMANOVA test for the model including 1990–2008. CV, estimates of components of variation

Source	df	SS	MS	F	P(perm)	CV
Year = Y	9	1.3862E5	15403	4.2148	− 0.001	911.87
Station = S	2	25600	12800	9.567	0.001	269.36
Y×S	18	65779	3654.4	2.7314	− 0.001	539.39
Residual	99	1.3246E5	1337.9			1337.9
Analyses based on Bray–Curtis similarity. No transformation. Number of perm. used: 999. Permutation of residuals under a reduced model.

Table II. Results of the SIMPER analysis showing the year-by-year average Bray Curtis similarity (Av. S.) among replicates of each station, and the relative contribution of the species (Sp. Contr.%). Cut off for low contributions: 50%

	107		005		007
	Av. S.	Sp. Contr.%	Av. S.	Sp. Contr.%	Av. S.	Sp. Contr.%
1990	68.39	P. auricoma 51.70	61.59	O. fusiformis 55.41	66.28	O. fusiformis 48.12
						P. auricoma 17.98
1991	61.38	N. conchylega 30.33	60.99	O. fusiformis 82.16	63.66	E. vittata 21.10
		E. vittata 22.71				O. fusiformis 20.16
						P. auricoma 16.02
1992	47.12	E. vittata 36.45	49.30	O. fusiformis 46.86	47.73	O. fusiformis 38.31
		N. conchylega 16.13		P. auricoma 18.42		E. vittata 13.93
1994	47.94	P. auricoma 30.08	55.35	P. auricoma 53.28	39.40	O. fusiformis 23.80
		O. fusiformis 11.00				N. conchylega 15.09
		A. gunneri 10.95				P. auricoma 12.56
2003	47.38	N. latericeus 71.50	59.30	N. latericeus 59.88	40.87	N. latericeus 40.58
						O. fusiformis 15.07
2004	37.20	N. latericeus 40.47	55.55	N. latericeus 34.79	45.71	O. fusiformis 31.37
		G. oculata 12.33		O. fusiformis 17.27		G. oculata 14.07
						A. grubii 10.13
2005	38.12	N. latericeus 34.44	47.33	N. latericeus 24.07	44.24	O. fusiformis 30.99
		P. auricoma 16.06		P. auricoma 23.23		G. oculata 14.74
				A. acutifrons 16.66		H. gracilis 7.20
2006	46.43	N. latericeus 30.43	56.46	A. acutifrons 34.71	48.77	O. fusiformis 34.53
		A. acutifrons 14.66		G. oculata 9.95		G. oculata 20.63
		S. scutata 11.89		O. fusiformis 9.45
2007	34.81	S. kroyeri 26.36	56.31	A. acutifrons 40.04	32.47	A. acutifrons 19.00
		N. latericeus 16.57		S. scutata 18.51		P. auricoma 17.41
		P. fauchaldi 11.13				O. fusiformis 15.93
2008	44.46	S. kroyeri 21.49	50.63	A. acutifrons 22.05	50.34	P. auricoma 38.58
		P. auricoma 15.85		P. auricoma 16.38		O. fusiformis 8.60
		N. latericeus 11.64		S. kroyeri 14.94		G. unicornis 7.51
		T. stroemi 11.49

Table III reports the results of DISTLM-forward analysis conducted between the abundance matrix of polychaetes at station 107 and environmental predictor variables. The sequential tests (i.e. fitting each variable one at a time, conditional on the variables that are already included in the model) revealed that all the considered variables were significantly correlated to the assemblage spatial–temporal distribution, altogether explaining 60.7% of the assemblage variation. The analysis ranked salinity first (explaining 17.77% of the total variation). The relative contribution of each of the remaining variables were due to pH (8.11%), SiO₄ (7.89%), PO₄ (6.88%), temperature (4.94%), chlorophyll a (4.53%), oxygen (3.96%), density (3.66%) and total inorganic nitrogen (2.98%). A separate DISTLM-forward analysis considering, as a predictor variable, only the lowest oxygen observed in each year before the macrofaunal sampling (oxygen values recorded from September to November depending on the year) revealed a significant (p < 0.01) correlation, with the polychaete spatial–temporal distribution explaining 10.35% of the assemblage variation.

Table III. Results of the DISTLM-forward analysis: Sequential tests. SS (trace), portion of sum of squares on the predictor variable analysed; pseudo-F, statistical; Prop., proportion of variability explained by the predictor variable; Cumul., sum of the variable individual contributions; O₂, oxygen (ml/l); temp, temperature (°C); sal, salinity (psu); dens, density (kg/m³); PO₄, orthophosphate (μmol/l); TIN, Total inorganic nitrogen (μmol/l); SiO₄, orthosilicate (μmol/l); chl a, chlorophyll a (μg/l)

Variable	SS (trace)	Pseudo-F	P	Prop.	Cumul.
sal	22424	8.8487	0.001	0.17751	0.17751
pH	10242	4.3741	0.001	8.1076E-2	0.25859
PO4	8695.9	3.9916	0.001	6.8837E-2	0.32742
SiO4	9972.7	5.0534	0.001	7.8943E-2	0.40637
O2	4999.4	2.6428	0.001	3.9575E-2	0.44594
temp	6236.2	3.5213	0.001	4.9366E-2	0.49531
density	4619.3	2.7339	0.001	3.6567E-2	0.53187
TIN	3766.9	2.3131	0.003	2.9818E-2	0.56169
chl a	5722.9	3.804	0.001	4.5303E-2	0.60700
Correlation among variables
	O2	temp	sal	dens	PO4	TIN	SiO4	pH	chl a
O2
temp	−0.576
sal	0.427	−0.137
dens	0.673	−0.895	0.562
PO4	−0.970	0.430	−0.499	−0.585
TIN	−0.767	0.286	−0.470	−0.458	0.850
SiO4	−0.987	0.495	−0.459	−0.622	0.987	0.842
pH	0.639	−0.447	0.382	0.552	−0.660	−0.747	−0.642
chl a	0.638	−0.809	0.498	0.893	−0.566	−0.416	−0.585	0.345

Discussion

Natural communities are characterized by stochastic changes of the structure of assemblages (Boero 1994; Benedetti-Cecchi 2009). Our results highlighted a large degree of yearly variation of the polychaete assemblages in the research area through the observed period. Differences through years varied with stations but, at all stations, visual inspection of nMDS plots did not suggest a pattern of random variations of assemblages through years but, rather, a global pattern of the interannual trajectories through the whole study period. Our results, therefore, could suggest that, despite spatial inconsistencies in the magnitude of assemblage variation through years, the polychaete assemblages continue to change after the anoxic event by not following a pattern of random through-years variation. The use of a 2-mm mesh size may have reduced the possibility of properly characterizing the polychaete assemblages, as some small opportunistic species may be overlooked. However, the Spearman correlation test showed high correlation between the polychaete assemblages' structure of the whole polychaete macrofauna (1- and 2-mm mesh size) and the subset of data (2-mm mesh size) in the most recent period. This reveals that, at least for the most recent period, our data subset may well represent the whole polychaete macrofauna assemblage patterns. Nonetheless, care must be taken when considering the actual changes in the biodiversity patterns.

A relevant component of this study is that it extends through more than two decades and that it is based on fine taxonomy (i.e. identification at species level), which is not common (Fraser et al. 2006; Schirosi et al. 2010). Long-term studies are needed to obtain reliable estimates of both natural variations in undisturbed systems as well to access the long-term effects of anthropogenic disturbance on complex systems (Rhoads & Germano 1986; Ibanez & Dauvin 1988). Understanding the pattern in biotic assemblages is important before sensible statements can be made about plausible causal processes. Since describing assemblages at species level is logistically difficult, requiring the counting of all organisms and their taxonomic identification, at present taxonomic sufficiency (TS) (Ellis 1985) is recommended as an ideal short-cut procedure in community perturbation studies (Ferraro & Cole 1995; Olsgard et al. 1997; Dauvin et al. 2010). According to Bertrand et al. (2006), TS means ‘identifying organisms only to a level of taxonomic resolution sufficient to satisfy the objective of a study’. However, due to the limited number of case studies, an array of possible biases (i.e. the relationships among TS, spatial scale, habitat features and assemblage structure) including the relationship between TS and temporal scale call for the further investigation to make generalizations (Terlizzi et al. 2003). In a recent study analysing the same polychaete assemblages (Musco et al. 2011), it has been demonstrated that using TS could lead to a hazardous loss of information and it is suggested that periodical analysis at fine taxonomic level should be alternated routinely with long-term monitoring based on TS. Our analysis of temporal pattern at a fine level of taxonomic resolution within a diversified group of benthic invertebrates might be of interest for the study of natural variability and polychaetes biodiversity both at a structural and a functional level. Such studies would be difficult using surrogates of the species (Musco et al. 2009).

The analyses highlighted the first period after the anoxia as characterized by tubiculous deposit-feeders (P. auricoma and O. fusiformis) possibly because of the large quantities of the decomposing organic matter lying in and on the sediment after the 1989 mass mortality. It could be hypothesized that the biomass derived from dead organisms was not immediately consumed by scavengers, which disappeared from the area with the majority of the macrofauna (Jaklin 2002). Beside the deposit-feeders, the following years (1992–1994) were also characterized by the motile carnivorous (E. vittata) and omnivorous (N. conchylega) species that might be interpreted as an increase of the functional complexity of the polychaete assemblages. The year 2003, in particular, as well as the years 2004–2005, were characterized by the burrowing deposit-feeder N. latericeus, whose importance declined with time. The second research period showed a renewed importance of deposit-feeders that were represented by a large number of species. This last condition might be interpreted as an increase of functional redundancy in the assemblages possibly implying increased resilience of the system (Boero 1994). In fact, instability and high variability characterized the very first post-anoxic period due to the abrupt increase and decline of abundance of just few species.

The post-anoxic period was characterized by the rapid growth and decline of the polychaete abundance. This trend can be explained by rapid population growth of a small number of opportunistic species, such as E. vittata, N. conchylega, O. fusiformis and P. auricoma. Afterwards, the population of these species declined rapidly, possibly as a result of competition that started with the arrival of other species, thus evidencing a kind of ecological succession. In fact, while the polychaete abundances declined after the initial bloom, the number of species generally continued to grow by the end of the first research period.

The analysis of dominant species highlighted that some acted as more tolerant (i.e. survived the anoxia, G. oculata), some as opportunistic showing some kind of population ‘flush and crash’ (N. conchylega and P. auricoma), while others showed both behaviours (O. fusiformis and E. vittata). Other species appeared to be sensitive, becoming abundant only in the most recent period (A. acutifrons). Simboura and Zenetos (2002) consider E. vittata, G. oculata, O. fusiformis and A. acutifrons to be tolerant species. Differences in the prevalence of dominant species were noticed in different research periods. This state, however, is not necessarily connected with the repeated assemblage perturbations, but could be the result of natural variations in the community characterized by multiple stable points. Most sediment-living communities show large changes in dominance through time, and cannot be encompassed within the formal stability theory asserting that in the stable populations the dominance should remain at a uniform level (Gray 1977). Apart from changes in dominance due to the variation of the physical environment, or those induced by anthropogenic perturbations, changes can be induced by different biological interactions such as the competition between species, predation on certain species, differential larval settlement and the activities of the animals themselves rendering the sediment unsuitable for their own or other species. According to Gray (1977), marine sediment communities typically have multiple stable points and show neighbourhood stability and correspond to the poly-climax theory of community development (Gleason 1926). Low oxygen concentrations can distract the equilibrium of the whole ecosystem of a certain area, as different groups of organisms show different tolerance levels (Hagerman & Szaniawska 1989). This is very well represented in the NA through the altered relations between the most important macrobenthic groups. Previous researches in the area showed the dominance of bivalves over the polychaetes and echinoderms in communities affected by the different stress conditions, such as dithering, trawling and lack of oxygen, while with the shift toward more favourable conditions the dominance of polychaetes and echinoderms became re-established (Crema et al. 1991; Jaklin 1992; Hrs-Brenko et al. 1994; Zavodnik et al. 1994, Šimunović et al. 1999). The results of the analyses of the whole macrobenthic fauna on our research stations from the first period (1989–1994) confirmed that benthic communities were in the unstable transition period, probably because of the constant environmental pressure of repeated low-oxygen episodes (Jaklin 2002). Independently from the differences between the stations, bivalves were the dominant macrobenthic group overall among the surviving organisms over the whole research area (over 40% in overall diversity and almost 60% in overall macrobenthos abundance). However, during the recovery, the contribution of bivalves to overall diversity and abundance gradually decreased. On the other hand, the importance of polychaetes in the community grew continuously over time. After about one and a half years post- anoxia, polychaetes became more or less equal to bivalves in the community in terms of their diversity and abundance. After two years, they became the dominant macrobenthic group considering both parameters, and retained that status until 1994 (four years after the anoxia) (Jaklin 2002).

Although not recorded at our research stations, anoxic events and mass mortalities were sporadically recorded in other areas of the NA during the last 20 years (Malej et al. 1991; Penna et al. 2004). According to Stachowitsch and Fuchs (1995), the recovery of the macrobenthic community after the mass mortality in another area of the NA was not reached even 10 years after the anoxic event, possibly due to repeated disturbances. In the NA, instability has been introduced by recurring perturbations involving hypoxia/anoxia along with intensive dredging and trawling activities (Kollmann & Stachowitsch 2001). Currently the frequency of such disturbances greatly exceeds the duration of recolonization processes (Riedel et al. 2008). The recurrence of bottom disturbances is probably too high in relation to the time required for community reconstruction and stabilization (Crema et al. 1991). The situation of NA has been described as ‘rapid death, slow recovery’ (Stachowitsch 1991). Similar sources of disturbances (trawling, dredging, hypoxias) may also influence the presently analysed polychaete assemblages, which, manifestly, still show high levels of variation. Low oxygen concentrations (< 2 ml/l) that are retained for longer times can cause disturbances in reproduction and feeding, decreases in abundance or disappearance of oxygen-sensitive species, and in extreme cases (< 0.5 ml/l) mass mortality of benthic organisms (Diaz & Rosenberg 1995). After the anoxic event in 1989, low oxygen concentrations at station 107 were recorded only in September and October 1991. This could have slowed down the recovery process of the polychaete assemblages during those years. Since station 107 is the most critical among the three investigated considering minimum oxygen levels, it can be supposed that the other two stations did not suffer from the oxygen deficiencies. Our results underline that changes in oxygen after the anoxic event did not relate to the pattern of assemblage variations. Salinity is an important variable in explaining the distribution patterns of benthic assemblages in estuarine or lagunar conditions (e.g. Schirosi et al. 2010). NA is highly influenced by the Po river waters, which in years of higher inflow markedly decrease the salinity and increase the concentrations of nutrients, stimulating algal blooms followed by a decrease of oxygen in the bottom layer (Degobbis et al. 2000). Our analyses, however, show very low correlation of salinity with the other variables, including oxygen. The variation of salinity through years after the anoxia (min. 37.88; max. 38.66) was also not large enough to likely suggest direct effects on the benthic fauna.

Explanatory models about the causes that determine changes in biodiversity are difficult to formulate over the long term. Correlative analyses of changes in assemblages in relation to abiotic factors are perhaps the only tools available to unravel this issue. There is global awareness on the ongoing changes in marine biodiversity as a result of a wide range of large-scale human disturbances, particularly climate change. Long-term monitoring programmes are essential for quantifying the extent of such effects. Unfortunately, although many programmes were implemented from the first half of the twentieth century, most were considered ‘simply monitoring’ by administrators and stopped during the 1980s. Our study is a clear example of how long-term analyses at fine taxonomic level based on well-designed monitoring programmes should be supported, as they can provide precious insights to the way we perceive changes in marine biodiversity over time as a result of large-scale natural changes and human impact.

Acknowledgements

We thank C. Arvanitidis, R. Sardá and the third anonymous reviewer for their helpful comments, which improved the manuscript. BM would like to thank to Elvis Zahtila (Public Institution ‘Natura Histrica’) and Andrej Jaklin (CMR, Ruđer Bošković Institute) for taxonomic support. The results of this investigation as well as the work of BM were carried out within the projects ‘Mechanisms of long-term changes in the Adriatic Sea ecosystem’ (project no. 0098111) and ‘Biodiversity of benthic communities in the Adriatic: natural and human impacts’ (project no. 098-0982705-2732), both funded by the Ministry of Science, Education and Sports of the Republic of Croatia. The work of LM, AG and AT has been carried out within the MARBEF Network of Excellence which is funded in the Community's Sixth Framework Programme (Contract no. GOCE-CT-2003-505446), the EU Integrated Project SESAME and the CMCC (Centro Euro-Mediterraneo per i Cambiamenti Climatici).