Coastal Ecosystem Effects of Increased Summer Temperature and Contamination by the Flame Retardant HBCDD

The combined effects of ocean warming and contaminants on marine ecosystems are poorly understood. In this study, we exposed model ecosystems comprising typical shallow coastal Baltic Sea communities to elevated temperature (+5 ◦C) and the flame retardant hexabromocyclododecane (HBCDD), both singly and in combination, for 13 days. Higher temperatures caused the release of PO4 from the sediment, which in turn stimulated the growth of the cyanobacteria Dolichospermum sp. This in turn led to an increase in the copepod Acartia bifilosa and other indirect effects in the plankton, interpreted as being caused by changes in predation, grazing, and competition. Elevated temperatures also stimulated benthic primary production and increased production of benthic mollusk larvae. Although increased temperature was the dominant driver of effects in these systems, HBCDD also appeared to have some effects, mainly in the zooplankton (both direct and indirect effects) and benthic meiofauna (an interactive effect with temperature). Although the study used model ecosystems, which are an approximation of field conditions, it highlights that interactive ecosystem effects between two stressors are possible and demonstrates the ecological and temporal complexity of such responses. Such unpredictable responses to warming and contaminants are a major challenge for ecosystem management to deal with multistressor situations in the Baltic Sea.


Introduction
Marine ecosystems are exposed to a wide range of anthropogenic stressors, such as eutrophication, global warming, overfishing, and chemical pollution.Substantial research has been carried out on the threats of each of these stressors on marine ecosystems, but the potential effects of these threats in combination are understudied and highly uncertain.There are strong possibilities for unforeseen interactive effects between multiple stressors [1], non-linear responses, feedback loops, and overall erosion of ecosystem resilience [2].During the last decade there has been an increased awareness of the potential interactions between multiple stressors caused by climate change (e.g., [3]) and specifically between climate change and contaminants [4][5][6], but knowledge is lacking, particularly at an ecosystem level.
Increased temperature may alter ecosystem effects of contaminants in a number of ways.Physicochemical properties of the contaminant, and thus contaminants' environmental behavior, bioavailability, and toxicity, may be changed.The rates and efficiency of metabolic processes (including homeostasis, growth, contaminant uptake, and elimination) in organisms may also be altered [4][5][6].Increased temperature may alter biogeochemical cycles, further impacting marine species through changes in, e.g., nutrient availability.Both factors may cause physiological stress but relative sensitivities of species differ, so exposure to these may change the relative fitness of species and indirectly alter the ecosystem or food web composition [7], in turn influencing exposure pathways.Species living at the edge of their range of physiological tolerance, such as many Baltic Sea species, may be particularly susceptible to such interactive effects [4,6].
Global warming and chemical pollution are recognized as two of the major threats to Baltic Sea ecosystems in general [8] and also to Baltic Sea red-listed biotopes [9].In the Baltic Sea, climate models have predicted increased water temperatures as a result of global warming.The exact increase is of course difficult to predict, and will vary spatially and with season, but estimates are generally in the range of 2-5 • C by the end of the century [10][11][12].In all sub-basins, the surface layer is projected to warm more than the deep-water [12].It should, however, be noted that although these models have a good resolution (down to 10-50 m; [11]), these models are at best regional rather than local and shallow coastal bays could be expected to have more extreme temperature fluctuations.The Helsinki Commission (HELCOM) [13] identified temperature as an anthropogenic pressure likely to affect several Baltic Sea "core indicators" such as zooplankton size and abundance and a number of benthic indicators.Temperature is acknowledged by HELCOM to have potentially large impacts on the abundance and distributions of species in Baltic ecosystems, but potential interactions with contaminant pollution are not mentioned.
The Baltic Sea also experiences a number of other environmental stresses.Its limited water exchange with the North Sea and short geological history have resulted in a permanent north-south salinity gradient (3-15 psu) and species-poor flora and fauna of both marine and freshwater origins.Parts of the Baltic coastline are spatially complex, with multiple islands (i.e., archipelagos) or a large number of shallow inlets, lagoons, and shallow soft bottom ecosystems (e.g., [14]).Coastal areas such as these are often areas with high levels of anthropogenic pressures, and several of these types of biotopes have been identified as "vulnerable" or "near-threatened" in a recent evaluation by HELCOM [9].In addition, the Baltic Sea has a large and heavily populated drainage basin which, in combination with the slow water exchange rate, leads to the accumulation of nutrients and contaminants.These include legacy contaminants such as polychlorinated biphenyls (PCBs), dioxins, dichlorodiphenyltrichloroethane (DDT), and metals, but also emerging contaminants such as pharmaceuticals, perfluoroalkyl substances, and flame retardants.Increasing our currently poor understanding of the potential effects of these emerging contaminants is essential for making robust assessments of environmental quality and regulatory decisions.
In the Baltic Sea, HELCOM's list of eleven substances of "specific concern" includes HBCDD (hexabromocyclododecane) and it is one of the substances in the set of HELCOM core indicators for hazardous substances [9], which are used to assess the state of the Baltic marine environment and the EU Marine Strategy Framework Directive (MSFD) qualitative descriptors ("contaminants in environment").HBCDD is a brominated flame retardant (BFR) used widely in building materials, polystyrene, textiles, and electronics.It is used additively, i.e., it is not chemically bound in the matrix, and is thus gradually released from the products to which it has been added or from waste.It may therefore enter the marine environment directly, via atmospheric fallout or runoff from land.It is classified as a Persistent, Bioaccumulative, and Toxic (PBT) substance and a "Substance of very high concern" [15].In 2013, HBCDD was added to Annex A of the Stockholm Convention on Persistent Organic Pollutants (POPs), a global treaty to protect human health and the environment from environmentally persistent chemicals by for example reducing production, import and export, use, and release to the environment.HBCDD is also listed in Annex XIV list of the European Union Regulation REACH (Registration, Evaluation, Authorisation and Restriction of Chemicals).
The main source of HBCDD to the Baltic is "leakage from products (e.g., insulation, electrical products, and textiles) via rivers, atmosphere, or wastewater treatment plants" [9].HBCDD concentrations are measured in a range of species as part of the Swedish national monitoring program and temporal trends vary spatially and temporally and with the organism measured [16].For example: concentrations in herring muscle and cod liver are generally higher in the Baltic Proper than on the Swedish West coast; concentrations in guillemot eggs have been decreasing over the past decade from their peak in the mid-2000s, but steadily increasing since 1980 in cod liver from SE Gotland.Concentrations may also vary widely between samples from the same species, site, and time point [16].
The effects of HBCDD used in its classification as a PBT were determined using standard ecotoxicological tests (for a summary see [15]).For aquatic species, reduced size and reproductive output was seen in Daphnia magna (21-day exposure, Lowest Observed Effect Concentration (LOEC) 5.6 µg L −1 ), and reduced growth rates were seen in three microalgae (72-h 50% Effect Concentration (EC 50 ) c. 40-50 µg L −1 ).No effects were seen in growth of the amphipod Hyalella azteca at concentrations up to 1000 mg•kgdw −1 but numbers of the worm Lumbriculus variegatus were reduced at 50 mg•kgdw −1 .No measureable effect was seen on eggs, larvae, or adult rainbow trout at concentrations up to 6.8 µg L −1 .Assessments are often complicated by the low solubility of HBCDD and the different solubilities and toxicities of the different diastereomers present in the technical mixture of HBCDD.A few experiments have also been done using Baltic species, applying HBCDD in a more realistic manner, associated to suspended organic matter; Smolarz and Berger reported cytotoxic and genotoxic effects (significant increases in nuclear and nucleolar abnormalities and frequency of dead cells) of HBCDD in the bivalve Limecola balthica [17].Pirzadeh et al. suggested that copepods may be a particularly sensitive part of the Baltic plankton [18].Bradshaw et al. used an 8-month mesocosm study to demonstrate that HBCDD caused effects on benthic population structure and ecosystem function and possible changes in benthic-pelagic coupling of nutrients with knock-on effects in the plankton [19].
The main aim of this study was to determine if increased temperature could modify the ecological response of a shallow coastal Baltic ecosystem to a contaminant, HBCDD.This was done with an experiment using model ecosystems ("cosms") to approximate the natural environment whilst maintaining control over experimental conditions.The cosms were exposed to the two stressors either singly or in combination.The study found that increased temperature had the largest effect, mainly through stimulating the release of P from the sediment into the water column.HBCDD had a less dominant, though detectable effect in some parts of the ecosystem.

