Surface (ca. 0.5 m) water was collected from three different sites along an inshore-offshore gradient in a coastal lagoon (Al Monseini), situated between King Abdullah University of Science and Technology (KAUST) and King Abdullah Economic City (KAEC), in the eastern coast of the central Red Sea (Figure 1). The three sites were dominated by different primary producers, two of them benthic macrophytes (one overlying Cymodocea serrulata seagrass meadows and the other in the vicinity of a patch of Avicennia marina mangroves) and the other one pelagic (phytoplankton). Hereinafter we will use the names Seagrass, Mangrove and Phytoplankton or their initials (S, M, and P) when referring to each site. The Seagrass (22°24.5′ N, 39°7.8′ E) and Mangrove (22°22.8′ N; 39°7.9′ E) sites were located inside the coastal lagoon, with depths of 1.5–2 m, and the Phytoplankton (22°19.1′ N, 39°3.2′ E, ca. 20 m) site was located offshore the lagoon entrance (Figure 1). Aiming at covering the seasonal cycle, sampling dates were staggered every 2–3 months. Actual samplings took place on 15 February 2016, 1 June 2016, 5 September 2016, and 4 December 2017. Temperature and salinity were measured immediately prior to sampling using a CTD (conductivity, temperature and depth) unit (Valeport MONITOR CTD + 500) and water sampling was performed by submerging two acid-washed polycarbonate carboys of 9 L per site.

For each site, experimental incubations of 4–6 days were performed following the seawater culture method (Ammerman et al., 1984). Two different treatments were used for monitoring bacterial and DOM dynamics at each site: unfiltered water was incubated in order to assess the dynamics of the entire microbial community (Community treatment) and water pre-filtered through pre-combusted Whatman GF/C filters (1.2 μm nominal size) was incubated to target the direct interaction between heterotrophic bacteria and DOM after removing protistan grazers and phytoplankton, virtually retained completely onto the filters (Filtered treatment). Incubations were performed in triplicates in 2 L acid-washed polycarbonate bottles in temperature-controlled incubators (Percival – I-22LLVL) with in situ light (12 h light:12 h dark cycle was set for all incubations, as in that Saudi Arabia location there are minor changes between light hours year-round) and in situ temperature conditions (which were adjusted for each site and season). The aim of the diel cycle of light and darkness in the Filtered treatment was to allow for photoheterotrophy of the natural bacterioplankton assemblages. Samples for inorganic nutrients, dissolved organic carbon (DOC) and nitrogen (DON) concentrations, DOM fluorescence (FDOM) components, bacterial abundance and bacterial physiological properties were collected twice per day to better capture the exponential phase of growth until the abundance of bacteria reached a plateau or started to drop (typically after 1 or 2 days). Incubations lasted up to 6 days, and were stopped when a stationary or decay phase was well-established. We focused on the initial responses and did not consider any secondary growth phases in our calculations, which were inconsistently detected.

Ambient size-fractionated chlorophyll a (Chl a) concentration was measured at each sampling site after sequential filtration of 200 ml samples through polycarbonate (Isopore™ membrane filters) filters of 20, 2, and 0.2 μm pore-size in order to respectively retrieve the micro-, nano-, and picoplankton size-fractions. Filters were kept at –80°C until analysis. Chl a was extracted in 90% acetone for 24 h in dark conditions at 4°C. A Turner Trilogy fluorometer was used to measure Chl a, using the acidification method (Holm-Hansen and Riemann, 1978) calibrated with a chlorophyll a standard (from Anacystis nidulans, Sigma Aldrich). Total Chl a was calculated by summing the three size-fractions.

Samples for nutrient analyses were collected in 15 mL sterile polypropilene Falcon tubes after being filtered through 0.2 μm Millipore^® polycarbonate filters and stored frozen at –20°C until analysis. A segmented flow analyzer from Seal Analytical^® was used to analyze nitrate (NO₃^–), nitrite (NO₂^–), and phosphate (PO₄^3–). Detection limits were 0.2, 0.06, and 0.01 μmol L^–1 for NO₃^–, NO₂^–, and PO₄^3–, respectively. Standards were prepared with a nutrient-free artificial seawater matrix in acid-washed glassware. Consumption and production rates of inorganic nutrients (μmol L^–1 d^–1) were estimated as the slope of the linear model of nutrient concentrations vs. time during the bacterioplankton exponential growth phases. Positive and negative values correspond to production and consumption, respectively.

