Microbial plankton responses to multiple environmental drivers in marine ecosystems with different phosphorus limitation degrees

natural microbial communities from Mediterranean Sea and Atlantic Ocean exposed to temperature, nutrients, CO 2 , and UVR interactions. At thedecadescale,wefoundthatmoreintenseandfrequent(andlongerlasting)Saharandustinputs (andmarine heatwaves) were only coupled with reduced phytoplankton biomass production. When microbial communities wereconcurrentlyexposedtofuturetemperature,CO 2 ,nutrient,andUVRconditions(i


H I G H L I G H T S
• Increased intensity and frequency (by 3-fold) of Saharan dust inputs since 1985 • Heatwaves are more frequent, intense, and lasting (up to 6-fold) since 1985.• Multi-driver effects were evaluated on phosphorus (P)-limited microbial plankton.• The ecosystem behaved as a C-source under a more P-limited multi-driver's scenario.• A more P-limited multi-driver's scenario weakened the algae-bacteria interaction.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Understanding how global change influences trophic interactions, and ultimately the carbon (C) budget of the ecosystem poses an ongoing challenge in ecology.Over recent decades, global-change models, time series, and manipulative experiments have tracked the potential ramifications of global environmental alterations.Such approaches have evidenced that shifts in environmental drivers (e.g.temperature, nutrients) are affecting the physiology of many organisms (Wiltshire and Manly, 2004;González-Olalla et al., 2017), their metabolism (Bopp et al., 2013;Gao et al., 2018), and their interactions between trophic levels (Joint et al., 2002;Carrillo et al., 2015a;Durán et al., 2016).
Temperature and nutrient availability are two dominant environmental drivers controlling the pathways and rates of materials and energy moving through the ecosystem (Kleiber, 1975;Kaspari, 2012).Nutrient uptake is mediated by imbalances between an organism's demand of resources and their relative availability in the environment (Sterner and Elser, 2002), whereas temperature is a master variable controlling all biological activity in ecosystems through its effect on metabolic rates (Gillooly et al., 2001).Advances in ecological stoichiometry (Sterner and Elser, 2002) and metabolic ecology (Brown et al., 2004;Schramski et al., 2016) have developed a framework that enables researchers to quantify how nutrients and temperature interact to control ecological dynamics (Allen and Gillooly, 2009;Billings and Ballantyne, 2013).Under this conceptual framework, a review by Cross et al. (2015) proposes that rising temperature accelerates the metabolic rates regardless of nutrient availability because the first driver masks or counteracts the effect of the second one.This proposal is supported by recent findings reflecting a strong stimulatory warming effect on key processes such as primary (Lee and Kang, 2020;Schubert et al., 2019) and bacterial (Scofield et al., 2015) production or on microbial communities' β-diversity (Ren et al., 2017) regardless of nutrient status.Contrary to the proposal of Cross et al. (2015), ocean-scale evidence supports the contention that nutrient availability controls the temperature dependence of metabolism and thus temperature may be far less decisive than expected under nutrient limitation (Marañón et al., 2014).In fact, results with different phylogenetic phytoplankton taxa have shown that the stimulatory effect of temperature on metabolic processes is reduced (or even canceled) as the nutrient limitation intensifies (Marañón et al., 2018), likely due to diffusion-limited nutrient uptake rather than a limited rate of internal enzymatic processes (Thingstad and Aksnes, 2019).The finding by Marañón et al. (2018) is consistent with results reported elsewhere on respiration rates (Garnier et al., 2017;Tabi et al., 2019) as well as on bacterial C demands to fuel their metabolism (Huete-Stauffer et al., 2018) in natural communities.In view of all these findings, two questions arise: Could the degree of natural nutrient limitation ultimately govern the magnitude and direction (from stimulation to inhibition) of the impacts of warming and nutrient pulses on ecosystems?Could such impacts be modified when additional global-change drivers than temperature and nutrients are altered?
The two above questions are ecologically relevant and timely because: (1) previous experimental findings with natural microbial plankton communities from the Atlantic Ocean evidenced that as the oligotrophy accentuates, bacterial production (BP) (and respiration) increased up to 700-800%, whereas primary production (PP) decreased ~60% (Marañón et al., 2010).The higher ability of heterotrophic bacteria to augment their productivity and C expenditure (i.e. its bacterial carbon demand, BCD) has also been observed in the Eastern Mediterranean Sea (Herut et al., 2005;Pulido-Villena et al., 2014), as well as in the subtropical North Atlantic Ocean (Moore et al., 2008).Therefore, considering a higher BCD and ability of bacteria (compared to phytoplankton) to acquire nutrients at lower ambient concentrations (Cotner and Biddanda, 2002),we predict a potential accentuation of the heterotrophic processes and, consequently, major shifts in the ecosystem C budget (i.e.primary production:community respiration ratio, PP:CR), and phytoplankton-bacteria commensalistic interaction [BCD:excreted organic carbon (EOC)].(2) Together with temperature and nutrient changes mentioned above, the alteration of a wide range of environmental drivers tends to occur simultaneously.For example, anthropogenic activities since the Industrial Revolution have influenced global climate through the release of CO 2 into the atmosphere, increasing its concentration by more than 40% to date (IPCC, 2013).Also, the ongoing stratification of the water column resulting from elevated sea-surface temperatures traps the organisms in shallower upper mixed layers, exposing them to higher ultraviolet radiation levels (UVR) (Gao et al., 2012;Kwiatkowski et al., 2020).Thus, recent research is demanding more complex studies to improve our understanding and predictions regarding how organisms and their interactions respond to multiple environmental drivers because such interactions trigger nonadditive effects, misguiding predictions based on single or doubledriver studies (Brennan and Collins, 2015;Boyd et al., 2016;Galic et al., 2018;Villar-Argaiz et al., 2018;Feng et al., 2020).
The ocean, the dominant biome on Earth, is limited primarily by nutrients, mainly P and N (Moore et al., 2013).Recent modeling results predict that an intensification of the nutrient limitation coupled with warmer, stratified, and more acidic environments may reduce global total primary production (PP) up ~9% by 2100 (Kwiatkowski et al., 2020).However, a recent review by Van de Waal and Litchman (2020) proposes that these complex environmental conditions may reduce general metabolic and photosynthesis costs, allowing energy reallocation to the acquisition of increasing limiting nutrients, thereby dampening the reduction of PP.Because available information shows contrasting evidence, we stress the need for experimental approaches that not only quantify how multiple interacting environmental drivers alter plankton metabolism and trophic interaction between phytoplankton and bacteria, but also determine whether such effects are modulated by in situ nutrient limitation degree.
To help fill this gap in knowledge, we experimentally manipulated UVR, pCO 2 , dust-derived nutrients, and temperatures to simulate conditions for the year 2100.For this we used a collapsed design, and tested the effect of these factors on C-budget and commensalistic phytoplankton-bacteria interaction in two marine microbial communities that differed in the degree of P limitation.By using this approach, we seek to achieve a more comprehensive understanding about how planktonic communities may respond to future complex environmental scenarios.

Remote-sensing approach
Daily, weekly, and monthly area-averaged aerosol index (AI), biomass of phytoplankton taxa (Bacillariophyceae, Chlorophyceae, Cyanophyceae, and Haptophyceae), chlorophyll a (Chl a), and vertical attenuation coefficient of photosynthetically active radiation with depth (as Kd PAR at 490 nm) data, respectively, for both marine study areas were downloaded from the Geospatial Interactive Online Visualization and Analysis Infrastructure (Giovanni v. 4.33) Gregg et al., 2017), whereas the data for the phytoplankton groups were validated against surfacelayer observations of phytoplankton group abundances (Global mean difference model-data: < 15%; Gregg and Casey, 2007) and phytoplankton pigments (~20 different types, and ~4000 observations) derived from high-performance liquid chromatography (global mean difference model data ≤ 6%; Hirata et al., 2011).
From these datasets, we calculated the mean seasonal value, i.e. for spring, summer, fall, and winter.For AI, we considered only values >0.5 because they constitute an intense event of dust deposition and are the main contributors to the total annual AI (Bullejos et al., 2010).We also estimated the annual frequency of these events as the total number of AI >0.5 per year.For phytoplankton taxa, we calculated their proportion (%) of the total community biomass as the ratio of the biomass of each specific taxon to the total phytoplankton community biomass.
The frequency, duration, mean, maximum, and cumulative intensity (in °C) of marine heatwaves (MHW) in the south-western Mediterranean Sea and north-eastern Atlantic Ocean, for the period 1985-2018, was examined using data downloaded from the Hadley Centre Sea Ice and Sea Surface Temperature v1.1 (Rayner et al., 2003).An MHW event was defined as the temperature that exceeded the 90th percentile of the seasonally varying threshold for at least 5 consecutive days (Hobday et al., 2016).Thus, the MHW frequency is the total number of events per year, and the cumulative intensity of MHW is the cumulative number of °C per year that temperature exceeded the seasonally varying threshold.The duration of this event is the number of days per year with MHW, and the severity is expressed as the difference between mean and maximum MHW intensity registered within each year.

