Polycyclic Aromatic Hydrocarbons Have Adverse Effects on Benthic Communities in the Baltic Sea: Implications for Environmental Status Assessment

Raymond, Caroline; Gorokhova, Elena; Karlson, Agnes ML

doi:10.3389/fenvs.2021.624658

ORIGINAL RESEARCH article

Front. Environ. Sci., 12 April 2021

Sec. Toxicology, Pollution and the Environment

Volume 9 - 2021 | https://doi.org/10.3389/fenvs.2021.624658

Polycyclic Aromatic Hydrocarbons Have Adverse Effects on Benthic Communities in the Baltic Sea: Implications for Environmental Status Assessment

$\r\nCaroline Raymond$ Caroline Raymond¹

Elena Gorokhova²

Agnes ML Karlson^1,3*

¹Department of Ecology, Environment and Plant Sciences (DEEP), Stockholm University, Stockholm, Sweden
²Department of Environmental Science, Stockholm University, Stockholm, Sweden
³Stockholm University Baltic Sea Centre, Stockholm, Sweden

Changes in benthic macrofaunal communities are indicative of environmental stressors, including eutrophication and hypoxia. However, some species are sensitive not only to hypoxia but also to various environmental contaminants. We tested which of the environmental predictors (sediment organic carbon, sediment concentrations of metals and polyaromatic hydrocarbons [PAHs], bottom water oxygen, salinity, temperature, and surface chlorophyll-a concentration) that best explained the following response variables: (1) macrofauna community composition, (2) abundance of a benthic sentinel species, the amphipod Monoporeia affinis; and (3) the Benthic Quality Index (BQI). All data originated from 29 reference monitoring stations in the Baltic Sea and the statistical tests included both uni- and multivariate analyses. The community composition and BQI were best explained by the same combination of salinity, depth, temperature and PAH concentrations. The abundance of M. affinis, which is sensitive to hypoxia and chemical exposure, was best explained by PAHs as a single predictor. Our findings suggest that benthic communities in the Baltic Sea are influenced by anthropogenic contaminants, which should be taken into account when benthos is used for eutrophication status assessment.

Introduction

Anthropogenic impacts, such as pollution with nutrients and hazardous substances, together with the progressing climate change threaten biodiversity and ecosystem integrity worldwide (Halpern et al., 2008; Cloern et al., 2016). Benthic communities in marine soft sediments, the world’s second largest habitat, are increasingly subjected to multiple stressors resulting from these pressures. Hypoxia as a consequence of eutrophication, and to some extent increasing temperatures, are considered the major threat (e.g., Breitburg et al., 2018; Ehrnsten et al., 2020). In coastal environments, benthic ecosystems are also commonly exposed to contaminant loading. Benthic communities, hereafter macrobenthos, respond relatively fast to environmental disturbances but with species-specific sensitivity/tolerance to stressors (Pearson and Rosenberg, 1978). The composition of macrobenthos integrates environmental conditions over longer periods (years to decades). Therefore, macrobenthos is commonly used as indicators of ecosystem status and health (Pearson and Rosenberg, 1978; Dauer, 1993; Borja et al., 2000; Muniz et al., 2005).

The Baltic Sea is strongly influenced by eutrophication and contaminants and is often considered one of the most polluted sea areas in the world (Elmgren, 1989; Karlson et al., 2002; Belkin, 2009; Carstensen et al., 2014). In addition, it has a naturally low species richness due to its low salinity (Elmgren, 1984), and it is one of the most studied ecosystems making it suitable for potential future scenarios for other sea areas (Elmgren and Hill, 1997; Reusch et al., 2018). The restricted water exchange from the permanent halocline and seasonal pycnoclines in coastal water results in a natural oxygen deficiency, which is exacerbated by eutrophication (Conley et al., 2009, 2011; Carstensen et al., 2014). Anoxia has increased several-fold since the early 1990s, leading to the elimination of macrobenthos below the halocline at 70-80 m depth in the Baltic Proper (Karlson et al., 2002; Conley et al., 2009, 2011), but also shallower areas can be exposed to seasonal or episodic hypoxia (Gammal et al., 2017). Increasing organic enrichment and hypoxia are likely to alter macrobenthic species composition via differential effects on sensitive and tolerant species (Pearson and Rosenberg, 1978). However, in the species-poor Baltic Sea, such changes in species composition due to variability in oxygen concentrations have not been found significant (e.g., Villnäs and Norkko, 2011; Rousi et al., 2013). Nevertheless, macrobenthos and various metrics based on the community structure are broadly used as indicators of ecological effects of eutrophication with low oxygen conditions being considered the primary cause for potential changes in species composition (HELCOM, 2009; Leonardsson et al., 2009).

