Viral activity in lake analogs of anoxic early Earth oceans

Abstract

Background

Meromictic lakes, with their stratified water columns, are modern analogs for ancient euxinic (anoxic and sulfidic) oceans, where anaerobic sulfur-oxidizing purple and green sulfur bacteria (PSB and GSB) dominated as primary producers. Recent studies suggest a potential role of viruses in the metabolisms and biosignatures of these bacteria, but conclusive evidence of viral replication and activity in such lakes is still lacking.

Results

Here, we investigate viral activity in the upper mixed layer (mixolimnion), the anoxic bottom (monimolimnion), and the microbial plate (a dense layer of phototrophic sulfur bacteria forming at the boundary between the oxygenated mixolimnion and the anoxic monimolimnion) of three meromictic lakes: Poison and Lime Blue Lakes (WA, USA) and Mahoney Lake (BC, CA). Geochemical profiles of two lakes, Mahoney and Poison, which are dominated by PSB, show a sharp chemocline, whereas Lime Blue displays a less steep chemical gradient and hosts a mixture of PSB and GSB. Viral gene transcription and epifluorescence microscopy revealed depth-dependent patterns in viral activity. The two strongly stratified, PSB-dominated lakes showed a significant decrease in the virus-to-microbe ratio (VMR) in their microbial plates, suggesting reduced viral particle production via lysis. Metatranscriptome data corroborated this trend by showing lower levels of viral gene expression in these microbial plates, higher expression of CRISPR defense and lysogeny-related genes, and relatively high expression of photosynthesis-related viral genes. Conversely, the third lake, which harbors a mix of PSB and GSB, exhibited low microbial density, high VMR, and high viral transcriptional activity. Viral transcription levels significantly correlated with VMR in the microbial plates and bottom layers, but this relationship was absent in low-density, oxic surface samples.

Conclusions

Here, two independent lines of evidence, abundances and gene expression, show reduced viral lytic production in microbial plates dominated by PSB in stratified lakes. This suggests that viral lysis may contribute less to bacterial community structuring in these high-density microbial plates. Rather, other viral-mediated mechanisms, such as lysogeny and the expression of auxiliary metabolic genes, may represent a more significant viral influence on bacterial physiology and geochemistry. These patterns in virus-bacteria interactions may be consequential for the interpretations of biosignatures left by these bacterial groups in the geologic record.

Video Abstract

Introduction

Meromictic lakes consist of highly stratified water columns considered analogs of ancient oceans [1]. This stratification typically comprises an oxygenated surface layer (mixolimnion), a microbial plate associated with sharp chemical gradients (chemocline), and an anoxic bottom layer (monimolimnion). The chemocline is marked by a rapid decline in oxygen and a corresponding rise in hydrogen sulfide concentrations, creating euxinic conditions. In these euxinic waters, light penetration supports the development of a unique microbial community dominated by anaerobic sulfur-oxidizing purple and green sulfur bacteria (PSB and GSB, respectively). These bacteria use sulfide for anoxygenic photosynthesis and can form syntrophism with sulfate reducers [2, 3]. Therefore, the microbial plate supports a microbial community that is distinct from the community in the mixolimnion and the monimolimnion [4,5,6,7,8]. Modern euxinic lakes offer valuable insights into the ecological conditions during periods of ocean anoxia in Earth’s history, where anoxygenic photosynthetic bacteria, such as PSB and GSB, are believed to have played important roles as primary producers [9, 10]. This hypothesis is supported mainly by the detection of biomarkers produced by PSBs and GSBs, such as carotenoids, in ancient sedimentary rocks [11, 12]. These biomarkers, along with distinct carbon and sulfur isotopic ratios resulting from their metabolic processes, serve as biosignatures of past or present microbial life with specific activities and environmental conditions, such as light availability, sulfide requirements, and oxygen sensitivity.

Recently, viruses have been proposed to significantly affect PSB and GSB community dynamics [13, 14]. Specifically, viruses could affect biosignatures by modulating community composition via lysis and interfering with the rates at which these bacteria make their metabolic products [13]. Virus-mediated changes in metabolic pathways could decouple microbial biosignatures from environmental factors such as light and sulfide concentration, and confound the interpretation of these markers. This may help explain inconsistencies observed between the abundance and activity of extant bacterial phototrophs and their geochemical signatures in meromictic lakes. For instance, a growing body of evidence suggests that the abundance of PSB does not directly explain the concentrations of okenone, a carotenoid pigment unique to PSB, nor the isotopic composition of carbon and sulfur that is expected to reflect their metabolic activities [15,16,17]. One such case is found in euxinic Fayetteville Green Lake (NY, USA), where even though GSB was the dominant group in the water column, okenone was the predominant biomarker observed [18]. In Lake Cadagno (Switzerland), one species of PSB, Chromatium okenii, made up 0.3% of the bacterial community and yet was responsible for 70% of the carbon uptake [16]. Recent work showing that viruses in euxinic lakes encode auxiliary metabolic genes (AMGs) related to carbon, sulfur, and pigment metabolisms supports the idea that these biogeochemical inconsistencies may be caused by viral infection [13]. Viruses commonly encode genes that can alter the metabolisms and behavior of the bacteria they infect [19,20,21,22,23,24]. In Lime Blue Lake (WA, USA), a meromictic lake with a chemocline harboring both PSB and GSB, viruses encode AMGs associated with the light reactions of photosynthesis (psbD, psbA), carbon fixation and recycling (CP12, rpe, rbsK, G6PD/zwf, talA/B, PGD/gnd/gntz), sulfur relay system (this, moeB, moaA, mec, iscS/NFS1, phsA/psrA, nrnA, cysH, cyst), and pigment production (crtF). However, it remains unclear whether these viruses are actively replicating or expressing metabolic genes in situ. Determining the rate at which viruses can manipulate PSB and GSB, thereby influencing their biosignatures, could have implications for interpreting Earth’s redox evolution and improve the accuracy of proxies for euxinic water columns. Additionally, viral predation and gene transfer may accelerate the rate of cellular evolution by facilitating adaptation to euxinic environments.