Experimental approach
Coastal waters from the Mediterranean Sea and Atlantic Ocean are transparent and oligotrophic (Prieto et al., 1999;Ramírez et al., 2005) (Fig. S1A, B; Table S1) ecosystems with low biomass production (Chl a range 0.20-2.85μg L −1 ; Fig. S1).These study areas are also subjected to recurrent nutrient inputs of inorganic phosphorus (P) [although also nitrogen (N) and trace metals] through Saharan dust-deposition events (Gallisai et al., 2014), and a rising frequency of heatwaves over recent decades (see the Results section).
Surface-water samples (0-10 m deep, upper mixed layers) were collected on 18 and 25 July 2018 in two stations, 36°38′13.92″N-4°12′ 19.80″W (Mediterranean Sea) and 36°35′16.8″N-6°27′27.72″W(Atlantic Ocean), respectively, aboard the R/V UCÁDIZ using an automated CTD-rosette equipped with 10-L Niskin bottles.The samples were prescreened through 180-μm Nitex mesh to remove mesozooplankton, placed in 8-L acid-washed opaque bottles (40-L in total) and transported in darkness in temperature controlled containers to the laboratory (between 1 and 4 h away from the sampling sites).At the laboratory, water samples were maintained in darkness and at the sampling temperature overnight before to being used in the experiments.The following day (19 and 26 July 2018, for the Mediterranean Sea and Atlantic Ocean, respectively; hereafter initial time), plankton communities were placed in microcosms (600-mL quartz vessels), submerged in 25-L rectangular glass aquaria (see below), and exposed for 6 days (19-25 July for the Mediterranean Sea and 26 July-1 August 2018 for the Atlantic Ocean) to the conditions described below.
Four abiotic factors (temperature, nutrients, UVR, and CO 2 ) have been identified as the main regulators controlling the metabolism of microbial food and the commensalistic phytoplankton-bacteria relationship in oligotrophic ecosystems (Cabrerizo et al., 2017;Gao et al., 2018;Goudie, 2019).From the remote-sensing approach, we contextualized the biotic/abiotic environment of our study areas, and identified temperature (using MHWs as a proxy) and nutrient inputs (using AI as a proxy) as major/dominant drivers.Hence we experimentally evaluated its impact on isolation and combined under current and future UVR and CO 2 environmental conditions.For this, we used a collapsed experimental design (in triplicate; Fig. 1), as proposed by Boyd et al. (2016), which relies on prior identification of the dominant driver (s) and then collapses (groups together) the non-dominant drivers.The collapsed design allows the measurement of the relative importance of the dominant driver(s), and the interacting effect of the all drivers considered, whether dominant or not.Thus, it enables the researcher to make powerful inferences concerning the manner in which future complex conditions may alter community responses, whereas such deductions could be unmanageable (and even impractical) when a full or fractional factorial design is used.
Considering that temperature and nutrients mediate physiological/ metabolic processes (Brown et al., 2004;Sterner and Elser, 2002), and that both drivers have been altered in our study areas over recent decades (see the Results section), we maintained temperature and nutrients separate (as the dominant drivers), whereas the other two, i.e.UVR and CO 2 , were grouped together (i.e.collapsed; Boyd et al., 2016) into a combined scenario (Fig. 1).For temperature, we had two levels for each ecosystem: 'in situ' = 20.5 °C (for the Mediterranean Sea) and 20 °C (for the Atlantic Ocean); and 'warming' = in situ + 3 °C.The increased temperature (T) simulates the worst-case scenario predicted by 2100 (IPCC, 2019).For nutrients, we also had two levels: 'ambient' (unmanipulated) = nutrient concentrations at the time of sampling; and 'enr' (enriched) = + 29.03 μM N and + 1 μM P respect in situ conditions.This nutrient enrichment with inorganic N (as NH 4 NO 3 ) and P (as NaH 2 PO 4 ) was intended to mimic (not to counteract) the N and P inputs after an intense Saharan dustdeposition event (≈ 1 μM P at a molar TN:TP ratio = 30) on surface waters in Mediterranean ecosystems (Morales-Baquero et al., 2006), these inputs having been used in previous experimental studies by our group with microbial plankton communities (Cabrerizo et al., 2016;González-Olalla et al., 2017).Finally, for the collapsed factor, we had two levels: 'Current' = UVR and CO 2 (unmodified levels of both factors) and 'Future' = UVR and CO 2 according to predictions for the year 2100 (doubled levels of both factors; Fig. 1).The experimental UVR intensification simulated a 2-fold increase in surface UVR irradiance with respect to the amount currently received by phytoplankton communities in surface waters.By so implementing this UVR increase, we sought to simulate future conditions due to the ongoing stratification of the water bodies and subsequent shallower upper mixed layers (Kwiatkowski et al., 2020).The future pCO 2 scenario was simulated by bubbling both communities with CO 2 -enriched air at 750 ppm (Air Liquide S. A., París, France), mimicking atmospheric concentrations expected by 2100 (IPCC, 2013).

