Sampling was carried out on board the R/V Aegaeo in April 2016, at the onset of seasonal thermal stratification. Seawater was collected between 09:00 and 12:00 am by Rosette-Niskin deployment from 4 stations in the Western Levantine basin (Eastern Mediterranean Sea), encoded as LV3 (35.0333oN, 23.4667oE), LV10 (34.6667oN, 24.3667oE), LV13 (34.2500oN, 25.4833oE), and LV18 (34.4333oN, 26.3833oE). Four depth layers were sampled in each station: 5, 50, 75, and 1,000 m. The collected volume was transferred from the Niskin bottles into acid-cleaned and deionized water-rinsed low-density polyethylene containers, and further processed for the various analyses.

Hydrographic profiles of salinity (S), temperature (T) and density (D) were obtained using a Seabird SBE 9plus CTD underwater unit connected to a Seabird SBE 11 deck unit (Velaoras et al., 2019). The methodologies for the determination of dissolved oxygen (DO), nitrate (NO₃⁻), nitrite (NO₂⁻), silicate (SiO₄), ammonium (NH₄), phosphate (PO₄³⁻), dissolved organic carbon (DOC), phosphorus (DOP) and nitrogen (DON), total nitrogen (TN) and phosphorus (TP) and chlorophyll a (Chl) concentrations are presented in the Supplementary material.

The abundance of virus-like particles (VLP) was determined based on the protocol of Brussaard (2004) and that of heterotrophic and autotrophic bacteria based on Marie et al. (1997) from glutaraldehyde-fixed samples (0.5% final concentration). Samples were stained with SYBR™ Green I nucleic acid stain (ThermoFisher Scientific) at 5 × 10⁻⁵ and 4 × 10⁻⁴ final dilution of the stock solution for VLP and heterotrophic bacteria, respectively, and then incubated for approx. 15 min at 80°C and for approx. 30 min in the dark for VLP and heterotrophic bacteria, respectively. Autotrophic bacteria were distinguished based on their auto-fluorescence signals. Yellow-green latex beads of 1 μm nominal size (Polysciences) were added and used as an internal standard of fluorescence. A FACSCalibur™ instrument (Becton Dickinson) was used at conventional air pressure, with an air-cooled laser at 488 nm and standard filter setup. Data were processed with the CellQuest™ Pro software (Becton Dickinson).

Viral production rates were estimated following the viral reduction approach (Winter et al., 2004; Winget et al., 2005). Two L of seawater were pre-filtered through 0.8 μm pore size polycarbonate membranes to exclude large particles and grazers. The 0.8 μm-filtrate was concentrated using a tangential flow filtration system (Sartorius, VF20P7, 0.2 μm MWCO). The 0.22-μm filtrate was back-flushed with the use of the reverse flow mode of a peristaltic pump (Masterflex, EW-07523-80) in order to produce an ultra-concentrated-in-bacteria seawater volume (100 mL). The remaining 0.22-μm filtrate was further filtered with a spiral-wound ultra-filtration cartridge (Sartorius, VF20P2, 30.000 MWCO) so to produce a virus-free ultra-filtered volume (200 mL). Bacterial-concentrated and virus-free water masses were mixed gently, and then equally distributed into 50-mL centrifuge tubes (triplicate incubations). Tubes were incubated at in situ temperature (as determined by CTD) and in the dark for 24 h. Every 0, 1, 3, 6, 12, and 24 h during the incubation, samples were fixed for the determination of viral and bacterial abundances as described above. Lytic viral production (lytic VP, expressed in virus-like particles mL⁻¹ h⁻¹) was estimated from the slopes of the relationships between net increases in viral abundance over the respective time period of the net increase.

Lytic VP was further corrected for the bacterial losses due to the filtration; the recovery percentage was determined by dividing the mean bacterial abundance at the onset of the incubations (B0) by the in situ bacterial abundance (B original), and this factor was multiplied to the lytic VP measurement.

The methodology for the assessment of heterotrophic bacterial production is presented in the Supplementary material.

