Depth-driven patterns in lytic viral diversity, auxiliary metabolic gene content, and productivity in offshore oligotrophic waters

Introduction Marine viruses regulate microbial population dynamics and biogeochemical cycling in the oceans. The ability of viruses to manipulate hosts’ metabolism through the expression of viral auxiliary metabolic genes (AMGs) was recently highlighted, having important implications in energy production and flow in various aquatic environments. Up to now, the presence and diversity of viral AMGs is studied using -omics data, and rarely using quantitative measures of viral activity alongside. Methods In the present study, four depth layers (5, 50, 75, and 1,000 m) with discrete hydrographic features were sampled in the Eastern Mediterranean Sea; we studied lytic viral community composition and AMG content through metagenomics, and lytic production rates through the viral reduction approach in the ultra-oligotrophic Levantine basin where knowledge regarding viral actions is rather limited. Results and Discussion Our results demonstrate depth-dependent patterns in viral diversity and AMG content, related to differences in temperature, nutrients availability, and host bacterial productivity and abundance. Although lytic viral production rates were similar along the water column, the virus-to-bacteria ratio was higher and the particular set of AMGs was more diverse in the bathypelagic (1,000 m) than the shallow epipelagic (5, 50, and 75 m) layers, revealing that the quantitative effect of viruses on their hosts may be the same along the water column through the intervention of different AMGs. In the resource- and energy-limited bathypelagic waters of the Eastern Mediterranean, the detected AMGs could divert hosts’ metabolism toward energy production, through a boost in gluconeogenesis, fatty-acid and glycan biosynthesis and metabolism, and sulfur relay. Near the deep-chlorophyll maximum depth, an exceptionally high percentage of AMGs related to photosynthesis was noticed. Taken together our findings suggest that the roles of viruses in the deep sea might be even more important than previously thought as they seem to orchestrate energy acquisition and microbial community dynamics, and thus, biogeochemical turnover in the oceans.