The Water Framework Directive (WFD) and Marine Strategy Framework Directive (MSFD) provide a general framework to evaluate the environmental quality of biological, hydromorphological, and physicochemical characteristics of water bodies in Europe (Council directive 2000/60/EC). Assessing the environmental quality of soft-bottom macrofauna in coastal and open sea areas differs between countries bordering the Baltic Sea and its sub-basins (HELCOM, 2009) since the strong salinity gradient and habitat types determine the species diversity (Elmgren, 1984; Leppäkoski and Bonsdorff, 1989; Rumohr et al., 1996). Thus, no single index has been applicable to capture these differences, and each country uses different benthic community metrics to assess the state of the environment. In the northern Baltic Proper and the Bothnian Sea, where Sweden and Finland are the bordering countries, two similar indices are used for assessment of benthic quality: Sweden uses Benthic Quality Index (BQI), and Finland uses Brackish-water Benthic Index (BBI) (Perus et al., 2007; HELCOM, 2009, 2018b; Leonardsson et al., 2009). Both indices build on abundance weighted proportion of sensitive to tolerant taxa and the number of species within the community, and good benthic quality is associated with a high proportion of sensitive taxa and a high number of species. The species sensitivity values are primarily based on tolerance to low oxygen levels and not to other stressors, such as environmental contaminants.

Organic contaminants and heavy metals are, however, considered major threats to the Baltic Sea (HELCOM, 2010). Many environmental contaminants remain in the Baltic Sea for a long time because of the slow water exchange (Meier, 2007), the persistent nature of many pollutants, and the association of hydrophobic chemicals to sediments (Axelman et al., 2001). In particular, polycyclic aromatic hydrocarbons (PAHs) and heavy metals cause various adverse effects on macrobenthos (Sundelin, 1983; Long et al., 1995; Landrum et al., 2003; Martins et al., 2013; Löf et al., 2016). The amphipods Monoporeia affinis and Pontoporeia femorata are key components in the benthic communities, ranked among the most hypoxia-sensitive species, but their reproduction and embryo development are also highly sensitive to contaminants (Sundelin and Eriksson, 1998). In the MSFD, embryo development of M. affinis is, therefore, used as a supplementary indicator of biological effects of contaminants (HELCOM, 2018a).

Here, we tested which of the environmental predictors (sediment organic carbon, sediment concentrations of metals and PAHs, bottom water oxygen, salinity, temperature, and surface chlorophyll-a concentration) that best explained the following response variables: (1) macrobenthic species composition, (2) abundance of a macrobenthic sentinel species, the semelparous amphipod Monoporeia affinis (life cycle of 2 years); and (3) the Benthic Quality Index (BQI). The data originated from 29 monitoring stations located along the Swedish Baltic Sea coast with no known point sources of pollution.

We specifically asked three questions.

(1) Can variability in the macrobenthic species composition be better explained by including contaminant data to environmental variables?

(2) Can variability in the abundance of the sentinel species Monoporeia affinis be better explained by including contaminant data?

(3) Does the benthic quality index (BQI), aimed for assessing eutrophication effects, also correlate with sediment contaminant levels?

Materials and Methods

Site Descriptions and Study Species

The Baltic Sea is a semi-enclosed non-tidal brackish water body, the largest in the world area-wise. The two largest basins are the Baltic Proper (BP) and the Bothnian Sea (BoS) (Figure 1). The Baltic Sea is characterized by a permanent salinity gradient that is decreasing northwards and a permanent halocline at about 80 m depth in the Baltic Proper. The low salinity with mean surface salinity of 6-7.5 and 4-5 in BP and BoS, respectively, and the short evolutionary time because of the recent geological history, have resulted in a uniquely low species richness of benthic macrofauna (Elmgren and Hill, 1997).

FIGURE 1

Figure 1. The 29 stations (color-coded) within regions along the Swedish Baltic Sea coast. Regions 1-3 in the Bothnian Sea (BoS) and regions 4-8 in the Baltic Proper (BP). The Swedish names for the regions (and number of stations included in each region) are as follows; BoS-1 Norrbyn (1 station), BoS-2 Höga kusten (4 stations), BoS-3 Söderhamn (5 stations), BP-4 Stockholm (7 stations), BP-5 Askö (5 stations), BP-6 Hävringe (3 stations), BP-7 Östergötland (3 stations), and BP-8 Utklippan (1 station).

