Viral potential to modulate microbial methane metabolism varies by habitat

Domain: earth_science
Published: 2024-02-29
Editorial outcome: major_revision

Abstract
Methane is a potent greenhouse gas contributing to global warming. Microorganisms largely drive the biogeochemical cycling of methane, yet little is known about viral contributions to methane metabolism (MM). We analyzed 982 publicly available metagenomes from host-associated and environmental habitats containing microbial MM genes, expanding the known MM auxiliary metabolic genes (AMGs) from three to 24, including seven genes exclusive to MM pathways. These AMGs are recovered on 911 viral contigs predicted to infect 14 prokaryotic phyla including Halobacteriota, Methanobacteriota, and Thermoproteota. Of those 24, most were encoded by viruses from rumen (16/24), with substantially fewer by viruses from environmental habitats (0–7/24). To search for additional MM AMGs from an environmental habitat, we generate metagenomes from methane-rich sediments in Vrana Lake, Croatia. Therein, we find diverse viral communities, with most viruses predicted to infect methanogens and methanotrophs and some encoding 13 AMGs that can modulate host metabolisms. However, none of these AMGs directly participate in MM pathways. Together these findings suggest that the extent to which viruses use AMGs to modulate host metabolic processes (e.g., MM) varies depending on the ecological properties of the habitat in which they dwell and is not always predictable by habitat biogeochemical properties.

Introduction
Earth is currently warming at an unprecedented speed over at least the last 2000 years1, partly owing to the increased concentration of greenhouse gases in the atmosphere2. Methane (CH4) is ranked second after carbon dioxide (CO2) in terms of the overall contribution to atmospheric warming and accounts for ~20% of the greenhouse gas-driven warming3–5. Approximately 50% of global methane emissions originate from aquatic ecosystems, of which freshwater lakes contribute up to 53%6. Wetlands are another important natural source of methane, whereas non-water-logged terrestrial ecosystems generally function as methane sinks7. Methane cycling is largely driven by microbes, with microbial methanogenesis (all mediated by archaea) producing ~69% of the total methane released to the atmosphere8. Among anthropogenic sources, about 30% of methane production is microbially mediated, almost all of which derives from ruminant livestock farming7. Understanding how cellular microbes and the viruses that infect them might impact methane metabolism (MM) across various habitats is therefore crucial to inform efforts to mitigate microbially driven methane emission and climate warming.

The impact of viruses on MM has only recently begun to be investigated9–11. Viruses are found ubiquitously in the environment and have important roles in ecological, biogeochemical, and evolutionary processes through cells lysis, horizontal gene transfer, and modulation of host metabolism (including carbon, sulfur, and nitrogen metabolism)12–15. Recently, some viruses have been found to encode pmoC and cofF as auxiliary metabolic genes (AMGs), with the potential to supplement the aerobic oxidation of methane by their bacterial host in freshwater lakes9. In addition, putative cofF and fae genes were found in viruses from deep‑sea hydrothermal vents11. Notably, MM AMGs were recently also found in a novel group of extrachromosomal elements called “Borgs”, which have been shown to encode mcr genes phylogenetically related to those of anaerobic methanotrophic archaea (ANME) in the genus Methanoperedens16. However, beyond these initial observations, little is known about other MM genes encoded by viruses or other extrachromosomal elements, or how they influence methane production.

Here, we sought to explore the potential effects of viruses on MM in habitats with microbially-derived MM. First, we analyzed 982 publicly available metagenomes from a range of environments, which are known from literature to potentially host methane-cycling microbial communities, and in which we were able to confirm the presence of microbial MM genes, to identify virus-encoded AMGs that could be involved in the MM pathway (MMP), including those participating in methane production (i.e., methanogenesis by archaea) and oxidation (either by aerobic methanotrophic bacteria or anaerobic methanotrophic archaea). We then generated an additional 11 metagenomes using lake sediments, in which methane emission has been detected, from Croatia’s largest freshwater lake (Vrana Lake)17 to sample bacterial/archaeal viruses and investigate their potential impacts on MM during infection.

Results
Some viruses encode genes that could modulate microbial methane metabolism
To discover new MM AMGs, 982 publicly available metagenomes from 15 environments (Supplementary Data 1; including rumen, marine water, marine sediment, lake water, lake sediment, river estuary sediment, wetlands sediment, and permafrost active layers, among others), were analyzed for microbial genes involved in MMP and viral genomes encoding MM AMGs. The assembled contigs excluding viral contigs, were used to identify microbial genes involved in MMP (based on their KEGG and PFAM annotations and the KEGG MM pathway modules18) and each of the environments contained 138–183 distinct microbial genes involved in MMP (in total 184 distinct genes from all environments after dereplication; Supplementary Data 2 & 3). Viral genomes were also identified from the assembled contigs of these metagenomes, using a combination of three tools: VirSorter19, DeepVirFinder20, and MARVEL21 (see Methods). In the identified viral genomes, we predicted and annotated viral genes and screened for putative virus-encoded MM AMGs using VIBRANT22 and manual curation. Particularly, MM AMGs were extracted based on their KEGG annotations and the MM pathway modules18. After rigorous inspection (see Methods), 911 viral contigs were identified to contain MM AMGs (Supplementary Data 4), resulting in the discovery of 24 distinct AMGs that potentially participate in 25 metabolic reactions in the MMP (Supplementary Data 5, Figs. S1 and S2). These 911 viral contigs originated from ~32% (316 of 982) of the here analyzed metagenomes (Supplementary Data 3). We compared the 911 viral contigs to the viral genomes/contigs from the NCBI RefSeq database (cultivable viral genomes) and IMG/VR database (uncultivated viral genomes from metagenomics)23 using a genome-based network approach (see Methods)24,25. About 34% of these viruses (308 of 911) could be assigned to taxonomy and all belonged to the class Caudoviricetes of the phylum Uroviricota, except one that belonged to an unclassified class of the phylum Nucleocytoviricota (Supplementary Data 4 and 6). About 28% (n = 257) of the 911 viruses were successfully linked to their microbial hosts (by iPHoP26) in four archaeal (Halobacteriota, Methanobacteriota, Thermoplasmatota, and Thermoproteota) and 10 bacterial (Actinobacteriota, Bacteroidota, Bdellovibrionota, Campylobacterota, Chloroflexota, Cyanobacteria, Firmicutes, Marinisomatota, Patescibacteria, and Proteobacteria) phyla (Supplementary Data 4). All their hosts contained genes (2 to 89 distinct genes) involved in MMP (Supplementary Data 7). About two-thirds of the host-linked viral contigs (163 of 257) encoded MM AMGs that were also detected in their hosts (Supplementary Data 4).

For each of the 24 MM AMGs, we selected one protein sequence from a highly confident viral contig (Supplementary Fig. S2) as an example to investigate the conserved domain and putative protein structure (see Methods). These in silico analyses revealed that all 24 AMGs exhibited the conserved functional domains and structural configurations (100% confidence for all the tested AMGs, except the fwdF gene with 99%; the confidence represents the probability that the match between the studied sequence and the template in the database is a true homology27) of their corresponding enzymes (Supplementary Data 5), suggesting that they likely encode functional AMGs. These results indicate that viruses could be largely underexplored players in ecosystem MM.

Investigating the metabolic roles of the 24 MM AMGs, we found that 17 of them could also participate in metabolic pathways other than MM, while the remaining seven AMGs (i.e., mtrA, pmoC, fwdF, fae, cofE, cofF, and frhB) exclusively participate in the MMP and thus had a high confidence in supporting direct viral modulations of microbial MM (Fig. 1; Supplementary Data 5, Figs. S1, S3, and S4). Among these seven AMGs, fwdF and fae were each detected on only one viral contig, while the others were identified on 3 to 25 viral contigs (see Supplementary Information for additional descriptions about viral contigs containing these seven AMGs). Functionally, the pmoC gene participates in the aerobic methane oxidation pathway of bacterial methanotrophs28, while the other six genes mtrA, fwdF, fae, cofE, cofF, and frhB are involved in the pathways of methanogenesis and/or anaerobic oxidation of methane (AOM)29 (Fig. 1A & Supplementary Fig. S1). The pmoC (methane monooxygenase subunit C) gene encodes a subunit of the particulate methane monooxygenase (pMMO) that catalyzes the aerobic oxidization of methane to methanol in bacteria (Fig. 1A & Supplementary Fig. S1)28. In the methanogenetic pathway, the mtrA (tetrahydromethanopterin S-methyltransferase subunit A) gene encodes a subunit of the membrane-associated multienzyme complex Mtr that transfers the methyl group of N5-methyltetrahydromethanopterin to coenzyme M (CoM) and produces Methyl-CoM30, which is an exergonic (ΔG°′ = −29 kJ/mole), sodium-ion-translocating step contributing to ion motive force in the methanogens’ energy metabolism. This energy conservation mechanism happens in all methanogens being able to produce methane from CO2 or acetate31. The gene product, Methyl-CoM, is essential for the final step of methanogenesis by methanogens32. The fwdF gene encodes an iron-sulfur protein as the subunit F of formylmethanofuran dehydrogenase, which can catalyze the reduction of methanofuran and CO2 to formylmethanofuran (Fig. 1A and Supplementary Fig. S1), in the first step of methanogenesis from CO233,34. The fae gene encodes the formaldehyde activating enzyme (Fae) catalyzing the condensation of formaldehyde with tetrahydromethanopterin (THMPT) to methylene-THMPT35, an intermediate in methanogenesis from CO2 (Fig. 1A & Supplementary Fig. S1). The remaining three genes cofE, cofF, and frhB are relevant to the synthesis of coenzyme F42036–38, which impacts the production of methylene-THMPT and 5-Methyl-THMPT, also intermediates for methanogenesis from CO2 (Fig. 1A and Supplementary Fig. S1).

Phylogenetic analyses of the above seven exclusive MM AMGs suggested that viruses have potentially acquired the mtrA genes (n = 3) from methanogens of the genus Methanobrevibacter (Euryarchaeota) (Fig. 1C & Supplementary Fig. S5A); the virus-encoded pmoC genes (n = 25) have potentially been transferred from methylotrophs belonging to several different genera within the phylum Proteobacteria, including Methylobacter, Methylomagnum, and Methylocystis (Supplementary Fig. S5B); and the fae gene (n = 1) might have been transferred from Methylophaga or Pseudomethylobacillus (Supplementary Fig. S5D). The remaining four genes (i.e., fwdF, cofE, cofF, and frhB) were more divergent from known and taxonomically characterized microbial genes, and thus could not be confidently linked to the potential gene transfer events from hosts to viruses (Supplementary Fig. S5C, E–G).

We assessed the habitats associated with each of the 24 MM AMGs, finding that host-associated samples (i.e., rumen) contained 16, whereas environmental habitats contained between one to seven MM AMGs, including marine water (7 AMGs), marine sediment (5), lake water (3), lake sediment (1), and hot spring sediment (2) (Fig. 2 & Supplementary Data 5). All 24 MM AMGs were also found on microbial contigs from the same environment where the AMGs were identified (Supplementary Data 2). Surprisingly, we did not find MM AMGs in some of the environmental habitats where we found 138–180 microbial MM genes, such as river estuary sediment, permafrost active layers, and wetlands sediment (Supplementary Data 2 and 3). Focusing only on AMGs involved in methane production, we identified 10 genes that can directly participate in or synthesize an intermediate for the pathway of methanogenesis from CO2 or acetate, including six that exclusively participate in MMP (i.e., mtrA, fwdF, fae, cofE, cofF, and frhB; Fig. 1A) and four that could also be involved in other metabolic pathways (ackA, pta, cooS, and glyA; Supplementary Fig. S1 & Data 5; see Supplementary Information for their potential roles in methane production pathways); nine of them came from host-associated rumen samples and only one to three were found in the environmental habitats including marine water, marine sediment, lake water, and lake sediment. Thus, despite representing less than 30% (286/982; Supplementary Data 1 and 3) of the metagenomes analyzed, host-associated samples (i.e., rumen) contained most of the identifiable MM AMGs (including those potentially participating in methane production), which were less common in environmental habitats (e.g., lake sediment, lake water, marine water, and marine sediment) where microbial MM genes (from 138 to 183 distinct genes) were also present. These results suggest that the extent to which viruses use AMGs to modulate host MM processes, including methane production, may vary depending on the habitats in which they dwell.