Introduction
Viruses are the most abundant biological entities in the oceans and have a large impact on biogeochemical processes.Indeed viruses modulate the cycling of nutrients, since organic matter is released from the lysed cells and returns to prokaryotes in readily-available dissolved form rather than being transferred to higher trophic levels (Suttle, 2007).Viral activity also results in the formation of gelatinous, aggregated particles from the products of lysis and through this process, organic material is diverted from surface layers to deeper ones (Zimmerman et al., 2020).
But even without lysing their hosts, viruses may reprogram hosts' metabolism and thus alter nutrients flow, through the expression of viral Auxiliary Metabolic Genes (AMGs).AMGs are expressed by viral genomes when the host is infected and still intact; their expression causes a "shift" in hosts' metabolism toward viral replication, for instance through enhancement of microbial nutrient uptake under low nutrients availability.This pattern was particularly studied under P limited conditions (Zeng and Chisholm, 2012;Kelly et al., 2013;Lin et al., 2016) and resulted in boosted host survival and viral replication too.Expression of AMGs has been described at the transcriptional level mostly in cultures (Lin et al., 2016;Bachy et al., 2018;Howard-Varona et al., 2020) and rarely in natural communities (Sieradzki et al., 2019;Heyerhoff et al., 2022) but also at the genome level in various ecosystems including the open and coastal ocean (Crummett et al., 2016;Luo et al., 2020;Tsiola et al., 2020;Jian et al., 2021).Special focus has been given on photosynthesis and carbonmetabolism related AMGs (Thompson et al., 2011) as well as on deciphering the wide variability of phage-host interactions that lead to totally different host metabolisms (Howard-Varona et al., 2018, 2020).
High-throughput exploration of viral populations is nowadays predominant (Guo et al., 2021;Nayfach et al., 2021;Pons et al., 2021;Gregory et al., 2022;Roux et al., 2023) in comparison to earlier studies that mostly applied flow cytometry (Marie et al., 1999) despite the fact that the usual nucleic acid stains may have underestimated viral counts.Viral community structure and diversity is now known across spatial and temporal scales (Mizuno et al., 2013;Hurwitz et al., 2015;Paez-Espino et al., 2016;Coutinho et al., 2019;Gregory et al., 2019;Luo et al., 2020;Tsiola et al., 2020;Wu et al., 2020).Global expeditions revealed that patterns of diversity are linked to temperature as well as oxygen and nutrients concentration and Prochlorococcus abundance (Brum et al., 2015;Coutinho et al., 2017;Gregory et al., 2019).Similarly to diversity, selection of AMGs was also found to be driven by environmental variables at the global and smaller scales, mainly by temperature and depth (Williamson et al., 2008).Temperature in specific seems to control not only diversity patterns (Steward et al., 2000;Hurwitz et al., 2013Hurwitz et al., , 2015) ) but also viral abundances and cytometric characteristics (Winter et al., 2012).
Next, to the so called "qualitative" implications, the most straightforward implications of viruses in the oceans are quantitative; viruses cause mortality, i.e., removal of standing stock and release of carbon that may reach approx.10 billion tons day −1 in the oceans (Suttle, 2007;Brussaard et al., 2008).Quantitative data, such as the frequency of lytically-infected and lysogenic cells, mostly originate by applying the viral reduction approach (Wilhelm et al., 2002;Winget et al., 2005).The dynamics of infections seem to be highly related to temperature (Mojica and Brussaard, 2014) and depth (Brum, 2005) but also to hosts' abundance and productivity (Brum et al., 2016) and organic matter and energy availability (Lara et al., 2017).
To link quantitative and qualitative implications, relevant data needs to be worked together.We hereby describe viral abundances and lysis rates, viral community structure, and AMG content in the lytic fraction (<0.2 μm) in relation to the environmental setting of a largely unexplored basin in the Eastern Mediterranean Sea (EMS).The EMS is a semi-enclosed, concentration-type basin where evaporation largely exceeds precipitation and river run-off.It has unique circulation patterns (Velaoras et al., 2019) and physicochemical features, including high bottom-water temperatures, salinity and transparency, extreme oligotrophy and limited microbial plankton growth associated mainly to low phosphate and inorganic nitrogen sources (Azov, 1991;Krom et al., 2004).We sampled the Levantine Basin, the easternmost part of the EMS, at 5, 50, 75, and 1,000 m depths in order to investigate water masses with discrete temperature, density, and levels of oxygen, nutrients and chlorophyll; these factors that are among the most frequent controlling the microbial communities, including the viral component (Techtmann et al., 2015).According to Velaoras et al. (2019) who performed the hydrographic study in this cruise, distinct water masses are encountered in the studied area from the surface to the seafloor.Down to 75 m, Atlantic Water and Levantine Surface Water masses are rather warm, saline and saturated in oxygen, with no major differences in current distributions across the vertical scale of this depth layer.Levantine Surface Water masses are mostly prominent in the easternmost side of the sampling area (LV13 and LV18 sampling stations, see section 2.1), while the westernmost side (LV3 and LV10 sampling stations, see section 2.1) is additionally affected by a branch of Mid-Mediterranean Jet that carries low-salinity Atlantic Water from the Ionian to the Levantine Sea.Deeper than 75 m, the intermediate layer extends down to ~165 m, and further below, three water masses can be identified: the Transitional Mediterranean Water between ~600 and 1,200 m, the old Cretan Deep Water between ~1,400 and 2,500 m, and the Eastern Mediterranean Deep Water below ~3,000 m.As part of a multidisciplinary sampling expedition with pre-defined sampling depths, we chose to collect water from 5 m (representing surface water), 75 m [representing the deep euphotic layer, approaching the deep-chlorophyll maximum (Livanou et al., 2019)] and 50 m (representing an additional depth within the mixed layer and the thermocline, with nearly equal primary productivity levels as the deep-chlorophyll maximum, not published data).All surface/ subsurface depths were chosen so they remain within the mixed layer (during the sampling time, stratification of the water column was just starting) and within the thermocline.Further, we chose to collect water from 1,000 m [representing a depth layer dominated by Transitional Mediterranean Water masses, rich in nutrients and of low salinity and oxygen (Velaoras et al., 2019)].We coupled traditional and modern ecological and genomic approaches to shed light on ecosystem processes that are affected by viral actions.We expected to find discrete viral populations in terms of activity and genomic signature in the coldest waters in the bathypelagic zone, dependent on temperature and energy.Sampling constraints and high analytical costs did not allow us to study a higher number of replicated sampling stations and seasons.Thus, we wish to point out that the presented findings should be generalized to other ecosystems/seasons with caution, and that further testing is needed to confirm our hypotheses.According to previous findings in the same project that revealed depth-specific eukaryotic community patterns (Santi et al., 2020) and other recent similar surveys (Coutinho et al., 2019(Coutinho et al., , 2023)), we hypothesized that a wider repertoire of viral AMGs characterize the cold bathypelagic waters, available to aid energy acquisition under resource-depleted conditions.

