Does Warming Enhance the Effects of Eutrophication in the Seagrass Posidonia oceanica?

Seagrass meadows are disappearing at rates comparable to those reported for mangroves, coral reefs, and tropical rainforests. One of the main causes of their decline is the so-called cultural eutrophication, i.e., the input of abnormal amounts of nutrients derived from human activities. Besides the impact of eutrophication at a local scale, the occurrence of additional stress factors such as global sea warming may create synergisms in detriment of seagrass meadows’ health. In the present study, we aimed to evaluate if plants undergoing chronic cultural eutrophication and plants growing in relatively pristine waters are more (or less) sensitive to heat stress, nutrient load and the combination of both stressors. To address this question, a mesocosm experiment was conducted using Posidonia oceanica collected from two environments with different nutrients load history. Plants were exposed in controlled conditions to high nutrient concentrations, increased temperature and their combination for 5 weeks, to assess the effect of the single stressors and their interaction. Our results revealed that plants experiencing chronic cultural eutrophication (EU) are more sensitive to further exposure to multiple stressors than plants growing in oligotrophic habitats (OL). OL and EU plants showed different morphological traits and physiological performances, which corroborates the role of local pressures in activating different strategies in response to global environmental changes. EU-plants appeared to be weaker during the treatments, showing the greatest percentage of mortality, particularly under increased temperature. Temperature and nutrient treatments showed opposite effects when tested individually and an offset response when combined. The activation of physiological strategies with high energetic expenses to cope with excess of nutrients and other stressors, could affect plants present and future persistence, particularly under eutrophic conditions. Our results represent a step forward in understanding the complex interactions that occur in natural environments. Moreover, unraveling intraspecific strategies and the role of local acclimation/adaptation in response to multiple stressors could be crucial for seagrass conservation strategies under a climate change scenario.