Here, we investigate viral activity in three euxinic lakes by comparing 33 metatranscriptomes and 43 epifluorescence microscopy counts across three water column layers (the mixolimnion, the microbial plate within the chemocline, and the monimolimnion) informed by the lakes’ geochemical profiles. We show that viral activity in these lakes significantly differs across depths and between PSB-dominated and GSB-dominated lakes.

Methods

Sample collection

The three study lakes were sampled from July 19 to July 25, 2023 (see Table S1 for exact sampling locations and dates). Water column temperature, dissolved oxygen (DO), specific conductivity, and pH were measured using a multiparameter digital water quality meter (ProDSS; YSI Incorporated, Yellow Springs, OH, USA). Photosynthetically active radiation (PAR) was measured using a spherical underwater quantum sensor and data logger (LI-193 and LI-1500; LI-COR Environmental, Lincoln, NE, USA). Both sensors were cast from the side of an inflatable boat. Measurements were recorded at every meter to detect the location of the mixolimnion, chemocline, and monimolimnion and to create a vertical profile of lake water chemistry. We then increased vertical measurements to every 10 cm when noticing a decrease in oxygen to obtain a higher resolution of chemical profiles across the chemocline. We also pumped small volumes of water and passed them through a 0.45-µm filter to observe color changes in the filter as an indication of the position of the microbial plate. Samples were collected by pumping water using a peristaltic pump (Electra V2.0; Proactive Environmental Products, Lakewood Ranch, FL, USA) with a custom-assembled acrylic plate to ensure water intake laterally from the measured depths. Once the position of the microbial plate was identified, sampling depths for geochemical and microbiological analyses were determined (mixolimnion, microbial plate, and monimolimnion, respectively): 3 m, 8.8 m, and 13 m for Mahoney Lake; 3 m, 6.7 m, and 7.5 m for Poison Lake; and 3 m, 12 m, and 13.7 m for Lime Blue Lake. For viral activity estimates, two approaches were employed: epifluorescence microscopy and metatranscriptomics. Samples were collected from each lake and depth in five replicates.

Metal extraction

After collection, water samples were filtered through a 0.45-µm PCTE filter (Cytiva, Marlborough, MA, USA) and immediately acidified with 2% HNO₃ and stored in 15 mL acid-cleaned tubes (cleaning solution with 8% HNO₃ and 2% HCl; Merck, Germany). Metal concentrations were determined using an ICP–MS (Thermo Fisher Scientific, XSERIES 2, Dreieich, Germany), equipped with Collision Cell Technology (CCT), a conical chamber, and a Meinhard® concentric nebulizer. The instrument was optimized following the manufacturer’s guidelines and calibrated with 10 mg/L multi-elemental solutions (High Purity Standards, Charleston, SC, USA). Certified reference material (SPS-SW2 from LGC Standards) was used for analytical quality control. Calculated recoveries ranged from 92 to 110%, and for uncertified elements, the variation among recoveries (N = 6) was below 10%, demonstrating the robustness and reliability of the adopted procedures.

Pigment analysis

Lipid extracts were obtained from freeze-dried POM filters (one entire 142 mm GFF per depth) using a modified Bligh and Dyer ultrasonic extraction [25, 26] with methanol (MeOH)/dichloromethane (DCM)/phosphate buffer (50 mM, pH 7.4), (2:1:0.8, v/v/v), respectively [27]. The total lipid extracts were reduced to dryness under a gentle stream of nitrogen, flushed with argon, sealed, and stored in the freezer (− 20 ℃) until analysis. Pigments were resolved by reverse-phase high-performance liquid chromatography (HPLC)—tandem atmospheric pressure chemical ionization (APCI) ion trap mass spectrometry in positive ion mode. A Dionex Ultimate 3000 HPLC system with a vacuum degasser, ternary gradient pump, refrigerated autosampler, and diode array detector connected to a Thermo Scientific LTQ XL ion trap mass spectrometer controlled by Xcalibur 4.1 software. Two coupled ES Industries Spherisep ODS2 (3 μm, 80 Å, 100 mm × 2.1 mm) columns protected by a guard column were used for chromatographic separation. A 32-min ternary gradient method was adapted from a previously published method [28] (for details, see Supplementary Data). Pigments were identified by retention time, UV–vis absorbance spectra, molecular ion m/z ratios, and MSⁿ fragmentation patterns. Pigments were quantified based on UV–visible light absorption peaks integrated at the maximum absorbance wavelength for carotenoids and at specific Q_y band maxima for bacteriochlorins [29]. Response factors (f_i) for each identified pigment were calculated in reference to the response factors for pure standards (f_std), using chlorophyll a (Chl a; 7.03 absorbance units (AU)/nmol, R² = 0.999) for tetrapyrroles and lutein (13.6 AU/nmol, R² = 0.999) and beta-carotene (8.66 AU/nmol, R² = 0.999) for carotenoids. Differences in molar extinction coefficients (ε) were corrected using Eq. 1. Molar extinction coefficients in acetone for the identified pigments and standards were obtained from literature values [30,31,32].

$${\mathrm f}_{\mathrm i}\;=\;\frac{\varepsilon_{\mathit i}}{\varepsilon_{\mathit s\mathit t\mathit d}}\;\times\;{\mathrm f}_{\mathrm{std}}$$

(1)

Epifluorescence microscopy

