Three 10-day cruises were conducted on board B/O Ramón Margalef in February (winter), April (spring), and August (summer) 2016 to cover a wide range of hydrographic and ecological conditions. Samples were taken at two locations in the eastern Atlantic Ocean, in the upwelling system near the Ría de Vigo: the Coastal station (88 m depth, 42.14° N, 8.95° W, sampled between 8:00 and 8:20 am) and the Offshore station (260 m depth, 42.10° N, 9.60° W, sample between 6:30 and 6:45pm) ( Figure 1A ). It is important to note that the different time of sampling may have biased the metatranscriptome comparisons between both sampling stations, as it has been previously shown that prokaryotic transcriptomes from natural samples may exhibit diel periodicity, particularly those of Cyanobacteria in open ocean waters (Ottesen et al., 2013; Ottesen et al., 2014).

Detailed information on the experimental setup and specific methods below are given in Supplementary Material .

Vertical profiles of temperature (°C), salinity, turbidity (NTU), total chlorophyll fluorescence, and photosynthetically active radiation (PAR) were obtained using a Seabird CTD rosette. Inorganic-nutrient determinations of collected water were done with a Bran + Luebbe segmented flow analyzer (Hansen and Grasshoff, 1983).

Chlorophyll a (Chl-a) concentration was measured by the non-acidification technique as a phytoplankton biomass proxy (Welschmeyer, 1994). Fluorescence was determined with a TD-700 Turner Designs fluorometer (absorption coefficient as 87.7 at 663 nm) (Lorenzen and Newton Downs, 1986). Prokaryotic abundance in seawater samples was calculated as described by Gasol and Del Giorgio (2000) using a Becton Dickinson FACSCalibur flow cytometer (488 nm light produced by argon‐ion laser).

Community composition was assessed from 24 samples by sequencing the 16S rRNA gene for prokaryotes and the 18S rRNA gene for eukaryotes (0.22 - 3µm size-fraction and > 3µm size-fraction, respectively). DNA was amplified using the universal primers 515F and 926R for prokaryotes (Parada et al., 2016) and TAReuk454FWD1 and TAReukREV3 for eukaryotes (Logares et al., 2014). Amplified regions were sequenced in an Illumina Miseq platform. The sequencing depth for amplicon 16S and 18S was 10,000 reads per sample The sequences obtained were analyzed with the software DADA2 for amplicon sequence variants (ASVs) (Callahan et al., 2016) using the SILVA reference database (v138.1) for taxonomic assignment of 16S (Quast et al., 2013), and the databases PR2 (Guillou et al., 2012) and the marine protist from the BioMarKs project (Massana et al., 2015) for taxonomic assignment of 18S ASVs.

RNA samples (< 3µm size-fraction) were taken during the experimental period at the surface in the coastal and oceanic stations (a total of 24 samples). Sampling for gene expression analyses was carried out on days 1, 3, 5, and 7 of each cruise. mRNA samples were extracted and sequenced using Illumina Miseq sequencing technology at the Science for Life Laboratory in Stockholm (SciLifeLab; www.scilifelab.se). The sequencing depth was 70,000 sequences. Illumina HiSeq 2500 raw paired-end reads (2x125bp) were quality checked with FastQC and MultiQC (Ewels et al., 2016), primer sequences were removed with cutadapt (Marcel, 2011), reads were trimmed with Sickle (Joshi and Fass, 2011) and filtered with ERNE (Del Fabbro et al., 2013) against a “contamination reference” database to filter reads designated as artifacts. Subsequently, forward and reverse reads were merged with PEAR (Zhang et al., 2014). All Merged reads were aligned with DIAMOND (Buchfink et al., 2014) against the NCBI RefSeq protein database (O’Leary et al., 2016). Afterward, all the sequences were compared to the 162 protein families (Pfams and TIGRfams) by running MEGAN (Huson et al., 2016) to annotate taxonomy and all possible transcripts to select marker genes relevant to a variety of ecological processes. The 162 genes found out of the 198 target protein profiles were grouped into nine categories and 44 subcategories ( Table S1 ). To make comparisons between samples, considering that the transcripts have different lengths, and the libraries differ in size, we used a Transcript Per Million (TPM) normalization.