Seawater was filtered through 0.2 μm polycarbonate membranes under low vacuum. The 0.2-μm filtrate was chemically treated with 1 mg L⁻¹ FeCl₃ to achieve viral particle flocculation within the following 6–12 h (John et al., 2011) and filtered again through 1 μm polycarbonate filters. The 1 μm filters were stored at 4°C, pending resuspension in ascorbic acid buffer. The solution of ascorbate-EDTA buffer was prepared daily (0.25 M ascorbic acid, 0.2 M Mg₂EDTA, pH 6–6.5), kept in the dark and added in the viral flocculate, followed by hand shaking and overnight rotation at 4°C. After resuspension, viral particles in liquid were retained from the 1 μm filter by ultra-centrifugation at 141.000 g (SorvallTM WX100 ultracentrifuge, ThermoFisher Scientific, Sorvall TH 641 swing out rotor). Viral DNA was extracted following a CTAB protocol (Winnepenninckx et al., 1993) and purified as described in details elsewhere (Tsiola et al., 2020). The extracted DNA was dissolved in ultrapure water and stored at −20°C. Quantification of viral dsDNA was done with the Qubit high sensitivity assay kit in a 3.0 Qubit™ fluorometer (ThermoFisher Scientific).

Viral DNA shearing was done at 300 bp using the standard protocol for Covaris™ focused ultra-sonicator system. An indexed library for Illumina sequencing was prepared using the NEBNext Ultra DNA Library Prep Kit for Illumina (New England BioLabs) following the manual instructions. Size selection was done using AMPure XP beads (Beckman Coulter). PCR cycles were 6 (according to the manufacturer advice with regards to the amount of DNA input). Metagenomic libraries were sequenced in the Illumina Hiseq 4000 platform available at KAUST Bioscience Core Lab using paired-end sequencing.

Viral metagenome reads in FASTQ format were imported to CLC Genomics Workbench v.7 (CLC Bio) and trimmed using a minimum phred score of 20, a minimum length of 50 bp, allowing no ambiguous nucleotides and trimming off Illumina sequencing adaptors if found. The trimmed metagenome reads were assembled using CLC’s de novo assembly algorithm, using a k-mer of 63 and a minimum scaffold length of 500 bp. Raw reads were deposited on the NCBI with reference PRJNA996089.¹ The assembled contigs were then analyzed using the iVirus pipeline (Bolduc et al., 2017) through the Cyverse platform (Goff et al., 2011) and VirFinder software (Ren et al., 2017) with criteria as described previously. Viral sequences with a VirFinder score ≥ 0.7 and p < 0.05 and VirSorter categories 1 and 2 were used for further analysis. The abundance and metabolic potential of auxiliary and metabolic potential of AMGs was determined through VIBRANT (v1.2.1; Kieft et al., 2019) using the virome option. The taxonomy of viral contigs was determined with the VPF-Class software which uses the IMG/VR v3 database for this purpose (Pons et al., 2021; Roux et al., 2021). Further, viral metagenomic sequence data were processed with the MetaPop multi-functional bioinformatic pipeline (Gregory et al., 2022) for macrodiversity analyses with the “vegan” R package. MetaPop was run with default parameters and cut-offs. Raw population abundances were normalized in order to mitigate sample-to-sample variation in the number of reads, via normalizing to the library with the highest number of reads. The abundance table of viral operational taxonomic units and then, it was used to generate a Bray-Curtis distance matrix. MetaPop macrodiversity outputs included population abundances, and alpha- (within community) and beta-diversity (between community) indices. Realizing the need to apply up-to-date pipelines for quality assessment of our viral metagenomes, we tested the clustering function of CheckV (Nayfach et al., 2021). However, no viruses were clustered indicating a lack of duplicates in our analyses. The raw number of viral metagenome reads, the quality of reads, the number of total contigs, the contig length (total, minimum and maximum) and the N50 contig length are summarized in Table 1.