Samples for DOC and total dissolved nitrogen (TDN) were filtered through 0.2 μm Millipore^® polycarbonate filters. To avoid contamination, filters were pre-washed with 50 mL of 1.2 mol L^–1 hydrochloric acid and rinsed with abundant MQ before sample filtration. The initial filtrate was discarded and the subsequent filtrate was used for sample storage and analysis. Each filtered sample was acidified with H₃PO₄ and kept at 4°C until further analysis by high temperature catalytic oxidation (HTCO) using a Shimadzu TOC-L. Reference materials of deep-sea carbon (42–45 μmol C L^–1 and 31–33 μmol N L^–1) and low carbon water (1–2 μmol C L^–1) were used to monitor the accuracy of DOC and TDN concentration measurements. DON concentrations were calculated after subtracting the dissolved inorganic nitrogen (DIN) to the TDN (DON = TDN – DIN), where DIN (μmol C L^–1) = [NO₃^–] + [NO₂^–].

Samples for DOM fluorescence analyses were filtered through 0.2 μm Millipore^® polycarbonate filters and measured immediately after filtration. UV-VIS fluorescence spectroscopy was measured using a HORIBA Jobin Yvon AquaLog spectrofluorometer with a 1 cm path length quartz cuvette. Three-dimensional fluorescence excitation emission matrices (EEMs) were recorded by scanning with an excitation range between 240 and 600 nm and an emission wavelength range of 250–600 nm, both at 3 nm increments and using an integration time of 8 s. To correct and calibrate the fluorescence spectra post-processing steps were followed according to Murphy et al. (2010), in which Raman-normalized Milli-Q blanks were subtracted to remove the Raman scattering signal (Stedmon et al., 2003). All fluorescence spectra were Raman area (RA) normalized by the subtraction of daily blanks performed using Ultra-Pure Milli-Q sealed water (Certified Reference, Starna Cells). Inner-filter correction (IFC) was also applied according to McKnight et al. (2001). MATLAB (version R2015b) was used for the RA normalization, blank subtraction, IFC and generation of EEMs. The EEMs obtained were subjected to PARAFAC modeling using drEEM Toolbox (Murphy et al., 2013). Before the analysis, Rayleigh scatter bands [first order at each wavelength pair where Ex = Em ± bandwidth; second order at each wavelength pair where Em = 2 Ex ± (2 × bandwidth)] were trimmed. Only one sample was identified as an outlier and the model was validated using half split validation and random initialization (Stedmon and Bro, 2008). The data array consisted of 568 samples that were split into two random halves of 284 EEMs each. The nested PARAFAC algorithm was then applied stepwise to both data arrays for 2–7 components. Ten iterations and a convergence criterion of 1e-6 were used for each model. The spectral properties of the components derived from each half were compared and found to be congruent according to the Tucker Coefficient described in Lorenzo-Seva and Ten Berge (2006) for a four-component model. The maximum fluorescence (Fmax) is reported in Raman units (RU) (Stedmon et al., 2003; Murphy et al., 2010). The four components from the validated model are: peak C1 at Ex₁(Ex₂)/Em 252(324)/397 nm (humic-like peak M, Coble, 2007), peak C2 at Ex₁(Ex₂)/Em 252(366)/467 nm (marine humic-like peak C, Coble, 2007), peak C3 at Ex/Em 303/337 [protein-like peak T (Coble, 2007) attributed to Tryptophan], and peak C4 at Ex/Em 270/312 nm [corresponds to protein-like peak B (Coble, 2007) and attributed to Tyrosine]. Consumption and production rates of each DOM fluorescent components (ΔRU/Δt, RU d^–1) were estimated as the slope of the fluorescence intensity vs. time during the exponential growth phases. Positive and negative values were considered as production and consumption, respectively. The percentage of protein-like compounds was calculated for each EEM as the sum of the fluorescence intensity of peaks C3 and C4 divided by the sum of fluorescence intensity of all peaks: (Fmax _C3 + Fmax _C4)/(Fmax _C1 + Fmax _C2 + Fmax _C3 + Fmax _C4) × 100.