General Experimental Design
This paper presents the results of a two-week experiment using 5 L model ecosystems (cosms) constructed from water, sediment, and naturally occurring species at relevant densities, from a coastal bay in the Baltic Proper.A crossed two-factorial design was used, with the factors of temperature (ambient and warm (+5 • C)) and contaminant (with or without HBCDD).The experiment was designed to be as environmentally relevant as possible, within the constraints of working with cosms.Thus, a concentration of HBCDD was used that may be found in the environment and at which effects have previously been seen in similar systems.The exposure pathway was spiked dead phytoplankton, such as might occur after a spring bloom settles out to the seabed.In addition, outdoor systems were used that reflected natural summer temperatures in a shallow coastal bay, as well as a realistic increased temperature (+5 • C) that might be expected with global warming in the Baltic region (see Section 1).An array of relevant structural and functional endpoints was measured in both the pelagic and benthic parts of the ecosystem.

Collection of Components of the Experimental System
Sediment, water, and organisms (the benthic bivalves L. balthica (previously Macoma balthica) and Cerastoderma glaucum, Hydrobiidae snails) were collected between March and May 2007 from c. 1 m water depth in a bay near Stockholm University's Askö Laboratory (58 • 49.4 N, 17 • 38.2 E) (see Figure 1 for the experimental time line).The sediment was sieved on a 1 mm mesh to remove all macrofauna, which was saved in separate sediment and kept aerated, cool, and dark until the start of the experiment.The sieved sediment and water was transferred to a large plastic tank (100 × 100 × 50 cm) which was placed outdoors to allow warming of the sediment and hatching of plankton from the sediment.After one week, the water was visibly cloudy with plankton.Extra zooplankton was also collected just before the start of the experiment by passing 450 L of water from the same bay through a 25 µm mesh plankton net.This sample was dominated by rotifers, with some cladocerans and copepods.the experiment.The sieved sediment and water was transferred to a large plastic tank (100 × 100 × 50 cm) which was placed outdoors to allow warming of the sediment and hatching of plankton from the sediment.After one week, the water was visibly cloudy with plankton.Extra zooplankton was also collected just before the start of the experiment by passing 450 L of water from the same bay through a 25 μm mesh plankton net.This sample was dominated by rotifers, with some cladocerans and copepods.

Spiking of Organic Material with HBCDD
Phytoplankton was collected in April 2007, partly during a monitoring program in the area, and partly by trawling from a small boat in the immediate vicinity of the Askö Laboratory, in both cases using a 90 μm plankton net.These samples were concentrated and stored at −20 °C until later spiking with HBCDD.At that time, they were thawed out, thoroughly mixed, and divided into two stocks, one that was spiked with a solution of HBCDD in an acetone carrier in an amount calculated to achieve a final nominal amount per cosm of 78.75 μg, the other unspiked but with an equal amount of acetone.These two batches were kept cold (4-5 °C) on magnetic stirrers for 1 week before being added in equal amounts to each of the experimental containers.

Experimental Set-Up and Cosm Construction
The experiment was conducted outdoors under an open-sided transparent plastic shelter that allowed sunlight through while preventing rain or debris entering the cosms.The cosms were randomly distributed between 4 troughs.Two of these troughs had cool water flow-through pumped directly from the sea.The other two were maintained at c. 5 °C above ambient using aquaria water heaters (2 × 300 W heaters per trough).Two water pumps ensured thorough mixing of the warmed water in the troughs.Photographs of the set-up can be found in Figure S1.Distilled water was added to the cosms five times during the experiment to compensate for evaporation and maintain constant salinity.
During the experiment, the temperature of the ambient treatments varied between c. 12-20 °C, depending on the time of day and the weather.The warm treatments varied between c. 17-26 °C.The maximum daily variation was 7.6 °C and 9.1 °C in the ambient and warm treatments,