Principal coordinates analysis was applied to coordinate the viral diversity and AMG data (Clarke and Ainsworth, 1993). Canonical analysis of principal coordinates (CAP) was applied to define clusters of samples based on “depth,” using the set of physical and chemical measurements (S, T, D, and the concentrations of PO₄³⁻, NO₃⁻, NO₂⁻, NH₄, SiO₄, TN, TP, DO, DOC, DOP, DON). The clustering was tested by permutational multivariate analysis of variance (PERMANOVA). Then, patterns in viral diversity and AMG content between the stations were tested for significance by applying one-factor PERMANOVA using the factor “depth.” The null hypothesis was that there are no differences. When 5, 50, and 75 m were not different to each other, they were named as “surface/subsurface samples” for brevity. Bray–Curtis dissimilarity matrices on square-root transformed metagenomics data were constructed (Clarke and Warwick, 1994) to avoid misleading interpretations due to the sample-to-sample variation in read counts (Table 1). Hypothesis testing was performed using 999 permutations and pairwise tests using a significant level of 0.05. A principal coordinate analysis (PCoA) plot of all Bray-Curtis distances was created via the MetaPop pipeline and is presented as a means of macrodiversity visualization of the viral metagenomic data. The list of explanatory variables was normalized and the normalized matrix was used to create a resemblance matrix using Euclidean distances. Statistical analyses were done with the software package PRIMER v6 (PRIMER-E Ltd, Plymouth Marine Laboratory, Natural Environmental Research Council, United Kingdom) with PERMANOVA + add-on software (Anderson et al., 2008; Somerfield, 2008).

One-way analysis of variance (ANOVA) was applied to check for significant differences among the different depth layers in: the percentage contribution of individual genera and families over total viral contigs, the percentage contribution of viral AMGs over total number of AMGs and over their assigned KEGG category, the abundances of bacteria and virus-like particles, and finally in the physical, chemical and biological variables. The significance of the differences was assessed with post hoc Tukey test. Homogeneity of variance was checked using Levene’s test. ANOVAs were performed using IBM SPSS statistics software v23.

None of the measured variables differed between stations at the horizontal scale (PERMANOVA), thus, stations LV3, LV10, LV13, and LV18 are considered as replicates. Supplementary Table 1 summarizes output of all PERMANOVAs testing for differences between the sampling stations. Output of all PERMANOVAs testing variation among sampling depths are presented in Table 2, and not in the main body text (degrees of freedom, mean square, Pseudo-F ratio, p value). The Canonical Analysis of Principal coordinates (CAP) of physical and chemical variables (S, T, D, and the concentrations of NOx, PO₄³⁻, DO, SiO₄, DO, DOC, DOP, DON, TN and TP) revealed four significantly different groups, which were confirmed by PERMANOVA (p < 0.001). Samples from the same depth layers were grouped together and distantly from the others (Figure 1A).

S averaged 38.98 ± 0.13 in the sampling stations with equally high values at surface/subsurface depths (one-way ANOVA, p > 0.05). Approx. 0.21 lower S was measured at 1,000 m compared to averaged surface/subsurface values (Table 3) but the difference was significant only when 1,000 m were compared to 75 m (one-way ANOVA, p < 0.05). T ranged between 13.80°C and 18.30°C (Table 3). Small variation was seen between 50 and 75 m (16.81°C ± 0.46°C). All other depth layers exhibited significantly different T between them (one-way ANOVAs, p < 0.05), with the coldest waters seen at 1,000 m (post hoc Tukey test, p < 0.05).

The concentrations of dissolved inorganic and organic nutrients significantly differed between surface/subsurface depths and 1,000 m (one-way ANOVAs, p < 0.05), being higher at 1,000 m (post hoc Tukey test, p < 0.05). DO concentration followed the opposite pattern (i.e., it was significantly lower at 1,000 m compared to surface/subsurface layers, post hoc Tukey test, p < 0.05). Measurements of these variables as well as of D, DOC, DON, DOP are presented in Table 3, and are discussed in the Supplementary material. Chl concentrations differed significantly between 5 and 75 m, and between 50 and 75 m (one-way ANOVAs, p < 0.05). Chl concentration was ~4x higher at 75 m than 5 m, and ~3x higher at 75 m than 50 m (post hoc Tukey tests, p < 0.05, Table 3).