Incubation system
The experiment was conducted in the laboratory of simulation and modeling of aquatic ecosystems (Aqualab) of the Institute of Water Research of the University of Granada.The laboratory has a glassaquarium system with independent lines for controlling temperature through refrigerators (Teco® SRL tank TK2000, Ravenna, Italy) and solar radiation through light-emitting diode lamps (nano-LED light v.2.0; BLAU aquaristic, Barcelona, Spain) for photosynthetically active radiation, and Q-Pannel UVA-340 (Miami, FL, USA) tubes, which simulate the UVR portion of the spectrum from 365 nm down to 295 nm (peak emission at 340 nm), for UVR.To elicit realistic planktonic responses with respect to those in the natural environment, communities were exposed to constant mean irradiances of 100 μmol photons m −2 s −1 for photosynthetically active radiation in current and future environmental scenarios, and 2.42 (for current scenario) and ~5 (for future scenario) W m −2 for UVR over diel light-dark (16 h:8 h) cycles.2.42 W m −2 mimics mean integrated (0-10 m depth) underwater UVR values measured in the study area at noon during summer (i.e.2.96 W m −2 ; Cabrerizo et al. unpub.data).

Physico-chemical and biological characterization of the water column
Vertical profiles of temperature, salinity, dissolved oxygen, Chl a, and underwater photosynthetically active radiation were done using a CTD Seabird SBE 25 (Seabird Scientific Inc., CA, USA).
The response variables presented below were quantified at the initial time (i.e.19 and 26 July 2018, for Mediterranean Sea and Atlantic Ocean, respectively) and at the end of the exposure period (i.e. 26 July and 1 August 2019, for Mediterranean Sea and Atlantic Ocean, respectively) (Fig. 1).