Non-parametric Mann-Whitney U Test analyses were performed to compare the salinity, nutrient concentrations, bacterial abundance and chlorophyll a concentration during the different seasons at both stations. T-tests analyses were performed to compare the contribution of each taxonomical group to the 16S dataset and to the metatranscriptome. To test the null hypothesis that data came from a normally distributed population, a Kolmogorov-Smirnov & Lillierors test was used.

Analysis of Non-parametric MultiDimensional Scaling (nMDS) and PERmutational Multivariate ANalysis Of VAriance (PERMANOVA) were used to test the grouping of metatranscriptomes by seasons and stations. ReDundancy Analyses (RDA) were applied to show the variability in the metatranscriptomes that is explained by environmental data and by the composition of phytoplankton and bacteria communities and PERMANOVA was used to test the significance of explanatory variables.

A Benjamini-Hochberg correction was applied to the Spearman rank correlation analyses.

The results of the hydrographic surveys showed that downwelling conditions prevailed during winter and spring ( Figures 1B ; S1 – S2 ). During this period, the important vertical mixing and reduced light availability prevented phytoplankton blooms and relatively high concentrations of inorganic nutrients were measured in surface waters (NO₃ ^-, NO₂ ^-and SiO₂ concentrations and the N:P ratio during winter and spring at both stations were significantly higher than those measured in summer, Mann-Whitney U Test, p-value < 0.05; Table S2 ; Figure S3 ). At the coastal station, a persistent surface halocline associated with rain and terrestrial runoff was observed during winter and spring (during winter and spring surface salinity was significantly lower than during summer, Mann-Whitney U Test, p-value < 0.05) ( Figures S1 , S3 ; Table S2 ). By contrast, during summer, the upwelling brought deep, salty, and nutrient-rich water to the surface. In summer, a massive phytoplankton bloom was observed at the coast and spread to the offshore station: at the surface significantly higher Chl-a concentrations accompanied by significantly higher bacterial abundance were measured compared to spring and winter (Mann-Whitney U Test, p-value < 0.05) (see satellite images, Figure S4 ). The observed very low nutrients concentrations during summer in surface waters at the coastal station are probably due to the massive consumption of nutrients by the phytoplankton bloom ( Figures 1B ; S1 – S3 ; Table S2 ). At the coastal station, pronounced changes in Chl-a concentration were observed at the end of the winter cruise ( Figures 1B ; S1 , S5 ) and even more so, at the end of the summer cruise (likely due to the rapid wind relaxation that interrupted the upwelling and disrupted the bloom ( Figure S1 ), see Joglar et al., 2020 (Joglar et al., 2020). Collectively, there was more daily variation in physicochemical and biological oceanographic conditions (e.g. salinity, nutrients, chlorophyll a and bacterial biomass) in the coastal compared to the offshore and in summer compared to winter ( Figures S1 – S3 , S5 ; Table S2 ). An interesting exception to this was temperature, which registered a similar large daily variation in the offshore and coastal stations in summer due to the intense upwelling of cold water that precisely promoted the changes in the rest of variables.

The eukaryotic community at both stations was dominated by dinoflagellates and diatoms (Dinophyceae, Thalassiosira, and Chaetoceros, particularly in summer) ( Figure S6A ) (Joglar et al., 2020). Among the bacteria, Pelagibacterales and Synechococcales showed higher relative proportions at the offshore station, particularly in winter, while other Alphaproteobacteria and Flavobacteria reached higher relative abundances in spring and summer on the coast. At both stations, Oceanospirillales were relatively more abundant in winter, whereas Cellvibrionales increased in summer ( Figure S6B ) (Joglar et al., 2020).

Our metatranscriptomics analysis showed a dominance by Alphaproteobacteria (mainly Pelagibacterales and Rhodobacterales), Gammaproteobacteria (including Alteromonadales and Cellvibrionales), and Bacteroidota (mainly Flavobacteriales) in both stations and in the three sampled seasons ( Figure 2A ). The taxonomic distribution of transcripts for functional genes was similar in winter and spring, except for higher levels of Cyanobacteria at the offshore station in winter. During the summer phytoplankton bloom at the coastal station there was an important increase in the relative transcription levels of Cellvibrionales ( Figure 2A ). A relatively higher representation in the metatranscriptomics than in the DNA-based 16S rRNA gene amplicon data was observed for both Pelagibacterales (offshore station) and Other Gammaproteobacteria and Cellvibrionales (coastal station) during the three different seasons (t-test, p-value < 0.05) ( Figure S7 ).

Analysis of Non-parametric MultiDimensional Scaling (nMDS) and PERmutational Multivariate ANalysis Of VAriance (PERMANOVA) of bacterioplankton metatranscriptomes revealed a grouping of samples based on seasons and sampling stations (PERMANOVA, p-value < 0.05) ( Figure S8 ). Redundancy Analysis (RDA) showed that functional profiles of winter and spring bacterial communities from both stations were significantly (PERMANOVA, p-value < 0.001) associated with high concentrations of specific nutrients (NO₃ ^- and SiO₂) and with fairly high relative abundances (estimated from DNA-based 16S rRNA gene amplicon data presented in Figure S6 ) of the eukaryotic phytoplankton Ciliophora, Cryptophyceae, and Dictyophiceae (PERMANOVA, p-value < 0.005) and the bacteria Oceanospirillales, Pelagibacterales, and Alphaproteobacteria (PERMANOVA, p-value < 0.001) ( Figures 2B–D ). During the summer bloom, the bacterial functional profiles showed higher variability and appeared to be linked to increased bacterial abundance, Chl-a concentration, PO³, and NH₄ ⁺ concentration (PERMANOVA, p-value < 0.001) as well as with Chaetoceros, Dinophyceae, and Fungi at both stations (PERMANOVA, p-value < 0.005) ( Figures 2B–D ). At the offshore station, bacterial functioning during summer was associated with high relative abundances of Amylibacter, Flavobacteriales, and Rhodobacterales (PERMANOVA, p-value < 0.001), while at the coastal station, it was associated with higher numbers of Cellvibrionales (PERMANOVA, p<0.001) ( Figures 2B–D ).

We next assessed the microbial functions related to nutrient fluxes (C, N, P, S, and Fe) and to vitamin synthesis, response to stress and stimuli, and phototrophy in this upwelling system. Different seasonal ( Figures 3 ; S9 – S17 ), daily ( Figures 4 ; S18 – S23 ) and spatial patterns were observed in the different functions studied.

Since eukaryotic phytoplankton are key components of upwelling ecosystems – potentially shaping the conditions for bacterial activities – we carried out Spearman rank correlation analyses (applying a Benjamini-Hochberg correction) between the relative abundance of the major eukaryotic taxa in the 18S rRNA gene amplicon data and the relative transcription levels of the studied prokaryotic functional categories ( Figure 6A ) and environmental variables ( Figure 6B ). Notably, bacterial biomass, temperature, and salinity displayed significantly positive correlations (p-value < 0.05 marked with asterisks in Figure 6 ) with categories related to bacterial uptake and utilization of resources (e.g., metabolism of carbohydrates, amino acids, phosphate and sulfate, and ammonium uptake) while the opposite pattern was observed for nutrients concentration (p-value < 0.05 Figure 6B ). Significant positive correlations were also found between the categories related to bacterial uptake and utilization of resources and the relative abundances of Stramenopiles (particularly Chaetoceros), Dinophyceae, marine stramenopiles (MAST), and Fungi (p-value < 0.05) ( Figure 6A ). A positive relationship was found between the transcription of DMSP assimilation genes by the prokaryotic community and the abundance of Dinophyceae (p-value = 0.052) and marine alveolate (MALV) (p-value < 0.05) ( Figure 6A ). In contrast, other metabolic traits, such as sulfite assimilation, iron storage, and photosynthesis genes, showed a significantly positive correlation (p-value < 0.05) with Chlorophyta (particularly Ostreococcus and Micromonas) and Cryptophyceae ( Figure 6A ) and with low chlorophyll a values (p-value < 0.05) ( Figure 6B ).

Our analyses provided novel information about the potential mechanisms governing microbial interactions, nutrient utilization, and energy and matter fluxes in a representative upwelling system. We reason that it is essential to identify and characterize the large diversity of microbial processes and patterns for two important reasons. First, it partly results from two-way interactions with the inherently diverse distribution patterns of phytoplankton and bacteria within the microbial communities. Second, because of the complex interdependency between different physicochemical conditions and the myriad of ecological processes driven by different microbes in the ocean. Thus, focusing on the transcriptional investment by any one taxon in any particular metabolic pathway or a specific season of choice, would unavoidably have resulted in underestimating the magnitude of dynamics and distribution patterns characteristic of coastal upwelling systems. In the following we will discuss these findings with the intention to understand potential linkages between physicochemical and biotic processes.

The present work is, to our knowledge, the first studying high resolution spatial and temporal changes in bacterial functioning by means of metatranscriptomics in an upwelling system. The coastal area of the NW Iberian Peninsula is a highly dynamic ecosystem, affected by the intermittent upwelling of cold and inorganic nutrient-rich water that promotes phytoplankton and bacterial growth (Figueiras and Ríos, 1993; Nogueira et al., 1997; Álvarez-Salgado et al., 2002). Our dataset suggests that the coast to offshore gradient and the temporal changes in environmental and biological conditions promoted by the upwelling are key drivers of the observed pronounced taxonomic and functional bacterial diversification in this system. Several examples extracted from our results may be given.

The present dataset exemplifies some specific functional mechanisms behind the important role of nitrogen (N) and phosphorus (P) in structuring the functioning of microbial communities previously suggested to exist in this ecosystem (Martínez-García et al., 2010; Teira et al., 2016). Regarding nitrogen, the observed high expression of the amt, glnAB, and icd genes by Gammaproteobacteria (particularly Cellvibrionales) during summer at the coastal station was remarkable. The enzymes encoded by these genes are linked in the process of nitrogen assimilation from ammonium in that amt encodes an ammonium transporter, glnAB encode glutamine synthase that synthesizes glutamine from ammonium and glutamate (glutamine can next be used for amino acid synthesis), and icd encodes the TCA cycle enzyme isocitrate dihydrogenase that produces 2-oxoglutarate, which is the precursor of glutamate. This may suggest that the competitive ability of Cellvibrionales during the phytoplankton bloom in part lies in their ability to gear their metabolism towards utilizing N from ammonium as a means to utilize the C that they process. This is in line with the emerging view that the dialogue between the metabolic N and C cycles is largely mediated at the intercept between the TCA-cycle and the synthesis of amino acids (Forchhammer et al., 2022). At the offshore station, on the other hand, Pelagibacterales showed a relatively higher transcription of ammonium transporters compared to other taxa, suggesting this group is more capable of assimilating ammonium at the low concentrations typically found in oceanic areas. We think analysis of adjustments of central metabolism between biosynthetic as compared to energy-yielding pathways (such as respiration) has a large potential to inform on how marine bacteria meet challenges in elemental stoichiometry of their resources (as imposed through e.g. C- or N-limitation). Especially if advances are made in characterizing gene systems like amino acid transporters and extracellular proteases involved in the utilization of dissolved organic nitrogen.

The high expression of the gene that encodes the transport system for orthophosphate (pstB) during summer, when phosphate concentration was low, is coherent with the suggested role of the Pst system (pstABCS) in inorganic phosphate (Pi) uptake only under Pi-limiting conditions (Harke and Gobler, 2013; Jimenez-Infante et al., 2017). On the other hand, the transcription of the pit gene, a low-affinity phosphate membrane transporter typically expressed under high phosphorus concentrations and post-bloom conditions (Alonso-Sáez et al., 2020), peaked at the coast in spring, coinciding with the highest ammonium and phosphate concentrations registered. Note here that, even though the phosphate concentration was similar in winter and spring at the coastal site, dissolved inorganic N, particularly ammonium, was considerably higher in spring, suggesting that higher expression of the pit gene in spring than in winter might be related to the higher ammonium to phosphate ratio. Overall, the different spatial-temporal patterns of variability observed in genes with different phosphate affinity in this upwelling system, regardless of the relatively limited variability in phosphate ambient concentration, reinforces the hypothesis of a complex role of phosphate in structuring bacterial functioning (Harke and Gobler, 2013; Satinsky et al., 2014).

Our results support the idea of a taxonomic diversification in iron (Fe) provision strategies in this upwelling system. Thus, while Cellvibrionales and Flavobacteriales are suggested to mostly use siderophore-associated Fe, Pelagibacterales probably use ferric ions and Synechococcales accumulate Fe as ferritin in this ecosystem (Ahlgren et al., 2020). On the other hand, a widespread genetic capacity for the uptake and utilization of sulfate in surface waters off NW Iberian Peninsula can be deduced from our data, as opposed to previous knowledge limiting sulfate use during phytoplankton blooms to a few bacterial groups (Zhou et al., 2020). Interestingly, our results also highlight a probable dominant role of Pelagibacterales in sulfite assimilation in the system, extending to upwelling areas previous findings in offshore epipelagic layers (Needham et al., 2017). The present dataset suggests for the first time the relevant role of organic sulfur (S) compounds as sources of reduced S for marine prokaryotes in this upwelling system, confirming previous studies in other areas (Kiene et al., 2000; Moran et al., 2003). Interestingly, the high relative abundance of dmdA transcripts at the offshore station during the summer upwelling period coincided with a dominance of the phytoplankton community by Dinophyceae. Dinoflagellates are known to have high dimethylsulfoniopropionate (DMSP) production rates (Zhao et al., 2021), and previous work has suggested that their DMSP-producing ability could affect the structure of their associated bacterial community (Lin et al., 2021). In this regard, the observed dominant role of Pelagibacterales in the transcription of DMSP utilization genes expands the work by Varaljay et al. (2012) in the North Pacific subtropical gyre to this upwelling system. Furthermore, DMSP has been previously reported to have a photoprotective function (Archer et al., 2009), which may be related to the high radiation levels and the strong thermal stratification registered in the offshore station during summer.

Finally, the present dataset suggests that different strategies for using light as a source of energy by bacteria occur in this upwelling system. Thus, the relative transcription of photosynthesis related genes by Synechococcales was more important in the offshore station during winter, where the abundance of this group was relatively higher (Joglar et al., 2020). It is important to note here that, as explained in the Methods section, the coastal and offshore stations were sampled at different time points in the daily cycle. Thus, the high relative transcription of photosynthesis-related genes by Synechococcales during winter in the offshore station may be not only related to the broad differences in physicochemical and biological conditions (higher relative abundance of Synechococcales during winter in the offshore station, see Figure S6 ) but also partly to an increase in the investment in light capture due to the sampling time. In this regard, previous studies have shown that prokaryotic transcriptomes from natural samples may exhibit diel periodicity, particularly those of Cyanobacteria in open ocean waters (Ottesen et al., 2013; Ottesen et al., 2014). Interestingly, the data on the transcription of the prd gene suggest that photrophy mediated though proteorhodopsin was relevant and widespread among different taxa, including Pelagibacterales, in accordance with the results by Olson et al. (2018) in oligotrophic waters and Gammaproteobacteria during the coastal phytoplankton bloom.

Previous studies have shown recurring daily patterns in the functioning of heterotrophic bacteria from relatively stable oligotrophic environments (Ottesen et al., 2014; Muratore et al., 2022). It is also well known that the metabolism of naturally occurring bacteria experience inter- and intra-seasonal changes associated to phytoplankton blooms (Teeling et al., 2012; Teeling et al., 2016). However, this is to our knowledge the first time that systematic changes in specific microbial functions (identified from gene expression analyses) associated to short-term oceanographic variability are shown. Our results evidence important short-term dynamics in bacterioplankton functioning superimposed on seasonal variability that coincided with changes in hydrography, nutrient availability and biotic variables in this system. Interestingly, this day-to-day changes in transcripts for specific functions related to nutrient fluxes and response to stimuli were even more pronounced when daily variability in oceanographic conditions (e.g. temperature, phytoplankton biomass, resource availability) was highest, that is: close to the coast during summer.

In summer, a wind relaxation event interrupted the ongoing upwelling, which modified the vertical density gradient at the coastal station, altered the hydrographic conditions and nutrient availability, and disrupted the phytoplankton bloom (see (Barbosa et al., 2001) and (Joglar et al., 2021b)). Concomitantly, a sharp decrease in the transcription of functions related with bacterial utilization of resources derived from the phytoplankton bloom was observed, including amino acid metabolism (glnAB), TCA cycle (icd), and ammonium uptake (amt). This suggests that bacteria quickly responded to the changes in resource supply and potentially also in the C/N ratio of resources, as previously suggested by others (see review by Forchhammer et al. (2022)). Also, the relative importance of high-affinity phosphate uptake genes (e.g., pstS and pstB, which are considered an indicator of microbial P limitation and that varies seasonally in relation to P availability (Dyhrman et al., 2007)) decreased quickly after bloom disruption (on day 5), this probably due to the cease of the competition with phytoplankton. In contrast, the transcription of low-affinity phosphate uptake gene pit followed the opposite pattern and relatively increased at the end of the summer sampling at the coast when phosphorus availability increased. Similarly, bacteria shifted from a strategy based on iron storage to an increase in the relative importance of iron uptake (as suggested by the high relative transcription of ftnA and fbpAB, respectively) when the phytoplankton bloom was interrupted. Interestingly, the important short-term changes observed in the transcription of DMSP assimilation genes coincided with day-to-day peaks in the relative abundance of dinoflagellates during the final phase of the summer cruise at both stations. Overall, these results suggested a close coupling between the physicochemical conditions and the bacterial and algal component of microbial communities from this upwelling system.

Overall, our results demonstrate that specific functional processes of the prokaryotic community identified from gene expression analyses) are associated with different phytoplankton groups in the upwelling system off the NW Iberian Peninsula. The presence of different phytoplankton and bacterioplankton populations during periods with different environmental conditions has been suggested in previous studies both in the area (Joglar et al., 2020; Costas-Selas et al., 2023) and in other marine ecosystems (Pinhassi et al., 2004; Needham et al., 2017; Bunse et al., 2019). Recent work in the NW Iberian Peninsula suggests a key role of biotic interactions as structuring factors of the eukaryotic community, mostly driven by positive associations between phytoplankton and bacteria (Costas-Selas et al., 2023). These findings may suggest some coevolution between the algae and the associated bacterial community that extends beyond the ability to transform phytoplankton-derived macromolecules by bacteria. The present dataset unravel the specific metabolic processes and bacterial and phytoplankton taxa that may be behind these associations through space and time. It is important to note here again that the interactions between phytoplankton and bacteria may change throughout the light cycle as it is known to affect the functioning of photosynthetic organisms (Ottesen et al., 2013; Ottesen et al., 2014). Therefore, the sampling at different times in both stations might influence to some extent the results discussed below.