Neither viral (4.6 × 10⁶ ± 4.3 × 10⁵ mL⁻¹) nor heterotrophic bacterial (3.2 × 10⁵ ± 8.4 × 10⁴ mL⁻¹) abundances differed within the surface/subsurface layers (one-way ANOVAs, p > 0.05, Figures 2A,B). Significant differences were seen between surface/subsurface and 1,000 m samples at both cases (one-way ANOVA, p < 0.05); abundances were significantly lower at 1,000 m than at the other depths (post hoc Tukey test, p < 0.05). The abundance of Synechococcus did not differ between the sampling depths, opposite to Prochlorococcus (one-way ANOVA, p < 0.05), which was significantly lower at 5 than 75 m, as well as at 50 than 75 m (post hoc Tukey tests, p < 0.05). Synechococcus overcame in abundance Prochlorococcus by approx. one order of magnitude at 5 m, and the opposite was noticed at 75 m (Figures 2C,D). The resulting virus-to-bacteria ratio (VBR) varied greatly between surface/subsurface (9–21) and 1,000 m (56–107) considering in the calculation total bacterial abundances (Figure 2E). Significantly higher VBR was measured at 1,000 m than surface/subsurface samples (post hoc Tukey test, p < 0.05).

Lytic VP did not differ between the sampling depths (one-way ANOVA, p > 0.05), and it was slightly higher at 5 m (1.9 × 10⁶ ± 1.3 × 10⁶ cells mL⁻¹ h⁻¹) than the other depths (8.3 × 10⁵ ± 3.6 × 10⁵ cells mL⁻¹ h⁻¹; Figure 2F). BP did not differ within the surface/subsurface samples (one-way ANOVA, p > 0.05) and ranged between 2.35 and 9.28 ng C L⁻¹ h⁻¹ (Figure 2G). There was a significant difference between surface/subsurface and 1,000 m samples (one-way ANOVA, p < 0.05); significantly lower BP levels were seen at 1,000 m than surface/subsurface (post hoc Tukey test, p < 0.05).

Viral community composition differed with the sampling depth at the family and genus level (PERMANOVAs, p < 0.01). Bray-Curtis distances were plotted using the principal coordinates analysis method (PCoA, Figure 1B). The two 1,000-m samples were more dissimilar to the rest and exhibited the highest species richness among all samples, while surface/subsurface samples were less dissimilar to each other and clustered together (Figure 1B). Alpha diversity indices (richness, Chao1, ACE, Shannon’s H, Simpsons, inverse Simpsons, Fisher, and Pielou’s J) for the samples are presented in Supplementary Table 1A. A percentage of 30 ± 7% of the contigs remained unassigned to the family level. Within the assigned contigs, the dominant families were Podoviridae (29% ± 3%) and Siphoviridae (23% ± 3%, Figure 3A). At 1,000 m, significantly lower contribution of Podoviridae (24% ± 2%) was seen compared to 5 m (27% ± 1%) and 50 m (30% ± 2%; one-way ANOVAs, p < 0.05). The same was noticed for Siphoviridae that contributed less at 1,000 m (19% ± 4%) compared to 5 m (22% ± 2%) and 50 m (24% ± 2%; one-way ANOVAs, p < 0.05). The third most abundant family was that of Myoviridae, exhibiting a decreasing trend (post hoc Tukey test, p < 0.05) from surface (23% ± 5% relative abundance) to deep waters (14% ± 0.3%, Figure 3A). Phycodnaviridae family was seen at all surface/subsurface stations with no significant differences in its contribution (0.7%–1.4%). The remaining taxa (<0.5%) were rare families (27 in total). Most of them were found in few stations each. The families Adenoviridae, Baculoviridae, Herpesviridae, Inoviridae, Iridoviridae, Marseilleviridae, Microviridae, Mimiviridae, Poxviridae were found in all stations (Figure 3B). The contribution of genera is presented in the Supplementary material (Supplementary text and Supplementary Table 2).

Information about the potential host assignment is presented in the Supplementary material.

