Global environmental changes negatively impact temperate seagrass ecosystems

The oceans are increasingly affected by multiple aspects of global change, with substantial impacts on ecosystem functioning and food-web dynamics. While the effects of single factors have been extensively studied, it has become increasingly evident that there is a need to unravel the complexities related to a multiple stressor environment. In a mesocosm experimental study, we exposed a simplified, multi-trophic seagrass ecosystem (composed of seagrass, two shrimp species, and two intermediate predatory fish species) to three global change factors consisting of simulated storm events (Storms), heat shocks (Heat), and ocean acidification (OA), and the combination of all three factors (All). The most striking result indicated that when all factors were combined, there was a negative influence at all trophic levels, while the treatments with individual factors revealed species-specific response patterns. It appeared, however, that single factors may drive the multi-stressor response. All single factors (i.e., Storms, Heat, and OA) had either negative, neutral, or positive effects on fish and shrimp, whereas no effect was recorded for any single stressor on seagrass plants. The findings demonstrate that when several global change factors appear simultaneously, they can have deleterious impacts on seagrass ecosystems, and that the nature of factors and food-web composition may determine the sensitivity level of the system. In a global change scenario, this may have serious and applicable implications for the future of temperate seagrass ecosystems.


INTRODUCTION
Anthropogenic influence on global ecosystems has reached unprecedented rates with nearly the entire ocean (97.7%) being affected by multiple factors such as fishing, climate change, and ocean and land-based nutrients and toxins (Halpern et al. 2015). The evidence of the deleterious effects of global changes in the world's oceans is mounting, with increasing numbers of studies showing a range of associated impacts deriving from sea ice melt, sea surface temperature increases, sea level rise, and increased incidences of severe storms and ocean acidification (OA; e.g., Perry et al. 2005, Harley et al. 2006, Harvey et al. 2013, Mueter and Litzow 2013, Milazzo et al. 2013, Verges et al. 2014. With elevated temperatures, the frequency and magnitude of storms increase and can create more precipitation as well as cause physical damage to coastal areas (IPCC 2014). Increased severe weather events can lead to increased wave action in normally unexposed coastal shallow-water areas. Despite the relatively short duration of extreme events, it is such events that cause the largest amount of damage to ecological systems (Jentsch et al. 2007, Harvey et al. 2013. Additionally, OA has been identified as one of the most influential consequences of atmospheric CO 2 rise, causing significant alterations to the physiology of organisms creating repercussions for entire food webs (Harvey et al. 2013). Under current conditions, without substantial reductions in CO 2 emissions, it is predicted that there will be serious changes to marine organisms as well as ecosystem services by the end of the century . While it is clear that individual global change factors are influencing the marine environment, it is still less well understood, despite much current effort, what kind of effect the impacts have when in combination. Marine global change research has recently begun to explore interactive effects in order to understand and manage ecosystems influenced by multiple factors (Duarte et al. 2014, Murray et al. 2014, Boyd et al. 2018. Thomson et al. (2015) argue, in order to accurately predict how the seagrass ecosystem responds to global changes, it is necessary to examine the effects of both trend-and event-driven climate variables. Understanding how species and ecosystems respond in the short-term versus more long-term physiological changes is also of importance as some factors may be more acute (though increasingly frequent) in nature such as storm or heat shock events, while others progress more slowly and present a more long-term stress such as OA. For instance, organisms may respond to stress initially by increasing their rate of respiration (Vinagre et al. 2012, Gr€ ans et al. 2014, while, when exposure to factors becomes long term, oxidative stress can occur causing physiological changes resulting in alterations of lipids, proteins, and even DNA (Abele and Puntarulo 2004, Heise et al. 2006, Hernroth et al. 2012.
While much considerable effort has been made in the last decade to understand the influence of global changes, much of the earlier work focused, for practical reasons, on single species and single factors (O'Brien et al. 2019). Through this research, it has been clear that species responses to various factors differ widely, with some factors such as OA and temperature rise, having positive and negative effects depending on scenarios and species studied (for a review, see Harvey et al. 2013, Vargas et al. 2017. Though valuable for developing knowledge of a complicated and relatively new problem, it is difficult to draw ecosystem-level conclusions from them, which is where natural effects are, and will be seen in the future. Given the interconnectedness of ecosystems, it can be expected that beneficial or deleterious effects at one trophic level, particularly lower trophic levels, can have enormous consequences for entire food webs, such has been predicted with krill in the Arctic and Antarctic regions (Hill et al. 2013, Mcbride et al. 2014. Understanding the compounding problems created by multiple factors applied to species and ecosystems is a critical step in the process of mitigating negative impacts given that multiple factors can act in a complex non-linear interactive way (Coelho et al. 2013, Salo and Pedersen 2014, Halpern et al. 2015. Though healthy ecosystems tend to have a natural resilience to acute stress and species have a natural degree of adaptability to climate factors (Woodson et al. 2019), resilience is lost as problems become more chronic (Lake 2013). Such compromised ecosystems are likely to cross a threshold of resilience and shift to a non-recoverable state (Loreau et al. 2001, Lundberg andMoberg 2003) if research and management do not work to understand, reduce, and eliminate the stresses and threats to the oceans.
Although coastal areas, including seagrass meadows, are relatively resilient to environmental changes such as temperature and salinity alterations, it is unclear how the habitats will manage to adapt to the unparalleled rates of change the oceans are currently experiencing. Seagrass ecosystems, an ecologically critical habitat (Green and Short 2003), have already shown global gradual decline (e.g., Cambridge and McComb 1984, Baden et al. 2003, Waycott et al. 2009, Sobocinski et al. 2013. The system can also exhibit largescale, abrupt diebacks, triggered by extreme weather events such as hurricanes (Wilson et al. 2019) and periods of heat shock (Waycott et al. 2009, Thomson et al. 2015. For instance, extreme weather may alter the physical structure of habitat-forming species, such as seagrasses, and coastal habitats may also suffer from changes in species composition as mobile species shift ranges due to temperature alterations (Albouy et al. 2014, Hattab et al. 2014. Habitat-forming foundation species susceptible to climate change alterations create a situation where the entire ecosystem is vulnerable to shifts or reorganization as a result of habitat complexity loss . For instance, Ekl€ of et al. (2012) found that under increased temperature scenarios, invertebrate mesograzers were unable to control algal growth and therefore the seagrass system shifted to an algal-dominated habitat.
As habitat-forming systems, seagrass meadows are of high ecological value  because they support many organisms and function as important nursery grounds (Heck et al. 2003, Perry et al. 2018) for economically important fisheries species (Beck et al. 2001, Nordlund et al. 2017; it is therefore of critical importance to understand this system's response to predicated ocean changes. The current study evaluated a simplified seagrass system via a mesocosm laboratory experiment under acidified conditions, a stress parameter following a gradual trend in reduced seawater pH (increased pCO 2 ), as well as event-driven heat shock and storm scenarios. In particular, we aimed to understand (1) whether global climate change-associated events influenced a seagrass ecosystem and, if so, (2) how different trophic levels, from plants up to fish, responded to projected future ocean surface temperature rises, acidification, and storm events (causing turbidity), and lastly (3) whether global climate change-associated factors more negatively influence the seagrass ecosystem when occurring simultaneously. We hypothesized that each driver, when evaluated as a single entity, would have positive, negative, or no effect depending on the species, however, that the cumulative effect of all factors combined would be deleterious at all trophic levels.

Mesocosm organisms
Simplified seagrass mesocosms were created using twenty glass aquaria (61 9 39 9 35.5 cm; 84.5 L). The mesocosms consisted of seagrass (Zostera marina L.) collected directly outside the laboratory at the Kristineberg Marine Innovation and Research Center, Sweden (58°15 0 00.55″ N, 11°26 0 45.66″ E). The seagrass was collected by hand, digging down in the sediment in order to collect shoots with associated rhizomes and roots. After collection, the seagrass was taken directly to the laboratory and transplanted along with the sediment, into the aquaria. No algae were intentionally included in the seagrass collection, and any filamentous algae were removed prior to transplantation into the tanks. In order to minimize disturbance to the seagrass shoots to try to ensure the maximum health of the plants at the start of the experiment, the leaves were not scraped of any grazers though grazers were not intentionally added to the system. Some grazers may have been placed in the tanks either via the sediment or those smaller grazers remaining on the leaves during the transplant process; however, the fish are likely to have predated on the few unintended tagalongs during the acclimation phase or in the early days of the experiment. While a small amount of algae and grazers may have been included with the naturally collected seagrass when transplanted into the tanks, all tanks were collected, replanted, and set up in the same way and therefore no large differences would be expected between the tanks. The sediment in each tank was approximately 5 cm deep, which was enough to completely cover all rhizomes and roots. Following density measurements recorded from the same location (Staveley et al. 2017),~28 (AE2) shoots were planted in each aquarium. Mobile fauna was collected from seagrass habitats located in the surrounding area using beach seine. Two species of fish were selected for the experiment-the three-spined stickleback (Gasterosteus aculeatus L.) and the broadnosed pipefish (Syngnathus typhle L.)-representing intermediate predators mainly eating invertebrates. Sticklebacks are intermittent visitors in Z. marina being especially abundant in the autumn, whereas pipefish are resident in Z. marina . Two species of shrimp, the European rock shrimp (Palaemon elegans Rathke) and the Baltic prawn (Palaemon adspersus Rathke), were selected. Both are abundant in Z. marina habitat along the Swedish west coast feeding on filamentous algae and ❖ www.esajournals.org invertebrates (Jephson et al. 2008). Whereas P. elegans is migrating between different shallow coastal habitats, including physically demanding habitats such as rock pools, P. adspersus is either a strict resident or found in areas close to Z. marina meadows (Berglund 1980, Baden and Pihl 1984, Pihl et al. 2006. Once caught, the animals were placed in buckets and transported back to the laboratory where they were weighed and randomly placed in different aquaria, ensuring no differences in mean weights of species between treatment groups. Each aquarium received five specimens of G. aculeatus and three specimens each of S. typhle, P. elegans, and P. adspersus, respectively, with densities mirroring those found naturally in Z. marina meadows in the area (Emanuelsson 2008). Once all organisms had been placed in the aquaria, the mesocosms were given 1-3 d to acclimate prior to the start of the experiment. The aquaria were supplemented with live prey (Artemia) daily.

Experimental manipulations and treatment groups
The experiment was run for 28 d in August 2016, in a semi-open laboratory equipped with a glass roof in order to enable ambient light conditions. Each aquarium was set up with air diffusers and flow-through seawater pumped in from the bay outside (mean flow, 1.45 L per min AE 0.20 standard error [SE]), and therefore, water conditions fluctuated naturally. Twenty aquaria were randomly assigned to a control group or one of four treatment groups, with each group consisting of four replicates. The treatment groups received either one of three global change manipulated factors (T 1 -storm, T 2 -heat shock, and T 3 -OA or pH), or a combination of all three (T 4 -All, i.e., combined effects of storms, heat shocks, and OA). For the Storm condition tanks, turbid conditions were created by manually mixing the sediment with the water column every two hours for around 5 min over a 24-h period. This was done twice over the course of the experiment so that there were two storm events during the 28-d experiment. For the Heat treatment group, heat shock events occurred for a 48-h period, where ambient water temperature (~19°C) was raised by 10°C (~29°C) simulating heat waves experienced by coastal environments. The OA treatment had pure CO 2 bubbled into the tanks using a feedback pH-stat computer system (Aqua Medic, Bissendorf, Germany) to maintain a reduced pH of 7.7 (average pH expected by 2100 as compared to 8.1 today at the test site; Dorey et al. 2013) throughout the experiment. For the All treatment, the tanks were exposed to the same storm conditions for two 24-h periods, the same heat shock conditions twice for 48 h each time, and a reduced pH for the entire 28-d experimental period. Storms were conducted on experimental days 8 and 22, while heat shock events were performed on days 5-6 and 19-20.
Daily observations of tank conditions were recorded to ensure all equipment was functioning properly and experimental organisms were accounted for. The mesocosms were cleaned weekly to remove excess algal buildup on the glass and ensure that daily observations could be made. To monitor the conditions throughout the experiment, a number of parameters were measured regularly. The pH of the tanks with acidified conditions was measured using a WTW Multi 3430 pH meter 2-3 times per week. The pH meter was calibrated weekly according to manual instructions using NBS buffers (i.e., standard buffers for seawater; Perez and Fraga 1987). During the two storm events, turbidity was measured (HACH 2100 Q turbidity meter) prior to the start of the storm and then three times throughout the course of the event. For the heat shock events, the water temperature was also recorded prior to the start of both 48-h periods and then 7-10 times over the course of the 48 h.
Additionally, water chemistry conditions were recorded twice a week. Temperature, salinity, and pH were recorded in the aquaria, after which water samples were collected from all tanks in order to determine tank alkalinity. All samples were filtered using a 0.45-lm filter and remeasured for total pH (pH T ) using a pH meter calibrated with TRIS (Tris/HCl) and AMP (2-aminopyridine/HCl) buffer solutions with a salinity of 32, after which point samples were placed in the alkalinity machine (TA05 plus/TW alpha plus, SI Analytics, Mainz, Germany) for total alkalinity (TA) readings. Carbonate system parameters were then calculated from pH T and TA (Table 1).
There were no significant differences between treatments or between replicates within treatments for temperature, alkalinity, and salinity, ❖ www.esajournals.org respectively. Among all the parameters of the carbonate system, the high pH treatment had a significant effect (ANOVA; P < 0.001), but there was no significant difference (P > 0.05) between replicates within the pH treatment. The high pH groups had a pH TS of 8.03 AE 0.01 (pCO 2 of 429 AE 10 latm), and the low pH treatment had a pH TS of 7.65 AE 0.04 (pCO 2 of 1287 AE 564 latm). All treatments were supersaturated for both calcite and aragonite.
The increase in temperature of 10°C was selected given that shallow-water habitats, such as the seagrass ecosystem studied, naturally experience large variations in temperatures, particularly during warm seasons (Heise et al. 2006, Roth et al. 2010, and as heat shock events are projected to increase in number and duration (Russo et al. 2014, 2015, Fischer and Knutti 2015, Diffenbaugh et al. 2017. Additionally, it is expected that severe storm events will also increase (discussed in Philippart et al. 2011) with increasing levels of atmospheric CO 2 , making increased wave action causing turbid conditions a necessity to study as seagrass meadows rely on photosynthesis and are therefore highly affected by suspended particulate matter in the water column.

Animal response parameters
Response parameters were selected to cover both acute and chronic factors in animals of two different trophic levels. A tradeoff between number of animals per species and response parameters possible to analyze was carefully considered.

Physiological measurements
Respirometry: fish and shrimp.-Respiration, measured as oxygen consumption rate (MO 2 ), was recorded for the fish and shrimp. MO 2 was recorded for individuals placed in the respirometry chambers for all five experimental groups at the end of the 28-d experimental period. Organisms were placed individually in intermittent flow-through cylindrical respirometry chambers based on their weight. Glass respirometry chambers were specially designed according to recommendations from Clark et al. (2013) based on the studied species' mass and the volume of chambers (0.236 L for the small chambers and 0.491 L for the large chambers). The chambers were submersed in surface seawater (same source as the mesocosms). Water was continuously pumped through the chambers via a recirculating peristaltic pump. After calibration, the oxygen concentration was continuously measured using a FireSting O 2 system (PyroScience, Aachen, Germany). Freshly oxygenated water from the surrounding seawater flow-through bath was flushed through the respirometers for 10 min out of every 20-min measurement cycle, a process ensuring that oxygen levels within the chambers never dropped below 80% air saturation. Resting MO 2 was calculated for each organism after the measurement period (a minimum of 12-h nighttime measurement period for the fish and a minimum of 6-h daytime measurement period for the shrimp) using the decline in oxygen concentration in the respirometers over the 10-min period, when the water was recirculating within the chambers between flush cycles. Outliers (measurements of >2 standard deviations below the mean of the lowest 10% of measured values) were excluded from the calculations of resting MO 2 .
Respirometry chambers were cleaned daily between uses with hot water and then an ethanol ❖ www.esajournals.org rinse to control for possible background respiration caused by microbial growth. For the experimental treatment groups receiving acidified conditions throughout the study (i.e., pH and All), the seawater bath in which the respirometry chambers were placed during oxygen consumption recordings was also fitted with a pH computer system to ensure the water was held at pH 7.7. Fish analyses of oxidative stress markers.-1. Tissue homogenates.-Livers of G. aculeatus from each of the treatments were dissected on dry ice and stored in À80°C until analyzing the markers, including total antioxidant capacity (TCA), protein carbonyls (PC), and lipid peroxidase (LPO). The mean weight of the livers was 37.7 AE 13.7 mg, and there were no significant differences between treatments (F = 0.894; P = 0.437). The livers were homogenized in 750-lL ice-cold phosphate buffer saline (10 mmol/L Na 2 HPO 4 , 1.8 mmol/L KH 2 PO 4 , 2.7 mmol/L KCL, 137 mmol/L NaCl; pH: 7.4) with protease inhibitors (Protease Inhibitor Cocktail, #P8340; Sigma-Aldrich, St. Louis, Missouri, USA) and butylated hydroxytoluene to prevent further oxidation. After centrifugation at 16,000 9 g for 10 min at 4°C (5415R, Eppendorf AG, Hamburg, Germany), the protein concentration of supernatants was determined using BCA assay (Thermo Scientific, Waltham, Massachusetts, USA) with bovine serum albumin as standard. Aliquots were stored at À80°C prior to the use in the assay. For the TCA and PC analyses, supernatants were first cleaned by an additional centrifugation through a Micron centrifugal filter (# 42413; Millipore, Burlington, Massachusetts, USA). A positive control used for all the oxidative stress analyses was obtained by heat treatment (50°C for 45 min) of liver tissue from an individual of G. aculeatus, not included in any of the experimental treatments. All samples were analyzed in duplicates.
2. Total antioxidant capacity.-The TCA assay (#STA-360, OxiSelect; Cell Biolabs, San Diego, California, USA) is based on reduction of copper (II) to copper (I) by antioxidants. The analyses were prepared according to the manufacturer 0 s instruction and with components provided in the analyses kit. Briefly, absorbance units, measured with a microplate reader at 490 nm (Labsystems iEMS Reader MF, Helsinki, Finland), was calculated by subtracting the initial absorbance from the final, measured five minutes after addition of Copper Ion Reagent and followed by Stop solution. Uric acid was used to perform the standard curve, and values were converted to copper reducing equivalents, presented as CRE/mg protein.
3. Protein carbonyls. -For the PC assay (#STA-310, OxiSelect Protein Carbonyl ELISA Kit; Cell Biolabs), samples were diluted to 10 lg/mL and analyzed according to the manufacturer's instructions. The assay is based on derivatization of protein carbonyls to dinitrophenylhydrazine (DNPH) followed by immunoblotting with an anti-DNP antibody and a horse reddish peroxidase-conjugated secondary antibody. 3,3 0 ,5,5 0 -Tetramethylbenzidine was used as substrate, and absorbance was measured with a microplate reader at 450 nm (Victor TM 2030 Multiple Reader; Perkin Elmer, H€ agersten, Sweden). Concentrations of PC (nmol/mg protein) were calculated in relation to a standard curve, prepared from predetermined reduced and oxidized BSA standards, all provided in the kit.
4. Lipid peroxidase.-The LPO assay (# 21044, BIOXYTECH â MDA-586; Oxis Research, Los Altos, California, USA) is based on the reaction of the chromogenic reagent N-methyl-2phenylindole with malondialdehyde (MDA) at 45°C. The assay is regarded MDA specific as 4hydroxyalkenals are excluded by carrying out the reaction in hydrochloric acid with the addition of Probucol. The manufacturer's protocol was followed with one exception: The volumes were modified to half of the recommended. The absorbance was measured using microcuvette in a spectrophotometer at 586 nm (BioPhotometer, Eppendorf AG, Hamburg, Germany) after centrifugation of the samples. The MDA concentrations (nmol MDA/mg protein) were calculated in relation to an MDA-standard (provided in the kit) curve.
Shrimp phagocytosis assay.-Phagocytosis assays were performed for both species of shrimp, P. adspersus and P. elegans, in order to determine whether the individuals were immune-compromised as a result of exposure to the global change factors. Following the experimental exposure period, the shrimp were carefully removed from their tank and 10 lL of hemolymph was removed ❖ www.esajournals.org using a Hamilton syringe through a puncture made between the carapace and abdomen along the dorsal segment join. Directly after this 10 lL of heat-killed Bacillus thuringiensis (concentratioñ 10 8 cells/mL) dissolved in insect, Ringer saline solution, labeled with FITC (methodology in Kurtz 2002), was injected into the animal via the same puncture hole. Animals were then placed into an isolation aquarium for a period of 2 h. After 2 h, 5 lL of hemolymph was removed from each individual using a melanized capillary tube and placed in a Lab-Tek chamber slide well together with 250 lL of Grace's insect medium, where they were kept on ice for 15 min, followed by 30 min in a wet chamber. The phagocytosis assay is described in full in Roth and Kurtz (2009) and Roth et al. (2010). Trypan blue was added to each well for 15 min to quench the non-phagocytosed fluorescent bacteria, and then washed away with Grace's insect medium until the slides were only slightly blue. Subsequently, a glycerine mountant with DAPI was used to mount a cover. The ratio of hemocytes phagocytosing (identified by the fluorescent bacteria content) to those nonphagocytosing allowed phagocytic activity to be quantified as a measure of the individual's immunocompetence. A Leica UV DIC digital microscope with epifluorescent light was used to count total and phagocytosing bacteria in five fields of vision per well (one well per individual). A ratio was calculated by dividing the number of phagocytosing hemocytes by the total number of hemocytes to determine the percentage of hemocytes containing endocytosed bacteria (Steenbergen et al. 1978). The same methodology was used for both shrimp species, with each species kept separate throughout the assay. Treatments were processed in a random order and the slides coded and read blind (i.e., with no knowledge of the species or treatment the sample was related to).

Seagrass, Zostera marina
Seagrass growth and photosynthesis were recorded. At the start of the experiment (day 1), two shoots from all 20 mesocosm aquaria were punctured on the lower part of the leaf sheath using a needle following the method by Short and Duarte (2001) and then marked with a zip tie for relocation. Growth (mm) and number of new leaves were then recorded for each punctured shoot on days 7, 14, 21, and 28. Using the growth data collected weekly, the seagrass turnover rate was calculated for the average growth per day as well as for the average number of new leaves per day. In addition to investigating growth as in indicator of seagrass health, the photosynthetic efficiency of seagrasses was followed throughout the experimental period. A pulse amplitude modulated fluorometer (PAM; Aquation, Woy Woy, Australia) was used simultaneously with photosynthetic active radiation (PAR) loggers (Odyssey Dataflow Systems, Christchurch, New Zealand) to measure the photosynthetic yield of photosystem II (Φ PSII ) and record irradiance (lmol photons Ám À2 Ás À1 ), respectively. Absorption factor (AF) values were also recorded on seagrass collected from the same location as the seagrass used in the mesocosms using a Diving PAM (Heinz Walz, Effeltrich, Germany) following the method described by Beer et al. (2014). Electron transport rates (ETR) were then calculated using Eq. (1).
The derived values of ETR were fitted against PAR values by a double exponential function (P/I curve) suggested by Platt et al. (1980), where the slope alpha (a) of the curve is the photosynthetic efficiency. Photosynthesis measurements were carried out at the start, middle, and end of the 28-d experimental period. Measurements were taken for each of the five treatment groups for a 24-h period. The sensors were then moved to another tank and run for another 24-h period (n = 3 for each treatment at the start, middle, and end of the experimental period). Additionally, during the climate manipulated events, three replicates were taken for the treatment tanks and three replicates for the control tanks during the storm and heat shock treatments (n = 3 9 2 storm events for control as well as Storm and All and n = 3 9 2 heat shock events for control as well as Heat and All).

Statistical analysis
Fish and shrimp respirometry, shrimp phagocytosis assay, and seagrass growth data were analyzed using STATISTICA version 6.4 (StatSoft, Uppsala, Sweden). Data were checked for normal distribution prior to analysis using the Shapiro-Wilk test and natural log-or squareroot-transformed when appropriate in order to ❖ www.esajournals.org validate statistical assumptions. Subsequently, the data were analyzed using general linear models based on nested analysis of variance (ANOVA), where tank was nested in treatment for seagrass growth, and a one-way ANOVA used for the respirometry and phagocytosis data. When statistically significant differences were seen between treatments, Fisher's least significant difference (LSD) post-hoc tests were run in order to determine which treatments differed from one another.
It should be noted that a large number of pipefish, S. typhle, died throughout the experiment and were therefore not used in the respirometry assessments or analyses of stress markers. However, the number of days alive for the pipefish was analyzed using a nested ANOVA (tank nested in treatment) and no differences were found between treatment groups (F 4,15 = 0.976, P = 0.449).
The number of G. aculeatus used for the analyses of TCA, PC, and LPO differed between treatments as well as between aquaria. Homogeneity of data, tested with the Shapiro-Wilk method, failed, and thus, the Kruskal-Wallis one-way analysis of variance on ranks was used for investigating potential aquarium effects. As no significant differences were found, all data from each of the treatments were regarded as true replicates: Storm, n = 12; Heat, n = 17; pH, n = 13; All, n = 14; and Control, n = 16. One-way ANOVA on ranks was used to explore differences between treatments, and a multiple comparisons test was used for multiple pairwise comparisons. These analyses were also performed using the STATISTICA software. A oneway ANOVA was used to analyze alpha slope data for the seagrass at the end of the experiment also using STATISTICA. Additionally, a two-way ANOVA was used to test for differences between all treatments and the different periods (start, middle, and end of the study) and a one-way ANOVA was used to test for differences within the Heat and Storm treatments throughout the study period (start: first heat wave or first storm; middle: second heat wave or second storm). Data were tested for equal variance and normality using a Shapiro-Wilk test. Seagrass global driver event analyses were performed with SigmaPlot v. 12.3 (Systat Software, San Jose, California, USA).

Fish-respirometry and oxidative stress markers
The oxygen consumption rate (MO 2 ) was significantly lower in three of the disturbance treatments (i.e., Storm, Heat, and All) compared to the control, while the pH treatment did not differ from the control (Fig. 1, Fisher's LSD post-hoc test). The control and pH treatments have añ 20% higher MO 2 than the other treatments. For the TCA results, the mean (AESE) concentration of CRE obtained from the different treatments was 530 AE 160 mg protein À1 . There was no significant difference between the median values of treatments (P = 0.221). For the PC results, the mean concentrations of PC (nmol/mg protein) AE SE obtained from the different treatments indicated that pH and the combined treatment were significantly different from all the other treatments (P < 0.05; Fig. 1). For the LPO analysis, the mean concentration of MDA (nmol/L per mg protein) AE SE obtained from the different Fig. 1. Bars indicate Gasterosteus aculeatus mean (AEstandard error [SE]) oxygen consumption rate (MO 2 ). Black diamonds show mean values (AESE) of protein carbonyls (PC) per mg protein oxidative stress marker measured in liver tissues of G. aculeatus. Letters indicate significant differences between treatments (having the same letter indicates no significance) based on Fisher's least significant difference post-hoc results with uppercase letters showing post-hoc results for the PC data and lowercase letters for the MO 2 data. treatments was 3.06 AE 0.62 without any significant differences between treatments.

Shrimp respirometry and phagocytosis
MO 2 rates differed between groups for P. adspersus (Fig. 2a). The control group differed significantly only from Storm, which had a reduced respiration rate. No other significant differences were seen between groups. Phagocytosis results for P. adspersus showed the lowest phagocytotic activity in the pH treatment, which was significantly lower than all other treatments, especially compared to the control. Additionally, the Storm and All treatments were also significantly immunocompromised when compared to the control (Fig. 2a).
Compared to P. adspersus, a different MO 2 pattern was seen for P. elegans (Fig. 2b). In this case, the lowest MO 2 was seen in the control group, which was clearly different from all other treatments. In turn, the Storm treatment was significantly lower than the other three treatments, that is, Heat, pH, and All (Fig. 2b). In terms of phagocytotic activity by P. elegans, all treatments showed significantly lower values compared to the control group (Fig. 2b); however, no differences between only the treatment groups were seen (Fig. 2b).

Seagrass growth
At the end of the experiment period (day 28), the longest seagrass growth (length in mm) from initial needle puncture was seen in the control and storm groups, while the lowest was in the All treatment. The control group was significantly different from the Heat, pH, and All groups.
The total number of new leaves grown (i.e., showing no hole from the needle puncture done on day 1 of the experiment) differed slightly The black diamonds on the left (a) show the phagocytotic activity as a proportion of phagocytosis (AESE) for P. adspersus, and the black diamonds on the right (b) show the phagocytotic activity as a proportion of phagocytosis (AESE) for P. elegans. Letters indicate significant differences between treatments (having the same letter indicates no significance) based on Fisher's least significant difference post-hoc results with uppercase letters showing post-hoc results for the phagocytosis data and lowercase for the MO 2 data. between the different treatment groups (Fig. 3). No significant differences between groups were, however, seen on day 28 (nested ANOVA F 4,13 = 2.256, P = 0.119).

Seagrass photosynthesis
For seagrass photosynthesis evaluating alpha slope, measured using the PAM, no differences were seen between treatments (one-way ANOVA; F 4,10 = 1.406, P = 0.301) at the end of the experiment (Fig. 4). Additionally, no statistical differences were seen when analyzing only measurements taken at the time of disturbance events, that is, when the storm or heat shock events were taking place (F 6,58 = 1.342, P = 0.253).

DISCUSSION
Our results suggest that temperate seagrass (Z. marina) ecosystems are sensitive to expected future global changes. Each species had different responses to simulated storms, simulated heat shock, and OA but when the three stress factors were combined all species were significantly affected (Table 2). Interestingly, the results show no clear interaction between the tested factors. This may be a result of a phenomenon seen by Brennan and Collins (2015) in which a species response to multiple factors tends to be swayed by a single dominant driver (factor). From a broader ecological perspective, if stress responses tend to be driven by single dominant factors, managing ecosystems for resilience may become slightly less challenging despite the increased risk of multi-stressor environments; however, the problem still remains that an ecosystems dominate stress factor is not universal but rather often unique for each species and would strongly depend on intensity and duration.

Fish response
The respiration results for the three-spined stickleback G. aculeatus showed significantly higher oxygen consumption rates in the control group compared to the Storm, Heat, and All treatment groups. Perhaps this lowered respiration in three of the treatment groups is a result of metabolic depression or something like the faint or freeze phenomenon seen in other animals under acute stress (Lang et al. 2000). During the  ) of the PI curve (Electron transport rates/photosynthetic active radiation). No significant differences were seen between treatments. faint or freeze phase of the stress response, an animal's heart rate typically slows and breathing decreases. While this behavior is typical of a response to a physical threat such as a predator, bradycardia (slowed heart rate) has been documented in fish exposed to stress as a result of hypoxic conditions (Axelsson et al. 2002, Monteiro et al. 2013. The pH treatment, however, had a similar oxygen consumption rate to that of the control. The treatments experiencing acidified conditions (pH and All) were exposed to lowered pH levels for the full duration of the 28d experiment, whereas storms and heat shock treatments experienced these stresses as events that the tanks were exposed to twice each throughout the experimental period. This difference may be an explanation for why there was no change in the pH treatment's respiration as it may indicate a level of habituation to this chronic exposure. Other experiments have found that longer-term stress exposure in fish can result in habituation (Lankford et al. 2005) meaning the fish develop an adaptation to no longer respond to the stress and in some cases have a reduced ability to respond (Barton et al. 1987). However, while the fish no longer show a respiratory response to the stress, internal physiological consequences of chronic stress occur (Barton 2002). This is exactly what we see in our results from the oxidative stress analysis. While respiration rates were unchanged in the pH group, the acidified condition still clearly increased the level of protein carbonyls when compared to the control group as well as the Storm and Heat treatments.
Here, it appears that OA may be the dominant stress factor. It is likely that underlying mechanisms of stress resilience can facilitate habituation of fish, and explain unaffected growth and survival of fish larvae, under the pressure of OA (Sswat et al. 2018).

Shrimp response
Despite both shrimp species being similarly sized, and exposed to the same conditions, basic respiration rate of the control group of Palaemon adspersus was, surprisingly, double that of Palaemon elegans controls. Notably, P. adspersus is typically considered a specialist inhabiting Zostera marina areas, while P. elegans is more of a generalist inhabiting many different shallow-water habitats (Berglund 1980). From the studies of Notes: The * symbol indicates a significant difference, while "ns" indicates no significance "not significant." Arrows indicate whether the difference from control is higher (↑) or lower (↓). Berglund (1980) and Berglund and Bengtsson (1981), P. elegans also inhabits physically demanding rock pools and is far more tolerant when exposed to large variations in salinity, temperature, and hypoxia than P. adspersus. If any explanation exists for this interspecific difference, it could be that P. elegans, via a different behavioral strategy, saves energy with the ability to move away from a factor and thereby increasing the respiration in the face of extreme habitats; a response shown here under the global factors presented in the experiment. Pihl et al. (2006) compared density and biomass of the two shrimp species in Z. marina meadows and areas previously covered with Z. marina. Both species were found in the two habitat types and especially P. adspersus showed to be more active during night than day compared to P. elegans, thus supporting the theory about differences in behavior of these common shrimp species. For P. adspersus, with the generally high respiration rate, the Storm treatment seemed to be the dominate stress factor and the only factor that resulted in a significantly changed and decreased respiration rate compared to all other treatments including the control.
The two shrimp species showed a marked difference in their immunocompetency response to the heat shock treatment, with P. elegans showing a much lower immune functioning after exposure to heat shock, as compared to all other treatments, whereas P. adspersus showed no immune deficiency due to heat exposure. A heat shock study conducted on another seagrass-associated species, Idotea baltica, showed a decrease in immunocompetency of~50%, a similar result to that shown for P. elegans (Roth et al. 2010). A similar response was also seen for the American lobster Homarus americanus, when exposed to increases in temperatures (Steenbergen et al. 1978). While P. elegans showed significantly decreased immunocompetency for all treatments and P. adspersus did not, both species were significantly affected by all factors combined, a result that has very negative potential consequences for these species in many Z. marina meadows in the future. Phagocytosis constitutes an important defense against pathogens in crustaceans. Therefore, decreased immunocompetency of both shrimp species may highly affect their defense against pathogens, which would have serious repercussions for survival of shrimp populations, especially considering the increasing abundance and distribution of pathogens in the wake of climate change (Hernroth and Baden 2018).
Overall, the results indicate that P. adspersus is more tolerant to global change stress factors than P. elegans. This was somewhat unexpected, since P. elegans lives both in seagrass meadows and more physically demanding environments with rapid changes in salinity, temperature, and oxygen, and has shown better survival in extreme situations than P. adspersus (Berglund 1980, Berglund andBengtsson 1981). An explanation could be that P. elegans cope with energy-demanding stresses by moving to less stressful environments, a so-called flight response, whereas P. adspersus are more stationary and exhibit a fight response. Exposing P. elegans to in vitro stress situations, where it is unable to move away, may be compromising their normal flight reaction causing a turn to energy-demanding respiration and resultant immune suppression.

Seagrass response
Strikingly, seagrass shoot growth (measured as length increase of existing leaves) in the All treatment was lower as compared to the control and Storm groups. The seagrass response to the different treatments indicated that growth was highest in the Storm treatment and control group. The slightly higher growth in the Storm treatment may be a result of the storms causing turbid conditions that shaded the seagrass and triggered a lengthening response in the plant. This type of elongation in the response to shading has been documented in seagrass (Ralph et al. 2007, Deyanova et al. 2017 and is the plants' way of growing up to the light. Another possibility is that increased velocity could also reduce the diffusive boundary layer around the leaves, which would result in easier CO 2 uptake, potentially stimulating growth of the plant (Lucas 1983). Conversely, while no significant differences were seen in number of new leaves, there was an indication that the fewest number of new leaves was seen in the Storm treatment, possibly attributable to increased energy reserves being devoted to lengthening instead of the plant devoting resources to new leaves. Interestingly, Sunday et al. (2017) demonstrate that seagrass shoot density is predicted to increase with OA, and in some cases, this is found to be true; however, in other instances, the shoot density showed either a decrease or no change at all. For our results, while not significant, there was an indication of increased turnover (number of new leaves) with decreased pH. This is in keeping with our predictions given that seagrass photosynthesis is limited indirectly by the availability of CO 2 (Borum et al. 2016), and with increased CO 2 , the plant has increased resources (Koch et al. 2013) to devote to becoming denser (Zimmerman et al. 2017).
A possible explanation as to why there was a negative effect seen in seagrass growth for the combined treatment while not for individual treatments is the oxygenation of the sediment. Seagrass respiration increases faster than photosynthesis in high temperatures (Larkum et al. 2006), which also leads to higher respiration in the belowground parts within sediment, thereby reducing the sediment oxygen content which promotes anaerobic respiration. Increased temperatures will also have a direct effect on the microbial community within the sediment, further decreasing the oxygen content of the sediment and eventually leading to toxic hydrogen sulfide production that will negatively affect seagrass productivity (Greve et al. 2003, George et al. 2018). The shading caused by the storm event in the experiment would have direct negative effects on the amount of oxygen produced in the leaves of the seagrass, thus reducing the amount of oxygen transported from the leaves to the roots. Typically, some of this oxygen leaks into the sediment, ensuring the prevention of sulfide formation around the root tip, which would otherwise cause sulfide damage to the plant (Greve et al. 2003). While the Storm treatment alone had no significant effect on the seagrass, possibly due to the high mixing, and therefore potentially no effect on the sediment properties, it is possible that when combined with the effects caused by the heat, changes to the sediment were so severe as to cause decreased plant growth.
Seagrass plants are, to some extent, resilient and adapted to variable and dynamic environments. As such, a temporary compensatory response could be the reason for maintaining the same photosynthetic efficiency for all treatments compared to control. Therefore, longer exposure to such unfavorable conditions could result in a decrease in photosynthetic efficiency (increase in alpha slope), although such assumptions must be further investigated.

Possible ecological implications
Global environmental changes create a multitude of alterations from temperature increases, OA, increased storm events, and more, and each one of these changes influences species and ecosystem differently (Harvey et al. 2013, Halpern et al. 2015. Climate change adds to stress from already existing anthropogenic factors such as fishing and ocean and land-based nutrients and toxins (Harvey et al. 2013, Halpern et al. 2015. A relevant example is higher frequency and intensity of rain due to climate change, which results in even more washout of nutrients to the coastal water, increasing eutrophication. Given this complexity, and the seemingly immeasurable number of interactions possible between abiotic and biotic factors, studying the effects of global changes is infinitely complex and yet it is absolutely critical in order to try to mitigate future effects (Boyd et al. 2018). This study is an attempt to evaluate a simple multi-trophic seagrass ecosystem in response to interacting global change factors. While data interpretation is complex, it is necessary to understand how the system responds not to a single factor but multiple factors as interactions take place. This is evidenced by the seagrass growth results from the current study where length measurements were most negatively affected by interacting climate factors. Using model simulations, a study by Zimmerman et al. (2015) indicated that photosynthesis in Z. marina meadows on the east coast of the United States was stimulated by OA, enough to offset the negative impacts of increased temperature. Bach et al. (2018) showed that pelagic plankton communities vary greatly in their response to global change scenarios in the laboratory as compared to under natural conditions and they therefore highlight the importance of performing research in conditions as close to the natural environmental complexity as possible.
The results from this experiment complemented with results from previous ecological studies give valuable insight into future scenarios and explanations to processes going on in Z. marina meadows along the Swedish west coast, an ecosystem already showing biodiversity changes as a result of anthropogenic stresses (reviewed by € Ostman et al. 2016). Zostera marina has changed in spatial extent and species composition the last three decades (Baden and Pihl 1984, Baden et al. 2003, Pihl et al. 2006. Filamentous algae, thriving from a surplus of nutrients, have gradually smothered the Z. marina, making light and nutrients less available for the plants. The areal extension of Z. marina has decreased by about 60% since 1980s (Baden et al. 2003, Nyqvist et al. 2009). Simultaneously, overfishing has left~10% of top predators in coastal areas, increasing the abundance of intermediate predatory fish (dominated by G. aculeatus) by as much as ten times (Sved€ ang and Bardon 2003. However, despite their opportunistic behavior and apparent adaptability, the current results show global change consequences and, therefore, indicate problems for this species. Despite our findings that P. adspersus may be more capable to cope with future climate change factors simulated in the experiment, Z. marina meadows have become a favorable habitat for P. elegans as well. Therefore, future changes will have consequences for both of these important crustacean species inhabiting seagrass meadows and surrounding habitats.
Pipefish are sensitive to changes in the environment (Pollom 2014) and strongly favor Z. marina without filamentous algae (Sundin et al. 2011). Additionally, research has shown that pipefish are influenced by turbidity and shifts in oceanic pH, with it altering their mating behavior (Sundin et al. 2010(Sundin et al. , 2013. The present and the expected exacerbated conditions in Z. marina with even more filamentous algae, turbid and turbulent warmer water with large changes in oxygen concentrations and OA, may be disadvantageous to the survival of stationary pipefish, and may explain the generally high mortality in this experiment. Most notably, the results indicate that, while individual factors may be either negative, positive, or cause no effect depending on the species, all trophic levels in our simulated seagrass ecosystem were deleteriously affected by the combination of all three factors. Hence, even though single factors may drive the response in a multifactor environment, all species show a heightened negative response compared to control when the global factors are combined. This has very applicable and potentially challenging implications for the future of the seagrass ecosystems on the Swedish west coast. As Pecl et al. (2017) pointed out "Although species may adapt to changing climates, either through phenotypic plasticity or natural selection (Valladares et al. 2014), all species have limits to their capacity for adaptive response to changing environments (Williams et al. 2008). . ." and in a future where oceans will be experiencing a multitude of changes, species and ecosystems are unlikely to be able to adapt quickly enough.

ACKNOWLEDGMENTS
The authors extend a big "thank you" to the interns Jacqueline Fuchs, Stefania Charisiadou, Felisa Friedman, and Johanna Hormann who assisted with the setup and running of the experiment. Our appreciation also to the staff at the Sven Lov en Centre for Marine Sciences-Kristineberg, Sweden, and in particular Bengt Lundve for all the logistical assistance and Alex Ventura for his guidance with the water chemistry equipment. We are very grateful to Jeroen Brijs for a tremendous amount of help with respirometry setup and equipment and data processing, as well as Michael Axelsson and Erik Sandblom for loan of the equipment. The assistance offered from Fredrik Jutfelt made the experiment possible and we are extremely appreciative of his help. Additionally, we are thankful for Lina Rasmusson's help with ideas for stress measurements in seagrass and thank Jakob Walve for the loan of the turbidity meter. Also, much thanks to the many, many others who assisted with either brains or brawn during planning, setup, or execution of the experiment. The stay at the Sven Lov en Centre for Marine Sciences was funded by a scholarship from Kungliga Vetenskapsakademien (