Two interesting examples of possible phytoplankton-bacteria linkages observed during the productive season were associated to the bacterial transcription of genes related to the metabolism of DMSP (see above) and to the synthesis of B-vitamins. The observed increase in the relative abundance of transcripts of vitamin-related genes during summer is likely related to the high phytoplankton and bacterial requirements characteristic of blooming populations (Croft et al., 2005; Koch, 2012; Sañudo-Wilhelmy et al., 2012; Koch and Trimborn, 2019). Interestingly, genes belonging to the aerobic pathway of B12 synthesis (e.g. cobG and cobT), which is mostly associated with phytoplankton blooms (Zhou et al., 2020), increased during summer, while transcription of cbiX, which belongs to the anaerobic pathway of B12 synthesis that has been previously suggested to be independent of phytoplankton activity (Grossman, 2016), did not follow a clear seasonal pattern. The relative increase in thiL and cobG transcription by Cellvibrionales and Flavobacteria during summer suggests that these groups provided B1 and B12 vitamins to blooming eukaryotic phytoplankton. Interestingly, this hypothesis is also supported by microcosms experiments stimulating phytoplankton growth and including B1 and B12 vitamins additions carried out in summer in the NW Iberian Peninsula upwelling system that showed relatively important contributions of Cellvibrionales and Flavobacteria to the transcription of genes related to vitamin metabolism (Joglar et al., 2021a).