Inorganic nutrients
Three hundred milliliters (in triplicate at the initial time and for each environmental scenario in both study areas) samples were used to determine total nitrogen (TN) and phosphorus (TP), and total dissolved nitrogen (TDN) and phosphorus (TDP) concentrations using standard protocols (N, Koroleff, 1977;P, APHA, 2017).

Chlorophyll a (Chl a)
For chlorophyll measurements, onetwo hundred fifty millilieters sample (in triplicate at the initial time and for each environmental scenario in both study areas) were immediately filtered onto GF/F glass fiber filters (Whatman Inc., UK) and stored at −20 °C until analysed.Before analysis, filters were thawed, placed in vials with 90% acetone and incubated in darkness at 4 °C for 24 h (Porra, 2002).Samples were measured at an excitation wavelength of 460 nm and emission at 670 nm, using a fluorometer (Perkin Elmer, model LD55, USA) routinely calibrated using a pure Chl a standard (Sigma-Aldrich, UK).

Primary production (PP)
PP was measured using the radio-labelled 14 C-incorporation technique (Steemann, 1952).Firstly, forty milliliters samples (triplicate clear bottles and one dark per scenario) were placed into Teflon FEP narrow-mouth bottles (Nalgene®, USA), inoculated with 5 μCi of labelled bicarbonate and incubated for 4 h in the same environmental scenario/aquarium from which they came.After this period, particulate PP was determined by filtering 35 mL of a subsample through GF/B glass fiber filters (Whatman Inc., UK).A subsample of 4 mL of the filtrate was also directly collected in scintillation vials to assess the excreted organic carbon by phytoplankton (EOC).Filters and filtrates were exposed to acid fumes for 24 h by adding 1 N HCl (2%) to eliminate the nonassimilated 14 C, and the samples were measured using a scintillation counter (Beckman LS-6000 TA, USA).Dark values were subtracted Fig. 1.Graphical scheme of the collapsed experimental design in which two microbial plankton communities were exposed to current and future scenarios of temperature, nutrients enrichment, carbon dioxide concentration (CO 2 ), and ultraviolet radiation (UVR).Note that UVR and CO 2 acted always together whereas that the other two drivers acted either in isolation or simultaneously in current and future scenarios (see material and methods for a detailed description).
from the corresponding light values.The total CO 2 in seawater samples was calculated from pH and alkalinity measurements performed with an automatic potentiometric titrator (Titrando 905, Metrohm Inc., USA).
Total PP was calculated as the sum of particulate PP and EOC.From this compartment, and together with the total community respiration (see description below), we calculated the C-budget of the microbial community (i.e.PP:R).

