Simon Roux

Jennifer R Brum

Bas E Dutilh

Shinichi Sunagawa

Melissa B Duhaime

Alexander Loy

Bonnie T Poulos

Natalie Solonenko

Elena Lara

Julie Poulain

Stéphane Pesant

Stefanie Kandels-Lewis

Céline Dimier

Marc Picheral

Sarah Searson

Corinne Cruaud

Adriana Alberti

Carlos M Duarte

Josep M Gasol

Dolors Vaqué

Tara Oceans

Silvia G Acinas

Patrick Wincker

Matthew B Sullivan

Ecogenomics and potential biogeochemical impacts of globally abundant ocean viruses

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

definition explanation). This expands ocean viral populations nearly 3-fold over the prior TOV dataset [START_REF] Brum | Patterns and ecological drivers of ocean viral communities[END_REF] , while also improving average contig lengths and genomic context 2.5-fold for TOV-known populations (Supplementary Table 2). Rarefaction analyses show that while mesopelagic viral communities remain undersampled, epipelagic viral communities now appear near-completely sampled (Extended Data Fig. 1A). Because bathypelagic communities were underrepresented due to cellular contamination, we focused the remaining analyses on 15,222 non-bathypelagic viral populations.

We first categorized viral populations into viral clusters, or VCs using shared gene content information and network analytics [START_REF] Lima-Mendez | Reticulate representation of evolutionary and functional relationships between phage genomes[END_REF] (see Supplementary Fig. 1 for VC definition schematic). This method starts from genome fragments (≥10kb) and results in VCs approximately equivalent to known viral genera [START_REF] Lima-Mendez | Reticulate representation of evolutionary and functional relationships between phage genomes[END_REF][START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] . Clustering of the 15,222 GOV viral populations with 15,929 publicly available bacterial and archaeal viruses revealed 1,259 VCs (see Supplementary Table 3, Supplementary Text & Extended Data Fig. 2 for comparison with alternative classification methods). Of these, 658 included exclusively GOV sequences, which approximately doubles known bacterial and archaeal virus genera [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] , and another 209 VCs contained at least one GOV sequence (Fig. 1B). As with viral populations, rarefaction analyses suggested that VC diversity was undersampled in mesopelagic waters, but near-completely sampled in epipelagic waters (Extended Data Fig. 1B).

We next identified the most abundant and widespread VCs based on read recruitment of VC members. In each sample, a fraction of the VCs were identified as abundant based on their cumulative contribution to sample diversity (estimated with Simpson Index, abundant VCs represent 80% of the total sample diversity, Extended Data Fig. 1C). By these criteria, only 38 of 867 observed VCs were abundant in two or more stations, and together recruited an average of 50% and 35% of reads from viral populations for epipelagic and mesopelagic samples, respectively (Supplementary Table 3). Four of these 38 abundant VCs were also relatively ubiquitous as they were abundant in more than 25 stations, and 62 of the 91 non-bathypelagic samples were dominated by 1 of these 4 VCs (Fig. 2 A &B). Among the 38 abundant VCs, only 2 corresponded to well-studied viruses, from the T4 superfamily [START_REF] Sullivan | Genomic analysis of oceanic cyanobacterial myoviruses compared with T4-like myoviruses from diverse hosts and environments[END_REF][START_REF] Zhao | Abundant SAR11 viruses in the ocean[END_REF] (VC_2, 1 of the 4 ubiquitous) and the T7virus genus [START_REF] Labrie | Genomes of marine cyanopodoviruses reveal multiple origins of diversity[END_REF] (VC_9). Eight represented known but unclassified viral isolates, 10 included viruses known only from environmental sequencing [START_REF] Mizuno | Expanding the marine virosphere using metagenomics[END_REF][START_REF] Chow | Combining genomic sequencing methods to explore viral diversity and reveal potential virus-host interactions[END_REF] , and the remaining 18 VCs were completely novel (Fig. 2C, Extended Data Fig. 3).

Given this global map of the dominant dsDNA viral types in the oceans, we next sought to identify the range of hosts these viruses infect. This is challenging, as culture-based methods insufficiently capture naturally-occurring diversity, whereas metagenomic approaches broadly survey viral diversity but often without host information. Fortunately, sequence-based approaches are emerging that examine similarities between (i) viral genomes and host CRISPR spacers [START_REF] Andersson | Virus population dynamics and acquired virus resistance in natural microbial communities[END_REF] , (ii) viral and microbial genomes due to integrated prophages or gene transfers [START_REF] Mizuno | Expanding the marine virosphere using metagenomics[END_REF] , and (iii) viral and host genome nucleotide signatures (here, tetranucleotide frequencies [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] , see Supplementary Table 4 and Supplementary Text for discussion of the accuracy/sensitivity of in silico host prediction methods). We applied all 3 methods to GOV to predict hosts at the phylum level, or class level for Proteobacteria (Supplementary Table 5), then summarized these results at the VC level. This led to host range predictions for 392 of 867 VCs -all with confidence assessed by comparison to a null model (Supplementary Fig. 2 and Supplementary Table 3).

The hosts of the 38 globally abundant VCs were largely restricted to abundant and widespread epipelagic-ocean microbes that were previously identified via mi Tag-based OTU counts in Tara Oceans microbial metagenomes [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] . Notably, the 4 ubiquitous and abundant VCs were predicted to infect 7 of the 8 globally abundant microbial groups (Actinobacteria, Alpha-, Delta-, and Gammaproteobacteria, Bacteroidetes, Cyanobacteria, Deferribacteres; Fig. 2C, Extended Data Fig. 4). The 8 th abundant microbial group, Euryarchaeota, was not linked to these 4 VCs, but was predicted as a host for 3 of the 34 other abundant VCs (VC_3, VC_27, and VC_63, Extended Data Fig. 3). Among the 38 abundant VCs, the number of VCs predicted to infect a given microbial host phylum (or class for Proteobacteria) was positively correlated with host global richness rather than relative abundance (Extended Data Fig. 4B). This suggests that, likely because ocean viruses appear globally distributed [START_REF] Brum | Patterns and ecological drivers of ocean viral communities[END_REF] , widespread and abundant hosts that are minimally diverse (e.g. Cyanobacteria) provide few viral niches, whereas more diverse host groups, even at lower abundance (e.g. Betaproteobacteria), provide more opportunity for viral niche differentiation. Hence, these host associations provide critically-needed empirical support for hypotheses derived from global virus-host network models [START_REF] Flores | Multi-scale structure and geographic drivers of crossinfection within marine bacteria and phages[END_REF] .

Having mapped viral diversity and predicted virus-host pairings, we next sought to identify virusencoded AMGs that might modify host metabolism during infection and likely impact biogeochemistry. To maximize AMG detection, all 298,383 viral contigs >1.5kb were examined, including small contigs not associated with a viral population. This revealed 243 putative AMGs (Supplementary Table 6). While 95 of these AMGs were known (reviewed in ref. [START_REF] Hurwitz | Depth-stratified functional and taxonomic niche specialization in the "core" and "flexible" Pacific Ocean Virome[END_REF] ), others offer insights into how viruses may directly manipulate microbial metabolisms. Here we focus on 4 (dsrC, soxYZ, P-II and amoC; see Extended Data Table 1, Supplementary Figs. 3456and Supplementary Text for functional affiliation of these AMGs) because of their putative roles in sulfur or nitrogen cycling. Three of these are not known in viruses, and one, dsrC, has only been observed in viruses from anoxic deep-sea environments [START_REF] Roux | Ecology and evolution of viruses infecting uncultivated SUP05 bacteria as revealed by single-cell-and meta-genomics[END_REF][START_REF] Anantharaman | Sulfur Oxidation Genes in Diverse Deep-Sea Viruses[END_REF] .

Sulfur oxidation in seawater involves two central microbial pathways -dissimilatory sulfur reductase (Dsr) and sulfur oxidation (Sox) [START_REF] Friedrich | Prokaryotic sulfur oxidation[END_REF] -and GOV AMG analyses revealed that epipelagic viruses encode key genes for each. First, 11 dsrC-like genes were identified in viral contigs (Extended Data Fig. 5). The Dsr operon is used by sulfate/sulfite-reducing microbes in anoxic environments, as well as sulfur-oxidizing bacteria in oxic and anoxic environments (Fig. 3A) [START_REF] Friedrich | Prokaryotic sulfur oxidation[END_REF] . DsrC, specifically, provides sulfur to DsrAB-sulfite reductase for processing through a conserved C-terminal motif (Cys B X 10 Cys A ), and dictates sulfur metabolism rates [START_REF] Santos | A protein trisulfide couples dissimilatory sulfate reduction to energy conservation[END_REF] . Other DsrC-like proteins (also known as TusE) lack Cys B and instead participate to tRNA modification [START_REF] Venceslau | The "bacterial heterodisulfide" DsrC is a key protein in dissimilatory sulfur metabolism[END_REF] . In GOV, four clades of DsrC-like sequences were similar to TusE (DsrC-1 to DsrC-4), whereas the fifth (DsrC-5) was similar to bona fide DsrC (Extended Data Fig. 5, Extended Data Table 1, Supplementary Fig. 3, and Supplementary text). Second, 4 soxYZ genes were identified on viral contigs (Extended Data Fig. 6). Like DsrC, SoxYZ is an important sulfur carrier harboring a conserved functional motif identified in all GOV SoxYZ proteins (Fig. 3A, Supplementary Fig. 4, and Supplementary text) [START_REF] Dahl | Sulfite oxidation in the purple sulfur bacterium Allochromatium vinosum: identification of SoeABC as a major player and relevance of SoxYZ in the process[END_REF] .