Finally, our results suggest strong functional linkages between bacteria and phytoplankton during non-blooming winter conditions. In this regard, the association between transcripts related to bacterial metabolic traits such as sulfite assimilation, iron uptake and storage, and photosynthesis, and Clorophyta and Cryptophyceae may reflect the close association between these phytoplankton groups and bacteria from the orders Pelagibacterales and Synechococcales, which played prominent roles in the relative transcription of genes associated with these metabolic functions. The close connection between these phytoplankton groups and Synechococcus sp., also described by Costas-Selas et al. (2023) in the same study area, were hypothesized to be linked with B12 provision. Our results point to an additional interaction mechanism between these groups, which could be related to the high capacity of Synechococcus to store Fe. As Clorophyta and Cryptophyceae may have mixotrophic behavior (Stoecker et al., 2017), predation on Synechococcus might be a source of Fe for them. A previous experimental study observed a sharp decrease in genes related to Fe uptake in a mixotrophic Haptophyta fed with bacteria, suggesting that bacterivory may provide Fe to the algae (Liu et al., 2015).

The close correlation between Chlorophyta (particularly Ostreococcus and Micromonas), Cryptophyceae, and Pelagibacterales suggests a link mediated by the sulfur metabolites, given the strong correlation between these phytoplankton groups and the relative abundance of sulfite assimilation genes. In this regard, a recent study suggests that adenylylsulfate reductase genes (aprAB) in Pelagibacterales clades could be related to taurine metabolism, oxidizing the sulfite derived from the conversion of taurine into acetyl-CoA (De Corte et al., 2021; Ruiz-Perez et al., 2021), a result consistent with the demonstrated ability of Ostreococcus and Micromonas to produce high quantities of taurine (Durham et al., 2019).