Bacterial production (BP)
BP was measured using the radio-labelled thymidine-incorporation technique (Fuhrman and Azam, 1982).Samples of one and a half mL (3 live replicates +2 blanks per scenario) were placed in acid-cleaned sterilized tubes, inoculated with 3 H-thymidine (Perkin Elmer, USA) to a final saturating concentration of 22 nM, and incubated for 1 h in darkness in the same environmental scenario/aquarium from which they came.After the incubation period, in a laboratory adjacent to the housing the aquarium system, samples were centrifuged at 16000 ×g, rinsed twice with cold trichloroacetic (5% final concentration), and measured using a scintillation counter (Beckman LS-6000 TA, USA).The data for blanks (killed controls) were subtracted from the corresponding values for live replicates.Finally, the correction factors, 1 × 10 18 cells mol −1 of thymidine and 2 × 10 −14 g C cell −1 (Bell, 1993;Lee and Fuhrman, 1987), were used to convert the incorporated radiotracer into BP expressed in terms of C.

Community (CR) and bacterial (BR) respiration
Raw water (for CR) and GF/F (Whatman Inc., UK) filtrated (for BR) forty milliliters samples were placed into Teflon FEP narrow-mouth bottles without bubbles (Nalgene®, USA).Each Teflon bottle equipped with an oxygen sensor-spot (SP-Pst3-NAU-D5-YPO, Presens GmbH, Germany) was incubated following the same procedure as mentioned above for PP and BP, and the oxygen concentration was measured during 24 h daily cycles using a Fibox Oxy-mini 10 optode-probe oxygen transmitter (Presens GmbH, Germany).Prior to experimentation, the transmitter was calibrated using a two-point procedure (0 and 100% oxygen concentration) together with data of atmospheric pressure and temperature.Respiration rates were calculated as the difference between oxygen concentrations measured at the beginning and at the end of the incubation period, and converted into C units assuming a respiratory quotient of 1 (del Giorgio and Cole, 1998).

Bacterial carbon demand (BCD)
The BCD was calculated as the sum of BP plus BR.

Commensalistic interaction
The coupling between phytoplankton and bacteria has been defined as the capacity of EOC by phytoplankton to support the BCD (Morán et al., 2002).As EOC is the C source preferentially used by bacteria to fuel their metabolism and ultimately their growth, we quantified the strength of the phytoplankton-bacteria coupling as the ratio (in percentage) between BCD and EOC.Thus, a BDC:EOC ratio > 100 (or 〈100) indicates a decoupling (or coupling) between the two trophic levels, i.e.BCD exceeds EOC (or EOC exceeds BCD).

Data and statistical analyses
The individual and interactive effects of the drivers considered on the C budget and commensalistic interaction were assessed using ln response ratios (LRR) as follows: where Current refers to the samples exposed to unmanipulated conditions of UVR, CO 2 , nutrients, and temperature; Current enr or Current warming and Future enr or Future warming refer to the samples exposed to unmanipulated (Current) and increased levels of UVR and CO 2 (Future) under enriched nutrient ( enr ) or warming ( warming ) conditions.Future enr+warming represent samples exposed to conditions in which all drivers were manipulated (i.e.increased).LRR values > or < 0 indicate a stimulation or inhibition of the processes assessed, respectively, and LRR = 0 signifies an absence of change.Commensalistic interaction values were multiplied by −1 to ease the biological comparison with the C budget.All data are reported as mean values and standard deviations or 95% confidence intervals, whereas a pooled standard deviation for equal sample size was used in the LRR of the C budget and commensalistic interaction.Linear and non-linear regression models were fitted for AI, Chl a and MHW data vs. time in each study area.A nested analysis of the variance (ANOVA), with the scenario (i.e.Current enr , Current warming , Future enr , Future warming , and Future enr+warming with respect to Current; See Supplementary Information, Figs.S3, S4) nested to the area (i.e. the Mediterranean Sea and the Atlantic Ocean), was used to test the effects of scenario on LRR C-budget and commensalistic interaction depending on the area.A least significant differences (LSD) post hoc test was used to determine differences between specific scenarios and areas.Student's t-tests were used to identify significant differences between TDN:TDP ratios (i.e.degree of P limitation) of the two marine areas.The statistical models were validated by checking the assumptions of normal distribution and homoscedasticity of residuals by Shapiro-Wilk's and Levene's tests, respectively.