Other AMGs suggest marine viruses may manipulate nitrogen cycling. First, 10 GOV contigs encoded P-II, a gene widespread across bacteria and archaea and central in nitrogen metabolism regulation (Fig. 3B) [START_REF] Huergo | P(II) signal transduction proteins: nitrogen regulation and beyond[END_REF] . Three AMG clades (P-II-1, P-II-2, and P-II-4) displayed both P-II conserved motifs and had predicted structures similar to bona fide P-II, whereas the fourth clade (P-II-3) is functionally ambiguous as it lacked a conserved motif (Supplementary Fig. 5, and Supplementary text). Second, two P-II AMG clades (P-II-1 and P-II-4) were proximal to an ammonium transporter gene, amt, in GOV contigs (Extended Data Fig. 7). In bacteria, such an arrangement is a signature of P-II-like genes that specifically activate alternative nitrogen production and ammonia uptake pathways during nitrogen starvation [START_REF] Huergo | P(II) signal transduction proteins: nitrogen regulation and beyond[END_REF] . Third, one GOV contig included amoC, encoding the subunit C of ammonia monooxygenase, suggesting a role in ammonia oxidation [START_REF] Stahl | Physiology and diversity of ammonia-oxidizing archaea[END_REF] . While functional annotation is challenging for these genes [START_REF] Stahl | Physiology and diversity of ammonia-oxidizing archaea[END_REF] , and functional motifs are not yet known, the translated AMG was 94% identical to functional AmoC from Thaumarchaeota -a level of identity only observed among expressed and functional AMGs (Extended Data Fig. 8, Supplementary Fig. 6, and Supplementary text).

Next, we investigated the origin, evolutionary history, and diversity of these AMGs in epipelagic viruses (see Supplementary Text for additional discussion about taxonomic affiliation and host prediction for AMG-containing GOV sequences). The 15 GOV contigs encoding dsrC or soxYZ genes, when affiliated, were all associated with members of the abundant and ubiquitous VC_2 (T4 superfamily, Extended Data Fig. 5 and6, Extended Data Table 1). Phylogenies suggested that these viruses obtained AMGs from S-oxidizing proteobacterial hosts, with likely a single transfer event for soxYZ and two for dsrC (Extended Data Fig. 5 and6). Among the latter, the bona fide S-oxidation DsrC-5 was most closely related to a clade of uncultivated S-oxidizing Gammaproteobacteria (MED13k09, Supplementary Fig. 7). These bacteria are widespread in the epipelagic ocean [START_REF] Loy | Reverse dissimilatory sulfite reductase as phylogenetic marker for a subgroup of sulfur-oxidizing prokaryotes[END_REF] and suspected to degrade dimethyl sulfide, a key reduced sulfur species involved in ocean-to-atmosphere sulfur transport and cloud formation. If confirmed, DsrC5-encoding viruses infecting these bacteria would impact critical sulfur cycling steps throughout surface waters. In contrast to sulfur AMGs, phylogenies suggest that P-II AMGs originated from diverse viruses (6 VCs including the abundant VC_2 and VC_12), and were acquired at least 4 times independently from Bacteroidetes, Proteobacteria, and possibly Verrucomicrobia (Extended Data Fig. 7, and Supplementary Text). Finally, while a single amoC AMG offers only preliminary evaluation of its evolutionary history, this amoCencoding contig appears to represent novel and rare archaeal dsDNA viruses (VC_623), predicted to infect ammonia-oxidizing Thaumarchaeota, known for their major role in global nitrification [START_REF] Pester | The Thaumarchaeota: an emerging view of their phylogeny and ecophysiology[END_REF] (Extended Data Fig. 8).

Finally, we investigated the ecology of viruses encoding these AMGs by mapping their distribution across GOV. Seven AMG clades were geographically restricted (DsrC-unc, DsrC-1, DsrC-2, DsrC-4, P-II-2, P-II-3, and amoC), and 5 were widespread throughout epipelagic (DsrC-3, DsrC-5, SoxYZ, P-II-1) or mesopelagic (P-II-4) waters (Fig. 3C). All widespread epipelagic AMGs were detected in waters of mid-range temperatures. In contrast, DsrC-5 and SoxYZ were predominantly detected in lownutrient conditions, while P-II-1 was predominantly detected in high-nutrient conditions (Fig. 3D, Extended Data Fig. 9). Thus, we hypothesize that viruses utilize DsrC-5 or SoxYZ to boost sulfur oxidation rates when infecting sulfur oxidizers in low-nutrient conditions, and P-II under high-nutrient conditions. The latter could be useful to viruses by activating expensive alternative N-producing pathways typically used only under N-starvation conditions [START_REF] Huergo | P(II) signal transduction proteins: nitrogen regulation and beyond[END_REF] . Consistent with this, metatranscriptomes from three low-nutrient stations (11_SRF in Mediterranean Sea, 39_DCM in Arabian Sea, and 151_SRF in Atlantic Ocean) revealed expression of viral homologs of dsrC and soxYZ but not of P-II (Extended Data Table 1).

Overall, this systematically collected and processed GOV dataset provides a critical resource for marine microbiology. This map of global dsDNA ocean viral diversity, at both the population and VC level, and viral-encoded AMGs brings global ecological context to abundant surface and deep ocean viruses. Both will also help interpret future (meta)genomic datasets and select experimental systems to develop. Together with recent experimental, informatic and theoretical advances [START_REF] Brum | Rising to the challenge: accelerated pace of discovery transforms marine virology[END_REF][START_REF] Dutilh | A highly abundant bacteriophage discovered in the unknown sequences of human faecal metagenomes[END_REF][START_REF] Weitz | A multitrophic model to quantify the effects of marine viruses on microbial food webs and ecosystem processes[END_REF] , this fundamental resource will accelerate the field towards understanding and dynamically predicting the roles and planetary impacts of viruses in nature.

Methods

Sample collection and processing Tara Oceans expedition

Ninety samples were collected between October 10, 2009, and December 12, 2011, at 45 locations throughout the world's oceans (Supplementary Table 1) through the Tara Oceans Expedition 32 . These included samples from a range of depths: surface, deep chlorophyll maximum, bottom of mixed layer when no deep chlorophyll maximum was observed (Station 123, 124, and 125), and mesopelagic samples. The sampling stations were located in 7 oceans and seas, 4 different biomes and 14 Longhurst oceanographic provinces (Supplementary Table 1). For TARA station 100, two different peaks of chlorophyll were observed, so two samples were taken at the shallow (100_DCM) and deep (100_dDCM) chlorophyll maximum. For each sample, 20 L of seawater were 0.22 µm-filtered and viruses were concentrated from the filtrate using iron chloride flocculation [START_REF] John | A simple and efficient method for concentration of ocean viruses by chemical flocculation[END_REF] followed by storage at 4ºC. After resuspension in ascorbic-EDTA buffer (0.1 M EDTA, 0.2 M Mg, 0.2 M ascorbic acid, pH 6.0), viral particles were concentrated using Amicon Ultra 100 kDa centrifugal devices (Millipore), treated with DNase I (100U/mL) followed by the addition of 0.1M EDTA and 0.1M EGTA to halt enzyme activity, and extracted as previously described [START_REF] Hurwitz | Evaluation of methods to concentrate and purify ocean virus communities through comparative, replicated metagenomics[END_REF] . Briefly, viral particle suspensions were treated with Wizard PCR Preps DNA Purification Resin (Promega, WI, USA) at a ratio of 0.5 mL sample to 1 mL resin, and eluted with TE buffer (10 mM Tris, pH 7.5, 1 mM EDTA) using Wizard Minicolumns. Extracted DNA was Covaris-sheared and size selected to 160-180 bp, followed by amplification and ligation per the standard Illumina protocol. Sequencing was done on a HiSeq 2000 system (101 bp, paired end reads) at the Genoscope facilities (Paris, France).

Temperature, salinity, and oxygen data were collected from each station using a CTD (Sea-Bird Electronics, Bellevue, WA, USA; SBE 911plus with Searam recorder) and dissolved oxygen sensor (Sea-Bird Electronics; SBE 43). Nutrient concentrations were determined using segmented flow analysis 35 and included nitrite, phosphate, nitrite plus nitrate, and silica. Nutrient concentrations below the detection limit (0.02 µmol kg -1 ) are reported as 0.02 µmol kg -1 . All data from the Tara Oceans expedition are available from ENA (for nucleotide) and from PANGAEA (for environmental, biogeochemical, taxonomic and morphological data) [START_REF] Oceans | Registry of all samples from the Tara Oceans Expedition[END_REF][START_REF] Oceans | Environmental context of all samples from the Tara Oceans Expedition[END_REF][START_REF]Biodiversity context of all samples from the Tara Oceans Expedition[END_REF] .