The soft-bottom communities in the Baltic Sea deeper than 20 m are thus represented only by a handful species, where the long-lived (up to 30 years; Segerstråle, 1960) tellenid clam Limecola balthica (previously known as Macoma balthica) often dominates the biomass, except for the northernmost parts. Two species of semelparous amphipods with a 2 year life-cycle, Monoporeia affinis (a glacial relict) and the morphologically similar Pontoporeia femorata (of marine origin), often dominates the abundance in the Baltic Proper although in lower numbers now than in the 1970s (Ankar and Elmgren, 1976; Rousi et al., 2013). In the Bothnian Sea, M. affinis dominates the deposit-feeding guild. In recent decades, the non-indigenous polychaetes Marenzelleria (three morphologically similar species; Bastrop and Blank, 2006) have successfully colonized the entire Baltic Sea, and since the early 2000s, it is a dominant genus (Kauppi et al., 2015). The most abundant predatory invertebrates in the Baltic Proper are the isopod Saduria entomon, the polychaete Bylgides sarsi, and the priapulid Halicryptus spinulosus; the latter two are of marine origin and less abundant in the Bothnian Sea.

Sampling Methods and Data Extraction

Macrobenthos was collected within the Swedish National and Regional monitoring program for benthic macrofauna with a yearly (May to early June) sampling along the Swedish coast. One van Veen grab sample (0.1 m²) is collected at each station, sieved onboard over a 1 mm mesh, and the retained material is fixed with 4% buffered formaldehyde. High spatial coverage instead of replication within stations is recommended for Baltic macrobenthic monitoring (SEPA, 2007); variation between replicate samples in this species-poor system is low, and species rarefaction curves are nearly saturated after a single sample (Villnäs and Norkko, 2011). The macrobenthic organisms are identified to the lowest possible taxonomical level in the lab, counted, and weighed (see species list in Supplementary Table 2). The sampling, sample processing, and taxonomic analysis follow the European standard (CSN EN ISO 16665, 2013).

Sediment samples were collected within the program for biological effect monitoring of contaminants, carried out since 1994 in the Swedish Baltic Sea coastal areas by assessing embryo development in the sentinel amphipod species M. affinis and P. femorata. The sampling is conducted in January, when gravid females can be obtained (Sundelin and Eriksson, 1998; Sundelin et al., 2008; HELCOM, 2018a) using a benthic sledge set to sample the upper 2-3 cm of the sediment. From the year 2012, the effect monitoring was expanded to cover a larger geographical region and improve the spatial overlap with the Swedish National and Regional monitoring program for benthic macrofauna. As a result, 30 matching stations within the benthic macrofauna monitoring program were established. One station (SR1A, BoS-3 region) was populated only by M. affinis at ten-fold higher densities than the median density; hence, it was considered an outlier and excluded from all statistical analyses. Thus, the study is restricted to 29 stations; 10 of them located in the Bothnian Sea and 19 stations in the Baltic Proper. For visualization, these stations are grouped into eight regions dependent on their latitude; Figure 1. All stations included in this study (Supplementary Table 1) are representing areas without any known point sources of pollution.

Salinity, temperature and oxygen concentrations are measured simultaneously in the bottom water when the macrobenthic sampling is carried out; however, these data represent only a snapshot of these abiotic variables. Moreover, there is no sampling of Chlorophyll-a (Chl-a) from the monitoring stations used in this study. We, therefore, extracted modelled environmental data from DAS (Data Assimilation System¹; Sokolov et al., 1997; Savchuk, 2018), which is a program that constructs a 3D grid of hydrographic and chemical monitoring data collected from all countries around the Baltic, stored in a database at the Stockholm University Baltic Sea Centre. Extrapolated data for our 29 stations and a relevant period (see below) are based on interpolations of data from varying spatial resolution across time (Supplementary Table 1 and Supplementary Figure 1). Temperature, salinity, and oxygen concentrations from August-October, the time of the year when oxygen minimum and temperature maximum could be expected to affect macrobenthos the following spring, were extracted from the bottom water for each site; note that the bottom depth varies across stations (Supplementary Table 1). In the upper part of the water column (0-10 m), the Chl-a concentrations in spring (average March-May) and summer (average June-August) was extracted from each site. As there was not enough Chl-a data in the Bothnian Sea during the springtime to allow sufficient resolution for site-specific values in this basin, the spring Chl-a was only used for analyses in the Baltic Proper (Supplementary Figure 1).

Chemical Analyses

The surface sediment samples from each station were immediately frozen (−20°C). Samples of 20 mg were analyzed for total organic carbon (TOC; g kg^–1 dry weight [dwt]) content after drying and homogenization. The TOC values were determined using a Flash 2000 (Thermo Scientific) elemental analyser at Stockholm University, Dept. Ecology, Environment and Plant Science, and expressed as a percentage per dry weight (dwt). Acidification is not necessary since the inorganic carbon content of Baltic sediment is negligible (J. Walve, Stockholm University, pers comm).