Samples of 1 mL were collected in cryovials and fixed with paraformaldehyde (2% final concentration) for 20 min on site and kept refrigerated for 3 to 5 h until further processing [33]. Depending on densities estimated from filter color and volume taken to clog filters, 20 to 1000 µL of the filtered sample were then filtered through a 0.02-µm Anodisc (Cytiva, Marlborough, MA, USA) using a vacuum hand pump, dried, and stored on dry ice until staining in the laboratory at the University of Miami. For samples where less than 1000 µL was used, molecular grade water was added to the sample for a final volume of 1000 µL before filtration onto the Anodisc (Table S1 for volumes). Nucleic acid-containing particles on the filters were stained with SYBR Gold Nucleic Acid Gel Stain (10X final concentration; Invitrogen, Waltham, MA, USA). Excess stain was removed by wetting the Anodisc from below with two drops of 100 µL molecular grade water and gently blotting the bottom of the filter with a lint-free wipe. Filters were air-dried for 30 min without exposure to light and mounted between a glass microscope slide and coverslip with 15 µL of mounting solution (0.02-µm-filtered 1X PBS, 0.1% ascorbic acid, and 50% glycerol). Slides were visualized under oil immersion at 630 × magnification using a ZEISS Axio Imager.A2 equipped with an Axiocam 506 mono camera and the X-Cite Mini light source (Excelitas Technologies, Pittsburgh, PA, USA) and stored at − 20 ℃. Ten images were captured per slide using the Zeiss software Zen. Cells and virus-like particles were size-delimited (> or < 0.22 µm) using Zen and counted manually to obtain bacterial cells and viral particle abundances per mL of sample after correcting for the volume of sample filtered onto the Anodisc.

Transcriptomics

Water samples were directly filtered through a 0.22-µm Sterivex filter (Millipore-Sigma, Burlington, MA, USA) coupled to the peristaltic pump output. Due to high cell densities in microbial plates and monimolimnion, sampling volumes varied between layers (see Table S1). One mL of RNAlater (Thermo Fisher Scientific, Waltham MA, USA) was immediately added to metatranscriptome samples, and filters were flash-frozen in liquid nitrogen on site. RNA extraction followed a Sterivex extraction protocol [34]. This protocol modifies a DNeasy PowerWater Sterivex kit (Qiagen, Germantown, MD, USA) by adding B-Marcaptoethanol and Invitrogen Turbo DNase. Twelve samples yielded too little RNA for sequencing. However, we had backup raw lake samples (frozen at -80 ℃) for some of these and were able to extract sufficient RNA from two of them using a manual Chloroform/TRIzol Reagent (Invitrogen, Waltham, MA, USA) extraction, leading to a total of 33 RNA samples. RNA was quantified using a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, USA), and RNA integrity was measured using the RNA Screen Tape on an Agilent 2200 TapeStation (Agilent Technologies, Palo Alto, CA, USA). To remove DNA contaminants, samples were treated with TURBO DNase (Thermo Fisher Scientific, Waltham, MA, USA). rRNA depletion was performed using QIAseq® FastSelect™ − rRNA 5S/16S/23S kit (Qiagen, Germantown, MD, USA) following the manufacturer’s protocol. RNA sequencing libraries were constructed with the NEBNext Ultra II RNA Library Preparation Kit for Illumina following the manufacturer’s recommendations. Briefly, enriched RNAs are fragmented for 15 min at 94 ℃. First-strand and second-strand cDNA were subsequently synthesized. cDNA fragments were end-repaired and adenylated at 3’ends, and universal adapters were ligated to cDNA fragments, followed by index addition and library enrichment with limited cycle PCR. Sequencing libraries were validated using the Agilent Tapestation 4200 (Agilent Technologies, Palo Alto, CA, USA) and quantified using Qubit 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) as well as by quantitative PCR (KAPA Biosystems, Wilmington, MA, USA).

The sequencing libraries were multiplexed and clustered onto a flow cell on the Illumina NovaSeq instrument according to the manufacturer’s instructions. The samples were sequenced using a 2 × 150 bp paired-end (PE) configuration at a target depth of 20 M reads per sample. The NovaSeq Control Software (NCS) conducted image analysis and base calling. Raw sequence data (.bcl files) generated from Illumina NovaSeq was converted into fastq files and de-multiplexed using Illumina bcl2fastq 2.20 software. One mismatch was allowed for index sequence identification.

Gene expression analyses

Gene expression analyses were adapted from Coclet et al. [35]. Raw reads were quality-filtered using BBDuk v.38 [36]. The adapters were trimmed, and reads annotated as ribosomal via k-mer matching to the RNA Silva 132_99 database were removed [35]. Reads containing four or more N bases with an average quality score of less than ten or a length of 51 bp or less were removed [35]. Reads were assembled by sample using Spades v3.15.5 at default with the –meta flag option, resulting in 9,999,946 contigs with an average length of 303.5 bp (min 55 bp, max 41,336 bp) [37]. Reads per kilobase million (rpkm) values and coverage statistics were calculated using BBMap v.38 [36]. Contigs with a minimum of 200 bp were annotated using MetaCerberus using all available HMMs [38]. Lysogenic and CRISPR-associated genes were identified from MetaCerberus annotations if they contained the keywords “CRISPR-Cas” or “CRISPR-associated protein” for CRISPR genes and “transposase,” “integrase,” “recombinase,” and “excisionase,” for lysogenic genes. If a gene contained both terms, CRISPR keywords took precedence. The rpkm values were normalized by rpkm per sample to allow comparisons of CRISPR and lysogeny-related genes across samples. Viral contigs were identified with geNomad v1.8.0 using end-to-end configuration [39]. Viral contigs were only considered active if reads mapped to at least 10% of the viral contig [35]. Viral contigs were dereplicated using CheckV anicalc.py and aniclust.py scripts at MIUViG of 95% ANI and a minimum of 85% coverage [40, 41]. Only viruses with no host gene contamination, as defined by CheckV, were used for downstream processes. Statistical analysis was carried out in RStudio version 2022.12.0 + 353. Rpkm values were clustered using pheatmap v1.0.12 with default parameters, which uses an agglomerative hierarchical clustering method with Euclidean distances. Viral diversity at the class level was calculated using the Shannon Index from the R package vegan v.2.6–4 based on the rpkm values at the class level. For statistical comparisons, the R package rstatix v 0.7.2 was used to test the normality of the distribution of data using Shapiro–Wilk test, and statical differences using one-way ANOVA and Tukey’s HSD post hoc test were performed using the R package stats v4.2.2.