Malaspina expedition

Thirteen bathypelagic samples and one mesopelagic sample were collected between April 19, 2011 and July 11, 2011 during the Malaspina 2010 global circumnavigation covering the Pacific and the North Atlantic Ocean. All samples were taken at 4,000 m depth except two samples from stations 81 and 82 collected at 3,500 and 2,150 m respectively (Supplementary Table 1). Additionally, Station M114 was sampled at the OMZ region at 294 m depth. For each sample, 80 L of seawater were 0.22 µm-filtered and viruses were concentrated from the filtrate using iron chloride flocculation [START_REF] John | A simple and efficient method for concentration of ocean viruses by chemical flocculation[END_REF] followed by storage at 4°C. More details about the sampling and additional variables used in the Malaspina expedition can be found in ref. [START_REF] Salazar | Global diversity and biogeography of deep-sea pelagic prokaryotes[END_REF] . Further processing was done as for the Tara Oceans samples, except that Illumina sequencing was done at DOE JGI Institute (151 bp, paired end reads).

Dataset generation Contigs assembly

An overview of the contigs generation process is provided in Supplementary Fig. 8. The first step consisted in the generation of a set of contigs using as many reads as possible from the 104 oceanic viromes, including 74 epipelagic and 16 mesopelagic samples from the Tara Oceans expedition [START_REF] Karsenti | A holistic approach to marine eco-systems biology[END_REF] , and 1 mesopelagic and 13 bathypelagic from the Malaspina expedition [START_REF] Duarte | Seafaring in the 21st Century : The Malaspina 2010 Circumnavigation Expedition[END_REF] . This set of contigs was generated through an iterative cross-assembly [START_REF] Dutilh | A highly abundant bacteriophage discovered in the unknown sequences of human faecal metagenomes[END_REF] (using MOCAT [START_REF] Kultima | MOCAT: a metagenomics assembly and gene prediction toolkit[END_REF] and Idba_ud 41 , Supplementary Fig. 8) as follows: (i) high-quality (HQ) reads were first assembled sample by sample with the MOCAT pipeline as described in [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] , (ii) all reads not mapping (Bowtie 2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] , options --sensitive, -X 2000, and --nondeterministic, other parameters at default) to a MOCAT contig (by which we denote 'scaftigs', that is, contigs that were extended and linked using the paired-end information of sequencing read [START_REF] Luo | SOAPdenovo2 : an empirically improved memory-efficient short-read de novo assembler[END_REF] ) were assembled sample by sample with Idba_ud (iterative k-mer assembly, with k-mer increasing from 20 to 100 by step of 20), (iii) all reads remaining unmapped to any contig were then pooled by Longhurst province (i.e. unmapped reads from samples corresponding to the same Longhurst province were gathered) and assembled with Idba_ud (with the same parameters as above), and (iv) all remaining reads unmapped from every samples were gathered for a final cross-assembly (using Idba_ud). This resulted in 10,845,515 contigs (Supplementary Fig. 8B).

Genome binning and re-assembly

The contigs assembled from the marine viral metagenomes might still contain redundant sequences derived from the same, or closely related populations. We set out to merge contigs derived from the same population into clusters representing population genomes. To this end, contig sequences were first clustered at 95% global average nucleotide identity (ANI) with cd-hit-est [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF] (options -c 0.95 -G 1 -n 10mask NX, Supplementary Fig. 8B), resulting in 10,578,271 non-redundant genome fragments. Next, we used co-abundance (i.e. correlation between abundance profiles estimated by reads mapping) and nucleotide usage profiles of the non-redundant contigs to further identify contigs derived from the same populations with Metabat [START_REF] Kang | an efficient tool for accurately reconstructing single genomes from complex microbial communities[END_REF] . Briefly, Metabat uses Pearson correlation between coverage profiles (determined from the mapping of HQ reads of each sample to the contigs with Bowtie 2 [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] , options -sensitive, -X 2000, and --non-deterministic, other parameters at default) and tetranucleotide frequencies to identify contigs originating from the same genome (Metabat parameters: 98% minimum correlation, mode "sensitive", see Supplementary Text for more detail about the selection of these parameters). The 8,744 bins generated, including 3,376,683 contigs, were further analyzed, alongside 623,665 contigs not included in any genome bin but ≥1.5kb.

In an attempt to better assemble these genome bins, two additional sets of contigs were generated for each genome bin (beyond the set of initial contigs binned by Metabat [START_REF] Kang | an efficient tool for accurately reconstructing single genomes from complex microbial communities[END_REF] ), based on the de novo assembly of (i) all reads mapping to the contigs in the genome bin, and (ii) only reads from the sample displaying the highest coverage for the genome bin (both assemblies with Idba_ud 41 , Supplementary Fig. 8C). The latter might be expected to lead to the "cleanest" genome assembly because it includes the minimum between-sample sequence variation, lowering the probability of generating chimeric contig [START_REF] Mavromatis | Use of simulated data sets to evaluate the fidelity of metagenomic processing methods[END_REF] . The former may be necessary if the virus is locally rare, so that sequences from multiple metagenomes are needed to achieve complete genome coverage. Thus, if the assembly from the single "highest coverage sample" was improved or equivalent to the initial assembly (longest contig in the new assembly representing ≥95% of the longest contig in the initial assembly), this set of contigs was selected as the sequence for this bin (n=6,423). This optimal single-sample assembly was thus privileged compared to a cross-assembly (either based on the initial contigs or on the re-assembly of all sequences aligned to that bin). Otherwise, the "all samples" bin re-assembly was selected if equivalent or better than the initial assembly (longest contig representing ≥95% of the longest initial contig, n=999). The assumption that cross-assembly would be needed for locally rare viruses, without a highcoverage sample, was confirmed by the comparison between the highest coverage of these two types of bins: on average, bins for which the "optimal" assembly were selected displayed a maximum coverage of 5.47 per Gb of metagenome, while the bins for which the "cross-assembly" was selected displayed a maximum coverage of 1.37 per Gb of metagenome (Supplementary Table 2). Finally, if both reassemblies yielded a longest contig smaller (<95%) than the one in the initial assembly, the bin was considered as a false positive (i.e. binning of contigs from multiple genomes, n=1,356), and contigs from the initial assembly were considered as "unbinned" (263,006 contigs, added to the 623,665 contigs ≥1.5kb initially retained as "unbinned").

Identification of viral contigs and delineation of viral populations

Despite efforts to completely remove cellular DNA during sample preparation, the resulting viral metagenomic datasets will only ever be enriched for viruses [START_REF] Roux | Assessment of viral community functional potential from viral metagenomes may be hampered by contamination with cellular sequences[END_REF] . Thus, assembled sequences in the GOV dataset were in silico filtered a posteriori to identify and remove clearly non-viral signal. In this way, our purification methods should have greatly enriched for viruses, but the in silico decontamination step served as a back-up for problematic samples. Together these two "filters" mean that virtually no known cellular signal should have been considered in our analyses. For the in silico cleaning step, VirSorter [START_REF] Roux | VirSorter: mining viral signal from microbial genomic data[END_REF] was used to identify and remove microbial contigs using the "virome decontamination" mode, with every contig ≥10kb and not identified as a viral contig being considered as a microbial contig. Sequences with a prophage predicted were manually curated to distinguish actual prophages (i.e. viral regions within a microbial contig) from contigs that belonged to a viral genome and were wrongly predicted as a prophage. Contigs originating from an eukaryotic virus were identified based on best BLAST hit affiliation of the contig predicted genes against NCBI RefseqVirus (see Supplementary Text).

The genome bins were affiliated as microbial (if 1 or more contigs were identified as microbial, n=1,763), eukaryotic virus (if contigs affiliated as eukaryotic virus comprised more than 10kb or more than 25% of the genome bin total length, n=962) or viral (i.e. archaeal and bacterial viruses, n=4,341), with the 356 remaining bins, lacking a contig long enough for an accurate affiliation, considered as "unknown" (see Supplementary Text).

Viral bins were then refined to evaluate if they corresponded to a single or a mix of viral population(s). To that end, the Pearson correlation and Euclidean distance between abundance profiles (i.e. profile of the contig average coverage depth across the 104 samples) of bin members and the bin seed (i.e. the largest contig) were computed, and a single-copy viral marker gene (TerL) was identified in binned contigs (Supplementary Fig. 8E). Thresholds were chosen to maximize the number of bins with exactly one TerL gene and minimize the number of bins with multiple TerL genes (Supplementary Fig. 8G). For each bin, contigs with a Pearson correlation coefficient to the bin seed <0.96 or a Euclidean distance to the seed >1.05 were removed from the bin, and added to the pool of unbinned contigs. Eventually, every bin still displaying multiple TerL genes after this refinement step were split, and all corresponding contigs added to the pool of "unbinned" contigs (Supplementary Fig. 8E).