About 50 g of the frozen sediment were used for metal analyses (As, Cd, Co, Cr, Cu, Hg, Ni, Pb, V and Zn) at ALS Scandinavia, Sweden, and the same amount of sediment for PAH analyses (11 priority congeners; SEPA, 1999; Josefsson, 2017) at the Swedish Environmental Research Institute (IVL), Gothenburg, Sweden. For three stations, the sediment from the year 2012 was lost during the sample processing; therefore, sediment from the year 2014 from the same station and same sampling month was used instead (Supplementary Table 1). Since the sediment is collected with a sledge, which integrates material from approximately 2-3 cm, the measured contaminants and organic carbon content represent the conditions in this top layer, regardless of the yearly representation within this short period of time (2 years). Metal concentrations are expressed as mg per gram dwt sediment and PAHs as μg per gram dwt sediment (Supplementary Tables 3, 4).

The measured concentrations of PAHs were compared to the Swedish sediment quality guidelines for PAHs in Baltic sediments (SEPA, 1999; Josefsson, 2017). PAHs are classified on a 5-grade scale, from insignificant to very high deviation established for both individual PAHs and the sum of 11 priority PAHs (SEPA, 1999; Josefsson, 2017). For metals, we used the risk quotient (RQ) approach. Here, the measured environmental concentration (MEC) is compared to the reference sediment values in Löf et al. (2016). The MEC-to-reference values ratio is the RQ, and RQ exceeding 1 for a metal suggests that it may pose a toxicity risk to the environment (EFSA, 2013; Arias-Andres et al., 2018). The sum of these RQ values was calculated for each station and used in statistical analyses since all metals correlated significantly and highly to RQ (Pearsons r ≥ 0.90 except for As, Hg and Co, which had values of 0.73-0.89; Supplementary Table 6) and to each other (mean r = 0.80, with only exception being As and Hg, which were not significantly correlated; Supplementary Table 6). Similarly, all individual PAHs correlated significantly and highly to sum PAHs (Pearsons r ≥ 0.86 except for anthracene with r = 0.70; Supplementary Table 5), and all PAHs were significantly correlated to each other (mean r = 0.80, except for anthracene and benso(k)fluoranthene, which were not significantly correlated; Supplementary Table 5). Threshold values based on benthic toxicity data, so-called Environmental Quality Standards for sediment (EQS_sediment; EC, 2011) are available for four compounds; the PAHs anthracene (24 μg kg^–1 dwt, standardized for 5% TOC) and fluoranthene (2,000 μg kg^–1 dwt, for 5% TOC) and the metals Pb and Cd (120 and 2.3 mg kg^–1 dwt, respectively, not TOC standardized). The EQS_sediment values were compared to our measured concentrations (standardized for TOC when relevant).

Macrobenthic Index Calculation

The Benthic Quality Index (BQI) based on the abundance-weighted proportion of sensitive to tolerant taxa and the number of species within the community was calculated for each station and used as a univariate response variable in statistical analyses as described in section “Statistical Analyses”. The BQI index is defined as:

\begin{matrix} B Q 1 = [\sum_{i = 1}^{{}^{S}c l a s s i f i e d} (\frac{N_{i}}{N_{c l a s s i f i e d}} * S e n s i t i v i t y v a l u e_{i})] * \\ \log_{10} {(S + 1)}^{*} (\frac{N}{N + 5}) \end{matrix}

where S_classified is the number of taxa having a sensitivity value, N_i is the number of individuals of taxon i, N_classified is the total number of individuals of taxa having a sensitivity value, the Sensitivity value_i is the sensitivity value for taxon i, S is the total number of taxa, and N is the total number of individuals in the sample (0.1 m²). The species are classified on a sensitivity-disturbance scale, based on the species oxygen tolerance reported in the literature and expert knowledge (Leonardsson et al., 2009).

Statistical Analyses

The macrobenthic data were restricted to the stations with available contaminant data (BP: 19 stations, and BoS: 10 stations). We carried out ordinations, permutational variance analyses, and regression analyses as described in the following steps.

Step 1: Establishing Gradients in Macrobenthic and Environmental Data