INTRODUCTION
Seagrasses are important aquatic angiosperms that form meadows of great ecological and economic value for marine and global ecosystems. The productivity of seagrass meadows and their ability to capture carbon are comparable to those of terrestrial forests (Fourqurean et al., 2012). These macrophytes are considered ecosystem engineers due to their habitatforming capacity, providing food and shelter for a range of organisms such as finfish, shellfish, waterfowl, and herbivorous mammals (Boudouresque et al., 2009). There is strong evidence that seagrasses are disappearing at rates comparable to those reported for mangroves, coral reefs, and tropical rainforests (Waycott et al., 2009). Such decline has been reported worldwide especially in populations occurring in sheltered embayments and lagoons, where water recirculation is low and prone to nutrient loading from neighboring human population, the socalled "cultural eutrophication" (Touchette and Burkholder, 2000). Besides the impact of local anthropic pressures such as cultural eutrophication, additional global stressors (e.g., sea warming) may generate synergisms, thus increasing loss rate of these precious ecosystems (Lloret et al., 2008). However, little is known about the potential interplay between multiple sources of stress and the response that meadows growing in different environmental conditions may have.
Eutrophication or over-enrichment of nitrogen and phosphorus in the water column has several indirect and direct effects on the physiology of seagrasses (Unsworth et al., 2015;Ceccherelli et al., 2018). Excessive inorganic nitrogen (N i , as NO 3 − and NH 4 + ) concentrations can indirectly inhibit seagrass growth and survival by reducing light availability due to stimulation of phytoplankton, macroalgae and epiphytic algae overgrowth (Touchette and Burkholder, 2000). Additionally, N i enrichment can directly affect growth in several seagrass species by altering their cellular function and generating a negative physiological response (Burkholder et al., 2007). Some studies suggest that feedback inhibitory mechanisms for N i uptake are missing in seagrasses as a result of an evolutionary adaptation to oligotrophic habitats (Burkholder et al., 1992;Touchette and Burkholder, 2000). Thus, as observed in the eelgrass Zostera marina (Burkholder et al., 1992), plants continuously uptake N i even during dark periods, when surrounding water is rich in nitrate. This uncontrolled uptake of N i generates metabolic imbalances because N assimilation is a highly energy-demanding process (Touchette et al., 2003). The toxic effect of ammonium addition was observed to be stronger during winter in Zostera noltii (Brun et al., 2002), suggesting interactions between different environmental drivers. In the latter, the addition of NH 4 + implied a mobilization of non-structural carbohydrates accompanied by a reduction in growth. Nitrate can be stored in vacuoles or reduced in the cytosol to nitrite and subsequently in the chloroplast to ammonium, which is then converted to glutamine (Touchette and Burkholder, 2000). The tissue content of carbon and nitrogen has provided good indication of the nutritional status of seagrasses, integrating the long-term nutrient exposure history of plants (Lee et al., 2004). Variations in N content can arise from specific environmental conditions. For instance, Mediterranean Posidonia oceanica growing in high CO 2 natural conditions showed a clear up-regulation of nitrate transporter genes, when fertilized with nutrients, suggesting a plastic response of plants for balancing C:N ratio (Ravaglioli et al., 2017). Although the C content shows low variability across species (Duarte, 1990), the rate of change in C:N ratio should shift from high to small as nutrient supply meets the plants' demands. However, inferences on a particular species should be taken carefully as substantial inter-and intra-specific variation is common (Touchette and Burkholder, 2000).
In addition to eutrophication, global increase of sea temperature also threatens seagrasses. Rapid sea warming is occurring at unprecedented alarming rates due to anthropogenic activities altering climatic conditions worldwide (Francour et al., 1994;Bianchi, 2007;Vargas-Yáñez et al., 2008;Coma et al., 2009;Lejeusne et al., 2010;Marbà et al., 2015). The Intergovernmental Panel on Climate Change (IPCC) predicts a global increase of 2.58 • C in mean sea surface temperatures (SST) by the end of the century (IPCC, 2019). In the short term, SST anomalies in the form of heat waves have also been observed, e.g., in 2010/11 in the southeast Indian Ocean, where a heatwave hit the Western Australian coast, experiencing record water temperatures of 2-4 • C above long-term averages for about 10 weeks during late summer (Pearce and Feng, 2013;Kendrick et al., 2019). In the Mediterranean Sea, warming due to heatwaves, particularly from late spring/early summer to early autumn, are expected to be stronger, lasting up to a whole month (Diaz-Almela et al., 2007;Hobday et al., 2016;Oliver et al., 2018). For instance, in 2003 the average temperature reached at 1 m depth was 25.6 ± 1.5 • C, with a maximum of 28.6 • C for the Gulf of Naples (Garrabou et al., 2009). The increasing frequency of those events raised the alarm on their impacts to marine biota, including seagrass meadows. For instance, mortality of P. oceanica has been observed in the Balearic Islands in 2003 and 2006 after summer heatwaves (Díaz-Almela et al., 2009;Marbà and Duarte, 2010). Moderate temperature increase can stimulate biochemical reactions such as photosynthesis and respiration (Olsen et al., 2012;Marín-Guirao et al., 2018), or even promote flowering Marín-Guirao et al., 2019). However, beyond certain thresholds, thermal stress affects the stability of photosystem II (PSII) reaction centers in the chloroplasts and reduces electron transport rates (York et al., 2013;Marín-Guirao et al., 2016;Repolho et al., 2017;Ruocco et al., 2019a). The acceleration of respiration over photosynthetic rates can also lead to carbon imbalances under heat stress (Collier and Waycott, 2014;Marín-Guirao et al., 2018). In any case, the effect of temperature is highly variable and the specificity of a particular species or population responses will be highly dependent on the natural regimes where plants grow (Collier et al., 2011;Winters et al., 2011;Marín-Guirao et al., 2018).
Interactions between multifactorial environmental conditions are indeed influencing the physiology of marine organisms. However, experiments assessing the interactions between temperature and nutrients are rare due to experimental and budgetary constraints. Some studies on seagrasses demonstrated that synergistic interactions occur when plants are exposed to a combination of stressors, such as osmotic, light, temperature and eutrophication (Touchette and Burkholder, 2000;Collier et al., 2011;Villazán et al., 2013Villazán et al., , 2015Ontoria et al., 2019b). However, the responses were highly specific for the species assessed and varied according to the life stage of the plant. For instance, in Enhalus acoroides seedlings, the effect of interaction between temperature and nutrient enrichment was weak (Artika et al., 2020). On the other hand, an antagonistic interaction was observed in Cymodocea nodosa adult plants, where acidification improved ammonium assimilation and buffered the enhanced respiration promoted by temperature (Egea et al., 2018).
The endemic seagrass P. oceanica forms extensive monospecific meadows on rocky and sandy bottoms, representing one of the most significant benthic species in the Mediterranean Sea (Serra and Mazzuca, 2011). Declines in P. oceanica populations over the last decades have been attributed, among many other causes, to both sea warming and eutrophication (Pergent et al., 1999;Seddon et al., 2000;Ruiz et al., 2001;Pergent-Martini et al., 2006;Díaz-Almela et al., 2009;Marbà and Duarte, 2010;Marbà et al., 2014;Moore et al., 2014;de los Santos et al., 2019). However, the consequences of the combined effect of both stressors on the plant physiology are virtually unknown. A recent microcosm study on P. oceanica examined the interactions between ammonium addition and high temperature (Ontoria et al., 2019a). In this case, synergism between factors caused a decrease of about 70% in photosynthetic performance indicating the harming potential of eventual heat waves in areas exposed to high anthropic pressures. Hence, the presence of meadows in highly impacted and nutrient enriched areas (e.g., the Gulf of Naples, Italy; Cianelli et al., 2012) is cause of concern and raises questions about the physiological strategies eventually adopted to respond to additional stress. For instance, Halophila ovalis plants growing in turbid waters were more sensitive to further shading in respect plants growing in clear waters that were able to decrease shoot density and increase photochemical efficiency (Fv/Fm; Yaakub et al., 2013).
In this study, we aimed to evaluate whether plants undergoing chronic cultural eutrophication and plants from relatively pristine sites are more (or less) sensitive to heat stress, nutrient load and the combination of both stressors. To address this question, a mesocosm experiment was conducted using P. oceanica plants from two environments with different nutrients load history. Plants were exposed in controlled conditions to high nutrient concentrations, increased temperature and their combination for 5 weeks, to assess the effect of single stressors and their interaction. At the end of the exposure period, we analyzed morphological, photophysiological and biochemical endpoints including the content of C and N to corroborate nutrient uptake.