We followed the methodology described in more detail in Gasol and Morán (2015) to characterize the physiological structure of heterotrophic bacteria at the three sites by considering five different single-cell groups (del Giorgio and Gasol, 2008): cells with low and high nucleic acid content (LNA and HNA, respectively), cells with healthy and damaged membranes (Live and Dead, respectively) and actively respiring cells (CTC+). LNA and HNA cells were differentiated according to relative nucleic acid content based on their green fluorescence signal after being stained with SybrGreen (Marie et al., 1997). Samples to distinguish these two groups were previously fixed with 1% paraformaldehyde + 0.5% glutaraldehyde (final concentration), deep-frozen in liquid nitrogen and stored at –80°C until analysis. Samples were thawed, stained and run in a flow cytometer within 1–3 months after collection. Membrane-intact (Live) and membrane-compromised (Dead) cells were assessed by combining two nucleic acid stains, SybrGreen (Molecular Probes) and propidium iodide (PI, Sigma Chemical Co), yielding green and red fluorescence signals, respectively. PI can only penetrate into cells with large holes in their membranes. Live and Dead cells were analyzed in vivo as described in Grégori et al. (2001). Actively respiring cells (CTC+), able to reduce the tetrazolium salt 5-cyano-2,3-ditolyl tetrazolium chloride or CTC, were distinguished by the red fluorescence signal resulting from the deposition of oxidized crystals of CTC (Sherr et al., 1999). CTC+ cells are Live cells representing typically a subset of the HNA group. This group was analyzed in vivo after 90 min incubation in the dark with the compound. All samples were analyzed in a BD FACSCanto II flow cytometer.

The abundances of the five physiological groups were calculated after gravimetric calibration of the flow cytometer flow rates. Total bacterial abundances correspond to the sum of LNA and HNA cells described in the previous section. Filtration by GF/C filters resulted in an overall mean loss (i.e., retention onto the filters) of 26 ± 10% (SD) of the heterotrophic cells present in the Community treatment. From a seasonal perspective, retention was minimum in September (14 ± 9%) and maximum in February (32 ± 5%). Latex fluorescent beads of 1 μm diameter (Molecular Probes, ref. F-13081) were added to each sample as an internal standard. Side scatter (SSC, light scatter at 90°) values relative to the fluorescent beads were converted into cell diameter following the empirical calibration of Calvo-Díaz and Morán (2006). Cell diameter was then converted to biovolume (Bv) and finally converted into bacterial biomass (BB) using the formula of Gundersen et al. (2002): fg C cell^–1 = 108.8 × [Bv]⁰.⁸⁹⁸ and then converted into μmol C L^–1.

The growth rates (μ) of four physiological groups (LNA, HNA, Live and CTC+ cells) were calculated as the slope of the natural logarithm of bacterial abundances vs. time during the corresponding exponential growth phases. Total heterotrophic bacterial growth rates correspond to the changes in abundance of total bacteria [LNA + HNA (reversed order for consistency)]. Viruses were unlikely to affect bacterial growth rates due to the short-time of the exponential growth periods (1–2 days, Guixa-Boixareu et al., 1999). Therefore, we will use the terms “specific growth rate” for the Filtered treatment and “net growth rate” for the Community treatment, containing the entire microbial food web (i.e., with protistan grazers). Carrying capacity is defined here as the maximum abundance recorded for each bacterial group at the plateau stage of the incubations. The ratio between the Filtered and the Community treatment carrying capacities was used as a proxy of top–down control.

