Arctic Ocean virus communities and their seasonality, bipolarity, and prokaryotic associations

Domain: earth_science
Published: 2025-07-11
Editorial outcome: major_revision

Abstract
Viruses of microbes play important roles in ocean environments as agents of mortality and genetic transfer, influencing ecology, evolution and biogeochemistry. However, we know little about the diversity, seasonality, and host interactions of viruses in polar waters. Here, we study dsDNA viruses in the Arctic Fram Strait across four years via 47 long-read metagenomes of the cellular size-fraction. Among 5662 vOTUs, 98% and 2% are Caudoviricetes and Megaviricetes, respectively. Viral coverage is, on average, 5-fold higher than cellular coverage, and 8-fold higher in summer. Viral community composition shows annual peaks in similarity and strongly correlates with prokaryotic community composition. Using network analysis, we identify putative virus-host interactions and six ecological modules associated with distinct environmental conditions. The network reveals putative novel cyanophages with time-lagged correlations to their hosts (in late summer) as well as diverse viruses correlated with Flavobacteriaceae, Pelagibacteraceae, and Nitrosopumilaceae. Via global metagenomes, we find that 42% of Fram Strait vOTUs peak in abundance in high latitude regions of both hemispheres, and encode proteins with biochemical signatures of cold adaptation. Our study reveals a rich diversity of polar viruses with pronounced seasonality, providing a foundation for understanding viral regulation and ecosystem impacts in changing polar oceans.

Introduction
Polar regions are subject to the strongest seasonal cycles on Earth, and experience intense pressure from climate change1,2. The functioning of polar ecosystems is critical to biogeochemical cycles3–5, and is under marked ecological and evolutionary constraints6. Such ecosystem processes are strongly driven by microorganisms and their interactions, including viral dynamics. In addition to cell death, viruses also drive evolution via gene exchange, frequency-dependent selection7,8, and transmission and expression of auxiliary metabolic genes9–11. Therefore, characterizing the diversity and dynamics of viruses is paramount for understanding the function and stability of polar ocean ecosystems.

Relative to more accessible temperate, subtropical, and tropical oceans, few studies have examined virus diversity and ecology in polar waters12,13. It is known, however, that polar viromes are distinct from their warm-water counterparts13–22 and are typically dominated by bacteriophages with a high level of diversity14. In addition to spatial structuring, polar viral communities also shift over time, for example in the Antarctic, following bloom dynamics and assembling into distinct communities across seasons23. One year-round study showed strong seasonal variation among virus communities, with sharp decreases in virus-to-prokaryote ratios at the onset of the spring bloom and highest ratios in winter22. Others have reported increases during spring-summer, and highest abundances of viruses in winter24,25. Some evidence suggests that different lifestyles of viruses (i.e., lytic vs. lysogenic) may exhibit distinct dynamics across seasons, with lytic infection more prevalent in the spring bloom26.

However, due to the challenges of continuous sampling in the polar regions across multiple years, the degree of seasonality among polar viruses—here, meaning annually repeating patterns of populations and communities across the same seasons of different years27—and the potential ecological implications remain to be discovered. Such seasonality has been previously observed in marine temperate environments28–30. In the Arctic, time-series have been critical to advancing an understanding of the biological carbon pump31, elucidating benthopelagic coupling and biotic interactions32, both in the water column33 and on sinking particles34. Furthermore, continuous sampling, over the long-term, can discern the impact of ‘Arctic Atlantification'35–37: the northward expansion of subarctic habitats through the Fram Strait—the major connection between Atlantic and Arctic Oceans38. Here, polar water outflowing the central Arctic Ocean via the East Greenland Current (EGC) meets the West Spitsbergen Current (WSC), transporting warmer Atlantic water into the Arctic Ocean39–41. In the WSC, prokaryotic communities exhibit pronounced seasonality, underpinned by changes in photosynthetically active radiation (PAR) and mixed layer depth42.

Given this, and the tight coupling of viruses and their hosts9–11, we aim to examine the degree that viral populations are seasonally structured in the Arctic. The strong seasonal gradients in the Arctic provide an opportunity to observe host-virus dynamics and how they relate to prevailing conceptual models such as the Piggyback-the-winner43–45, Constant-Diversity46,47, and Red-Queen29,48 hypotheses. We examine the diversity and seasonality of dsDNA virus communities in the Arctic Fram Strait over four complete annual cycles at roughly monthly resolution. This time-series of virus communities in the Arctic demonstrates considerable seasonality in viral diversity, their association with environmental conditions and potential microbial hosts, and their distribution across the global ocean.

Results
We examined 47 long-read metagenomes from samples collected at near-monthly resolution over a four-year period (Sep 2016–Jul 2020) in the WSC (Fig. 1a, and Supplementary Data 1), from an average depth of 29 m. As samples originate from the cellular size-fraction (>0.2 µm), our study characterizes the diversity of actively infecting (i.e., intracellular) viruses, free-living viruses with a size of >0.2 µm, viruses attached to particles, and/or integrated viruses.

The total number of reads per sample were 196,489 ± 18,358, with an average length of 5435 ± 405 bp. Viral sequences were predicted based on both mapping to assembled contigs and raw reads (see Virus Prediction, Methods). The overall number of predicted viral sequences was concordant between the two approaches (Supplementary Data 1). On average, each sample harbored 9815 ± 965 viral reads based on raw read counts (5.3 ± 0.5%), and 8825 ± 1020 reads mapping to vOTUs (viral Operational Taxonomic Units) (6.8% ± 1.8%). Notably, despite not pre-filtering to remove eukaryotes, the large majority of cellular reads in the dataset were prokaryotic (90.2% ± 1.4%), with eukaryotic reads making up a smaller proportion (7.2% ± 1.0%). Hereafter, we focus on the contig-based predictions because of their improved genomic context.

Diversity of Fram Strait viruses over time
We first investigated how the community of Fram Strait viruses is structured over time. Through contig-based predictions, we identified 5662 vOTUs (95% nucleotide identity clusters, see Methods, Supplementary Data 2) which were derived from 7775 viral contigs that were greater than 10 Kb in length. These vOTUs predominantly represented the classes Caudoviricetes and Megaviricetes, with Caudoviricetes being the largest group, comprising 5531 vOTUs.

To assess the temporal dynamics of vOTUs, we normalized their abundance based on the estimated number of cells sequenced in each metagenome—a metric we term coverage-based Virus to Cell Ratio (cVCR). The cVCR approach not only accounts for differences in sequencing depth but also captures shifts in virus to host ratios across samples. Using this metric, we observed clear seasonal structuring in the cVCR of viral communities and taxa. Community cVCR values reached annual maxima of 20–35 and were four-fold higher during July–September (average cVCR 8.6) than during October–June (Fig. 1b). These patterns were consistent also when we sub-sampled the total reads to 100,000 reads per sample, as a test of the robustness of this metric to sequencing depth, as well as when we calculated coverage per Gigabase (Supplementary Fig. 1), thus we focus the remaining analyses on cVCR from the full dataset. This seasonal variation in cVCR was generally consistent across Caudoviricetes. Notably, the most abundant vOTU, a member of the Caudoviricetes, was also the most persistent, occurring in 45 of the 47 sampled time points. Similarly, the most persistent vOTUs overall, present in half of the samples, were predominantly members of the Caudoviricetes. In contrast, Megaviricetes were found less frequently, but were prevalent in late summer (July–September) (Fig. 1b, and Supplementary Fig. 1a), when eukaryotic phytoplankton, their presumed hosts, are most abundant49. Overall, lytic lifestyle was predicted to be predominant, with an average of 94.0 ± 0.4% of viruses, a pattern that was invariable across seasons (Fig. 1c). The Southern Ocean featured higher rates of lysogeny during periods of low production26. As for the latter point, it is important to point out that some lysogens may not be recognized due to missing integrases or other factors, and, vOTUs may not represent the full prophage, which could then be misinterpreted as lytic rather than integrated. Additionally, lytic viruses may outnumber lysogens in terms of copies within the cellular metagenome, making them easier to detect and assemble, which could bias the recovery or obscure lysogens.

We next investigated how the overall diversity of viruses changes across seasons. As viral diversity was not saturated at any sequencing depth (Fig. 2a), we compared the diversity of vOTUs across months after subsampling to a normalized depth. The diversity of viruses reached a maximum between late summer-autumn (Aug–Oct) (Fig. 2b, Supplementary Fig. 2), with lowest richness values typically in May. These patterns were consistent both by estimates from the slopes of rarefaction curves to the sub-sampled depth (Fig. 2b) as well as the extrapolated richness (Supplementary Fig. 2). Notably, prokaryotic community richness was also lowest in May, though the timing of highest prokaryotic richness was different than for viruses, with the highest prokaryotic richness observed in winter (vs. late summer-autumn for viruses)42 (Supplementary Data 1). Inter-annual variability in richness was also observed, for example, with May 2019 showing elevated richness relative to other years (Fig. 2c). The dynamics in richness calls for further sequencing efforts, with increased depth of sequencing to further unravel the viral diversity that persists, especially, through late polar night.

In addition to the variations in overall relative abundance (cVCR) and diversity, we also observed strong seasonality in the composition of the major bacteriophage groups (Fig. 3a). The composition of vOTUs within both Caudoviricetes and Megaviricetes showed a sinusoidal-like pattern over time, with peaks in Bray-Curtis similarity at 12, 24, and 36 month intervals (Fig. 3a). Furthermore, in some cases, similarity between samples from opposing seasons was zero, which is likely a combination of strong seasonality, low relative abundance (and hence detection) of vOTUs in winter, as well as variable sequencing depth. Notably, the seasonality observed in virus community composition was consistent regardless of minimum coverage value utilized for detected viral presence, though similarity was less at higher thresholds (Supplementary Fig. 3).

To further contextualize the dynamics at the vOTU level, we examined the underlying sequence diversity of vOTUs in two ways. First, we explicitly compared all non-clustered viral contigs >10,000 bp across years. Out of 7360 such contigs, we found 21 clusters of contigs, encompassing 55 contigs total that were 100% matches between at least two years. This suggests that detection of viruses with no genetic changes across years is rare, though there is a limitation of sequencing depth. Second, based on analysis of single nucleotide variations in reads to vOTUs (via InStrain, see Methods), we find that there tends to be microdiversity underlying the vOTU populations, though the amount varies between vOTUs ranging from likely insignificant and potentially related to sequencing error, to high differentiation (e.g., 97% similarity) (Supplementary Fig. 4). These analyses highlight that, in general, vOTUs are made up of populations of highly similar viruses, rather than being clonal. This has also been observed elsewhere, including the Arctic14.

Virus-host and environmental relationships
To explore the temporal structuring of viruses and their association with prokaryotic hosts, we employed community- and taxon-level analyses in the context of eight physicochemical parameters and prokaryotic community composition data, comprising 3748 prokaryotic amplicon sequence variants (ASV).

At the whole community level, Mantel tests demonstrated the strongest correlation to prokaryotic community composition (Mantel ρ = 0.632, p = 0.001, n = 46) followed by oxygen (Mantel ρ = 0.371, p = 0.001, n = 46) and mixed layer depth (MLD) (Mantel ρ = 0.259, p = 0.001, n = 46) (Supplementary Data 3, Supplementary Fig. 5). To further elucidate physicochemical drivers, we used canonical correspondence analysis (CCA) to reveal a seasonal clustering of samples, with temperature (summer), MLD (winter), photosynthetically active radiation (PAR; late spring-early summer), and polar water fraction (late summer) accounting for 15.4% of the total variance (Fig. 3b). These results suggest that viral communities are primarily correlated with prokaryotic community composition which is in turn driven by environmental conditions plus biological interactions—leading to a complex network of interdependencies and physicochemical linkages, similar to dynamics in temperate environments7,50.

Given the coupling between viral and prokaryotic communities, we next explored potential virus-host associations over time at the individual vOTU and 16S rRNA gene ASV level. We constructed a Convergent Cross Mapping (CCM) network based on co-occurrences, similar to the approach described in Oldenburg et al.49, which includes information about the direction of associations (i.e., predicting causal relations) between vOTUs and ASVs. The CCM network comprised 5136 vOTUs and 850 prokaryotic ASVs (Supplementary Fig. 6a), with directional associations (correlations >0.7) where vOTUs dynamics are ‘following’, or are ‘caused’ by, prokaryotic ASV dynamics (akin to Lotka-Volterra dynamics51) totaling 15,930. Louvain clustering resulted in six distinct modules, comprising vOTU-ASV associations occurring during the same temporal period. We named the six modules based on their seasonal abundance patterns (Fig. 4a); for example, M1 peaks in late spring, followed by M2. The temporal distinction between modules was confirmed by their correlations with specific physicochemical conditions (Fig. 4b). M1 showed the strongest positive correlation with PAR, while M2 had the strongest positive correlation with temperature (Fig. 4b). Interestingly, these two modules also have the most diverse virus communities, with M1 containing 1448 and M2 containing 1334 vOTUs, respectively (Fig. 4c). M3 and M4 peaked in late summer, with M3 having the strongest positive correlation with polar water fraction, and M4 having the strongest negative correlation with PAR. M5 and M6 both negatively correlated with temperature, cVCR, and positively with MLD. However, among these two, only M5 negatively correlated with PAR, being particularly abundant during winter (December–April), while M6 was persistent throughout the year. Generally, the number of vOTUs in a given module outnumbered the number of prokaryotic ASVs, except M5 where ASVs outnumbered vOTUs. Notably, M5 also contained the fewest vOTUs, 132 (Fig. 4c). The overall trend of more vOTUs than ASVs in modules corresponds to the overall larger number of vOTUs vs. prokaryotic ASVs examined.

The modules were consistently dominated by Caudoviricetes (Supplementary Fig. 6b), and the proportion of lytic to lysogenic viruses was rather invariable (Supplementary Fig. 6c). Megaviricetes were primarily present in the spring and early summer modules M1 and M2 (Supplementary Fig. 6b) where they constituted 6.3% and 1.4% of the vOTUs, respectively. Together, the major differences between modules, especially relating to bacteriophages, are hence at the vOTU level.

In terms of prokaryotic membership, the modules varied considerably. M1 was dominated by taxa typically associated with copiotrophic conditions or phytoplankton blooms, including vOTUs with Flavobacteriaceae, Rhodobacteraceae, and Porticoccaceae (Fig. 4d). In contrast, M5, which contained the largest number of ASVs, comprised diverse prokaryotic taxa, including Pelagibacteraceae, Nitrosopumilaceae, GCA−002718135 (aka HIMB59), Nitrospinaceae, SAR86, Pirellulaceae, and Planctomycetaceae (Fig. 4d). The smaller modules M2, M3, M4, and M6 comprised a variety of taxa, with the most prevalent taxon within each being UBA1611 (Marinimicrobia), and Pelagibacteraceae, respectively (Fig. 4d).

In order to evaluate host association patterns, we utilized both the CCM network correlations and host predictions via iPHoP52. In total, 42.5% of vOTUs received a host prediction by iPHoP. Among bacterial families with a high number of predicted vOTUs, Flavobacteriaceae and Pelagibacteraceae were the most frequently predicted. In M1, Flavobacteriaceae and Akkermansiaceae (Verrucomicrobia) dominated, while Pelagibacteraceae was less common. Conversely, in M5, Pelagibacteraceae was the most predominant, with fewer Flavobacteriaceae and Akkermansiaceae. Meanwhile, vOTUs with Cyanobiaceae as their predicted host were more common in M2–M4 (Supplementary Fig. 6d). Among the vOTUs correlated with an ASV in the network, the majority of host-virus pairwise correlations were to diverse prokaryotes present in M5, with Pelagibacteraceae (20.3%) and Nitrosopumilaceae (6.0%) being the most common (Supplementary Fig. 6e). M1 had the second most bacteria-to-virus connections with 37.5% of these being to Flavobacteriaceae. Among the families with most pairwise correlations (i.e., over 250), the overall pattern of host prediction of the viruses and module membership was similar, with predicted hosts of Flavobacteriaceae and Pelagibacteraceae being the dominant host predictions for M1 and M5, respectively (Supplementary Fig. 6f). However for individual vOTUs and ASVs, only 3.7% of the predictions overlapped at the family level between CCM network correlations and iPHoP (i.e., 122 out of 3,280 correlations with associated host predictions). Among these, again, Flavobacteriaceae and Pelagibacteraceae made up 46.7% and 6.6%. For Flavobacteriaceae, these vOTUs were diverse, coming from 14 different novel orders (via VConTACT353, see Methods), while for Pelagibacteraceae they all came from a single novel order (Supplementary Data 4). However, overall, the low number of matching CCM correlations and host predictions at the vOTU-ASV level demonstrate the challenge to discern host-virus relationships in diverse and complex ecosystems where a variety of technical and biological challenges may complicate such analyses (see Discussion). Nonetheless, the results provide valuable indications for particular lineages at an overall module level.

Diversity and dynamics of cyanophages and Nitrosopumilaceae-associated viruses
From the predicted interactions between vOTUs and ASVs, we further explored the diversity and dynamics of putative cyanophages—as cyanobacteria are of particular interest with respect to changing conditions in Arctic ecosystems5,54,55. We identified putative cyanophage vOTUs based on phylogenetic analysis of psbA genes, classification based on VPF-Class56 and host prediction via iPHoP52. The phylogenetic reconstruction of psbA revealed that putative cyanophages are distinct from cultivated relatives from more temperate locations (Fig. 5a). The putative cyanophage vOTUs represented some of the largest viral contigs recovered, with half of the assembled cyanophages being medium-to high-quality, and a quarter >80% complete (Fig. 5b, and Supplementary Data 2). The cyanophage vOTU abundance peaked during August–September 2017 (Fig. 5b), complementing previous observations of Synechococcus in the WSC42,54. Although the seasonal dynamics were consistent across the cyanophage vOTUs, their maximal cVCR varied from 0.003 to 0.06 (Fig. 5b). In addition, all of the cyanophage vOTUs exhibited lower abundances in 2018–2020 compared to 2017, indicating interannual variation of these viruses and their presumed hosts (Fig. 5b). Overall, the abundances of the cyanophage vOTUs with Synechococcus ASVs were correlated when considering no time-lag (Spearman’s ρ = 0.48, p = 0.0008, df = 44), but a stronger correlation (Spearman’s ρ = 0.51, p = 0.0008, df = 44) occurred with a vOTU time-lag of one time-point (average interval of 32 days). Thus, putative cyanophages increased in abundance ~1 month after Synechococcus peaks. For individual ASVs and vOTUs, the strongest correlations were without time delay (n = 53), followed by one (n = 44) and two (n = 32) time-points, indicating some temporal variability between individual cyanophage vOTUs and their putative hosts compared to the groups at-large (Supplementary Data 5). Our results reveal that both cyanobacteria and their viruses are present in the Arctic, highlighting a need to further understand their interactions for the future Arctic Ocean, as they are expected to increase due to Atlantification.

