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

Comparative metagenomics of phytoplankton blooms after nutrient enrichment of oligotrophic marine waters John K. Pearman1,2, Laura Casas2,3, Craig Michell2,4, Naroa Aldanondo2,5, Nazia Mojib2,6, Karie Holtermann2,7, Ioannis Georgakakis2, Joao Curdia2, Susana Carvalho2, Amr Gusti2, Xabier Irigoien2,5,8 1 Cawthron Institute, Coastal and Freshwater Group, 98 Halifax Street East, Nelson, 7010, New Zealand 2 King Abdullah University of Science and Technology (KAUST), Red Sea Research Center, Thuwal, 23955-6900 Saudi Arabia 3 Institute of Marine Research (IIM-CSIC), Vigo, Spain 4 University of Eastern Finland, Department of Environmental and Biological Sciences, Yliopistokatu 7, Joensuu, 80101, Finland 5 AZTI – BRTA, Marine Research, Herrera Kaia, Portualdea z/g – 20110 Pasaia (Gipuzkoa), Spain 6 Spelman College, 350 Spelman Lane, S.W., Atlanta, GA 30314-4399, USA 7 City of Unalaska, P.O. Box 610, Unalaska, AK 99685, USA 8 IKERBASQUE, Basque Foundation for Science, Bilbao, Spain Corresponding author: John K. Pearman (john.pear[email protected]g.nz) Academic editor: Uwe John | Received 11 December 2021 | Accepted 30 March 2022 | Published 12 April 2022 Abstract 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 understanding of the functional and taxonomic shifts in the community during blooms is essential to reverse or mitigate human impacts on coastal environments. Key Words algal blooms, eutrophication, metagenomics, nutrient enrichment Introduction Ecological and biogeochemical processes in the marine system are driven by a diverse assortment of microbial organisms including members from the kingdoms Archaea, Bacteria and Eukarya. Microbial plankton are responsible for a variety of processes, including carbon fixation (Strom 2008; Jardillier et al. 2010), biogeochemical cycling (Falkowski et al. 2008; Strom 2008) and the transfer of energy between trophic levels (Weisse 1989). The biogeochemical cycling is dependent on the functional capacities of the community to cycle different elements (Louca et al. 2016). Copyright John K. Pearman et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Metabarcoding and Metagenomics 6: 101– 115 DOI 10.3897/mbmg.6.79208 Research Article https://mbmg.pensoft.net John K. Pearman et al.: Metagenomics algal blooms102 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. 2015; Pearman et al. 2016b, 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 (Sunagawa et al. 2015) 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 phytoplankton 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. (2016a, 2016b, 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. Methods Experimental design Mesocosm bags of 8000 l (depth: 2.5 m) were placed in the harbour of King Abdullah University of Science and Figure 1. Framework showing the experimental set-up used to evaluate taxomic and functional changes of the microbial community throughout a nutrient-enriched mesocosm phytoplankton bloom. Metabarcoding and Metagenomics 6: e79208 https://mbmg.pensoft.net 103 Technology (KAUST), Thuwal, Saudi Arabia (22.304°N, 39.103°E). Sampling took place at solar noon over 20 days between the 27 January and 15 February 2013 (Suppl. material 6: Table S1). Mesocosm bags (8000 l) were subjected to five different treatments (2 replicates per treatment): 1) a single addition of nitrate (NaNO3; 16 µM) and phosphate (H2NaPO4.H2O; 1 µM), referred to as the NP treatment; 2) a single addition of nitrate (16 µM) and phosphate (1 µM) and silicate (Na2SiO3.9H2O; 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 SqueezeMeta 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 https://mbmg.pensoft.net John K. Pearman et al.: Metagenomics algal blooms104 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 VennDiagram 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. Results 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 Metabarcoding and Metagenomics 6: e79208 https://mbmg.pensoft.net 105 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 Figure 2. The relative abundance (%) of taxonomically-classified ORFs (only those with a kingdom level classification are depicted) at the phylum level by treatment and day of the experiment. Other includes taxa not classified at the phylum level and lower abundance phyla. NP = the single nitrate and phosphate addition; NP.cont = the continuous nitrate and phosphate addition; NPS = the single nitrate, phosphate and silicate addition; and NPS.cont = the continuous nitrate, phosphate and silicate (NPS.cont) addition. https://mbmg.pensoft.net John K. Pearman et al.: Metagenomics algal blooms106 (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 Figure 3. Canonical correspondence analysis, based on the ordination of the KO IDs and constrained by environmental variables (A) and some of the outermost KO IDs related to energy metabolism (B). Points are coloured by the bloom period with shapes indicating different nutrient treatments. NH4 = Ammonium; NO2; = Nitrite; SiO4 = Silicate, PO4 = Phosphate; Syns = Flow cytometry counts of Synechococcus; Pico = Flow cytometry counts of picoeukaryotes; Nano = Flow cytometry counts of nanoeukaryotes and Bact = Flow cytometry counts of total bacteria. NP = the single nitrate and phosphate addition; NP.cont = the continuous nitrate and phosphate addition; NPS = the single nitrate, phosphate and silicate addition; and NPS.cont = the continuous nitrate, phosphate and silicate (NPS.cont) addition. Metabarcoding and Metagenomics 6: e79208 https://mbmg.pensoft.net 107 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 Figure 4. Enriched pathways in the nutrient treatments, based on significantly differentially-abundant KO IDs between the control and A) the continuous nitrate and phosphate (NP.cont) addition; B) the single nitrate, phosphate and silicate addition (NPS); and C) the continuous nitrate, phosphate and silicate (NPS.cont) addition. https://mbmg.pensoft.net John K. Pearman et al.: Metagenomics algal blooms108 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). Discussion 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 Figure 5. The average similarity based on 1 – Bray Curtis dissimilarity between the KO ID composition of the nutrient treatments compared to the control on each sampling day. Treatments are depicted in different colours. NP = single nitrate and phosphate addition; NPS = single nitrate, phosphate and silicate addition; NP.cont = continuous nitrate and phosphate addition; NPS.cont = continuous nitrate, phosphate and silicate addition. Metabarcoding and Metagenomics 6: e79208 https://mbmg.pensoft.net 109 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 Figure 6. Mean reads per million (RPM) of selected pathways during the pre-bloom, bloom and post-bloom periods. NP = the single nitrate and phosphate addition; NP.cont = the continuous nitrate and phosphate addition; NPS = the single nitrate, phosphate and silicate addition; and NPS.cont = the continuous nitrate, phosphate and silicate (NPS.cont) addition.