First, the macrobenthos data from 2011-2014 were used to test whether the year 2012 was a representative year regarding the community composition. The year 2012 was of special interest for macrobenthos because of the high resolution Chl-a and physicochemical data available from summer 2011 (lagged effect on benthos) and benthic contaminant data available from 2012 (see section “Sampling Methods and Data Extraction”). To visualize the variability in the macrobenthos abundance related to the region and sampling year, we used non-metric multidimensional scaling (nMDS; Clarke and Gorley, 2006), based on the Bray-Curtis similarity index. The abundance data were non-transformed prior to creating the resemblance matrix since subtle changes in dominant species are relevant. However, for comparison, all analyses are also performed on fourth-root transformed data downweighing dominant species (Supplementary Figure 2). Species composition (based on Bray-Curtis similarity index) was tested for the effect of Year (four levels; 2011, 2012, 2013, and 2014) and Region (eight levels; see map in Figure 1) and their interaction (Year x Region) using Permutational analyses of variance, PERMANOVA (PERMANOVA + in PRIMER v6; Anderson, 2001). As a complement, to test for differences in beta-diversity among Regions, Jaccard similarity index on presence/absence data was also calculated and used in PERMDISP (Anderson et al., 2006).

Thereafter, multidimensional scaling using principal coordinate analysis (PCO) was performed on the species composition (based on the Bray-Curtis dissimilarity of abundance data) for the year 2012 only (since the analysis above confirmed that it was a representative year) and, separately on the environmental/chemical data (Euclidean Distance). We considered species and environmental data with a Spearman correlation greater than 0.3 with any of the first two ordination axes as significantly contributing to the differences between the stations. For all data types (abundances, environmental, and chemical data), the stations were grouped and visualized by Region as shown in Figure 1.

Step 2: Linking Macrobenthic Data to Environmental Data

Distance-based linear models (DistLM in PERMANOVA + for PRIMER 6; Anderson et al., 2008) is a multiple regression modeling, where the resemblance matrix (here the Bray-Curtis distances for species composition) is regressed against a set of explanatory variables. As shown in step 1, the year 2012 was representative for macrobenthic community development and, therefore, the abundance data for this year were used as response variable in DistLM analyses with environmental data for 2011 that were extracted from DAS (section “Sampling Methods and Data Extraction”) and contaminant levels (section “Chemical Analyses”) as explanatory variables. The time spans of the environmental and contaminant data were considered as ecologically relevant for the benthic community data from the monitoring survey 2012. The univariate responses for the abundance of Monoporeia affinis (recorded at all stations; resemblance based on the Bray-Curtis dissimilarity) and the BQI index (resemblance based on the Euclidean distance) were also analyzed by the DistLM.

The skewness of the explanatory variables was inspected using pair-wise Draftsman plots prior to analysis. The explanatory variables were generally not strongly correlated to each other (Anderson et al., 2008), and distributions were not strongly skewed (Supplementary Table 7). To validate the robustness of the results, both non-transformed and log-transformed contaminant concentrations were tested as predictors. Similarly, we also run all models with fourth-root transformed community composition data, which downweigh the importance of dominant species as stated in Step 1. Finally, as no spring Chl-a data were available for the BoS, a model with spring Chl-a data as an additional explanatory variable was analyzed for the Baltic Proper stations (n = 19).

Marginal DistLM was used to determine the contribution of each predictor to the macrobenthic species composition, BQI and M. affinis abundance when other predictors are ignored. The variables included in the final DistLM output were selected using the stepwise selection procedure, which utilizes all possible combinations of the explanatory variables to determine the combination accounting for the best improvement of the Akaike information criterion corrected for a small sample size, AICc. To confirm that results were robust to the criteria used, we also performed the DistLM analyses with forward and backward selection and adjusted R² as a stop criterion instead of AICc. The multivariate model (macrobenthic community) was visualized using a db-RDA (McArdle and Anderson, 2001). Finally, since organic carbon content and contaminant levels were highly correlated (Supplementary Table 7), we performed distLM analyses for the three models (macrobenthic species composition, M. affinis and BQI) with the variables PAH and metals that were force-excluded to compare the results between the concurrent models using the percentage of the overall variance explained.

Step 3. Predictive Modelling

To determine whether there was a significant relationship between the species composition and each environmental variable, a canonical analysis of principal coordinates, or CAP (Anderson and Robinson, 2003; Anderson and Willis, 2003) was applied. This allows a constrained ordination on the basis of, e.g., Bray–Curtis similarity, determining the axes that best discriminates an environmental gradient. A separate CAP analysis was done for each dependent variable (Ellis et al., 2015).

Results

Biotic Gradients

In total, 17 species were recorded from the 29 stations (Table 1), with 3 to 11 species per station (Supplementary Table 2). Densities of the three most abundant macrobenthic species, i.e., the sentinel amphipod Monoporeia affinis, the tellenid bivalve Limecola balthica, and the non-indigenous polychaetes Marenzelleria spp., ranged from 10-30 ind. m^–2 to 2,500-3,900 per m^–2, with species-specific median densities of 825 (M. affinis), 687 (L. balthica), and 1271 (Marenzelleria spp.). Species composition with details on covariation among species is shown in Figure 2A.