Sampling
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 waterrinsed low-density polyethylene containers, and further processed for the various analyses.

Assessment of viral and bacterial abundances
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 −5 and 4 × 10 −4 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.Yellowgreen 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).

Assessment of viral and bacterial production rates
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.000MWCO) 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 −1 h −1 ) was estimated from the slopes of the relationships between net increases in viral abundance over the respective time period of the net increase.

VP VLPmax VLPmin tmax tmin
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.

Assessment of viral metagenomic content
Seawater was filtered through 0.2 μm polycarbonate membranes under low vacuum.The 0.2-μm filtrate was chemically treated with 1 mg L −1 FeCl 3 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 2 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 ultracentrifugation 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. 1  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-tosample 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.

Statistics
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 4 3− , NO 3 − , NO 2 − , NH 4 , SiO 4 , 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.

Physical, chemical and biological features of the sampling area
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 (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).

Viral community composition
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 sum of nitrate and nitrite (NOx), silicate (SiO 4 ), total nitrogen (TN) and phosphorus (TP), dissolved oxygen (DO), dissolved organic carbon (DOC), phosphorus (DOP) and nitrogen (DON) and chlorophyll a (Chl) at the sampling stations.Abundances of virus-like particles (A: VLPs), heterotrophic bacteria (B), Synechococcus (C), Prochlorococcus (D), the ratio between virus-like particles and total bacteria (E: VBR), lytic viral production rate (F: Lytic VP) and heterotrophic bacterial production (G: BP) at the sampling area represented as average values and standard deviations that derive from the four sampling stations (LV3, LV10, LV13, and LV18) at four sampling depths (5, 50, 75, and 1,000 m).
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.

Viral AMG content
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.

Amino acid metabolism AMGs
"Amino acid metabolism" AMGs were mostly related to "cysteine and methionine metabolism" and "arginine and proline metabolism" contributing 59% ± 8% and 32% ± 7% within the category, respectively (Supplementary Table 3).The "cysteine and methionine metabolism" AMGs contributed significantly less within the category at 5 m than  the rest depths (one-way ANOVA, p < 0.05, post hoc Tukey test).The remaining percentage at 5 m was attributed to "glycine, serine and threonine metabolism" AMGs (Supplementary Table 3), involving glyA, GATM and serA.Within "cysteine and methionine metabolism, " the involved genes were DNMT1, DNMT3A, mtnN, mtn and pfs at all depths, while at 50 and 75 m the gene yrrT was also detected.The "arginine and proline metabolism" AMGs contributed significantly more within the category at 5, 50, and 75 m than 1,000 m (one-way ANOVA, p < 0.05, post hoc Tukey tests).