Results
AI intensity and frequency, a proxy of Saharan dust-deposition events, significantly increased in both areas, being up to 3-fold more frequent over the 1985-2019 period (Fig. 2A, B; Table S2).Over the same period, the frequency and duration of MHWs significantly increased by 4-fold, cumulative intensity by 6-fold, and severity by 2-fold (Fig. 2C-F, Table S2).The seasonal Chl a and Kd PAR values varied greatly over the last two decades (Fig. S1A, B), registering maximum values during spring/summer (2.5-3 μg L −1 for Chl a; 0.3-0.5 m −1 for Kd PAR ) and minimums during winter (< 0.5 μg L −1 for Chl a; < 0.10 m −1 for Kd PAR ) in both areas over the last 35 years (Fig. S1A, B).The variations in Chl a concentrations were coupled with shifts in the phytoplankton community structure, from a dominance by Cyano-/ Chlorophyceae during summer to a dominance by Bacillariophyceae in winter (Fig. S2A, B).
The increases registered in AI intensity and in the MHW duration over the 1985-2019 period were associated with reductions in Chl a (Fig. 3).We found a significant negative relationship between Chl a and AI intensity in both areas (Mediterranean Sea: R 2 = 0.32, F 7.14 , p < 0.01; Atlantic Ocean: R 2 = 0.30, F 5.13 , p = 0.05; Fig. 3A), and a similar response pattern when Chl a was compared with MHWs duration, although it was significant only for the Mediterranean Sea (R 2 = 0.22, F 3.85 , p < 0.05; Fig. 3B).
Chl a concentrations were < 1 μg L −1 in both study areas at the sampling moment, and the communities were dominated by diatoms (mainly Nitzschia sp.).Dissolved nutrients reached similar concentrations in both areas for TDN, ranging between 34.88 (±0.08) and 30.12 (±0.03) μM for the Mediterranean Sea and Atlantic Ocean, respectively.By contrast, they differed by 2-fold for TDP concentration (mean values of 0.47 ± 0.04 and 0.23 ± 0.03 μM, respectively).These differences in the N and P concentrations translated into significantly weaker P limitation (in terms of the TDN:TDP ratio) in the Mediterranean as compared to the Atlantic (t-test = −4.88,p < 0.01; Table S1).The experimental nutrient addition reduced, but did not eliminate, the in situ Plimitation; that is, TDN:TDP ratios changed from 74.86 (±6.46) and 133.74 (±4.56) (before the nutrient pulse) to 43.55 (±1.17) and 48.23 (±1.18) (after the nutrient pulse) for the Mediterranean Sea and the Atlantic Ocean, respectively.Fig. 4 shows the effect size, calculated as LRR from C-budget (i.e.PP: R ratio) and commensalistic interaction (i.e.BCD:EOC ratio) values (Figs.S3, S4), of all environmental scenarios tested in two microbial plankton communities from the Mediterranean Sea and Atlantic Ocean areas.Current enr and Current warming had a stimulatory effect on C-budget (i.e.increasing autotrophy; by 80-100%) and commensalistic interaction (i.e.strengthening the interaction; by ~180%) at lower N:P ratio (i.e.Mediterranean Sea).This stimulatory effect was maintained or attenuated with increasing P limitation (i.e.high TN:TP ratios; Atlantic Ocean) under Current warming but became inhibitory (−10 and − 180%) under Current enr (Fig. 4A, C; Table S3; LSD post hoc test, p < 0.01).Under future conditions, the above-reported stimulatory effect under Current warming conditions turned strongly inhibitory (−40 and − 90%) when the P limitation was accentuated (i.e.Atlantic Ocean; Fig. 4B, D; Table S3; LSD post hoc test, p < 0.01).For Future enr conditions, we found a response pattern similar to that described under Current ones, although the magnitude of the impact was diminished on the C budget (only under lower P limitation) and on commensalistic interaction (Fig. 4B, D).Finally, on comparing the Future enr+warming scenario with Current conditions, we found that a stronger P limitation decreased the C budget of the microbial community (i.e.shifts towards heterotrophy, Figs. 4, S4) and weakened the commensalistic phytoplankton-bacteria interaction (i.e.BCD was not satisfied by EOC, Figs. 4, S4); however, such effects were lower than those found when the main drivers, i.e. nutrient enrichment or warming, acted separately.