Vrana Lake sediment comprises mostly novel viral genera
Given that MM AMGs were apparently less common in publicly available metagenome datasets from environmental habitats, we adopted a targeted approach to look for additional MM AMGs from methane-rich environmental samples and further explore the impact of viruses on MM via host infection, by generating metagenomes from the methane-rich sediment of Vrana Lake in Zadar County, Croatia. Six pairs of bulk metagenomes and viromes were constructed for the Vrana Lake sediment (VLS) samples, recovered from two sediment cores at 50, 100, and 225 cm deep below the lake sediment surface. The cores were obtained from two sites within the lake: a muddy site consisting of organic-rich sediments within a concave depression (pockmark) of the sediment surface with fluid and gas (e.g., methane) efflux; and a sandy site in an area of sandy sediments and no visible pockmark depressions (Supplementary Fig. S6 & Data 8).

We recovered 3,260 viral contigs from the above VLS metagenomes. These contigs were clustered into vOTUs if they shared ≥95% nucleotide identity across 80% of their lengths39, resulting in 3,146 vOTUs (≥5 kb), including 1,050 “long” (≥10 kb) vOTUs (Supplementary Data 8). Taxonomic analyses, by comparing VLS viruses to viral genomes in both the NCBI RefSeq and IMG/VR databases (see Methods), revealed that most of the VLS long vOTUs (911 of 1,050) could not be taxonomically classified, indicating a high degree of novelty among VLS viruses. The remaining 139 vOTUs were assigned to Caudoviricetes, Faserviricetes, and Megaviricetes (Fig. 3A; Supplementary Fig. S7 & Data 9).

Viral communities differ between sediment sites and across sediment depths
The cellular microbial communities, investigated based on relative abundances of the 99 bacterial/archaeal metagenome-assembled genomes (MAGs) recovered from VLS metagenomes (Supplementary Data 10; see Methods), were distinct between muddy and sandy sites (Supplementary Fig. S8), which had very different physicochemical conditions (e.g., the total nitrogen and dissolved organic carbon were 8.9 and 2.3 times higher, respectively, in the muddy vs. sandy sediment; Supplementary Data 8). Similarly, the muddy and sandy sampling sites comprised mostly different viruses, with only 4.2% (131 of 3,146) of VLS vOTUs shared between sites (Fig. 3B). Ordination analysis, using the relative abundance data of vOTUs (Supplementary Data 11), confirmed that viral communities were significantly (p = 0.015) different between sites (Fig. 3C). Viral communities also varied with depth (i.e., 50, 100, and 225 cm deep), with only 5.7% (181 of 3,146) of vOTUs detected in samples from all three depths and the majority (75.0%; 2,360 of 3,146) of vOTUs being unique to a single depth (Supplementary Fig. S9A). Additionally, a comparison between bulk metagenomes and viromes found that 97.2% of VLS vOTUs (3,058 of 3,146) were retrieved exclusively from bulk metagenomes, suggesting that most of the recoverable VLS viruses might be within the cellular fraction captured by bulk metagenomes rather than in the viral particle fraction, though our data could not eliminate the possibility that the viromes might have only captured a subset of VLS extracellular viruses (e.g., some extracellular viruses might have been adsorbed to the sediment particles which were removed from viromes via filtering, but captured in bulk metagenomes) (Supplementary Information; Supplementary Fig. S9B & Data 11). For maximizing the virus recovery, we combined viruses identified from both bulk metagenomes and viromes for all further analyses.

Abundant viruses likely infect dominant microbes of the Thermoproteota and Chloroflexi to impact the sediment ecosystems
To explore the potential viral impacts on VLS ecosystems, we investigated virus-host linkages as reported previously (e.g., in soil and seawater13,40), via the iVirus tool VirMatcher41 that aggregates four different methods for host predictions (see Methods). Using the 99 VLS bacterial/archaeal MAGs as the host database (Supplementary Data 10; See Methods), we could link 2,167 of the 3,146 vOTUs (68.9%) to microbial hosts belonging to 17 different phyla (Fig. 4A; Supplementary Data 12). The VLS microbial communities were dominated by Thermoproteota (relative abundance: average 24.7% and range 12.7–41.4%; archaea) and Chloroflexi (average 23.5% and range 17.9–29.4%; bacteria) (Supplementary Fig. S10). We then calculated lineage-specific virus/host abundance ratios to assess viral infections for specific phyla and found that the most abundant VLS viruses were predicted to infect the above two most dominant microbial phyla, Thermoproteota and Chloroflexi (Fig. 4B; Supplementary Fig. S10). A substantial portion of VLS vOTUs were also linked to Desulfobacterota (Fig. 4A), which, however, was present at low levels or absent across the samples in terms of the identifiable MAG relative abundances (Supplementary Fig. S10). Some MAGs of Thermoproteota contained genes encoding for key steps of MM42–45, while some members of Chloroflexi are aerobic methanotrophs46,47 or are able to reduce sulfate to benefit methanogens/methanotrophs via syntrophic interactions48,49. In our data, we recovered 23 VLS MAGs belonging to the Thermoproteota, and each of them contained 24–102 (average 65) genes involved in MMP (Supplementary Data 13). Overall, these findings showed that the VLS viruses infected dominant microbial phyla, including ones involved in MM, and thus likely had an important impact on the sediment ecosystems.

Some Vrana Lake sediment viruses encode genes to modulate host metabolisms
The 99 VLS MAGs contained a total of 5,503 genes (136 distinct genes after dereplication) involved in MMP, including genes that can impact the key steps of MM such as the genes encoding methylenetetrahydromethanopterin dehydrogenase (Mtd), 5,10-methylenetetrahydromethanopterin reductase (Mer), tetrahydromethanopterin S-methyltransferase (Mtr), and methyl-coenzyme M reductase (Mcr) (Supplementary Data 13). To assess if VLS viruses also encode MM AMGs and thus could be modulating the hosts’ MM, we annotated genes for all the VLS vOTUs (Supplementary Data 14) and screened them for putative virus-encoded AMGs, including MM AMGs. After rigorous inspection (see Methods), we identified 13 putative AMGs from 11 VLS vOTUs (Supplementary Data 15). Interestingly, none of these AMGs were predicted to be directly involved in the MMP, which would agree with the inference of our analysis of publicly available metagenomes that the extent to which viruses modulate hosts’ MM may vary by habitats, and that MM AMGs seem to be less common in environmental habitats including lake sediments (<2% of the publicly available lake-sediment metagenomes had ≥1 MM AMG).

However, the fact that we identified 13 AMGs suggests that VLS viruses do still have the potential to modulate host metabolism in energy, carbohydrates, amino acids, nucleotides, cofactors, and vitamins (Supplementary Data 15). Particularly, we identified a vOTU that was predicted to infect a putative methanogenic Bathyarchaeia (a class of the phylum Thermoproteota) and which encoded a Thermoproteota-derived AMG for bacterioferritin (bfr) that oxidizes Fe2+ to Fe3+50 (Fig. 4C, D; Supplementary Fig. S11; Supplementary Information). Evolutionary pressure assessments within species and across lineages found that this virus-encoded Bfr was likely functional and under purification selection (pN/pS = 0; average dN/dS = 0.114; Supplementary Data 15 and 16). Iron is essential for numerous metabolic processes51, including microbial MM52–54. These results suggest that this virus might have the potential to modulate iron metabolism of a methanogenic host and thus indirectly impact MM in VLS (see Supplementary Information for additional discussion).

In silico analyses suggested that the 13 AMGs detected are likely functional. All of them had conserved functional domains (Supplementary Data 15), and when their protein sequences were structurally modeled using Phyre227, they had 100% confidence scores to their closest template proteins (Supplementary Data 15 and S10). Furthermore, microdiversity analyzes found that the pN/pS values, a proxy for gene selection pressure55,56, were <1 for all the testable VLS AMGs, suggesting that they were under purifying selection (Supplementary Data 15). While no AMGs that can directly participate in MMP were detected from these methane-rich lake sediments, these results indicate that VLS viruses encode functional AMGs that likely alter microbial metabolisms in the Varana Lake sediments, including an AMG that have an indirect influence on MM through manipulating a putative methanogen’s iron metabolism.

Methods
Published metagenome analyses
To investigate how viruses might modulate hosts’ metabolic processing to participate in methane cycling, we analyzed 982 publicly metagenomes from both host-associated and environmental habitats that contained 138–183 genes involved in MMP (Supplementary Data 1, 2, and 3; including rumen, marine water, marine sediment, lake water, lake sediment, river estuary sediment, hot spring, and permafrost active layers, among the 15 habitats). Depending on data availability (indicated in Supplementary Data 1), these metagenomes were analyzed by assembling contigs, identifying viral genomes, annotating viral genes, and/or rigorously screening them for putative virus-encoded MM AMGs, as described in below method sections.

VLS site characterization and field sampling
Two deep sediment cores were collected in 2015 from two sites of the Vrana Lake in Zadar, Croatia: One core was obtained within a pockmark depression in muddy sediment area (muddy site), and a second core was sampled from a sandy sediment area with no visible pockmarks (sandy site; Supplementary Fig. S6 & Data 8). Three sediment samples were collected from each core, at 50, 100, and 225 cm, respectively, below the lake sediment surface of each site. These six sediment samples were frozen at −20°C once sampled in the field, and then were transported to the laboratory, where they were stored at −20°C for further analyses, including filtration and DNA extraction.

Sample processing and genomic DNA isolation
Each sample (0.5 g sediment) was used for bulk DNA extraction with a DNeasy PowerSoil Isolation Kit (Cat No. 12888-100, QIAGEN) according to the manufacturer’s instructions. In addition, the extracellular viruses were extracted from each sample (0.9 g sediment) by suspending the sediment using AKC buffer (1% potassium citrate, 1% PBS, and 150 mM MgSO4) by horizontally shaking at 400 rpm for 15 min at 4°C, according to a previously established protocol65. The liquid suspension (about 12 mL) was then passed through a polycarbonate 0.22-μm-pore-size filter (Cat No. GTTP02500, Isopore) to remove cells and particles >0.22 μm. Samples were incubated with 100 U DNase I per 1 mL of sample (Roche) with DNase I reaction buffer (final 10 mM Tric-HCl, 2.5 mM MgCl2, 0.5 mM CaCl2, pH 7.6) at 4°C for 48 hr. DNase was inactivated by addition of EDTA and EGTA to a final concentration of 100 mM. The virus-like particles in the filtrate were concentrated to 0.5 mL using 100 kDa Amicon Ultra Concentrators (EMD Millipore, Darmstadt, Germany) and preserved at 4°C until DNA extraction (within 2 hours). Genomic DNA from viral concentrates was isolated using the same protocol as isolating the bulk DNA above. Both bulk and viral DNA were preserved at −20°C until further processing.

Metagenomic sequencing
Theoretically, bulk DNA was able to capture all viruses (both intracellular and extra-cellular viruses) from the sediments, while the viral DNA extracted from the filtrates specifically captured the extra-cellular viruses. To maximize viral discovery and gain insight into the proportion of the extra-cellular viruses in VLS, this study analyzed both bulk and viral DNA, which were subjected for bulk and viral metagenome (virome) sequencing, respectively. All metagenomes (i.e., six bulk metagenomes and six viromes) were sequenced at the Joint Genome Institute (JGI), Department of Energy, USA. Briefly, the DNA libraries were prepared using the Nextera® XT Library Prep Kit (Cat No. 15032354, Illumina) and sequenced on the Illumina NovaSeq platform (2×150 bp). Sequencing failed for one bulk metagenome sample (i.e., M50_M), which was collected from 50 cm sediment deep of the muddy site core (within a pockmark); thus this sample only had a virome (i.e., M50_V) for further analyzes (Supplementary Data 8).