Metabolism of cofactors and vitamins AMGs
Only at 1,000 m, AMGs related to "metabolism of cofactors and vitamins" included "biotin metabolism" (fabG, fabF, fabZ, Figure 4).At surface/subsurface depths, the detected AMGs were involved Bubble chart reporting the percentage contribution of the specific AMG metabolic pathways (as defined by KEGG using the VIBRANT method) over the total number of AMGs identified at the four sampling depths.For 5, 50, and 75 m: averaged values derive from the four sampling stations (LV3, LV10, LV13, and LV18).For 1,000 m, the values of the two deep stations are presented separately (dark red: LV10 and light red: LV18).The dimension of each bubble is proportional to the relative percentage contribution of the specific AMG categories.

Lipid metabolism AMGs
The percentage of "lipid metabolism" AMGs significantly differed between surface/subsurface depths and 1,000 m (one-way ANOVA, p < 0.05) being higher at 1,000 m (post hoc Tukey test; Table 4).All AMGs in the category were related to "fatty acid biosynthesis" at 1,000 m (Supplementary Table 3; Figure 4) involving fabG, fabF, and fabZ.In the surface/subsurface depths, no "fatty acid biosynthesis" AMGs were found and the great majority were related to "biosynthesis of unsaturated fatty acids" (desC; Figure 4).

Same lysis rate at all depths-different taxa and AMGs involved
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)  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.

Change of amino-acid frequency aids adaptation to the cold in the deep
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).

Surface (5 m) vs. subsurface (50 and 75 m) patterns in AMG content
The most noticeable difference in the AMG content between between 5, 50, and 75 m was that of the DCM; at 75 m, photosynthesisrelated 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 10.3389/fmicb.2023.1271535Frontiers in Microbiology 13 frontiersin.org(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.

Conclusion
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 fieldwork 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.
FIGURE1(A) Canonical analysis of principal coordinates (CAP) of the seawater physicochemical variables (blue solid-lined circle), depicting the sample grouping based on "depth" and Pearson correlations of the physicochemical variables with the canonical axes.Explanatory variables were normalized and the normalized matrix was used to create a resemblance matrix using Euclidean distances.S, salinity; T, temperature; D, density, concentrations of NOx, dissolved inorganic nitrogen; PO4, phosphate; DO, dissolved oxygen; SiO4, silicate; DON, dissolved organic nitrogen; DOC, dissolved organic carbon; TN, total nitrogen; and TP, total phosphorus.The length and direction of the vectors indicate their relative strength and direction of relationship, respectively.(B) Principal coordinates analysis (PCoA) ordination plot of all Bray-Curtis distances.An abundance table of viral operational taxonomic units was generated via MetaPop and then, the abundance table was used to generate a Bray-Curtis distance matrix to visualize viral community distributions using a PCoA ordination method.The color of each circle represents the species richness within each sample.
FIGURE 3(A,B) Percentage contribution of the major viral families to total viral taxa at the sampling area represented as average values and standard deviations that derive from the four sampling stations (LV3, LV10, LV13, and LV18) at four sampling depths (5, 50, 75, and 1,000 m)."Other" families contribute <0.5% to the total reads and are expanded in the bottom left plot (at panel B: the 100% x axis refers to the "Other" fraction of panel A).

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 4

TABLE 1
Raw numbers of viral metagenome reads, number of reads after quality control, number of total contigs, total contig length, the N50 contig length, and Shannon index of the viral metagenomes at the sampling stations.
Samples from 1,000 m from stations LV3 and LV13 were not considered in this study since the raw number of reads was too low.

TABLE 2
Output of permutational analysis of variance (PERMANOVA) testing variation among sampling depths.

TABLE 4
Percentage contribution of the general AMG metabolic pathways (as defined by KEGG using the VIBRANT method) identified at the sampling stations.From left to right the full category description is: 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, Xenobiotics biodegradation and metabolism.
seem to drive differences in viral genomes from the community to the protein level.While the total