Comparative metagenomics of phytoplankton blooms after nutrient enrichment of oligotrophic marine waters

Increasing anthropogenic pressures on the coastal marine environments impact these ecosystems via a variety of mechanisms including nutrient loading, leading to eutrophication and increases in algal blooms. Here, we use a metagenomics approach to assess the taxonomic and functional changes of the microbial community throughout a nutrient enriched mesocosm phytoplankton bloom. We tested four different nutrient treatments consisting of either nitrate and phosphate or nitrate, phosphate and silicate, administered on the first day or continuously for the first two weeks of the experiment. Our results show a shift in the taxonomic composition of the community over time that is dependent on the nutrient addition regime. Significant differences in the functional potential of the communities were detected, with an interaction between bloom period (pre-bloom, bloom and post-bloom) and nutrient treatment (p = 0.004). A sharp drop in functional similarity was observed in the first week in all treatments and after 20 days had not returned to pre-bloom levels. Changes within energy metabolism pathways showed a remarkable enrichment of the dissimilatory nitrate reduction pathway in the post-bloom period. Eukaryotic oxidative phosphorylation and photosynthetic antenna proteins were more abundant during the bloom, especially in the continuous treatment with silicate. Our results suggest that continuous (i.e. chronic) nutrient enrichment has a larger effect on the functioning of marine systems compared to a single (i.e acute) addition. A deep under standing of the functional and taxonomic shifts in the community during blooms is essential to reverse or mitigate human impacts on coastal environments.

Despite the importance of microbial organisms to the functioning of the marine system, the inherent difficulty to understand what controls a complex system with multiple connections (Huisman and Weissing 1999) has limited our understanding of the mechanisms organising microbial plankton communities. With the advent of high throughput sequencing technologies, there has been an increased understanding of the marine microbial community. Studies have primarily looked at changes in the taxonomy of the organisms contributing to the community and relating them to changes in the environment (Fuller et al. 2006;Gilbert et al. 2012;Massana et al. 2015, e.g. de Vargas et al. 2015Pearman et al. 2016bPearman et al. , 2017. However, more recently, investigations involving deep sequencing the genes of these marine communities (functional potential) have increased with studies revealing a vast wealth of novel genes (Yooseph et al. 2007;Sunagawa et al. 2015;Duarte et al. 2020). Global expeditions, such as the Global Ocean Survey (Yooseph et al. 2007), Tara Oceans  and Malaspina 2010 Circumnavigation (Acinas et al. 2021), have provided insight into spatial patterns in the metabolism of marine microbes on a worldwide scale, while seasonal patterns have been shown during temporal studies (Biller et al. 2018;Yoshitake et al. 2021).
Human population growth has resulted in the increased development of coastal regions leading to greater anthropogenic pressures on the coastal marine environments. One of the consequences of the increased development of the coastal zone is the higher discharge of anthropogenically related material into the marine environment. This can have an impact on the ecosystem via a variety of mechanisms including increased nutrient loading (Malone and Newton 2020) and consequent eutrophic processes. For instance, higher nutrient loading can cause increased rates of primary production, changes in algal biomass and composition, reductions in transparency and loss of biodiversity (Smith 2006). Transient increases in the abundance of phytoplank-ton are referred to as blooms and are usually dominated by one or a few species (Irigoien et al. 2004). These blooms are increasing in frequency, with elevated nutrient availability being one of the principal drivers behind this pattern (Lancelot et al. 1987;Paerl 1988;Huppert et al. 2002;Beman et al. 2005). An understanding of how the functional potential of a microbial community changes during nutrient enrichment and the resulting blooms has received less attention, although differences have been shown between bloom microcosms and controls (Rinta-Kanto et al. 2012) and with deep water nutrient enrichment experiments (Shi et al. 2012). Further, metatranscriptomic approaches have also been used to highlight changes in the metabolic pathways related to dinoflagellate blooms (Gong et al. 2017). Pearman et al. (2016aPearman et al. ( , 2016bPearman et al. ( , 2017, using a metabarcoding approach, showed that nutrient additions led to changes in the composition of the microbial community within the Red Sea microbial plankton communities. We undertook shotgun metagenomic sequencing at various time points throughout a nutrient-enriched mesocosm phytoplankton bloom to investigate the taxonomic and functional changes of the microbial community ( Fig. 1). We hypothesised that: i) the metagenomic sequencing will show changes in the taxonomic profile of the microbial community amongst treatments and over time; ii) compared with the control, the treatments will have a different functional potential; and that iii) the functional potential will change during the bloom with the post-bloom functional potential returning to a similar state to that of the pre-bloom community.

Experimental design
Mesocosm bags of 8000 l (depth: 2.5 m) were placed in the harbour of King Abdullah University of Science and   Table S1). Mesocosm bags (8000 l) were subjected to five different treatments (2 replicates per treatment): 1) a single addition of nitrate (NaNO 3 ; 16 µM) and phosphate (H 2 NaPO 4 .H 2 O; 1 µM), referred to as the NP treatment; 2) a single addition of nitrate (16 µM) and phosphate (1 µM) and silicate (Na 2 SiO 3 .9H 2 O; 39 µM), referred to as the NPS treatment; 3) a continuous addition of nitrate (2 µM each day) and phosphate (0.12 µM each day) for the first two weeks, referred to as the NP.cont treatment; 4) a continuous addition of nitrate (2 µM each day), phosphate (0.12 µM each day) and silicate (3.75 µM each day) for the first two weeks, referred to as the NPS.cont treatment; and 5) control without nutrient additions. The ratios were adapted from those published by Wyman et al. (2000) and confirmed by a small scale microcosm experiment to show that the concentrations were able to induce a phytoplankton bloom (according to microscope counts).