A list of 143 AMGs was counted and associated with a metabolic pathway as defined by the Kyoto Encyclopedia of Gene (KEGG) using the VIBRANT method. The general pathways are shown in Table 4 and include: “amino acid metabolism,” “carbohydrate metabolism,” “metabolism of cofactors and vitamins,” “energy metabolism,” “lipid metabolism,” “glycan biosynthesis and metabolism,” “nucleotide metabolism,” “biosynthesis of secondary metabolites,” “folding, sorting and degradation,” “metabolism of other amino acids,” “metabolism of terpenoids and polyketides,” and “xenobiotics biodegradation and metabolism.” The specific metabolic pathways within the general ones are shown in Figure 4 and Supplementary Table 3. Statistical analysis for differences between the sampling depths was performed (1) in percentage contribution of specific AMG pathways within their general metabolic category, as well as (2) in percentage contribution of specific AMG pathways over total AMG reads.

The content of AMGs differed with the sampling depth when considering both general and specific metabolic pathways (PERMANOVAs, p < 0.01). AMG content differed between 1,000 m and surface/subsurface samples. In terms of percentage contribution, the most frequently-detected general metabolic pathway was “amino acid metabolism” (34% ± 6%), followed by “carbohydrate metabolism” (18% ± 3%), “metabolism of cofactors and vitamins” (14% ± 4%), “energy metabolism” (12% ± 5%), and then “lipid metabolism” and “glycan biosynthesis and metabolism” (Table 4). The remaining categories are presented in the Supplementary material (Supplementary text and Supplementary Table 3).

We studied oceanic viral populations taking into consideration the lytic efficiency and the possible mechanisms that viruses employ to alter their hosts’ metabolism using their own metagenome footprint. Our results demonstrate clear depth-dependent patterns (surface/subsurface shallow vs. deep waters) in viral community composition and AMG content during the early stratification season in the Eastern Mediterranean Sea. The vertical environmental variability of this oligotrophic basin was considered in order to uncover, and propose viral actions that lead to energy acquisition and lysis sustenance. The depth-dependent patterns were linked to differences in temperature, nutrient availability (total and dissolved nitrogen and phosphorus), host productivity and host population density. While viral productivity levels did not change significantly along the water column, the numbers of free viruses per bacteria (VBR) were significantly higher at 1,000 m, resembling a similar peak in VBR in shallower mesopelagic depths in the same study area (Magiopoulos and Pitta, 2012). Interestingly, the diversity of viral taxa and AMGs was also substantially different at this depth. The presented findings suggest that a similar quantitative effect of viruses on their hosts (lytic activity) is achieved by different viral populations carrying different auxiliary metabolic genes.

Differences in viral communities with temperature were highlighted in the Mediterranean and globally (Coutinho et al., 2019; Gregory et al., 2019) and differences with depth were also reported for the hadal zone (Jian et al., 2021) and the Pacific Ocean using either high-throughput sequencing technology (Hurwitz et al., 2015) or genome fingerprinting approaches (Steward et al., 2000; Brum, 2005). Coutinho et al. (2019) have reported a peak in the Shannon diversity index in deep water of the western-Mediterranean. Lower contribution of the dominant families in deep rather than surface waters was seen in the hadal zone of Challenger Deep (Gao et al., 2022). Similarly, the highest species richness was recorded at 1,000 m, with the community encompassing less Podoviridae and Siphoviridae members but more unassigned reads, which is reasonable as a more complex community may sustain and optimize ecosystem functioning under low-oxygen and low-nutrients conditions.

Additionally, at 1,000 m we found that the AMG content was different than at the other depths. The main distinctive features of the deep layers in comparison to the shallower ones were (a) a higher contribution of genes involved in the sulfur relay system, the synthesis and metabolism of glycan and the metabolism of fructose, mannose, fatty acids, cofactors and vitamins and secondary metabolites, (b) a lower contribution of genes involved in energy metabolic pathways, and the metabolism of arginine, proline and galactose, (c) different AMG content within the categories of metabolism of amino sugar and sugar nucleotide, fructose and mannose, lipopolysaccharide, sulfur and methane, lipids and purine, and finally (d) the absence of genes related to pentose phosphate pathway, and the unique finding of genes related to pyruvate and biotin metabolism, and prodigiosin biosynthesis. At this point, we urge the readers to take into consideration the low number of replicated samples when deepening to our findings. Consideration of seasonal and inter-annual changes in irradiance, temperature and nutrient levels with depth would have also improved our exploration, since it is widely known that these changes largely impact viral activity (Brum et al., 2016; Coutinho et al., 2017; Gainer et al., 2017; Puxty et al., 2018).