The final set of contigs was formed by compiling (i) all contigs belonging to a viral bin, (ii) "unbinned" viral contigs (i.e. contigs affiliated to archaeal and bacterial virus and not part of any genome bin), and (iii) viral contigs identified in microbial or eukaryote virus bins (considered as "unbinned" contigs, Supplementary Fig. 8F). Within this set of contigs, all viral bins were considered as viral populations, as well as every unbinned viral contig ≥10kb, leading to a total of 15,222 epi-and mesopelagic populations, and 58 bathypelagic populations (Supplementary Fig. 1, Supplementary Table 2, and Supplementary Text). In this study, we focus only on the 15,222 epi-and mesopelagic populations, totaling 24,353 contigs. For the detection of AMGs, we added to these populations all short epi-and mesopelagic unbinned viral contigs (<10kb), adding up to a total of 298,383 contigs.

Sequence clustering and annotations

Dataset of publicly available viral genomes and genome fragments

Genomes of viruses associated with a bacterial or archaeal host were downloaded from NCBI RefSeq (1,680 sequences, v70, 05-26-2015). To complete this dataset of reference genomes, viral genomes and genome fragments available in Genbank but not in RefSeq were downloaded (July 2015) and manually curated to select only bacterial and archaeal viruses (1,017 sequences). These included viral genomes not yet added to RefSeq, as well as genome fragments from fosmid libraries generated from seawater samples [START_REF] Mizuno | Expanding the marine virosphere using metagenomics[END_REF][START_REF] Chow | Combining genomic sequencing methods to explore viral diversity and reveal potential virus-host interactions[END_REF] . Mycophage sequences (available from http://phagesdb.org [START_REF] Pope | Whole genome comparison of a large collection of mycobacteriophages reveals a continuum of phage genetic diversity[END_REF] ) were downloaded (July 2015) and included as well if not already in RefSeq (734 sequences). Finally, 12,498 viral genome fragments from the VirSorter Curated Dataset, identified in publicly available microbial genome sequencing projects, were added to the database [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] .

Genome (fragments) clustering through gene-content based network analysis

Proteins predicted from 14,650 large GOV contigs (≥10kb and ≥10 genes), were added to all proteins from the publicly available viral genomes and genomes fragments gathered, and compared through all-vs-all blastp, with a threshold of 10 -5 on e-value and 50 on bit score. Protein clusters were then defined using MCL (using default parameters for clustering of proteins, similarity scores as log-transformed e-value, and 2 for MCL inflation [START_REF] Enright | An efficient algorithm for large-scale detection of protein families[END_REF] ). vContact (https://bitbucket.org/MAVERICLab/ vcontact) was then used to calculate a similarity score between every pair of genome and/or contigs based on the number shared of PCs between the two sequences (as in [START_REF] Lima-Mendez | Reticulate representation of evolutionary and functional relationships between phage genomes[END_REF][START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] ), and then compute a MCL clustering of the genomes/contigs based on these similarity scores (thresholds of 1 on similarity score, MCL inflation of 2). The resulting viral clusters (or VCs, clusters including ≥2 contigs and/or genomes), consistent with a clustering based on whole-genome BLAST comparison, corresponded to approximately genus-level taxonomy, with rare cases closer to subfamily-level taxonomy (Extended Data Fig. 2 and Supplementary Text). A total of 1,259 viral clusters were obtained, with 867 including at least one GOV sequence. Notably, however, automatically defined VCs merely serve as a starting place for assigning viral taxonomy. Current ICTV convention for formal taxonomic consideration of these VCs would require manual comparison of genomes and genome fragments to identify signature genes, comparison of phylogenetic signals, and ideally observation of morphological features of corresponding viruses, although this process is currently being reviewed as advanced computational analytics and genome datasets, such as those presented here, are being developed.

Viral contigs annotation

A functional annotation of all GOV predicted proteins was based on a comparison to the PFAM domain database (v27 [START_REF] Finn | Pfam: the protein families database[END_REF] ) with HmmSearch [START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] (threshold of 30 on bit score and 1e-3 on e-value), and additional putative structural proteins were identified through a BLAST comparison to protein clusters detected in viral metaproteomics dataset [START_REF] Brum | Illuminating structural proteins in viral "dark matter" with metaproteomics[END_REF] . This metaproteomics dataset led to the annotation of 13,547 hypothetical proteins lacking a PFAM annotation. A taxonomic annotation was performed based on a blastp of the predicted proteins against proteins from archaeal and bacterial viruses from NCBI RefSeq and Genbank (threshold of 50 on bit score and 10 -3 on e-value).

VCs were affiliated based on isolate genome members, when available. When multiple isolates were included in the VC, the VC was affiliated to the corresponding subfamily or genus of these isolates (excluding all "unclassified" cases). This was the case for VC_2 (T4 subfamily [START_REF] Sullivan | Genomic analysis of oceanic cyanobacterial myoviruses compared with T4-like myoviruses from diverse hosts and environments[END_REF][START_REF] Zhao | Abundant SAR11 viruses in the ocean[END_REF] ), and VC_9 (T7virus [START_REF] Labrie | Genomes of marine cyanopodoviruses reveal multiple origins of diversity[END_REF] ). When only one or a handful of affiliated isolate genomes were included in the VC and lacked genus-level classification, a candidate name was derived from the isolate (if several isolates, from the first one isolated). This was the case for VC_5 (Cbaphi381virus 54 ), VC_12 (P12024virus [START_REF] Kang | Complete genome sequences of two Persicivirga bacteriophages, P12024S and P12024[END_REF] ), VC_14 (MED4-117virus), VC_19 (HMO-2011virus [START_REF] Kang | Genome of a SAR116 bacteriophage shows the prevalence of this phage type in the oceans[END_REF] ), VC_31 (RM378virus 57 ), VC_36 (GBK2virus 58 ), VC_47 (Cbaphi142virus 54 ) , and VC_277 (vB_RglS_P106Bvirus [START_REF] Halmillawewa | Genomic and phenotypic characterization of Rhizobium gallicum phage vB_RglS_P106[END_REF] ). Otherwise, VCs were considered as "new VCs".

"Phage proteomic tree" (i.e. "whole-genome comparison tree") computation and visualization

All publicly available complete genomes (see above), all complete (circular) and near-complete (extrachromosomal genome fragment >50kb with a terminase) from the VirSorter Curated Dataset, and all complete and near-complete GOV contigs were compared to generate a phage proteomic tree, as previously described [START_REF] Mizuno | Expanding the marine virosphere using metagenomics[END_REF][START_REF] Rohwer | The Phage Proteomic Tree : a Genome-Based Taxonomy for Phage[END_REF] . Briefly, a proteomic similarity score was calculated for each pair of genome based on a all-vs-all tblastx similarity as the sum of bit scores of significant hits between two genomes (e-value ≤ 0.001, bit score ≥30, identity percentage ≥ 30). To normalize for different genome sizes, each genome was also compared to itself to generate a self-score, and the distance between two different genomes was calculated as a Dice coefficient (as in [START_REF] Mizuno | Expanding the marine virosphere using metagenomics[END_REF] ), i.e. for two genomes A and B with a proteomic similarity score of AB, the corresponding distance d would be 1-(2*AB)/(AA+BB), with AA and BB being the self-score of genomes A and B respectively. For clarity, the tree displayed in Extended Data Fig. 2 only include non-GOV sequences found in a VC with GOV sequence(s) or within a distance <0.5 to a GOV sequence, adding for a total of 1,522 reference sequences. iTOL [START_REF] Letunic | Interactive Tree Of Life (iTOL): an online tool for phylogenetic tree display and annotation[END_REF][START_REF] Letunic | Interactive Tree Of Life v2: online annotation and display of phylogenetic trees made easy[END_REF] was used to visualize and display the tree.

Distribution and relative abundance of viral populations and VCs Detection and estimation of abundance for viral contigs and populations

The presence and relative abundance of a viral contig in a sample were determined based on the mapping of HQ reads to the contig sequences, computed with Bowtie 2 (options --sensitive, -X 2000, and --non-deterministic, default parameters otherwise [START_REF] Langmead | Fast gapped-read alignment with Bowtie 2[END_REF] ), as previously described [START_REF] Brum | Patterns and ecological drivers of ocean viral communities[END_REF] . A contig was considered as detected in a metagenome if more than 75% of its length was covered by aligned reads derived from the corresponding sample. A normalized coverage for the contig was then computed as the average contig coverage (i.e. number of nucleotides mapped to the contig divided by the contig length) normalized by the total number of bp sequenced in this sample. The detection and relative abundance of a viral population was based on the coverage of its contigs: a population was considered as detected in a sample if more than 75% of its cumulated length was covered, and its normalized coverage was computed as the average normalized coverage of its contigs.

Relative abundance of VCs

The relative abundance of VCs was calculated based on the coverage of its members within the 15,222 viral populations identified. If a population included contigs all linked to the same VC, or linked to a single VC except for unclustered (because too short) contigs, this population coverage was added to the total of the corresponding VC. In the rare cases where the link between population and VC was ambiguous because different contigs within a population pointed toward different VCs (n=475, i.e. 3.1% of the populations), the population coverage was equally split between these VCs. Finally, if no contig in the population belonged to any VC (n=2,605, 17% of the populations), the population coverage was added to the "unclustered" category. Eventually, for each sample, the cumulated coverage of a VC was normalized by the total coverage of all populations to calculate a relative abundance of the VC among viral populations.

The selection of abundant VCs within a sample was based on the contribution of the VC to the sample diversity as measured by the Simpson index. For each sample, the overall Simpson index was first calculated with all VCs. Then, VCs were sorted by decreasing relative abundance and progressively added to a new calculation of the Simpson index. VCs considered as abundant were the ones which, once cumulated, represented 80% of the sample diversity (i.e. a Simpson index greater or equal to 80% of the sample total Simpson index, Extended Data Fig. 1C). The 38 VCs identified as abundant in at least 2 different stations were selected as "recurrently abundant VCs in the GOV dataset" (Fig. 2 and Extended Data Fig. 3).

Host prediction and diversity

Three different approaches were used to link viral contigs and putative host genomes: blastn similarity, CRISPR spacer similarity, and tetranucleotide frequencies similarities. An overview of the contigs generation process is provided in Supplementary Fig. 8, and an extended discussion about the efficiency and raw results of these host prediction methods is provided in Supplementary Text, Supplementary Table 4, and ref. [START_REF] Edwards | Computational approaches to predict bacteriophage-host relationships[END_REF] . A list of all host predictions by viral sequence is available in Supplementary Table 5.

Generation of host database

A genome database of putative hosts for the epi-and mesopelagic GOV viruses was generated, including all archaea and bacteria genomes annotated as "marine" from NCBI RefSeq and WGS (both times only sequences ≥5kb, 184,663 sequences from 4,452 genomes, downloaded in August 2015), and all contigs ≥5kb from the 139 Tara Oceans microbial metagenomes corresponding to the bacteria and archaea size fraction (791,373 sequences) [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] . For these microbial metagenomic contigs, a first blastn was computed to compare them to all GOV contigs, and exclude from the putative host dataset all metagenomic contigs with a significant similarity to a viral GOV sequence (thresholds of 50 on bit score, 0.001 on e-value, and 70% on identity percentage) on ≥90% of their length, as these are likely sequences of viral origin sequenced in the bacteria and archaea size fraction (these represented 2.2% of the contigs in the assembled microbial metagenomes). The taxonomic affiliation of NCBI genomes was taken from the NCBI taxonomy. For Tara Oceans contigs, a last common ancestor (LCA) affiliation was generated for each contig based on genes affiliation [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] , if 3 genes or more on the contig were affiliated.

BLAST-based identification of sequence similarity between viral contigs and host genome

All GOV viral contigs were compared to all archaeal and bacterial genomes and genome fragments with a blastn (threshold of 50 on bit score and 0.001 on e-value), to identify regions of similarity between a viral contig and a microbial genome, indicative of a prophage integration or horizontal gene transfer [START_REF] Edwards | Computational approaches to predict bacteriophage-host relationships[END_REF] . A host prediction was made when (i) a NCBI genomes displayed a region similar to a GOV viral contig ≥5kb at ≥70% id, or (ii) when a Tara Oceans microbial metagenomic contig (≥5kb) displayed a region similar to a GOV viral contig ≥2.5kb at ≥70% id.

Matches between GOV viral contigs and CRISPR spacers.

CRISPR arrays were predicted for all putative host genomes and genome fragments (NCBI microbial genomes and Tara Oceans microbial metagenomic contigs) with MetaCRT [START_REF] Bland | CRISPR recognition tool (CRT): a tool for automatic detection of clustered regularly interspaced palindromic repeats[END_REF][START_REF] Rho | Diverse CRISPRs Evolving in Human Microbiomes[END_REF] . CRISPR spacers were extracted, and all spacers with ambiguous bases or low complexity (i.e. consisting of 4 to 6 bp repeat motifs) were removed. All remaining spacers were matched to viral contigs with fuzznuc [START_REF] Rice | EMBOSS: the European Molecular Biology Open Software Suite[END_REF] , with no mismatches allowed, which although rarely observed yields highly accurate host predictions [START_REF] Edwards | Computational approaches to predict bacteriophage-host relationships[END_REF] ( Supplementary Table 4).

Nucleotide composition similarity: comparison of tetranucleotide frequency

Bacterial and archaeal viruses tend to have a genome composition close to the genome composition of their host, a signal that can be used to predict viral-host pairs [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF][START_REF] Edwards | Computational approaches to predict bacteriophage-host relationships[END_REF][START_REF] Ogilvie | Genome signature-based dissection of human gut metagenomes to extract subliminal viral sequences[END_REF] . Here, canonical tetranucleotide frequencies were observed for all viral and host sequences using Jellyfish [START_REF] Marçais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF] , and mean absolute error (i.e. average of absolute differences) between tetranucleotide frequency vectors were computed with in-house Perl and Python scripts for each pair of viral and host sequence as in ref. [START_REF] Roux | Viral dark matter and virus-host interactions resolved from publicly available microbial genomes[END_REF] . A GOV viral contig was then assigned to the closest sequence (i.e. lowest distance d) from the pool of NCBI genomes if d<0.001 (because both the tetranucleotide frequency signal and the taxonomic affiliation of these complete genomes are more robust than for metagenomic contigs), and otherwise assigned to the closest (i.e. lowest distance) Tara Oceans microbial contig if d<0.001.

Summarizing host prediction at the VC level

Overall, 3,675 GOV contigs could be linked to a putative host group among the 24,353 GOV contigs associated with an epi-or mesopelagic viral population. To summarize these affiliations at the VC level, a Poisson distribution was used to estimate the number of expected false positive associations for each VC -host group combination based on (i) the global probability of obtaining a host prediction across all pairs of viral and host sequences tested and for all methods (5.8x10 -08 ), (ii) the number of potential predictions generated for the VC, corresponding to 3 times the number of sequences in the VC (to take into account the three methods), and (iii) the number of sequences from the host group in the database (Supplementary Figure 2). By comparing the number of links observed between a VC and a host group to this expected value, which takes into account the bias in database (i.e. some host groups will be over-or under-represented in our set of archaeal and bacterial genomes and genome fragments) and the bias linked to the variable number of sequences in VCs, we can determine if the number of associations observed for any VC -host group combination is likely to be due to chance alone (and calculate the associated p-value).

Microbial community diversity and richness indexes

Diversity and richness indexes for putative host populations were based on the OTU abundance matrix generated from the analysis of mi TAGs in Tara Oceans microbial metagenomes [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] . These indexes were computed for each host group at the same taxonomic level as the host prediction, i.e. the phylum level except for Proteobacteria where the class level is used. The R package vegan [START_REF] Oksanen | The vegan Package[END_REF] was used to estimate for each group (i) a global Chao index (i.e. including all OTUs from all samples) through the function estaccumR, (ii) a sample-by-sample Chao index with the function estimateR, and (iii) Sorensen indexes between all pairs of samples with the function betadiver. Diversity indexes presented in Extended Data Fig 4 are based on epipelagic samples only, as the 38 VCs identified as abundant were mostly retrieved in epipelagic samples. Candidate division OP1 was excluded from this analysis because no OTU affiliated to this phylum was identified.

Identification and annotation of putative AMGs Detection of AMGs

Predicted proteins from all GOV viral contigs were compared to the PFAM domain database (hmmsearch [START_REF] Eddy | Accelerated Profile HMM Searches[END_REF] , threshold of 40 on bit score and 0.001 on e-value), and all PFAM domains detected were classified into 8 categories: "structural", "DNA replication, recombination, repair, nucleotide metabolism", "transcription, translation, protein synthesis", "lysis", "membrane transport, membraneassociated", "metabolism", "other", and "unknown" (as in ref. [START_REF] Hurwitz | Depth-stratified functional and taxonomic niche specialization in the "core" and "flexible" Pacific Ocean Virome[END_REF] ). Four AMGs (i.e. similar to a domain from the "metabolism" category) were then selected for further study because of their central role in sulfur (dsrC and soxYZ) or nitrogen (P-II, amoC) cycle, and the fact that these had never been detected in a surface ocean viral genome so far (dsrC/tusE-like genes have been detected in deep water viruses [START_REF] Roux | Ecology and evolution of viruses infecting uncultivated SUP05 bacteria as revealed by single-cell-and meta-genomics[END_REF][START_REF] Anantharaman | Sulfur Oxidation Genes in Diverse Deep-Sea Viruses[END_REF] ). To evaluate if an AMG was "known", a list of PFAM domain detected in NCBI RefSeqVirus and Environmental Phages was computed based on a similar hmmsearch comparison (threshold of 40 on bit score and 0.001 on e-value), and augmented by manual annotation of AMGs from [START_REF] Hurwitz | Depth-stratified functional and taxonomic niche specialization in the "core" and "flexible" Pacific Ocean Virome[END_REF][START_REF] Sharon | Comparative metagenomics of microbial traits within oceanic viral communities[END_REF] . These corresponded for the most part to photosynthesis and carbon metabolism AMGs previously described in cyanophages [START_REF] Thompson | Phage auxiliary metabolic genes and the redirection of cyanobacterial host carbon metabolism[END_REF][START_REF] Dammeyer | Efficient phage-mediated pigment biosynthesis in oceanic cyanobacteria[END_REF][START_REF] Lindell | Photosynthesis genes in marine viruses yield proteins during host infection[END_REF][START_REF] Lindell | Genome-wide expression dynamics of a marine virus and host reveal features of co-evolution[END_REF][START_REF] Sullivan | Prevalence and Evolution of Core Photosystem II Genes in Marine Cyanobacterial Viruses and Their Hosts[END_REF] . The complete list of PFAM domains detected in GOV viral contigs is available as Supplementary Table 6.

Phylogenetic tree generation and contigs map comparison

Sequences similar to these AMGs were recruited from the Tara Oceans microbial metagenomes 18 based on a blastp of all predicted proteins from microbial metagenome to the viral AMGs identified (threshold of 100 on bit score, 10 -5 on e-value, except for P-II where a threshold of 170 on bit score was used because of the high number of sequences recruited). The viral AMG sequences were also compared to NCBI nr database (blastp, threshold of 50 on bit score and 10 -3 on e-value) to recruit relevant reference sequences (up to 20 for each viral AMG sequence). These sets of viral AMGs and related protein sequences were then aligned with Muscle [START_REF] Edgar | MUSCLE: a multiple sequence alignment method with reduced time and space complexity[END_REF] , the alignment manually curated to remove poorly aligned positions with Jalview 77 , and two trees were computed from the same curated alignment: a maximum-likelihood tree with FastTree (v2.7.1, model WAG, other parameters set to default [START_REF] Price | FastTree 2--approximately maximum-likelihood trees for large alignments[END_REF] ) and a bayesian tree with MrBayes (v3.2.5, mixed evolution models, other parameters set to default, 2 MCMC chains were run until the average standard deviation of split frequencies was <0.015, relative burn-in of 25% used to generate the consensus tree [START_REF] Huelsenbeck | MRBAYES: Bayesian inference of phylogenetic trees[END_REF] ). In all cases except AmoC, the mixed model used by MrBayes was 100% WAG, confirming that this model was well suited for archaeal and bacterial virus protein trees. Manual inspection revealed only minor differences between each pair of trees, so an SH test was used to determine which tree best fitted the sequence alignment, using the R library phangorn [START_REF] Schliep | phangorn: phylogenetic analysis in R[END_REF] . Itol 61 was used to visualize and display these trees, in which branches with supports <40% were collapsed. Annotated interactive trees are available online at http://itol.embl.de/shared/Siroux. Contigs map comparison were generated with Easyfig 81 , following the same method as for the VCs (see Supplementary Information).

Functional characterization of putative AMGs

Conserved motifs were identified on the different AMGs based on the literature: dsrC conserved motifs were obtained from ref. [START_REF] Venceslau | The "bacterial heterodisulfide" DsrC is a key protein in dissimilatory sulfur metabolism[END_REF] , soxYZ conserved residues were identified from the PFAM domains PF13501 and PF08770, and P-II conserved motifs from PROSITE documentation PDOC00439. A 3D structure could also be predicted for P-II AMGs by I-TASSER [START_REF] Roy | I-TASSER: a unified platform for automated protein structure and function prediction[END_REF] (default parameters), the quality of these predictions being confirmed with ProSA web server [START_REF] Wiederstein | ProSA-web: interactive web service for the recognition of errors in three-dimensional structures of proteins[END_REF] . To further confirm the functionality of these genes, selective constraint on these AMGs was evaluated through pN/pS calculation, as in ref. [START_REF] Schloissnig | Genomic variation landscape of the human gut microbiome[END_REF] . Briefly, synonymous and non-synonymous SNPs were observed in each AMG, and compared to expected ratio of synonymous and non-synonymous SNPs under a neutral evolution model for this genes. The interpretation of pN/pS is similar as for dN/dS analyses, with the operation of purifying selection leading to pN/pS values < 1. Finally, AMG transcripts were searched in metatranscriptomic datasets generated through the Tara Oceans consortium (ENA Id ERS1092158, ERS488920, and ERS494518). For generating these metatranscriptomes, bacterial rRNA depletion was carried out on 240-500 ng total RNA using Ribo-Zero Magnetic Kit for Bacteria (Epicentre, Madison, WI) for 0.2-1.6 and 0.22-3µm filters. The Ribo-Zero depletion protocol was modified to be adapted to low RNA input amounts [START_REF] Alberti | Comparison of library preparation methods reveals their impact on interpretation of metatranscriptomic data[END_REF] . Depleted RNA was used to synthetize cDNA with SMARTer Stranded RNA-Seq Kit (Clontech, Mountain View, CA) [START_REF] Alberti | Comparison of library preparation methods reveals their impact on interpretation of metatranscriptomic data[END_REF] . Metatranscriptomic libraries were quantified by qPCR using the KAPA Library Quantification Kit for Illumina Libraries (KapaBiosystems, Wilmington, MA) and library profiles were assessed using the DNA High Sensitivity LabChip kit on an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA). Libraries were sequenced on Illumina HiSeq2000 instrument (Illumina, San Diego,CA) using 100 base-length read chemistry in a paired-end mode. High quality reads were then mapped to viral contigs containing dsrC, soxYZ, P-II, or amoC genes with SOAPdenovo2 43 within MOCAT [START_REF] Kultima | MOCAT: a metagenomics assembly and gene prediction toolkit[END_REF] (options screen and filter with length and identity cutoffs of 45 and 95%, respectively, and paired-end filtering set to yes), and coverage was defined for each gene as the number of bp mapped divided by gene length (including only reads mapped to the predicted coding strand).