Discussion
This work evidences that aeolian dust storms and MHWs are more frequent, more intense, and longer lasting in coastal Mediterranean and Atlantic Ocean waters, and that both phenomena are related to lower Chl a production.On the one hand, the negative effect of aeolian dust on phytoplankton agrees with previous experimental (Marañón  , 2010;González-Olalla et al., 2017) and modeling (Gallisai et al., 2014) results reported in Mediterranean and Atlantic ecosystems.Such effects have been attributed either to an inhibition of phytoplankton growth by trace metals also present in dust inputs (Paytan et al., 2009;Jordi et al., 2012), or to harmful effects of UVR on PP when radiation interacts with dust (or P) inputs (Carrillo et al., 2015b;González-Olalla et al., 2017).On the other hand, the reductions found in Chl a production in the Mediterranean Sea as MHWs duration increased are consistent with the weakened phytoplankton blooms (i.e.biomass production) reported by Hayashida et al. (2020) in Fig. 4. Mean (± pooled SD) ln response (LRR) of carbon budget (i.e.production:respiration ratio, PP:R) and microbial trophic interaction (i.e.bacterial carbon demand:excreted organic carbon BCD:EOC) in planktonic communities from an study area with low (Mediterranean Sea) and high (Atlantic Ocean) nutrients limitation exposed to current conditions of ultraviolet radiation (UVR), carbon dioxide concentration (CO 2 ), nutrients and temperature, Future conditions of UVR and CO 2 under a nutrient enrichment (Future enr ) or warming (Future warming ) alone, and under Future enr+warming conditions.Note that positive and negative values represent a stimulatory and inhibitory effect, respectively.nutrient-poor ecosystems subjected to MHWs.This negative effect of both drivers, i.e. nutrients (linked to AI) and temperature (linked to MHWs), over long-term scales contrast with a recent meta-analysis which have shown that neither experimental increases in P inputs nor warming significantly affected Chl a production in coastal regions, estuaries or the open ocean (Wu et al., 2021).Based on our long-term results concerning the effects of Saharan dust-deposition events and MHWs, and contrasting findings in the literature (nutrients and warming) over mid-term scales, we cannot disregard the possibility that additional environmental drivers (e.g.UVR, CO 2 ) and even biotic interactions, e.g. with potential competitors such as bacteria, may also be modulating the functioning of the ecosystem and communities.
Because nutrients and temperature are two major environmental drivers rapidly altered over the last four decades  in our study areas, two ecosystems with different nutrient limitation (i.e.Mediterranean Sea and Atlantic Ocean), and available evidences in the literature presented above show a highly variable effect of them: What are the nutrient×temperature effects under future globalchange conditions on microbial plankton communities subjected to different degrees of in situ P limitation?
The response to the above question is that the magnitude and direction of the cumulative impact of UVR, CO 2 , nutrients, and temperature on marine microbial plankton communities is strongly influenced by in situ P limitation.This result of nutrient limitation has already been stressed in laboratory experiments indicating more intense photosystem II photoinhibition by UVR (Litchman et al., 2002), decreased C fixation and maximal relative electron transport rates under high CO 2 (Zhang and Gao, 2021), or increased phytoplankton sensitivity to warming, undergoing reductions of between 3 and 6 °C in its optimum temperature for growth (Thomas et al., 2017).We provide evidence for the first time in a microbial plankton community, no effect (stimulatory or inhibitory) of Future enr+warming on PP with respect to current conditions in either of the two study areas.These results contrast with reductions found in PP (and growth) by Feng et al. (2020) when E. huxleyi populations were exposed to light, nutrients, CO 2 , and warming interactions.These authors reported that the reductions in PP and growth were due to a downregulation (to the extent of being virtually inhibited) of the RuBisCO and carbonic anhydrase activity expression.Although these activities were not explicitly measured in our study, we can rule out this possibility because we found that PP and Chl a (excepting Chl a in Atlantic Ocean) were stimulated in both areas after the acclimation to the individual and interactive effect of the drivers tested compared with the response measured before this period (i.e.Initial time; Fig. S3).
In contrast to PP, EOC underwent an opposite effect under Future enr +warming conditions in both marine study areas compared with current ones.That is, the former conditions boosted the release of EOC by two-fold in the Mediterranean Sea, and depressed it in the Atlantic Ocean (i.e. the area with greater P limitation).A higher EOC entails a higher total PP (i.e.PP + EOC), which, joined to the lower CR rates reported, supports the finding of autotrophy (i.e.PP:R > 1) in the Mediterranean Sea.Based on these results, we cannot discard the possibility that the higher EOC by phytoplankton in the Mediterranean Sea is a mechanism controlling bacteria.It is well known, both in freshwater and marine ecosystems, that bacteria exhibit a selective preference for this organic C source when available (Medina-Sánchez et al., 2002;Lagaria et al., 2011;Seymour et al., 2017), and hence it could have enabled phytoplankton to outcompete the bacteria in this nutrientlimited environment.A potentially higher competitive advantage of phytoplankton by limiting nutrients on bacteria in the Mediterranean Sea could support the reductions found both in BP and BR, and ultimately those on BCD, explaining the commensalistic coupling found (i.e.BCD:EOC < 100).This phytoplankton-bacteria coupling has been consistently reported in previous observational studies in other Mediterranean coastal waters (Fouilland et al., 2014), as well as in experimental (Carrillo et al., 2015a;Cabrerizo et al., 2019) and observational (González-Olalla et al., 2018) approaches in oligotrophic freshwater ecosystems.Thus, in agreement with the proposal by González-Benítez et al. (2019), the metabolic balance of the ecosystem (i.e.autotrophy vs. heterotrophy) appears to determine the source of C through which bacterial metabolism and communities are supported (i.e.EOC vs. another C source).Although it was not explicitly quantified in this study because it was not our aim, we speculate that a top-down control of bacteria by microzooplankton grazers (including mixotrophic species) could be a plausible explanation to the aforementioned reductions in BP and BR.In fact, in mesocosm experiments in the eastern Mediterranean Sea, Ptacnik et al. (2016) found a consistently strong negative relationship between heterotrophic bacteria net growth rate (i.e. a proxy of bacterial production) and bacterivorous phototrophic eukaryote abundance in microbial communities during the summertime.
For the Atlantic Ocean, we found reductions of 50% in EOC but similar CR rates under Future enr+warming (and Future warming ) that those measured under the current scenario.These similar CR rates denote persistent stressful conditions for phytoplankton communities in this ecosystem during summer.The predominance of a metabolic equilibrium (i.e.PP:R ~1; Future enr+warming or Future warming ) and net heterotrophy (PP:R < 1; Future enr ) in the Atlantic Ocean agrees with the results recently reported by Gazeau et al. (2021) for plankton communities from the ultraoligotrophic eastern Mediterranean Sea when subjected to dust×warming×CO 2 interactions; however, these authors suggested that this response pattern could be more dependent on the initial metabolic state in the microbial communities (i.e.auto-or heterotrophy) than on their own taxonomic composition (mainly nanoeukaryotes), or the global-change conditions tested or the in situ environment.The attenuated C-sink capacity evidenced in the Atlantic Ocean matched with an incipient decoupled commensalistic relationship (BCD:EOC > 100).This trend towards decoupling may be expected when organic C sources other than EOC become available for bacteria.Two additional explanations that are not mutually exclusive could be the following: (1) The changes in the composition (Viviani and Church, 2017) and lability (Cabrera-Brufau et al., 2021) of photosynthetically produced EOC.This varies with the growing phase (Myklestad, 2000) and taxonomic composition (Sarmento and Gasol, 2012) of phytoplankton.We did not quantify the molecular EOC composition over our experiments but we showed that two different groups Chlorophyceae (i.e.Atlantic Ocean) and Cyanophyceae (i.e.Mediterranean Sea) dominated in the studied areas during summer.Following this reasoning, we suggest that such communities could potentially release organic C with a different nutritional/energetic content and lability.(2) The effect of the grazing pressure also mentioned above.Previous results by Teira et al. (2015) with coastal microbial communities from the Atlantic Ocean suggested that a strong top-down (grazing) control could explain short-term disconnections between C supplies (EOC) and demands (BCD).Under this scenario, bacteria would essentially rely on DOC released as a consequence of predation.This decoupling phytoplanktonbacteria would imply reductions in C flux between the two trophic groups as the P-limitation accentuates (as occurred in Atlantic Ocean respect Mediterranean Sea), and thus a weakened microbial loop.
Our findings reinforce the importance that the nutrient background has in modulating the effects of multiple interacting globalchange drivers on microbial plankton communities, this providing diagnostic help to design future multi-driver experimental approaches.We also recognize that microcosm experiments cannot fully replicate natural environmental complexity (Hillebrand et al., 2012;Flombaum et al., 2020), and hence we do not envisage their outcomes to be completely representative of the impacts in natural systems.Nevertheless, it should be considered that a fuller understanding the effects of multiple global-change drivers on natural communities using systems of intermediate complexity is a mandatory step in building empirically grounded models that predict such impacts on microbial diversity dynamics.
In summary, these results of our multiple interacting driver manipulation constitute a starting point to understand how marine plankton communities will respond to complex future global climate.If we consider that the global ocean is already nutrient limited (Moore et al., 2013), the potential effects of UVR × CO 2 × temperature×nutrients on marine productivity and trophic interactions could be accentuated as the limitation exacerbates.Further studies going beyond acclimation processes, and considering ocean-basin scale approaches will determine whether microbial plankton adaptation will allow them to palliate (or not) the effects of rapidly changing ocean conditions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 2 .
Fig.2.Long-term response patterns of mean aerosol index (AI, relative units) intensity (A) and frequency (B), and marine heatwaves (MHW) frequency (C), duration (D), cumulative intensity (E) and severity (F) over 1985-2019 period in the Southwestern Mediterranean Sea and Northeastern Atlantic Ocean.Blue and red solid and lines represent the fit regression, and dotted thick lines 95% confidence intervals for Atlantic Ocean and Mediterranean Sea, respectively.

- 1 Fig. 3 .
Fig. 3. Mean chlorophyll a (Chl a) and aerosol index (AI, relative units; A) or marine heatwaves duration (MHW, B) relationship over 1985-2019 period in the Southwestern Mediterranean Sea and Northeastern Atlantic Ocean.Blue and red solid and lines represent the fit regression, and dotted thick lines 95% confidence intervals for Atlantic Ocean and Mediterranean Sea, respectively.