Bacterial growth efficiencies (BGE, %) were estimated by following changes in DOC concentrations (ΔDOC/Δt) and bacterial biomass (ΔBB/Δt), during the exponential growth phase, as shown in the following formula: BGE (%) = [(ΔBB/Δt)/(ΔDOC/Δt)] × 100, where ΔBB/Δt represents the increase of total bacterial biomass (μmol C L^–1 d^–1), and ΔDOC/Δt is the carbon consumption (μmol C L^–1 d^–1) estimated from the decrease in total DOC.

Five liter of each site were filtered through 0.2 μm Sterivex filters for each sampling period (12 samples in total), frozen and extracted with the Power Soil DNAeasy kit (Qiagen). The regions V4-V5 of the 16S rRNA gene were amplified using the primers 515F-Y and 926R (Parada et al., 2016). Sequencing of 2 × 300 bp was performed on a MiSeq platform, where all steps performed and reagents used followed the MiSeq Illumina protocol. After trimming of indexes with cutadapat 2.3 (Martin, 2011), the clean sequences were analyzed with the DADA2 pipeline (Callahan et al., 2016, 2017) and the taxonomy was assigned from the reference database SILVA 132 at 80% bootstrapping. In total, 5950 final amplicon sequence variants (ASVs) were produced after removal of chimeras, singletons, and sequences assigned to mitochondria, chloroplasts or wrongly assigned to phyla different than Archaea or Bacteria. Samples were normalized to the minimum number of reads (21423) for all subsequent analyses.

JMP PRO 14 and R (version 3.6.1) software were used for statistical analyses. Repeated measurements analyses of variance (ANOVA), considering the month nested to station and post hoc comparisons between levels were performed using R “psycho” package (Makowski, 2018). Two-way ANOVAs were used to check for differences between bacterial physiological growth rates. Pearson r coefficients are given for correlation analyses. Sequence processing was performed with the dada2 package (Callahan et al., 2016) while diversity estimates and multidimensional scaling were calculated with vegan (Oksanen et al., 2019).

Warm temperatures were found year-round in the lagoon (Seagrass and Mangrove) and offshore water (Phytoplankton) sites, though showing a clear seasonality from 22.4°C in February to 33.3°C in September (Figures 2A–C). Salinity at the Phytoplankton site was relatively stable (39.3–39.5) while both the Seagrass and Mangrove sites showed higher variability (37.1–41.1 and 38.5–41.0, respectively). Total chlorophyll a concentrations were low to moderate (0.10–1.35 μg L^–1), with the maximum observed at the Seagrass site in June and the lowest values systematically detected at the Phytoplankton site (Table 1). The contribution of picophytoplankton (0.2–2 μm) ranged from 4 to 79%, while that of nanophytoplankton (2–20 μm) ranged from 10 to 62% (data not shown). Microphytoplankton (>20 μm) contribution to total chlorophyll a was the lowest of the three size-fractions in 58% of the samples (with a minimum value of 3.0%), but it made up to 74 and 87% inside the lagoon (Mangrove and Seagrass sites, respectively) in June.

The minimum and maximum inorganic nutrient concentrations at the onset of the incubations were similar for both nitrate and phosphate, with minima in June at Seagrass and maxima in December at Mangrove (Table 1). Contrary to inorganic nutrients, the Seagrass site showed much higher DOM concentrations than the other sites. DOC concentration averaged 121.6 ± 4.0 μmol C L^–1 there, with the Phytoplankton site consistently characterized by the lowest values (repeated measurements ANOVA, p < 0.001, n = 30). Overall, DOC concentration ranged between 77.2 and 137.3 μmol C L^–1 (Figures 2D–F) and it was positively correlated with temperature inside the lagoon (r = 0.89, p < 0.0001 at the Seagrass site; r = 0.74, p = 0.006 at the Mangrove site, n = 10). The concentration of DON (3.4–8.9 μmol N L^–1, Figures 2G–I) followed a similar pattern, with the Seagrass site showing the highest values and the Phytoplankton site the lowest. The corresponding C:N ratios ranged from 14 to 23 without any clear geographical difference (Table 1). Positive correlations between DON concentration and temperature were observed within each of the three sites (r = 0.68, p = 0.02 in Phytoplankton; r = 0.95, p < 0.0001 in Mangrove; r = 0.87, p < 0.001 in Seagrass, n = 10).