TABLE 1

Table 1. Species list of the 29 stations with species dominance and frequency of occurrence (number of stations with recorded species).

FIGURE 2

Figure 2. Principal coordinate analysis PCO for year 2012 on (A) biological data (Bray-Curtis similarity on abundance data, untransformed), and (B) environmental and chemical data (Euclidean distance). See Table 1 for a complete list of abundance and occurrence of all species and Supplementary Table 2 for the station-specific list.

There was no significant effect of neither Year nor Year x Region interaction for macrobenthic species composition (Supplementary Figure 2 and Supplementary Table 8). By contrast, the effect of Region on the species composition was significant; differences between the regions were found both within and across basins (Supplementary Figure 2 and Supplementary Table 9). The Region effect on the beta diversity was not significant (Permdisp F₇,₂₁: 4.2753, p < 0.09).

The BQI values at all stations except two (both in the southern Bothnian Sea, Figure 3) reached the threshold of good environmental status (assessment of a water body should however be based on at least five stations, see Leonardsson et al., 2009). The BQI index was significantly positively correlated with the abundance of M. affinis, H. spinulosus and P. femorata (Supplementary Table 10).

FIGURE 3

Figure 3. BQI values for each station for the year 2012, error bars show max and min values for the period 2011-2014. The horizontal dashed line is the threshold for good ecological quality assessment (added for comparison but please note that assessment of a water body should be based on at least five stations, see Leonardsson et al., 2009).

Environmental Gradients

Salinity differed between the basins, with average values 5.5 (range 5.2-5.8) and 6.5 (5.8-7.1) at the BoS and BP stations, respectively. The corresponding values for temperature were 5.2 (2.3-8.8°C) and 8.0°C (4.8-11.8°C) in BoS and BP, respectively. The other variables were not significantly different with regard to basin factor, thus, showing no clear latitudinal gradients (Supplementary Table 1). To summarise, the range in oxygen concentrations for the BoS and BP basins were 6.7-8.1 and 6.7-7.4 mg O₂ L^–1; organic carbon content in sediment was 0.3-2.7 and 0.3-6.1% of dwt and summer surface water Chl-a 1.6-5.0 and 2.0-6.3 mg L^–1, respectively. The overall correlations between the environmental variables can be seen in Figure 2B and univariate correlations in Supplementary Table 7. Org C% was significantly positively correlated to metal RQ and sum PAHs, while the temperature was significantly negatively correlated to metal RQ.

PAH concentrations were classified as medium to high values according to national guidelines at 4 stations (Josefsson, 2017); in addition, 8 stations exceeded the pristine level (Supplementary Table 1). Metals exceeded the total RQ (of 10) at 7 stations. The effect-based threshold values (EQS_sediment) for fluoranthene was never exceeded, foranthracene it was exceeded at the southernmost station (BP-8.1). Eight stations in the Baltic Proper and one station in the Bothnian Sea had similar or higher values than EQS_sediment for bensofluor(b and k)anthen and the EQS_sediment for benso(a)pyrene was exceeded at 3 stations in the Baltic proper. Concentrations of Pb and Cd were always below the EQS_sediment.

Macrobenthic Data Explained by Environmental Data

The best-fit models generated by the DistLM analyses for 1) multivariate macrobenthic species composition data, 2) abundance of M. affinis, and 3) the BQI index, are summarized in Table 2; marginal tests are reported in Supplementary Table 11. The same explanatory variables (salinity, depth, temperature, and PAH concentrations), albeit in different order of importance, were included in the species composition model and the BQI model; they all were significant, except for temperature in the BQI model. The M. affinis abundance was best explained by the sumPAHs as a single predictor. In the models with PAHs and metals force-excluded, organic carbon replaced PAH as a fourth (community data) and second (BQI) most influential predictors and the only predictor for M. affinis. However, the total variance explained was lower in the two models (35% for the community, and 15% for M. affinis; Table 3) compared to the models allowing PAH as a possible predictor. The BQI model was not affected (58 instead of 57%) by adding a fifth variable instead of four variables when the PAHs variable was included (Supplementary Table 12).

TABLE 2

Table 2. DistLM output from sequential tests for the models for community composition, BQI and M. affinis abundance based on stepwise selection and AICc as stop criteria (see text).

TABLE 3

Table 3. Summary of the cumulative proportion explained for each DistLM model (community composition, BQI and M. affinis); when abundance data was untransformed and fourth-root transformed, when all environmental variables initially included were untransformed and when PAHs and RQmetals were log(x + 1)-transformed, and, finally when the variable PAHs was force-excluded (due to its high correlation with orgC%).