Distribution of AMGs and association with geochemical metadata

The distribution and relative abundance of AMGs was based on the read mapping and normalized coverage of the contig including the AMG. To get a range of temperature and nutrient concentrations for the widespread AMGs (detected in >5 stations) that takes into account both the samples in which these AMGs were detected and the differences in normalized coverage, a set of samples was selected through a weighted random drawing replacement, with the weight of each sample corresponding to the AMG's normalized coverage. That way, a range of temperature or nutrient concentration values associated with the AMG's distribution and abundance could be generated for each AMG and each environmental parameter tested. The number of samples randomly selected for each AMG was the same as the total number of samples for which a value of this parameter was available.

Code and data availability

Scripts used in this manuscript are available on the Sullivan lab bitbucket under project "GOV_Ecogenomics" (http://bitbucket.org/MAVERICLab/gov_ecogenomics/overview). Scripts used in the assessment of microbial diversity are gathered in the directory "Host_diversity", the ones used by the AMG coverage across these samples (see Extended Data Fig. 9 for underlying coverage data). Distributions significantly different from the "All Samples" distribution (two-sided KS-test) are indicated with stars. Boxes represent the first and third quartiles around the median.

Extended Data Figure 1: Accumulation curves of populations (A) and viral clusters (VCs, B) and identification of abundant VCs in GOV samples (C). A & B.

Accumulation curves were computed from 50 randomly shuffled samples (blue dots), with all, epipelagic, mesopelagic, or bathypelagic subsets of the data. For each curve, the average of 50 iterations is displayed with red dots. C. Schematic of the selection process of abundant VCs. For each sample, VCs accounting for (up to) 80% of the sample diversity (as assessed by Simpson index) were considered as abundant (example for sample 125_MIX on the left). VCs detected as abundant in at least two different stations were included in the 38 VCs described in Fig. 2 and Extended Data Fig. 3. Extended Data Figure 2: Comparison of VCs with other classification methods: phage proteomic tree and percentage of shared genes. The phage proteomic tree includes the 756 GOV complete and near-complete genomes from epi-and mesopelagic samples, and closest references from RefSeq and Environmental phages (d<0.5 to a GOV sequence or found in the same VC as a GOV sequence). Branches of monophyletic clades including more than 3 GOV and/or uncultivated marine sequences with no isolate reference are highlighted in blue. All VCs with more than 8 representatives in the tree or part of the 38 abundant VCs are indicated with coloring of the outer ring. The name and affiliation (if available) of the 38 abundant VCs are indicated next to the VC on the colored ring. VCs whose members were gathered in a single monophyletic clades are indicated with a solid black outline, while VCs for which all but one members were gathered in a single monophyletic clades are highlighted with a dashed black outline. Inset: distribution of number of shared genes estimated based on the number of shared PCs (protein clusters) for viral genome/contigs pairs either between different VCs or within VCs. On average, 73% and 39% of sequences within a VC shared more than 20% and 40% of their genes, respectively, which represent the current thresholds currently accepted for sub-family and genus designations. Similarly, 83% of sequences within a VC were consistently affiliated in the phage proteomic tree as they formed a monophyletic group including only members of the particular VC. Thus all three classification methods are largely consistent for the GOV dataset (see Supplementary Text). Extended Data Figure 3: Summary of 34 of the 38 abundant viral clusters (VCs, the 4 other abundant VCs being the ubiquitous ones presented in Fig. 2). Predicted genome size is based on the set of isolates and circular contigs in the VC (NA corresponds to VCs without any circular contigs, or for which the relative standard deviation of estimated genome size across the different isolate(s) and/or circular contigs is greater than 15%). Host association values are based on the number of cluster members associated with each host group, the statistical significance of this number of predictions being evaluated by comparison with an expected number of associations calculated from a Poisson distribution. Host associations based on known isolates are indicated with a star (for associations based on cultivated isolates) or a dot (for associations based on the detection of a cluster member in a microbial genome from the VirSorter Curated Dataset). The abundant epipelagic microbial groups (representing >1% of the microbial OTUs abundance of epipelagic samples) are highlighted in bold. Distribution and relative abundance of VCs are based on the cumulated coverage of VC members among sample viral populations. The main oceanic basins are indicated for each set of sample, Med. Sea-Mediterranean Sea. Extended Data Figure 4: Association between abundant viral clusters (VCs) and host group abundance and diversity A. Abundance and diversity of bacterial and archaeal host groups associated with the 38 abundant VCs (see Fig. 2A). For each host group (phylum level, except for Proteobacteria where the class level is used), the different panels display from top to bottom (i) the number of VCs associated with this host group, (ii) the global relative abundance of this group estimated from the microbial metagenomic OTU counts, (iii) the global diversity of this group based on a Chao index computation including all Tara Oceans microbial metagenome samples (i.e. including both Alpha and Beta diversity), (iv) the distribution of Chao indexes by sample for this group (Alpha diversity), and (v) the average Sorensen index between pairs of samples including at least one OTU of this group (Beta diversity). OTU counts were derived from the 109 epipelagic microbial metagenomes described in [START_REF] Sunagawa | Structure and function of the global ocean microbiome[END_REF] . B. Pearson correlations between host group relative abundance or diversity indexes (Global Chao, Average Chao across samples, and Average Sorensen across samples) and the number of VCs. Extended Data Figure 5: Diversity, distribution, and genome context of dsrC genes in GOV contigs. A. Maximum-likelihood tree (from an amino-acid alignment) including the 11 viral DsrC and microbial sequences from microbial metagenomes and NCBI nr database. The presence of conserved C residues (named Cys-A & Cys-B, as in ref. [START_REF] Venceslau | The "bacterial heterodisulfide" DsrC is a key protein in dissimilatory sulfur metabolism[END_REF] ) is indicated with color circles next to each sequence or clade, and the corresponding type of DsrC-like protein is indicated by coloring the branch or clade. The microbial metagenomic contigs affiliated to uncultivated, marine sulfur-oxidizing Gammaproteobacteria (as confirmed by complementary phylogenetic analysis of DsrAB, Supplementary Fig. 7) are indicated with a star next to the sequence or clade. Viral AMG sequences are highlighted in blue, internal nodes SH-like supports are represented by proportional circles (all nodes with support < 0.40 were collapsed). Each dsrC AMG is associated with an abundance profile (on the right) displaying the relative abundance of the contig across the 91 epi-and mesopelagic samples (based on normalized coverage, i.e. contig coverage / Gb of metagenome). B. Comparison of dsrCcontaining contigs maps. T4-like marker gene (T4 baseplate) is indicated on the maps, alongside putative AMGs (Fe-S biosyn for Iron-sulfur cluster biosynthesis, and Amt for Ammonia transporter). Extended Data Figure 6: Diversity, distribution, and genome context of soxYZ genes in GOV contigs. A. Bayesian tree (from an amino-acid alignment) including the 4 viral SoxYZ and microbial sequences from microbial metagenomes and NCBI nr database. The affiliation of microbial clades (either from the NCBI reference or from the LCA affiliation of metagenomic contigs) is indicated by coloring of the grouped clades or with a colored square next to the sequence. Viral AMG sequences are highlighted in blue, posterior probabilities are represented by proportional circles (all nodes with posterior probability < 0.40 were collapsed). Clades including sulfur-oxidation proteobacteria are indicated on the tree. Each soxYZ AMG is associated with an abundance profile (on the right) displaying the relative abundance of the contig across the 91 epi-and mesopelagic samples (based on normalized coverage, i.e. contig coverage / Gb of metagenome). B. Comparison of soxYZ-containing contigs maps. For contig GOV_bin_4310_contig-100_0, the second largest contig from the same bin (GOV_bin_4310_contig-100_1) is displayed. T4-like marker genes (Gp23 and T4 baseplate) are indicated on the maps, alongside putative AMGs (Fe-S biosyn: Iron-sulfur cluster biosynthesis). Extended Data Figure 7: Diversity, distribution, and genome context of P-II genes in GOV contigs. A. Maximum-likelihood tree (from an amino-acid alignment) including the 10 viral P-II and microbial sequences from microbial metagenomes and NCBI nr database. The affiliation of microbial clades (either from the NCBI reference or from the LCA affiliation of metagenomic contigs) is indicated by coloring of the grouped clades or with a colored square next to the sequence. The sequences lacking the conserved uridylation site of P-II (Supplementary Fig. 5) are highlighted with a star next to the sequence name or clade. Viral AMG sequences are highlighted in blue, internal nodes SH-like supports are represented by proportional circles (all nodes with support < 0.40 were collapsed). Each P-II AMG is associated with an abundance profile (on the right) displaying the relative abundance of the contig across the 91 epi-and mesopelagic samples (based on normalized coverage, i.e. contig coverage / Gb of metagenome). B. Comparison of P-II-containing contigs maps. Ammonia transporter genes linked to P-II are indicated on the map (Amm Transp, dark red). When available, the VC affiliation of each contig is indicated next to the contig name. Contig GOV_bin_5834_contig-100_7 is too short to be clustered based on a shared PC network, however the seed contig of its population was clustered (in VC_12, Siphoviridae -P12024virus), hence this seed contig affiliation is indicated. Extended Data Figure 8: Diversity, distribution, and genome context of amoC gene in GOV contigs. A. Maximum-likelihood tree (from an amino-acid alignment) including the GOV amoC AMG and microbial sequences from microbial metagenomes and NCBI nr database. The affiliation of microbial clades (either from the NCBI reference or from the LCA affiliation of metagenomic contigs) is indicated by coloring of the grouped clades or with a colored square next to the sequence. Viral AMG sequence is highlighted in blue, internal nodes SH-like supports are represented by proportional circles (all nodes with support < 0.40 were collapsed). B. Abundance profile displaying the relative abundance of the contig across the 91 epi-and mesopelagic samples (based on normalized coverage, i.e. contig coverage / Gb of metagenome). C. Map of the amoC-containing contig. Extended Data Figure 9: Normalized coverage of contigs harboring AMG as function of the temperature and nutrient concentrations (NO 2 , NO 3 , PO 4 ) of the corresponding samples. AMGs are grouped by clade based on the phylogeny (see Extended Data Fig. 5-6-7), and coverages are cumulated when a clade included multiple contigs. Plots display the cumulated normalized coverage of a clade (y-axis) as function of the temperature or nutrient concentration (x-axis) across all epipelagic samples (mesopelagic samples were excluded from the analysis since the AMG signal was detected in epipelagic samples), only for clades not geographically restricted (i.e. found in >5 samples, see Fig. 3C). Samples are color-coded according to their ocean and sea region (Supplementary Table 1). The calculated preferential range of temperature or nutrient concentration is displayed below each plot for the epipelagic AMGs (P-II-4 distribution could not be linked to specific environmental conditions, but this AMG is the only one consistently retrieved in mesopelagic samples). Extended Data Table 1: Summary of genes and contigs characteristics for new viral DsrC, SoxYZ, and P-II AMGs. Each gene is linked to its contig, and when available, to the corresponding viral cluster and predicted host (from BLAST hit, CRISPR spacer similarity, or nucleotide composition similarity, Alphaprot.-Alphaproteobacteria, Gammaprot.-Gammaproteobacteria). Widespread and abundant VCs are highlighted in bold. In addition, the calculated pN/pS of each gene is indicated (measuring the strength of selection pressure occurring for this gene, the gene with a pN/pS not representing a strong purifying selection is highlighted in red), as well as the coverage of these genes and other genes in the contigs in 3 metatranscriptomic samples from 3 open ocean Tara stations (cases where the AMG coverage is >0.5 and associated with the coverage of other genes from the same viral contig are highlighted in green).

Acknowledgements

We thank Joshua Weitz for advice on statistics, Claus Pelikan for help with the DsrAB phylogenetic tree, Christiane Dahl for valuable discussion regarding DsrC function, and members of the Sullivan and Rich labs for suggestions and comments on this manuscript. We acknowledge support from UA high-

for host predictions are in "Host_prediction", and the scripts used to identify abuntant VCs are in "Virus_clusters_prevalence". All raw reads are available through ENA (Tara Oceans) or JGI (Malaspina) using the dataset identifiers listed in Supplementary Table 1. Processed data are available through iVirus (http://mirrors.iplantcollaborative.org/browse/iplant/home/shared/ivirus/GOV/), including all sequences from assembled contigs, list of viral populations and associated annotated sequences as genbank files, viral clusters composition and characteristics, map comparisons of genomes and contigs of the 38 abundant VCs, and host predictions for viral contigs.


Joanna Daaboul

Catherine Da Cunha

Differentiation and Customer Decoupling Points: Key Value Enablers for Mass Customization

Keywords: Mass customization, value, decoupling point, differentiation points

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Many companies in different fields (Adidas, Nike, DELL, Woonio for customized furniture, Spreadshirt for customized shirts, Louis Vuitton, Motorola, BMW etc.) implemented Mass Customization (MC). It enables to benefit both from costs reduction (mass production) and good response to customers' requirements (customization). Even though MC was the center of many research projects (trying to better identify this strategy, understand its enablers and how it can be implemented successfully), many questions remained unanswered. There are two main decisions when implementing MC:

1. What customization to offer to the customer: which components of the product will be standard and which will be customized, thus where to position the product differentiation points (PDP).

2. How to produce a mass customized product: which processes will be Make-to-Stock (MTS) and which will be Make-to-order (MTO), thus where to position the customer-order-decoupling-point (CODP).

Until now those two questions were often answered separately. This paper advocates that PDP and CODP should be considered simultaneously when defining the best MC customization strategy for a company. The paper is organized as follows: section 2 presents a definition and a literature review for the three main concepts: value, PDP, and CODP. Section three describes the methodology used. Section 4 presents a case study and finally a conclusion in section 5 on the opportunities for MC.

2

Literature review

Value

Based on a literature review combined with an empirical study, [START_REF] Daaboul | Value of personalized products: modelling the customer perception[END_REF] concluded that value can be identified as "the amount of satisfaction created, by fulfilling a certain physical, biological, or psychological need of a beneficiary party". Many criteria (such as cost, delay, perceived quality, and perceived price) influence it.

It can be objective or subjective. It depends on circumstances and is related to the specific goals of the beneficiary party. Different performance indicators are used as its measurement. Value has different beneficiary parties. Those can be the customer, suppliers, enterprise, stakeholders, etc. All those form a value network which is a group of partners collaborating together in order to generate value. A value network, unlike a value chain, is not a linear sequential order of activities transforming materials into products. It allows the consideration of the interactions between the different activities, and relies on the definition of value being multicriteria, subjective and with different beneficiary parties [START_REF] Daaboul | Value network modelling and simulation for strategic analysis: A discrete event simulation approach[END_REF]. For the enterprise, the financial value (its objective value) is measured by its profit. Its subjective value is measured by its image, customer loyalty, ranking among completion, etc. For the customer, value is purely subjective. It is measured by dividing the perceived quality (of the product and related services) to the perceived price.

Product Differentiation Point (PDP)

The product differentiation point denominates an operation that transforms a product common to all the products of the same family to a customized or personalized product (finished or not). This transformation can be done by the adjunction of specialized components and/or by the action of a special process. PDP can be multiple, as many attributes can differentiate products among a family. Figure 1 schematizes the multiple product differentiation points that enable obtaining a final product from common inventories. PDPs are not limited to technological products. For instance in the agro-food industry, the large variety of formats and/or packaging, which creates customization, compensates the relative simplicity of the product (van Donk 2001).

The choice of PDP encompasses the 2 dimensions of MC optimization: product and process. Product differentiation is either studied on the product portfolio angle (how many different options should be offered) or on the position in the chain angle (should the PDP be moved to the customer: postponement). The differentiation may occur in any attribute of the product including price and localization. A good localization of PDP is closely linked not only to the process used, but also to the information describing the demand [START_REF] Zhang | Postponed product differentiation with demand information update[END_REF]. 2013). The CODP position should be determined based on the overall value generated by the system for each beneficiary party (e.g. customers, shareholders, company). In fact, the CODP position influences the performance of the network via its impact on costs (cost of resources, cost of materials, inventory cost, production cost, transportation cost), delays (setup time, processing time, delivery delays, order delay time, etc.), and quality (conforming commands, customer perceived quality of the product and the service, etc.). Its position should be optimized by taking into consideration not only capacity constraints, cost and time, but also the product differentiation points and the customization possibilities offered to the customer.

Coupled analysis

The question of PDP and CODP positioning were often considered as exclusive alternatives. Zinn and Bowersox (1988) addressed both questions; naming the point of differentiation "form postponement" and the customer-order-decoupling-point "time postponement". They concluded that product price is the variable that justifies form postponement while uncertainty is determinant for time postponement. Nevertheless in their model form postponement can only occur after the customer order is received. Su et al (2005) considered form and time postponement as exclusive alternatives.

They have concluded that once the number of products increases above some threshold, the time postponement structure is preferred under both performance metrics (cost and waiting time). According to them, form postponement structure is more favorable in the case of higher arrival time and process time variations. Time postponement is more favorable in the case of higher interest rates. Hegde et al ( 2005) developed a framework to address what degrees of customization to offer. Their results showed that product conformance gradually decreased when the degree of customization exceeded the manufacturing threshold. This paper advocates that the decision on the customization offer should be made simultaneously with the decision on the manufacturing system configuration, and in particular with the CODP determination. [START_REF] Fujita | Mathematical model for simultaneous design of module commonalization and supply chain configuration toward global product family[END_REF] propose a mathematical model for strategy-level simultaneous design of module commonalization and supply chain configuration. They consider necessary to address both decisions simultaneously.

Methodology

The developed decision aid system for value networks consists of an evaluation of a strategic decision (such as moving to MC) considering its impact on generated value for all partners. It is based on a what-if analysis and comparison of different possible scenarios. In step 3 of the methodology shown in Figure 3, possible TO-BE scenarios are defined. In our case, a TO-BE scenario consists of a combination of CODP position and PDPs position. This is realized by identifying possible PDP positions, possible CODP positions and then by identifying all possible combination of these positions.

Step 4 analyses the results.

A specific library was developed by [START_REF] Daaboul | Value networks: pulling the triggers. A combined approach of modelling and simulation for performance evaluation[END_REF] for discrete event simulation software: Arena (Rockwell Automation). It is intended to model and simulate easily and in short time, large and complex value networks. It concentrates on value evaluation. It also permits easy modeling and simulation of different PDPs and CODPs. It is formed of seven modules: order generator, decision activity, trigger, execution activity, partner, variant, and physical flow. The library permits evaluating customer perceived value and value for enterprise. A shoe is formed essentially of an upper, outsole, insole, midsole, eyelets, buckle, laces, and accessories. After analysis and discussions with Alpina, two possible PDPs were identified. The first PDP is at the laced shoe stage, meaning that the customer can only customize the color of the laces. The second PDP is at the upper stage, meaning that the customer may customize the color of its shoe, the type and color of the outsole, and the laces. The shoe production process essentially consists of seven main processes: outsole production, insole and midsole production, cutting upper parts, stitching upper parts, lasting upper and assembling it to the soles, cleaning the shoes, and adding laces. This process leads to five possible CODPs (CODP 1 before outsole production, CODP2 before cutting upper parts, CODP3 before stitching upper parts, CODP4 before lasting upper, and CODP5 before adding laces. Figure 4 shows both the process and the shoe parts.

Simulation of value network of Alpina

Needed data (concerning products, processes, supply chain, distributor, suppliers, etc.) were collected from the enterprise using different templates and many interviews with different departments of Alpina. An empirical study realized by Alpina permitted obtaining the perceived quality attributes (fitting, aesthetics, material, thermal comfort, flexibility, weight, stability, breathability, durability and waterproofing, brand name, ecological level, shop assistance, customization time and number of proposed product variants). Thus all needed data for building the value model was provided by the real world. The performance indicators measuring customer value were chosen in coherence with the firm's strategy. Concerning the service related attributes, they are chosen according to (Daaboul et al., 2011). The value network with PDP 1 and CODP 1 stands for the AS-IS situation. After collecting all needed data, this scenario was simulated and validated by Alpina. After validating the AS-IS scenario (scenario 1), 8 other scenarios were identified. The scenario formed of PDP 2 and CODP 5 is not feasible since it is not possible to have such customization of the shoe and have the customer order decoupling point positioned at the last stage.

Results

The obtained results were validated by Alpina. The results were coherent with other analysis and studies realized at Alpina. As shown in Figure 5, the value for Alpina is higher in the case of PDP at position 1. This is due to the fact that Alpina was not ready to make the necessary changes and investments such as changing its agreements with its suppliers in order to offer MC shoes. Whereas the customer perceived value has a less predictable behavior. It is impacted by many indicators such as product perceived quality, offered customization, and order delay. The results (Figure 5) show that the value for Alpina is higher for PDP1 than PDP2 no matter the position of CODP. The value for Alpina is highest when the CODP is further in the chain (CODP5). In addition it reduces dramatically for CODP 1, 2, 3 & 4 compared to CODP5. This is due to:

 An increase in cost of raw material caused by the increase in the required stock for each material.  Reduced economies of scale as the size of the production batch in the case of MC is equal to 1.  An increase in the cost of usage of resources, since the execution time of an activity increased with increasing set-up time.  A slight increase in transportation costs caused by increased late orders that require faster delivery therefore more expensive.  A slight increase in turnover. As part of the MC, Alpina estimated a 5% increase in sales. However, this increase in sales (an increase of 37,559 Euros) was not sufficient to cover the additional costs.

The Pareto Front for Alpina's results was built. It shows the dominating combinations of CODP and PDP. (Mason-Jones and Towill, 1999) stressed that to maximize competitive opportunity CODP position should be considered jointly with information decoupling point (when information from forecast and market meet). The early usage of distorted information

Fig. 1 . 3 CODPFig. 2 .

 132 Fig. 1. Multiple PDPs

Fig. 3 .Fig. 4 .

 34 Fig. 3. : Methodology for evaluation of value network evolution

Fig. 5 .

 5 Fig. 5. Scenarios' performance and Pareto front

could improve capacity management. Hence, future work is needed in order to investigate the joint impact of CODP, PDP and the information decoupling point (IDP).