Results

Water column stratification

Three meromictic lakes, Lime Blue and Poison Lakes in Washington, USA, and Mahoney Lake in British Columbia, Canada, were sampled to determine viral activity (Fig. 1a). The mixolimnion, chemocline, and monimolimnion were detected through geochemical vertical profiling of dissolved oxygen (DO), salinity, pH, turbidity, temperature, and light. A sharp decrease in DO was observed in Mahoney Lake from the maximum at 4 m (248.0%) to 1% at 9.5 m and at Poison Lake (from 205.5% at 5 m to 1.4% at 7.1 m), whereas DO in Lime Blue Lake decreased more slowly from approximately 5 m (135.8%) to 12.5 m (1.3%) (Fig. 1b). Mahoney Lake had the highest difference in DO, with a 248-fold decrease in the chemocline (Fig. 1b, Table S2), followed by Poison Lake with a 146-fold reduction, albeit over a shorter distance (3.5 m vs. 2.1 m, respectively). The lake with less sharp stratification, Lime Blue Lake, displayed a 103-fold decrease in DO over 7.50 m. Salinity varied widely in Mahoney and Poison Lake with a total increase from 11,298 µS/cm on the surface to 60,090 µS/cm at the bottom in Mahoney and from 6077 to 43,823 µS/cm in Poison Lake (Fig. 1b). Along the depth profiles, all lakes showed a drop in pH, the most pronounced in Poison Lake where pH dropped from 8.9 at the surface (3.00 m) to 7.2 at 9 m (Fig. S1a). Lime Blue Lake was mainly stratified by temperature, ranging from 25.5 to 12.0 ℃ (Fig. S1b), while Mahoney and Poison Lake exhibited strong salinity stratification. The microbial plate in each lake was detected through changes in 0.45-µm filter coloration. For Mahoney and Poison Lakes, the location of the microbial plate, just at the end or below the chemocline, was indicated by a bright pink coloration of the filter due to the dominance of PSBs (at 8.8 m and 7.5 m depth, respectively; Fig. 1c). The mixolimnion was generally light in color, the microbial plate pink to purple, and the monimolimnion brown. In Lime Blue Lake, no obvious microbial plate could be identified, except for a slight increase in the green–brown coloration of the filters at a depth of 12 m, concurrent with an increase in turbidity, where the sample identified as the microbial plate was chosen for subsequent analyses.

Fig. 1

Sampling sites and euxinic lake profiles. a Geographical location of lakes sampled in this study: Mahoney Lake (49°17′19″ N, 119° 34′58″ W), Lime Blue Lake (48°34′01″ N, 119° 36′45″ W), and Poison Lake (48°32′05″ N, 119° 34′36″ W). b A vertical profile of the three lakes, including dissolved oxygen (DO), specific conductivity, and Photosynthetically Active Radiation (PAR). Horizontal lines indicate the sampling depth for the mixolimnion, microbial plate, and monimolimnion at each lake. c Visual coloration of filters from the mixolimnion, microbial plate, and monimolimnion samples from Mahoney Lake

Stratification of photosynthetic pigments

Mahoney Lake and Poison Lake both had maximum pigment concentrations, including mostly okenone and bacteriochlorophyll a (Bchl a) derivatives, in the microbial plate (Table 1). This confirms that PSBs are the predominant photosynthetic organisms in this sample. In Mahoney Lake, Bchl a was the main tetrapyrrole pigment, and there was a lower concentration of its oxidized derivative hydroxy-Bchl a. In contrast, hydroxy-Bchl a was dominant in Poison Lake, and Bchl a was below the level of quantification in the microbial plate sample. The Lime Blue Lake microbial plate sample also included okenone, but isorenieratene from GSB was more abundant. Its occurrence extended downward below the chemocline and was more concentrated in the monimolimnion. This apparent stratification of PSB above GSB in meromictic lakes is common, including in Fayetteville Green Lake [29]. The GSB tetrapyrrole pigments Bchl d and Bchl e, as well as bacteriopheophytins c, d, and e, were only detected in Lime Blue Lake, and their concentrations were greater in the monimolimnion compared with the microbial plate. Our reported concentrations consist of the sum of all structural homologs for each pigment type. The mixolimnion was omitted from Table 1, as none of the pigments was detected in these samples.

Table 1 Pigment concentrations (µg/L). Absent values indicate that the pigment concentration was below detection thresholds

Low virus-to-microbe ratio in microbial plates of strongly stratified lakes

In Mahoney and Poison Lakes, both cell (Table S3 and Fig. 2a) and viral abundances (Fig. 2b) varied significantly between each of the three layers (p < 0.01, Tukey HSD test), increasing with depth. In contrast, cell abundance in Lime Blue Lake was lower at the microbial plate relative to the other lakes and did not vary significantly with depth (p > 0.05, Tukey HSD test). Viral abundance in Lime Blue Lake’s mixolimnion was significantly lower than abundances at and below the microbial plate (p < 0.01, Tukey HSD test).