The fluorescent properties of DOM showed clear seasonal variations, with the percent contribution of protein-like substances higher in February and December and lower during the warmest months, especially at the Mangrove and Phytoplankton sites (Figures 2J–L). Total FDOM (i.e., the sum of the four fluorescent components detected) also showed higher values at the Seagrass site, followed by the Mangrove site (except in December, where the Mangrove site showed 1.7-fold higher fluorescence values due to a huge increase in the C3 peak, which dominated the fluorescence signal, Figure 3A). The lowest FDOM values were consistently observed at the Phytoplankton site regardless of the season (Figure 3A). In the warmest months sampled (June and September), the protein-like fluorescent components C3 (Tryptophan-like) and C4 (Tyrosine-like) were less abundant than the humic-like components (C1 and C2) at the 3 sites, accounting for less than 40% (Figures 2J–L, 3B). At the Seagrass site, a low contribution of protein-like compounds was observed throughout the year (repeated measurements ANOVA, p < 0.005, n = 30). A remarkable increase in the relative fluorescence of the C3 component was observed in December, reaching 62 and 74% at the Phytoplankton and Mangrove sites, respectively (Figure 3B). The humic-like fluorescent component C1, corresponding to peak A in Coble (2007), showed consistently higher fluorescence intensities than the other components at the Seagrass site. At the Mangrove and Phytoplankton sites, the fluorescent components showed higher variability (Figure 3B).

At the onset of the incubations, total heterotrophic bacterial abundance (i.e., the sum of the LNA and HNA groups) ranged from 2.3 to 7.7 × 10⁵ cells mL^–1 (Phytoplankton and Seagrass sites, respectively, Figures 4A–C). The Seagrass site showed significantly higher abundances than the other two sites (repeated measurements ANOVA, p < 0.001, n = 30), while the Phytoplankton site never exceeded 4.5 × 10⁵ cells mL^–1. HNA cells dominated the community at all sites most of the time, except at the Phytoplankton site in February and December, when %HNA dropped to 47%. A clear seasonal pattern in %HNA emerged at the three sites, with minima in February and December and maxima close to 80% in September (Figures 4D–F). %HNA increased significantly with temperature at the three sites (Supplementary Table 1). The Mangrove site showed a significantly higher contribution of HNA cells (repeated measurements ANOVA, p < 0.05, n = 30), exceeding on average by 4 and 10% the values of the Seagrass and Phytoplankton sites, respectively. Likewise total abundance, the Seagrass site also showed consistently higher abundance of Live cells (maximum of 5.14 × 10⁵ cells mL^–1), while the lowest abundances of membrane-intact cells were typically found at the Phytoplankton site (minimum of 1.91 × 10⁵ cells mL^–1). The contribution of Live cells (%Live) ranged from 84 to 98% (Figures 4D–F) and no significant differences were found between sites. %Live cells decreased significantly with increasing temperature only at the Seagrass site (Supplementary Table 1). The abundance of actively respiring cells (CTC+) also peaked at the Seagrass site (maximum 1.93 × 10⁵ cells mL^–1 in September) with one order of magnitude lower minimum values found at the Phytoplankton site in December (1.50 × 10⁴ cells mL^–1). The contribution of CTC+ cells followed the same geographical pattern, ranging between 6.6 and 28.0% (Figures 4D–F).

