Cyanobacterial taxonomy and diversity
Cyanobacteria were detected in all metagenome and RNA-seq sediment samples at all stations (n = 3 per station for both DNA and RNA data). The number of classified reads belonging to cyanobacteria was on average 3.3% (range 2.5–4.5%) in the DNA data, and showed a decreasing trend along the water depth (i.e. oxic station A to anoxic station F) in the RNA data from 7.9 to 4.2% (Additional File 1: Fig. S1). Both the cyanobacterial DNA and RNA datasets were dominated by populations belonging to the oxygenic photosynthetic class Oxyphotobacteria, with especially two taxonomic orders: 1) Synechococcales (up to 65% of the cyanobacteria in the RNA data) and 2) Nostocales (up to 44% of the cyanobacteria in the RNA data; Fig. 1). In the DNA data there was an increasing pattern of reads classified to Nostocales, from 13 to 33% of all cyanobacteria, along the water depth gradient (i.e. station A to F; Fig. 1). Looking closer at the cyanobacteria in the hypoxic-anoxic sediment (station D, E and F), most of the classified reads were attributed to 1) the genera Cyanobium and Synechococcus both belonging to the order Synechococcales (both in the DNA and RNA data), and 2) the genus Anabaena belonging to the order Nostocales (15–20% of all cyanobacteria; mainly in the RNA data; Fig. 2). Cyanobium was attributed more DNA and RNA-seq classified reads at the oxic station A compared to the hypoxic-anoxic stations D, E and F, with DNA 51% and RNA 14% of all cyanobacteria compared to DNA 14% and RNA 5%, respectively. Synechococcus showed the opposite pattern with more numbers of reads at station D, E, and F (DNA 19% and RNA 16% of all cyanobacteria compared to station A with DNA 12% and RNA 11%; Fig. 2). These results show that cyanobacteria, even bloom forming oxygenic photosynthetic taxa, are present and possibly active (as indicated in the RNA transcripts data) in the dead zone sediment. A full list of read counts and sequence classifications against cyanobacterial genomes in the RefSeq database is available in Additional File 2: Data S1.
Cyanobacterial alpha diversity (Shannon’s H) in the sediment was significantly higher in the oxic station A (Shannon’s H 4.17 ± 0.02, one standard deviation shown) compared to the hypoxic-anoxic stations D, E, and F (Shannon’s H 3.92–3.94; One-Way ANOVA, F(3,8) = 129.76, P = 0.0000004; and Tukey post hoc tests P < 0.00001; Fig. 3).
Metabolism of cyanobacteria in the dead zone sediments
We detected more RNA transcripts attributed to cyanobacteria at the oxic station A (363 ± 45 GeTMM) compared to the hypoxic-anoxic stations D, E, and F (170 ± 37, 158 ± 42, 126 ± 21, respectively; F = 24.9, P < 0.001, One-Way ANOVA with Tukey Post hoc test; Fig. 4A). RNA transcripts for the KEGG MODULES fatty acids biosynthesis and glycolysis were prevalent at all stations, indicating that anaerobic carbon metabolism such as fermentation was likely occurring (with only “Fatty acid biosynthesis, elongation” being significantly higher at station A compared to D; F = 4.9, P < 0.05). Other transcribed processes included e.g. KEGG MODULES for the reductive citrate cycle (Arnon-Buchanan cycle), cystine amino acid synthesis (no difference between stations), and the photosystems which were indicated in the RNA transcript data at all stations (Fig. 4A). In more detail, The Arnon-Buchanan cycle was higher at station A compared to E and F (F = 4.1, P < 0.05). The Photosystem I had more RNA transcripts at station A (F = 85.7, P < 0.001). Interestingly, the Photosystem II had a different number of RNA transcripts for each station (when each station was tested between each other (F = 31.3, P < 0.05; Fig. 4A). Nitrate assimilation was prominent only at station A (F = 30.8, P < 0.001; Fig. 4A). Because not all proteins could be classified into KEGG MODULES, we also decided to inspect the UniProtKB/Swiss-Prot classified protein data manually. A full list of proteins affiliated with Cyanobacteria is available in Additional File 2: Data S2, KEGG KOs and Modules in Additional File 2: Data S3, while a list of proteins with a significant difference between stations are available in Additional File 2: Data S4.
RNA polymerase was present at all stations in the cyanobacterial RNA data (Additional File 2: Data S2), also the oxic station A had more RNA transcripts classified to ribosomal proteins (30S and 50S) when compared to the other stations in hypoxic and anoxic conditions (edgeR analysis, FDR < 0.05). Furthermore, the number of significantly different proteins between station A and D, E, and F increased with water depth, with: A vs D 58 proteins, A vs E 69 proteins, and A vs F 90 proteins (Additional File 2: Data S4). This shows that with decreasing oxygen concentrations, and further away the oxic station A, the difference in cyanobacterial metabolism increases.
Cyanobacterial RNA transcripts for N2 fixation (genes nifDEHKNU) were present at all stations (Fig. 4B). The number of nifU RNA transcripts was higher at station A compared to D and F (FDR < 0.05) and NifH RNA transcripts was higher in station E compared to A (FDR < 0.05). Genes related to assimilatory nitrate reduction and uptake (nitrate reductase narB and nitrate transporters nrtABCDP) had together more RNA transcripts at the oxic station A than in the hypoxic anoxic stations (20–30 GeTMM at station A compared to < 1 GeTMM at the other stations; Fig. 4B). Of these, RNA transcripts attributed to genes narB and nrtP were significantly higher at station A (FDR < 0.05) than in station D, E, and F (Additional File 2: Data S2).
RNA transcripts attributed to genes affiliated with the photosystems decreased with water depth and O2 concentrations (photosystem I: ycf3, ycf4, ycf12, btpA, psaABCDEFIJKLM, and photosystem II: psbA1A2A3BDHIJKLMOTUXYZ, psb27, psb28; Fig. 4B). There was a variation between stations with e.g. psbD and psbA2 having more RNA transcripts at station A compared to the others (FDR < 0.05), and psaD showing the opposite pattern with lower RNA transcripts in station A (FDR < 0.05). Furthermore, the psbA123 genes together had the highest GeTMM values in the mRNA dataset except at station F, with an average of 86 ± 4% of photosynthesis GeTMM-values at A, D, E and 27 ± 4% at F (One-Way ANOVA, F = 349.4, P < 0.001). In addition to the photosystem, cyanobacterial phytochromes used to detect light were attributed RNA transcripts at all stations with no significant difference (phyAB; Fig. 4B).
The hypoxic-anoxic stations D, E, and F had a higher amount of RNA transcripts attributed to gas vesicle proteins (range 52–96 GeTMM) when compared to the oxic station A (< 3 GeTMM, genes gvpAJN; Fig. 4B and Table 2), with RNA transcripts for gvpA being dominant of the gvp proteins and significantly higher at stations D, E and F (FDR < 0.05; Additional File 2: Data S2).
Table 2
RNA transcripts attributed to proteins with a significant higher difference (edgeR analysis, FDR < 0.05) in the hypoxic/anoxic sediments of station D, E, and F when compared to station A. The values in the table shows the log fold change difference for stations D, E, and F compared to station A.
Protein name | D | E | F |
3',5'-cyclic adenosine monophosphate phosphodiesterase CpdA | | 4.1 | 3.7 |
30S ribosomal protein S18 | 5.9 | | 5.6 |
3-oxoacyl-[acyl-carrier-protein] reductase | | | 5.2 |
Acetyl-coenzyme A synthetase (AcCoA synthetase) (Acs) | | 3.0 | |
Acyl carrier protein (ACP) | | 3.1 | 3.1 |
Adenylate cyclase | | 1.9 | |
ATP synthase gamma chain | 5.5 | 4.5 | 4.4 |
ATP-dependent DNA helicase RecG | 2.0 | | |
Carbon dioxide-concentrating mechanism protein CcmK homolog 2 | | | 2.1 |
Carbonic anhydrase | 3.9 | 3.8 | 3.5 |
Chaperone protein ClpB 1 | | | 4.4 |
C-phycoerythrin class 1 subunit alpha | | | 5.1 |
Demethyl-4-deoxygadusol synthase (DDGS) | 5.4 | 5.1 | 4.9 |
DNA-binding protein HU | | | 5.9 |
Elongation factor Tu (EF-Tu) | | | 5.4 |
Gas vesicle structural protein (GVP), gvpA | 6.2 | 5.4 | 5.9 |
Glucokinase | | 2.6 | |
Glutathione synthetase | 1.8 | 1.9 | |
HTH-type transcriptional activator CmpR | | 1.5 | |
Iron uptake protein A2 (Iron deficiency-induced protein A) (PP2) | | | 5.2 |
Metalloprotease slr0863 | 3.6 | | 3.9 |
Nitrogenase iron protein | | 2.1 | |
Photosystem I reaction center subunit II (PSI-D) | 5.9 | 5.1 | 5.7 |
Phytochrome-like protein cph1 | | 1.6 | |
Probable transposase for insertion sequence element IS701 | 1.4 | 1.4 | 1.4 |
Protein GvpF/L | 5.9 | | 6.6 |
Protein RcaC | | 7.3 | |
Protein Smf | 2.4 | 2.0 | 2.1 |
Putative biopolymer transport protein ExbB-like 1 | 5.4 | | |
Putative dihydroflavonol 4-reductase (DFR) | | | 2.0 |
Putative pyruvate-flavodoxin oxidoreductase | | | 1.4 |
Putative RNA-binding protein RbpA | 2.3 | 2.1 | 1.6 |
Pyruvate kinase 1 (PK 1) | 6.3 | | 5.8 |
Pyruvate kinase 2 (PK 2) | | 5.3 | 5.6 |
Ribonuclease 3 | 1.7 | 2.7 | 2.4 |
Ribosome hibernation promotion factor (HPF) (Dark-response protein 21) | 2.8 | 2.8 | 2.2 |
Ribosome hibernation promotion factor (HPF) (Light-repressed protein A) | | | 2.7 |
RNA 3'-terminal phosphate cyclase (RNA cyclase) | | | 5.5 |
SsrA-binding protein (Small protein B) | | | 5.6 |
Thiamine-phosphate synthase (TP synthase) (TPS) | 3.8 | 3.9 | 3.5 |
TPR repeat-containing protein SYNPCC7002 A0425 | | | 5.8 |
UDP-N-acetylmuramoyl-tripeptide–D-alanyl-D-alanine ligase | 4.8 | 3.5 | 3.5 |
Uncharacterized hydrolase sll0100 | 3.7 | | |
Uncharacterized methyltransferase slr0309 | | | 1.7 |
Uncharacterized phycocyanin operon protein Z (ORF Z) | | | 2.1 |
Uncharacterized protein sll0925 | | 6.4 | |
Uncharacterized protein sll1119 | | | 5.4 |
Uncharacterized protein sll1483 | 1.9 | 1.6 | 1.8 |
Uncharacterized protein slr0039 | | | 1.9 |
Uncharacterized serine protease syc0938 d | | | 2.5 |
In addition to these metabolic processes, compared to the oxic station A the deepest anoxic station F had more RNA transcripts attributed to proteins used in e.g. iron deficiency (Iron uptake protein A2), stress related protein Chaperone protein ClpB, Pyruvate kinase proteins used in glycolysis (PK1 and PK2), and ribosome hibernation to stop protein synthesis (dark-response and light-repressed hibernation promoting factors) (all FDR < 0.05; Table 2).
Cyanophage community composition and diversity
A total of 90 partial cyanophage genomes (i.e. contigs) were assembled from the co-assembly metagenome data with an average genome size of 5076 bp (range 2000–47,157 bp). These cyanophages were used in network analyses conducted using vConTACT [42], which clusters genome-wide proteins based on the MCL algorithm, and showed that the cyanophages formed three clusters (Fig. 5). These clusters were: x similar to Siphoviridae cyanophages, y similar to Podoviridae cyanophages, and z similar to Myoviridae cyanophages (according to the cyanophage taxonomic classifications; Fig. 5). See Additional File 1: Fig. S2 for the whole vConTACT network analysis including all reference viruses. Of these 90 identified cyanophages 77 had a mapped read coverage of at least 75% and were therefore considering to be present in the sediment at the stations (Fig. 6 and Additional File 2: Data S6 for a full list of these 77 cyanophages). Cyanophage alpha diversity was 2.37 ± 0.11 (station A, Shannon’s H), 3.07 ± 0.11 (D), 3.51 ± 0.43 (E), and 3.11 ± 0.39 (F). The alpha diversity was significantly higher at station E when compared to station A (P = 0.007, One-Way ANOVA with Tukey post hoc test, whole model F(3,8) = 7.39, P = 0.011; Fig. S3). The community composition of the cyanophages was different in oxic sediment at Station A when compared to the hypoxic-anoxic sediment at stations D, E, and F (Bray-Curtis beta diversity, PERMANOVA (9999 permutations), F(3,8) = 5.71, P = 0.0006; Additional File 1: Fig. S4). These differences in cyanophage alpha and beta diversity is also visible in Fig. 6, and interestingly most cyanophages with a high relative abundance in the oxic sediment at station A belonged to cluster x, while cyanophages in the hypoxic-anoxic sediment mainly belonged to cluster y and z (Fig. 6).
Cyanobacteria and associated cyanophages in the dead zone sediments
The community composition of the cyanophages in the sediment were found to largely cluster with stations D and E. Similarly the composition of cyanobacteria formed three clusters: stations D + E, A, and F, which was in relation to the oxygen concentrations at the stations (RNA data; Fig. 7A). The cyanophages distributed in the sediments also showed an association with the cyanobacteria proteins. Here, cyanophages x16, x27, x57, z7, z79, z81 occurring at the oxic station A (Fig. 6) clustered towards the beta diversity of the UniProtKB/Swiss-Prot classified cyanobacteria proteins (mRNA data) in the sediment of station A (Fig. 7B). Majority of the cyanophages and cyanobacteria proteins clustered with the hypoxic-anoxic stations (right side in Fig. 7B). BVSTEP analysis showed that the combination of the three cyanophages y2, x16, and z58 best explained the beta diversity of the cyanobacteria proteins (rho = 0.80). That the cyanophages had an effect on the cyanobacterial community in the hypoxic-anoxic sediment was also indicated by several cyanophages (z5, z10, z12, z13, z17, z19, z45, z46, y48, z58, z74, z76, z78, z82) correlating positively with the number of mapped RNA transcripts attributed to cyanobacterial protein CRISPR-associated endoribonuclease Cas2 3 (gene cas2-3; rho = 0.60–0.86, n = 12 for each cyanophage, all P < 0.05; Fig. 7C). Finally, the protein classification of the cyanophage genomes revealed that they were carrying a large number of gene copies related to hypothetical proteins, but also e.g. glycosyl transferase related genes, heat-shock protein IbpA, p-starvation inducible protein (gene phoH), and photosystem II D1 (psbA genes) (Additional File 2: Data S5). Cyanophage genes did not explain the presence of cyanobacteria genes in the metagenome data as indicated by the lower cyanophage gene count when compared to cyanobacteria. For example psbA genes (cyanophage 2.84 ± 1.05, cyanobacteria 8.92 ± 3.07 GeTMM, One-Way ANOVA, F = 42.3, P < 0.001), and phoH (cyanophage 1.00 ± 0.41, cyanobacteria 20.33 ± 6.66 GeTMM, F = 100.9, P < 0.001; data from all stations n = 12) (Additional File 2: Data S2 and Data S5).