Sample collection
A CTD profiler (Valeport Monitor CTD Profiler with an attached chlorophyll sensor) was used to measure daily profiles for temperature, salinity and fluorescence. The CTD was calibrated as per Millard and Yang (1993).
Samples for metagenomic analysis were collected in 20 l Niskin bottles and directly transferred to carboys for immediate filtration. To avoid destruction of delicate cells during filtration, peristaltic pumps at low speed (70 rpm) filtered approximately 4 l of seawater through a 0.2 µm hollow fibre CellTrap (Mem-Teq, UK). The concentrated cells within the CellTrap were eluted using 2 ml of filtered seawater (from the same sample) and the eluted cells were immediately frozen in liquid nitrogen and stored at -80 °C for later analysis.
Nutrient samples were collected as per the protocols of Brügmann and Kremling (1999) and Kremling and Brügmann (1999). In brief, water samples were filtered through 0.45 µm filters into acid washed (10% HCl acid overnight) sample bottles. Nitrate, nitrite and orthophosphate concentrations were assessed by means of Autoanalyser AAIII pentacanal BRAN-LUEBBE (with 10 nm optical path) using the protocols of Cauwet (1999) at AZTI Tecnalia. Ammonia was analysed by a SEAL Analytical AutoAnalyzer 3 at KAUST using the manufacturer's protocols.
Counts of total bacteria, the photosynthetic bacterial genus Synechococcus, photosynthetic picoeukaryotes and photosynthetic nanoeukaryotes were determined using a FACSVerse flow cytometer (Becton Dickinson, Belgium). Samples filtered through a 0.4 µm mesh were fixed in glutaraldehyde (2.5% final conc.) and flash frozen in liquid nitrogen. Cells were excited with a blue laser (488 nm wavelength) and red fluorescence emissions measured to discriminate the eukaryotic phytoplankton. Beads of size 1.002 µm (Polysciences, Europe) were added for size verification.

DNA extraction and metagenomic sequencing
DNA from the concentrated cells was extracted using a phenol:chloroform:isopropanol methodology combined with bead-beating, as described in Pearman et al. (2016a). Based on results described in Pearman et al. (2016a), samples targeting specific aspects of the bloom were selected for metagenomic sequencing. The pre-bloom sample was from day 1, while the peak bloom samples were taken on days 5, 6 and 8. Post-bloom samples were taken on days 13, 16 and 20 (note that samples for the control on day 16 are missing as they did not sequence). Paired end sequencing libraries (100 × 2 bp) were prepared following the manufacturer's protocols of the NEBNext Ultra DNA kit (#E7370L). After library construction, six samples per lane were sequenced on an Illumina HighSeq 2000 sequencer at the KAUST Biosciences Corelab (BCL). Raw sequences were stored in the National Centre for Biotechnology Information (NCBI) short read archive under the accession number: PRJNA395437.