Spiking of Organic Material with HBCDD
Phytoplankton was collected in April 2007, partly during a monitoring program in the area, and partly by trawling from a small boat in the immediate vicinity of the Askö Laboratory, in both cases using a 90 µm plankton net.These samples were concentrated and stored at −20 • C until later spiking with HBCDD.At that time, they were thawed out, thoroughly mixed, and divided into two stocks, one that was spiked with a solution of HBCDD in an acetone carrier in an amount calculated to achieve a final nominal amount per cosm of 78.75 µg, the other unspiked but with an equal amount of acetone.These two batches were kept cold (4-5 • C) on magnetic stirrers for 1 week before being added in equal amounts to each of the experimental containers.

Experimental Set-Up and Cosm Construction
The experiment was conducted outdoors under an open-sided transparent plastic shelter that allowed sunlight through while preventing rain or debris entering the cosms.The cosms were randomly distributed between 4 troughs.Two of these troughs had cool water flow-through pumped directly from the sea.The other two were maintained at c. 5 • C above ambient using aquaria water heaters (2 × 300 W heaters per trough).Two water pumps ensured thorough mixing of the warmed water in the troughs.Photographs of the set-up can be found in Figure S1.Distilled water was added to the cosms five times during the experiment to compensate for evaporation and maintain constant salinity.
During the experiment, the temperature of the ambient treatments varied between c. 12-20 • C, depending on the time of day and the weather.The warm treatments varied between c. 17-26 • C. The maximum daily variation was 7.6 • C and 9.1 • C in the ambient and warm treatments, respectively.The salinity in the cosms was c. 7.5 psu and oxygen water concentrations were consistently >8 mg•L −1 .
The cosms were prepared as follows.Water from the hatching tank was drained off into holding containers.The sediment was then thoroughly mixed and divided equally (610 mL wet sediment each) into 58 separate 5 L buckets (cosms).These were factory-new, made of polypropylene, and approved for food use.Two of these cosms were taken immediately for analysis of CNP and HBCDD in the water and CNP content of the sediment (see Section 2.5.4).Hatching tank water was then equally distributed (3.75 L each) between the remaining 56 buckets, which were then placed in the troughs and allowed to stand for 3 days.After this time the sediment had settled out to form a layer 2 cm thick.The number of buckets allowed for 4-6 replicates of each of the four treatments (ambient control (A−), warm control (W−), ambient with HBCDD (A+), and warm with HBCDD (W+)), and for destructive sampling at three time points (1, 6, 13 days after HBCDD addition; t1, t6 and t13).
The previously collected organisms were sieved out of their holding sediment and four L. balthica (11.5 ± 1.2 mm), two C. glaucum (shell length c. 16 mm and c. 8.5 mm), and 10 Hydrobiidae snails were added to each of the treatments.This represented natural field densities of 187, 93, and 468 individuals per m 2 , respectively.The L. balthica burrowed quickly into the sediment.
Finally, 4.5 mL HBCDD-spiked plankton suspension was added to half of the cosms (A+ and W+ treatments) to give a nominal amount in each cosm of 78.75 µg HBCDD.The same amount of unspiked suspension was added to the A− and W− cosms to ensure equal amounts of extra organic material were added to all cosms (63 mg dw or 21.4 mg TOC (total organic carbon) per cosm).
The experiment was run for 13 days.During this time the weather was in general warm and sunny.

Sampling
Non-destructive sampling was used for measurements of whole system metabolism, and pelagic primary production.Destructive sampling (i.e., disassembling the contents of the cosms) was done for the rest of the endpoints (pelagic chlorophyll a, phyto-and zooplankton community composition), and sediment fauna community composition (macrofauna, meiofauna), CNP, and HBCDD analysis.This was done three times during the experiment, at 1, 6, and 13 days after the HBCDD was added (t1, t6, t13) (see Figure 1).

Whole System Metabolism
The whole system metabolism was determined five times during the experiment from O 2 measurements using a WTW Multi 350i (WTW GmbH, Weilheim, Germany), though only data from t1, t6, and t13 is used here.Each time, O 2 concentrations in the water were measured twice in the dark (starting c. 22:00) and twice in the light (starting c. 10:30).The difference between the dark measurements was used to calculate the whole system respiration (R, in mg O 2 h −1 ) and the difference between the light measurements was used to calculate the net primary production (NPP, in mg O 2 h −1 ).Gross primary production (GPP, in mg O 2 h −1 ) was calculated as NPP-R and the production/respiration ratio (P/R) as GPP/R.

Pelagic Endpoints
Pelagic primary production was measured five times during the experiment (t1, t3, t6, t10, t13), though only data from t1, t6, and t13 is used in this paper.Two 10 mL water samples were taken (in the morning) from each microcosm and placed in two 20 mL scintillation vials.One hundred µL 14 C was added as buffered NaH 14 CO 3 (specific activity 185 kBq mL −1 ) and the samples were incubated at ambient temperature for a few hours, one of the samples in the light and the other in the dark (covered in black plastic).Three drops of 10% HCl was then added to stop primary production, excess 14 C bubbled off, and 10 mL Lumagel Safe scintillation cocktail (Perkin Elmer, Waltham, MA, USA) was added. 14C activity in the samples was measured in a 1214 Rackbeta scintillation counter (LKB Wallac, Turku, Finland) and decays per minute (dpm) converted to mg 14 C L −1 h −1 .Net primary production was calculated as the difference between 14 C uptake in the light minus that in the dark.
All other pelagic endpoints were measured at t1, t6, and t13.In order to measure pelagic chlorophyll a, phyto-and zooplankton community composition, different fractions of the overlying water were obtained.Firstly, the overlying water was siphoned off into a bucket from which a well-mixed water sample of known volume was filtered through a 90 µm sieve.The fraction retained on the sieve (zooplankton) and the fraction that had passed through the sieve (largely phytoplankton) were fixed separately in Lugol solution and saved in the dark until quantification and identification to the lowest possible taxonomic level.A further 1-2 L was filtered through a 90 µm sieve to remove the zooplankton and the filtrate was further filtered on 0.7 µm cellulose-nitrate filters and stored at −20 • C until the analysis of chlorophyll a, a measure of phytoplankton production/biomass.At that time, frozen filters were thawed and extracted in 5 mL acetone; the extract was analyzed with a U-2000 spectrophotometer (Hitachi, Tokyo, Japan) according to [21].