When spring Chl-a was considered as a possible predictor in addition to the summer Chl-a and the other environmental variables (only possible for Baltic Proper, see methods), the DistLM analysis resulted in a final model with depth and sumPAH as predictors (Supplementary Table 13).

The correlation between PAH and the response variables (BQI and M. affinis) were negative, while salinity was positively correlated to BQI. A db-RDA is illustrating the relationships between species and environmental variables for the macrobenthic species composition model (Figure 4).

FIGURE 4

Figure 4. db-RDA with the first two axes explaining 31.8% of the total variation. Environmental data selected in the best distLM model (Table 2) are shown as vectors (blue) and species characterizing the stations (gray) are overlaid to facilitate the interpretation.

In the CAP analyses, a significant community change was observed along the PAH gradient (Figure 5, R² = 0.34, p < 0.04); however, correlations with salinity (R² = 0.70, p < 0.001), temperature (R² = 0.49, p < 0.015) and depth (R² = 0.52, p < 0.001) were stronger and more significant. Correlations with the other variables were not significant (org C: R² = 0.29, p > 0.08, oxygen: R² = 0.23, p > 0.08, Chl-a: R² = 0.07, p > 0.6, metals: R² = 0.39, p > 0.08).

FIGURE 5

Figure 5. CAP axis 1 correlated with PAH concentrations (sum of the 11 priority PAHs). Horizontal dashed lines represent the classification for PAH contamination in relation to pristine status and based on the Swedish EPA guidelines (SEPA, 1999; Josefsson, 2017); the thresholds are set to 170 and 440 representing thresholds between very low to low levels and low to medium high levels, respectively.

Discussion

This study shows that at reference stations within the Swedish macrobenthic monitoring program along the Baltic Sea coast, concentrations of 11 priority congeners of PAHs in sediment often exceed the level “low” in the national guidelines (Josefsson, 2017). Our analysis suggests that PAHs, together with non-chemical factors, such as salinity, temperature and depth, contribute significantly to variation in the macrobenthic species composition. Moreover, the same combinations of predictors explained the Benthic Quality Index (BQI) variation, indicating that BQI adequately represents the macrobenthic community state. Interestingly, the sum of PAHs alone was the best negative predictor of Monoporeia affinis abundance. The importance of PAHs in explaining changes in macrobenthic species composition and abundance of M. affinis suggests that populations and communities in sea areas with no known point sources of pollution are adversely affected by these environmental contaminants.

The positive relationship between PAHs and organic carbon (orgC%) in sediment was expected because these compounds, especially high-molecular weight (hmw) PAHs, are lipophilic and sorb to organic particles associated with the sediment. However, after the force-exclusion of PAHs and replacing it with orgC, the models for both community composition and M. affinis abundance as dependent variables were weaker. Thus, PAHs had greater explanatory capacity than organic carbon. The environmental quality standards (EQS_sediment) threshold value was exceeded for anthracene in one station (the southernmost). Furthermore, high concentrations relative to the indicative EQS values of other PAHs listed as priority substances under the EU Water Framework Directive, WFD (Council directive 2000/60/EC, 2000) – benzo[b and k]fluoranthene and benzo[a]pyrene – were found at several stations in the Baltic Proper. High sediment concentrations of these PAHs have previously been reported in other Baltic areas without known point sources (Pikkarainen, 2004b); however, other PAHs, e.g., chrysene, pyrene, benzo[a]anthracene and fluoranthene bioaccumulated most in the clam L. balthica inhabiting these sediments (Pikkarainen, 2004a). Thus, the sediment concentrations do not necessarily translate into body concentrations in macrobenthos because of the variability in bioavailability (driven by sediment properties, grain size distribution, and organic content, etc.) and physiological/biochemical properties of the organisms.

PAHs mainly affect invertebrates via narcosis, i.e., non-specific toxicity (Douben, 2003), endocrine disruption, and genotoxicity (White et al., 1997; Douben, 2003; Zhang et al., 2016; Gorokhova et al., 2020). A majority of PAHs (19 out of 24) have shown positive correlations to various reproductive disorders in M. affinis, especially to the occurrence of dead broods (Löf et al., 2016). Macrobenthic species can, however, respond to PAH exposure in different ways. Following the large Tsesis oil spill in 1977 that took place in the Baltic Sea archipelago, close to the Askö area, where several of the highest concentrations were also measured in this study, the bivalve L. balthica population was relatively unaffected despite the considerable contamination by oil hydrocarbons. Conversely, the amphipods M. affinis and P. femorata showed drastically reduced population densities for at least three years after the accident, as well as increased malformed embryos (Elmgren et al., 1983). The macrobenthic species composition may, therefore, at least partly, reflect contaminant history in the area.