Total carrying capacity (i.e., the maximum abundance of LNA + HNA cells reached during the incubations) in the Filtered treatment without protistan grazers ranged one order of magnitude, from 3.23 to 10.8 × 10⁵ cells mL^–1 (Phytoplankton and Seagrass, respectively, data not shown but see Supplementary Table 2 for the individual physiological groups). Minima were observed in December at the three sites, while maximum values were reached in September within the lagoon and in June in offshore waters. Except in June, the Phytoplankton site showed consistently lower carrying capacities than the other two sites, but only significantly lower than the Mangrove site (repeated measurements ANOVA, p < 0.05, n = 36). The carrying capacity values were generally higher in the Filtered than in the Community treatment, although with some exceptions (Seagrass in June and December and Phytoplankton in September). The corresponding ratio of the Filtered to Community treatment carrying capacities (a surrogate for protistan grazers top–down control of heterotrophic bacterioplankton standing stocks) therefore ranged from 0.86 to 1.64, with the Phytoplankton site showing a different seasonal pattern to both lagoon sites (Figures 5A–C), which seemed to share a common seasonality in grazing pressure likely caused by their proximity (Figure 1B).

The growth rates of total heterotrophic bacteria (i.e., LNA + HNA cells) in the Filtered treatments (equivalent to specific growth rates) showed consistently higher values (paired t-test, p = 0.0013, n = 12) than the Community ones (equivalent to net growth rates), similarly to carrying capacities. In the Community treatment we observed mostly positive net growth rates (ranging from 0.06 to 0.85 d^–1, data not shown), except at the Seagrass site in February, where a decrease rather than an increase in abundance was recorded. In the Filtered treatment, specific growth rates ranged from 0.35 to 1.44 d^–1 and differed seasonally (Figures 5D–F), with the Mangrove site showing consistently and significantly higher values than the other sites (repeated measurements ANOVA, p < 0.01, n = 36). Although clear peaks were observed in June for the Mangrove and Phytoplankton sites, no significant correlations were found between specific growth rates and in situ temperatures either within individual sites or with all data pooled. Total specific growth rates showed, however, consistently negative relationships with the C:N ratio of DOM, which became significant at the Phytoplankton and Seagrass sites (r = –0.69, p = 0.0132 and r = –0.75, p = 0.0054, n = 12, respectively, Figure 6).

When it comes to single-cell physiological groups, the specific growth rates of HNA cells ranged from 0.34 (Seagrass in September) to 1.82 d^–1 (Mangrove in June, Supplementary Figure 1A), while for LNA cells minimum values were observed at the Seagrass site (0.15 d^–1) and maximum at the Phytoplankton site (1.05 d^–1), both in June (Supplementary Figure 1B). Similarly to HNA cells, Live bacteria also showed minimum and maximum specific growth rates at the Seagrass (0.33 d^–1) and Mangrove (1.33 d^–1) sites, respectively (Supplementary Figure 1C). The specific growth rates of both groups shared a similar seasonality (r = 0.76, p < 0.0001, n = 36). Lack of sufficient CTC+ specific growth rate measurements precluded us from apprehending any general pattern besides the fact than values were generally low, never reaching 1 d^–1 (Supplementary Figure 1D).

Since production and consumption of DOC was concurrently occurring in the Community treatment, drawing any sound conclusion about its dynamics was precluded. On the contrary, with phytoplankton and protistan grazers excluded, DOC dynamics generally showed a decrease with time in the Filtered treatment (Supplementary Figure 2), therefore allowing us to calculate DOC consumption rates. ΔDOC/Δt values ranged from 0.92 to 9.23 μmol C L^–1 d^–1 (Table 2). Hereinafter, only results from the Filtered treatment will be presented. The percentage of DOC consumed daily (%DOC) ranged from 0.8 to 9.0% (Seagrass and Mangrove site, respectively) with the Mangrove site showing consistently higher %DOC values than the other sites, but only significantly higher than the Seagrass site (repeated measurements ANOVA, p < 0.001, n = 27).