Knowledge of the ecophysiology and ecology of Cellvibrionales is still scarce, but they appear to reach highest relative abundances in surface waters and coastal areas (Chénard et al., 2019; Pajares et al., 2020; Reji et al., 2020). Although Cellvibrionales is not the most abundant bacterial group in this system, it appears to play a relevant and unexpected role during the summer phytoplankton bloom formation and development.

First, our field data suggest a prominent role of Cellvibrionales in degradation of phytoplankton-derived organic matter. Thus, during the Dinophyceae and Chaetoceros summer bloom, Flavobacteriales and Cellvibrionales dominated the relative expression of genes related to resource utilization and remineralization. It is well known that the primary role of Flavobacteriales in the degradation of polymeric organic matter (Teeling et al., 2012; Fernández-Gómez et al., 2013) is especially important during diatom blooms in upwelling ecosystems (Klindworth et al., 2014). However, the high transcription by Cellvibrionales of carbohydrate and TCA cycle-related genes as well as ammonium utilization and amino acid synthesis genes during the phytoplankton bloom presented here is remarkable. This behavior may be related to efficient energy production (Wang et al., 2022) and to adjustments to C or N-limitation by this group when resources are available during summer phytoplankton blooms (see above). Thus, the increase in the transcription of C metabolism-related genes by Cellvibrionales is consistent with recent findings in concurrent mesocosms experiments (performed in the same sampling area during the present study) (Pontiller et al., 2022) and in other regions (Poretsky et al., 2010; Liu et al., 2020), suggesting a quick response of Cellvibrionales to phytoplankton-derived organic matter. In surface waters of coastal systems like the NW Iberian Peninsula, the Gullmar Fjord and the German Bight, Cellvibrionales have been shown to importantly invest energy in the expression of different polymer degrading enzyme systems including for example laminarases that allow bacteria to decompose the algal storage glucan laminarin (Teeling et al., 2012; Pontiller et al., 2021; Pontiller et al., 2022) emphasizing their significant role in the turnover of phytoplankton-derived organic matter in these systems. Similarly, the relatively abundant Cellvibrionales transcription of motility-related genes during the summer bloom is in accordance with their description as a copiotroph group, known to quickly respond to the increase of DOM surrounding the microenvironment of phytoplankton cells in productive environments (Smriga et al., 2016). In addition, the high transcription of the proteorhodopsin gene (i.e. prd) by Cellvibrionales during summer may indicate that this group uses proton pumps as a supplemental energy source to uptake recently released complex carbon compounds during phytoplankton blooms (Pinhassi et al., 2016).