Metals were not identified as significant predictors despite the high RQ values of metals at several stations, suggesting that their effects were low compared to PAHs. Heavy metals should not, however, be dismissed as a potential toxicity source to the benthos. Copper, zinc, and lead at levels below sediment quality guidelines in New Zealand and Australia have been linked to altered benthic community composition and a decline in large species, i.e., infaunal bivalves (Hewitt et al., 2009). In shallow Baltic harbors with zinc and copper concentrations comparable to those measured in our study, gastropods suffered sublethal effects as well as increased mortality (Bighiu et al., 2017a, b). In the EU WFD, mercury, cadmium and lead are considered priority substances (Council directive 2000/60/EC, 2000); these metals are particularly toxic to crustaceans (e.g., Connor, 1972; Greenwood and Fielder, 1983; Giudici and Guarino, 1989; Reddy et al., 1997; Itow et al., 1998). The measured values in our study was below the EQS_sediment values for cadmium and lead, however, cadmium concentrations comparable to those in our study were found to increase embryo aberrations in M. affinis (Sundelin, 1983; Löf et al., 2016).

Despite the fact that the Baltic Sea is considered one of the most contaminated sea areas in the world (Reusch et al., 2018), sediment contaminants have not previously been included in models of macrobenthic species composition from sites that are not directly affected by pollution. Depauperate macrobenthic communities in heavily contaminated areas compared to reference sites have, however, been reported in the southern Baltic Sea (Bettinetti et al., 2009). Exposure to multiple stressors, e.g., contaminated sediments and hypoxia, can lead to a decreased capacity to cope with these stressors, as shown for M. affinis (Gorokhova et al., 2010, 2013). In our study, none of the environmental variables linked to the eutrophication, i.e., oxygen concentrations in the bottom water, the organic carbon content of sediment and water column chlorophyll-a during summer and spring, were identified as significant predictors for any of the three response variables: species composition, the BQI index and abundance of the sentinel species M. affinis. This raises concerns about using BQI (or benthic community composition in general) as eutrophication indicators and highlights the problem of evaluating anthropogenic pressures in isolation. Based on the correlation between BQI and contaminants (PAHs) found in this study, we believe that the species sensitivity values in the BQI calculation need an update to include ecotoxicological sensitivity to account for other types of disturbances than eutrophication only.

Our results suggest that Baltic macrobenthic communities, even in the relatively unpolluted areas, are impacted by anthropogenic contaminants. This impact may need to be accounted for when macrobenthos is used as an indicator of eutrophication. However, long-term contaminant data in sediments are very limited as well as the spatial coverage for macrobenthos community data. To investigate the contaminants in the open sea, the national Swedish status and trend monitoring programme for contaminants in marine sediments was launched in 2003, with 16 offshore stations (12 within the Baltic Sea) and sampling every 5-6 years (Apler and Josefsson, 2016). Unfortunately, six of the 8 stations in the Baltic Proper in this monitoring are located at depths where macrobenthos is absent due to anoxic conditions; hence, the data on the contaminant concentrations have relatively low spatial overlap with the biological data collected within Swedish benthic monitoring programs. Better coordination between the monitoring programs would facilitate the contribution of PAHs (and, probably, other contaminants) in combination with eutrophication to community responses at a larger scale if we are to use changes in the macrobenthic species composition, dynamics of the sentinel amphipod M. affinis, and, ultimately, BQI for environmental quality assessment.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

Author Contributions

AK designed the study and together with CR compiled the data and drafted the manuscript. CR performed all the final statistical analyses. EG provided data on PAHs and contributed to the manuscript writing. All authors contributed to the article and approved the submitted version.

Funding

Metal analyses was funded by The King Carl XVI 50-year Foundation for Science, Technology and Environment to AK. The collection of sediments and PAH analyses were funded via project ReproIND (FORMAS) and the national program for biological effect monitoring of contaminants and INSERT project (Swedish Environmental Protection Agency) to EG. The Swedish national and regional benthic monitoring program is funded by the Swedish Agency for Marine and Water Management (Havs och Vattenmyndigheten) and the county administrative board (Länsstyrelsen).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Matias Ledesma (Stockholm University) for compiling the environmental data from DAS. Jan Albertsson (Umeå University) and Mats Blomqvist (Hafok AB) provided benthic data from the Bothnian Sea. Jonas Gunnarsson is acknowledged for comments on a previous draft.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenvs.2021.624658/full#supplementary-material

Footnotes

^ http://nest.su.se/das

References

Anderson, M. J. (2001). A new method for non-parametric multivariate analysis of variance. Austral. Ecol. 26, 32–46. doi: 10.1111/j.1442-9993.2001.01070.pp.x