Metagenomic read processing and viral identification
Metagenomic data analyses were supported by the Ohio Supercomputer Center, unless stated otherwise. Sequencing reads were filtered for quality by JGI using their previously established standard pipeline66, generating a total of 9.5 × 1010 bases of sequencing data (range 0.3–1.5 × 1010 bases, average 8.6 × 109 bases per library; Supplementary Data 8). Then the metagenomic sequence data was assembled to contigs by metaSPAdes67, using a previously established pipeline for assembling pre-amplified metagenomes (parameters: read deduplication + read error correction + --sc + -k 21,33,55,77,99,127)68. The assembled contigs (length ≥5 kb or circular contigs with length 1.5–5.0 kb) from all metagenomes were used for identifying viruses following previously described methods69, as also described below. Three tools VirSorter v1.1.019, DeepVirFinder v1.020, and MARVEL v0.221 were used for predicting viruses. Contigs were classified as viruses if they met one of the following four criteria: (i) Categories 1, 2, 4, or 5 of VirSorter v1.1.0; (ii) DeepVirFinder score ≥0.9 and p < 0.05; (iii) MARVEL probability score ≥90%; or (iv) DeepVirFinder score ≥0.7 and p < 0.05 and MARVEL probability score ≥70%. Viral contigs identified by the above methods were combined for further analyses.

Viral contigs were first inspected and filtered for potential contaminants by comparing them to viral genomes considered as putative laboratory contaminants (e.g., phages cultivated in our laboratory: Synechococcus, Cellulophaga, and Pseudoalteromonas phages) using Blastn. The remaining contigs were clustered into vOTUs (~species-level taxonomic unit) if they shared ≥95% nucleotide identity across 80% of their lengths39. The longest contig within each vOTU was selected as the seed sequence to represent that vOTU. These efforts generated a total of 3,260 viral contigs, that were clustered into 3,146 vOTUs, including 1,050 “long” vOTUs with length ≥10 kb. The coverages of vOTUs (≥5 kb) were generated using the iVirus’ Bowtie Batch and Read2RefMapper tools, by mapping quality-filtered reads to vOTUs, and the resulting coverage depths were normalized by library size to “coverage per gigabase of virome” to assess the viral communities in VLS41,70.

Taxonomy and ecology analyses
Because viruses lack any single, universally shared gene, we established taxonomy using gene-sharing network analysis from viral sequences ≥10 kb in length using vConTACTv224. Briefly, this analysis compared the 1,050 “long” VLS vOTUs and the 911 public datasets-originated viruses that contained MM AMGs to viral genomes in the National Center for Biotechnology Information (NCBI) RefSeq database (release v201) and the IMG/VR v4 database and generated viral clusters approximately equivalent to known viral genera13,24,71. Principal coordinate analyzes (PCoA) were performed using Bray Curtis distance matrices based on the coverage of each vOTU. PERMANOVA (Permutational Multivariate Analysis of Variance; permutations=999) tests72 were used to calculate the statistical differences in communities between both sampling sites and metagenome types.