In the energy-limited bathypelagic zone, it is expected that viruses, as their hosts (Konstantinidis et al., 2009), need to have higher metabolic flexibility and thus, exploit a variety of AMGs to support their reproduction in comparison to the more limited resources. Recent reports consider that viral genomes are under the selective pressure of the environment (Heyerhoff et al., 2022); in natural populations, the availability of nitrogen (Jian et al., 2021) and phosphorus (Kelly et al., 2013) and the variable temperature conditions (Coutinho et al., 2019) seem to drive differences in viral genomes from the community to the protein level. While the total amount of AMGs entailed in energy acquisition was lower at 1,000 m compared to surface/subsurface waters, the particular genes involved were rather more variable and seemed to ensure that diverse substrates and energy sources can be used in the EMS deep environment where light, organic matter and nutrients are limited and the growth of bacteria is carbon-limited (Rahav et al., 2019). In the bathypelagic realm of the oligotrophic EMS, nutrients are scattered in time and space, thus their high demand may have caused relevant viral genome adaptions that are discussed in chapters 4.2. and 4.3.

A higher contribution of AMGs related to energy production via gluconeogesis (pps, ppsA) was seen in the deepest samples, suggesting that the use of alternative carbon sources for the synthesis of glucose and polysaccharides is positively selected at 1,000 m. In deep waters, hexoses are preferentially removed from high-molecular weight DOC by bacteria (Garel et al., 2021). The number of hexose metabolism genes was higher in the deep, and this fact led us to think that viruses may provide a vital mechanism to build the precursors for ATP synthesis and aid viral proliferation, thus supporting the high microbial respiration needs in the bathypelagic zone (Acinas et al., 2021). AMGs involved in the pentose phosphate pathway (PPP) was proposed to divert host metabolism to nucleotide synthesis and energy production (Thompson et al., 2011). PPP is probably not boosted in the bathypelagic zone in our study area but this hypothesis needs confirmation using a larger set of samples and replicates. A range of PPP genes (rpi, gnd, zwf) was seen in the surface/subsurface samples (but not in the deep) where light and energy availability is not a limiting step for host and viral productivity, similar to recent observations for surface coastal waters of the Mediterranean (Tsiola et al., 2020).

A higher contribution of AMGs related to fatty-acid biosynthesis was seen at 1,000 m, suggesting that the elongation of fatty acids is promoted through viral AMGs. fabG encodes for the NADPH-dependent reduction of beta-ketoacyl-ACP substrates to beta-hydroxyacyl-ACP products, and has been found in viromes in the hadal zone (Jian et al., 2021). Alongside, the reaction of the enzyme of fabZ was a bottleneck in the fatty-acid production in numerous microbial hosts (Jeon et al., 2012). The presence of fabG and fabZ in the deep may be related to the need of fatty acids as electron donors in several bacteria. Sulfate-reducing bacteria for example are optimized to fully oxidize fatty acids; interestingly, their presence was indicated by the concomitant detection of CysH gene at 1,000 m.

The higher contribution of AMGs involved in the glycan biosynthesis in the deep supports the scenario of accelerated viral folding and structuring in these systems, something that totally agrees with the high VBR and VP lytic measurements in our study. The identification of genes involved in the amino sugar and nucleotide sugar metabolism along with genes involved in glycoprotein synthesis and metabolism (waaF, rfaE; now known as hldE, having key roles in the biosynthetic process of lipopolysaccharides, Graves et al., 2001; Weynberg et al., 2017) is a unique coupled finding for deep waters, that confirms the magnitude of hosts-gene hijacking in natural populations. In addition, the different repertoire of AMGs at the four depths involved in glycan metabolism, thus in correct host recognition, attachment and entry to the bacterial host cell (Kortright et al., 2020). indicates that attachment and penetration to hosts may happen by exploiting different strategies in deep vs. shallow waters, ultimately leading to equally high lytic activity at all depth layers.