Our data suggest that the association between blooming phytoplankton and Cellvibrionales may be more complex than just the ability to transform phytoplankton derived macromolecules. We found higher levels of expression of genes for vitamins B1 and B12 synthesis by these bacteria associated with phytoplankton groups dominant during the blooms (i.e., Chaetoceros and Dinophyceae). This finding suggests that Cellvibrionales have an important role as suppliers of micronutrients needed for phytoplankton growth and bloom development. Similarly, the increase in the expression of siderophore-uptake related genes by Cellvibrionales associated with bloom-forming phytoplankton may indicate that these bacteria compete with Chaetoceros and Dinophyceae algae for organically-bound iron during the bloom season (Amin et al., 2009). Therefore, our data suggest that Cellvibrionalles develop a mutualistic relationship with blooming phytoplankton by producing growth-promoting compounds like vitamin B12 and facilitating phytoplankton functioning. They also show, as explained above, that Cellvibrionalles seem to act as competitors of blooming phytoplankton for iron provision in this productive upwelling system.

In conclusion, the present work highlights the mosaic of microbial processes governing energy and matter transformations in this representative upwelling system. Uncovering the linkages between bacterial activity and spatial-temporal variability in environmental conditions represents a relevant goal for microbial oceanography The systematic patterns that emerge from this work are encouraging in the context of future modeling attempts aiming to predict biogeochemical processes in such highly productive ecosystems.