Benthic Endpoints
Directly after removing the water for the plankton and pelagic chlorophyll a measurements, a 1 mL (4 mm diameter) cut-off plastic pipette was then used to take 4 small sediment cores from the sediment.The top 2 mm of each mini-core (c. 25 mm 3 ) was saved and all four pooled to a single sample in a Winkler flask.20 mL filtered seawater and 100 µL buffered NaH 14 CO 3 solution (specific activity 185 kBq mL −1 ) was then added to each.Double samples of this type were taken from each cosm for light and dark incubations and the flasks were incubated for several hours in the troughs at the temperature from which they had come.Twenty-five mL filtered seawater were then added and 10 mL samples of well-shaken suspension were taken.Three drops 10% HCl were added to stop all production, excess 14 C bubbled off, and 10 mL Lumagel Safe (Perkin Elmer, Waltham, MA, USA) was added. 14C activity was measured in a 1214 Rackbeta scintillation counter (LKB Wallac, Turku, Finland) and dpm counts per sample were converted to mg 14 Macrofauna were sieved out of the remaining sediment on 1 mm and 5 mm sieves.L. balthica and C. glaucum shell lengths were measured and all species were counted, allowed to clear their guts for several hours, and were then frozen at −20 • C for later HBCDD analysis.The number of dead individuals was noted.At t13, 100 mL of the top 1 cm sediment was scraped off using a spatula and preserved with 40 mL 95% ethanol for identification and quantification of meiofauna.This was done by sieving the sediment on a 63 µm sieve and quantifying the organisms to the lowest possible taxon level.

CNP and HBCDD Sampling and Analysis
Water samples for the analysis of C, N, and P were taken at t1, t6, and t13; C was analyzed in unfiltered water and N and P in water filtered through a 0.45 µm filter.Nitrogen (NO 2 − , NO 3 − , and NH 4 + ) was analyzed with a Leco-CHNS-932 analyzer (Isomass Scientific Inc., Calgary, AB, Canada) (EDTA and acetanilide as standards) and phosphorus (PO 4 3− ) with a segmented flow system (ALPKEM Flow Solution IV) following combustion at 500 • C and acidic persulfate oxidation, at the accredited laboratory at the Department of Ecology, Environment, and Plant Sciences, Stockholm University, Sweden.
To determine the fate of the HBCDD in the experimental systems, samples of the spiked organic material were taken, as well as water and sediment samples at t1 and t13, from both A+ and W+ treatments.The soft tissues of L. balthica and C. glaucum collected at t13 (see above) were also analysed.All samples were frozen at −20 • C until analysis using gas chromatography-mass spectrometry (GC-MS/MS, Agilent Technologies, Santa Clara, CA, USA) and liquid chromatography-mass spectrometry (LC-MS/MS, Shimadzu Corp., Kyoto, Japan) at the Department of Environmental Chemistry, Stockholm University.The relative amounts of the three diastereomers (α, β, γ) were also determined.Full details of these chemical analyses are available in Bradshaw et al. [20] (Supplementary Information).

Data Handling and Statistics
Single and interactive effects of warming (ambient and warm), HBCDD addition (with HBCDD and without HBCDD), abiotic variables, and phyto-and zooplankton species abundances, were assessed using constrained analysis of principle coordinates (CAP), based on Bray-Curtis dissimilarities.Prior to testing the effects of temperature and HBCDD addition and abiotic variables on species abundances of both phyto-, and zooplankton, species matrices were fourth root transformed to account for homoscedasticity and to downweigh the importance of taxa with very high abundances.Additionally, the abiotic variables were tested for autocorrelation between predictor variables and factors and no autocorrelation was found.Furthermore, the effect of the blocks (i.e., troughs) was tested, but since it did not have an effect, it was removed from the analyses.All multivariate analyses were performed in R version 3.2.[22].
The model for phytoplankton included both factors of warming and HBCDD addition as single and interactive effects, as well as all abiotic variables chlorophyll a, phosphate (PO 4 ), total nitrogen, P/R ratios, and pelagic primary production, with the holding times (t6 and t13) constant.Significance of the model, terms, and axis were assessed using a permutation procedure with 999 permutations.The minimal adequate model (containing only significant parameters α = 0.05) was identified through step-wise deletion of non-significant parameters [23].The analyses were performed using the "capscale" function from the {vegan} R package [24].A similarity percentage test (SIMPER) was used to identify which phytoplankton species drove the community-level effects, using significant factor(s).For the zooplankton species matrix, both the abiotic variables (same as above) as well as all phytoplankton species were used as predictors.The analysis was then performed the same way as described for the phytoplankton species abundance.
Three-way ANOVA (factors: time, temperature, HBCDD) was used to test for: differences in pelagic and benthic primary production (standardized to average A− at each time point to remove the effect of differences in light intensities during incubations at the different time points), total N, and PO 4 − concentrations in the water.Homogeneity of variance was tested using Cochran's test and non-homogeneous data was Box-Cox transformed where necessary.The analyses were performed in Statistica version 12.5 (Dell Inc., Tulsa, OK, USA).Since too few benthic community data were available for a meaningful benthic CAP analysis, and several benthic endpoints were only measured at one time point, benthic parameters were analyzed using correlation analysis and ANOVA.Two-way ANOVA (factors: temperature, HBCDD) was used to test for differences in: abundance of each meiofauna taxa at t13; abundance of Hydrobiidae at t6; L. balthica biomass per cosm at t6 and t13.Homogeneity of variance was tested using Cochran's test and non-homogeneous data were square-root transformed where necessary.
To investigate potential connections between the benthic and pelagic components of the ecosystem, two multiple correlations were performed (t6 and t13 separately) on the following data: abundances of the 5 phytoplankton and 6 zooplankton groups identified by SIMPER (see above) to cumulatively contribute to 70% of the dissimilarity between treatments in the pelagic analysis; pelagic and benthic primary production; chlorophyll a, total N, and PO 4 concentrations.For t6, the number of Hydrobiidae was also available.For the t13 data correlation, the following variables were also available: total number of meiofauna, numbers of Turbellaria and Microarthridion littorale, and total biomass of L. balthica.Although the focus of this analysis was benthic-pelagic linkages, it was also used to help explain the relationships between plankton species.The significance of the correlation coefficients r was tested using a two-tailed t test.The analyses were performed in Statistica v.12.5.
A significance level of 0.05 was used in all statistical tests.