Fig. 2

Microbial densities across depths at each lake. a Cell abundances. b Viral-like particle (VLP) abundances. c Virus-to-microbe ratios (VMR). Comparisons were made using the Tukey HSD test, where * indicates p < 0.05, **p < 0.01, ***p < 0.001, and ns indicates not significant

Virus-to-microbe ratios (VMR) were calculated by dividing the viral abundance by the number of microbial cells. A high VMR is a proxy for lytic viral production, whereas a low VMR could be associated with lysogeny [42,43,44,45]. Mahoney and Poison Lakes showed a significant decrease in VMR in the microbial plate, indicating decreased viral lytic activity relative to the mixolimnion (Fig. 2c). Monimolimnion samples also had significantly lower VMR than the surface water samples. This pattern was the opposite for Lime Blue Lake, where VMR was lowest at the mixolimnion and significantly higher at and below the microbial plate.

Viral gene transcription

A total of 19,152 assembled contigs were identified as viral genomes or genome fragments (Table S4), with 11,953 representative sequences after dereplication with MIUViG standards (95% sequence identity and 85% coverage). Of these viral contigs, 83.8% were predicted to belong to the viral realm Duplodnaviria, which was dominated (99.8%) by the viral class Caudoviricetes, commonly known as tailed bacteriophages. Other abundant contigs belonged to Riboviria (6.52%), followed by Varidnaviria (1.93%), and Monodnaviria (0.05%), while 8.15% of contigs could not be classified at the realm level. Across all lakes and depths, Caudoviricetes was the viral group with the highest normalized transcript abundance (reads per kilobase million, rpkm), ranging from 68.2% of transcript reads (SE = 0.64) in the mixolimnion of Lime Blue Lake to 87.3% (SE = 1.90) in the Mahoney Lake microbial plate (Fig. 3a, Fig. S2, Table S4). The transcript reads mapping to the class Megaviricites were enriched in mixolimnion samples, with the highest representation in Lime Blue Lake (7.9%, SE = 0.38). Compared to other microbial plates, Mahoney Lake’s microbial plate was enriched in Chrymotiviricetes, making it the second most active viral group (1.8%, SE = 0.28). In Lime Blue Lake’s chemocline, the second most abundant viral class was Pisonvirecetes (1.93%, SE = 0.93), whereas in Poison Lake, it was Duplopiviricetes (1.94%, SE = 0.44). In all lakes, the mixolimnion had the highest average viral diversity, indicated by their Shannon Index, which was calculated using base e (Fig. S3, Table S5). Poison Lake’s mixolimnion had the highest diversity (H' = 1.47, SD = 0.08), followed by the monimolimnion (H' = 0.881, SD = 0.21, ANOVA, p = 0.01) and the microbial plate (H' = 0.749, SD = 0.45, ANOVA, p = 0.001). Viral diversity was lower in the microbial plate and monimolimnion for all lakes.

Fig. 3

Viral activity levels by depth and lake. a Taxonomic annotation of viral genomes and genome fragments across three sampling depths and lakes. The relative abundances of viral genomes were averaged across biological replicates. b Hierarchical clustering of viral genomes by layer and lake. In blue, viral activity is shown in log₁₀(x + 1) transformed reads per kilobase million (rpkm). The upper color bar indicates the lake depth, where Mixo. = mixolimnion, Plate = microbial plate, and Moni. = monimolimnion. The second color bar indicates the lake, abbreviated as MAH = Mahoney Lake, POI = Poison Lake, and LIMB = Lime Blue Lake. The dendrogram is based on hierarchical clustering of log₁₀(x + 1) rpkm values. c Number of active viruses based on transcriptional activity in each lake across layers. Significant differences from Tukey’s HSD test are displayed with stars, where * indicates p < 0.05, **p < 0.01, and ***p < 0.001

Hierarchical clustering of total viral transcript abundance using Euclidean distance measures of rpkm values by sample revealed that viral activity clustered primarily by the depth and then by the lake (Fig. 3b). Mahoney and Poison Lakes’ microbial plates clustered together, while four of five Lime Blue Lake microbial plate samples formed a separate cluster from all other samples. We quantified the number of active viruses in each sample to compare how viral activity varied across depths within lakes. An active virus was defined as those recruiting transcripts covering at least 10% of the viral contig length [35]. To obtain the proportion of active viruses, the number of active viral contigs was normalized by the total number of contigs per sample and compared with Tukey’s HSD test (Fig. 3c, Table S7). Lime Blue Lake’s microbial plate viruses were significantly more active than those in the mix- or monimolimnion (p < 0.01 and p < 0.001, respectively). Conversely, Mahoney and Poison Lake showed significant decreases in viral activity at the microbial plate (p < 0.001). There was also a slight increase in the proportion of active viruses in Poison Lake’s monimolimnion relative to its microbial plate (p < 0.05).

CRISPR defense systems and lysogeny signatures

To investigate whether the decrease in viral activity corresponded with lysogeny or CRISPR defense systems, we investigated the relative abundance (sample-normalized rpkm) of CRISPR-associated or lysogeny-associated genes. Expression of these genes displayed significant differences across the three depths in all lakes examined (ANOVA, p < 0.001; Fig. 4). In the Mahoney Lake microbial plate, CRISPR defense gene expression increased on average by 156-fold (Tukey’s, p < 0.001), and genes associated with lysogeny increased by 15.7-fold (Tukey’s, p < 0.001) compared to the mixolimnion (Fig. 4a). Similarly, the Poison Lake microbial plate exhibited a 121-fold increase in CRISPR-associated gene expression (Tukey’s, p < 0.001) and a 19-fold increase in expression of genes associated with lysogeny (Tukey’s, p < 0.001; Fig. 4b). Between the mixolimnion and monimolimnion in both lakes, CRISPR and lysogeny genes also showed significant increases in expression, although at lower magnitudes (30-fold for CRISPR and 9.77-fold for lysogeny; Tukey’s, p < 0.001). Lime Blue Lake followed the same trend, with CRISPR gene expression increasing by 32.7-fold and lysogeny genes by 31.02-fold (Tukey’s, p < 0.001) at the microbial plate compared to the mixolimnion (Fig. 4c). However, the CRISPR defense system expression in Lime Blue Lake’s microbial plate was significantly lower compared to Poison Lake (Tukey’s, p = 0.021) and Mahoney Lake (Tukey’s, p = 0.004). Additionally, the average proportion of lysogeny-associated genes in Lime Blue Lake was lower than in Mahoney Lake and Poison Lake, although these differences were not statistically significant (Tukey’s, p = 0.072, p = 0.505, respectively).

Fig. 4

Expression of CRISPR- and lysogeny-related genes across depths. a The relative rpkm of genes associated with CRISPR defense and lysogeny with depth in Mahoney Lake, panel b in Poison Lake, and c Lime Blue Lake. *** indicates significant differences (p < 0.001) based on Tukeys HSD test

Relationship between gene expression, virus-to-microbe ratios, and chemical and physical variables

To investigate the relationship between VMR and viral activity, we compared the epifluorescence microscopy counts to the proportion of active viruses in each sample (Fig. 5). Viral-like particle (VLP) abundance displayed a weak negative relationship with viral activity in the microbial plate (Linear regression R² = 0.34, p = 0.045, slope = − 0.158; Fig. 5a), while cell abundance displayed a strong and significant negative relationship with viral activity (Linear regression R² = 0.82, p < 0.001, slope = − 0.532, Fig. 5a). VMR was a strong predictor of viral activity (Linear regression R² = 0.66, p = 0.001, slope = 0.358; R² = 0.73, p = 0.029, slope = 0.666, respectively) in both the microbial plate and monimolimnion (Fig. 5c). These trends were driven by the large density differences between the depth layers in the highly stratified lakes (Mahoney and Poison Lakes) and the lesser stratified Lime Blue Lake.

Fig. 5

Relationship between densities and viral activity from transcriptomics. a Linear regression between viral-like particle (VLP) abundance and the proportion of active viruses (mixolimnion R² = 0.30, p = 0.265; microbial plate R² = 0.34, p = 0.045; monimolimnion R² = 0.05, p = 0.442). b Linear regression between cell abundance and the proportion of active viruses (mixolimnion R² = 0.33, p = 0.235; microbial plate R² = 0.82, p < 0.001; monimolimnion R² = 0.06, p = 0.407). c Linear regression between virus-to-microbe ratio (VMR) and the proportion of active viruses (mixolimnion R² = 0.22, p = 0.106; microbial plate R² = 0.66, p = 0.001; monimolimnion R² = 0.73, p = 0.029). MAH = Mahoney Lake, POI = Poison Lake, LIMB = Lime Blue Lake, Mixo = mixolimnion, Plate = microbial plate, Moni = monimolimnion

We examined the relationship between viral activity and environmental variables using linear regression for normally distributed data (Fig. S4). Variables that met the normality threshold (Shapiro–Wilk test, p > 0.05) included molybdenum concentration, specific conductivity, turbidity, PAR, bacteriopheophytin a, and pH. Among these, specific conductivity and turbidity showed significant negative correlations with viral activity (specific conductivity: R² = 0.45, p = 0.040; turbidity: R² = 0.55, p = 0.020). Turbidity also displayed a significant negative linear relationship with viral Shannon diversity (R² = 0.44, p = 0.043). For non-normally distributed variables with at least five data entries, temperature (℃), okenone (µg/L) dissolved oxygen (%), nickel (ppb), iron (ppb), and manganese (ppb), we utilized Spearman’s rank correlation; however, no significant correlations with viral activity were observed (temperature: p = 0.793, okenone: p = 0.075, dissolved oxygen: p = 0.755, nickel: p = 0.882, iron: p = 0.556, manganese: p = 0.793). None of the above variables correlated significantly with VMR.

Highly expressed viral genes

In Lime Blue Lake and Mahoney Lake, major capsid (head) proteins and other structural genes were the most expressed viral genes (320.47 to 1029.14 rpkm) in most samples across all depths, with many serving as strong indicators for specific sample types (Fig. 6, Table S8, Table S9). Poison Lake was unique because viruses in the anaerobic layers (microbial plate and monimolimnion) primarily expressed genes associated with host metabolism. In the microbial plate of Poison Lake, the most active viral contig carried a gene annotated as a D1-like protein of the photosystem II reaction center. This gene was a significant indicator of Poison Lake’s microbial plate (indicator analysis, p = 0.048, Table S9). The most expressed viral gene in the monimolimnion encoded a gas vesicle structural protein used by prokaryotes to control buoyancy [46]. The second most expressed gene, a nitric oxide reductase NorQ protein, is a key enzyme for the assembly of catalyst cNOR for denitrification [47].

Fig. 6

Highly expressed viral genes at each lake and layer. The ten genes of viral contigs with the highest expression levels in each lake and layer (sorted by the mean across the replicates) are displayed in decreasing order, with mixolimnion samples at the top, chemocline in the center, and monimolimnion samples at the bottom. Each row displays one of three lakes, except for Mahoney Lake’s monimolimnion, which was lost during sample processing. Rpkm values are based on the contig where the gene was identified

A principal coordinate analysis (PCoA) based on rpkm abundances using Bray–Curtis dissimilarity explained 34.75% of the variance between genes expressed across sample types (Fig. S4). A PERMANOVA with 9999 permutations revealed significant differences in viral gene composition among sample types (p < 0.001). The PCoA shows Poison Lake’s and Mahoney Lake’s microbial plates clustering with Poison Lake’s monimolimnion, while all mixolimnion samples formed a cluster, with some overlap with Lime Blue Lake’s monimolimnion and microbial plate. A pairwise PERMANOVA revealed that most gene profiles of sample types were significantly different from each other (p < 0.05); however, some gene profiles were similar: Mahoney Lake and Poison Lake microbial plates (p = 0.056), Lime Blue Lake monimolimnion vs. its microbial plate (p = 0.38), and Mahoney Lake microbial plate vs. Lime Blue Lake monimolimnion (p = 0.063).

The ten genes with the highest expression levels (Fig. 6) in each sample group also contributed to the largest average dissimilarity between clusters, i.e., major head proteins differentiated the Lime Blue Lake microbial plate from the other two lakes (SIMPER, 22.9% of the dissimilarity with Poison Lake and 17.7 with Mahoney Lake, p = 0.002 and 0.043, respectively). Although contributing less to the dissimilarity, expression of an 8.4 kDa Cro protein was significantly higher in the microbial plate of Lime Blue Lake than in Mahoney Lake (SIMPER, p = 0.027, average dissimilarity = 0.071%) and Poison Lake (SIMPER, p = 0.013, average dissimilarity = 0.080%) microbial plates. This may indicate that there are more temperate viruses in the microbial plate of Lime Blue Lake undergoing the lytic cycle, as the expression of Cro is responsible for inducing prophages in the model phage lambda [48].

Discussion

Dense microbial plates display low virus-to-microbe ratios

The two highly stratified, PSB-dominated lakes studied here, Mahoney and Poison Lakes, exhibited lower VMR at and below the microbial plate compared to the mixolimnion (Fig. 2). This was also observed in a study investigating VMR using flow cytometry in Lake Cadagno, a highly stratified meromictic lake in Switzerland [14]. These results contrast with the pattern observed in Lime Blue Lake, where VMR was higher at the microbial plate and monimolimnion compared to the mixed layer. Though difficult to fully decouple from other environmental factors, this pattern could be partially attributed to variations in cell density among these lakes. Lime Blue Lake had the lowest cell densities among the three study lakes, and the cell densities did not differ significantly across the water column (Fig. 2). In contrast, the cell densities increased by one order of magnitude between the surface and microbial plate of Mahoney and Poison Lakes. The decrease in VMR in these lakes indicates a decrease in the frequency of lytic viral infections. One possible mechanism explaining the reduced lysis is the increase in lysogenic phage infections by temperate phages, causing an overall decrease in VMR at the community level [49,50,51]. During lysogeny, the viral genome integrates into the host genome and is replicated alongside the host as it grows, as opposed to the viral particle replication and cell death caused by lytic infections. Temperate phages have been detected in meromictic lakes, and in one case, GSB populations were associated with two prophages consistently infecting over 10 years [52]. The Piggyback-the-Winner framework has proposed that the prevalence of lysogeny increases in high-density microbial communities, driven by an increase in phage coinfections in these conditions of high encounter rates [49,50,51, 53]. Higher coinfections increase the likelihood of lysogeny by increasing the concentration of repressors of the lytic cycle and activating integrase expression [51, 54]. This explanation would fit the observations of the present study, where higher densities were observed in Mahoney Lake and Poison Lake, concurrent with the lowest VMRs. A non-mutually exclusive hypothesis to explain the decrease in VMR is that defense mechanisms such as CRISPR systems arise in these dense microbial plates. Previous research has shown that CRISPR systems are prevalent in PSB [55], which dominate in both Poison and Mahoney Lakes [6, 8, 13]. These defense systems could contribute to a reduction in detectable viral activity by inhibiting viruses’ lytic cycles. CRISPR-Cas systems are abundant in environments with high viral abundance and low viral diversity [56], which is the case in the microbial plates of this study (Fig. S3).

Four lines of evidence indicate that both CRISPR defense and lysogeny contribute to the difference in observed VMR and viral activity. (1) High transcription levels of lysogenic and CRISPR-associated genes in the microbial plates compared to the mixolimnion (Fig. 4) and higher in Poison Lake and Mahoney Lake relative to Lime Blue Lake. (2) Higher expression of the prophage-inducing cro in Lime Blue Lake’s microbial plate compared to the other two lake’s microbial plates. (3) High expression of structural genes as a driver of the dissimilarity between Lime Blue Lake samples compared to other lakes (Fig. 6). (4) The negative relationship between specific conductivity, a proxy for salinity, and viral activity (Fig. S4). Salinity has been positively associated with lysogeny in high-altitude lakes [57] and estuaries [58]. Although the mechanism behind this remains to be explored, salinity can impact phage viability and particle stability [59, 60], which could lead to a preferential lifestyle of lysogeny. This agrees with our finding, where the two saline lakes, Mahoney and Poison Lake, also have the highest estimated frequency of lysogeny.

Transcribed viral genes show evidence for high viral activity in Lime Blue Lake and lower activity in PSB-dominated microbial plates

The most transcriptionally active viral contigs encoded genes related to viral particle production during the lytic cycle, specifically major capsid and tail proteins. Lime Blue Lake’s microbial plate sample had the highest viral transcription levels (Fig. 6), with SIMPER analysis identifying the major head protein as a significant contributor to the dissimilarity between Lime Blue Lake and the other lakes, which corroborates the inference of higher lytic activity inferred from high VMR (Fig. 2c) and the proportion of active viruses (Fig. 3c). Conversely, the viral transcripts with the highest expression levels in Poison Lake’s microbial plate and monimolimnion were associated with host metabolism, followed by structural genes. This result indicates that these viruses may be primarily influencing the metabolic processes of the host organisms through either lytic or lysogenic infections. The ability of phages to manipulate host metabolism has been observed in marine environments [18, 19, 44, 61]. The most active virus in the microbial plate of Poison Lake encoded a gene with significant similarity to the D1 protein of the Photosystem II in Cyanobacteria (Fig. 6). The L and M subunits of the PSB photosystem reaction center belong to the type II reaction center proteins and are homologs of the D proteins in the reaction centers of Cyanobacteria [62]. Therefore, the highly expressed viral gene observed here in the PSB-dominated Poison Lake likely encodes for the L subunit of the PSB reaction center. The gene observed here has only 63.73% amino acid sequence identity with the closest relative Photosystem II D1 protein from NCBI’s Refseq protein database (WP_224340432.1), indicating significant divergence. The expression of the L protein of the PSB reaction center observed here is likely analogous to the expression of phage-encoded psbA genes for the D1 protein during infection of Cyanobacteria. This expression keeps the reaction centers running during infection, powering the light reactions of photosynthesis and providing energy in the form of ATP for viral particle production [63,64,65].

The Poison Lake microbial plate displayed a higher concentration of hydroxy-Bchl a, a product of Bchl a oxidation, compared to the microbial plate in Mahoney Lake (Table 1), although the microbial plates from both lakes are dominated by PSB. This observation coincides with the high expression of viral-encoded metabolic genes, including photosystem genes, in Poison Lake compared to Mahoney Lake, as discussed above (Fig. 6), and leads us to speculate where viral gene expression is associated with the change in pigment profile. Cyanophages have been shown to affect pigment production. For instance, they can downregulate host genes involved in chlorophyll biosynthesis and other light-harvesting molecules [66]. Inversely, cyanophages can catalyze the production of phycocyanobilin by skipping an intermediary step in the biosynthesis pathway by encoding genes such as pebS [67]. A similar mechanism using a different gene has also been observed in phages infecting Alphaproteobacteria [68]. Changes in pigment profiles have also been reported in Serratia marcescens, where lytic infection increases the pigment prodigiosin [69]. A eukaryotic virus infecting the marine prymnesiophyte Emiliania huxleyi causes structural changes to the photosynthetic apparatus and increased hydroxylation of Chl a, culminating in the decreased quantum efficiency of Photosystem II and reduced chlorophyll fluorescence [70, 71]. An analogous process could explain the higher concentration of hydroxy-Bchl a in the Poison Lake microbial plate observed here. If that process is occurring, light and dark reactions of photosynthesis could be decoupled, resulting in a mismatch between the sulfur oxidation and carbon fixation rates, ultimately observed as mismatched sulfur and carbon isotopic signatures of PSB [13]. Because the light reactions of photosynthesis in PSB are also dependent on carotenoid pigments such as okenone, an increase in the light reactions of photosynthesis could also indirectly increase the demand for okenone synthesis, further decoupling the organic and inorganic geochemical signatures of PSB.

The metabolic effect of viruses on their hosts would have important implications for interpreting geological evidence used to infer the biogeochemical dynamics of ancient oceans. The stratified meromictic lakes examined in this study are modern analogs for ancient euxinic oceans. Viral activity could be decisive in shaping microbial metabolisms, particularly through lysogenic infection and the expression of auxiliary metabolic genes (AMGs), and significantly influence the production and preservation of metabolic outputs, which are used as biosignatures to reconstruct these ecosystems in the geological record. The observed decoupling in biosignatures may help explain inconsistencies in the different biosignatures of PSB in the sediments of meromictic lakes [17, 18] and improve the interpretation of the geologic record of anoxic water columns. Physiological studies quantifying the coupled rates of sulfur oxidation, carbon fixation, and pigment production in infected and uninfected PSB are necessary to test this hypothesis.

Viral-encoded metabolic genes were absent in the ten most abundant viral contigs in Mahoney Lake, which may be due to the difficulty in assembling long viral genome sequences solely through metatranscriptomic data. Highly transcribed AMGs from integrated prophages are more likely to be missed with our method, especially if they have a high sequence identity with the host version of the gene, leading to misidentification as host (not viral) gene transcription. Viruses in these lakes contain genes that could disrupt key metabolic pathways such as sulfur oxidation, carbon fixation, and pigment production [13]. Because these biosignatures are used to interpret the geological record, incorporating the role of viruses into models of ancient ocean biogeochemistry is crucial for refining our understanding of Earth’s redox evolution.

Conclusions

Here, we show that the viral life strategies and activity are strongly stratified in lakes dominated by photosynthetic sulfur-oxidizing bacteria. In all lakes, the viral community composition differed across depths, and viral diversity decreased in anoxic layers compared to the mixolimnion. The highly stratified lakes with high densities of PSB, Mahoney and Poison Lakes, exhibited lower VMRs at the microbial plate, due to decreased lytic replication, likely caused by lysogeny and defense systems. In contrast, Lime Blue Lake, which displays a lower-density microbial plate with a mix of GSB and PSB, showed higher VMR at the microbial plate and the highest viral transcriptional activity, including genes for viral particle production, and indicating higher lytic activity than in the microbial plates of the other two lakes. Transcriptional analysis of viral genes indicates active viral influence on host metabolism, particularly in Poison Lake, although the detection of AMGs remains challenging. These findings underscore the need for further metagenomic approaches to fully understand viral dynamics in stratified lakes and potential impacts on early Earth oceans, highlighting the complex interactions between viruses, their hosts, and environmental factors.