Nevertheless, given the scarcity of data on viruses during the polar night, we also examined the microbial dynamics of viruses related to prokaryotes dominant in winter. To do so, we focused on viruses associated with Nitrosopumilaceae due to their importance for wintertime nitrogen and carbon cycling57,58, their prevalence in the dataset with the top five Nitrosopumilaceae ASVs peaking from October–May (average 4.2%) vs. lower abundances from June–September (average 0.5%), and being among the taxa with the most correlations to viruses (Supplementary Fig. 7). We identified associations between 58 vOTUs and five Nitrosopumilaceae ASVs (Fig. 6), which were primarily associated with the winter module M5. In total, 53 of the 58 vOTUs were classified as Caudoviricetes (the other five were unclassified), and came from 14 different novel orders, based on vConTACT3 (Supplementary Data 6). Notably, among vOTUs with correlations to hosts in the correlation network, only one had a predicted host (based on iPHoP) of Nitrosopumilaceae. This vOTU was correlated with an unclassified Gammaproteobacterium (no family-level prediction) within the correlation network, highlighting again the difficulty comparing directly these two complex and independent approaches. The persistence of these vOTUs despite low virus cVCR demonstrates the ability to track pronounced seasonality even amongst low abundance viruses. However, more focused investigations on the temporal linkages of these vOTUs and their potential hosts are needed to better understand the environmental impacts of such associations.

Bipolarity of Fram Strait viruses
Considering the long-standing discussion on latitudinal microbial diversity gradients59–62 and the endemicity of Arctic and Antarctic microbiomes63,64, we assessed the distribution of Fram Strait viruses across the global ocean through metagenomic datasets from various large-scale sampling campaigns, such as Malaspina, Tara Oceans and Bio-GO-SHIP14,65–69 (Supplementary Data 7).

Overall, the abundance of Fram Strait viruses peaked in the epipelagic (<200 m) around 60–70°N, with decreasing abundances towards 50°N, and typically no detection in subtropical and tropical epipelagic waters (Fig. 7a). A similar pattern occurred in southern hemisphere epipelagic waters, with peaks in abundance around 45–55°S (Fig. 7a). At individual vOTU level, 42% of vOTUs displayed bimodal peaks in abundance with separate maxima in each hemisphere (i.e., the two highest peaks in abundance occurred in separate hemispheres). The average peak latitude of vOTUs was 61°N and 51°S, respectively (Fig. 7b). Fram Strait viruses were also commonly detected in the deep ocean (here, operationally defined as sampling depth of >200 m). In particular, they were detected in 87% (338 of 390) of deep global samples. Like the epipelagic waters, deep-water viruses were more prevalent in northern and southern higher latitudes (Fig. 7b), though in deep samples the viruses were more commonly detected. Overall, however, when detected, their relative abundances were lower in deep samples than in shallow samples. Notably, there was significant variation in community structure of the detected Fram Strait viruses between epipelagic and deep samples (PERMANOVA: F(1, 767) = 99.16, p = 0.001; R² = 0.1145, 95% CI [0.102–0.131]) which aligns with established differences in microbial communities in the deep vs. epipelagic ocean, and past results of virus biogeography14. Future work could consider the differences in the total virus communities between epipelagic and deep, as the current analysis utilizes only mapping to the Fram Strait viruses.

Expanding on these observations, we investigated the distributional patterns of the six Fram Strait modules on the global scale. All modules were more prevalent in the northern than in the southern hemisphere, but their relative distributions varied, with M1 and M2 being the most abundant in northern epipelagic waters, while M6 being more prominent in higher southern latitudes (Supplementary Fig. 8a), indicating a relatively restricted distribution of M1 and M2. In deeper waters, all modules were less prevalent than in epipelagic waters, with higher prevalence in the northern hemisphere. In deep waters, M6 was the most prevalent (Supplementary Fig. 8b).

Distinctive amino acid signatures of polar viruses
To expand on the global perspective, we assessed potential adaptations of viruses to polar waters by examining properties that have been attributed to cold adaptation in prokaryotes, in particular amino acid signatures that increase protein flexibility in cold environments. To do this, we examined the amino acid features of proteins from the Fram Strait viruses and from the GOV 2.0 dataset, which spans both polar and non-polar sampling locations14. To ensure fair comparisons between the datasets, we focused only on GOV samples originating from ≤35 m water depth (i.e., the average depth of Fram Strait sampling) and first analysed the environmental context, which shows strong clustering of GOV 2.0 polar samples to the Fram Strait, based on their similarity in distance to equator, oxygen, and temperature (Fig. 8a). Linking those environmental data to the examined protein features for all viruses, revealed that the aliphatic index and the nitrogen usage score were positively correlated with temperature (Fig. 8b), whereas other traits of potential cold adaptation, in particular polar charged and uncharged amino acids, showed significant correlation with one or more of the other environmental parameters, but not directly with temperature. This was also the case for Caudoviricetes, while Megaviricetes amino acid traits were more generally positively correlated with oxygen, temperature, and chlorophyll, and other vOTUs (non-Caudoviricetes, non-Megaviricetes) were generally positively correlated with salinity and oxygen (Supplementary Fig. 9). In the latter cases, these groups are more poorly sampled due to their low abundance resulting in weaker correlations overall, necessitating further study. Examining the average values for these protein parameters according to sample and temperature, across all viruses, reflected the patterns observed via other statistical analyses (Supplementary Fig. 10), reinforcing our observations.

Beyond the amino acid level, we found that 17.0% of viral protein annotations were significantly enriched in high latitudes (253 of 1490), while another 7.1% were enriched in lower latitudes (106 of 1490) (Supplementary Data 8). Notably, among a variety of annotations enriched at high latitudes, most significant were chaperone proteins and stress response such as Cold (CSD) and Heat Shock Proteins (HSP70), DnaJ (Supplementary Fig. 11), as well as genes common in cyanophages in the lower latitudes including a photosystem gene (Photo_RC) and Transaldolase/Fructose-6-phosphate aldolase (TAL_FC)70. Together, our evidence suggests that biochemical and biological traits distinguish polar viruses from their tropical and subtropical counterparts, linked to amino acid sequences.