HBCDD Partitioning and Exposure of the Model Ecosystem
After only one day, >90% of the HBCDD was detected in the sediment [20] and this increased to >99% by t13.Concentrations in the water were very low (Table 1), with >80% of the HBCDD in the water being present in the particulate fraction at t1, decreasing to 66.5%-78.7%by the end of the experiment, depending on the treatment (Table 1).Concentrations in both water fractions decreased throughout the experiment, with the warmer treatment (W+) having lower HBCDD concentrations than the ambient treatment (W−).No significant effect of either temperature or time was seen on sediment concentrations, although HBCDD concentrations tended to be lower at higher temperatures.There was no obvious difference between treatments for HBCDD body concentrations in L. balthica at the end of the experiment.For full details of partitioning, including that of HBCDD diastereomers, see [20].
Table 1.Measured concentrations of hexabromocyclododecane (HBCDD) in the various compartments and treatments one day after the start and at the end of the experiment.For more details, see [20].

Pelagic Part of the Ecosystem
Pelagic primary production was consistently higher in the warm treatments (W− and W+) (3-way ANOVA, significant main effect of temperature, F 1,60 , p < 0.0001), but there was no effect of HBCDD (Figure 2, Table S12).
water being present in the particulate fraction at t1, decreasing to 66.5%-78.7%by the end of the experiment, depending on the treatment (Table 1).Concentrations in both water fractions decreased throughout the experiment, with the warmer treatment (W+) having lower HBCDD concentrations than the ambient treatment (W−).No significant effect of either temperature or time was seen on sediment concentrations, although HBCDD concentrations tended to be lower at higher temperatures.There was no obvious difference between treatments for HBCDD body concentrations in L. balthica at the end of the experiment.For full details of partitioning, including that of HBCDD diastereomers, see [20].
Table 1.Measured concentrations of hexabromocyclododecane (HBCDD) in the various compartments and treatments one day after the start and at the end of the experiment.For more details, see [20].

Pelagic Part of the Ecosystem
Pelagic primary production was consistently higher in the warm treatments (W− and W+) (3-way ANOVA, significant main effect of temperature, F1,60, p < 0.0001), but there was no effect of HBCDD (Figure 2, Table S12).The total nitrogen concentration (sum of NO2 − , NO3 − and NH4 + ) decreased in all four treatments at t6 compared to t1 (3-way ANOVA, significant main effect of time, F2,47, p = 0.0001).The decrease was more pronounced in the warm treatments compared to ambient treatments, though this was not statistically significant (3-way ANOVA, interaction of time and temperature, F2,47, p = 0.072) (Figure 3a, The total nitrogen concentration (sum of NO 2 − , NO 3 − and NH 4 + ) decreased in all four treatments at t6 compared to t1 (3-way ANOVA, significant main effect of time, F 2,47 , p = 0.0001).The decrease was more pronounced in the warm treatments compared to ambient treatments, though this was not statistically significant (3-way ANOVA, interaction of time and temperature, F 2,47 , p = 0.072) (Figure 3a, Table S13).At t13, the total N concentrations further decreased in the ambient temperature treatments, while a slight increase could be observed in the warm treatments.Phosphate concentration slightly increased at t6 in the ambient treatments, while the increase in the warm treatments was significantly more pronounced (Figure 3b), as was the decrease in PO 4 by t13 (3-way ANOVA, interaction of time and temperature, F 2,47 , p < 0.0001; Table S14).
J. Mar.Sci.Eng.2017, 5, 18 9 of 21 Table S13).At t13, the total N concentrations further decreased in the ambient temperature treatments, while a slight increase could be observed in the warm treatments.Phosphate concentration slightly increased at t6 in the ambient treatments, while the increase in the warm treatments was significantly more pronounced (Figure 3b), as was the decrease in PO4 by t13 (3-way ANOVA, interaction of time and temperature, F2,47, p < 0.0001; Table S14).