A small but variable part of the amount of DOC consumed by heterotrophic bacteria was used to build up their biomass. This increase in heterotrophic bacterial biomass was minimum in December at the Phytoplankton site (0.08 μmol C L^–1 d^–1) and maximum in June in the Mangrove site (0.44 μmol C L^–1 d^–1, Table 2), reflecting largely the changes in specific growth rates (Figures 5D–F). Bacterial growth efficiency (BGE), which represents the amount of consumed DOC used for bacterial BP, ranged from 3.1 to 18.0% (Table 2). We observed a significant correlation between the initial DOC concentrations and BGE at the Mangrove (r = 0.77, p = 0.0081, n = 9) and the Seagrass sites (r = 0.92, p = 0.0012, n = 8), but not at the Phytoplankton site. BGE was also significantly and negatively correlated with temperature at the Seagrass site (r = –0.77, p = 0.0253, n = 9) and with all data pooled (r = −0.60, p = 0.0011, n = 27).

The DOM fluorescent components C2 (humic-like) and C4 (protein [Tyrosine]-like) showed consistent consumption during the exponential phase of bacterial growth while C1 (humic-like) and less clearly C3 (protein [Tryptophan]-like) were mostly produced rather than taken up (Supplementary Figure 3). Overall, FDOM components ranged from a consumption rate of −0.0078 R.U. d^–1 to a production rate of 0.0112 R.U. d^–1 (C4 and C3, respectively, Supplementary Figure 4A). FDOM dynamics were clearer at the Seagrass site, likely because of its higher initial DOC and FDOM components concentrations regardless of the sampling period (Figures 2D–F and Supplementary Figure 3A). Interestingly, C4 consumption rates increased significantly with nanophytoplankton (chlorophyll a size-class) there (r = 0.82, p = 0.0012, n = 9) and at the Mangrove site (r = 0.81, p = 0.0013, n = 9), a relationship that held with all data pooled (r = 0.60, p = 0.0001, n = 26). Bacterial specific growth rate and BGE were both positively correlated to the production rate of C1 at the Seagrass site (r = 0.55, p = 0.06 and r = 0.77, p = 0.02, respectively, n = 9). BGE was also significantly and positively correlated to the consumption rate of C4 (r = 0.71, p = 0.05, n = 9) there. Consumption of inorganic nutrients was also observed in most incubations (Supplementary Figure 4B). The dynamics of nitrate ranged from a consumption rate of −0.14 μmol L^–1 d^–1 to a production rate of 0.32 μmol L^–1 d^–1 (Mangrove and Phytoplankton sites, respectively). At the Seagrass site no consumption of nitrate was observed. Changes in nitrite concentration were usually lower than for nitrate, except in February at the Seagrass site where a marked consumption was detected. Phosphate was consistently consumed at the Phytoplankton and Mangrove sites, peaking at −0.054 μmol L^–1 d^–1 in September (Mangrove). Phosphate consumption was only observed in September at the Seagrass site, with a maximum production of 0.075 μmol L^–1 d^–1 in December.

The analysis of the 16S rRNA gene confirmed that Archaea in surface coastal waters represented a very low fraction of the total abundance of prokaryotes (0.5% ± 0.2%). Overall, the most abundant phyla were Proteobacteria (49%), Cyanobacteria (24%, of which 94% were Synechococcus spp.) and Bacteroidetes (19%) (Figure 7A). Although their relative contribution differed between sites, differences at the phylum level were not statistically significant. The Seagrass and Mangrove sites showed fewer Cyanobacteria than the Phytoplankton site (ANOVA, marginally significant differences p = 0.055, n = 12) and different Proteobacterial groups, mainly in the Alphaproteobacteria class. The Seagrass and Mangrove sites had a higher percentage of the order Rhodospirilalles and Rhodobacteralles, while the Phytoplankton site had higher percentages of the SAR11 clade and Puniceispirillales (Supplementary Figure 5). Overall, the number of ASVs detected at each site was not significantly different (406 ± 33, n = 12), nor was the Shannon index (3.9 ± 0.1, n = 12). The distribution of samples in a Principal Coordinates Analysis showed the existence of 2 clusters along the PCoA 1 axis (Figure 7B), one gathering Mangrove and Seagrass samples (except June of the latter one) and another separating all Phytoplankton samples, representing the inshore-offshore gradient in the lagoon. Of all the environmental and biological variables considered, only chlorophyll a had a significant effect explaining the distribution of the samples.