The higher contribution of the mec gene at 1,000 m seems to assist nutrient acquisition in this deep ecosystem, since the gene encodes for a sulfur carrier protein that aids the biosynthesis of cysteine (Hesketh-Best et al., 2022) and possibly the degradation of sulfated organic matter (Han et al., 2023), thus viral particle formation too. Despite the low understanding of the role of mec in the viral genomes, it seems that recent investigations on viral implications on the oceanic sulfur cycling are reinforced and extended by our findings (Kieft et al., 2021); viral genes related to sulfur relay were found at 1,000 m, thus possibly contributing to sulfur’s budget manipulation.

We propose that the various above-mentioned mechanisms may have been adapted by the viral populations in order to more efficiently use the scattered energy sources in the oligotrophic deep part of the EMS water column. In such case, the selective advantage of viruses can reasonably justify the high VBR and lytic production in all depth layers, despite the low nutritional status of the bathypelagic zone.

Arginine and proline metabolism AMGs contributed less in the deepest samples, confirming that viruses “select” for certain amino acids for tolerance to cold temperatures (Alarcón-Schumacher et al., 2021) and optimal capsid formation (Weynberg et al., 2017).

The most noticeable difference in the AMG content between between 5, 50, and 75 m was that of the DCM; at 75 m, photosynthesis-related psbA and psbD outnumbered the energy-metabolism AMGs (>90%), while at 5 and 50 m these genes were significantly less. Indeed, comparing with in situ primary production measurements (Livanou et al., 2019; Psarra unpubl. data), primary production maxima at subsurface layers were consistently recorded just above the DCM layer. It was found that the psb genes are expressed during the viral latent period, thus their overrepresentation in the DCM seems supportive of extensive photosynthesis at this layer (Sharon et al., 2007; Sieradzki et al., 2019). At 75 m, no energy-production related to sulfur seemed to occur through viral AMGs, and photosynthesis governed. On the contrary, at 5 and 50 m, a repertoire of viral AMGs related to sulfur exploitation either for ATP production or for sulfur assimilation and biosynthetic processes was found; dissimilatory (e.g., dsrA) and assimilatory (e.g., cysC) sulfur metabolism genes support the auxiliary metabolism of sulfur in the host cells (Kieft et al., 2021). The dsrA, encoding for dissimilatory sulfur reductase within the pathway of sulfur oxidation was noted in the latter study, while Roux et al. (2016) and references therein have mentioned only the dissimilatory sulfur reductase subunit C (dsrC) in various pelagic and benthic samples. The authors proposed that viruses assist their sulfur oxidizer hosts through the dsrC-like genes under nutrient limiting conditions (Roux et al., 2016) that could be valid in the oligotrophic conditions of the studied area. Our findings are in agreement with previous observations that viral diversity does not differ substantially between surface and deep-chlorophyll maximum (DCM) layers except when considering protein clusters (Brum et al., 2015). An additional indication that excessive energy at DCM was diverted toward viral replication is the peak of purine metabolism genes (Puxty et al., 2015) at 75 m.

A comprehensive description of viral dynamics in surface (5 m), subsurface (50 and 75 m) and deep (1,000 m) off-shore waters of the EMS was achieved by combining whole viral metagenome sequencing, flow cytometric and productivity analyses, assuming that such a combination is the foundation of meaningful viral ecology discussions. The findings confirm our hypothesis that a more diverse repertoire of AMGs is present in bathypelagic than shallow EMS waters, diverting hosts’ metabolism toward energy production and thus, viral replication too. Further work is needed to generalize our implications to other ecosystems and seasons with regards to the diel, seasonal and annual variations in irradiance, temperature and other physicochemical characteristics of the water masses, since the sampling power in this study was limited by field-work constraints and analytical costs. AMG content was substantially different between 1,000 m and the shallower layers, revealing novel implications for the energy-production mechanisms in the resource-depleted deep waters of the EMS. The lysis rate was, however, similar along the vertical profile. Altogether, our findings advocate that viral interference with hosts’ metabolism is affected by the particular environmental setting, and must be carefully considered for deciphering the biogeochemical turnover in the oceans.