Meiofauna
Total meiofauna (Figure 9a) abundances at the end of the experiment were dominated by the abundance of nematodes.Although total nematode numbers were lower in the warm treatments than in the ambient ones (Figure 9b), this was not significant (2-way ANOVA, F1,15, p = 0.12; Table Figure 8. Benthic primary production, standardized to the A− treatment to enable valid comparisons between time points (mean ± SE).

Meiofauna
Total meiofauna (Figure 9a) abundances at the end of the experiment were dominated by the abundance of nematodes.Although total nematode numbers were lower in the warm treatments than in the ambient ones (Figure 9b), this was not significant (2-way ANOVA, F 1,15 , p = 0.12; Table S16) and there was no significant effect of HBCDD (2-way ANOVA, F 1,15 , p = 0.88), nor any significant interaction between the two factors (2-way ANOVA, F 1,15 , p = 0.44).There was no overall effect of HBCDD or temperature on Turbellaria or M. littorale, but there was a significant interaction between these stressors (2-way ANOVA, F 1,15 , p < 0.05) for both species (Figure 9c,d, Tables S17 and S18).

Benthic-Pelagic Coupling
Benthic and pelagic variables seemed in general to be only weakly correlated (Tables S10 and  S11).We looked especially carefully at the potential influence of L. balthica biomass, hydrobid numbers, and benthic primary production in the systems and at the two zooplankton taxa that are the larvae of benthic species (gastropod and bivalve larvae).

Benthic-Pelagic Coupling
Benthic and pelagic variables seemed in general to be only weakly correlated (Tables S10 and S11).We looked especially carefully at the potential influence of L. balthica biomass, hydrobid numbers, and benthic primary production in the systems and at the two zooplankton taxa that are the larvae of benthic species (gastropod and bivalve larvae).
Benthic primary production was significantly correlated to PO 4 , pelagic primary production, and Dolichospermum and Melosira sp.abundance at t6 (Table 2).At t13 it was not correlated to any of the variables.L. balthica biomass was only significantly correlated (r = 0.48, p = 0.05) with one of the other 19 variables, Melosira sp., (Table 2 and Table S11) a result that is most likely not causative.In addition, L. balthica biomass per microcosm showed no significant difference between treatments, either at t6 or t13 (2-way ANOVA, F 1,15 , p > 0.05; Tables S19 and S20), though there was a tendency for lower biomass in HBCDD treatments at t13.The number of Hydrobiidae increased from 10 at the start of the experiment to 15-20 by t6.There was a tendency for more Hydrobiidae in the HBCDD treatments at t6 (Figure S2), but this was not significant (Table S21).The number of Hydrobiidae at t6 only correlated significantly (r = 0.45, p = 0.05) with total N concentrations in the water (Table 2).Abundances of gastropod larvae were significantly positively correlated with PO 4 at both time points, positively correlated with pelagic chlorophyll a and Dolichospermum heterocysts at t6, but negatively at t13.There was also a significant correlation (r = 0.62, p = 0.008) between the number of gastropod larvae and the total number of meiofauna, but this result is probably coincidental.

Benthic-Pelagic Coupling: Nutrient Release from Sediment and Effects on the Plankton
In previous work, with larger cosms run over eight months [19], we showed that HBCDD negatively affected the biomass of "large" (>8 mm) L. balthica in the same size class as in this experiment.This in turn affected benthic-pelagic coupling, reducing the recirculation of nutrients from the sediment to the water.In the current study, no changes in L. balthica biomass were observed, probably due to the short duration of this study.Trends in benthic and pelagic parameters were not linked to differences in L. balthica biomass between treatments, but to overall effects of temperature, and to a lesser degree HBCDD (see Section 4.3).
It is well-known from freshwater systems that increased temperature causes the release of PO 4 from sediments [25][26][27].The few studies that have been performed in marine ecosystems suggest the same (e.g., [28][29][30]).This PO 4 release is thought to occur due to increased microbial activity increasing the decomposition of organic matter at the sediment surface [31][32][33] and by consuming oxygen at the sediment-water interface [25], leading to a redox-mediated release of P [27].Studies of anoxic marine benthic ecosystems also confirm that low oxygen conditions facilitate fluxes of nutrients from the sediment to the water [28,30,34].Although overall oxygen levels in the water column remained high in this experiment, it is possible that local anoxia/hypoxia at the sediment surface occurred in the warm treatments.We did not measure benthic bacterial activity specifically in this study, but measurements of benthic primary productivity support the idea that benthic processes were higher at higher temperatures, particularly at t6.By t13, P levels in the sediment were probably quite low, leading to the reduction in benthic primary productivity at this time point.It should be noted that by adding organic matter to the cosms at t0, we also added a source of P, which potentially could have increased water PO 4 concentrations.However, we estimate that this addition was negligible compared to the P pool in the sediment.Calculations using measured TOC and C:P ratios in the added suspended organic matter and the sediment suggest that the organic matter provided c. 0.35 mgP per cosm compared c. 1200 mgP in the sediment.
This nutrient regime in the warm treatments therefore has similarities with the Baltic Sea nutrient regime after the spring bloom, where N and P are depleted and the sedimentation of organic material leads to increased microbial activity, leading to anoxia, which in turn results in an increased P release by the anoxic sediment [35].While our microcosms did not turn anoxic, filamentous nitrogen-fixing cyanobacteria Dolichospermum sp.increased from relatively low numbers at t6 to the most abundant species in the warm treatments at t13.As P is the primary limiting nutrient for N-fixation, the increase in phosphorus concentrations, through PO 4 release from the sediment, triggered a rapid increase in Dolichospermum sp.abundances [36] and in the number of nitrogen-fixing heterocysts present.While cyanobacteria blooms are generally considered to be harmful, recent studies indicate that the N fixed and released as NH 4 by cyanobacteria might actually be beneficial for other primary producers (i.e., [35]).Stable isotope analysis and nanoscale secondary ion mass spectrometry (nanoSIMS) have shown that N derived from cyanobacterial N-fixation can be taken up by diatoms in the form of NH 4 .Furthermore, recent studies indicate that cyanobacteria can be a direct food source for mesoand microzooplankton.Stable isotope analyses show that cyanobacteria are a major food source of copepods species such as Acartia spp.[37].This could explain the positive relationship of Acartia with increasing chlorophyll a concentrations, as this increase was most likely driven by Dolichospermum filaments and therefore represents an increase in food items.The increase in rotifers was most likely due to an increase in abundance of a potential food species, Melosira sp., in ambient temperature treatments, or could be explained by predator-prey release.Copepods are known to not only graze on phytoplankton but also prey on other zooplankton species including rotifers [38].Due to an increase in PO 4 and a consecutive cyanobacteria increase in the warm treatments, Acartia increased in the warm treatments, while abundances in the ambient treatments stayed relatively low, which in turn resulted in lower predation pressure of rotifers in the ambient treatments.
The overwhelming effect of temperature-induced PO 4 release, and the resulting changes in the phytoplankton and pelagic primary production, was probably responsible for the significant correlations observed at t6 between the benthic primary production and PO 4 , chlorophyll a, and some phytoplankton species.It is unlikely that these correlations are a cause-effect relationship; rather, they are probably reflecting the common increase in rates of ecological processes in the warm treatments at that time point (Figure 10).
predation pressure of rotifers in the ambient treatments.
The overwhelming effect of temperature-induced PO4 release, and the resulting changes in the phytoplankton and pelagic primary production, was probably responsible for the significant correlations observed at t6 between the benthic primary production and PO4, chlorophyll a, and some phytoplankton species.It is unlikely that these correlations are a cause-effect relationship; rather, they are probably reflecting the common increase in rates of ecological processes in the warm treatments at that time point (Figure 10).

Benthic-Pelagic Coupling: Benthic Species' Larvae in the Zooplankton
Benthic species were also present in the water column as larvae; bivalve and gastropod larvae were identified as two of the planktonic taxa that were responsible for differences in the

Benthic-Pelagic Coupling: Benthic Species' Larvae in the Zooplankton
Benthic species were also present in the water column as larvae; bivalve and gastropod larvae were identified as two of the planktonic taxa that were responsible for differences in the zooplankton community structure.Although bivalve larvae may survive in the water column for many weeks, their general in numbers in the warm treatments throughout the experiment suggests that these were released from adults in the cosms, rather than being present in the water column from the start.In this case, they were most probably L. balthica (10 adults present in each cosm), or possibly C. glaucum (2 per cosm), both of which are known to spawn at the time of year that the experiment was performed (e.g., [39,40]).
There were no clear effects on adult L. balthica or C. glaucum in terms of mortality or size, so the clear effects on their larvae must be either due to direct effects on adult spawning or larval survival, or indirect effects on the larvae related to foodweb dynamics.Drent showed that L. balthica larvae grew faster and developed more quickly at higher temperatures (experimental range 10-20 • C) [41], and several other studies have shown that spawning is triggered by the temperature rising to a critical level [42,43], albeit at lower temperatures than in this experiment.Benthic primary production at t6 was also positively correlated with bivalve larval abundance at t13; adult L. balthica may have been able to spawn more successfully in the warmer treatments by t13 due to increased food levels on the sediment surface.L. balthica larvae are thought to feed on small unicellular phytoplankton, bacteria, flagellates, and organic matter [44][45][46][47], and will likely compete for these food items with other small zooplankton, such as other larvae and rotifers.They are also likely to be preyed on by larger zooplankton, such as copepods [38].Relationships between bivalve larvae and these other taxa could therefore be expected.However, multiple correlations of bivalve and gastropod larval abundances with all possible prey species, competitors, and predators revealed no significant correlations with the abundances of any of these species.The possibility of indirect effects should, however, not be ruled out.It is difficult in studies such as these to capture the fine scale temporal changes with a limited number of "snapshot" sampling times.Processes or events that have happened in between sampling points may well be the most important in determining the effects that are subsequently measured.
Most of the common gastropod species present in shallow Baltic bays (e.g., Hydrobia ventrosa, H. neglecta, Potamopyrgus antipodarum, Theodoxus fluviatilis, Bithynia tentaculata, Radix spp.) either produce live young, or lay eggs that hatch into juveniles.The only common species that has a larval stage, hatching from eggs laid on the seabed, is Hydrobia ulvae [48], which was most probably the species introduced to the cosms ("Hydrobiidae") at the start of the experiment.Adult Hydrobiidae increased in number from 10 at the start to c. 25 at t6 and >35 (not quantified) at t13, with no significant differences between treatments at t6.However, the larval abundance showed a clear treatment difference, with warm treatments having higher numbers at t6 and ambient treatments having higher numbers at t13.Given that the time to hatching takes at least 5 days, it seems probable that the larvae seen in this experiment hatched from eggs already present in the sediment at the start, a process that was stimulated (or occurred earlier) in the warm treatments.Fish and Fish [49] showed that H. ulvae from a Welsh had a clear optimum temperature with respect to embryonic development and survival as well as the time taken for eggs to hatch, as is the case for most invertebrates [3].
H. ulvae are generalist feeders that graze on organic matter and biofilm on the sediment surface [50].They also consume more biomass at higher temperatures [51].They might be expected to take up a considerable amount of HBCDD from the spiked algal material added to the cosms, particularly at increased temperature, but if this was the case, it did not seem to have adversely affected adult survival in any way.One might also expect to see a correlation between adult or larval hydrobid numbers and benthic primary production (i.e., food supply), but this was not the case.Again, this may due to comparing these parameters at the same time point, whereas in reality there would be a lag between the increased food supply and the increase in population size.

HBCDD Effects
Increased temperature was the factor that had the largest effect on these systems.However, HBCDD was also partly responsible for changes in the zooplankton community: higher numbers of gastropod larvae, rotifers, bivalve larvae, and A. bifilosa at t6, and higher numbers of the latter two at t13; fewer nauplii, rotifers, and Eurytemora affinis at t13 (and rotifers at t6).The decrease in numbers of some taxa may be direct negative effects of the contaminant, but the increase in numbers of some taxa in the presence of HBCDD are likely indirect effects caused by the release of these taxa from predation pressure or competition, or the presence of more food species in these treatments.However, exploration of the potential correlations between competitors, food species, and predators, both within and between time points, revealed little.This is in contrast to Pirzadeh et al. [18], who found copepods (all life stages) to be the most sensitive to HBCDD, resulting in increases in their prey taxa, rotifers, in HBCDD treatments.However, this apparent discrepancy may be explained by the time and frequency of sampling, HBCDD dosage, and other experimental factors.These types of experimental systems are highly dynamic and the temporal resolution of sampling is not sufficient to capture the fine-scale and complex temporal changes in ecosystem structure and function that occur.
It is surprising that more effects of HBCDD were not seen in the benthos, since this is where the majority of the contaminant was present (see Section 3.1).Benthic primary production, L. balthica biomass, and the abundance of most meiofauna species was unaffected by HBCDD.The harpacticoid copepod M. littorale feeds preferentially on diatom biofilm [52] and fresh organic matter settling on the sediment surface, e.g., after phytoplankton blooms [53,54].In this experiment, it is therefore likely that it ingested a considerable amount of HBCDD as the spiked algae settled out quickly on to the sediment.However, HBCDD negatively affected M. littorale's abundance only in the ambient treatments; in warm treatments, abundances were higher in the presence of HBCDD.It is possible that this is due to release from predation by Turbellaria, whose abundances showed approximately the opposite trend to M. littorale.Turbellaria are hard to identify and little is known about their ecology [55], but some species are known to predators on other meiofauna.Turbellaria were least abundant in the treatments A− and W+ where L. balthica were the largest (although these size differences were non-significant).L. balthica's deposit-feeding activity at the sediment surface is known to disturb the feeding activity of some meiofauna groups [56] and this may have been the case for Turbellaria in this study, although there was no correlation between these variables.

Conclusions
This model ecosystem study experimentally demonstrated the complex nature of shallow Baltic ecosystem responses to a combination of two stressors; the increase in temperature and the presence of a contaminant (HBCDD).Increased temperature was the dominant driver of ecosystem changes, mainly through stimulating the release of PO 4 from the sediment to the water column, resulting in a cyanobacteria bloom and knock-on effects in the plankton community.Higher temperature also enhanced benthic primary production and the production of benthic mollusk larvae.HBCDD effects were detectable in the zooplankton and meiofauna but were more complex, in some cases being different depending on the temperature (i.e., an interactive effect; Figure 10).
Both stressors were applied at levels that are not unreasonable to expect in nature.Climate models for the Baltic Sea predict warming in the range of 2-5 • C by the end of the century, and this is also within the natural fluctuation of a shallow bay in the summer.HBCDD was added at a realistic concentration and applied to the system in an ecologically relevant way (associated with suspended organic matter).Although the study was based on model ecosystems, which are an approximation of field conditions, it has shown that interactive ecosystem effects between two stressors are indeed possible in shallow Baltic ecosystems, as well as demonstrating the ecological and temporal complexity of such responses.This is only a first step to increasing our understanding of such complex interactions, using only two stressors.In reality, ecosystems may be exposed to multiple chemicals, temperature stress, salinity stress, overfishing, and eutrophication.A major challenge for future management of Baltic Sea ecosystems is how to ensure "good ecological status" [13] in such multistressor situations.

Supplementary Materials:
The following are available online at www.mdpi.com/2077-1312/5/2/18/s1, Figure S1: Photographs of the experimental set-up, Figure S2: Number of Hydrobiidae per cosm at t6, Tables S1-S9: Results of the similarity percentage (SIMPER) analyses for phytoplankton and zooplankton communities at different time points, Tables S10 and S11: Full results for multiple correlations at t6 and t13, Tables S12-S22: Full results of ANOVA analyses.

Figure 2 .
Figure 2. Pelagic primary production over time standardized to A− treatment.Values are means and SE of each treatment.

Figure 2 .
Figure 2. Pelagic primary production over time standardized to A− treatment.Values are means and SE of each treatment.

Figure 4 .
Figure 4. CAP (constrained analysis of principle coordinates) biplot, based on the phytoplankton community matrix, and constrained by time, where each symbol represents each cosm with the individual treatment.The centroids are ambient (A) and warm (W) temperature treatments, also indicated by different colors.The arrows represent the different predictor variables and point into the direction of increasing value.The species shown represent the taxa that contribute 70% to the dissimilarity identified in the similarity percentages (SIMPER) analysis (TablesS1-S3).

Figure 4 .
Figure 4. CAP (constrained analysis of principle coordinates) biplot, based on the phytoplankton community matrix, and constrained by time, where each symbol represents each cosm with the individual treatment.The centroids are ambient (A) and warm (W) temperature treatments, also indicated by different colors.The arrows represent the different predictor variables and point into the direction of increasing value.The species shown represent the taxa that contribute 70% to the dissimilarity identified in the similarity percentages (SIMPER) analysis (TablesS1-S3).

Figure 6 .
Figure 6.CAP biplot based on the zooplankton community matrix, and constrained by time, where each symbol represents each cosm with the individual treatment.The arrows represent the different predictor variables and point in the direction of increasing value of each variable and their lengths indicate their relative importance.The centroids displayed are ambient (A) and warm (W)temperature treatments, which were based on a higher significance of the temperature effect compared to the HBCDD effect.The species shown represent the taxa that contribute 70% to the dissimilarity identified in the SIMPER analysis (TablesS4-S9).

Figure 6 . 21 Figure 7 .
Figure6.CAP biplot based on the zooplankton community matrix, and constrained by time, where each symbol represents each cosm with the individual treatment.The arrows represent the different predictor variables and point in the direction of increasing value of each variable and their lengths indicate their relative importance.The centroids displayed are ambient (A) and warm (W) temperature treatments, which were based on a higher significance of the temperature effect compared to the HBCDD effect.The species shown represent the taxa that contribute 70% to the dissimilarity identified in the SIMPER analysis (TablesS4-S9).J. Mar.Sci.Eng.2017, 5,18 13 of 21

Figure 8 .
Figure 8. Benthic primary production, standardized to the A− treatment to enable valid comparisons between time points (mean ± SE).

Figure 10 .
Figure 10.Conceptual diagram of direct and indirect effects of increased temperature and HBCDD addition in model ecosystems.

Figure 10 .
Figure 10.Conceptual diagram of direct and indirect effects of increased temperature HBCDD addition in model ecosystems.

Table 2 .
Summary of significant results from the multiple correlations of the benthic-pelagic parameters.r is the correlation coefficient, p is the level of significance.