Microbial genomic analyses
For microbial genomic analyses, quality-controlled reads of the bulk metagenomes were co-assembled using metaSPAdes v3.11.167. The assembled contigs (≥1.5 kb) were then used to bin microbial metagenome-assembled genomes (MAGs), by MetaBat2 v2.12.173 using each present binning strategy with and without contig coverage profiles74. A total of 99 MAGs, with medium to high quality (completeness ≥40% and contamination ≤10%, via CheckM v1.1.1075), were generated and then were assigned to a taxonomy using GTDB-Tk v1.3.076,77. Assembly, binning, and quality estimation of prokaryotic MAGs was performed at the Life Science Compute Cluster (https://lisc.univie.ac.at) at the University of Vienna. These MAGs were dereplicated to 83 MAG populations using dRep v1.0.0 with default parameters (sharing ≥95% nucleotide identity across ≥10% of their length)78. Metagenomic reads were mapped to MAG populations to characterize their relative abundances using CoverM v0.3.2 with default parameters (https://github.com/wwood/CoverM).

Viral host prediction
The putative virus-host linkages were predicted in silico using the iVirus tool VirMatcher41, which aggregates four different methods to provide a statistical confidence score for each host prediction and these methods are based on: (i) tRNA match, (ii) nucleotide sequence composition, (iii) nucleotide sequence similarity, and (iv) CRISPR spacer match. Since viral host prediction benefits from the database that contains microbial genomes from the same ecosystems as viruses40, we used the microbial MAGs (n = 99), that were recovered from the VLS bulk metagenomes described above, as the microbial database for linking the VLS viruses to their hosts. A summary of the host predictions is available in Supplementary Data 12. The lineage-specific virus/host abundance ratios at phylum level were assessed by comparing the relative abundances of microbial phylum and viruses infecting each phylum40.

Virus-encoded AMG identification
The putative AMGs were identified and evaluated for viruses recovered from both 982 publicly published metagenomes from a range of environments where microbial MM genes were detected (see next paragraph) (Supplementary Data 1, 2, and 3) and the 11 VLS metagenomes originally constructed in this study, according to our previously established methods79. Specifically, once viral contigs were recovered from metagenomes, they were processed with VIBRANT to obtain gene functional annotations against the KEGG and PFAM databases and identify putative AMGs by the default parameters22. To obtain high-quality and rule out false-positive AMGs from microbial contamination, CheckV (with default parameters, v0.3.0) and manual inspection were then used to assess host-virus boundaries and remove the potential host fraction of the viral contigs80. Only AMGs that were surrounded by phage genes, did not contain transposon regions, and had consistent annotations between the KEGG and PFAM databases were included for further analyzes. Metabolism categories of AMGs, including those participating in MMP, were summarized based on KEGG annotations and the pathway modules18.

For the 982 published metagenomes, we used their pre-existing viral contigs for AMG recovery if the data was publicly available, and otherwise de novo assembled the metagenomes (n = 265), recovered viral contigs, and/or identified putative AMGs (the data type used for each of the 982 metagenomes is indicated in Supplementary Data 1), using the methods described in the preceding sections. The KEGG IDs of AMGs were used for extracting the genes that could be involved in the MMP, resulting in a discovery of 24 distinct AMGs (on a total of 911 viral contigs containing ≥1 MM AMG) that could participate in 25 steps of the MMP, including seven genes that exclusively participate in MMP (Supplementary Data 5 & Fig. S1). Hosts of the 911 viral contigs containing ≥1 MM AMG were predicted by iPHoP26 (confidence score ≥90%), resulting successfully virus-host linkages for 257 viral contigs (Supplementary Data 4). Of the 24 AMGs, 17 were detected from only one environment type, while the remaining seven were found in two to four environment types; similarly, six AMGs were identified on only one viral contig and the other 18 AMGs were found on two to 60 viral contigs, except one AMG glyA that presented on 642 viral contigs (Supplementary Data 5). For each of the 24 AMGs, we used one viral contig/genome (a contig with a highest viral quality score was selected as the representative, via CheckV’s assessment, if the AMG was presented on more than one contig) to illustrate the viral genomic context and AMG position (Supplementary Data 5 & Fig. S2). In addition, the AMGs of the selected viral contigs were further used as examples for analyzing conserved domains by comparing to the domains in the public Conserved Domain Database (CDD v3.20) via NCBI CD-Search81 and predicting three-dimensional protein structures by Phyre227 (Supplementary Data 5). The assembled contigs, excluding viral contigs, were annotated by DRAM against the KEGG and PFAM databases by the default parameters and further used for recovering microbial MM genes based on their KEGG and PFAM annotations (with consistent annotation in the two databases) and the KEGG MM pathway modules18.

Visualization of the genome maps for the viruses was performed using Easyfig v2.2.582. Phage genes, hallmark genes, and potential cellular genes were identified by VIBRANT, CheckV, and VirSorter19,22,80,83. Protein sequences from the AMGs were structurally modeled using Phyre227 in normal modeling mode to confirm and further resolve functional predictions. These seven exclusive MM AMGs (mtrA, pmoC, fwdF, fae, cofE, cofF, and frhB) and the VLS AMG bfr were subjected to phylogenetic analyses to infer its evolutionary history. DIAMOND (v2.0.15) BLASTP84 was used to query the gene’s amino acid sequence against the NCBI RefSeq database (release v214) in a sensitive mode with default settings, to obtain the top 40 hits (top 20 hits if an AMG was identified on more than two viral contigs) as the reference sequences. In addition, microbe-encoded bfr genes were extracted from the VLS microbial metagenomes to study possible gene transfers between viruses and their microbial hosts. Multiple sequence alignment was performed using MAFFT (v. 7.017)85 with the E-INS-I strategy for 1000 iterations. The aligned sequences were then trimmed using TrimAl86 with the flag gappyout. The substitution model was selected by ModelFinder87 for accurate phylogenetic analysis. Phylogenies were generated using IQ-TREE88 with ultrafast 1,000 bootstrap replicates, and then visualized in iTOL (v5)89. Potential recombination among genes was evaluated using nine programs: RDP90, GENECONV91, BootScan92, MaxChi93, Chimaera94, SiScan95, LARD96, Phylpro97, and 3Seq98 within RDP5 (v5.23)99. A Bonferroni correction with a p value cut-off of 0.05 was applied in each of the tests. A sequence was considered as a true recombinant if supported by at least four of the nine programs. The selection pressure (pN/pS) of VLS AMGs were calculated by recruiting VLS metagenomic reads to the AMG-containing vOTUs and identifying the SNPs on AMGs, using the tool MetaPop v1.0 through default parameters56. For the VLS AMG bfr, branch and site selection pressure (dN/dS) analysis across lineages was carried out using codon models with maximum likelihood estimated with the codeml package in PAML (v4.9)100 (Supplementary Data 16).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Zhong et al. describe a comprehensive computational search for viral genes involved in methane metabolism as well as an equivalent computational search in a new metagenomic dataset from methane-rich lake sediments. This is an interesting scientific question. The methods the authors use to assemble and identify viral genomes in the metagenomic assemblies are robust. Fundamentally, their primary conclusion, the report of 24 viral AMGs involved in methane metabolism, is most likely accurate. However, crucial analyses of these viral genomes are missing, and there are some organizational issues in the manuscript. Specifically:\n1. The authors need to provide an analysis of what the likely hosts are of the viruses encoding viral AMGs. What general classes of phage are these viruses? What are their likely or possible hosts? Do their hosts have genes involved in methane metabolism? Can the authors use a broad host-prediction tool like iPhoP to aid in these predictions?\n2. Some conclusions by the authors consist of comparing the distribution of viral AMGs across environments expected to have microbial methane metabolism. But can they quantify the presence of prokaryotic methane pathways in these environments to shed light on this? How do viruses encoding methane metabolism genes correlate with prokaryotes that do so? Can a virus potentially provide a gene otherwise entirely missing from a host, or are they augmenting existing host metabolism? These are the pressing questions that are left unanswered by the manuscript.\nSome specific comments:\nAbstract: the abstract should focus a bit more on what the viruses encoding the MM AMGs are: what types of viruses encode these genes? What organisms might they infect? This information is missing from the abstract.\nLine 84: publicly available metagenomes from environments highly likely to have active biological MM\ncan the authors explain more about how they evaluated whether a metagenome fit this criteria?\nLine 87: \"isolating viral genomes\" can the authors adjust terminology (e.g. \"identifying\" viral genomes) and briefly mention the methodology used here?\nLine 88-90: More details should be in the results describing how MM genes were identified. What qualifies as an MM gene? How were they identified? Just the most basic description of the methods should go here so that a reader does not have to skip to the methods to understand what is actually being reported.\nLine 95: Please specify what the \"confidence\" reported here actually is and what it implies.\nLine 138-140: \"Where we expect active biological MM...\" So, were prokaryotic MM genes present in these environments? How do viral MM genes correlate with the presence of prokaryotic MM genes, and/or prokaryotes involved in methane metabolism? This crucial comparison is missing from the work.\nFor example, in a later section, it is noted that \"The VLS microbial communities were dominated by Thermoproteota\" and \"Many Thermoproteota members are archaeal methanogens\" ... but from the MAGs can the authors identify genes for methanogenesis? What about the methanotrophs?", "reviewer_2": "Microorganisms are main players in the cycling of the greenhouse gas methane, but recent discoveries of extrachromosomal elements associated with these microbes suggest that they influence the cycling of CH4, too. The authors Zhong et al. now analyze 982 publicly available metagenomes and generate 11 new metagenomes from a freshwater lake in Croatia, to identity viruses associated with methane metabolizing microorganisms. They find viruses that seem to infect methanogens and aerobic methanotrophs in the existing and their own datasets and then analyze the AMGs encoded on these viruses/phages. This is a generally very interesting topic and the workflows described seem well-chosen. Yet some of the data presented seems a bit disconnected and preliminary for publication. For example, it is not surprising that viral communities differ between sediment sites because the microorganisms differ between those sites. Thus it would be important to resolve the microbial community of the sites or at least present those microorganisms that are involved in methane cycling in these sites (and present this data), and then present the viral MAGs found in those sites. It is also not clear how novel the viral contigs are since they were only compared to cultivable viruses and a selection of metagenomic datasets. There is also no clear boundary between methanogenic archaea and methanotrophic bacteria. A reader that is no expert in this would likely be very confused as everything is presented together. This manuscript unfortunately requires a major revision before it can be published.\nMajor comments:\n- Headers of sections are imprecise. E.g. “Viruses encode 24 AMGs modulating…” suggests this is a general finding. It would be better to rephrase to “Some viruses encode genes that could augment/modulate microbial methane metabolism”. Generally, abbreviations should be avoided in results sections because the reader will be lost..\n- L: 100-110: The aerobic methane oxidation should be more clearly separated form (reverse) methanogenesis, since these are very distinct processes mediated by completely different microbes\n- L111f: what is cool about the finding that mtrA is encoded on a viral contig is that the Mtr complex is directly involved in energy conservation. This part could be really improved by diving in a bit deeper into the function and stating this!\n- L169: why were the vOTUs only compared to cultivable viruses? This of course will lead to the (likely false) inferral that they are all novel… There is a huge database of viruses/phages derived from metagenomics (and were thus not cultivated) and the viral contigs should be compared to those? The authors then do compare them to “selected published environmental metagenomes”(384) but how and why were these selected?\n- L:254: how do the viruses alter the ecophysiology of the sediments? That seems like a really huge claim with no evidence to support this\n- L259-275: This reads like a summary or even abstract and not like a discussion.\n- L284-287: where does this hypothesis come from and how does the data presented in this paper connect to it?\n- Figure 1A: The representation is extremely hard to follow. It also suggests that pmoC is in one cell (virus ?) with the other enzymes? Gene names are to be written in italics. This Figure should be revised.\n- Figure 1C: Is it necessary to show this tree in such a collapsed way?It seems to provide little information. As a sidenote: protein names should be capitalized.\n- Figure 2: Is it necessary to show this as a main figure? I am more confused by looking at it than having read the section in the text, unfortunately.\n- The Figures generally look a bit unpolished… e.g., Figure 3 (A) and (B) are outside the boxes, where (C) and (D) are inside. There are lots of pie charts and circles in several figures.. Figure 4A why does the green arrow restate in green font what is written in the left pie-chart? Were all 2,167 vOTUs assigned to either bacteria or archaea? This seems suprising that taxonomic classification/linkage is possible for every vOTU found..\nMinor comments:\n- L28: “methanogenic sediments” seems very colloquial, maybe rephrase to “sediments where methane emissions/consumption has been detected”(or whatever is applicable)\n- L90: resulting in the discovery..\n- L93L if only one protein sequence was investigated, how can the “conserved domains” be investigated ?\n- L104: on only one viral contig (not “from”)\n- L106: it should be stated that pmoC is part of the aerobic methane oxidation pw of bacteria\n- L107: how do the authors know that these genes are not involved in AOM ? Unless there is evidence for this (e.g., the metagenomic samples these viral contigs originate from, only comprise methanogens, and not ANME). If this was not investigated, it should be added “pathway of methanogenesis or reverse methanogenesis”\n- L168: how were the vOTUs established?\n- L 211: how were vOTUs linked to microbial hosts? This should be stated in brief in the results.\n- L221: I don’t think it is accurate to write “many” Thermoproteota members are methanogens. Also not sure why “archaeal” is required at this point, because all methanogens (we know so far!) are archaeal (this would fit in the beginning, but not here)\n- L224: had an important\n- L223-5: This seems pretty disconnected. What is the evidence for this or why bring this up here?\n- L242: “conserved functional domains of their corresponding enzymes” What does this mean?\n- L:251: Why was this pN/pS analysis done?\n- L282: “significant” may not be the best term here\n- L293: CO is also a substrate for methanogens\n- L:294: It is a bit odd to have (n=17) in the discussion section as such specific values usually belong to the results\n- L297ff: it should be stated that mcr genes were recently found in “Borgs” – archaeal ECEs that are neither typical viruses nor plasmids"}, "author_response": "Response #1: We appreciate your kind words in assessing our methods and primary conclusion and thank you for the constructive criticisms on the original manuscript’s shortcomings, which greatly helped us improve our work during revision. In the revised version, we performed additional analyses as requested:\n(1) We identified and added information on the putative hosts for these viruses. Please see Response #2 for further details.\n(2) We assigned the taxonomy of these viral genomes by comparing them to viruses from the NCBI RefSeq and IMG/VR databases. Please see Response #2 for further details.\n(3) We assessed the microbial genes involved in MMP in the same environments where we searched for virus-encoded MM AMGs. This assessment helps shed light on the later work about identifying virus-encoded MM AMGs in these environments, and thus greatly improves the contextualization of our data. Please see Response #3 for further details.\nPlease note that all the line numbers cited in our responses are based on the revised manuscript with tracked changes.\n1. The authors need to provide an analysis of what the likely hosts are of the viruses encoding viral AMGs. What general classes of phage are these viruses? What are their likely or possible hosts? Do their hosts have genes involved in methane metabolism? Can the authors use a broad host-prediction tool like iPhoP to aid in these predictions?\nResponse #2: Thank you for these thoughtful suggestions. In the revised manuscript, we used iPHoP to predict hosts for the 911 viral contigs that we detected to encoded MM AMGs. The results are provided in Table S4 (Columns Y–AD; a general summary of viral hosts) and Table S7 (full gene annotations of viral hosts). Briefly, 257 of the 911 viruses were successfully linked to their hosts. Notably, all the predicted hosts contained genes (between 2 and 89) involved in MMP, and about two-thirds of the host-linked viral contigs (163 of 257) encoded MM AMGs that were also detected in their hosts (Column AA in Table S4). To contextualize these results, the below sentences were added to the revised main text (Lines 121-128):\n“About 28% (n=257) of the 911 viruses were successfully linked to their microbial hosts (by iPHoP25) in four archaeal (Halobacteriota, Methanobacteriota, Thermoplasmatota, and Thermoproteota) and 10 bacterial (Actinobacteriota, Bacteroidota, Bdellovibrionota, Campylobacterota, Chloroflexota, Cyanobacteria, Firmicutes, Marinisomatota, Patescibacteria, and Proteobacteria) phyla (Supplementary Table S4). All their hosts contained genes (2 to 89 distinct genes) involved in MMP (Supplementary Table S7). About two-thirds of the host-linked viral contigs (163 of 257) encoded MM AMGs that were also detected in their hosts (Supplementary Table S4).”\nRegarding the comment on “general classes of phages for these viruses”, we assigned the taxonomy of viruses that encoded MM AMGs by comparing their genomes to the viral genomes in NCBI RefSeq and IMG/VR databases via the tool vConTACT v2. About 34% of these viruses (308 of 911) could be assigned a taxonomy and all of them belonged to the class Caudoviricetes of the phylum Uroviricota, except one viral contig that belonged to an unclassified class of the phylum Nucleocytoviricota. These results were provided in Table S4 (Column X) and Table S6, as well as described in the main text (Lines 115-121):\n“We compared the 911 viral contigs to the viral genomes/contigs from the NCBI RefSeq database (cultivable viral genomes) and IMG/VR database (uncultivated virus genomes from metagenomics)22 using a genome-based network approach (see Methods)23,24. About 34% of these viruses (308 of 911) could be assigned to taxonomy and all belonged to the class Caudoviricetes of the phylum Uroviricota, except one that belonged to an unclassified class of the phylum Nucleocytoviricota (Supplementary Table S4 & Table S6)”\n2. Some conclusions by the authors consist of comparing the distribution of viral AMGs across environments expected to have microbial methane metabolism. But can they quantify the presence of prokaryotic methane pathways in these environments to shed light on this? How do viruses encoding methane metabolism genes correlate with prokaryotes that do so? Can a virus potentially provide a gene otherwise entirely missing from a host, or are they augmenting existing host metabolism? These are the pressing questions that are left unanswered by the manuscript.\nResponse #3: Thanks for your thoughtful comments. We now also annotated the prokaryotic contigs of the publicly available 982 metagenomes, where we searched for viral MM AMGs, and explored the genes involved in MMP. These metagenomes were from 15 environments, in which the prokaryotes contained 138–183 distinct genes (i.e., after dereplication) involved in MMP (see Table S2 & Table S3). All the 24 viral MM AMGs were also detected from prokaryotes, in the same environment where the viral and prokaryotic genes originated. However, the number of viral and prokaryotic MM genes were uncorrelated across these environments (see Table S3). These results were discussed in the revised main text:\nLines 96-104: “To discover new MM AMGs, 982 publicly available metagenomes from 15 environments (Supplementary Table S1; including rumen, marine water, marine sediment, lake water, lake sediment, river estuary sediment, wetland sediment, and permafrost active layers, among others), were analyzed for microbial genes involved in MMP and viral genomes encoding MM AMGs. The assembled contigs excluding viral contigs, were used to identify microbial genes involved in MMP (based on their KEGG and PFAM annotations and the KEGG MM pathway modules17) and each of the 982 metagenomes contained 138–183 distinct microbial genes involved in MMP (in total 184 distinct genes from all environments after dereplication; Supplementary Table S2 & Table S3).”\nLines 178-186: “We assessed the habitats associated with each of the 24 MM AMGs, finding that host-associated samples (i.e., rumen) contained 16, whereas environmental habitats contained between one to seven MM AMGs, including marine water (7 AMGs), marine sediment (5), lake water (3), lake sediment (1), and hot spring sediment (2) (Fig. 2 & Supplementary Table S5). All 24 MM AMGs were also found on microbial contigs from the same environment where the AMGs were identified (Supplementary Table S2). Surprisingly, we did not find MM AMGs in some of the environmental habitats where we found 138–180 microbial MM genes, such as river estuary sediment, permafrost active layers, and wetland sediment (Supplementary Table S2 & Table S3).”\nSome specific comments:\nAbstract: the abstract should focus a bit more on what the viruses encoding the MM AMGs are: what types of viruses encode these genes? What organisms might they infect? This information is missing from the abstract.\nResponse #4: We have revised the abstract to include the information on taxonomy and putative hosts of these viruses (Lines 34-40):\n“By analyzing viruses … we expand this to 24 viral MM AMGs identified on 911 viral contigs, including seven genes that are exclusive to MM pathways. About 34% of the 911 viruses could be taxonomically assigned and almost all belonged to the Caudoviricetes; ~28% were successfully linked to their hosts, in four archaeal and 10 bacterial phyla, including Halobacteriota, Methanobacteriota, and Thermoproteota.”\nLine 84: publicly available metagenomes from environments highly likely to have active biological MM\ncan the authors explain more about how they evaluated whether a metagenome fit this criteria?\nResponse #5: Yes. In the revised version, we added the analyses of microbial genes in these environments (see Response #3). The results suggested that each of these environments contained 138 to 183 distinct microbial genes involved in MMP and they were thus used to search for viral MM AMGs. Overall, we initially selected metagenomes from environments, based on literature knowledge, known to potentially host methane-cycling microbial communities. Further, we annotated their microbial genes and confirmed the presence of microbial MM genes in these environments. We have modified the below text to explain this:\nLines 83-86: “First, we analyzed 982 publicly available metagenomes from a range of environments, which are known from literature to potentially host methane-cycling microbial communities, and in which we were able to confirm the presence of microbial MM genes, to identify virus-encoded AMGs that could be involved in the MM pathway (MMP), including …”\nLines 96-104: “To discover new MM AMGs, 982 publicly available metagenomes from 15 environments (Supplementary Table S1; including rumen, marine water, marine sediment, lake water, lake sediment, river estuary sediment, wetland sediment, and permafrost active layers, among others), were analyzed for microbial genes involved in MMP and viral genomes encoding MM AMGs. The assembled contigs excluding viral contigs, were used to identify microbial genes involved in MMP (based on their KEGG and PFAM annotations and the KEGG MM pathway modules17) and each of the 982 metagenomes contained 138–183 distinct microbial genes involved in MMP (in total 184 distinct genes from all environments after dereplication; Supplementary Table S2 & Table S3).”\nLine 87: \"isolating viral genomes\" can the authors adjust terminology (e.g. \"identifying\" viral genomes) and briefly mention the methodology used here?\nResponse #6: We have adjusted the terminology and added a brief description of methods for identifying viral genomes. The modified text reads as (Lines 108-110):\n“Viral genomes were also identified from the assembled contigs of these metagenomes, using a combination of three tools: VirSorter18, DeepVirFinder19, and MARVEL20 (see Methods).”\nLine 88-90: More details should be in the results describing how MM genes were identified. What qualifies as an MM gene? How were they identified? Just the most basic description of the methods should go here so that a reader does not have to skip to the methods to understand what is actually being reported.\nResponse #7: We have modified these sentences to include a description of the methods for identifying MM genes. They read as (Lines 110-116):\n“In the identified viral genomes, we predicted and annotated viral genes and screened for putative virus-encoded MM AMGs using VIBRANT21 and manual curation. Particularly, MM AMGs were extracted based on their KEGG annotations and the MM pathway modules17. After rigorous inspection (see Methods), 911 viral contigs were identified to contain MM AMGs (Supplementary Table S4), resulting in the discovery of 24 distinct AMGs that potentially participate in 25 metabolic reactions in the MMP (Supplementary Table S5, Fig. S1, & Fig. S2).”\nLine 95: Please specify what the \"confidence\" reported here actually is and what it implies.\nResponse #8: Thank you for pointing out this omission. The confidence represents the probability (from 0% to 100%) that the match between the studied sequence and the template in the database is a true homology. A higher confidence implies a higher accuracy of the protein structure prediction. We have modified the text to include an explanation of the “confidence”, it reads as (Lines 132-136):\n“These in silico analyses … and structural configurations (100% confidence for all the tested AMGs, except the fwdF gene with 99%; the confidence represents the probability that the match between the studied sequence and the template in the database is a true homology26) …”\nLine 138-140: \"Where we expect active biological MM...\" So, were prokaryotic MM genes present in these environments? How do viral MM genes correlate with the presence of prokaryotic MM genes, and/or prokaryotes involved in methane metabolism? This crucial comparison is missing from the work.\nFor example, in a later section, it is noted that \"The VLS microbial communities were dominated by Thermoproteota\" and \"Many Thermoproteota members are archaeal methanogens\" ... but from the MAGs can the authors identify genes for methanogenesis? What about the methanotrophs?\nResponse #9: Please refer to “Response #3” for our arguments about the prokaryotic MM genes present in these environments. We have modified this sentence as (Lines 184-187):\n“Surprisingly, we did not find MM AMGs in some of the environmental habitats where we found 138–180 microbial MM genes, such as river estuary sediment, permafrost active layers, and wetland sediment (Supplementary Table S2 & Table S3).”\nFor the comment on “Thermoproteota members”: We recovered 23 MAGs belonging to Thermoproteota from VLS, and each of them contained 24–102 (average 65) genes involved in MMP. Although mcr genes were previously found from MAGs of the uncultivated Thermoproteota lineages (e.g., Evans et al., 2015. Science; Qu et al., 2022. ISME J), they were not detected on any of the above 23 MAGs. However, these VLS MAGs contained some other MM genes that can impact the key steps of methanogenesis and/or reverse methanogenesis (i.e., anaerobic methane oxidation), such as the genes encoding heterodisulfide reductase (hdrA1, hdrA2, hdrB2, hdrC2, and hdrD), methenyltetrahydromethanopterin cyclohydrolase (mch), methylenetetrahydromethanopterin dehydrogenase (mtd, mtdB), 5,10-methylenetetrahydromethanopterin reductase (mer), tetrahydromethanopterin S-methyltransferase (mtrA and mtrH), dimethylamine corrinoid protein (mtbC), methylamine corrinoid protein Co-methyltransferase (mtmB), and trimethylamine corrinoid protein (mttBC). These genomic data suggested that at least some of the VLS Thermoproteota members were able to participate in MM, though the available data is unable to determine whether they are methanogens or methanotrophs, or both.\nIn the revised version, we provided the gene annotations of the above 23 MAGs, along with the other 76 VLS MAGs, and highlighted the genes involved in MMP (see Table S13). In addition, we have revised the text as (Lines 281-287):\n“Some MAGs of Thermoproteota contained genes encoding for key steps of MM41-44, … In our data, we recovered 23 VLS MAGs belonging to the Thermoproteota, and each of them contained 24–102 (average 65) genes involved in MMP (Supplementary Table S13).”\nEvans PN, Parks DH, Chadwick GL, Robbins SJ, Orphan VJ, Golding SD & Tyson GW (2015) Methane metabolism in the archaeal phylum Bathyarchaeota revealed by genome-centric metagenomics. Science 350: 434-438. DOI: https://doi.org/10.1126/science.aac7745\nOu YF, Dong HP, McIlroy SJ, et al. (2022) Expanding the phylogenetic distribution of cytochrome b-containing methanogenic archaea sheds light on the evolution of methanogenesis. ISME J 16: 2373-2387. DOI: https://doi.org/10.1038/s41396-022-01281-0\nReviewer #2 (Remarks to the Author):\nMicroorganisms are main players in the cycling of the greenhouse gas methane, but recent discoveries of extrachromosomal elements associated with these microbes suggest that they influence the cycling of CH4, too. The authors Zhong et al. now analyze 982 publicly available metagenomes and generate 11 new metagenomes from a freshwater lake in Croatia, to identity viruses associated with methane metabolizing microorganisms. They find viruses that seem to infect methanogens and aerobic methanotrophs in the existing and their own datasets and then analyze the AMGs encoded on these viruses/phages. This is a generally very interesting topic and the workflows described seem well-chosen. Yet some of the data presented seems a bit disconnected and preliminary for publication. For example, it is not surprising that viral communities differ between sediment sites because the microorganisms differ between those sites. Thus it would be important to resolve the microbial community of the sites or at least present those microorganisms that are involved in methane cycling in these sites (and present this data), and then present the viral MAGs found in those sites. It is also not clear how novel the viral contigs are since they were only compared to cultivable viruses and a selection of metagenomic datasets. There is also no clear boundary between methanogenic archaea and methanotrophic bacteria. A reader that is no expert in this would likely be very confused as everything is presented together. This manuscript unfortunately requires a major revision before it can be published.\nResponse #10: We appreciate your kind evaluation of our research topic and workflows, and thank you for highlighting the weakness about the data presented to help us improve the work. We have substantially revised the manuscript based on the above comments, which should more deeply explain and contextualize the findings in this study. Please note that all the line numbers cited in our responses are based on the revised manuscript with tracked changes. The revisions include:\n(1) We now introduced Vrana Lake sediment (VLS) microbial communities/functions first, before we discussed VLS viral communities/functions. We also added the gene annotations of all the 99 MAGs obtained from VLS and highlighted all genes involved in MMP. Please see Table S13 and the below text:\n- Lines 238-247: “The cellular microbial communities, investigated based on relative abundances of the 99 bacterial/archaeal metagenome-assembled genomes (MAGs) recovered from VLS metagenomes (Supplementary Table S10; see Methods), were distinct between muddy and sandy sites (Supplementary Fig. S8), which had very different physicochemical conditions ... Similarly, the muddy and sandy sampling sites comprised mostly different viruses, with only 4.2% (131 of 3,146) of VLS vOTUs shared between sites (Fig. 3B). Ordination analysis, using the relative abundance data of vOTUs (Supplementary Table S11), confirmed that viral communities were significantly (p = 0.015) different between sites (Fig. 3C).”\n- Lines 267-275: “To explore the potential viral impacts on VLS ecosystems, we investigated virus-host linkages as reported previously (e.g., in soil and seawater11,39), via the iVirus tool VirMatcher40 that aggregates four different methods for host predictions (see Methods). Using the 99 VLS bacterial/archaeal MAGs as the host database (Supplementary Table S10; See Methods), we could link 2,167 of the 3,146 vOTUs (68.9%) to microbial hosts belonging to 17 different phyla (Fig. 4A; Supplementary Table S12). The VLS microbial communities were dominated by Thermoproteota (relative abundance: average 24.7% and range 12.7–41.4%; archaea) and Chloroflexi (average 23.5% and range 17.9–29.4%; bacteria) (Supplementary Fig. S10).”\n- Lines 294-302: “The 99 VLS MAGs contained a total of 5,503 genes (136 distinct genes after dereplication) involved in MMP, including genes that can impact the key steps of MM such as the genes encoding methylenetetrahydromethanopterin dehydrogenase (Mtd), 5,10-methylenetetrahydromethanopterin reductase (Mer), tetrahydromethanopterin S-methyltransferase (Mtr), and methyl-coenzyme M reductase (Mcr) (Supplementary Table S13). To assess if VLS viruses also encode MM AMGs and thus could be modulating the hosts’ MM, we annotated genes for all the VLS vOTUs (Supplementary Table S14) and screened them for putative virus-encoded AMGs, including MM AMGs.”\nWe would like to point out that a more detailed analyses of the microbial data for VLS samples will be used in a different manuscript that is in preparation. In the current paper, we specifically focused on viruses and only presented limited microbial data (e.g., MAG-based community, taxonomy, and genes involved in MMP) that helped reveal virus-microbe interactions in the sediments.\nApart from the revisions for the VLS datasets, we have also conducted additional analyses of microbial MM genes for the public metagenomic datasets where we searched for viral MM AMGs, as per the request of the other reviewer. Then we reorganized the text to introduce the microbial MM genes first, before we documented the viral MM genes in the environments where the public metagenomic datasets originated. Please refer to Response #3 for further details about this revision.\n(2) In the revised manuscript, we included both cultivated (NCBI RefSeq database) and uncultivated viral genomes (IMG/VR v4 database) for assessing viral novelty. Please see Response #14 for further details.\n(3) We have revised the text to clarify whether we refer to methanogenic archaea, methanotrophic bacteria, or anaerobic methanotrophic archaea:\nLines 72-75: “Recently, some viruses have been found to encode pmoC and cofF as auxiliary metabolic genes (AMGs), with the potential to supplement the aerobic oxidation of methane by their bacterial host in freshwater lakes7.”\nLines 86-88: “including those participating in methane production (i.e., methanogenesis by archaea) and oxidation (either by aerobic methanotrophic bacteria or anaerobic methanotrophic archaea).”\nLines 146-152: “Functionally, the pmoC gene participates in the aerobic methane oxidation pathway of bacterial methanotrophs27, while the other six genes mtrA, fwdF, fae, cofE, cofF, and frhB are involved in the pathways of methanogenesis and/or anaerobic oxidation of methane (AOM)28 (Fig. 1A & Supplementary Fig. S1). The pmoC (methane monooxygenase subunit C) gene encodes a subunit of the particulate methane monooxygenase (pMMO) that catalyzes the aerobic oxidization of methane to methanol in bacteria (Fig. 1A & Supplementary Fig. S1)27.”\n(4) We have carefully improved all the main figures. Please see Response #21 for further details.\nWe hope that these revisions have better contextualized this work for a publication.\nMajor comments:\n- Headers of sections are imprecise. E.g. “Viruses encode 24 AMGs modulating…” suggests this is a general finding. It would be better to rephrase to “Some viruses encode genes that could augment/modulate microbial methane metabolism”. Generally, abbreviations should be avoided in results sections because the reader will be lost..\nResponse #11: Thank you for these comments. We have changed this section header to “Some viruses encode genes that could modulate microbial methane metabolism” based on your suggestion (Lines 94-95). In addition, we modified the below four section headers:\nLine 204: Modified from “Vrana lake sediment comprises mostly novel viral genera and species” to “Vrana Lake sediment comprises mostly novel viral genera”.\nLine 237: Modified from “VLS viral communities differ between sediment sites” to “Viral communities differ between sediment sites and across sediment depths”.\nLines 265-266: Modified from “VLS viruses impact ecosystem via infecting dominant microbes” to “Abundant viruses likely infect dominant microbes of the Thermoproteota and Chloroflexi to impact the sediment ecosystems”.\nLine 292-293: Modified from “VLS viruses encode AMGs to modulate hosts’ metabolism” to “Some Vrana Lake sediment viruses encode genes to modulate host metabolisms”.\nWe removed all the abbreviations in the results headers, but kept the abbreviations in the results sections as we respectively disagree with your opinion that “abbreviations should be avoided in results sections because the reader will be lost”.\n- L: 100-110: The aerobic methane oxidation should be more clearly separated form (reverse) methanogenesis, since these are very distinct processes mediated by completely different microbes\nResponse #12: Thank you for pointing out this omission. We have modified the text to clearly clarify that the pmoC gene is involved in aerobic oxidation of methane in bacterial methanotrophs (Lines 146-152):\n“Functionally, the pmoC gene participates in the aerobic methane oxidation pathway of bacterial methanotrophs27, while the other six genes mtrA, fwdF, fae, cofE, cofF, and frhB are involved in the pathways of methanogenesis and/or anaerobic oxidation of methane (AOM)28 (Fig. 1A & Supplementary Fig. S1). The pmoC (methane monooxygenase subunit C) gene encodes a subunit of the particulate methane monooxygenase (pMMO) that catalyzes the aerobic oxidization of methane to methanol in bacteria (Fig. 1A & Supplementary Fig. S1)27”.\n- L111f: what is cool about the finding that mtrA is encoded on a viral contig is that the Mtr complex is directly involved in energy conservation. This part could be really improved by diving in a bit deeper into the function and stating this!\nResponse #13: Thanks for your suggestion to enhance our work. We agree and have modified the text to document that the enzyme complex Mtr catalyzes an energy conserving step in the MMP (Lines 152-158):\n“In the methanogenetic pathway, the mtrA (tetrahydromethanopterin S-methyltransferase subunit A) gene encodes a subunit of the membrane-associated multienzyme complex Mtr that transfers the methyl group of N5-methyltetrahydromethanopterin to coenzyme M (CoM) and produces Methyl-CoM29, which is an exergonic (ΔG°′ = −29 kJ/mole), sodium-ion-translocating step contributing to ion motive force in the methanogens’ energy metabolism. This energy conservation mechanism happens in all methanogens being able to produce methane from CO2 or acetate30.”\n- L169: why were the vOTUs only compared to cultivable viruses? This of course will lead to the (likely false) inferral that they are all novel… There is a huge database of viruses/phages derived from metagenomics (and were thus not cultivated) and the viral contigs should be compared to those? The authors then do compare them to “selected published environmental metagenomes”(384) but how and why were these selected?\nResponse #14: We agree that the viruses should be compared to both cultivable and uncultivable viral genomes in the databases when evaluating their novelty. In the revised version, we removed the analyses of comparing to the “selected published environmental metagenomes” (384). Instead, we included the viruses from both NCBI RefSeq and IMG/VR databases for taxonomic assignments. The later database contained uncultivated viral genomes from metagenomics. As a result, 139 vOTUs (instead of 3 in the original analyses that only included cultivable viruses for comparison) were assigned a taxonomy. We have modified the main text to reflect these changes:\nLines 221-226: “Taxonomic analyses, by comparing VLS viruses to viral genomes in both the NCBI RefSeq and IMG/VR databases (see Methods), revealed that most of the VLS long vOTUs (911 of 1,050) could not be taxonomically classified, indicating a high degree of novelty among VLS viruses. The remaining 139 vOTUs were assigned to Caudoviricetes, Faserviricetes, and Megaviricetes (Fig. 3A; Supplementary Fig. S7 & Table S9).”\nLines 486-491: “Because viruses lack any single, universally shared gene, we established taxonomy using gene-sharing network analysis from viral sequences ≥10 kb in length using vConTACT v223. Briefly, this analysis compared the 1,050 “long” VLS vOTUs and the 911 public datasets-originated viruses that contained MM AMGs to viral genomes in the National Center for Biotechnology Information (NCBI) RefSeq database (release v201) and the IMG/VR v4 database, and generated viral clusters approximately equivalent to known viral genera11,23,73.”\nIn addition, Fig. 3A was modified to accommodate the changes in taxonomy assignments, and provided below for review convenience:\n- L:254: how do the viruses alter the ecophysiology of the sediments? That seems like a really huge claim with no evidence to support this\nResponse #15: We apologize for this overstated claim, as our data did not sufficiently support it. This sentence has been modified as (Lines 325-329):\n“While no AMGs that can directly participate in MMP were detected from these methane-rich lake sediments, these results indicate that VLS viruses encode functional AMGs that likely alter microbial metabolisms in the Varana Lake sediments, including an AMG that have an indirect influence on MM through manipulating a putative methanogen’s iron metabolism.”\n- L259-275: This reads like a summary or even abstract and not like a discussion.\nResponse #16: We agree and have shortened this paragraph to ¼ length in the revised version. Now it reads as (Lines 332-354):\n“After carbon dioxide, methane is the second largest contributor to warming, accounting for approximately 20% of greenhouse gas-driven warming3,56,57. While it is widely accepted that bacteria and archaea are major players in the global methane cycling, little was known about how viruses might impact MM. This study identified 24 virus-encoded MM AMGs in 911 viral contigs by analyzing 982 published metagenomes from environments where microbial MM is known to occur, and where microbial genes involved in MMP were detected. We found that the extent to which viruses use MM AMGs to modulate host MMP may vary depending on the ecological properties of the habitat in which they dwell. This finding is consistent with previous reports of the habitat-specific association of AMGs in the environments13,59.”\n- L284-287: where does this hypothesis come from and how does the data presented in this paper connect to it?\nResponse #17: We apologize for this mistake. We realized that this is an overspeculation and our data presented cannot mechanistically support this hypothesis. In the revised manuscript, we removed this hypothesis and added the suggestion that future studies are necessary to test the presence pattern of MM AMGs and its mechanism in the environments (Lines 396-399):\n“Future studies are necessary to experimentally validate … and to further test the presence pattern of MM AMGs and its mechanism as more host-associated and environmental metagenomes become available.”\n- Figure 1A: The representation is extremely hard to follow. It also suggests that pmoC is in one cell (virus ?) with the other enzymes? Gene names are to be written in italics. This Figure should be revised.\nResponse #18: We apologize for the confusion, and thanks for your comments to help us improve the figure. The pmoC should not be visualized in one cell with the other enzymes. We have modified this plot by removing the cell cartoon, italicizing the gene names, and reorganizing the representation of MMP.\nFor review convenience, we provide the revised plot and the figure legend below:\nFigure 1. Characterization of exclusive MM AMGs. (A) Schematic for viral participations in key MMP steps via encoding seven AMGs that exclusively participate in MMP. Viruses encoded seven AMGs (fwdF, fae, frhB, cofE, cofF, mtrA and pmoC; as colored in purple text) to impact the key steps in both methane production and oxidation. The methanogenesis pathway from CO to methane is indicated by orange arrows. More information for 17 additional AMGs that could potentially participate in both MMP and other types of metabolism pathways is provided in Supplementary Fig. S1, Fig. S2, & Table S5. (B) Genome maps … MMP, methane metabolism pathway.\n- Figure 1C: Is it necessary to show this tree in such a collapsed way? It seems to provide little information. As a sidenote: protein names should be capitalized.\nResponse #19: Thanks for your comments. We agree that this tree should not be visualized in such a collapsed way. It has been revised by removing some collapses and modifying the protein name “MtrA” as capitalized.\nFor review convenience, the updated figure is provided below:\n- Figure 2: Is it necessary to show this as a main figure? I am more confused by looking at it than having read the section in the text, unfortunately.\nResponse #20: We apologize for the confusion. We prefer to keep it as a main figure, as these Venn plots visualize the variation of MM AMGs’s presence/absence by habitat, with an emphasis on the finding that rumen, the host-associated habitat, contained the majority of virus-encoded MM AMGs by comparing to the environmental habitats. We choose Venn plots as they can easily visualize both the total number of AMGs in each environment (this is described in the text) and the number of AMGs shared between/among environments (not described in the text).\nTo address a comment on “There are lots of pie charts and circles in several figures” in Response #21, we have tried using UpSet plots, instead of the Venn plots that contained lots of circles, to visualize the AMG numbers:\nHowever, we found that the UpSet plots did not convey the findings better than Venn plots. Thus, we decided to use the Venn plots, but not the UpSet plots. In the revised version, we have modified the Venn plots and the figure legend to clearly convey the habitat-associated AMG distribution pattern as mentioned above.\nFor review convenience, we have included the revised figure and its legend below:\nFigure 2. Venn diagrams visualizing the number of MM AMGs from different habitats. (A) All the 24 MM AMGs. We identified 24 distinct MM AMGs from six habitats: rumen (16 AMGs), marine water (7), marine sediment (5), lake water (3), lake sediment (1), and hot spring sediment (2). Seven of these genes were identified in 2–4 habitats, and the remaining 17 were found exclusively in one of these habitats. (B) Ten AMGs involved in methanogenesis pathway. Of the 24 MM AMGs, 10 genes (i.e., mtrA, fwdF, cofE, cofF, frhB, ackA, pta, cooS, glyA, and fae) can directly participate in or synthesize an intermediate for the pathway of methanogenesis from CO or acetate (Supplementary Fig. S1). Nine of these 10 AMGs were found in rumen, while only one to three were found from other detectable environmental habitats including marine water, marine sediment, lake water, and lake sediment. Three of these genes were identified in 2–4 habitats, and the remaining six and one were found exclusively in rumen and marine sediment, respectively.\n- The Figures generally look a bit unpolished… e.g., Figure 3 (A) and (B) are outside the boxes, where (C) and (D) are inside. There are lots of pie charts and circles in several figures.. Figure 4A why does the green arrow restate in green font what is written in the left pie-chart? Were all 2,167 vOTUs assigned to either bacteria or archaea? This seems suprising that taxonomic classification/linkage is possible for every vOTU found..\nResponse #21: Thank you for your valuable comments to help us improve the figures. We have made significant revisions to all four main figures, including the modifications you highlighted in the comments. Among the revisions, while we removed one of the pie chart plots (the original Fig. 3B), we have retained the other pie charts and Venn plots (i.e., circles), as we consider them essential and appropriate for conveying key findings. We believe that the new figures represent a considerable improvement.\nFor the comment “Were all 2,167 vOTUs assigned to either bacteria or archaea? This seems suprising that taxonomic classification/linkage is possible for every vOTU found”, we apologize for the confusion. We identified a total of 3,146 vOTUs, and 2,167 of them (not all of them) could be linked to their hosts. They were all predicted to infect either bacteria or archaea, as this study only identified bacterial and archaeal viruses and used VLS bacterial/archaeal genomes as the host database for host prediction. To clarify this, we modified the figure (i.e., Fig. 4A, see below), the main text (Lines 89-92 and 270-274), and figure legend (Lines 916-918):\nLines 88-92: “We then generated an additional 11 metagenomes using … sample bacterial/archaeal viruses and investigate their potential impacts on MM during infection.”\nLines 269-273: “Using the 99 VLS bacterial/archaeal MAGs as the host database (Supplementary Table S10; See Methods), we could link 2,167 of the 3,146 vOTUs (68.9%) to microbial hosts belonging to 17 different phyla (Fig. 4A; Supplementary Table S12).”\nLines 946-948: “Of the 3,146 vOTUs, 2,167 were linked to putative hosts, with 1,110 and 1,057 vOTUs infecting bacteria and archaea, respectively.”\nFor review convenience, the updated version of Fig. 1, Fig. 3, and Fig. 4 are included below (Fig. 2 was provided in Response #20):\nFigure 1. Characterization of exclusive MM AMGs. (A) Schematic for viral participations in key MMP steps via encoding seven AMGs that exclusively participate in MMP. Viruses encoded seven AMGs (fwdF, fae, frhB, cofE, cofF, mtrA and pmoC; as colored in purple text) to impact the key steps in both methane production and oxidation. The methanogenesis pathway from CO to methane is indicated by orange arrows. More information for 17 additional AMGs that could potentially participate in both MMP and other types of metabolism pathways is provided in Supplementary Fig. S1, Fig. S2, & Table S5. (B) Genome maps of three viral contigs carrying the AMG mtrA gene. The three viral contigs belonged to the same viral population (with 97.4–97.8% genomic identities among each other) and carried an identical mtrA gene. CheckV was used to assess host-virus boundaries and remove potential host fractions on the viral contig. Genes were marked by five colors to illustrate AMGs (purple), phage genes (orange), phage hallmark genes (blue), potential cellular genes (green), and hypothetical protein genes (grey). (C) Phylogenetic tree of the viral and microbial mtrA genes. The tree was inferred using maximum likelihood method with protein sequences. Parametric bootstrap values (expressed as percentages of 1,000 replications) are shown at branching points. The scale bar indicates a distance of 0.1 substitutions per position in the alignment. The viral and microbial MtrA sequences are indicated in red and black, respectively. The numbers in parentheses indicate the number of protein sequences assigned to each group. The full phylogenetic tree (without collapsed groups) is provided in Supplementary Fig. S5A. The genomic maps and phylogenetic trees for the other six exclusive MM AMGs (pmoC, fwdF, fae, cofE, cofF, and frhB) are provided in Supplementary Fig. S3 and Fig. S5B–G. MM, methane metabolism; MMP, methane metabolism pathway.\nFigure 3. Viral communities of Vrana Lake sediments (VLS). (A) Taxonomic assignments of VLS vOTUs. The left chart shows the fraction of “long” vOTUs (length ≥10 kb) among all VLS vOTUs (n = 3,146). The right chart shows the taxonomy of VLS “long” vOTUs, when compared to viral genomes in the NCBI RefSeq and IMG/VR databases. Further details of taxonomic results are listed in Supplementary Table S9. (B) Shared and unique vOTUs between the two sediment sites (Muddy vs. Sandy sites) in Vrana Lake as shown in Supplementary Fig. S6. Only 4.2% of the 3,146 vOTUs were presented in both sites, while the remaining 95.8% were unique in either site. (C) PCoA plot of VLS samples based on the relative abundances of vOTUs. Samples are marked by the two sediment sites (Muddy and Sandy sites in green and orange, respectively) and the two metagenome types (bulk metagenome and virome as triangles and circles, respectively). PERMANOVA (permutations = 999) shows a statistically significant difference in viral communities between both sampling sites and metagenome types. The p values <0.05 are indicated in red.\nFigure 4. VLS virus-host interactions. (A) Predicted phylum-level hosts of VLS viruses. Of the 3,146 vOTUs, 2,167 were linked to putative hosts, with 1,110 and 1,057 vOTUs infecting bacteria and archaea, respectively. The number of vOTUs putatively infecting each phylum is indicated in parentheses after the phylum name. (B) Relative abundances of the two most abundant VLS microbial phyla and their predicted viruses. The two phyla are Thermoproteota (left panel of box plots) and Chloroflexi (right panel of box plots). Relative abundances of microbes and vOTUs were obtained based on their coverages generated by read mapping to MAG populations and vOTUs. Box-plot elements: center line, median; x symbol, mean; box range, upper and lower quartiles; whiskers, data range. (C) Genome map of the vOTU S225_M_175_17980bp encoding the AMG bfr. The genome content was characterized by the same methods described in Fig. 1B. (D) Phylogenetic tree of the viral and microbial bfr genes. The tree was inferred using maximum likelihood method with protein sequences (see Methods). Parametric bootstrap values (expressed as percentages of 1,000 replications) ≥40 are shown at branching points. The scale bar indicates a distance of 0.2 substitutions per position in the alignment. The viral Bfr sequence is indicated in red and other sequences are indicated in black. The full phylogenetic tree (without collapsed groups) is provided in Supplementary Fig. S11. VLS, Vrana Lake sediment.\nMinor comments:\n- L28: “methanogenic sediments” seems very colloquial, maybe rephrase to “sediments where methane emissions/consumption has been detected”(or whatever is applicable)\nResponse #22: The text has been rephrased based on your suggestion (Lines 88-90):\n“We then generated an additional 11 metagenomes using lake sediments, in which methane emission has been detected, from Croatia’s largest freshwater lake …”\n- L90: resulting in the discovery..\nResponse #23: It was modified as “resulting in the discovery” in the revised main text (Line 115).\n- L93L if only one protein sequence was investigated, how can the “conserved domains” be investigated ?\nResponse #24: We apologize for the confusion. The conserved domains were analyzed by comparing to the sequences in the domain database. To clarify, we modified the method sentence in the main text as (Lines 548-550):\n“In addition, the AMGs of the selected viral contigs were further used as examples for analyzing conserved domains by comparing to the domains in the public Conserved Domain Database (CDD) via NCBI CD-Search84 and …”\n- L104: on only one viral contig (not “from”)\nResponse #25: Thank you for identifying this mistake. We changed “from viral contigs” to “on viral contigs” through the text:\nLines 143-144: “Among these seven AMGs, fwdF and fae were each detected on only one viral contig, while the others were identified on 3 to 25 viral contigs …”.\nLines 538-539: “a discovery of 24 distinct AMGs (on a total of 911 viral contigs containing ≥1 MM AMG) …”.\nLine 546: “the AMG was presented on more than one contig”.\nLine 562: “if an AMG was identified on more than two viral contigs …”.\n- L106: it should be stated that pmoC is part of the aerobic methane oxidation pw of bacteria\nResponse #26: This was stated in the revised text (Lines 146-152):\n“Functionally, the pmoC gene participates in the aerobic methane oxidation pathway of bacterial methanotrophs27, while the other six genes mtrA, fwdF, fae, cofE, cofF, and frhB are involved in the pathways of methanogenesis and/or anaerobic oxidation of methane (AOM)28 (Fig. 1A & Supplementary Fig. S1). The pmoC (methane monooxygenase subunit C) gene encodes a subunit of the particulate methane monooxygenase (pMMO) that catalyzes the aerobic oxidization of methane to methanol in bacteria (Fig. 1A & Supplementary Fig. S1) 27.”\n- L107: how do the authors know that these genes are not involved in AOM ? Unless there is evidence for this (e.g., the metagenomic samples these viral contigs originate from, only comprise methanogens, and not ANME). If this was not investigated, it should be added “pathway of methanogenesis or reverse methanogenesis”\nResponse #27: Thank you for pointing out this mistake. There is no evidence to support that these genes are not involved in AOM. We have revised text as (Lines 147-149):\n“while the other six genes mtrA, fwdF, fae, cofE, cofF, and frhB are involved in the pathways of methanogenesis and/or anaerobic oxidation of methane (AOM)28 (Fig. 1A & Supplementary Fig. S1).”\n- L168: how were the vOTUs established?\nResponse #28: Viral contigs were clustered into vOTUs if they shared ≥95% nucleotide identity across 80% of their length. We have modified this sentence to include the methods (Lines 216-219):\n“We recovered 3,260 viral contigs from these metagenomes. These contigs were clustered into vOTUs if they shared ≥95% nucleotide identity across 80% of their lengths38, resulting in 3,146 vOTUs (≥5 kb), including 1,050 “long” (≥10 kb) vOTUs (Supplementary Table S8).”\nThe method was also described in the Methods section of the original main text (new Lines 467-469):\n“The remaining contigs were clustered into vOTUs (~species-level taxonomic unit) if they shared ≥95% nucleotide identity across 80% of their lengths as described previously38”\n- L 211: how were vOTUs linked to microbial hosts? This should be stated in brief in the results.\nResponse #29: We have modified this sentence to briefly describe how the vOTUs were linked to microbial hosts (Lines 267-273):\n“To explore the potential viral impacts on VLS ecosystems, we investigated virus-host linkages as reported previously (e.g., in soil and seawater11,39), via the iVirus tool VirMatcher40 that aggregates four different methods for host predictions (see Methods). Using the 99 VLS bacterial/archaeal MAGs as the host database (Supplementary Table S10; See Methods), we could link 2,167 of the 3,146 vOTUs (68.9%) to microbial hosts belonging to 17 different phyla (Fig. 4A; Supplementary Table S12).”\n- L221: I don’t think it is accurate to write “many” Thermoproteota members are methanogens. Also not sure why “archaeal” is required at this point, because all methanogens (we know so far!) are archaeal (this would fit in the beginning, but not here)\nResponse #30: We agree and have modified this sentence as (Lines 281-283):\n“Some MAGs of Thermoproteota contained genes encoding for key steps of MM41-44 …”\nWe also modified the text in the introduction to state that microbial methane production is mediated solely by archaea (Lines 62-64):\n“Methane cycling is largely driven by microbes, with microbial methanogenesis (all mediated by archaea) producing ~69% of the total methane released to the atmosphere6.”\nPer the other reviewer’s suggestion, we also provided gene annotations for all the Thermoproteota MAGs (n = 23) we have recovered from VLS metagenomes and highlighted their genes involved in MMP. Each of the Thermoproteota MAGs contained 24–102 (average 65) genes involved in MMP. Please see Response #9 for further details.\n- L224: had an important\nResponse #31: The phrase “had important” was modified to “had an important” in the revised text (Line 289).\n- L223-5: This seems pretty disconnected. What is the evidence for this or why bring this up here?\nResponse #32: Thank you for bringing this to our attention. We agree and have removed the text “Based on the classical Kill-the-Winner model” from this sentence. The modified text reads as (Lines 287-289):\n“Overall, these findings showed that the VLS viruses infected dominant microbial phyla, including ones involved in MM, and thus likely had an important impact on the sediment ecosystems.”\n- L242: “conserved functional domains of their corresponding enzymes” What does this mean?\nResponse #33: We apologize for the confusion. Our intention was to convey that all the AMGs possessed conserved functional domains. We have modified the text as (Lines 319-320):\n“All of them had the conserved functional domains (Supplementary Table S15)”\n- L:251: Why was this pN/pS analysis done?\nResponse #34: We calculated the pN/pS values for these AMGs, because the published standards methodology for studying virus-encoded AMGs (Pratama et al., 2021. PeerJ) recommends investigating gene’s selection pressure (e.g., pN/pS value) which indirectly indicates AMG’s evolution and functionality information.\nPratama AA, Bolduc B, Zayed AA, Zhong ZP, Guo J, Vik DR, Gazitua MC, Wainaina JM, Roux S & Sullivan MB (2021) Expanding standards in viromics: in silico evaluation of dsDNA viral genome identification, classification, and auxiliary metabolic gene curation. PeerJ 9: e11447. DOI: https://doi.org/10.7717/peerj.11447.\n- L282: “significant” may not be the best term here\nResponse #35: We have removed these sentences in the revised version based on another comment (see Response #16).\n- L293: CO is also a substrate for methanogens\nResponse #36: Thanks for identifying this omission. We have modified the sentence to clarify that CO is also a substrate for methanogens (Lines 365-370):\n“For example, the carbon monoxide (CO) dehydrogenase gene (cooS; identified from a rumen metagenome; Supplementary Table S5) catalyzes the oxidation of CO to CO 60, which is the substrate of a methanogenesis pathway from CO in ruminants61. In addition, the oxidation of CO to CO in itself could be a step of methanogenesis using CO as the substrate62 and an energy generating metabolic reaction63.”\n- L:294: It is a bit odd to have (n=17) in the discussion section as such specific values usually belong to the results\nResponse #37: The phrase “(n=17)” was removed from the sentence, now it reads as (Lines 371-373):\n“Thus, while many of the identified AMGs have the potential to participate in MM, without further verification, their actual functions remain hypothetical.”\n- L297ff: it should be stated that mcr genes were recently found in “Borgs” – archaeal ECEs that are neither typical viruses nor plasmids\nResponse #38: We agree and have added a statement to the text (Lines 376-379):\n“Interestingly, while not yet found in viral genomes, the mcr genes were recently discovered in a novel group of extrachromosomal elements called “Borgs”, that are associated with ANME in the genus Methanoperedens15”"}, {"round": 2, "reviewer_comments": {"reviewer_1": "This manuscript has improved substantially during revisions, and the authors have responded to most of the points raised prior. I have some comments that I still think could improve the MS further:\n1. There should be a figure describing the predicted iPhoP hosts of the viruses with MM AMGs. This could be added to Figure 2.\n2. There is no evidence here that the viruses actually do modulate microbial methane production. I would encourage the authors to reframe their title in terms of their primary finding, which is that viral genomes encode MM AMGs.\n3. Throughout the manuscript \"KEGG annotations\" are referred to without properly specifying the methods used to annotate KEGG orthologs. Were HMMs used? KOFams or METABOLIC? What cutoffs were used? Please be specific and cite, as all annotations are not equal. Apologies if I missed it.\n4. The authors still need to put their own data in the context of their overall analysis in the discussion. On average in their global analysis, what % of metagenomics from these environments had at least one viral MM AMG? With that number in mind, is it unusual or expected that they did not recover one from their own lake samples?\nMinor points:\n1. iPhoP is not mentioned in the methods - please add it.\n2. The assembler used should be stated directly in the methods.\n3. In the methods there is a dn/ds analysis described that I do not think appears in the text.\n4. I think it is poor communication to state that Methane is a \"more potent\" greenhouse gas than CO2, as this can be misleading. I think the more appropriate statistic is the cited \"After carbon dioxide, methane is the second largest contributor to warming, accounting for approximately 20% of greenhouse gas-driven warming\" in the discussion, and this fact should be used in the intro instead of the \"potent\" comparison.\n5. Please state how deeply samples were sequenced (Gbp) in the methods.", "reviewer_2": "Upon careful review of the revised manuscript by Zhong et al., I am pleased to note that the authors have adequately addressed my previous concerns and suggestions. The revisions have significantly improved the clarity and quality of the work. Consequently, I can now recommend this manuscript for publication. I believe it makes a valuable contribution to the field and will be of interest to the readership of Nature Communications."}, "author_response": "Response #1: We appreciate your kind words in assessing our revisions and thank you for the additional comments, which we have implemented into this revision as described below. Please note that all the line numbers cited in our responses are based on the revised manuscript with tracked changes.\n1. There should be a figure describing the predicted iPhoP hosts of the viruses with MM AMGs. This could be added to Figure 2.\nResponse #2: We agree and have modified Figure 2 to include the information on host prediction as Fig. 2A. The original Fig. 2A and 2B were changed to Fig. 2B and 2C, respectively.\nFor review convenience, we provide Fig. 2 and the figure legend below:\nFigure 2. Predicted hosts of viruses encoding MM AMGs and habitat association of MM AMGs. (A) Phylum-level host predictions of 257 viruses that encoded MM AMGs. Of the 911 viral contigs encoding MM AMGs, 257 were successfully linked to hosts that belonged to four archaeal and 10 bacterial phyla. Additional information about the predicted hosts is provided in Supplementary Table S4. (B–C) Habitat association of all the 24 MM AMGs (B) and the 10 MM AMGs involved in methanogenesis pathway (C). We identified 24 distinct MM AMGs from six habitats: rumen (16 AMGs), marine water (7), marine sediment (5), lake water (3), lake sediment (1), and hot spring sediment (2). Seven of these genes were identified in 2–4 habitats, and the remaining 17 were found exclusively in one of these habitats. Of the 24 MM AMGs, 10 genes (i.e., mtrA, fwdF, cofE, cofF, frhB, ackA, pta, cooS, glyA, and fae) can directly participate in or synthesize an intermediate for the pathway of methanogenesis from CO2 or acetate (Supplementary Fig. S1). Nine of these 10 AMGs were found in rumen, while only one to three were found from other detectable environmental habitats including marine water, marine sediment, lake water, and lake sediment. Three of these genes were identified in 2–4 habitats, and the remaining six and one were found exclusively in rumen and marine sediment, respectively.\n2. There is no evidence here that the viruses actually do modulate microbial methane production. I would encourage the authors to reframe their title in terms of their primary finding, which is that viral genomes encode MM AMGs.\nResponse #3: We have modified the title from “Viral modulation of microbial methane production and oxidation varies by habitat” to “Viral potential to modulate microbial methane metabolism varies by habitat”.\n3. Throughout the manuscript \"KEGG annotations\" are referred to without properly specifying the methods used to annotate KEGG orthologs. Were HMMs used? KOFams or METABOLIC? What cutoffs were used? Please be specific and cite, as all annotations are not equal. Apologies if I missed it.\nResponse #4: We apologize for this omission. We used published tools to annotate genes against the KEGG database. The viral and microbial contigs were annotated by the tools VIBRIANT and DRAM, respectively (Kieft et al., 2020 Microbiome; Shaffer et al., 2020 Nucleic Acids Res), using their default parameters. We added the below sentences to the Methods section to clarify the gene annotation methods:\nLines 484-486: “Specifically, once viral contigs were recovered from metagenomes, they were processed with VIBRANT to obtain gene functional annotations against the KEGG and PFAM databases and identify putative AMGs by the default parameters22.”\nLines 514-517: “The assembled contigs, excluding viral contigs, were annotated by DRAM against the KEGG and PFAM databases by the default parameters and further used for recovering microbial MM genes based on their KEGG and PFAM annotations (with consistent annotation in the two databases) and the KEGG MM pathway modules18.”\nKieft K, Zhou Z & Anantharaman K (2020) VIBRANT: automated recovery, annotation and curation of microbial viruses, and evaluation of viral community function from genomic sequences. Microbiome 8: 90.\nShaffer M, Borton MA, McGivern BB, Zayed AA, La Rosa SL, Solden LM, Liu P, Narrowe AB, Rodriguez-Ramos J, Bolduc B, Gazitua MC, Daly RA, Smith GJ, Vik DR, Pope PB, Sullivan MB, Roux S, & Wrighton KC (2020) DRAM for distilling microbial metabolism to automate the curation of microbiome function. Nucleic Acids Res 48: 8883-8900.\n4. The authors still need to put their own data in the context of their overall analysis in the discussion. On average in their global analysis, what % of metagenomics from these environments had at least one viral MM AMG? With that number in mind, is it unusual or expected that they did not recover one from their own lake samples?\nResponse #5: On average, 32.2% of the metagenomes in the global analysis had at least one viral MM AMGs, whereas only 1.7% of the lake-sediment metagenomes contained at least one viral MM AMGs. Thus, the finding that we did not recover one MM AMG from the 11 Vrana Lake sediment metagenomes (our own data) would agree with the analyses of the public data. The below sentences, including one in the Discussion section, have been modified to contextualize the above arguments:\nLines 118-119: “These 911 viral contigs originated from ~32% (316 of 982) of the here analyzed metagenomes (Supplementary Table S3).”\nLines 283-288: “Interestingly, none of these AMGs were predicted to be directly involved in the MMP, which would agree with the inference of our analysis of publicly available metagenomes that the extent to which viruses modulate hosts’ MM may vary by habitats, and that MM AMGs seem to be less common in environmental habitats including lake sediments (<2% of the publicly available lake-sediment metagenomes had ≥1 MM AMG).”\nLines 318-324 (in the Discussion section): “We found that the extent to which viruses use MM AMGs to modulate host MMP may vary depending on the ecological properties of the habitat in which they dwell. Specifically in lake sediments, less than 2% of the publicly available metagenomes contained ≥1 MM AMG and no MM AMG was identified from the 11 metagenomes of Vrana Lake sediments, in which methane emission has been detected. This finding is consistent with previous reports of the habitat-specific association of AMGs in the environments15,57.”\nMinor points:\n1. iPhoP is not mentioned in the methods - please add it.\nResponse #6: We apologize for this mistake. We have added it to the Methods section of the revised manuscript (Lines 501-504):\n“Hosts of the 911 viral contigs containing ≥1 MM AMG were predicted by iPHoP26 (confidence score ≥90%), resulting successfully virus-host linkages for 257 viral contigs (Supplementary Table S4).”\n2. The assembler used should be stated directly in the methods.\nResponse #7: We added the assembler information to the Methods section in the revised version (Lines 415-420).\n“Sequencing reads were filtered for quality … Then the metagenomic sequence data was assembled to contigs by metaSPAdes67, using a previously established pipeline for assembling pre-amplified metagenomes (parameters: read deduplication + read error correction + --sc + -k 21,33,55,77,99,127)68.”\n3. In the methods there is a dn/ds analysis described that I do not think appears in the text.\nResponse #8: Thanks for pointing out this omission. The dS/dS result was described in the Supplementary Information in the last version. In the revised version, we also included this result in the main text (Lines 295-297):\n“Evolutionary pressure assessments within species and across lineages found that this virus-encoded Bfr was likely functional and under purification selection (pN/pS = 0; average dN/dS = 0.114; Supplementary Table S15 & Table S16).”\n4. I think it is poor communication to state that Methane is a \"more potent\" greenhouse gas than CO2, as this can be misleading. I think the more appropriate statistic is the cited \"After carbon dioxide, methane is the second largest contributor to warming, accounting for approximately 20% of greenhouse gas-driven warming\" in the discussion, and this fact should be used in the intro instead of the \"potent\" comparison.\nResponse #9: We have modified the sentence that has the phrase “more potent” in the introduction based on your suggestion. Now it reads as (Line 62-64):\n“Methane (CH ) is ranked second after carbon dioxide (CO ) in terms of the overall contribution to atmospheric warming and accounts for ~20% of the greenhouse gas-driven warming3-5.”\n5. Please state how deeply samples were sequenced (Gbp) in the methods.\nResponse #10: The sequencing depth of each sample was provided in Supplementary Table S8 in the last version, and now we also included this information in the Methods section (Lines 415-417):\n“Sequencing reads were filtered for quality … generating a total of 9.5 × 1010 bases of sequencing data (range 0.3–1.5 × 1010 bases, average 8.6 × 109 bases per library; Supplementary Table S8).”\nReviewer #2 (Remarks to the Author):\nUpon careful review of the revised manuscript by Zhong et al., I am pleased to note that the authors have adequately addressed my previous concerns and suggestions. The revisions have significantly improved the clarity and quality of the work. Consequently, I can now recommend this manuscript for publication. I believe it makes a valuable contribution to the field and will be of interest to the readership of Nature Communications.\nResponse #11: We appreciate your kind evaluations of our revisions and thank you for recommending our work for publication at Nature Communications."}]

Decision letter