Methods
Seawater collection and eDNA sequencing
Research and sampling complied with relevant ethical and international regulations. Moored Remote Access Samplers (RAS; McLane) autonomously collected and fixed seawater from an average depth of 29 m in the eastern and western Fram Strait (Fig. 1) at weekly to fortnightly intervals between 2016–202042,58. Sampling occurred in the framework of the FRAM/HAUSGARTEN Observatory. The resulting eDNA was used to sequence 16S rRNA gene fragments using primers 515F–926R101, processed into ASVs using DADA2 as described under (https://github.com/matthiaswietz/FRAM_eDNA). Subsequently, we only considered ASVs with ≥3 reads in ≥3 samples, corresponding to a total of 3748 prokaryotic ASVs. To complement iPHoP host assignment predictions (below), we assigned taxonomy of prokaryotic ASVs via Greengenes2102 using the classify-consensus-vsearch method103 in the q2-feature-classifier104 of QIIME2105.

DNA extracts from selected timepoints were additionally used to generate PacBio HiFi metagenomes at the Max Planck Genome Center, Cologne, Germany. Further details about the molecular analyses are described in previous reports42,58. In total, we herein analyze 94 amplicon samples and 47 metagenomes from the WSC, and 9 metagenomes from the EGC (Supplementary Data 1).

Environmental parameters
Attached to the RAS were Seabird SBE37-ODO CTD sensors that measured temperature, depth, salinity, and oxygen concentration. Sensor measurements were averaged over 4 h around each seawater sampling event. Physical sensors were manufacturer-calibrated and processed in accordance with https://epic.awi.de/id/eprint/43137. Employing multiple CTD sensors along the mooring depths enabled the determination of the minimum MLD at each sampling timepoint. Chlorophyll concentrations were measured via Wetlab Ecotriplet sensors. Surface water PAR data, with a 4 km grid resolution, was obtained from AQUA-MODIS (Level-3 mapped; SeaWiFS, NASA) and extracted in QGIS v3.14.16 (http://www.qgis.org). Polar water fraction, which is the proportion of polar water in the mixing of Atlantic and polar water masses, was calculated based on salinity and temperature33. Figure 1a map was created in QGIS v3.14.16 (http://www.qgis.org) using publicly available bathymetry data obtained from GEBCO. RAS illustration in Fig. 1a was generated via Inkscape (https://inkscape.org).

Virus assembly, prediction, classification, and host prediction
Long-reads from each sample were assembled individually using hifiasm-meta v.0.13-r308106 with default settings. Viral sequences were predicted from both assembled contigs and the long-reads themselves using a combination of tools, based on a VirSorter2-based Standard Operating Procedure (SOP, https://www.protocols.io/view/viral-sequence-identification-sop-with-virsorter2-5qpvoyqebg4o/v3). First, in VirSorter2 v.2.2.3107 we included all possible viral groups: dsDNA phages, RNA viruses, ssDNA viruses, nucleocytoplasmic large DNA viruses (NCLDV), and Lavidaviridae. CheckV v.1.0.1108 was then used to quality check the long-reads and contigs identified as viruses by VirSorter2107 as well as to trim potential host regions from identified proviruses. Contigs identified as viruses by VirSorter2107 (--min-score 0.5) with at least one viral gene (predicted by CheckV108) were further screened using DeepMicroClass v.1.0.3109 (not included in the original SOP, but intended to further reduce false positives). Contigs that were classified as either eukaryotic or prokaryotic by DeepMicroClass109 were discarded. Contigs that were identified as RNA viruses by VirSorter2107 were not analyzed. We clustered the remaining viral contigs into vOTUs, as defined previously110, using CD-HIT v.4.8.1111,112 with the following parameters: cd-hit-est -M 100000 -c 0.95 -d 100 -g 1 -aS 0.85.

For classification, we used the genomad annotate function (default settings) of geNomad v.1.8.1113 which uses taxon-specific marker proteins to assign taxonomy. To further resolve taxonomy to the subfamily level (Supplementary Data 2), we used vConTACT3 v.3.0.0b65 (https://bitbucket.org/MAVERICLab/vcontact3/src/master/). To predict hosts, we employed iPHoP v.1.3.352,114–117 with default settings, along with a custom database by adding previously assembled bacterial and archaeal MAGs from the Fram Strait (PRJEB67368) and removed MAGs from one oceanic study that was not manually curated, due to their high degree of potential virus contamination (PRJNA385857), which could result in false positives. When multiple host predictions were available, we selected the prediction with the highest confidence score.

Mapping and normalization of viral abundance
Metagenomic long-reads from all samples were mapped to vOTUs (>10 Kb, from assembly- and read-based approaches) using minimap2 v.2.28-r1209 (parameter: -x asm5)118, and only the primary alignments were considered in subsequent analyses.

To quantify the abundances of vOTUs across metagenome samples, we employed a two-step normalization process to account for differences in read lengths and sequencing depth. First, we determined the mean coverage across the viral sequence length from mapped read counts (with a threshold that >25% of the viral sequence must be covered). Second, we divided the mean coverage by the estimated number of cellular genomes in each metagenome, as predicted based on the mean coverage across 16 universal single-copy ribosomal proteins119,120. We call this metric coverage-based virus to cell ratio (cVCR). This approach is similar to that used to derive virus to microbial ratios on virus to prokaryotic marker genes elsewhere75, but considers the full viral coverage. In addition to this, we also calculated the coverage of each virus per gigabase pair (cVGB) of metagenome for a given sample and demonstrated that the two metrics yield comparable values (Supplementary Fig. 1). To gain an insight into the composition of the metagenome samples in terms of prokaryotic and eukaryotic reads, we employed Tiara v.1.0.3, a machine-learning tool, to assign domain-level classifications to the raw reads121. This revealed that the large majority of reads across the metagenomes, about 90%, were of prokaryotic origin.

Estimation of viral diversity
Given that the vOTUs were derived from metagenomes with different sequencing depths, we employed an iterative subsampling approach to quantify and compare the alpha diversity of viruses across samples. For this, we applied 100 iterations of subsampling the vOTU count table at a range of different depths, spanning from 50 up to 30,000 counts at 50 count intervals. For this step, all vOTUs that exceeded the 25% breadth of coverage cutoff when using all the data for a given sample were utilized. During each iteration, richness, evenness and Shannon diversity were calculated, with mean values being determined for each 50-count interval. The functions rrarefy, specnumber and diversity from the vegan package122 were used to perform subsampling and alpha diversity calculations. The mean values were used to generate a rarefaction curve of alpha diversity using the ggplot2 package123. To assess shifts in diversity over time, we compared the mean richness values obtained from subsampling at the 1000 count interval. In addition, we computed and compared the slopes of the sample rarefaction curves up until the subsampled depths, which represents the rate of vOTU discovery with increasing viral read counts. Comparing the slopes across samples enabled an assessment of whether the sample richness patterns observed at the subsampled depth would change if the viral read count was increased. To further explore this, we also employed the iNEXT v.3.0.1 package124,125 to estimate the vOTU richness based on extrapolations to 30,000 viral read counts.

Co-occurrence network
To identify temporal co-occurrence patterns, cVCR values per vOTU and relative abundance of prokaryotic ASVs were converted into temporal profiles by Fourier transformation. Temporal profiles were constructed based on 16 Fourier coefficients, which capture the majority of observed vOTU and ASV peaks within the four-year period. Pairwise correlations between individual temporal profiles were then computed between all vOTUs and ASVs. Higher Pearson correlation values indicated similar temporal profiles. For network construction, we first calculated Pearson correlations for all pairs resulting in an undirected graph, from which we only considered correlations >0.7 after multiple testing corrections using the Benjamini-Hochberg procedure. This threshold was determined based on previous findings on similar data49,58 as well as preliminary analyses conducted on our dataset, which indicated that interactions exhibiting a correlation below this value often lacked biological relevance. To delineate strongly connected components representing co-occurring taxa, the Louvain community detection algorithm was applied to the entire graph126. Next, the putative associations were further evaluated based on Convergent Cross Mapping (CCM) in order to discern causal relationships between taxa in time-series data, as outlined by ref. 127. CCM enables the prediction of a species’ time-series based on the knowledge of another species’ time-series. Given the nature of lytic virus-host systems, where interactions can be short-lived, we acknowledge that this method must be carefully interpreted within the context of our sampling strategy (i.e., biweekly to monthly on average). Initially, we constructed a CCM network encompassing all pairwise combinations. Subsequently, we extracted the in- and outgoing edges between nodes that were also connected in the co-occurrence network, utilizing resources from (https://gitlab.com/qtb-hhu/marine/publications/framphages2024). To quantify the strength of relationships, we employed Normalized Mutual Information (NMI) to account for non-linear relations49. A permutation approach was employed to compute significance values for edge weights, with the objective of determining whether the NMI values exceeded those expected for random edges. Further details on the construction and validation of the CCM network can be found in ref. 49. The CCM network was visualized in Cytoscape v.3.10.1128. Additionally, to identify potential time-lagged correlations, we employed extended local similarity analysis129 with a maximum delay of two sampling points on the original non-Fourier transformed data.

Statistical analyses
Bray-Curtis similarity, Mantel tests, CCA, and Spearman correlation were performed in R v.4.3.1130 using the vegan122 and Hmisc131 packages. We included only vOTUs with a breadth of coverage >0.25.

Cyanophage phylogenetics and comparative genomics
Putative cyanophages were predicted via VPF-class v.0.1.256 and iPHoP52 host predictions. Subsequently, we focused on the most confident predictions by selecting cyanophages harboring the psbA core photosystem, identified by hmmscan132 of viral predicted proteins (prodigal133) PFAM for Photo_RC (PF00124). Viral and prokaryotic reference psbA were extracted from the vConTACT2 viral reference database (ViralRefSeq-prokaryotes-v21153) and GTDB prokaryotic reference database (version 214134) by hmmscan132 with PFAM Photo_RC. Eukaryotic reference psbA genes were extracted via Uniprot (gene_exact:psbA) AND (taxonomy_id:2759, reviewed Swiss-Prot, n=157)135. A preliminary tree was constructed from these sequences to differentiate between D1 (psbA) and D2 (psbD) photosystem II (PSII) reaction center proteins, by alignment using MAFFT136 followed by trimming with trimAl (-gt 0.2)137, and FastTree138 for phylogenetic reconstruction. psbD were identified manually and excluded. The remaining sequences were considered, and phylogenetic reconstruction performed a second time with the same alignment and trimming steps (but with -gt 0.8); at this stage, tree building was performed with IQ-TREE (settings: -B 1000, -alrt 1000, -m MFP)139.

Mapping to global metagenomes
Metagenomes from ten oceanic datasets14,65–69,140–145, including Tara Oceans, Malaspina, and Bio-GO-SHIP (Supplementary Data 7), were downloaded from NCBI, and quality-filtered with Trimmomatic (parameters: LEADING:3 TRAILING:3 SLIDINGWINDOW:50:30 MINLEN:50)146. All data utilized is previously published. Only DNA samples from discrete depths (e.g., not composites of the mixed layers or similar) and the prokaryotic size-fraction (0.2–3 µm) were used to allow comparison with our dataset. Multiple other size fractions are available, but the 0.2–3 µm was utilized because of our focus on prokaryotic dsDNA viruses, its comparability to the present Fram Strait virus dataset, and because it is among the most common range sampled (i.e., many samples available). Reads were mapped to the clustered contigs using bbmap.sh from the BBTools package (version 39.01, parameters: minid=95 idfilt=95) (https://sourceforge.net/projects/bbmap/). Trimmed mean coverage was calculated via CoverM147 on contigs with minimum covered fraction of >25% for each sample. Trimmed mean coverage values were then divided by total Gbp per quality filtered paired reads. Maps were generated via publicly available NOAA bathymetry data via getNOAA.bathy and the maps package in R. For determining abundance trends in the total dataset and at the module level, we computed a Generalized Additive Module (GAM) between cVCR and latitude using the gam function from the mgcv package148, with ten degrees of freedom (k=10). For plotting the predicted GAM, negative values were converted to 0. For calculating the maximum latitude that viruses occurred in each hemisphere, we used the GAM model, computing where peaks occur and then determining which two peaks had the highest abundance. For evaluating statistical differences in community similarity between epipelagic (<200 m) and deep (>200 m) samples, we first removed samples that had less than 0.0001 cVGB (n=769 samples remaining, including 435 and 334, shallow and deep samples, respectively) and vOTUs that had less than 0.0001 summed cVGB across all samples (n=4,693 vOTUs remaining). Then, we ran adonis2122,149 with 999 permutations based on the epipelagic and deep grouping, with 95% confidence intervals for the R² statistic estimated using 1000 bootstrap replicates.

Microdiversity of viral contigs
To examine microdiversity among vOTUs, we used two separate approaches. First, to examine the extent that exact, long virus sequences re-occur across years, we used CD-HIT v.4.8.1111 to cluster viral contigs >10 Kb at a 100% sequence identity threshold (cd-hit-est -M 100000 -c 0.95 -d 100 -g 1 -aS 0.85). We then identified clusters that appeared across more than one year throughout the four-year sampling period. Second, we used inStrain v.1.9.0150, which compares reads mapped to vOTUs to characterize virus populations based on nucleotide diversity, total counts of single nucleotide substitution (SNS) and single nucleotide variation (SNV), and other related metrics. We applied a 98% minimum ANI, the most stringent threshold recommended by the developers, and used a minimum coverage of five to call a variant. To confirm a SNV, we required a minimum SNV frequency of 0.05 and a false discovery rate of 1e-06. Additionally, we filtered for contigs where more than 25% of their bases were covered by at least one read, i.e., breadth >0.25. Only vOTUs that exceeded the coverage minimum more than three times throughout the sampling period were included in the analysis.

Calculation of amino acid traits
Proteins were predicted from GOV 2.0 and Fram Strait viruses via Prodigal v.2.6.3 (parameter: -p meta)133. Then, amino acid composition and various biochemical properties of the predicted protein sequences were assessed using custom R scripts available at https://github.com/alyzzabc/fram_strait_viruses_2016-2020. Briefly, protein sequences containing the ambiguous amino acid “X” were removed, followed by counting the occurrences of each proteinogenic amino acid. The frequencies of single amino acids (G, S, P, C) and of amino acid groups (acidic, polar uncharged, polar charged, aromatic) were calculated relative to the protein length. For the amino acids arginine (R) and lysine (K), the ratio R/K was calculated instead of occurrence. The molecular weight was calculated according to 151 and for calculation of the aliphatic index we followed 152. The pI function from the Peptides package153 was used to calculate the pI of the predicted proteins using the EMBOSS scale. The nitrogen usage score (NUS) was calculated using the following formula:

NUS = (4F(R) + 3F(H) + 2F(KNQW) + F(DESTGPCAVILMFY)) / n(protein) (1)

where F(X) corresponds to the absolute count of each amino acid and n(protein) to length of the protein.

To analyze the variation of the viral protein features in relation to the environmental conditions, we omitted samples collected from a water depth >35 m followed by calculation of the Bray-Curtis distance similarity matrix and generation of NMDS coordinates using the metaMDS function from the vegan package122.

To calculate pairwise correlations, the vegan package122 was used: the distances for each protein feature (Bray-Curtis) and environmental condition (Euclidean) were calculated using the vegdist and dist function, respectively. Pairwise Mantel tests of each environmental parameter vs. each amino acid trait were done using the Spearman method with 9999 permutations via the mantel function. The NMDS plot and the correlation heatmap were prepared using the ggplot2 package123.

Calculation of protein family enrichments by latitude
For determination of proteins enriched in polar regions, we used the prodigal-predicted proteins for GOV 2.0 and the Fram Strait virus dataset. These predicted proteins were annotated with pfam154 via hmmer132 using the protein-specific cutoffs for significant (cut_ga). Then, individual protein families were summed within a given sample. Proteins that had prevalences of less than 10% were removed (n=1,490 remaining). Then, ANCOM-BC155 was used to determine which annotations were enriched at latitudes greater than 60° N or S, or less than 60° N or S.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Calayag et al. present a well-conducted and comprehensive study of dsDNA viral ecology in the Fram Strait, Arctic, addressing a crucial gap in our understanding of virus dynamics in polar regions. The study demonstrates the seasonal structuring of viral communities over a four-year period, with distinct annual peaks in the virus-to-prokaryote ratio during the summer months. Additionally, the discovery that a significant proportion of Arctic viruses exhibit bipolar distributions further contributes to our knowledge of their biogeography. The use of convergent cross mapping network to identify virus-host interactions within ecological modules enhances our understanding of DNA viruses' ecological roles and interactions. Moreover, the analysis of amino acid signatures in viral proteins provides valuable insights into how viruses adapt to such extreme environments.\nOverall, this study is valuable to environmental virologists, offering both important analysis results and data on the ecology of viruses in polar environments. The use of long-read metagenomics strengthens the findings by capturing larger genomic fragments. Additionally, the use of cVPR and CCM provides novel methods for understanding virus-host interactions.\nHowever, before this manuscript is considered for publication, I have some suggestions and comments listed below.\nGeneral Comments:\n# 1. In Fig. 5b, many vOTUs do not encode structural or packaging-related genes, raising my concerns about the completeness of these cyanophage genomes. I wonder if this is an issue of labeling in this plot or if it reflects the nature of these genomes. It would be helpful if the authors could provide CheckV results along with basic statistics for the viral genomes to clarify genome completeness and quality.\nAdditionally, the bootstrap support in Fig. 5a for the basal topology appears low, with most nodes having bootstrap values < 80. To improve the reliability of this cyanophage phylogeny, I suggest the authors consider using alternative genes or incorporating cellular psbA genes into the tree.\n# 2. While Convergent Cross Mapping (CCM) appears to be a robust method for detecting causality in complex systems, I have some concerns about the application to lytic virus-host systems (high proportion in this study), where the short duration of interactions (with latent periods ranging from few hours to days) may complicate the inference of directional causality under a resolution of monthly sampling. This could potentially explain the limited overlap (4.3%; LN 221) between CCM associations and iPHoP predictions. I recommend that the authors provide further justification for their choice to use CCM in this context, potentially addressing these challenges in the methods section or adding a caveat in the discussion.\nAdditionally, the use of certain parameters, such as a corrected Pearson correlation threshold (>0.7), may overlook biologically meaningful interactions, as viruses and their hosts do not always perfectly co-occur and gain a high correlation coefficient. A more detailed explanation of these parameter choices would be helpful to people who want to follow the study. Alternatively, the authors could present a naive network (without edge filtration) to compare against their CCM-based inferences, particularly regarding module clustering and mutual pairs identified by iPHoP.\n# 3. The amino acid composition analysis is presented as pooled data across all viral genomes, making it difficult for me to determine whether the observed amino acid signatures are driven by environmental adaptation or reflect community composition differences, especially given the strong viral endemism in polar regions. I suggest that authors consider analyzing amino acid traits in relation to phylogenetic distance, or constraining the analysis by same viral taxonomic groups, to better distinguish evolutionary lineage from environmental adaptation.\n# 4. Authors have highlighted amino acid traits, but it would also be valuable to consider the role of other mechanisms, such as gene gain, in viral adaptation. Have authors examined whether there are gene orthologs specifically present or absent in polar regions? Or an additional discussion on how mechanisms may influence viral adaptations in polar environments could help readers understand virus adaptation to environments.\nMinor Comments:\nLN135-139: The color palette in Fig. 2 makes it challenging to check the left panel. It would be helpful to use two distinct color sets between the two groups (Jan-May vs June-Sep).\nLN142: The lower diversity and slower accumulation observed in summer could be a result of higher overall biomass during that season. Few dominant bacteria and viruses may account for the majority of reads at the current sequencing depth (10Gb), so the detection of rare viruses has become relatively difficult. Please consider adding this as a potential limitation in the text.\nFig. 3b: Using different shapes for the dots to represent different years could make it easier to understand temporal trends across four years.\nLN168: I’m uncertain whether the phrase \"shaped by their hosts\" is appropriate, as it's difficult to attribute the observed viral patterns to direct host influence, given the potential roles of other prokaryotes. Additionally, I am a bit confused about which pair in Supplementary Data 1 represents the structure of viruses in the main text. Authors could add the values for each variable in the main text.\nLN179: Please add the number of ‘following’ and ‘caused by’.\nLN192-193 & Fig. 4d: The ratio of vOTUs/ASVs in M5 is significantly lower than in other modules, and the value is somehow counterintuitive. Could this be due to that most winter viral nodes being discarded during network preprocessing? A bit further discussion on potential reasons behind this observation would be helpful.\nLN265: Are the 58 vOTUs mentioned evolutionarily related (similar taxa)?\nFig. 4C: Please ensure the p-values have been corrected for multiple testing.\nSupplementary Data: It is not easy to understand the headers in supplementary data 1 and 2. Please consider using one spreadsheet to include descriptions of the column headers for all tables.", "reviewer_2": "Review Nature Communications. NCOMMS-24-55589\nArctic Ocean virus communities: seasonality, bipolarity, and prokaryotic Interactions\nThis study investigates the seasonality, abundance, diversity, associations with prokaryotic hosts, and environmental factors of Artic viruses in the Fram Strait over a four-years period. The authors used 47 long-read metagenomes obtained from the cellular size-fraction, which were automatically collected. It was found that viruses exhibited annual peaks, and that the composition of the viral community was strongly correlated with prokaryotic community composition. Additionally, it was analyzed presence of the polar virus in metagenome from different latitudes, concluding that virus abundance is higher in high-latitude regions of both poles, highlighting adaptations for cold environments.\nFirst, I would like to positively emphasize the rigorous mathematical and statistical methods employed by the authors in this work. While seasonal cyclic variability in viral communities has been observed in prior studies, the extended sampling period here provides a more robust dataset than previous research. The constructed metagenomes alone represent a valuable resource database. However, some sections lack sufficient biological interpretation and related discussion. Further clarification from the authors on the following specific points is recommended to enhance the overall quality of the manuscript.\nIn general, more quantitative data should be included in the main text and as Supplementary Data. Table of Supplementary Data 4, in particular, could be expanded to include specific information such as the number of viral sequences in each collection, the depth (or range) of the water column from which each sample was obtained, the maximum length of the viral raw reads obtained (L239), the estimate of the number of prokaryotic genomes in each metagenome, and, if possible, parameters indicating viral (L142) and prokaryotic diversity. Additionally, a global overview of prokaryotic taxonomy retrieved by ASVs and their periodicity is needed in the Supplementary Data to improve the interpretation of the seasonality observed in the viruses (L142, L80). Details regarding the sizes and number or replicas of the 16S rRNA libraries should also be provided as Supplementary Data). Furthermore, data on the ASV-vOTUs involved in the six modules, as well as those identified by iPHoP (L181 & L214), must also be included.\nTo better grasp the magnitudes involved (as the authors provide statistical data but almost no “raw numbers”), it would be helpful to include precise information in the main text. For example, it would be useful to know how many viral sequences were identified among the nearly 200 000 sequences obtained per metagenome (L96); how many viral sequences represent the 5662 vOTUs identified by the authors; which vOTU was the most abundant and how it varied along time and how many vOTUS are involved in each module (L173-195). Additionally, the 3748 ASVs mentioned in L391 should be cited earlier, at the beginning of the Results sections.\nSpecific questions:\n1.-In my opinion, the title or the article is inadequate. First, not all viruses form the community were considered, as RNA viruses were excluded from the analysis (L422), and the authors focused primarily on dsDNA viruses. Importantly, the authors examined viruses found exclusively in the cellular fraction; i.e., those actively infecting cells at the sampling points or integrated into cellular genomes as lysogens. Additionally, there is limited information on prokaryotic “interactions” per se, as host identification relies exclusively on computational predictions. I believe a title emphasizing the seasonality and the dynamics of the viral community would better reflect the study´s scope and attract a broader audience. Moreover, only viruses present in the cellular fraction are considered in the analysis, rather than the complete community- an important detail that should be more clearly stated throughout the manuscript.\n2.-The introduction provides quotes related to the ´Artic Atlantification`, but this information is not directly relevant to the goals or findings of this manuscript. In addition to the hypothesis proposed by the authors, L79, which is not entirely novel (reference 22, 25), I would encourage the authors to present hypotheses in the introduction regarding how polar viral communities would fit in the ecological models already described for virus-host populations., e.g., ‘Piggyback-the-Winner’ theory, Red-Queen dynamics, Constant Diversity Model, etc. (F.R. Valera et al. 2009, B. Knowles et al. 2016, F. H. Coutinho et al. 2017, H. Alrasheed et al. 2019, J.C. Ignacio-Espinoza et al. 2020). This would help frame the discussion to better align with fit the Lotka-Volterra dynamics observed.\n3. There are two points that need to be clarified regarding the methodology used:\n(i) Importantly, the authors claim that the lytic lifestyle was predominant over lysogeny. Although I find value in the analysis conducted by the authors, I have a primary concern regarding the number of lysogens detected. Beyond bioinformatic predictions, the authors lack direct evidence to fully support this conclusion. Many lysogens do not have recognizable integrases, which is a general challenge that needs to be addressed (out in the context of this work). In my view, viruses replicating inside the cells (i.e., lytic viruses) would outnumber lysogens in term of copies within the cellular metagenome, making them easier to assemble with the long-reads used, compared to lysogens. Therefore, a strong bias could obscure the results.\n(ii) Cut-off used to define the abundances of vOTUs. The authors use only 25% of the viral genome length, whereas typical values range from 50 to 70% at a minimum.\n4. The authors focus on Cyanobacteria and the dynamics of putative cyanophages (L232-257) due to their relevance to changing conditions in Arctic ecosystems (but it is unknown in which module are founded), along with another nice example of the Nitrososphaerales. In my opinion, there is a lack of detailed information regarding the largest group, which is represented in several modules and with multiple identified virus-host pairs, the Pelagibacteraceae. Providing information on how the Pelagibacteraceae ASVs are distributed across the different modules, along with the assignment vOTUs, would better complement this section of the manuscript.\nOther details:\nJust as a suggestion, did the authors consider whether Supplementary Figure 3c could appear as a main figure (perhaps in place of Figure 4a)?\nL63. Please, clarify. Sandaa et al 2018, (reference 22) found that (cited textually) “Viral diversity and virus-to-prokaryote ratios (VPRs) dropped sharply at the commencement of the spring bloom but increased across the season, ultimately achieving the highest levels during the winter season. These findings suggest that viral lysis may be an important process during the polar winter, when productivity is low”. Therefore, this contradicts what is stated in L59, 60.\nL93-94. Did the authors retrieve eukaryotic DNA in their sequencing? Since no pre-filtering was applied, there is a high possibility of this occurring, yet it is not mentioned in the manuscript. If present, could this affect the cVPR calculations?\nL98-99. The program CheckV is mainly used to confirm the nature of the viral sequences, not to make predictions per se. Therefore, remove from this form the line, and then the two approaches mentioned in L99 that refers the authors, will be correct.\nL107. “95% sequence clusters”. Please modify to indicate that this refers to 95% ID.\nL109. Please note that it is recommended to follow the updated nomenclature for phages. The bacteriophage families Myoviridae, Podoviridae, and Siphoviridae have been reclassified under a new taxonomic system (International Committee on Taxonomy of Viruses (ICTV), 2021). These are now part of the order Caudoviricetes, which encompasses tailed bacteriophages. Myoviridae members are often reclassified in families like Myoviricetes, while Podoviridae members generally fall under Autographiviridae, and Siphoviridae have been split across families like Drexlerviridae and Demerecviridae, among others.\nL111. Phycodnaviridae continues to be cited in in L122, L199 and Figure 3, so please, update this sentence accordingly.\nL166. Please define `PAR’ as it is the first time it appears in the text. Also, what parameter exactly is the “polar water fraction” in Figure 3b? Additionally, “iceDist” and “icePast” are not mentioned in the manuscript (also in Figure 4).\nL287. Is this in disagreement with the statement in L280? First, authors must clarify whether the vOTUs detected in surface waters are different from those detected in deep waters. If Fram Strait viruses are also abundant in deep ocean samples, does this imply that similar hosts are more abundant at these depths? Additionally, which depth range is being referred to? This information needs to be either depicted in the figure or stated in the Supplementary Data 3. It is generally understood that deep and surface prokaryotic communities differ… then, could mixed waters be influencing this observation (although the authors specifically mention that these “mixed” samples were not considered in the analysis)? The authors should provide a more detailed discussion on this point (L291-298).\nL310-L317. Please, provide more details on the changes observed in the amino acid composition of the viruses from Fram and GOV2.0 datasets versus those found in in northern or southern lower latitudes. The authors provide statistics, but there is no percentage or other parameter that can be quickly appreciated.\nL331. The quantity of free viruses cited here were obtained in 2014 during in partial transect compared to the one followed in this work. Although due to the periodicity observed, it is likely that those numbers were similar, the authors could support this by using metaviromes and metagenomes recovered from the same sample and check the vOTUs abundance recovered in both.\nL438. RaFAH results are not shown anywhere in the text. The iPHoP results incorporate RaFAH, which helps improve phage-host association predictions. Please, correct the sentence.\nL452. Is reference 88 correctly place here? Please cite the original source where it is described.\nL536. The metagenomes chosen to map the viruses were those obtained using a prokaryotic size fraction of 0.2-3 microns. However, authors do not cite to use any pre-filtering, so, collections with larger sizes could also be suitable. Please, revise the sentence accordingly, as these metagenomes were selected because the authors focused on prokaryotic dsDNA virus relationships.\nFigure 4a. Please, modify this figure to better differentiate viruses from ASVs.\nFigure 5a. psbA gene in italics.\nPlease, enlarge the figures for improved readability. At their current size, it is challenging to distinguish finer details.\nFigure S3c. May the authors order the microbial family names accordingly to a taxonomic relationship? It would make more sense than organizing them alphabetically.", "reviewer_3": "The authors present a well-written paper about a novel dataset that describes a rare 4-year time series of viral metagenomes form the Arctic. They introduce the broad seasonal patterns of viruses they observed and add extra analysis by comparing the viruses they found with other publicly available datasets of that type. They use a fairly straight forward analysis with appropriate and common pipelines and tools.\nHowever, there is one noticeable exception: the authors look at a mix of morphotype classification (podo-, myo- and sipho-viridae) for bacteriophages and genetic based taxonomic groups for eukaryote viruses (Phycodnaviridae, Mimiviridae). This seems like an outdated approach since the ICTV reported that they are not monophyletic, phylogenetically solid groups (See e.g. “Taxonomy of prokaryotic viruses: 2018-2019 update from the ICTV Bacterial and Archaeal Viruses Subcommittee”). A direct comparison between those different types of classification groups not ideal, and there is no clear reason why the more recent classification has not been chosen. The main focus of the paper is not on that aspect, so it is not an unfixable problem. But please update the analysis according to the most recent ICTV guidelines, or at least justify why the chosen approach is still valid (I will leave that up to the editor). The authors should at minimum acknowledge this issue in the text and explain why this comparison was chosen, rather than the more up-to-date classification established in recent years.\nIt is crucial that the authors carefully go through the entire manuscript and make absolutely sure that they provide a citation for every software tool, package and database and whatever else they use. Studies like this one would not be possible without all the public databases and free software tools developed by other scientists, it is only good practice to give everyone credit for their work.\nMinor comments:\n22-23: Perhaps specify that the paper is about viruses of microbes.\n26-28: Wording could be clearer.\n62-63: A new paper has been released after the authors would have submitted their manuscript, but there is a viral metagenomic time-series that might be relevant to consider here. https://academic.oup.com/ismej/advance-article/doi/10.1093/ismejo/wrae216/7833430\n78-80: It is known that microbes in the polar regions are highly seasonal and, obviously, virus abundance is directly linked to host abundance, a metric the authors use as well. Given that seasonality of viral communities has been shown in other regions and has been speculated for polar viruses as well, this is not a good hypothesis in my opinion. The authors present a very important and interesting dataset, which as far as I know is the first virome time series from the Arctic Ocean. They should emphasize the novelty of this dataset, rather than the general concept of viral seasonality.\nAn additional side point: the manuscript only provides data from the Fram Strait, so the hypothesis should be about Arctic seasonality, not polar seasonality, in any case. It’s reasonable to speculate similar patterns for the Southern Ocean, but the seasonal dataset only really covers the Arctic.\n90-91: What was the depths these samples were collected from? Please clarify.\n117: s missing for ratios?\n166: Please introduce the meaning of abbreviation PAR.\n162-171: Please provide a value or a figure for this. The stats are buried in a difficult-to-read supplementary table, but this is a fairly important point and the following paragraph is building on that point. it should be made clear to the reader how significant that correlation is compared to oxygen and MLD. The authors provide a main figure (3) and a percentage to show the results for the less important physiochemical drivers, but not for the correlation with their hosts. The text and figures should reflect the more important aspects of results first.\n241-242: This is not surprising since the authors group the viruses based on the outdated morphological classification, which is polyphyletic.\n287: Define what depth is considered deep ocean. Based on context, > 200m depth.\n467: Please cite the authors of the vegan and ggplot2 packages and acknowledge the use of R.\n501: Please reference the original Cytoscape paper.\n518: Please cite the original vConTACT2 paper.\n519: Please cite the hmmer paper.\n520: Please add citation for UniRef.\n532: Please add citations for Tara Oceans, Malaspina Expedition and Bio-GO-SHIP. They are hidden in a supplementary table as a link only, but should be properly cited in the text itself too. I will stop checking if the authors have referenced every tool, database and package, since they tend to only add a citation upon first mention, if at all, so it’s tedious and easy to miss."}, "author_response": "REVIEWER COMMENTS\nResponses in blue\nReviewer #1 (Remarks to the Author):\nCalayag et al. present a well-conducted and comprehensive study of dsDNA viral\necology in the Fram Strait, Arctic, addressing a crucial gap in our understanding of virus\ndynamics in polar regions. The study demonstrates the seasonal structuring of viral\ncommunities over a four-year period, with distinct annual peaks in the virus-to-\nprokaryote ratio during the summer months. Additionally, the discovery that a significant\nproportion of Arctic viruses exhibit bipolar distributions further contributes to our\nknowledge of their biogeography. The use of convergent cross mapping network to\nidentify virus-host interactions within ecological modules enhances our understanding of\nDNA viruses' ecological roles and interactions. Moreover, the analysis of amino acid\nsignatures in viral proteins provides valuable insights into how viruses adapt to such\nextreme environments.\nOverall, this study is valuable to environmental virologists, offering both important\nanalysis results and data on the ecology of viruses in polar environments. The use of\nlong-read metagenomics strengthens the findings by capturing larger genomic\nfragments. Additionally, the use of cVPR and CCM provides novel methods for\nunderstanding virus-host interactions. However, before this manuscript is considered for\npublication, I have some suggestions and comments listed below.\nWe thank the reviewer for their positive assessment regarding the value of our paper,\nand appreciate their comments as to the novelty of the methods and metrics we used.\nBelow we address their suggestions and comments which we believe have improved\nthe paper.\nGeneral Comments:\n# 1. In Fig. 5b, many vOTUs do not encode structural or packaging-related genes,\nraising my concerns about the completeness of these cyanophage genomes. I wonder if\nthis is an issue of labeling in this plot or if it reflects the nature of these genomes. It\nwould be helpful if the authors could provide CheckV results along with basic statistics\nfor the viral genomes to clarify genome completeness and quality.\nWe thank the reviewer for this constructive comment. Indeed, we only sporadically\nidentified structural and other marker genes within the cyanophages. The goal was to\nhighlight some genes to give the reader a sense of ‘what is where’, but not\nsystematically label as it would be too cluttered. Given that it is not really possible to do\nlabel exhaustively and the figure was already cluttered, we have removed this part, and\nreplaced it with the much more quickly understandable CheckV (completion) results, as\na barplot next to the heatmap to show the estimated completion of each respective\nvOTU.\nWe also now provide CheckV for all vOTUs as part of Supplementary Data 2.\nAdditionally, the bootstrap support in Fig. 5a for the basal topology appears low, with\nmost nodes having bootstrap values < 80. To improve the reliability of this cyanophage\nphylogeny, I suggest the authors consider using alternative genes or incorporating\ncellular psbA genes into the tree.\nWe thank the reviewer for pointing this out. The psbA gene was utilized as it enabled\nestablishing a phylogenetic relationship for all of the predicted cyanophage vOTUs, and\nalso, as the reviewer points out, cellular organisms. We have now included cellular\npsbA and improved the bootstrap support for the tree, as well as a second form of\nsupport via Shimodaira-Hasegawa-like approximate likelihood ratio test. We now\nindicate where the cellular psbA are located in the tree in the main figure. The bootstrap\nsupport is still sometimes lower than ideal, but upon reviewing the literature, other\nversions of the psbA phylogeny often suffers from the same issue. For us, it was just\nimportant, mainly to somehow orient/cluster the most similar vOTUs together to give a\nsense of how closely related the cyanophage vOTUs are to one another. We now\nprovide an update to the figure legend.\nFigure legend:\nLines 1951–1957: Fig. 5. Cyanophages. a, psbA gene phylogeny of cyanophage\nvOTUs and reference genomes. Support values are reported as aLRT (approximate\nlikelihood ratio test) / BS (maximum likelihood bootstrap). For visualization purposes, “*”\nindicates branches from this node are scaled 33% of their actual length for visualization\npurposes. b, Dynamics of correlated Synechococcus ASVs (top) and cyanophage\nvOTUs (bottom), along with completion estimates based on CheckV108. ** indicates no\ncompletion prediction as the vOTU was classified as a provirus.\nWe update the methods accordingly.\nLines 926–957: Viral and prokaryotic reference psbA were extracted from the\nvConTACT2 viral reference database (ViralRefSeq-prokaryotes-v21152) and GTDB\nprokaryotic reference database (version 214132) by hmmscan133 with PFAM Photo_RC.\nEukaryotic reference psbA genes were extracted via Uniprot (gene_exact:psbA) AND\n(taxonomy_id:2759, reviewed Swiss-Prot, n=157)134. A preliminary tree was constructed\nfrom these sequences to differentiate between D1 (psbA) and D2 (psbD) photosystem II\n(PSII) reaction centre proteins, by alignment using MAFFT135 followed by trimming with\ntrimAl (-gt 0.2)136, and FastTree137 for phylogenetic reconstruction. PsbD were identified\nmanually and excluded. The remaining sequences were considered PsbA and\nphylogenetic reconstruction performed a second time with the same alignment and\ntrimming steps (but with -gt 0.8); at this stage, tree building was performed with IQ-\nTREE (settings: -B 1000, -alrt 1000, -m MFP)138.\n# 2. While Convergent Cross Mapping (CCM) appears to be a robust method for\ndetecting causality in complex systems, I have some concerns about the application to\nlytic virus-host systems (high proportion in this study), where the short duration of\ninteractions (with latent periods ranging from few hours to days) may complicate the\ninference of directional causality under a resolution of monthly sampling. This could\npotentially explain the limited overlap (4.3%; LN 221) between CCM associations and\niPHoP predictions. I recommend that the authors provide further justification for their\nchoice to use CCM in this context, potentially addressing these challenges in the\nmethods section or adding a caveat in the discussion.\nWe acknowledge that the short duration of interactions in lytic virus-host systems, with\nlatent periods ranging from a few hours to days, poses a challenge in inferring\ndirectional causality with monthly sampling resolution. We have chosen CCM because it\noffers a robust framework for detecting causality in complex ecosystems where\ntraditional correlation methods may fall short. CCM allows us to infer potential causal\nrelationships based on time-series data by analyzing the dynamics of both viruses and\ntheir prokaryotic hosts. The limited overlap (4.3%) between CCM associations and\niPHoP predictions has been clarified in our revised manuscript. We hypothesize that this\nlow overlap may be due to the temporal resolution of our data and the inherent\nbiological variability in virus-host interactions. In the revised Discussion section, we\nelaborate on how different methodologies capture different aspects of viral dynamics,\nand we emphasize the need for further investigation to reconcile these findings.\nUpdated Discussion:\nLines 652–660: The use of CCM in studying virus-host interactions presents challenges\ndue to their rapid temporal dynamics. Viral latent periods range from hours to days,\ncomplicating interpretations based on monthly sampling intervals. The low overlap\n(3.7%) between CCM associations and iPHoP predictions may reflect limitations of our\nsampling resolution. Integrating multiple methodologies is critical for understanding the\ncomplexities of microbial interactions. Combining CCM findings with naive network\nanalysis and iPHoP results provides a more comprehensive view of microbial\ncommunity dynamics, with visual comparisons highlighting the biological significance of\nidentified interactions.\nAdditionally, the use of certain parameters, such as a corrected Pearson correlation\nthreshold (>0.7), may overlook biologically meaningful interactions, as viruses and their\nhosts do not always perfectly co-occur and gain a high correlation coefficient. A more\ndetailed explanation of these parameter choices would be helpful to people who want to\nfollow the study. Alternatively, the authors could present a naive network (without edge\nfiltration) to compare against their CCM-based inferences, particularly regarding module\nclustering and mutual pairs identified by iPHoP.\nRegarding the choice of parameters, particularly the corrected Pearson correlation\nthreshold of >0.7, we recognize that high correlation coefficients may not always\nrepresent biologically meaningful interactions. We have revised the Methods section to\ninclude a detailed rationale for our parameter choices, explaining how these thresholds\nwere determined based on preliminary analyses and ecological relevance.\nMethods, Lines 881–894: This threshold was determined based on previous findings on\nsimilar data58,124, as well as preliminary analyses conducted on our dataset, which\nindicated that interactions exhibiting a correlation below this value often lacked\nbiological relevance.\n# 3. The amino acid composition analysis is presented as pooled data across all viral\ngenomes, making it difficult for me to determine whether the observed amino acid\nsignatures are driven by environmental adaptation or reflect community composition\ndifferences, especially given the strong viral endemism in polar regions. I suggest that\nauthors consider analyzing amino acid traits in relation to phylogenetic distance, or\nconstraining the analysis by same viral taxonomic groups, to better distinguish\nevolutionary lineage from environmental adaptation.\nWe appreciate the reviewer’s point that the observed amino acid signatures could be\ndriven by environmental adaptation or community composition differences.\nThese factors may be difficult to tease apart, but to improve our analysis, at the\nreviewer's latter suggestion, we have now re-performed separate analyses on the\nCaudoviricetes, Megaviricietes, and ‘other’. Roughly 95% of the vOTUs are\nCaudoviricetes, but the patterns are the same across all groupings.\nWe update the text throughout to reflect the new taxonomy.\n# 4. Authors have highlighted amino acid traits, but it would also be valuable to consider\nthe role of other mechanisms, such as gene gain, in viral adaptation. Have authors\nexamined whether there are gene orthologs specifically present or absent in polar\nregions? Or an additional discussion on how mechanisms may influence viral\nadaptations in polar environments could help readers understand virus adaptation to\nenvironments.\nWe appreciate the reviewer’s point. We now include a table and a supplementary figure\nindicating the genes that are most enriched in high (greater than > 60° latitude N or S)\nand lower latitudes (< 60° degrees latitude).\nThis revealed some interesting patterns, which we highlight in the text:\nLines 588–595: Beyond the amino acid level, we found that 17.0% of viral protein\nannotations were significantly enriched in high latitudes (253 of 1,490), while another\n7.1% were enriched in lower latitudes (106 of 1,490) (Supplementary Data 8). Notably,\namong a variety of annotations enriched at high latitudes, among the most significant\nwere chaperone proteins and stress response such as Cold (CSD) and Heat Shock\nProteins (HSP70), DnaJ (Supplementary Fig. 9), as well as genes common in\ncyanophages in the lower latitudes including a photosystem gene (Photo_RC) and\nTransaldolase/Fructose-6-phosphate aldolase (TAL_FC)70.\nTogether, the annotation- and amino-acid-based enrichment analyses suggest that\nthere are biochemical and biological traits distinguishing the polar viruses from the\ntropical and sub-tropical counterparts encoded in their amino acids.\nCorrespondingly, we add to the methods:\nLines 1035-1043: Calculation of protein family enrichments by latitude\nFor determination of proteins enriched in polar regions, we used the prodigal-predicted\nproteins for GOV2.0 and the Fram Strait virus dataset. These predicted proteins were\nannotated with pfam152 via hmmer133 using the protein-specific cutoffs for significant\n(cut_ga). Then, individual protein families were summed within a given sample. Proteins\nthat had prevalences of less than 10% were removed (n=1,490 remaining). Then,\nANCOM-BC153 was used to determine which annotations were enriched at latitudes\ngreater than 60° N or S, or less than 60° N or S.\nWe also add a figure legend:\nLines 2060–2063: Supplementary Fig. 10: ANCOM-BC results showing protein families\nthat are significantly enriched (top 10) at sampling locations of greater than 60 latitude\n(p < 0.001) or less than 60 latitude. The remaining list with p-values and significance\nscores as Supplementary Data 8.\nMinor Comments:\nLN135-139: The color palette in Fig. 2 makes it challenging to check the left panel. It\nwould be helpful to use two distinct color sets between the two groups (Jan-May vs\nJune-Sep).\nWe now use two distinct color sets for the two groups of months that have high and low\ndiversity. The group of months with low diversity spans March and May, while the group\nwith higher diversity spans August and February.\nHowever, it should be noted that since the original submission we have re-run the alpha\ndiversity analysis because we realized that it is better to rarify the data and re-calculate\nthe statistics. This analysis revealed slightly different results for the alpha diversity\nanalysis, resulting in higher richness of vOTUs in late summer-autumn (Aug - Nov), and\nlower diversity in March and May. However, large variations in richness values\nobserved within months suggests variable viral dynamics.\nWe update the manuscript accordingly in the results (as long as changes to the\nmethods).\nLines 194–207: As viral diversity was not saturated at any sequencing depth (Fig. 2a),\nwe compared the diversity of vOTUs across months after subsampling to a normalized\ndepth. The diversity of viruses increased following the spring phytoplankton bloom,\nreaching a maximum between late summer-autumn (Aug–Nov) (Fig. 2b), with lowest\nrichness values in March and May. Notably, prokaryotic community richness was also\nlowest in May, though the timing of highest prokaryotic richness was different than for\nviruses, with the highest prokaryotic richness observed in winter (vs. late summer-\nautumn for viruses)41 (Supplementary Data 1). For the virus richness measurements,\nalthough these patterns were derived from subsampled counts, we observed that the\nfraction of viral diversity captured in winter was much lower than in late summer, as\nindicated by the slopes of rarefaction curves (Supplementary Fig. 2). These differences,\nin part, reflect variations in sequencing depth across samples and suggest that further\nsequencing efforts are required to unravel the viral diversity that persists through late\npolar night.\nLN142: The lower diversity and slower accumulation observed in summer could be a\nresult of higher overall biomass during that season. Few dominant bacteria and viruses\nmay account for the majority of reads at the current sequencing depth (10Gb), so the\ndetection of rare viruses has become relatively difficult. Please consider adding this as a\npotential limitation in the text.\nAfter rarefaction, summer now has higher diversity, along with winter, and it is spring\nthat has the lowest diversity in general. As the reviewer points out, this could be due to\nthe dominance of a few abundant viruses, and is challenged by the relatively low\nsequencing depth, and mention this as a caveat to these aspects.\nLines 194–207: We next investigated how the overall diversity of viruses changes\nacross seasons. As viral diversity was not saturated at any sequencing depth (Fig. 2a),\nwe compared the diversity of vOTUs across months after subsampling to a normalized\ndepth. The diversity of viruses increased following the spring phytoplankton bloom,\nreaching a maximum between late summer-autumn (Aug–Nov) (Fig. 2b), with lowest\nrichness values in March and May. Notably, prokaryotic community richness was also\nlowest in May, though the timing of highest prokaryotic richness was different than for\nviruses, with the highest prokaryotic richness observed in winter (vs. late summer-\nautumn for viruses)41 (Supplementary Data 1). For the virus richness measurements,\nalthough these patterns were derived from subsampled counts, we observed that the\nfraction of viral diversity captured in winter was much lower than in late summer, as\nindicated by the slopes of rarefaction curves (Supplementary Fig. 2). These differences,\nin part, reflect variations in sequencing depth across samples and suggest that further\nsequencing efforts are required to unravel the viral diversity that persists through late\npolar night.\nFig. 3b: Using different shapes for the dots to represent different years could make it\neasier to understand temporal trends across four years.\nFigure 3b has been revised to distinguish years based on different point shapes.\nLN168: I’m uncertain whether the phrase \"shaped by their hosts\" is appropriate, as it's\ndifficult to attribute the observed viral patterns to direct host influence, given the\npotential roles of other prokaryotes.\nLines 275–279: These results suggest that viral communities are primarily correlated with\nprokaryotic community composition which is in turn driven by environmental conditions\nplus biological interactions – leading to a complex network of interdependencies and\nphysicochemical linkages, similar to dynamics in temperate environments7,49.\nAdditionally, I am a bit confused about which pair in Supplementary Data 1 represents\nthe structure of viruses in the main text. Authors could add the values for each variable\nin the main text.\nThe reason we performed multiple sets of Mantel tests was because not all viruses are\npredicted to be viruses of prokaryotes, yet we correlated the entire virus community with\nthe prokaryotic community because it is the most comprehensive analysis to perform.\nBut recognizing some of the viruses may be viruses of eukaryotes, we wanted to see\nalso what the correlations were if we excluded likely viruses of eukaryotes. Thus, we\nalso correlated only the Caudoviricetes with the bacterial community, and the other\nenvironmental parameters for good measure. This tests showed the differences were\nnegligible, thus we focus our main analysis on the entire virus community (most\ncomprehensive), which is column 1 (All vOTUs), and thus we didn’t want to confuse the\nmain text by including so many values, that don’t add much to the main take-away\nmessage (i.e., suitable for the supplement).\nWe now provide a more thorough legend for Supplementary Data 1 (now\nSupplementary Data 3) to explain this result, along with environmental parameters. We\nalso embolden and mark with ** the columns used in the main text.\nAdditionally, we have added a new figure to highlight one the full virus community\ncomposition to bacterial community, along with correlating it with various other\nparameters as a network. Here, we also cite Supplementary Data 3 as supporting\ninformation.\nLines 268–271: At the whole community level, Mantel tests demonstrated the strongest\ncorrelation to prokaryotic community composition (Mantel ρ = 0.632, p = 0.001) followed\nby oxygen (Mantel ρ = 0.371, p = 0.001) and mixed layer depth (MLD) (Mantel ρ =\n0.259, p = 0.001) (Supplementary Data 3, Supplementary Fig. 4).\nLN179: Please add the number of ‘following’ and ‘caused by’.\nWe clarify:\nLines 286–325: The CCM network comprised 5,136 vOTUs and 850 prokaryotic ASVs\nSupplementary Fig. 5a), with directional associations (correlations >0.7) where vOTUs\ndynamics are ‘following’, or are ‘caused’ by, prokaryotic ASV dynamics (akin to Lotka-\nVolterra dynamics50) totalling 15,930.\nLN192-193 & Fig. 4d: The ratio of vOTUs/ASVs in M5 is significantly lower than in other\nmodules, and the value is somehow counterintuitive. Could this be due to that most\nwinter viral nodes being discarded during network preprocessing? A bit further\ndiscussion on potential reasons behind this observation would be helpful.\nThe reviewer makes an astute observation that vOTU to prokaryotic ASVs is\nsignificantly lower in M5. This is the module which contains the most prokaryotes. Yet,\nthis period of low number of vOTU also corresponds to the period at which overall cVCR\n(i.e., the coverage-based virus to cell ratio) is lowest. Thus, yes, it could be due to lower\nsensitivity to viruses during this period. However, it would be difficult to quantify if\nviruses that otherwise might appear in M5 were removed from the network analysis due\nto preprocessing, as we can’t predict which viruses would be in which module without\nrunning new networks with different settings etc, which may be less robust.\nWe now provide a bit further discussion on potential reasons for this observation in the\ndiscussion:\nLines 560-566: It should be considered that factors such as diversity and abundance\ncan impact the total number of viruses detected in modules, and co-occurrence patterns\noverall. M5, for example, which has the lowest vOTU to prokaryotic ASV ratio,\ncorresponds to the period during which vOTUs were the least relatively abundant.\nAdditionally, during this period, vOTU richness was relatively high. The combination of\nthese two factors challenges vOTU detection, and could result in vOTUs being excluded\nduring pre-processing or missed altogether.\nLN265: Are the 58 vOTUs mentioned evolutionarily related (similar taxa)?\nLines 486-490: In total, 53 of the 58 vOTUs were classified as Caudoviricetes (the other\nfive were unclassified), and came from 14 different novel orders (Supplementary Data\n6). Notably, in the entire dataset of ASV to vOTU correlations, only one had a predicted\nhost of Nitrosopumilaceae, which was correlated with an unclassified\nGammaproteobacterium (no family-level prediction).\nFig. 4C: Please ensure the p-values have been corrected for multiple testing.\nThe p-values have been adjusted for multiple testing using the Bonferroni correction\nmethod.\nSupplementary Data: It is not easy to understand the headers in supplementary data 1\nand 2. Please consider using one spreadsheet to include descriptions of the column\nheaders for all tables.\nWe have added descriptions of the column headers below relevant Supplementary Data\ntables to improve clarity.\nReviewer #2 (Remarks to the Author):\nReview Nature Communications. NCOMMS-24-55589\nArctic Ocean virus communities: seasonality, bipolarity, and prokaryotic Interactions\nThis study investigates the seasonality, abundance, diversity, associations with\nprokaryotic hosts, and environmental factors of Artic viruses in the Fram Strait over a\nfour-years period. The authors used 47 long-read metagenomes obtained from the\ncellular size-fraction, which were automatically collected. It was found that viruses\nexhibited annual peaks, and that the composition of the viral community was strongly\ncorrelated with prokaryotic community composition. Additionally, it was analyzed\npresence of the polar virus in metagenome from different latitudes, concluding that virus\nabundance is higher in high-latitude regions of both poles, highlighting adaptations for\ncold environments.\nFirst, I would like to positively emphasize the rigorous mathematical and statistical\nmethods employed by the authors in this work. While seasonal cyclic variability in viral\ncommunities has been observed in prior studies, the extended sampling period here\nprovides a more robust dataset than previous research. The constructed metagenomes\nalone represent a valuable resource database. However, some sections lack sufficient\nbiological interpretation and related discussion. Further clarification from the authors on\nthe following specific points is recommended to enhance the overall quality of the\nmanuscript.\nWe thank the reviewer for their positive evaluation of our statistical methods and\nappreciation of the novelty of the analysis. We appreciate the opportunity to clarify the\npoints below.\nIn general, more quantitative data should be included in the main text and as\nSupplementary Data.\nWe now have incorporated additional quantitative data in the main text (Lines 141–142,\nLines 153, Lines 179–180, Lines 331–333, and Lines 340–341), and have expanded\nSupplementary Data 1 to include comprehensive metadata. This metadata now covers\nenvironmental parameters, project identifiers (bioproject and accession codes), and\ntaxonomy-based read counts and proportions across the domains Archaea, Bacteria,\nand Eukarya.\nTable of Supplementary Data 4, in particular, could be expanded to include specific\ninformation such as the number of viral sequences in each collection, the depth (or\nrange) of the water column from which each sample was obtained, the maximum length\nof the viral raw reads obtained (L239), the estimate of the number of prokaryotic\ngenomes in each metagenome, and, if possible, parameters indicating viral (L142) and\nprokaryotic diversity.\nWe have expanded Supplementary Data 4 (now Supplementary Data 1) to include the\ndepth at which each water sample was collected, the Bioprojects, Accession numbers,\nas well as per-sample data on the number of viral contigs, the length of the longest viral\ncontig, the number of ASVs, and the read counts and proportions of archaea,\neukaryotes, prokaryotes, viruses, and organelles as predicted by Tiara v. 1.0.3.\nAdditionally, a global overview of prokaryotic taxonomy retrieved by ASVs and their\nperiodicity is needed in the Supplementary Data to improve the interpretation of the\nseasonality observed in the viruses (L142, L80).\nWe now add the rarified statistics for the prokaryotic diversity referenced in L142 to our\nSupplementary Data 1. This Prokaryotic diversity from the same time-frame is reported\nin (Priest et al. 2024). We now cite Priest et al. 2024, as well as additional references\nand point directly to the results in these papers, and point out that those results also\nshow that richness is lowest during the spring.\nLines 199–202: Notably, prokaryotic community richness was also lowest in May,\nthough the timing of highest prokaryotic richness was different than for viruses, with the\nhighest prokaryotic richness observed in winter (vs. late summer-autumn for\nviruses)41(Supplementary Data 1).\nDetails regarding the sizes and number or replicas of the 16S rRNA libraries should also\nbe provided as Supplementary Data).\nWe now provide the number of raw sequences and number of replicates (in all cases\nn=1), and accession numbers for the 16S rRNA libraries that accompany the\nmetagenomic samples as part of Supplementary Data 1.\nFurthermore, data on the ASV-vOTUs involved in the six modules, as well as those\nidentified by iPHoP (L181 & L214), must also be included.\nWe have now provided more information about the ASV-vOTUs involved in the six\nmodules, comparing and contrasting the results more explicitly between host predictions\n(iPHoP) and correlations between vOTUs and ASVs (CCM). We add two additional\nfigures, one additional table showing the cases in which they agree, as well as more in-\nline descriptions of the results:\nLines 363–413: Among bacterial families with a high number of predicted vOTUs,\nFlavobacteriaceae and Pelagibacteraceae were the most frequently predicted. In M1,\nFlavobacteriaceae and Akkermansiaceae (Verrucomicrobia) dominated, while\nPelagibacteraceae was less common. Conversely, in M5, Pelagibacteraceae was the\nmost predominant, with fewer Flavobacteriaceae and Akkermansiaceae. Meanwhile,\nvOTUs with Cyanobiaceae as their predicted host were more common in M2–M4\n(Supplementary Fig. 5d). Among the vOTUs correlated with an ASV in the network, the\nmajority of host-virus pairwise correlations were to diverse prokaryotes present in M5,\nwith Pelagibacteraceae (20.3%) and Nitrosopumilaceae (6.0%) being the most common\n(Supplementary Fig. 5e). M1 had the second most bacteria-to-virus connections with\n37.5% of these being to Flavobacteriaceae. Among the families with most pairwise\ncorrelations (i.e., over 250), the overall pattern of host prediction of the viruses and\nmodule membership was similar, with vOTUs with predicted hosts of Flavobacteriaceae\nand Pelagibacteraceae being the dominant host predictions for M1 and M5, respectively\n(Supplementary Fig. 5f). However for individual vOTUs and ASVs, only 3.7% of the\npredictions overlapped at the family level between CCM network correlations and iPHoP\n(i.e., 122 out of 3,280 correlations with associated host predictions). Among these,\nagain, Flavobacteriaceae and Pelagibacteraceae made up 46.7% and 6.6%. For\nFlavobacteriaceae, these vOTUs were diverse, coming from 14 different novel orders\n(via VConTACT52, see Methods), while for Pelagibacteraceae they all came from a\nsingle novel order (Supplementary Data 4).\nTo better grasp the magnitudes involved (as the authors provide statistical data but almost no “raw numbers”), it would be helpful to include precise information in the main text. For example, it would be useful to know how many viral sequences were identified among the nearly 200 000 sequences obtained per metagenome (L96); how many viral sequences represent the 5662 vOTUs identified by the authors; which vOTU was the most abundant and how it varied along time and how many vOTUS are involved in each module (L173-195). Additionally, the 3748 ASVs mentioned in L391 should be cited earlier, at the beginning of the Results sections.\nWe have incorporated the following information, in response to the aforementioned\nquestions, along with other specific details where relevant:\nLines 141—142: “...10,215 ± 1,116 viral reads and 180 ± 29 viral contigs.”\nLines 153: “...which were derived from 7,775 viral contigs that were greater than 10 Kb in\nlength.”\nLines 179–180: “...the most abundant vOTU, a member of the Caudoviricetes, was also\nthe most persistent, occurring in 45 of the 47 sampled time points.”\nLines 331–333: “...M1 containing 1,448 and M2 containing 1,334 vOTUs, respectively.”\nLines 340–341: “...M5 also contained the fewest vOTUs, with only 132 identified.”\nWe now also mention the 3,748 ASVs at the beginning of the Virus-host and\nenvironmental relationships subsection of the Results.\nLines 265–266: “...comprising 3,748 prokaryotic amplicon sequence variants (ASV).”\nSpecific questions:\n1.-In my opinion, the title or the article is inadequate. First, not all viruses form the\ncommunity were considered, as RNA viruses were excluded from the analysis (L422),\nand the authors focused primarily on dsDNA viruses. Importantly, the authors examined\nviruses found exclusively in the cellular fraction; i.e., those actively infecting cells at the\nsampling points or integrated into cellular genomes as lysogens. Additionally, there is\nlimited information on prokaryotic “interactions” per se, as host identification relies\nexclusively on computational predictions. I believe a title emphasizing the seasonality\nand the dynamics of the viral community would better reflect the study´s scope and\nattract a broader audience. Moreover, only viruses present in the cellular fraction are\nconsidered in the analysis, rather than the complete community- an important detail that\nshould be more clearly stated throughout the manuscript.\nWe thank the reviewer for feedback on the title. We are concerned that dsDNA may be\na less accessible qualification to the title. We mention dsDNA in the third sentence of\nthe abstract. We change the title to ‘prokaryotic associations’ rather than prokaryotic\ninteractions, since perhaps associations implies less depth of analysis. We update the\ntitle to:\nArctic Ocean virus communities: seasonality, bipolarity, and prokaryotic\nassociations\nLines 634–637: A further consideration is that our study focused on dsDNA viruses,\nthus not investigating the RNA virome81–85. Investigating RNA viruses, their seasonality\nand impacts on the prokaryotic community is an important future direction.\n2.-The introduction provides quotes related to the ´Artic Atlantification`, but this\ninformation is not directly relevant to the goals or findings of this manuscript. In addition\nto the hypothesis proposed by the authors, L79, which is not entirely novel (reference\n22, 25), I would encourage the authors to present hypotheses in the introduction\nregarding how polar viral communities would fit in the ecological models already\ndescribed for virus-host populations., e.g., ‘Piggyback-the-Winner’ theory, Red-Queen\ndynamics, Constant Diversity Model, etc. (F.R. Valera et al. 2009, B. Knowles et al.\n2016, F. H. Coutinho et al. 2017, H. Alrasheed et al. 2019, J.C. Ignacio-Espinoza et al.\n2020). This would help frame the discussion to better align with fit the Lotka-Volterra\ndynamics observed.\nWe thank the reviewer for helping with the framing.\nWe clarify in the section about cyanophages that cyanobacteria are expected to\nincrease in the Arctic due to Atlantification, so that these observations are related to\nthis:\nLines 84-116: Furthermore, continuous sampling, over the long-term, can discern the\nimpact of 'Arctic Atlantification'34–36: the northward expansion of subarctic habitats\nthrough the Fram Strait — the major connection between Atlantic and Arctic Oceans37.\nHere, polar water outflowing the central Arctic Ocean via the East Greenland Current\n(EGC) meets the West Spitsbergen Current (WSC), transporting warmer Atlantic water\ninto the Arctic Ocean38–40. In the WSC, prokaryotic communities exhibit pronounced\nseasonality, underpinned by changes in photosynthetically active radiation and mixed\nlayer depth41.\nLines 423-425: From the predicted interactions between vOTUs and ASVs, we further\nexplored the diversity and dynamics of putative cyanophages – as Cyanobacteria are of\nparticular interest with respect to changing conditions in Arctic ecosystems5,53,54.\nLines 473–476: Our results reveal that both cyanobacteria and their viruses are present\nin the Arctic, highlighting a need to further understand their interactions for the future\nArctic Ocean, as they are expected to increase due to Atlantification.\nFurthermore, we address the concerns raised about the hypothesis and mention the\nquestion of how polar viruses would fit ecological models described for virus-host\npopulations described elsewhere:\nLines 118–127: Given this, and the tight coupling of viruses and their hosts9–11, we aimed\nto examine the degree that viral populations are seasonally structured in the Arctic. The\nstrong seasonal gradients in the Arctic provide an opportunity to observe host-virus\ndynamics and how they relate to prevailing conceptual models such as the Piggyback-\nthe-winner42–44, Constant-Diversity45,46, and Red-Queen28,47 hypotheses. We examined\nthe diversity and seasonality of dsDNA virus communities in the Arctic Fram Strait over\nfour complete annual cycles at roughly monthly resolution. This unprecedented time-\nseries of virus communities in the Arctic demonstrates considerable seasonality in viral\ndiversity, their association with environmental conditions and potential microbial hosts,\nand their distribution across the global ocean.\nWe also frame the results within the context of these concepts:\nLines 698–715: Our results indicated a low percentage of predicted prophages without a\nseasonal enrichment, even in the cellular size fraction studied. This suggests that the\n“Piggyback-the-winner” scenario is neither very prevalent, nor seasonally variable,\nthough again this may be currently limited by computational predictions, and other biases\nlike lytic viruses occurring within cells in more than one copy, and thus we may\nunderestimate prophage prevalence. We found that most viruses are part of seasonally\nvariable modules, and that often there is no similarity between opposing seasons (e.g.,\nsix months apart), at least at the depth sequenced at in the investigated cellular size\nfraction. Hence rather than having a ‘Constant-Diversity’45,46, Arctic virus communities are\nundergoing substantial seasonal change. Likewise, the prokaryotic community is also\nhighly seasonal, thus the viruses are generally following and strongly correlated with host\nabundances. Thus, unlike how viral ecological models are sometimes idealized in a\n‘chemostat-like’ setting100, viral communities' strong seasonality impact how such models\ncan be of use in the Arctic Ocean (such as in the Red-Queen hypothesis28,47). Arctic virus\ncommunities, with their strong seasonality where the many microorganisms are not\npresent year-round (or at least often not detected), may be impacted by how these\ndynamics play out elsewhere, or it’s likely these dynamics occur primarily within seasons.\nIn each case, these dynamics could be further evaluated via more spatially-resolved\ntemporal sampling.\n3. There are two points that need to be clarified regarding the methodology used:\n(i) Importantly, the authors claim that the lytic lifestyle was predominant over lysogeny.\nAlthough I find value in the analysis conducted by the authors, I have a primary concern\nregarding the number of lysogens detected. Beyond bioinformatic predictions, the\nauthors lack direct evidence to fully support this conclusion. Many lysogens do not have\nrecognizable integrases, which is a general challenge that needs to be addressed (out\nin the context of this work). In my view, viruses replicating inside the cells (i.e., lytic\nviruses) would outnumber lysogens in term of copies within the cellular metagenome,\nmaking them easier to assemble with the long-reads used, compared to lysogens.\nTherefore, a strong bias could obscure the results.\nWe appreciate the points about lysogens potentially not being recognized due to\nrecognition limitations and about copy numbers of lytic infections vs. lysogens.\nState-of-the-art lysogen/prophage predictors (VIBRANT, CheckV, and geNomad) use a\nvariety of techniques to identify lysogens, most typically looking for areas of contigs that\nare enriched in viral genes surrounded by proteins that appear to be of cellular origin.\nOriginally, we used VIBRANT to predict lysogens, but upon critically evaluating the point\nby this reviewer, we realized our workflow for predicting lysogens was flawed because\nwe applied VIBRANT to contigs from which CheckV had already been used to ‘trim off’\nhost regions of prophages. We have now rectified this issue and replotted the data for\nFig. 1c, and Fig. 5c, though the main result, that lysogens/prophages are rare, remains\nthe same.\nWe have added the following sentences to aid in the interpretation and describe the\nlimitations.\nLines 184–192: Overall, lytic lifestyle was predicted to be predominant, with an average\nof 94.0 ± 0.4% of viruses, a pattern that was invariable across seasons (Fig. 1c). The\nSouthern Ocean featured higher rates of lysogeny during periods of low production25.\nAs for the latter point, it is important to point out that some lysogens may not be\nrecognized due to missing integrases or other factors, and, vOTUs may not represent\nthe full prophage, which could then be misinterpreted as lytic rather than integrated.\nAdditionally, lytic viruses may outnumber lysogens in terms of copies within the cellular\nmetagenome, making them easier to detect and assemble, which could bias the\nrecovery or obscure lysogens.\n(ii) Cut-off used to define the abundances of vOTUs. The authors use only 25% of the\nviral genome length, whereas typical values range from 50 to 70% at a minimum.\nWe appreciate this point by the reviewer, and we had discussed this amongst co-\nauthors previously. We evaluated different cut-offs and decided for 25% as a trade-off\nbetween confident genome coverage and sensitivity, and based our decision as follows.\nFirst, the reads we are using are long-reads with an average length of 5,435 ± 405.\nThus, this can be considered when defining the minimum coverage, especially as it\nrelates to the more common short-read metagenomes. For short reads, unlike long-\nreads, the 150 bp reads could map to multiple conserved regions across many genes of\na contig, which may be more problematic and predictable in terms considering that most\ngenes have highly conserved motifs and less conserved motifs (circumvented by long-\nreads). Additionally, considering that the shortest contig that we decided to use was\n10kb (which was a conservative cut-off, itself), on average, the minimum coverage\nwould be 50%, though some may be less than that.\nNevertheless, we considered how different cut-offs can influence our results, especially\nthe main result of seasonality within the dsDNA virus community. We found that the\nseasonality is similar across a variety of cut-offs from 0, 0.25, 0.5, and 0.75. The\ndifference is that the higher the coverage required, the less resolution for the number of\nviruses that can be observed in any given sample. Given the fact that we had long-\nreads and that the main result of seasonality still strongly held at various cutoffs, we felt\ncomfortable with 0.25 which would enable more resolution of viruses across months,\neven where they may be more rare.\nWe have updated the manuscript with a supplementary figure showing how the various\ncut-offs do not impact the result of seasonality.\nLines 257–259: Notably, the seasonality observed in virus community composition was\nconsistent regardless of minimum coverage value utilized for detected viral presence,\nthough similarity was less at higher thresholds (Supplementary Fig. 3).\n4. The authors focus on Cyanobacteria and the dynamics of putative cyanophages\n(L232-257) due to their relevance to changing conditions in Arctic ecosystems (but it is\nunknown in which module are founded), along with another nice example of the\nNitrososphaerales. In my opinion, there is a lack of detailed information regarding the\nlargest group, which is represented in several modules and with multiple identified virus-\nhost pairs, the Pelagibacteraceae. Providing information on how the Pelagibacteraceae\nASVs are distributed across the different modules, along with the assignment vOTUs,\nwould better complement this section of the manuscript.\nWe have now gone into more detail on the Pelagibacteraceae dynamics, as well as\nother dominant hosts, the Flavobacteriaceae; two groups which dominated the CCM\nnetworks. We provide a new chord diagram that indicates that most of the ASV-to-vOTU\ncorrelations occur to prokaryotes in M1 and M5; these patterns are similar to the host\npredictions of the vOTUs in these modules. We also provide information specifically\nabout the vOTUs that are both correlated with and have host predictions for\nPelagibacteraceae (Supplementary Data 4).\nLines 362–418: In total, 42.5.2% of vOTUs received a host prediction by iPHoP. Among\nbacterial families with a high number of predicted vOTUs, Flavobacteriaceae and\nPelagibacteraceae were the most frequently predicted. In M1, Flavobacteriaceae and\nAkkermansiaceae (Verrucomicrobia) dominated, while Pelagibacteraceae was less\ncommon. Conversely, in M5, Pelagibacteraceae was the most predominant, with fewer\nFlavobacteriaceae and Akkermansiaceae. Meanwhile, vOTUs with Cyanobiaceae as\ntheir predicted host were more common in M2–M4 (Supplementary Fig. 5d). Among the\nvOTUs correlated with an ASV in the network, the majority of host-virus pairwise\ncorrelations were to diverse prokaryotes present in M5, with Pelagibacteraceae (20.3%)\nand Nitrosopumilaceae (6.0%) being the most common (Supplementary Fig. 5e). M1\nhad the second most bacteria-to-virus connections with 37.5% of these being to\nFlavobacteriaceae. Among the families with most pairwise correlations (i.e., over 250),\nthe overall pattern of host prediction of the viruses and module membership was similar,\nwith vOTUs with predicted hosts of Flavobacteriaceae and Pelagibacteraceae being the\ndominant host predictions for M1 and M5, respectively (Supplementary Fig. 5f).\nHowever for individual vOTUs and ASVs, only 3.7% of the predictions overlapped at the\nfamily level between CCM network correlations and iPHoP (i.e., 122 out of 3,280\ncorrelations with associated host predictions). Among these, again, Flavobacteriaceae\nand Pelagibacteraceae made up 46.7% and 6.6%. For Flavobacteriaceae, these vOTUs\nwere diverse, coming from 14 different novel orders (via VConTACT52, see Methods),\nwhile for Pelagibacteraceae they all came from a single novel order (Supplementary\nData 4). However, overall, the low number of matching CCM correlations and host\npredictions at the vOTU-ASV level demonstrates the challenge to discern host-virus\nrelationships in diverse and complex ecosystems where a variety of technical and\nbiological challenges may complicate such analyses (see Discussion). Nonetheless, the\nresults provide valuable indications for particular lineages at an overall module level.\nOther details:\nJust as a suggestion, did the authors consider whether Supplementary Figure 3c could\nappear as a main figure (perhaps in place of Figure 4a)?\nWe have now moved it to the main text and moved the original Figure 4a to the\nsupplement.\nL63. Please, clarify. Sandaa et al 2018, (reference 22) found that (cited textually) “Viral\ndiversity and virus-to-prokaryote ratios (VPRs) dropped sharply at the commencement\nof the spring bloom but increased across the season, ultimately achieving the highest\nlevels during the winter season. These findings suggest that viral lysis may be an\nimportant process during the polar winter, when productivity is low”. Therefore, this\ncontradicts what is stated in L59, 60.\nWe reviewed these papers again, and have revised the manuscript accordingly:\nLines 69–71: One year-round study showed strong seasonal variation among virus\ncommunities, with sharp decreases in virus-to-prokaryote ratios at the onset of the\nspring bloom and highest ratios in winter22 Others have reported increases during\nspring-summer, and highest abundances of viruses in winter23,24.\nWe also add a line to the discussion:\nLines 611–615: Our observation of strong virus seasonality in the cellular size fraction\ncontrasts with some of those of free viruses via fluorescence microscopy. This seeming\ncontradiction may indicate a difference in the dynamics of infection and host-\nassociation, or loss factors of free-living viruses71,72, for example photo-degradation of\nfree viral particles73,74.\nL93-94. Did the authors retrieve eukaryotic DNA in their sequencing? Since no pre-\nfiltering was applied, there is a high possibility of this occurring, yet it is not mentioned in\nthe manuscript. If present, could this affect the cVPR calculations?\nWe have now evaluated the amount of eukaryotic DNA in the samples, and find that\neukaryotic reads made up 1.5% ± 0.53% of reads in the metagenomic samples, vs\n72.2% ± 15.3% of prokaryotes, indicating that eukaryotic sequences made up a small\nminority of the samples, even without pre-filtration.\nWe update the Results:\nLines 142–145: Notably, despite not pre-filtering to remove eukaryotes, the large\nmajority of reads in the dataset were prokaryotic (72.2% ± 15.3%), with eukaryotes and\nviruses making up smaller proportions (1.5% ± 0.53%, and 0.4% ± 0.5.%, respectively).\nDiscussion:\nLines 609–611: This increase was consistent regardless of the normalization used,\nnamely either coverage-based virus to cell ratio (our main metric), or also coverage per\ngigabase pair (which would also account for eukaryotic DNA).\nAnd Methods:\nLines 849–852: To gain an insight into the composition of the metagenome samples\noverall in terms of prokaryotic and eukaryotic reads, we employed Tiara v.1.0.3121, a\nmachine-learning tool, to assign domain-level classifications to the raw reads. This\nrevealed that the large majority of reads across the metagenomes, >95%, were of\nprokaryotic origin.\nThe reviewer has made a good point in asking how the cVPR can affect the cVPR\nvalues. In fact, upon further consideration, we realized that our cVPR metric actually\ndoes include eukaryotes within the calculation, as it calculates the average coverage of\n16 universal ribosomal proteins found across all three domains.\nWe have edited the methods to clarify this point throughout the manuscript, and\nrenamed the metric cVCR (coverage-based Virus to Cell Ratio), from cVPR (coverage-\nbased Virus to Prokaryote Ratio). However, again, the eukaryotes are making up a\nsmall amount of the ratio itself since they only make up ~1.4% of total reads.\nWe update the methods:\nLines 842–844: Second, we divided the mean coverage by the estimated number of\ncellular genomes in each metagenome, as predicted based on the mean coverage across\n16 universal single-copy ribosomal proteins119,120.\nL98-99. The program CheckV is mainly used to confirm the nature of the viral sequences,\nnot to make predictions per se. Therefore, remove from this form the line, and then the\ntwo approaches mentioned in L99 that refers the authors, will be correct.\nWe appreciate this point by the reviewer and agree. We were following a published\nVirSorter2-based standard-operating-procedure for virus prediction that includes\nCheckV. We now specify this more fully and cite the SOP that we based our analysis\non: https://www.protocols.io/view/viral-sequence-identification-sop-with-virsorter2-\n5qpvoyqebg4o/v3\nTo avoid confusion, we remove the tools from the results and cite the methods where\nwe can describe this more fully.\nLines 139–141: Viral sequences were predicted on both assembled contigs and raw\nreads (see Virus Prediction, Methods) and the overall number of predicted viral\nsequences was concordant between the two approaches (Supplementary Fig. 1).\nL107. “95% sequence clusters”. Please modify to indicate that this refers to 95% ID.\nWe edited, Line 152: “95% nucleotide identity clusters”\nL109. Please note that it is recommended to follow the updated nomenclature for\nphages. The bacteriophage families Myoviridae, Podoviridae, and Siphoviridae have\nbeen reclassified under a new taxonomic system (International Committee on\nTaxonomy of Viruses (ICTV), 2021). These are now part of the order Caudoviricetes,\nwhich encompasses tailed bacteriophages. Myoviridae members are often reclassified\nin families like Myoviricetes, while Podoviridae members generally fall under\nAutographiviridae, and Siphoviridae have been split across families like Drexlerviridae\nand Demerecviridae, among others.\nIndeed the ICTV is ratified/accepted, but the taxonomy and tools and strategies for\nemploying it are still being developed, especially for uncultivated viruses. This includes\nongoing work within the Caudoviricetes, the Class of 98% of the viruses in our\nmanuscript. There is only one tool, geNomad, that provides high-quality classifications\nbased on the ICTV at this point, but this tool did not more finely resolve most of our\nCaudoviricetes sequences beyond the level of Caudoviricetes. vConTACT3 could\nprovide more information regarding protein network structure (and theoretically,\ntaxonomy), but the results were not clearly easy to incorporate, with 91.3% coming from\nnovel orders and only 0.4% receiving a taxonomy at the order level (Crassvirales).\nNevertheless, we provide the vConTACT3 taxonomies as part of Supplementary Data\n2, and cite them in cases where we talk about specific vOTUs.\nWe now apply geNomad to update our taxonomy to the ICTV taxonomy, which gives\noverall very similar results to our initial submission, i.e., if one were to\n‘add-up’ the different out-dated morphological families (Myo-,Podo-,Sipho-) to\n‘Caudoviricetes’. If we had focused on differences within Caudoviricetes (at a finer\ntaxonomy), then this update would potentially have caused more issues, but as we did\nnot, it's not such a big problem (as also mentioned by the reviewer).\nWe update the methods, and recalculate values for class level variability we discussed,\nas well as regenerate figures to include the new taxonomy.\nOverall, the main conclusions and patterns remain unchanged.\nL111. Phycodnaviridae continues to be cited in in L122, L199 and Figure 3, so please,\nupdate this sentence accordingly.\nWe have removed the part of the sentence about not focusing on Phycodnaviridae (now\nincluded in ‘Megaviricetes’ in the paper).\nL166. Please define `PAR’ as it is the first time it appears in the text. Also, what\nparameter exactly is the “polar water fraction” in Figure 3b? Additionally, “iceDist” and\n“icePast” are not mentioned in the manuscript (also in Figure 4).\nWe have defined “PAR” (Photosynthetically Active Radiation) (Lines 273–274) as\nsuggested.\nPolar water fraction is the proportion of polar water in the sample and was derived from\nsalinity and temperature. We now describe this in the Methods (Lines 763–765).\nAdditionally, we have removed iceDist and icePast from the manuscript, as these\nparameters are less relevant in the context of the West Spitsbergen Current (WSC).\nL287. Is this in disagreement with the statement in L280? First, authors must clarify\nwhether the vOTUs detected in surface waters are different from those detected in deep\nwaters. If Fram Strait viruses are also abundant in deep ocean samples, does this imply\nthat similar hosts are more abundant at these depths? Additionally, which depth range\nis being referred to? This information needs to be either depicted in the figure or stated\nin the Supplementary Data 3. It is generally understood that deep and surface\nprokaryotic communities differ… then, could mixed waters be influencing this\nobservation (although the authors specifically mention that these “mixed” samples were\nnot considered in the analysis)? The authors should provide a more detailed discussion\non this point (L291-298).\nWe have revised this section for clarity. We didn’t mean to imply that the Fram Strait\nviruses were more abundant in the deep samples, only that the proportion of samples in\nwhich they were detected was higher in deep samples. In terms of abundance, Fram\nStrait viruses are more abundant in the surface waters, but often completely absent\nfrom surface ocean samples (for example in the low latitudes < 30 N and S).\nWe now also clarify that the communities of viruses present in surface vs. deep do differ\nsignificantly in terms of ‘who is where’ which we now point out via an adonis test which\ncan detect differences in community structure between two groups.\nWe clarify these points as follows:\nLines 537–548: Fram Strait viruses were also common in the deep ocean (here,\noperationally defined as sampling depth of >200m). In particular, they were detected in\n87% (338 of 390) of deep global samples. Like the surface waters, deep-water viruses\nwere more prevalent in northern and southern higher latitudes (Fig. 7b), though in deep\nsamples the viruses were more commonly detected. Overall, however, when detected,\ntheir relative abundances were less in deep than in shallow samples. Notably, there was\nsignificant variation in community structure of the detected Fram Strait viruses between\nsurface and deep samples (adonis2, Sum of Squares = 38.68, R² = 0.11448, p < 0.001),\nwhich aligns with established differences in microbial communities in the deep vs.\nsurface ocean, aligning with past results of virus biogeography14. Future work could\nconsider the differences in the total virus communities between surface and deep, as\nthe current analysis utilizes only mapping to the Fram Strait viruses.\nWe also add to the methods to describe how the adonis2 results were calculated:\nLines 980–985: For evaluating statistical differences in community similarity between\nsurface (<200m) and deep (>200m) samples, we first removed samples that had less\nthan 0.0001 cVGB (n=769 samples remaining, including 435 and 334, shallow and deep\nsamples, respectively) and vOTUs that had less than 0.0001 summed cVGB across all\nsamples (n=4,693 vOTUs remaining). Then, we ran adonis2122,148 with 999\npermutations based on the surface and deep grouping.\nL310-L317. Please, provide more details on the changes observed in the amino acid\ncomposition of the viruses from Fram and GOV2.0 datasets versus those found in in\nnorthern or southern lower latitudes. The authors provide statistics, but there is no\npercentage or other parameter that can be quickly appreciated.\nWe now add an additional multi-panel supplementary figure that shows each protein\nparameter’s relationship to temperature, which is, we believe, more relevant than\nlatitude itself for this type of analysis.\nThese panels reflect the correlations or lack of correlation for the various parameters\nshown in the main text and can be quickly appreciated, as requested by the reviewer,\nand we cite the figure as follows in the main text:\nLines 584–586: Furthermore, examining the average values for these protein\nparameters according to sample and temperature reflected the patterns observed via\nother statistical analyses (Supplementary Fig. 8, Supplementary Fig. 9), reinforcing our\nobservations.\nWe also add the script used for generation of these plots to the github repository.\nWe add the following Figure legend:\nLines 2050–2053: S. Fig. 8. NMDS plots (left) showing the relationships between\nsamples and environmental parameters and Spearman correlations (right) of\nenvironmental parameters to amino acid traits for a, Caudoviricetes; b, Megaviricetes;\nand c, Others (non-Caudoviricetes, non-Megaviricetes).\nLines 2055–2059: S. Fig. 9. Correlation of protein parameters with temperature. Each\npanel shows the relationship of temperature and the mean of the specified protein\nparameters calculated for each sample from this study for both the Fram Strait (blue)\nand GOV2.0 (red) datasets. The blue line represents a trend line using LOESS curve\nfitting with a 95% confidence interval (grey).\nL331. The quantity of free viruses cited here were obtained in 2014 during in partial\ntransect compared to the one followed in this work. Although due to the periodicity\nobserved, it is likely that those numbers were similar, the authors could support this by\nusing metaviromes and metagenomes recovered from the same sample and check the\nvOTUs abundance recovered in both.\nWe appreciate the reviewer’s idea here. The very nice 2014 study focused on amplicon\nsequencing of g23 and mcp proteins and cell and virus enumeration in free living size\nfractions. It is a very informative study, but does not allow for easy comparison to our\nstudy (or calculating cVCR), because there are no metagenomes (to the best of our\nknowledge). We searched for updated metagenomes from this study, but did not find\nany. Nevertheless, we have clarified that our cVCR results do not strictly correlate with\ntheir count data, as mentioned elsewhere, and again that our alpha diversity\nmeasurements correspond to theirs (which we further clarify is strongly correlated with\nprokaryotic alpha diversity).\nIn reference to this paper from 2014, we added new text at Lines 69–71 as indicated\nabove.\nL438. RaFAH results are not shown anywhere in the text. The iPHoP results incorporate\nRaFAH, which helps improve phage-host association predictions. Please, correct the\nsentence.\nThis sentence has been corrected as suggested:\nLines 807–808: When multiple host predictions were available, we selected the\nprediction with the highest confidence score.\nL452. Is reference 88 correctly place here? Please cite the original source where it is\ndescribed.\nIn reference 88 the method we utilized is described more thoroughly, and is the first\ndocumentation of that method. We now also additionally cite the paper from which the\nuniversal single copy genes (in this case, a set of 16 ribosomal proteins) were\ndeveloped (Hug et al., 2016).\nL536. The metagenomes chosen to map the viruses were those obtained using a\nprokaryotic size fraction of 0.2-3 microns. However, authors do not cite to use any pre-\nfiltering, so, collections with larger sizes could also be suitable. Please, revise the\nsentence accordingly, as these metagenomes were selected because the authors\nfocused on prokaryotic dsDNA virus relationships.\nLines 964–969: Only DNA samples from discrete depths (e.g., not composites of the\nmixed layers or similar) and the prokaryotic size-fraction (0.2–3 µm) were used to allow\ncomparison with our dataset. Multiple other size fractions are available, but the 0.2-3\nµm was utilized because of our focus on prokaryotic dsDNA viruses, its comparability to\nthe present Fram Strait virus dataset, because it is among the most common range\nsampled (i.e., many samples available).\nFigure 4a. Please, modify this figure to better differentiate viruses from ASVs.\nAs suggested by Reviewer 2, we have moved this figure to the supplement where it is\nnow approximately 3x as large as previous. We distinguish vOTUs from ASVs better by\naltering the border width around their symbols. Together, we hope that this helps\ndifferentiate vOTUs and ASVs.\nFigure 5a. psbA gene in italics.\nChanged throughout.\nPlease, enlarge the figures for improved readability. At their current size, it is\nchallenging to distinguish finer details.\nWe have increased figure size wherever possible to help with distinguishing fine details,\nand increased the font sizes to improve readability throughout.\nFigure S3c. May the authors order the microbial family names accordingly to a\ntaxonomic relationship? It would make more sense than organizing them alphabetically.\nWe have modified the figure accordingly.\nReviewer #2 (Remarks on code availability):\nWe thank the reviewer for checking the code. We now upload relevant inputs to our\ngithub: https://github.com/alyzzabc/fram_strait_viruses_2016-\n2020/tree/main/input_tables\nReviewer #3 (Remarks to the Author):\nThe authors present a well-written paper about a novel dataset that describes a rare 4-\nyear time series of viral metagenomes form the Arctic. They introduce the broad\nseasonal patterns of viruses they observed and add extra analysis by comparing the\nviruses they found with other publicly available datasets of that type. They use a fairly\nstraight forward analysis with appropriate and common pipelines and tools.\nWe appreciate the positive evaluation of our submission and recognition of the value of\nthe dataset.\nHowever, there is one noticeable exception: the authors look at a mix of morphotype\nclassification (podo-, myo- and sipho-viridae) for bacteriophages and genetic based\ntaxonomic groups for eukaryote viruses (Phycodnaviridae, Mimiviridae). This seems like\nan outdated approach since the ICTV reported that they are not monophyletic,\nphylogenetically solid groups (See e.g. “Taxonomy of prokaryotic viruses: 2018-2019\nupdate from the ICTV Bacterial and Archaeal Viruses Subcommittee”). A direct\ncomparison between those different types of classification groups not ideal, and there is\nno clear reason why the more recent classification has not been chosen. The main\nfocus of the paper is not on that aspect, so it is not an unfixable problem. But please\nupdate the analysis according to the most recent ICTV guidelines, or at least justify why\nthe chosen approach is still valid (I will leave that up to the editor). The authors should\nat minimum acknowledge this issue in the text and explain why this comparison was\nchosen, rather than the more up-to-date classification established in recent years.\nWe have replied to this point at length in response to Reviewer 2, and update the\nmanuscript accordingly throughout. Overall, it required re-generating multiple figures\nand text values, but has not significantly changed the interpretation of our paper since\nas most of our main results were framed by the taxonomy as a helpful descriptor, but\noverall were not a major differentiating factor (i.e., within Caudoviricetes).\nIt is crucial that the authors carefully go through the entire manuscript and make\nabsolutely sure that they provide a citation for every software tool, package and\ndatabase and whatever else they use. Studies like this one would not be possible\nwithout all the public databases and free software tools developed by other scientists, it\nis only good practice to give everyone credit for their work.\nWe agree with the reviewer on this, and apologize for missing some references\n(highlighted below). In the original version, yes, we had cited many softwares at only the\nfirst mention. We now cite also at any repeated mentions, as well as those highlighted\nbelow that were missed.\nMinor comments:\n22-23: Perhaps specify that the paper is about viruses of microbes.\nAmended as suggested.\n26-28: Wording could be clearer.\nWe simplify:\nLines 27–28: Among 5,662 vOTUs, 98% and 2% were Caudoviricetes and\nMegaviricetes, respectively.\n62-63: A new paper has been released after the authors would have submitted their\nmanuscript, but there is a viral metagenomic time-series that might be relevant to\nconsider here. https://academic.oup.com/ismej/advance-\narticle/doi/10.1093/ismejo/wrae216/7833430\nWe now cite this paper in the Introduction and Discussion where we compare our\nresults with those from more temperature environments:\nLines 77–81: However, due to the challenges of continuous sampling in the polar\nregions across multiple years, the degree of seasonality among polar viruses – here,\nmeaning annually repeating patterns of populations and communities across the same\nse\nseasons of different years26 – and the potential ecological implications remain to be\ndiscovered. Such seasonality has been previously observed in marine temperate\nenvironments27–29.\nLines 641–642: However, we observed low similarity in opposing seasons compared to\nstudies from more temperate areas27-29.\n78-80: It is known that microbes in the polar regions are highly seasonal and, obviously,\nvirus abundance is directly linked to host abundance, a metric the authors use as well.\nGiven that seasonality of viral communities has been shown in other regions and has\nbeen speculated for polar viruses as well, this is not a good hypothesis in my opinion.\nThe authors present a very important and interesting dataset, which as far as I know is\nthe first virome time series from the Arctic Ocean. They should emphasize the novelty of\nthis dataset, rather than the general concept of viral seasonality.\nWe address this point immediately below, editing the text in accordance with the\nsuggestion.\nAn additional side point: the manuscript only provides data from the Fram Strait, so the\nhypothesis should be about Arctic seasonality, not polar seasonality, in any case. It’s\nreasonable to speculate similar patterns for the Southern Ocean, but the seasonal\ndataset only really covers the Arctic.\nWe edit this paragraph accordingly.\nLines 118–127: Given this, and the tight coupling of viruses and their hosts9–11, we aimed\nto examine the degree that viral populations are seasonally structured in the Arctic. The\nstrong seasonal gradients in the Arctic provide an opportunity to observe host-virus\ndynamics and how they relate to prevailing conceptual models such as the Piggyback-\nthe-winner42–44, Constant-Diversity45,46, and Red-Queen28,47 hypotheses. We examined\nthe diversity and seasonality of dsDNA virus communities in the Arctic Fram Strait over\nfour complete annual cycles at roughly monthly resolution. This unprecedented time-\nseries of virus communities in the Arctic demonstrates considerable seasonality in viral\ndiversity, their association with environmental conditions and potential microbial hosts,\nand their distribution across the global ocean.\n90-91: What was the depths these samples were collected from? Please clarify.\nLines 131–133: We examined 47 long-read metagenomes from samples collected at\nnear-monthly resolution over a four-year period (Aug 2016 – Aug 2020) in the WSC\n(Fig. 1a, Supplementary Data 1), from an average depth of 29 m.\nLines 737–739: Moored Remote Access Samplers (RAS; McLane) autonomously\ncollected and fixed seawater from an average depth of 29 m in the eastern and western\nFram Strait (Fig. 1) at weekly to fortnightly intervals between 2016–202041,58.\n117: s missing for ratios?\nFixed.\n166: Please introduce the meaning of abbreviation PAR.\nThe abbreviation PAR has been defined as suggested (Lines 273–274).\n162-171: Please provide a value or a figure for this. The stats are buried in a difficult-to-\nread supplementary table, but this is a fairly important point and the following paragraph\nis building on that point. it should be made clear to the reader how significant that\ncorrelation is compared to oxygen and MLD. The authors provide a main figure (3) and\na percentage to show the results for the less important physiochemical drivers, but not\nfor the correlation with their hosts. The text and figures should reflect the more\nimportant aspects of results first.\nWe have included the Mantel ρ statistics in the text Lines 268-271, and added\nSupplementary Fig. 4 which shows the complete Spearman correlations together with\nthe Mantel test.\nWe cite this figure in the main text as follows:\nLines 268–271: At the whole community level, Mantel tests demonstrated the strongest\ncorrelation to prokaryotic community composition (Mantel ρ = 0.632, p = 0.001) followed\nby oxygen (Mantel ρ = 0.371, p = 0.001) and mixed layer depth (MLD) (Mantel ρ =\n0.259, p = 0.001) (Supplementary Data 3, Supplementary Fig. 4).\n241-242: This is not surprising since the authors group the viruses based on the\noutdated morphological classification, which is polyphyletic.\nTo improve clarity (as mentioned elsewhere) we have removed this figure panel and the\nassociated text, and replaced it with a simpler figure.\n287: Define what depth is considered deep ocean. Based on context, > 200m depth.\nWe clarify:\nLines 537–538: (here, operationally, sampling depth of greater than than 200 m)\n467: Please cite the authors of the vegan and ggplot2 packages and acknowledge the\nuse of R.\nThe citations have been added as requested for the use of R (Lines 917–918), vegan\n(Line 862, Line 918, Line 1026, Line 1029), and ggplot2 (Line 864, Line 1033).\n501: Please reference the original Cytoscape paper.\nFor the original Cytoscape paper, Line 911.\n518: Please cite the original vConTACT2 paper.\nFor the original vConTACT2 paper, Line 927.\n519: Please cite the hmmer paper.\nFor hmmer, Line 928.\n520: Please add citation for UniRef.\nFor UniRef, Line 930.\n532: Please add citations for Tara Oceans, Malaspina Expedition and Bio-GO-SHIP.\nThey are hidden in a supplementary table as a link only, but should be properly cited in\nthe text itself too. I will stop checking if the authors have referenced every tool,\ndatabase and package, since they tend to only add a citation upon first mention, if at all,\nso it’s tedious and easy to miss.\nWe now cite these main large-scale studies in the main text, and the rest we used data\nfrom are now referenced in the Methods:\nLines 499–502: ...we assessed the distribution of Fram Strait viruses across the global\noceans through metagenomic datasets from various large-scale sampling campaigns,\nsuch as Malaspina, Tara Oceans and Bio-GO-SHIP14,65–69 (Supplementary Data 7).\nLines 961–962: Metagenomes from ten oceanic datasets 14,65–69,139–144, including Tara\nOceans, Malaspina, and Bio-GO-SHIP (Supplementary Data 7), were downloaded from\nNCBI,....\nReviewer #3 (Remarks on code availability):\nThe codes seem fine, there is enough instructions and explanations.\nThank you for the positive assessment of our code availability.\n###### Additional edits #######\nWe also updated a few methodological aspects related to the network analysis: We\nadjusted the Git Repository link (https://gitlab.com/qtb-hhu/marine/publications/framphages2024), cited an implemented GUI-based framework\nfor community dynamics using network analysis called Otter (DOI\n10.5281/zenodo.13840702 (https://zenodo.org/records/13840807), and a recent paper\ndescribing the method (Beyond blooms: the winter ecosystem reset determines\nmicroeukaryotic community dynamics in the Fram Strait.\nhttps://www.nature.com/articles/s43247-024-01782-0)"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I’d like to thank the authors for their efforts in addressing my comments and improving the manuscript. After reading the rebuttal letter and the revised manuscript, I find that most of the responses are adequate, and the revisions are well-explained and feasible. The improvements, such as clarifying genome completeness and expanding discussions on arctic adaptations at gene level, are appreciated. Overall, the manuscript has largely improved.\nI still have a few remaining suggestions to enhance robustness, transparency, and clarity of this paper.\n# 1. First of all, I have to say sorry that I made a mistake in my first revision about the total sequencing depth. I originally thought it was 10 Gb per sample, but, in fact, it is only 1 Gb (196,489 reads × 5,435 bp per read). This low sequencing depth makes the quantification of viruses more challenging than I expected, especially given that viral diversity does not reach saturation in the dataset.\nThe updated results indicate that viral diversity peaks in late summer to autumn, which is the somehow opposite of what was stated in the initial submission. This is a big shift in the overall viral seasonal dynamics.\nMy question is:\nThe error bars in Figure 2B are very large, making it difficult to compare diversity across months with confidence. Would it be possible to restrict the analysis to only samples with adequate viral reads when calculating diversity? Alternatively, it might be helpful to present a fully expanded diversity/richness plot for each year, similar to what was done for cVCR in Fig.1 and Fig.S1.\nAlso, I have a question regarding the subsampling approach. As described, subsampling was performed with “100 iterations of subsampling the vOTU count table at a range of different depths, spanning from 25 up to 16,000 counts at 50-count intervals.”\nDoes the 25% mapping breadth threshold still apply in this subsampling process? Specifically, if removing just one read causes the mapping coverage to drop below 25%, would that vOTU still be considered? The method here is important, particularly to richness estimation, where rare viruses may be affected. I suggest that the authors provide more details on how subsampling was performed with the mapping threshold.\nIn addition, the usage of slope and direct subsampling ponits to compare diversity/richness across sequencing depths might be improved by using non-linear functions. For example, combining rarefaction and extrapolation using the iNEXT R-package could provide estimates of diversity/richness accounting for undetected viruses beyond observed data.\nGiven the big change in diversity results, I wonder whether another key metric, cVCR, and its associated results are also influenced by sequencing depth and normalization (subsampling). Since cVCR relies on the coverage of both viral and cellular genomes, low and varied sequencing depth could potentially bias cVCR values as well. So subsampling may potentially introduce some shifts in results, particularly in seasonal trends presented in Fig 1.\nI believe the observed cVCR-based abundance seasonal pattern is biologically meaningful, but I would prefer authors to reinforce the conclusions by addressing these potential biases.\nMy suggestion is to justify the use of cVCR by testing whether cVCR values remain stable when reads are subsampled to the same sequencing depth.\n# 2. As far as I understand, the parameter \"-B 1000\" used in the phylogenetic analysis, as written in the Methods section, refers to UFBoot. This should be explicitly stated in the Fig.5 legend, because different from traditional values, UFBoot values are only considered reliable if they are ≥ 95% (http://www.iqtree.org/doc/Frequently-Asked-Questions).\nAdditionally, it would be greatly appreciated if the authors could provide the full psb phylogenetic tree and alignment as supplementary files to allow readers to evaluate evolutionary relationships.\n#3. The paper [https://pubmed.ncbi.nlm.nih.gov/39448562/] provides a background on polar viruses in a time-series context. The authors may consider referencing it to strengthen the introduction and discussion.\n#4. I need some clarification for the newly added values in LN 145 because “On average, each sample harbored 10,215 ± 1,116 viral reads” (LN 142) and “The total number of reads per sample was 196,489 ± 18,358” (LN 138).\nIf the total reads per sample are ~196,489, how is the virus read proportion only 0.4% ± 0.5%? Additionally, when I checked Supplementary Table 1, the average number of viral reads per sample is ~800 (Column AN), which does not match the reported 10,215 ± 1,116.\nCould the authors clarify whether there is an inconsistency in these values?\n# 5. Fig. 2 needs clearer descriptions. For example, what does the error bar mean in the Fig.2B? And the of Fig2B seems to be “mean of richness values in subsampling” as I understand, but the legends said “b, Slope of richness rarefaction curves across sampling months.”, which is overlaped with the Fig S2.\nOther legends could also be improved by providing more details. For example in Fig8A, what input data was used to calculate NMDS?\n# 6. Viral taxa should be written in italics throughout the manuscript.", "reviewer_2": "Review Nature Communications. NCOMMS-24-55589\nArctic Ocean virus communities: seasonality, bipolarity, and prokaryotic associations\nI would like to sincerely thank the authors for their thorough and thoughtful responses to not only my comments but also those of the other reviewers. Their efforts to address the feedback provided have significantly improved the quality of the manuscript, and I appreciate the time and attention they dedicated to refining their work.\nHowever, I would just like to highlight the following points:\nL 488-490. Please clarify the following statement: “Notably, in the entire dataset of ASV to vOTU correlations, only one had a predicted host of Nitrosopumilaceae, which was correlated with an unclassified Gammaproteobacterium (no family-level prediction)”.\nDoes this mean that one virus has two putative host assignments, one being an archaeon (Nitrosopumilaceae) and the other a Gammaproteobacterium?\nThe discussion about the finding that many of Fram Strait viruses are also abundant in deep samples could include also the possibility that ocean currents in specific areas play an important role. The reason why these viruses are found at such depths remains puzzling.\nAdditionally, the term 'epipelagic' may be more appropriate than 'surface waters' when referring to depths of less than 200 m.\nRegarding the new discussion paragraphs on how well these findings fit into host-virus evolutionary models, I completely agree with the authors that ecological models are often idealized in a chemostat-like setting. These dynamics primarily occur within seasons and should be evaluated through more spatially resolved temporal sampling.\nHowever, I would have like the authors to conduct a simple analysis: checking whether the exact same virus appears in consecutive summers, springs, or winters (not at 95%, if not at 100% ID). Are there any viral sequences with no genetic changes or SNPs observed across the temporal samples? The “Red Queen model” hypothesizes that 'fast' co-evolution occurs between hosts and viruses. If no genetic changes are detected, this hypothesis could be ruled out, at least within the analyzed period and possibly for a particular type of viruses. Instead, this would suggest the presence of a “viral seed bank” that persists in the environment between seasons but remains undetected due to the low abundance of its host at the change of season, influencing biodiversity in a way perhaps more aligned with the constant-diversity model. However, further analysis would be needed to confirm this hypothesis.\nI have no other questions for this manuscript.", "reviewer_3": "The authors have addressed all comments to a satisfactory degree. I have no further concerns or suggestions and recommend this manuscript for publication."}, "author_response": "Thank you for the help in improving our manuscript.\nI still have a few remaining suggestions to enhance robustness, transparency, and clarity of this paper.\n# 1. First of all, I have to say sorry that I made a mistake in my first revision about the total sequencing depth. I originally thought it was 10 Gb per sample, but, in fact, it is only 1 Gb (196,489 reads × 5,435 bp per read). This low sequencing depth makes the quantification of viruses more challenging than I expected, especially given that viral diversity does not reach saturation in the dataset.\nThe updated results indicate that viral diversity peaks in late summer to autumn, which is the somehow opposite of what was stated in the initial submission. This is a big shift in the overall viral seasonal dynamics.\nMy question is:\nThe error bars in Figure 2B are very large, making it difficult to compare diversity across months with confidence. Would it be possible to restrict the analysis to only samples with adequate viral reads when calculating diversity? Alternatively, it might be helpful to present a fully expanded diversity/richness plot for each year, similar to what was done for cVCR in Fig.1 and Fig.S1.\nWe appreciate the reviewers point here, we have now replaced the barchart with error bars with a fully expanded diversity/richness plot as was done for cVCR in Fig. S1, as requested. This shows that there is intra- and inter-annual variability in richness, which was contributing to the large error bars, and this version indeed provides for better context and interpretation.\nNotably – and this applies to all points related to the alpha diversity assessments – immediately below, as well as Point #4. We realized that when we computed the total viral read numbers we did so incorrectly (just a downstream coding misinterpretation of a mapping table). After re-calculating this value, which does not affect the cVCR or cVGB (calculations which were done independently of this explicitly calculated value), the total mapped reads are actually 8,825 ± 1,020 virus reads per sample. The main impact on the alpha diversity aspect was it allowed us to set higher sub-sampling (see below) which provided an overall more robust assessment of the alpha diversity, though of course still limited by relatively low-throughput Pac Bio relative to Illumina libraries.\nWe discuss more where relevant, below.\nIn addition to the alpha diversity assessments, it also impacted Supplementary Table (Tiara results), and a line in the text. See Point #4 below for more information about this.\nL119-123: On average, each sample harbored 9,815 ± 965 viral reads based on raw read counts (5.3 ± 0.5%), and 8,825 ± 1,020 reads mapping to vOTUs (6.8% ± 1.8%). Notably, despite not pre-filtering to remove eukaryotes, the large majority of reads in the dataset were prokaryotic (90.2% ± 1.4%), with eukaryotic reads making up a smaller proportion (7.2% ± 1.0%).\nAlso, I have a question regarding the subsampling approach. As described, subsampling was performed with “100 iterations of subsampling the vOTU count table at a range of different depths, spanning from 25 up to 16,000 counts at 50-count intervals.”\nDoes the 25% mapping breadth threshold still apply in this subsampling process? Specifically, if removing just one read causes the mapping coverage to drop below 25%, would that vOTU still be considered? The method here is important, particularly to richness estimation, where rare viruses may be affected. I suggest that the authors provide more details on how subsampling was performed with the mapping threshold.\nWe appreciate the reviewer’s question. Only contigs that had 25% breadth of coverage considering all the data available for a given sample were considered for the analysis. From that point, we did the sub-sampling (which was able to go to a higher value since we corrected the total sequence numbers, as described above).\nWe now provide more details about how subsampling was performed in the methods.\nL762-776: For this, we applied 100 iterations of subsampling the vOTU count table at a range of different depths, spanning from 50 up to 30,000 counts at 50 count intervals. For this step, all vOTUs that exceeded the 25% breadth of coverage cutoff when using all the data for a given sample were utilized. During each iteration, richness, evenness and Shannon diversity was calculated, with mean values being determined for each 50-count interval. The functions rrarefy, specnumber and diversity from the vegan package122 were used to perform subsampling and alpha diversity calculations. The mean values were used to generate a rarefaction curve of alpha diversity using the ggplot2 package123. To assess shifts in diversity over time, we compared the mean richness values obtained from subsampling at the 1000 count interval. In addition, we computed and compared the slopes of the sample rarefaction curves up until the subsampled depths, which represents the rate of vOTU discovery with increasing viral read counts. Comparing the slopes across samples enabled an assessment of whether the sample richness patterns observed at the subsampled depth would change if the viral read count was increased.\nIn addition, the usage of slope and direct subsampling points to compare diversity/richness across sequencing depths might be improved by using non-linear functions. For example, combining rarefaction and extrapolation using the iNEXT R-package could provide estimates of diversity/richness accounting for undetected viruses beyond observed data.\nWe appreciate the reviewers point here as well. Given the coverage of the long-reads, and that it is not yet approaching saturation, we had preferred not to extrapolate the curves. Also for this reason, we think that showing the accumulation curves is valuable. However, we did recalculate the linear modelled values, after reviewing the reviewers point we realized it would be better to calculate it only until 1000 reads (the subsampled depth used in Fig 2c. In addition to this, we applied iNEXT which shows results very similar to the trends we observed from the linear approach. Both approaches show lowest diversity estimates in the spring, and higher diversity in the summer. We clarify this mainly in Fig 2. Legend, along with adding the iNEXT results as Supplementary Fig 2. Also we update the methods to state how the iNEXT was employed.\nL1584-1594: Intra-and inter-annual viral richness variability. a, To assess differences in sequencing depth, we iteratively subsampled viral read counts from 50 up to 30,000 at 50 count intervals and at each interval determined the mean richness from 100 iterations. The mean richness across subsampled intervals was visualized in a rarefaction-style curve. b, Boxplots illustrating the slope of the vOTU richness rarefaction curves up until the chosen subsampling depth (1000 viral reads). The slopes represent the trajectory of vOTU discovery, and thus indicate how richness may appear if a higher sequencing depth was achieved. c, vOTU richness at subsampled depth across sampling months. Viral read counts were subsampled to 1000 per sample and vOTU richness calculating through 100 iterations. vOTU richness and standard deviation from the 100 iterations are illustrated through the bar and error bars, respectively.\niNext Methods, L776-778: To further explore this, we also employed the iNEXT v3.0.1 package125,126 to estimate the vOTU richness based on extrapolations to 30,000 viral read counts.\nGiven the big change in diversity results, I wonder whether another key metric, cVCR, and its associated results are also influenced by sequencing depth and normalization (subsampling). Since cVCR relies on the coverage of both viral and cellular genomes, low and varied sequencing depth could potentially bias cVCR values as well. So subsampling may potentially introduce some shifts in results, particularly in seasonal trends presented in Fig 1.\nI believe the observed cVCR-based abundance seasonal pattern is biologically meaningful, but I would prefer authors to reinforce the conclusions by addressing these potential biases.\nMy suggestion is to justify the use of cVCR by testing whether cVCR values remain stable when reads are subsampled to the same sequencing depth.\nWe appreciate the reviewer’s concern here, and to address it, we randomly sub-sampled all samples to 100,000 reads and re-performed the cVCR analysis. This excluded seven samples with fewer than 100,000 total reads. The results are consistent with those from the non-sub-sampled dataset in Fig. 1 and Fig. S1. There are some month-to-month differences in the cVCR vs the full dataset, which would be expected to stochastic variation and sometimes low variation in viruses in these sub-samples. Encouragingly, the overall magnitude of cVCR and the overall seasonal trends with this sub-sampling approach are highly similar to when using the entire dataset. We now provide a new panel to Fig. S1 which shows this results, and cite it in the results, adding also a detail about an alternative method we calculated, coverage per Gigabase pair.\nL158-161: These patterns were consistent also when we sub-sampled the total reads to 100,000 reads per sample, as a test of the robustness of this metric to sequencing depth, as well as when we calculated coverage per Gigabase (Supplementary Fig. 1); given these consistencies, we focus the remaining analyses on cVCR from the full dataset.\n# 2. As far as I understand, the parameter \"-B 1000\" used in the phylogenetic analysis, as written in the Methods section, refers to UFBoot. This should be explicitly stated in the Fig.5 legend, because different from traditional values, UFBoot values are only considered reliable if they are ≥ 95% (http://www.iqtree.org/doc/Frequently-Asked-Questions).\nWe have now labeled the branches in accordance with the reference that the reviewer suggested (SH-aLRT > 95% and Ultrafast Bootstraps >95), as shown in the figure.\nAdditionally, it would be greatly appreciated if the authors could provide the full psb phylogenetic tree and alignment as supplementary files to allow readers to evaluate evolutionary relationships.\nAlso we have uploaded the full sequences, alignment, trimmed alignments, and treefile to figshare and add it to the data availability statement: 10.6084/m9.figshare.28398173\nL1030-1032: psbA alignment, trimmed alignment, treefile, and labels for cyanophage references are available via figshare (DOI 10.6084/m9.figshare.28398173).\nWe also corrected the heatmap legend -- in editing the figure for the previous re-submission we had mistakenly made the values 100x too high. It is now consistent with the original submission, and the text, which was correct in each version.\n#3. The paper [https://pubmed.ncbi.nlm.nih.gov/39448562/] provides a background on polar viruses in a time-series context. The authors may consider referencing it to strengthen the introduction and discussion.\nWe add to the Introduction: L58-60: In addition to spatial structuring, polar viral communities also shift over time, for example in the Antarctic, following bloom dynamics and assembling into distinct communities across seasons23.\nWe add to the Discussion: L513-515: In terms of total virus prevalence, we found strong increases in the ratio of viruses to prokaryotes in the late summer months, mirroring findings from the Antarctic where viruses peaked in mid to late summer23.\n#4. I need some clarification for the newly added values in LN 145 because “On average, each sample harbored 10,215 ± 1,116 viral reads” (LN 142) and “The total number of reads per sample was 196,489 ± 18,358” (LN 138).\nIf the total reads per sample are ~196,489, how is the virus read proportion only 0.4% ± 0.5%? Additionally, when I checked Supplementary Table 1, the average number of viral reads per sample is ~800 (Column AN), which does not match the reported 10,215 ± 1,116.\nCould the authors clarify whether there is an inconsistency in these values?\nAs noted above for point #1, we appreciate the reviewers point here which made us realize that indeed there was an inconsistency in the values between the text and Supplementary Data 1. The virus read values in Supplementary Data 1 were incorrect (just a misinterpretation based on a downstream column of mapping data). We have re-calculated the aspects of the paper in which this was relevant, namely the line in the text referenced by the reviewer, Supplementary Data 1, and the Tiara results (in which we remove viral reads before calculating, since Tiara has no explicit viral prediction), and viral richness maximum rarefaction values. Fortunately, no other aspects were impacted by this inconsistency.\nL119-123: On average, each sample harbored 9,815 ± 965 viral reads based on raw read counts (5.3 ± 0.5%), and 8,825 ± 1,020 reads mapping to vOTUs (6.8% ± 1.8%). Notably, despite not pre-filtering to remove eukaryotes, the large majority of reads in the dataset were prokaryotic (90.2% ± 1.4%), with eukaryotic reads making up a smaller proportion (7.2% ± 1.0%).\n#5. Fig. 2 needs clearer descriptions. For example, what does the error bar mean in the Fig.2B? And the of Fig2B seems to be “mean of richness values in subsampling” as I understand, but the legends said “b, Slope of richness rarefaction curves across sampling months.”, which is overlaped with the Fig S2.\nAs previously mentioned above, we have revised the Figure 2 legend for clarity and added another panel to Fig. 2. The updated Fig. 2 now includes (a) observed richness of vOTUs shown as rarefaction curves, (b) summary of the slopes of the rarefaction curves as estimates of overall richness and (c) vOTU richness per sampling time point presented as bar plots at a depth of 1000 reads per sample. We believe that the revised figure is more clear and robust over the previous versions, and thank the reviewer for their help in improving it.\nFull Fig. 2 legend (Line 1584-1594) is listed above as part of point #1\nOther legends could also be improved by providing more details. For example in Fig8A, what input data was used to calculate NMDS?\nWe now update Fig. 8A stating what input data was used to calculate NMDS, as well as providing some clarification in the results section related to 8a and other related figures, which we believe clarify the interpretation.\nL1650-1662: Fig. 8. Viral amino acid traits across environmental gradients. a, Nonmetric multidimensional scaling (NMDS) plot showing the environmental diversity of the sampled Fram Strait and GOV2.0 ecosystems. The Bray-Curtis distance similarity matrix was calculated based on the available environmental parameters of 38 Fram Strait and 69 GOV2.0 samples and used to generate NMDS coordinates of each sample. Point shapes represent dataset origin, and colors distinguish different temperature ranges. Vectors show correlations with environmental variables. This ordination provides context for the correlations to environmental parameters of amino acid traits shown in part b. b, Heatmap plot illustrating the Spearman correlation coefficients of environmental parameters to amino acid traits, with p-values represented by asterisks as indicated. The Spearman coefficients were calculated using a Mantel test using pairwise distances of each environmental parameter (Euclidean distance) vs. each amino acid trait (Bray-Curtis distance).\nWe also amended the results section to clarify that the NMDS is based only on environmental data which is meant to provide the environmental overview/context and a sense of the diversity of the samples surveyed, in the correlations shown in the Figure 8b. Also, we more directly cited the complementary Figures in this section, which breakdown these correlations by family in the supplement.\nL460-476: To ensure fair comparisons between the datasets, we focused only on GOV samples originating from ≤35 m water depth (i.e. the average depth of Fram Strait sampling) and first analysed the environmental context, which shows strong clustering of GOV2.0 polar samples to the Fram Strait, based on their similarity in distance to equator, oxygen, and temperature (Fig. 8a). Linking those environmental data to the examined protein features for all viruses, revealed that the aliphatic index and the nitrogen usage score were positively correlated with temperature (Fig. 8b), whereas other traits of potential cold adaptation, in particular polar charged and uncharged amino acids, showed significant correlation with one or more of the other environmental parameters, but not directly with temperature. This was also the case for Caudoviricetes, while Megaviricetes amino acid traits were more generally positively correlated with oxygen, temperature, and chlorophyll, and ‘other’ vOTUs were generally positively correlated with salinity and oxygen (Supplementary Fig. 9). In the latter cases, these groups are more poorly sampled due to their low abundance resulting in weaker correlations overall, and especially necessitating further study. Examining the average values for these protein parameters according to sample and temperature, across all viruses, reflected the patterns observed via other statistical analyses (Supplementary Fig. 10), reinforcing our observations.\n# 6. Viral taxa should be written in italics throughout the manuscript.\nAll mentions of viral taxa have been italicized.\nWe have now also italicized prokaryotic Families, Orders, etc (in accordance to recent Nat. Comm. papers and Parker CT, Tindall BJ, Garrity GM. International code of nomenclature of prokaryotes. Prokaryotic code (2008 Revision) International Journal of Systematic and Evolutionary Microbiology. 2019;69(1A):S1–S111. doi: 10.1099/ijsem.0.000778.\nReviewer #2 (Remarks to the Author):\nReview Nature Communications. NCOMMS-24-55589\nArctic Ocean virus communities: seasonality, bipolarity, and prokaryotic associations\nI would like to sincerely thank the authors for their thorough and thoughtful responses to not only my comments but also those of the other reviewers. Their efforts to address the feedback provided have significantly improved the quality of the manuscript, and I appreciate the time and attention they dedicated to refining their work.\nThank you for your help in improving our manuscript.\nHowever, I would just like to highlight the following points:\nL 488-490. Please clarify the following statement: “Notably, in the entire dataset of ASV to vOTU correlations, only one had a predicted host of Nitrosopumilaceae, which was correlated with an unclassified Gammaproteobacterium (no family-level prediction)”.\nDoes this mean that one virus has two putative host assignments, one being an archaeon (Nitrosopumilaceae) and the other a Gammaproteobacterium?\nWe have tried to clarify this:\nL390-395: Notably, among vOTUs that had correlations to hosts in the correlation network, only one had a predicted host (based on iPHoP) of Nitrosopumilaceae. This vOTU was correlated with an unclassified Gammaproteobacterium (no family-level prediction) within the correlation network, highlighting again the difficulty comparing directly these two complex and independent approaches.\nThe discussion about the finding that many of Fram Strait viruses are also abundant in deep samples could include also the possibility that ocean currents in specific areas play an important role. The reason why these viruses are found at such depths remains puzzling.\nWe changed one line of the results to be more precise, i.e., from ‘detected commonly’ vs ‘common’ which may imply more strongly that they are abundant, vs just observed.\nWe also add to the Discussion:\nL584-598: Additionally, in general, the detection of the Fram Strait viruses in the deep ocean is notable, and is potentially related to export from the surface and deep ocean currents, some of which originate at the Arctic. Further inspection into the specific vOTUs and from other datasets at the global level could help further evaluate these ideas. Part of the explanation may be that it is vOTUs we primarily observed in the deep waters are from M6, the ‘persistently present’ Fram Strait module, thus, they may be associated also with generally persistent/ubiquitous host lineages and/or have broad host ranges, and furthermore the relative lack of seasonality in the deep sea may help with their consistent detection. In the discussion, we add:\nL581-598: Further investigations to unravel the diversity and dynamics of viruses in Arctic and Antarctic are crucial, given that climate shifts are causing major biological and physicochemical perturbations in these regions. Additionally, in general, the detection of the Fram Strait viruses in the deep ocean is puzzling, and is potentially related to export from the surface and deep ocean currents, some of which originate at the Arctic. Further inspection into the specific vOTUs and from other datasets at the global level could help further evaluate these ideas. Part of the explanation may be that it is vOTUs we primarily observed in the deep waters are from M6, the ‘persistently present’ Fram Strait module, thus, they may be associated also with generally persistent/ubiquitous host lineages and/or have broad host ranges, and furthermore the relative lack of seasonality in the deep sea may help with their consistent detection.\nAdditionally, the term 'epipelagic' may be more appropriate than 'surface waters' when referring to depths of less than 200 m.\nWe have changed surface waters to epipelagic throughout the manuscript.\nRegarding the new discussion paragraphs on how well these findings fit into host-virus evolutionary models, I completely agree with the authors that ecological models are often idealized in a chemostat-like setting. These dynamics primarily occur within seasons and should be evaluated through more spatially resolved temporal sampling.\nHowever, I would have like the authors to conduct a simple analysis: checking whether the exact same virus appears in consecutive summers, springs, or winters (not at 95%, if not at 100% ID). Are there any viral sequences with no genetic changes or SNPs observed across the temporal samples? The “Red Queen model” hypothesizes that 'fast' co-evolution occurs between hosts and viruses. If no genetic changes are detected, this hypothesis could be ruled out, at least within the analyzed period and possibly for a particular type of viruses. Instead, this would suggest the presence of a “viral seed bank” that persists in the environment between seasons but remains undetected due to the low abundance of its host at the change of season, influencing biodiversity in a way perhaps more aligned with the constant-diversity model. However, further analysis would be needed to confirm this hypothesis.\nWe thank the reviewer’s thoughtful consideration of the discussion we wrote in response to their previous suggestions, and are glad that they appreciate it.\nIn regards to the new analyses proposed, while a comprehensive analysis in this direction would be somewhat outside the scope of our already quite long paper, we have considered and addressed it in two ways. First, we explicitly compared viral contigs >10,000 bp across years. Out of 7,360, we found 21 clusters of contigs, encompassing 55 contigs total, that were 100% matches (akin to 100% vOTUs) between at least two years. This suggests that detection of viruses with no genetic changes across years is rare, though of course there is a limitation of sequencing depth. Second, we employed a state-of-the-art tool, inStrain (https://www.nature.com/articles/s41587-020-00797-0), which compares reads mapped to vOTUs, and examines for SNPs. We used 98% minimum ANI, the most stringent threshold that the developers recommend, and a minimum coverage of 5 (default setting) to call a variant. Then we filtered for contigs where more than 25% of their bases were covered by at least one read, i.e. breadth > 0.25. To confirm a SNV, we required a minimum SNV frequency of 0.05 and a false discovery rate of 1e-06 (default settings). Among the 202 comparisons (i.e., enabled by coverage exceeding minima required to detect SNPs), we find that there tends to be an underlying diversity within the vOTU populations as reflected by inStrain’s consensus SNP metric, though the amount varies between values that are likely insignificant (in accordance with the documentation (>0.9999) to about 0.97. This highlights that, in general, vOTUs are made up of populations of highly similar viruses, rather than being clonal. This has also been observed elsewhere, including the arctic.\nTherefore, we cannot rule out the “Red Queen Model”, and we can also not rule out the constant-diversity model without more detailed and comprehensive temporal and spatial data.\nAs this analysis is outside the original scope of our paper and is thus not a comprehensive interrogation of the strain diversity which would indeed be a valuable future direction, we add two sentences in the alpha diversity section of the manuscript related to the new analysis which also enable us to add a point to the discussion. Also we add the inStrain results as a supplementary figure.\nResults, L211-222: To further contextualize the dynamics we observed at the vOTU level, we examined the underlying sequence diversity based on examining exact contigs across years and single nucleotide variation. First, we explicitly compared viral contigs >10,000 bp across years. Out of 7,360 contigs, we found 21 clusters of contigs, encompassing 55 contigs total, that were 100% matches between at least two years. This suggests that detection of viruses with no genetic changes across years is rare, though there is a limitation of sequencing depth. Second, based on analysis of SNPs in reads to vOTUs,we find that there tends to be an underlying diversity within the vOTU populations, though the amount varies between vOTUs ranging from likely insignificant variation related to sequencing error, to high differentiation (e.g., 97% similarity) (Supplementary Fig. #). These analyses highlight that, in general, vOTUs are made up of populations of highly similar viruses, rather than being clonal. This has also been observed elsewhere, including the Arctic14.\nDiscussion, L623-629: Our analysis of microdiversity of the Arctic virus communities demonstrated both instances of identical viruses across years, as well as a general tendency for vOTU to be underlaid by sequences with very high similarity but not exact matches to the vOTUs, akin to that seen elsewhere in the ocean29,30,88, including the Arctic Ocean14. Thus, based on these results we cannot rule out either the Red-Queen hypothesis or Constant-diversity.\nI have no other questions for this manuscript.\nReviewer #2 (Remarks on code availability):\nThe code is well-organized and easy to follow, which makes it a valuable resource for the community.\nThank you.\nReviewer #3 (Remarks to the Author):\nThe authors have addressed all comments to a satisfactory degree. I have no further concerns or suggestions and recommend this manuscript for publication.\nThank you."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I appreciate the authors for the thoughtful replies and the updates. I have no further concerns.", "reviewer_2": "Thank you for the detailed and thoughtful responses to all my comments and suggestions, as well as for the additional analyses and clarifications provided throughout the revised manuscript. I appreciate the effort made to address both methodological and conceptual aspects, including the clarification of host-virus associations, the discussion on deep-sea viral persistence, and the additional inStrain analysis on viral microdiversity. I have no further concerns and fully support the acceptance of this version of the manuscript. Just for the record, I would like to highlight my appreciation for the clarity of the revised discussion and the transparency of the supplementary data. Thank you again."}, "author_response": ""}]

Decision letter