Sampling Sites
Fragments of P. oceanica consisting of a horizontal rhizome bearing 10-20 vertical shoots, were collected by SCUBA diving on May 15 -16th 2019, from shallow-water meadows (7-9 m depth) from two locations with different history of nutrient loads, located at 8 nautical miles apart ( Figure 1A). Plants were collected from a meadow in the vicinity of Spiaggia del Poggio (Bacoli) in the Gulf of Pozzuoli (Italy, 40 • 47.930 N; 14 • 05.141 E) and from a meadow in the surrounding area of Castello Aragonese in the Island of Ischia (Italy,40 • 44.114 N;13 • 57.866 E). Because of their proximity, both locations experienced similar Sea Surface Temperature (SST) regimes. The island of Ischia is a marine protected area since 2007, whereas the Gulf of Pozzuoli is heavily impacted by human activities (Tornero and Ribera d'Alcalà, 2014;Chiarore et al., 2019). The Bacoli sampling site was close to the city of Bacoli and to the one of Baia, an important commercial harbor since Roman age. Sediments in the area of the Bacoli sampling site are classified as sandy with mud being an abundant component (Celia Magno et al., 2012), whereas P. oceanica meadows around Ischia are settled on a "matte, " a dense mixture of rhizomes, roots and accumulated sediment (personal observation). The marine area around Baia/Bacoli is particularly impacted through dense urban settlements, intense maritime traffic and mussel aquaculture (Appolloni et al., 2018). Additionally, along the region of the Gulf of Pozzuoli, untreated or not appropriately depurated sewage discharges are still present (Margiotta et al., 2020). The different exposure to anthropogenic pressures determines the current conditions of P. oceanica meadows at both locations. According to the percentage of N leaf content, a value that integrates the eutrophic status of coastal areas over time (Yang et al., 2018), the site where the Bacoli plants grow appears more eutrophic. The N leaf content value, in fact, was almost twice in Bacoli (%N leaves = 1.89 % ± 0.2; C/N ratio = 16.7 ± 0.9) than in Ischia (%N leaves = 0.97% ± 0.2; C/N ratio = 33.2 ± 2.4) (Supplementary Table 1 (Buia et al., 2004). Meadow percentage cover is almost 30% higher in Ischia than in Bacoli. Additionally, in the Bacoli site the shoot density (122.6 ± 70.3 number of shoots per m −2 ) is half than the minimum value defining a "highly disturbed meadow" (see Supplementary Table 1). Rhizomes of Ischia plants were strong and thick, while rhizomes from Bacoli plants were brittle and broke off easily (personal observation). Epiphyte cover on leaves was on average almost 4 times higher in Bacoli compared to Ischia plants (Helber et al., unpublished data). Based on all the above observation, we consider plants collected in Bacoli as more exposed to eutrophic conditions, compared to plants collected in Ischia. Hereafter, we refer to plants collected in Bacoli as relatively eutrophic (EU plants), and plants collected in Ischia as relatively oligotrophic (OL plants). From a genetic point of view, comparative analysis of P. oceanica populations at the two sites using 9 microsatellite markers (SSRs) (Jahnke et al., 2015), revealed they are not fully distinct, sharing most of the alleles (F ST : 0.092) and possess a very similar genotypic richness (Ischia, R = 0.30; Bacoli, R = 0.37) (data not shown).

Experimental Design
After collection in the field, plant material was kept in darkened cooler containers filled with ambient seawater and rapidly transported to the indoor mesocosm facility at Stazione Zoologica Anton Dohrn (SZN, Naples, Italy). A detailed description of the experimental mesocosm system can be found in Ruocco et al. (2019b). Two plant fragments consisting of 15-20 connected shoots for each EU-and OL-plants were individually attached to the bottom of one plastic net basket (base 34 × 24 cm, height 10 cm) filled with coarse sediment. Twelve baskets for each site were arranged randomly in 12 glass aquaria (500 L) filled with natural seawater from a nearby-unpolluted area and following an orthogonal experimental design ( Figure 1B). In this way, plants from both locations were placed in the same tank and in triplicates per experimental treatment/condition. The four treatments were: Control (C), Temperature (T), Nutrients (N) and Nutrients plus Temperature (NT). Stressful environmental factors for our experiment were set according to previous mesocosm experiments and environmental features of sampling sites. Experimental temperature level was selected from the SST data recordings obtained from a MEDA (Monitoring and Environmental Data Unit) type buoy placed in Bagnoli, Gulf of Naples (IRM, RIMAR department, SZN) about three miles far from the sampling site, within the Gulf of Pozzuoli. In years 2017 and 2018, the maximum temperatures reached in August were 30 and 28 • C, respectively, whereas the average temperatures measured for the same month were 24.6 • C in 2017 and 25.3 • C in 2018 (Supplementary Figure 1). Thus, it was decided to expose seagrasses to a temperature of 30 • C, which is 4-5 degrees above the seasonal average, a sub-lethal temperature level for the species in the mid-term (Traboni et al., 2018). There was an initial pre-acclimation (1 week) into the mesocosms under the environmental conditions found at the time of plant collection (21 • C). Then temperature was gradually increased (0.5 • C day −1 ) in all tanks to reach the experimental control temperature of 24 • C, at warming rates naturally occurring in late spring/early summer according to MEDA recordings. Plants were allowed to acclimate to the new temperature for an additional week. Then, in T treatment tanks, temperature was gradually increased (0.5 • C day −1 ) until reaching 30 • C in a 2-week period ( Figure 1C). Light was provided by two LED lamps (max. noon irradiance = ca. 400 µmol photons m −2 s −1 above the canopy; 14:10 h light:dark photoperiod) in accordance to the natural levels measured in the field during plant sampling; for a detailed description of the LED system and settings see Ruocco et al. (2019b). Seawater was circulated using a 10,000 l/h self-priming pump, allowing continuous seawater circulation and filtration as well as a continuous replacement of water in the system. Additionally, within each aquarium, incoming seawater was spread through a diffuser in order to favor water mixing and to create a homogenous movement of water. Temperature in aquaria was controlled and maintained by cooler/heaters (Teko TK 2000). Salinity was measured daily in each aquarium using a WTW Cond 3310 portable conductivity meter and kept constant (within the range 37.3-37.7) by adding freshwater to compensate for evaporation. Aquaria were cleaned every week, and a 50% (filtered 45 µm and UV irradiated) seawater was also renewed to maintain seawater quality in the system.
According to the literature, NO 3 − , PO 4 −3 and NH 4 + levels follow different trends among seagrass meadows with low concentrations at canopy level (0.1 -4 µM) and higher concentration into sediments, where NH 4 + can also reach 180 µM Burkholder, 2000, 2002). Dynamics of nutrient concentration in our mesocosm system and frequency of nutrient addition was previously determined in a trial, where nutrient concentration was observed to decrease starting 3 days from the initial spike (data not shown). Hence, we decided to add nutrients weekly, to maintain the value almost constant during the experiment (Figure 1). To avoid a sudden increase, the stock solution (170 mM total nitrogen) was added during two consecutive days at the beginning of each week, for a 5-week experimental period (Figure 1). We aimed to reach a nominal nutrient concentration of 100 µM (total nitrogen) in the nutrient treatment tanks. Since nutrient concentration was assessed in the tanks at the end of each week, measured values reflected the weekly nutrient consumption in the experimental system. The stock nutrient solution was prepared using Osmocote R Pro (6 months release: 19% N -3.9% P -8.3% K, ICL Specialty Fertilizers) fertilizer pellets (Ravaglioli et al., 2017. Water analysis of dissolved organic nitrogen (DIN) confirmed nutrient enrichment in N and NT treatments (26.8 ± 4.0 µM) and normal levels in C ant T (1.7 ± 1.1 µM). According to the European Environmental Agency for the Mediterranean Sea, total nitrogen levels lower than 2.28 µM are considered "low" (class boundaries concentration is determined by the 80/20 percentiles of the DIN dataset for the years 2007 to 2012), while values higher than 26.0 µM are considered "high" and are indicative of eutrophic waters 1 .

Plant Traits
Plant samples for biochemical analysis were collected at the end of 5 weeks and photo-physiology measurements were performed every week to assess the status of plants. In order to reduce within-shoot and within-leaf variability of P. oceanica responses to warming (Ruocco et al., 2019a), parameters were measured in the middle portion of young fully developed leaves (second-and third-rank leaves). The aquarium was our true experimental unit, hence measurements performed on shoots of the same aquarium (i.e., 'pseudoreplicates') were averaged to obtain an independent replicated value. Therefore, the number of replicates used in all statistical tests was n = 3.

Carbon and Nitrogen Content
The leaf total carbon and nitrogen content was measured on a 6-cm segment collected starting 14 cm above the ligule of the third-rank leaf of four randomly sampled shoots per treatment. The analysis was carried out on epiphyte-free dried and homogenized tissue using an automatic elemental analyser FlashEA 1112 (Thermo Fisher Scientific, Waltham, MA, United States). Acetanilide was used as standard. Carbon and nitrogen content were expressed as a percentage of dry weight and the values were used to calculate the C:N ratio.

Photochemical Efficiency and Chlorophyll Content
Photochemical efficiency was estimated on the basal healthy and epiphytes-free part of leaves (10-15 cm from ligule) of the same shoots selected for morphological and physiological analysis. Chlorophyll a fluorescence measurements were conducted with a diving-PAM portable fluorometer (Walz, Germany). Rapid light curves (RLC) were performed after 4 h of illumination to calculate maximum electron transport rates (ETR m ) as an estimate of photosynthetic efficiency. During RLC, a saturation pulse was applied after 20 s incubation at increasing irradiance levels to determine the basal (F) and maximal light-adapted fluorescence (F m ). These parameters allowed for the calculation of the effective quantum yield (Y). Relative electron transport rate (rETR) at each irradiance level of RLC was calculated according to the formula: rETR = Y × E i × 0.84 × 0.5 (Beer et al., 2014), where Y is the effective quantum yield in the light-adapted state, E i is the incident irradiance (µmol photons m −2 s −1 ), 0.84 is the default average light absorptance of leaves (Ralph et al., 1998) and 0.5 is the fraction of photons absorbed by PSII in plants (Beer et al., 2014). rETR data (µmol electrons m −2 s −1 ) were plotted against incident irradiance and the resulting curve was fitted to the model described by Eilers and Peeters (1988) to obtain ETR m values.
Following light-adapted chlorophyll a fluorescence measurements, selected leaves were immediately frozen at −80 • C for subsequent pigment analysis. Leaf segments (1 cm 2 ) were homogenized using 80% acetone solution buffered with MgCO 3 to prevent acidification of the extract. Then, extracts were stored in the dark at 4 • C for 24 h and finally centrifuged (1000 g for 10 min at 4 • C). The absorbance of extracted pigments was read spectrophotometrically at 470 nm, 646 nm, 663 nm and 725 nm, using 3 mL quartz cuvettes. Total chlorophyll content (chlorophyll a and b concentrations) was calculated using the equations described by Lichtenthaler and Wellburn (1983) and expressed as µg cm −2 .

Leaf Morphology and Relative Growth Rate
Morphometric measurements (leaf length and width) were performed on the same shoots used for photo-physiological analyses and carbon/nitrogen content, using a ruler. The necrotic surface area of leaves was measured considering the portion of leaves with necrotic marks for each shoot and reported as percentage. The remaining surface of healthy tissue was measured and estimated as total photosynthetic surface (TPS, cm 2 shoot −1 ). Two weeks before the end of the experiment, five randomly selected shoots from both EU and OL were marked using the method described by Zieman (1974), i.e., piercing the boundary limit between the leaf and the ligule with a fine needle. Shoots were collected at the end of the experiment to calculate the leaf area of newly growth tissue. Relative growth rate (RGR) was estimated by the ratio between the newly produced leaf area (width × length) and the initial area measured before marking (cm 2 cm −2 day −1 ).

Carbohydrate Content
The content of carbohydrates was analyzed in leaf and rhizome tissues following the method described by Dubois et al. (1951). A portion of 6 cm from the third leaf (epiphytes removed) and a 3 cm of clean rhizome apex were collected from all treatments and plants locations (EU and OL). Leaf and rhizome samples were dried at 50 • C for 24 h and the extraction of soluble sugars and starch was performed using a phenol-sulphuric acid reaction. Soluble sugars (sucrose, fructose and glucose) were suspended from 50 mg of the ground tissue using subsequently heated 80% ethanol reactions (80 • C for 15 min). Samples were centrifuged (3000 rpm for 10 min) and then soluble sugars were extracted from the supernatant using 3% phenol and sulfuric acid. The absorbance of extract was read spectrophotometrically at 750 nm and 490 nm using a 1 mL quartz-glass cuvette. The analysis of starch content was performed adding 3 mL of NaOH to the remaining solid pellet after the extraction of soluble sugars and stored at 4 • C for 24 h. Total non-structural carbohydrates (TNC) were calculated according to Sørensen et al. (2018) and expressed on a dry weight basis (mg g −1 DW).

Plant mortality
Plant mortality was assessed as change in shoot number between the beginning (T1) and the end (T4) of the experiment. All shoots were counted and the net change in shoot number was calculated for OL and EU. For graphical representation, average shoot number of the three replicates/treatment was considered and expressed as percentage.

Statistical Analysis
Morphological data (TPS, necrotic area), physiological data (carbon and nitrogen content, photochemical efficiency and chlorophyll content, carbohydrate content), relative growth rate and the shoot mortality were analysed with 3-way analysis of variance (ANOVA) to detect significant differences between the responses of EU and OL plants to experimental treatments (see Supplementary Table 2). The model included plants (P, with two levels: EU and OL), temperature (T, with two levels: control and high) and nutrients (N, with two levels: control and high) as fixed factors. When significant differences were found for the factor plants, EU and OL plants were analyzed individually with 2-way ANOVA. In this case, the model was similar to the previous one, including temperature and nutrient as fixed factors (Table 1).
Before carrying out ANOVA analyses, normality and homoscedasticity were checked using the Shapiro-Wilk and Levene's tests and data subsequently transformed where necessary. Student-Newman-Keuls (SNK) post hoc test was used whenever significant differences (P < 0.05) among treatments were detected using the statistical package STATISTICA (StatSoft, Inc., v. 10). All variables measured for EU and OL were plotted in a Principal Component Analysis (PCA) to analyze differences in similarity patterns among treatments and plants with the software PAST v.3.03 (Hammer et al., 2001). A correlation matrix was performed to standardize measurement scales of different variables.

Carbon and Nitrogen Content and Ratio
Leaf nitrogen content was significantly higher in treatments with addition of nutrients (N and the combined NT) compared to control (C) and high temperature (T) treatment in both EU and OL (P < 0.01) (Figure 2A and Table 1). In contrast, carbon content was similar among treatments. Post hoc comparisons indicated significant differences in leaf nitrogen content among treatments, where NT (2.84%) and N (2.37%) were significantly greater than C in EU (NT vs. N, P < 0.05), but not in T treatment, which was similar to C conditions in both plant groups (T vs. C, P > 0.05). This variability in terms of nitrogen content was reflected in the Carbon:Nitrogen (C:N) ratio ( Figure 2B). C:N ratio decreased with the increase of nitrogen in treatments enriched with nutrients (−70% in N and −45% in NT). This pattern was similar for EU and OL, although EU showed a lower C:N ratio under NT treatment (12%) compared to OL (15%; P < 0.01). In the latter, C:N ratio in the NT treatment was similar to the N treatment.

Photochemical Efficiency and Chlorophyll Content
The maximum electron transport rate (ETR m ) differed between OL and EU plants (3-way ANOVA; Supplementary Table 2). ETR m was higher in EU than in OL plants (35.3 ± 2.2 vs. 24 ± 1.9 µmol electrons m −2 s −1 ; P < 0.01) (Figure 3A). In both OL and EU plants, photochemical efficiency decreased in N and NT treatments compared to controls (N vs. C −23%, NT vs. C −14% in EU; N vs. C −48%, NT vs. C −24% in OL), whereas the presence of high temperature alone did not influence photosynthetic performances. Total chlorophyll (Chl a + Chl b) content was significantly higher in EU than OL (+12%; P < 0.05). In contrast to OL plants, which maintained constant chlorophyll levels in all treatments, EU plants increased their content under nutrient addition (N, +30%) and temperature increases (T, +9%), whereas NT treatments showed similar levels to controls (Figure 3A).

Morphology and Relative Growth Rate
The highest percentages of necrotic area was recorded in treatments enriched with nutrients (i.e., N and NT) in both TABLE 1 | Results of two-way ANOVA analysis in OL and EU for factors "Nutrients" (N) and "Temperature" (T) for nitrogen content (%N), carbon and nitrogen ratio (C:N Ratio), Total chlorophyll, maximum electron transport rate (ETR m ), total non-structural carbohydrates in leaf (TNC leaf) and rhizome (TNC rhizome), % necrotic surface, total photosynthetic surface, relative growth rate and shoot mortality. Values in bold indicate significant differences (P < 0.05).
OL (33 and 27%, respectively) and EU plants (28 and 32%, respectively; P < 0.01) (Figure 3B). In contrast, the percentage of necrosis was lower in OL plants under T treatment, in respect to the control (−14%; P < 0.01). Total photosynthetic surface (TPS) measured in OL displayed an inverted pattern in respect to necrotic area, whereas EU followed a different trend especially for T and NT treatments, where total photosynthetic surface remained similar to control conditions even if necrosis increased (+23 and +13%, respectively; Figure 3B). The relative growth rate (RGR) showed a decrease under N (−58%) and T (−30%) in OL plants, although not statistically significant. Growth was affected by the combined treatment (NT) only in EU plants (Table 1), even if the post hoc comparison was not significant. Statistical analysis revealed no significant differences between plants and treatments ( Figure 3D and Supplementary Table 2).

Carbohydrate Content
Total non-structural carbohydrates (TNC) in leaves showed a different pattern compared to rhizomes ( Figure 3C). In leaves, TNC measured from OL plants decreased in the N (−32%) as well as in the NT treatment (−24%) compared to control (P < 0.01). In contrast, TNC decreased in EU plants for all treatments especially under N treatment (−70 %). At rhizome level, EU plants showed the greatest reduction in N treatment (−38%) followed by NT (−10%) ( Table 1). This is in contrast with OL plants, where no clear trend was observed. TNC were higher in rhizomes compared to leaves in both OL and EU plants (+29 and +45%, respectively).

Shoot Mortality
The net change in shoots measured in EU and OL plants significantly changed during the experiment (P < 0.01, 3-way ANOVA; Figure 4 and Table 1). Temperature was the main factor that affected shoot mortality along the experiment (P < 0.05). The highest mortality (−41.6%) was observed in EU plants under high temperature conditions (T treatment), followed by NT (−27.9%) and N (−23%). OL plants showed an overall lower mortality. As for EU plants, T induced the highest mortality (−15.7%), followed by NT (−15.6%) and N (−7%). Along the experiment, the highest mortality was observed between T2 and T4 (Supplementary Figure 2).

Principal Component Analysis
Considering EU plants, the first PCA axis (PC1), which explains 43.5% of the total variance, clearly placed treatments where nutrients were involved (N and NT) on the negative side of the axis, away from control (C) and temperature (T) treatments that are on the positive side ( Figure 5A). The second axis (PC2; 25.6% of total variance) segregated nutrient treatment from the combined stress treatment, being positively related with TNC measured in the rhizome and negatively related with the total chlorophylls. The PCA performed on OL plants ( Figure 5B) showed a different distribution pattern compared to EU plants.
The first axis (PC1; 40.8% of total variance) clearly separates temperature from control and the combined treatment. This axis was negatively related with N content (%N) and positively related with TNC measured in the leaf. The second axis (PC2; 20.1% of total variance) segregated nutrient treatment from the combined one even if in this case treatments appeared more widely distributed, with a positive relation with ETR m and a negative relation with TNC measured in the rhizome.

DISCUSSION
Being sessile organisms, seagrasses must cope with possible changes of local environmental conditions, modifying their morphological and mechanical traits, diverting energy for the maintenance of photosynthesis and changing their physiological traits (La Nafie et al., 2012Nafie et al., , 2013de los Santos et al., 2013;Marín-Guirao et al., 2016;Roca et al., 2016). Our results revealed that plants undergoing chronic cultural eutrophication (EU) are more sensitive to a further increase of nutrients, particularly when in presence of a temperature increase, than plants growing in oligotrophic water conditions (OL). OL and EU plants showed different morphological and physiological performances, which corroborates the role of local conditions in activating different strategies in response to environmental changes. We also found that temperature and nutrient treatments showed opposite effects when tested individually and an offset response when combined (Figure 6). We found antagonistic interactions resulting from the combination of chronic nutrient increase and temperature rise for chlorophyll content and total non-structural carbohydrates in rhizomes, only in plants growing in eutrophic conditions. Intraspecific variations in terms of resistance and performance are widely documented in seagrasses (e.g., Dattolo et al., 2014;Sandoval-Gil et al., 2014;Zhang et al., 2014;Marín-Guirao et al., 2016;Procaccini et al., 2017;Beca-Carretero et al., 2019;Tuya et al., 2019) The differential responses showed in our study further confirm that P. oceanica can have different levels of plasticity in response to changes of environmental conditions. According to our results, plants from nutrient enriched conditions, which are very fragile and FIGURE 3 | Endpoints measured in OL (oligotrophic) and EU (eutrophic) plants after 5 weeks of the exposure to high temperature and nutrient enrichment. (A) Chlorophyll content and maximum electron transport rate (ETR m ); (B) Necrotic area and photosynthetic surface; (C) Total non-structural carbohydrates (TNC) measured in leaves and rhizomes; (D) Relative growth rate. Asterisks indicate differences between OL and EU plants (ANOVA 3-way); (+) indicates factors where significant differences were found between levels of factor after ANOVA 2-way and post hoc test. Different lowercase letters indicate the groups of homogeneous means obtained in the post hoc (P < 0.05).
FIGURE 4 | Shoot mortality measured as percentage comparing the mortality between the beginning (T1) and the end of the experiment (T4) in OL and EU plants exposed to nutrient (N), temperature (T) and their combination (NT). Asterisks indicate differences between OL and EU plants (ANOVA 3-way); (+) indicates factors where significant differences were found between levels of factor after ANOVA 2-way and post hoc test. Different lowercase letters indicate the groups of homogeneous means obtained in the post hoc (P < 0.05).
break off easily, may be at higher risk, compared to plants from more oligotrophic conditions. It has also been observed by studies on Z. noltei that nutrient enrichment is able to reduce the strength of leaves (La Nafie et al., 2012). Although plants undergoing a particular stress factor could acclimate to the environment or even improve their resistance to a subsequent unfavorable pressure (Alexieva et al., 2003), individuals growing under cultural eutrophication conditions could already be at the limit of their biological capacity. Thus, the high energetic costs for nutrient excess acclimation might have compromised their tolerance to a further nutrient enrichment and to thermal stress. OL plants showed higher C:N values than EU plants, suggesting the presence of specific nutrient-balancing strategies, that may be related to the different nutrient levels chronically experienced by two populations in natural conditions. The activation of contrasting acclimation mechanisms to face high nutrient loads has been already demonstrated in P. oceanica exposed to chronic or pulse enrichments in a long-term field experiment . In our experiment, plants living in eutrophic conditions (EU plants) treated with high nutrient concentrations, regardless of temperature, increased the percentage of nitrogen content followed by a decrease of TNC reserves in leaves and rhizomes. A similar response was observed in OL plants, but only for leaves. The reduction of carbohydrate reserves is thought to be related to the assimilation of NH 4 + , which needs C skeletons to create new amino acids, an energetically costly process (Alcoverro et al., 2000). This strategy is more evident in EU plants, where TNC in rhizomes were higher than in leaves, suggesting translocation processes as a response to nutrient enrichment. A possible explanation for the observed TNC reduction only in leaves of OL plants, could be the lack of an active mechanism to cope with nutrient excess due to their low natural nutrient exposure. From a different perspective, since carbohydrate reserves allow P. oceanica to maintain leaf growth and to assimilate C and N from the environment (Invers et al., 2004), the absence of carbohydrate reallocation could also be considered as an alternative strategy: as OL plants appeared healthier than the EU ones and their rhizome stronger, probably OL plants reduced carbon reserves in the leaves to cope with nutrient excess, maintaining constant reserves at rhizome level. In contrast, EU rhizomes were weaker due to chronic cultural eutrophication at their natural environment and used more carbohydrates to deal with nutrient assimilation.
Carbohydrate reserves in temperate species are fundamental for survival under low-light conditions, especially during winter (Alcoverro et al., 2001;Soissons et al., 2018). A proportion of the carbohydrates used to produce amino acids during nutrient assimilation will certainly go to the production of other compounds such as chlorophylls. In EU plants under only nutrient treatment, total chlorophyll content increased significantly compared to control and the combined treatment, suggesting that the accumulated amino acids, such as glutamate, were used toward the production of chlorophyll molecules. Similar effects have been reported for seagrasses under low pH conditions and nutrient additions, where the presence of nutrients implied higher concentrations of chlorophylls compared to control (Roca et al., 2016;Ravaglioli et al., 2017). The change in chlorophyll concentrations was not observed in OL plants, indicating once again the lack of response compared to EU plants, which may produce more chlorophylls as a strategy to capture more light in eutrophic (reduced) light environments. Despite a greater concentration of chlorophylls due to the increased nutrient treatment, the photochemical performance (ETR m ) decreased in both EU and OL plants. However, as suggested above, the reduction in light availability may be responsible for this effect, as PSII modulates quickly according to light changes (Beer et al., 2014) and this may not be reflected in biochemical variables. It is important to highlight that ETR m values in EU plants were overall higher than OL plants.
There is evidence that P. oceanica responses to warming vary among plants experiencing different thermal conditions in their natural range and some are able to recover after short-term warming events (Marín-Guirao et al., 2016. The high temperature treatment in the present study could have been below the thermal tolerance threshold of the sampled individuals, as not many physiological parameters responded to this factor. Optimum growth temperature for this temperate species range from 11.5 to 26 • C (Olsen et al., 2012). However, this optimal range probably depends on other co-occurring environmental factors, such as nutrient availability, leaf senescence and nutrient partitioning within plant tissues (Ontoria et al., 2019a). Indeed, photochemical efficiency and carbohydrate reserves of OL plants were not affected by warming in this study. This is in contrast with EU plants, where high temperature led to the reduction of TNC, which reflects the metabolic adjustment of plants accelerating leaf senescence. Different to other parameters measured after 5-week exposure, the shoot mortality under high temperature in EU plants was higher (−41.6%) than OL plants (−15.7%). A possible explanation is that plants exposed to high temperature tried to acclimate to stressful conditions diverting  N), total chlorophyll, maximum electron transport rate (ETR m ), total non-structural carbohydrates in leaf (TNC leaf) and rhizome (TNC rhizome), % necrotic area, total photosynthetic surface (TPS) and relative growth rate (RGR). Different colored circles and symbols refer to different treatments (• = Control; = Nutrients; = Temperature; = Nutrients + Temperature).
energy for the maintenance of high metabolic rates (Olsen et al., 2012). Thus, plant growth decreased due to high-energy costs . Accordingly, we found lower relative growth rates (RGR) in EU plants and reductions of RGR in nutrient and temperature treatments for OL plants. This is in agreement with previous studies performed on P. oceanica, where shoot survival decreased after exposure to high temperature  and high nutrient (Ceccherelli et al., 2018) concentrations, respectively. While the physiological changes may be subject to a rapid regulation, the observed mortality may indicate that temperature could be a determinant factor in the long-term survival of the whole ramet.
Few multi-factorial experiments have been performed on P. oceanica, including some assessing nutrient (ammonium) input and temperature (Gera et al., 2013;Ontoria et al., 2019a). In these studies, different levels of synergism have been found between stress factors, influencing negatively the survival, photosynthesis and plant growth. In contrast, we found antagonistic interactions resulting from the combination of chronic nutrient increase and temperature rise for pigment content and TNC in EU and for growth in OL. Indeed, both EU and OL plants exposed to the combined stressors showed similar leaf nitrogen contents as plants that were only exposed to increased nutrients, with no change in TNC at rhizome level. This suggests that despite nutrient excess negatively affected plant performance, the simultaneous exposure to increased temperature accelerated metabolic rates to cope with the impact induced by nitrogen assimilation. These results are in agreement with those obtained by Egea et al. (2018) on Cymodocea nodosa plants, where the simultaneous exposure to NH 4 + enrichment and the high temperature had an antagonistic response. In contrast, Guerrero-Meseguer et al. (2020), found a synergistic effect of temperature increase in P. oceanica seedlings development when occurring concomitantly with other stressors, such as seed burial and grazing. In our experiment, plants were adversely affected by the combined treatment although the overall evidence showed a less negative effect in comparison to single nutrient and temperature treatments for the physiological parameters. Similarly, in Zostera capensis it was shown that nutrient enrichment and warming had a limited interaction, and that eutrophication was a stronger stressor (Mvungi and Pillay, 2019). It would be interesting to assess if P. oceanica seedlings would show similar response since Artika et al. (2020) found that in E. acoroides seedlings rely more on seed nutrient resources and are less affected by experimental nutrient addition.
Overall, our results suggest that ongoing eutrophic conditions in the Gulf of Pozzuoli are weakening the local P. oceanica meadow, since these plants showed the greatest percentage of mortality, particularly under increased temperature. Our findings do not seem to be a general rule for seagrasses, since Connolly et al. (2018) found the opposite trend in Z. muelleri. In that study, plants growing in more disturbed sites were more resilient to further disturbance, although with a lower genotypic diversity, suggesting a genotypic selection. In our analysis, the two studied meadows are only slightly distinct and did not show differences in genetic or genotypic variability, indicating that factors other than genotypic selection can explain our results. The data we presented indicate that eutrophication is likely inducing severe effects on the local seagrass meadows, which have activated physiological strategies to cope with excess of nutrients. However, the lack of replication of relative EU and OL populations of P. oceanica in our experiment precludes the extrapolation of our findings to other than the studied populations. The strategy adopted by these plants probably implied large energetic costs affecting their present and future persistence. Considering the high level of cultural eutrophication of a large part of the coastlines, and taking into account the predicted increase of SST and heat waves frequency, our results highlight the intraspecific vulnerability to environmental changes of seagrass meadows and the importance of performing studies at a local scale. Unraveling intraspecific strategies and the role of local acclimation/adaptation in the response to multiple stressors could be crucial for seagrass conservation strategies under a climate change scenario (Duarte et al., 2018). It becomes fundamental to perform multifactorial experiments on seagrass populations from different environments, in order to understand their resilience to future local and global environmental changes and to support management strategies.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JP, AS-S, LM-G, and GP conceived and designed the experiments. SH performed initial field assessment. MR contributed to the analysis of carbohydrates and assessment of plant morphology. JP and AS-S performed the experiment, all the other analyses and drafted the manuscript, and all authors reviewed it critically.

ACKNOWLEDGMENTS
Special thanks are given to: Emanuela Dattolo, Hung Manh Nguyen, and Ludovica Pedicini for their help during the periodical sample collection and the mesocosm system maintenance; Giovanni De Martino for helping in the mesocosm maintenance. We are grateful to the IRM and MAA units of the RIMAR Department (SZN) for the sampling of seagrass ramets and for carbon and nitrogen analyses. We also thank the "Progetto infrastrutturale EMSO-MedIT, data management of the MEDA-A Bagnoli (IRM unit of the RIMAR Department SZN).