Bioinformatics
PhiX sequences were removed from the samples using the bbduk.sh script (http://jgi.doe.gov/data-and-tools/ bb-tools/) before processing using the software Squeeze-Meta v.1.3.0 (Tamames and Puente-Sánchez 2019) with the co-assembly option. Firstly, sequences were trimmed with Trimmomatic v.0.39 (Bolger et al. 2014) to remove leading and trailing bases below a quality score of 20 and reads with an average per base quality of less than 30 over a 4 bp window. Assembly was undertaken using MEGAHIT (Li et al. 2015) with short contigs (< 500 bp) removed with prinseq (Schmieder and Edwards 2011). Reads for each sample were mapped to the assembly using Bowtie2 (Langmead and Salzberg 2012). Open reading frames (ORFs) were predicted using Prodigal (Hyatt et al. 2010) and similarity searches undertaken against the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto 2000) databases using DIAMOND with the options of evalue = 0.001 and id = 30 (Buchfink et al. 2015). For comparison across samples, the ORFs were normalised to reads per million (RPM). The RPM is calculated as reads per kilobase (RPK) / scaling factor, where the scaling factor is the sum of all RPK values divided by 10^6 reads. DIAMOND was used for taxonomic classification of the ORFs against the NCBI nr database with the SqueezeMeta default options evalue = 0.001 and id = 40.
To assess the taxonomic communities, open reading frames that had a phylogenetic classification to at least the kingdom level was used. The taxonomic composition was visualised in two ways. Firstly, to give an overview of the taxonomy in each treatment, ORFs were amalgamated at the class level and the percentage abundance calculated for each sample. The mean percentage abundance for each class was then calculated for each treatment. Secondly, after calculating the percentage abundance for each class, the mean was calculated for each day in each treatment. The results were visualised in in ggplot2 v.3.3.5 (Wickham 2016).
The functional potential of the metagenomes was assessed at the level of KEGG Orthology ID (KO ID) with ORFs being amalgamated at this level. To assess the shared and unique KO IDs amongst treatments, all replicates within a treatment were merged together and then the presence/absences of KO IDs in the different treatments was assessed and visualised in the R package Ven-nDiagram v.1.6.0 (Chen and Boutros 2011).
To test multivariate differences on the relative abundance (Bray Curtis) of KO IDs, a PERMANOVA was undertaken in the R package vegan v.2.7 (Oksanen et al. 2019) using the factors bloom period (3 levels; prebloom, bloom and post-bloom) and treatment (5 levels; Control, NP, NPS, NP.cont, NPS.cont) with an interaction term. To assess the association of environmental variables with the functional potential, a canonical correspondence analysis (CCA) was undertaken using the variables: ammonia, nitrate, phosphate, silicate concentrations and the counts of total bacteria, the photosynthetic prokaryote genus Synechococcus, picoeukaryotes and nano eukaryotes.
Raw read numbers were used to calculate differentially-abundant KO IDs using the DESeq2 package v.1.32.0 (Love et al. 2014) in R. KO IDs were considered differentially abundant if they had an adjusted (FDR) p-value < 0.05 in pairwise comparisons between the control and a treatment. Enriched metabolic pathways were identified by filtering the differentially-abundant KO IDs to only KEGG pathways involved in metabolism. These differentially-abundant KO IDs were used as the input for enrichment pathway analysis using clusterProfilier (Wu et al. 2021) with a pvalue cutoff of 0.05 and a minimum gene set size of 5. Pathway enrichment analysis was also undertaken, based on differentially-abundant KO IDs from pairwise comparisons between the control and the treatments for the bloom period, as well as the post-bloom period.
Temporal differences in functional differences were calculated, based on the average similarity (1 -Bray Curtis) of the relative abundance of the KO IDs between the control treatment and the nutrient additions for each sampling day. The temporal response of select metabolic pathways was investigated. The KEGG IDs that were used for this analysis are given in Suppl. material 7: Table S2. Average RPMs per pathway were calculated by summing all KEGG IDs within the pathway of interest for each sample and then finding the average either per day or per period.

Reads, sequencing depth, contigs
A total 1,583,978,528 reads passed filtering and were used for SqueezeMeta. On average, 70% of a sample reads mapped to the contigs with an N50 of 1,606 bp. The number of reads mapped per sample is detailed in Suppl. material 8: Table S3. Overall, 2,979,235 ORFS were predicted using Prodigal (sample mean = 799,662, sd = 231,854; Suppl. material 8: Table S3 shows the numbers per sample). Of the almost three million ORFs, 40.4% could be annotated using the KEGG Orthology database and were attributed to 13,526 KO IDs (Suppl. material 8: Table S3 for the number of KEGG annotations per sample).

Taxonomic composition in treatments and over time
Proteobacteria dominated the metagenomic reads in all conditions (Suppl. material 1: Fig. S1) made up predominantly of Alphaproteobacteria (mean across treatments = 39.3%, sd = 9.3%) and Gammaproteobacteria (mean across treatments = 12.6%, sd = 6.4%). Cyanobacteria (mean across treatments = 14.9%, sd = 6.0%) and Flavobacteriia (mean across treatments = 5.6%, sd = 1.3%) also contributed substantially to the planktonic community. Eukaryotic taxa accounted for a higher proportion of the taxonomically-classified ORFs in the treatments with Bacillariophyceae accounting for a relative abundance of 7.0% in the NPS.cont treatment. Viruses of the class Caudoviricetes increased in relative abundance in the single addition treatments accounting for 4.7% in NP and 7.3% in NPS (Suppl. material 1: Fig. S1).
In the continuous treatments, relative abundances of Alphaproteobacteria were substantially higher during the bloom than in the pre-bloom period reaching an average of 59.4 ± 6.4% and 65.2 ± 7.0% in the NP.cont and NPS. cont, respectively (Fig. 2). The photosynthetic bacterial class Cyanobacteriia had the highest relative abundance in the pre-bloom period of the experiment, displaying a decrease in relative abundance during the bloom period, especially in the NP.cont and NPS.cont treatments where the class accounted for an average of only 1.5 ± 0.7% and 3.8 ± 3.0%, respectively (Fig. 2). Further declines were noted in the post-bloom period in the nutrient treatments with Cyanobacteriia accounting for less than 1% of the community in the NP.cont and NPS.cont treatments. The eukaryotic phytoplankton class Bacillariophyceae accounted for < 1% of the relative abundance, except during the bloom period of the NPS.cont treatment where it accounted for an average of 16.6 ± 3.4% of the community (Fig. 2). The bacteriophage Caudoviricetes increased in both the NP and NPS treatments during the bloom period with an average of 7.9 ± 2.5% and 17.7 ± 13.1% of the community, respectively before declining again towards the end of the experiment (Fig. 2). Proportional abundances per replicate are shown in Suppl. material 9: Table S4.

Functional composition
The majority (~ 92.5%) of the KO IDs were shared amongst all treatments (Suppl. material 2: Fig. S2). A total of 526 KO IDs were present in the nutrient enriched treatments, but absent from the control. The continuous treatments had the most unique sequences with thirty-one observed in the NP.cont treatment and thirty in the NPS.cont treatment. Mesocosms supplemented with silicate (NPS and NPS. cont) had a total of 75 unique KO IDs, while those mesocosms that were only enriched with nitrate and phosphate (NP and NP.cont) had 50 unique KO IDs.
Multivariate analysis indicated that there was a significant interaction between the variables treatment and bloom period in assessing the variation in the functional composition of the mesocosms (PERMANOVA; F = 2.14; p = 0.004). Canonical correspondence analysis ( Fig. 3A) showed that there was a cluster of samples from all treatments representing the start of the experiment associated with higher counts of Synechococcus. The end of the experiment was associated with higher concentrations of ammonia and total bacterial counts, while counts of nanoeukaryotes were correlated with the NPS.cont treatment. KO IDs related to photosynthesis (K08901, K08913) and carbon fixation (K18209, K01100) were shown to be associated with the bloom period of the NPS. cont treatment. A couple of sulphur cycles (K16951 and K17994) and nitrogen metabolism (K02591 and K15876) were associated with the post-bloom period of the NPS. cont treatment (Fig. 3B).

Differential abundance
In total, there were 5,789 KO IDs that had a differential abundance (padj < 0.05) between the control and at least one of the treatments (Fig. 4). In general, the differentially-abundant KO IDs had a higher abundance in the enriched treatments than the control. Of these 1,683 were related to metabolic pathways within KEGG. Comparisons between the control and the treatments showed that the continuous treatments had more enriched pathways compared to the single addition treatments with the NP treatment having no enriched pathways (Fig. 3). Biosynthesis of secondary metabolites was observed to be enriched in the NPS, NP.cont and NPS.cont treatments with nitrogen metabolism enriched in the continuous treatments (NP. cont and NPS.cont). Within the nitrogen metabolism pathway, KO IDs relating to the dissimilatory nitrate reduction to ammonia (DNRA) enzymes (NapAB and NrfAH) were differentially abundant in the NPS, NP.cont and NPS.cont treatments. Both the enzymes NarGHI and NirBD were differentially abundant in the NPS.cont treatment with the latter also significantly different in the NP.cont treatment. The first step, the conversion of nitrate to nitrite, in the assimilatory nitrate reduction was also significantly differentially abundant between the control and the NPS, NP.cont and NPS.cont treatments, although the second step showed no difference between the treatments and the control. The antenna protein pathway was enriched in the NPS.cont treatment compared to the control.

Temporal functional differences
The highest similarity between the nutrient treatments and the control was observed in the pre-bloom period with average similarity values ranging from 0.9 to 0.95. A decrease in similarity was observed in the bloom period with the largest declines observed in the treatments containing silicate with the NPS.cont treatment having a similarity value of 0.62 on day 6 and the NPS treatment having reached a low similarity value of 0.51 on day 8. The NPS.cont had an increase in similarity on day 8 before gradually declining during the post-bloom period, reaching an average similarity value compared to the control of 0.70 at the end of the experiment. For the NPS treatment, the post-bloom period showed an increase in similarity compared to the control with the average similarity being 0.72 on day 20. Regarding the NP.cont treatment, in general, similarity declined throughout the experiment, reaching an average similarity value of 0.39 compared to the control at the end of the experiment. The NP treatment, after an initial decline and then recovery in similarity compared with the control in the bloom period, showed a shallow decline in similarity compared with the control over the post-bloom period, ending up with a similar average similarity (0.71) to the NPS and NPS.cont treatments (Fig. 5).
Focussing on changes within energy metabolism pathways differential abundance analysis followed by pathway enrichment analysis was undertaken comparing the bloom period in the control to the bloom period in the nutrient treatments. It showed that oxidative phosphorylation was enriched in the NPS.cont treatment, while nitrogen metabolism was enriched in the NP.cont treatment.
Comparisons in the post-bloom period showed that the nitrogen metabolism pathway was enriched for all treatments, except the NP treatment.
Investigating the KO IDs in the nitrate reduction pathways showed a switch between assimilatory and dissimilatory nitrate reduction during the experiment, especially in the conversion of nitrite to ammonia ( Fig. 6; Suppl. material 3: Fig. S3). At the beginning of the experiment there was a rapid drop off in nitrite reductase in the assimilatory pathway (Fig. 6 and Suppl. material 3: Fig. S3), while this enzyme increased over time in the NPS, NP. cont and NPS.cont treatments in the dissimilatory pathway ( Fig. 6 and Suppl. material 3: Fig. S3).
Temporal differences were observed in the antenna proteins over the course of the experiment. Cyanobacterial antenna proteins (Allophycocyanin, Phycocyanin, Phycoerythrocyanin and Phycoerythrin) declined rapidly over the first 5-8 days, while the eukaryotic chlorophyll light harvesting complex increased, especially in the NPS. cont treatment, reaching a peak in the NPS.cont during the bloom on day 5 before declining again ( Fig. 6; Suppl. material 4: Fig. S4).
Temporal dynamics within the oxidative phosphorylation pathway showed that the F type ATPase (Bacteria) was dominant in all treatments across time with no KO IDs differentially abundant for this ATPase ( Fig. 6; Suppl. material 5: Fig. S5). KO IDs were differentially abundant in the NPS.cont for both the F type (Eukaryotes) and V type (Eukaryotes) with peaks on day 5 before declining in line with the increase in Bacillariophyta in the treatment during around day 5 ( Fig. 6; Suppl. material 5: Fig. S5).

Taxonomic changes
The development of coastal environments is having a substantial impact on marine communities. The addition of nutrients into the environment has been shown to affect the microbial community and has led to phytoplankton blooms (Lancelot et al. 1987;Paerl 1988). Metabarcoding results of the mesocosm samples used in this study showed that there was a distinct shift in the microbial communities in the continuous treatments compared with the control (Pearman et al. 2016a). This was especially evident in the NPS.cont treatment which saw a large diatom bloom. The taxonomic classification of the open reading frames, produced by the metagenomic sequencing, was in agreement in general with these results with Bacillariophyceae being more abundant in the NPS.cont treatment and especially during the bloom period. Proteobacteria were the dominant component in all treatments and had a substantial relative abundance throughout the experiment in agreement with the metabarcoding results (Pearman et al. 2016a) and is consistent with their distribution in marine systems (Eiler et al. 2009;West et al. 2016;Pearman et al. 2017). Cyanobacteria were substantial contributors at the beginning of the experiment, but their relative abundance and taxonomic composition changed over time with populations possibly controlled in a density dependent way by viruses which were especially prevalent in single addition treatments, as described for this dataset by Coello-Camba et al. (2020).
Metagenomic investigations are able to give a more holistic taxonomic view of the community incorporating taxa from viruses to eukaryotes without the potential of primer biases. This potentially gives a more robust understanding of the taxonomic composition of the community and, thus, providing better information for management decisions which rely on knowledge of taxonomic shifts. However, while the current sequencing effort was sufficient to detect the Bacillariophyta bloom, especially in the NPS.cont treatment, changes in other eukaryotic groups were hard to determine. While Pearman et al. (2016a) were able to show differences in the protists using 18S rRNA metabarcoding, these were not detected in the metagenomics. Eukaryotes tend to be orders of magnitude less abundant than prokaryotes in the marine system (Jardillier et al. 2010) and, indeed, this was the case in the current experiment (Pearman et al. 2016a; Suppl. material 6: Table S1). Therefore, with a relatively low sequencing depth, these eukaryotic taxa would not be as comprehensively sequenced, limiting the ability for contig assembly and, thus, analysis. To gain a fuller understanding of the eukaryotic dynamics during the bloom, a more substantial sequencing effort would be required.

Functional potential across treatments
Of the KO IDs that were annotated, over 90% of them were shared amongst all treatments. As a large proportion of the KO IDs are required for the cellular activity, the high proportion of shared KO IDs would be expected. While care should be taken due to the low number of replicates and, thus, lower statistical power, multivariate analysis showed that the addition of nutrients to the mesocosms and the progress of the bloom had a significant effect on the functional potential of the microbial community. This is further evidenced with approximately four percent of the KO IDs not being shared with the control, while 43% of the KO IDs were differentially abundant. The change in functional potential is in agreement with previous studies that have used metagenomic or transcriptomic approaches to assess phytoplankton blooms (Rinta-Kanto et al. 2012;Shi et al. 2012;Gong et al. 2017). Distinct differences were observed when comparing the effects of the treatments. The continuous treatments, which simulated a low-level chronic addition to the environment, showed more differentially-abundant KO IDs compared with the single addition treatments, which was akin to a high-level acute addition to the environment. This may indicate that the transformation of the additional nutrients into biomass was more efficient when persistent low concentrations of nutrients were added over time. This is supported by the fact that a higher number of enriched metabolic pathways in the continuous treatments were present.
The addition of silicate led to a more distinct community compared with the nutrient treatments containing only nitrogen and silicate. This can be explained by the greater increase in relative abundance of Bacillariophyta, which require silicate for growth (Egge and Aksnes 1992;Escaravage and Prins 2002), in these treatments and especially the NPS.cont treatment. This led to the differential abundance of eukaryotic genes, such as those involved in oxidative phosphorylation and the eukaryotic light harvesting complex which were not evident in the other treatments.

Temporal changes in functional potential
None of the treatments, including the control, had returned to a similar functional potential to that observed at the beginning of the experiment. However, with the exception of the NP.cont treatment, after initial drops in similarity compared to the control, there had been some recovery in similarity and stabilisation towards the end of the experiment. Different factors could explain why the community did not return to the same functional profile as before the bloom. Firstly, there may not have been sufficient time for the community to recover and this is especially true for the continuous treatments, which had nutrients supplied for the first couple of weeks and which brought about a larger and longer bloom. This could result in sufficient nutrients being recycled within the microbial loop (Azam et al. 1983;Fenchel 2008) preventing the return of the mesocosm to its pre-nutrient addition oligotrophic state. Secondly, although linked, is the fact that the mesocosms are semi-closed systems (only open to the air), so there is no water flow. Thus, organic matter or ammonia, for instance, are not removed from the system resulting in changes in the environmental niches present in the mesocosm. This could explain the almost total dominance of heterotrophic bacteria in the continuous and NPS treatments as they can rapidly utilise the organic matter (Obernosterer et al. 2008). Furthermore, as the mesocosms were open to the air, processes, such as evaporation (increasing salinity) and the input of dust (through sand storms), would also alter the environment within the mesocosm and could explain why the control also differed temporally. Finally, there is also the possibility that, once knocked out of its equilibrium, the ecosystem will not necessarily return to its previous state, at least in terms of species composition in the short term (Benincà et al. 2008;Dakos et al. 2009).
Distinct temporal differences were noted in the temporal patterns of nitrate reduction. Assimilatory nitrate reduction was shown to decline early in the experiment, showing a pattern similar to that of the contribution of cyanobacteria to the community. Cyanobacteria have previously been shown to undertake assimilatory nitrate reduction (Flores et al. 2005). In addition, the relationship between the reduction of assimilatory nitrate reduction and cyanobacteria proportional abundance could suggest that these taxa are responsible for the majority of assimilatory nitrate reduction in the surface waters of the Red Sea. Furthermore, assimilatory nitrate reduction is repressed by ammonia (Mohan and Cole 2007), which increased during the experiment. In contrast, dissimilatory nitrate reduction is not repressed by the increase in ammonia and was observed to increase throughout the experiment with KO IDs related to this pathway observed to be differentially abundant in the continuous nutrient treatments and especially the NPS.cont treatment. This process reduces nitrate to ammonia, recycling the fixed N in the system and is in competition with denitrification for free nitrate (van den Berg et al. 2015;Broman et al. 2021). Dissimilatory nitrate reduction requires anoxic conditions and has previously been described to happen in cyanobacterial aggregates (Klawonn et al. 2015). This could explain the increase of this process towards the end of the current experiment as the breakdown of the bloom creates anoxic micro-niches where dissimilatory nitrate reduction can occur and recycle the N within the system. Broman et al. (2021) suggested that approximately a fifth of the nitrate pool is recycled by dissimilatory nitrate reduction and has the potential to sustain algal proliferation and enhance eutrophication.
Differences were also observed in the temporal dynamics of the antenna proteins present in the mesocosms and highlights a shift in the taxa contributing to primary production with nutrient enrichment. The prokaryotic antenna proteins (allophycocyanin, phycoerythrocyanin, phycocyanin and phycoerythrin) declined in line with the drop in the proportion of the community related to cyanobacteria. The eukaryotic light harvesting complexes increased especially in the NPS.cont treatment in line with the increase in Bacillariophyta within the community, reinforcing the results showing the shifts in the photosynthetic community with additional nutrients. Similarly, changes in the oxidative phosphorylation enzymes were noted with the shift in community with an increase in the eukaryote F and V type ATPases being observed with the increase in Bacillariophyta in the community in the NPS.cont treatment.
The sampling methodology undertaken in this study allowed for a temporal investigation of the functions occurring before, during and after a nutrient-stimulated phytoplankton bloom. Understanding, the processes occurring in each stage of the bloom (e.g. changes in nitrate reduction) is vital to understand the mechanisms of a bloom progression. This knowledge could then be used in future monitoring or remediation processes to mitigate the impacts of anthropogenic impacts (Techtmann and Hazen 2016).

Conclusions
With increasing anthropogenic impact on coastal waters, it is vital to understand how the microbial community responds to environmental changes, such as nutrient addition. In this study, we show that the addition of nutrients led to both a taxonomic and functional response in the planktonic community with continuous (i.e. chronic) nutrient enrichment having a larger effect compared to a single (i.e. acute) addition. Temporal dynamics showed that, in the time frame of the experiment, communities had not returned to comparable similarity levels to those obtained in the control, with the exception of the NP treatment. Temporal dynamics in energy metabolism pathways and especially in nitrate reduction and photosynthetic antenna proteins were evident as the bloom progressed and the nutrient conditions in the bloom progressed. The demonstration that there were functional shifts in the community as well as taxonomic shifts can provide management agencies a better understanding of how to respond to anthropogenic impacts. Understanding the mechanisms of the processes occurring in the community has the potential to improve decisions to reverse or mitigate the impacts from human impacts of coastal environments (Techtmann and Hazen 2016).