{'fraction_non_alphanumeric': 0.032044376189057235, 'fraction_numerical': 0.020308524537557018, 'mean_word_length': 5.161276655839571, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Plankton networks driving carbon export in the oligotrophic ocean

The biological carbon pump is the process by which CO2 is transformed to organic carbon via photosynthesis, exported through sinking particles, and finally sequestered in the deep ocean. While the intensity of the pump correlates with plankton community composition, the underlying ecosystem structure driving the process remains largely uncharacterised. Here we use environmental and metagenomic data gathered during the Tara Oceans expedition to improve our understanding of carbon export in the oligotrophic ocean. We show that specific plankton communities, from the surface and deep chlorophyll maximum, correlate with carbon export at 150 m and highlight unexpected taxa such as Radiolaria, alveolate parasites, as well as Synechococcus and their phages, as lineages most strongly associated with carbon export in the subtropical, nutrient-depleted, oligotrophic ocean. Additionally, we show that the relative abundance of just a few bacterial and viral genes can predict most of the variability in carbon export in these regions.

Marine planktonic photosynthetic organisms are responsible for approximately fifty percent of Earth's primary production and fuel the global ocean biological carbon pump 1 . The intensity of the pump is correlated to plankton community composition 2,3 , and controlled by the relative rates of primary production and carbon remineralisation 4 . About 10% of this newly produced organic carbon in the surface ocean is exported through gravitational sinking of particles. Finally, after multiple transformations, a fraction of the exported material reaches the deep ocean where it is sequestered over thousand-year timescales 5 .
Like most biological systems, marine ecosystems in the sunlit upper layer of the ocean (denoted the euphotic zone) are complex 6,7 , characterised by a wide range of biotic and abiotic interactions [8][9][10] and in constant balance between carbon production, transfer to higher trophic levels, remineralisation, and export to the deep layers 11 . The marine ecosystem structure and its taxonomic and functional composition likely evolved to comply with this loss of energy by modifying organism turnover times and by the establishment of complex feedbacks between them 6 and the substrates they can exploit for metabolism 12 . Decades of groundbreaking research have focused on identifying independently the key players involved in the biological carbon pump. Among autotrophs, diatoms are commonly attributed to being important in carbon flux because of their large size and fast sinking rates [13][14][15] while small autotrophic picoplankton may contribute directly through subduction of surface water 16 or indirectly by aggregating with larger settling particles or consumption by organisms at higher trophic levels 17 . Among heterotrophs, zooplankton such as crustaceans impact carbon flux via production of fast-sinking fecal pellets while migrating hundreds of meters in the water-column 18,19 . These observations, focusing on just a few components of the marine ecosystem, highlight that carbon export results from multiple biotic interactions and that a better understanding of the mechanisms involved in its regulation will require an analysis of the entire planktonic ecosystem.
Advanced sequencing technologies offer the opportunity to simultaneously survey whole planktonic communities and associated molecular functions in unprecedented detail. Such a holistic approach may allow the identification of community-or gene-based biomarkers that could be used to monitor and predict ecosystem functions, e.g., related to the biogeochemistry of the ocean [20][21][22] . Here, we leverage global-scale ocean genomics datasets from the euphotic zone 10,[23][24][25] and associated environmental data to assess the coupling between ecosystem structure, functional repertoire, and carbon export at 150 m.
should keep in mind that fluxes are calculated from images of particles. These estimates are derived from an approximation of Stokes' law relating the equivalent spherical diameter of particles to carbon flux (see Methods). This exponential approximation is reasonable assuming similar particle composition across all sizes, as highlighted by the standard deviations of parameters in Eq. 5 (see Methods). Furthermore, because of instrument and method limitations, particles <250 μm were not used, which may underestimate total carbon fluxes. Finally, these fluxes are instantaneous because they do not integrate space and time as sediment traps would. However, the approach allowed us to assemble the largest homogeneous carbon export dataset during a single expedition, corresponding to more than 600 profiles over 150 stations. This dataset is of similar magnitude to the body of historical data available in the literature that includes the 134 deep sediment trap-based carbon flux time-series 27 from the JGOFS program and the 419 thorium-derived particulate organic carbon (POC) export measurements 28 .
From 68 globally distributed sites, a total of 7.2 Tb of metagenomics data, representing ~40 million non-redundant genes, around 35,000 Operational Taxonomic Units (OTUs) of prokaryotes (Bacteria and Archaea) and numerous mainly uncharacterized viruses and picoeukaryotes, have been described recently 23,25 . In addition, a set of 2.3 million eukaryotic 18S rDNA ribotypes was generated from a subset of 47 sampling sites corresponding to approximately 130,000 OTUs 24 (Fig. 1a). Finally, 5,476 viral "populations" were identified at 43 sites from viral metagenomic contigs, only 39 (<0.1%) of which had been previously observed 25 (see Methods). These genomics data combined across all domains of life and viruses together with carbon export estimates and other environmental parameters were used to explore the relationships between marine biogeochemistry and euphotic plankton communities (see Methods) in the top 150 m of the oligotrophic open ocean. Our study did not include high latitude areas due to the current lack of available molecular data and results should not be extrapolated to deeper depths.
Using a method for regression-based modeling of high multidimensional data in biology (specifically a sparse Partial Least Square analysis -sPLS 29 , Extended data Fig. 1), we detected several plankton lineages for which relative sequence abundance correlated with carbon export and other environmental parameters, most notably with NPP, as expected (Fig.  1b and see Supplementary Table 1). These included diatoms, dinoflagellates and metazoa (zooplankton), lineages classically identified as key contributors to carbon export. most associated with carbon export rather than predicting organisms associated with sinking particles.
In brief, the WGCNA approach builds a network in which nodes are features (in this case plankton lineages or gene functions) and links are evaluated by the robustness of cooccurrence scores. WGCNA then clusters the network into modules (hereafter denoted subnetworks) that can be examined to find significant subnetwork-trait relationships. We then filtered each subnetwork using a Partial Least Square (PLS) analysis that emphasizes key nodes (based on the Variable Importance in Projection (VIP) scores; see Methods and Extended data Fig. 1). These particular nodes are mandatory to summarize a subnetwork (or community) related to carbon export. In particular, they are of interest for evaluating (i) subnetwork robustness and (ii) predictive power for a given trait (see Methods and Extended data Fig. 1).
We applied WGCNA to the relative abundance tables of eukaryotic, prokaryotic and viral lineages [23][24][25] and identified unique subnetworks significantly associated with carbon export within each dataset (see Methods and Supplementary Tables 2, 3 Table 2). Surprisingly, this small subnetwork's structure correlates very strongly to carbon export (Pearson cor. = 0.87, p = 5e −16 , Extended data Fig. 2d) and it predicts as much as 69% (Leave-One-Out Cross-Validated (LOOCV), R 2 = 0.69) of the variability in carbon export (Extended data Fig. 2g). Only ~6% of the subnetwork nodes correspond to diatoms and they show lower VIP scores than dinoflagellates (Supplementary Table 2). This is likely because our samples are not from silicate replete conditions where diatoms were blooming. Furthermore, our analysis did not incorporate data from high latitudes, where diatoms are known to be particularly important for carbon export, so this result suggests that dinoflagellates have a heretofore unrecognized role in carbon export processes in subtropical oligotrophic 'type' ecosystems. More precisely four of the five highest VIP scoring eukaryotic lineages that correlated with carbon export at 150 m were heterotrophs such as Metazoa (copepods), non-photosynthetic Dinophyceae, and Rhizaria (Fig. 2a and  Supplementary Table 2). These results corroborate recent metagenomics analysis of microbial communities from sediment traps in the oligotrophic North Pacific subtropical gyre 34 . Consistently, in situ imaging surveys have revealed Rhizarian lineages, made up of large fragile organisms such as the Collodaria, to represent an until now under-appreciated component of global plankton biomass 35 , which here also appear to be of relevance for carbon export. Another 14% of lineages from the subnetwork correspond to parasitic organisms, a largely under-explored component of planktonic ecosystems when studying carbon export.
The prokaryotic subnetwork that associated most significantly with carbon export at 150 m (subnetwork-trait relationship to carbon export, Pearson cor. = 0.32, p = 9e −03 ) contained 109 OTUs (Extended data Fig. 2b and Supplementary Table 3), its structure correlated well to carbon export (Pearson cor. = 0.47, p = 5e −06 , Extended data Fig. 2e) and it could predict as much as 60% of the carbon export (LOOCV, R 2 = 0.60) (Extended data Fig. 2h). By far the highest VIP score within this community was assigned to Synechococcus, followed by Cobetia, Pseudoalteromonas and Idiomarina, as well as Vibrio and Arcobacter (Fig. 2b and  Supplementary Table 3). Noteworthy, Prochlorococcus genera and SAR11 clade fall out of this community, while the significance of Synechococcus for carbon export could be validated using absolute cell counts estimated by flow cytometry (Pearson cor. = 0.64, p = 4e −10 , Extended data Fig. 2k). Moreover, Prochlorococcus cell counts did not correlate with carbon export (Pearson cor. = −0.13, p = 0.27, Extended data Fig. 2j) whereas the Synechococcus to Prochlorococcus cell count ratio correlated positively and significantly (Pearson cor. = 0.54, p = 4e −07 , Extended data Fig. 2l), suggesting the relevance of Synechococcus, rather than Prochlorococcus, to carbon export. Interestingly, Pseudoalteromonas, Idiomarina, Vibrio and Arcobacter (of which several species are known to be associated with eukaryotes 36 ) have also been observed in live and poisoned sediment traps 34 and display very high VIP scores in the subnetwork associated with carbon export. Additional genera reported as being enriched in poisoned traps (also known as being associated with eukaryotes) include Enterovibrio and Campylobacter, and are present as well in the carbon export associated subnetwork.  Table 4), which represented a 16-fold enrichment (Fisher's exact test p = 6.4e −09 ). Its structure also correlated with carbon export (Pearson cor. = 0.88, p = 6e −93 , Extended data Fig. 2f) and could predict up to 89% of the variability of carbon export (LOOCV, R 2 = 0.89) (Extended data Fig. 2i). The significance of these convergent results is reinforced by the fact that sequences from these datasets are derived from organisms collected on independent size filters (see Methods), and further implicates the importance of top-down processes in carbon export.
With the aim of integrating eukaryotic, prokaryotic, and viral communities in the euphotic zone with carbon export at 150 m, we synthesized their respective subnetworks using a single global co-occurrence network established previously 10 . The resulting network focused on key lineages and their predicted co-occurrences (Fig. 3). Lineages with high VIP values (such as Synechococcus) are revealed as hubs of the co-occurrence network 10 , illustrating the potentially strategic key roles within the integrated network of lineages underappreciated by conventional methods to study carbon export. Associations between the hub lineages are mostly mutually exclusive which may explain the relatively weak correlation of some of these lineages with carbon export when using standard correlation analyses as shown in Fig. 1b.

Gene functions associated with carbon export
Given the potential importance of prokaryotic processes influencing the biological carbon pump 22 , we used the same analytical approaches to examine the prokaryotic genomic functions associated with carbon export at 150 m in the annotated Ocean Microbial Reference Gene Catalogue from Tara Oceans 23 . We built a global co-occurrence network for functions (i.e., Orthologous Groups of genes or OGs) from the euphotic zone and identified two subnetworks of functions that are significantly associated with carbon export (light and dark green subnetworks; FNET1 and FNET2, respectively, see Extended data Fig. 3a, 3b and 3c).
The majority of functions in FNET1 and FNET2 correlate well with carbon export (FNET1: mean Pearson cor. = 0.45, s.d. 0.09 and FNET2: mean Pearson cor. = 0.34, s.d. 0.10). Interestingly, FNET2 functions (n=220) encode mostly (83%) core functions (i.e., functions observed in all euphotic samples, see Methods) while the majority of FNET1 functions (n=441) are non-core (85%) (see Supplementary Tables 5, 6), highlighting both essential and adaptive ecological functions associated with carbon export. Top VIP scoring functions in the FNET1 subnetwork are membrane proteins such as ABC-type sugar transporters (Extended data Fig. 3c). This subnetwork also contains many functions specific to the Synechococcus accessory photosynthetic apparatus (e.g., relating to phycobilisomes, phycocyanin and phycoerythrin; see Supplementary Table 5), which is consistent with the major role of this genus for carbon export inferred from the prokaryotic subnetwork (Fig.  2b). In addition, functions related to carbohydrates, inorganic ion transport and metabolism, as well as transcription, are also well represented (Fig. 4), suggesting overall a subnetwork of functions dedicated to photosynthesis and growth.
The FNET2 subnetwork contains several functions encoded by genes taxonomically assigned to Candidatus pelagibacter and Prochlorococcus, known as occupying similar oceanic regions as Synechococcus, but overall most of its relative abundance (74%) is taxonomically unclassified (Extended data Fig. 3e). Top VIP scoring functions in FNET2 are also membrane proteins and ABC-type sugar transporters, as well as functions involved in carbohydrate breakdown such as a chitinase (Extended data Fig. 3c). These features highlight the potential roles of bacteria in the formation and degradation of marine aggregates 37 . Strikingly, 77% and 58%, of OGs with a VIP score > 1 in FNET1 and FNET2, respectively, are functionally uncharacterized 38,39 (Fig. 4), pointing to the strong need for future molecular work to explore these functions (see Supplementary Tables 5, 6).
The relevance of the identified bacterial functions to predict carbon export was also confirmed by PLS regression (Extended data Fig. 3d). As proposed for plankton communities, the functional subnetworks predict 41% and 48% of carbon export variability (LOOCV, R 2 = 0.41 and 0.48 for FNET1 and FNET2, respectively) with a minimal number of functions (Fig. 4, 123 and 54 functions with a VIP score > 1 for FNET1 and FNET2, respectively). Finally, higher predictive power was obtained using subnetworks of viral protein clusters (Extended data Fig. 4a, 4b and 4c), predicting 55% and 89% of carbon export variability (LOOCV R 2 = 0.55 and 0.89 for VNET1 and VNET2, respectively; Extended data Fig. 4d, Supplementary Tables 7, 8), suggesting the key role, of not only bacteria, but also their phages in processes sustaining carbon export at a global level.

Discussion
In this report we reveal the potential contribution of unexpected components of plankton communities, and confirm the importance of prokaryotes and viruses in the correlating with carbon export in the nutrient-depleted oligotrophic ocean. Carbon export at 150 m has been estimated from particle size distribution in a global dataset, but should be taken with caution, as the estimates do not account for particle composition. In addition, these export estimates evaluate how much carbon leaves the euphotic zone, but they are not related and should not be extrapolated to sequestration, which occurs after remineralisation, deeper in the water column, and over longer timescales. Nonetheless, the use of the UVP was the only realistic method to evaluate carbon flux over the 3 years expedition because deployment of sediment traps at all stations would have been impossible. While our findings are consistent with the numerous previous studies that have highlighted the central role of copepods and diatoms in carbon export 14,15,[17][18][19] , they place them in an ecosystem context and reveal hypothetical processes correlating with the intensity of export, such as parasitism and predation. For example, while viruses are commonly assumed to lyse cells and maintain fixed organic carbon in surface waters, thereby reducing the intensity of the biological carbon pump 40 , there are hints that viral lysis may increase carbon export through the production of colloidal particles and aggregate formation 41 . Our current study suggests that these latter roles may be more ubiquitous than currently appreciated. The importance of aggregation and cell stickiness as inferred from gene network analysis should be further explored mechanistically to investigate the biological significance of these findings.
The future evolution of the oceanic carbon sink remains uncertain because of poorly constrained processes, particularly those associated with the biological pump. With current trends in climate change, the size and biodiversity of phytoplankton are predicted to decrease globally 42,43 . Furthermore, in spite of the potential importance of viruses revealed in this study, they have largely been ignored because of limitations in sampling technologies. Consequently, as oligotrophic gyres expand and global mean NPP decreases 44 , the field is currently unable to predict the consequences for carbon export from the ocean's euphotic zone. By pinpointing key lineages and key microbial functions that correlate with carbon export at 150 m in these areas, this study provides a framework to address this critical bottleneck. However, the associations presented do not necessarily suggest a causal effect on carbon export, which will require further investigation.
One of the grand challenges in the life sciences is to link genes to ecosystems 45 , based on the posit that genes can have predictable ecological footprints at community and ecosystem levels [46][47][48] . The Tara Oceans data sets have allowed us to predict as much as 89% of the variability in carbon export from the oligotrophic surface ocean with just a small number of genes, largely with unknown functions, encoded by prokaryotes and viruses. These findings can be used as a basis to include biological complexity and guide experimental work designed to inform climate modeling of the global carbon cycle. Such statistical analyses, scaling from gene-to-ecosystems, may open the way to the development of a new conceptual and methodological framework to better understand the mechanisms underpinning key ecological processes.

Environmental data collection
From 2009-2013, environmental data (Supplementary Table 9) were collected across all major oceanic provinces in the context of the Tara Oceans expeditions 20 . Sampling stations were selected to represent distinct marine ecosystems at a global scale 49 . Note that Southern Ocean stations were not examined herein because they were ranked as outliers due to their exceptional environmental characteristics and biota 23,24 . Environmental data were obtained from vertical profiles of a sampling package 50,51 . It consisted of conductivity and temperature sensors, chlorophyll and CDOM fluorometers, light transmissometer (Wetlabs C-star 25cm), a backscatter sensor (WetLabs ECO BB), a nitrate sensor (SATLANTIC ISUS) and a Hydroptic Underwater Vision Profiler (UVP; Hydroptics 52 . Nitrate and fluorescence to chlorophyll concentrations as well as salinity were calibrated from water samples collected with Niskin bottle 50 . Net Primary Production (NPP) data were extracted from 8 day composites of the Vertically Generalized Production Model (VGPM 53 ) at the week of sampling 54 . Carbon fluxes and carbon export, corresponding to the carbon flux at 150 m, were estimated based on particle concentration and size distributions obtained from the UVP 51 and details are presented below.

From particle size distribution to carbon export estimation
Previous research has shown that the distribution of particle size follows a power law over the μm to the mm size range 3,55,56 . This Junge-type distribution translates into the following mathematical equation, whose parameters can be retrieved from UVP images: (eq. 1) where d is the particle diameter, and exponent k is defined as the slope of the number spectrum when equation (2) is log transformed. This slope is commonly used as a descriptor of the shape of the aggregate size distribution. The carbon-based particle size approach relies on the assumption that the total carbon flux of particles (F) corresponds to the flux spectrum integrated over all particle sizes: where n(d) is the particle size spectrum, i.e., equation (1), and m(d) is the mass (here carbon content) of a spherical particle described as: where α = πρ/6, ρ is the average density of the particle, and w(d) is the settling rate calculated using Stokes Law: (eq. 4) where β =g(ρ -ρ 0 )(18νρ 0 ) −1 , g is the gravitational acceleration, ρ 0 the fluid density, and ν the kinematic viscosity.
In addition, mass and settling rates of particles, m(d) and w(d), respectively, are often described as power law functions of their diameter obtained by fitting observed data, m(d). w(d) = Ad B . The particles carbon flux can then be estimated using an approximation of Eq. 2 over a finite number (x) of small logarithmic intervals for diameter d spanning from 250 μm to 1.5 mm (particles <250 μm and >1.5 mm are not considered, consistent with the method presented by Guidi et al., [2008] 57 ) such as (eq. 5) where A=12.5±3.40 and B=3.81 ± 0.70 have been estimated using a global dataset that compared particle fluxes in sediment traps and particle size distributions from the UVP images.

Genomic data collection
For the sake of consistency between all available datasets from the Tara Oceans expeditions, we considered subsets of the data recently published in Science [23][24][25] . In brief, one sample corresponds to data collected at one depth (surface (SRF) or Deep Cholorophyll Maximum (DCM) determined from the profile of chlorophyll fluorometer) and at one station. To study the eukaryotic community in our current manuscript, we selected stations at which we had environmental data and carbon export estimated at 150 m with the UVP and all size fractions. Consequently a subset of 33 stations (corresponding to 56 Supplementary Table 10). In addition a detailed table is provided summarizing which samples (depth and station) are available for each domain (Supplementary Table 11).

Eukaryotic taxa profiling
Photic-zone eukaryotic plankton diversity has been investigated through millions of environmental Illumina reads. Sequences of the 18S ribosomal RNA gene V9 region were obtained by PCR amplification and a stringent quality-check pipeline has been applied to remove potential chimera or rare sequences (details on data cleaning in de Vargas et al. [2015] 24 ). For 47 stations, and if possible at two depths (SRF and DCM), eukaryotic communities were sampled in the piconano-(0.8-5 μm), micro-(20-180 μm) and mesoplankton (180-2000 μm) fractions (a detailed list of these samples is given in Supplementary  Table 12). In the framework of the carbon export study, sequences from all size fractions were pooled in order to get the most accurate and statistically reliable dataset of the eukaryotic community. The 2.3 million eukaryotic ribotypes were assigned to known eukaryotic taxonomic entities by global alignment to a curated database 24 . To get the most accurate vision of the eukaryotic community, sequences showing less than 97% identity with reference sequences were excluded. The final eukaryotic relative abundance matrix used in our analyses included 1,750 lineages (taxonomic assignation has been performed using a last common ancestor methodology, and had thus been performed down to species level when possible) in 56 samples from 33 stations. Pooled abundance (number of V9 sequences) of each lineage has been normalized by the total sum of sequences in each sample.

Prokaryotic taxa profiling
To investigate the prokaryotic lineages, communities were sampled in the pico-plankton. Both filter sizes have been used along the Tara Oceans transect: up to station #52, prokaryotic fractions correspond to a 0.22-1.6 μm size fraction, and from station #56, prokaryotic fractions correspond to a 0.22-3 μm size fraction. Prokaryotic taxonomic profiling was performed using 16S rRNA gene tags directly identified in Illumina-sequenced metagenomes ( mi tags) as described in Logares 23 ). The photic-zone prokaryotic relative abundance matrix used in our analyses included 3,253,962 mi tags corresponding to 1,328 genera in 104 samples from 62 stations.

Prokaryotic functional profiling
For each prokaryotic sample, gene relative abundance profiles were generated by mapping reads to the OM-RGC using the MOCAT pipeline 61 . The relative abundance of each reference gene was calculated as gene length-normalized base counts. And functional abundances were calculated as the sum of the relative abundances of these reference genes, annotated to OG functional groups. In our analyses, we used the subset of the OM-RGC that was annotated to Bacteria or Archaea (24.4 M genes). Using a rarefied (to 33 M inserts) gene count table, an OG was considered to be part of the ocean microbial core if at least one insert from each sample was mapped to a gene annotated to that OG. For further details on the prokaryotic profiling please refer to Sunagawa et al. [2015] 23 . The final prokaryotic functional relative abundance matrix used in our analyses included 37,832 OGs or functions in 104 samples from 62 stations. Genes from functions of FNET1 and FNET2 subnetworks were taxonomically annotated using a modified dual BLAST-based last common ancestor (2bLCA) approach 62 . We used RAPsearch2 63 rather than BLAST to efficiently process the large data volume and a database of non-redundant protein sequences from UniProt (version: UniRef_2013_07) and eukaryotic transcriptome data not represented in UniRef (see Supplementary Table 5, 6, for full annotations).

Enumeration of prokaryotes by flow cytometry
For prokaryote enumeration by flow cytometry, three aliquots of 1 ml of seawater (prefiltered by 200-μm mesh) were collected from both SRF and DCM. The samples were fixed immediately using cold 25% glutaraldehyde (final concentration 0.125%), left in the dark for 10 min at room temperature, flash-frozen and kept in liquid nitrogen on board and then stored at −80°C on land. Two subsamples were taken to separate counts of heterotrophic prokaryotes (not shown herein) and phototrophic picoplankton. For heterotrophic prokaryote determination, 400 μl of sample was added to a diluted SYTO-13 (Molecular Probes Inc., Eugene, OR, USA) stock (10:1) at 2.5 μmol l −1 final concentration, left for about 10 min in the dark to complete the staining and run in the flow cytometer. We used a FacsCalibur (Becton & Dickinson) flow cytometer equipped with a 15 mW Argon-ion laser (488 nm emission). At least 30,000 events were acquired for each subsample (usually 100,000 events). Fluorescent beads (1 μm, Fluoresbrite carboxylate microspheres, Polysciences Inc., Warrington, PA) were added at a known density as internal standards. The bead standard concentration was determined by epifluorescence microscopy. For phototrophic picoplankton, we used the same procedure as for heterotrophic prokaryote, but without addition of SYTO-13. Data analysis was performed with FlowJo software (Tree Star, Inc.).

Profiling of viral populations
In order to associate viruses to carbon export we used viral populations as defined in Brum et al.
[2015] 25 using a set of 43 Tara Oceans viromes. Briefly, viral populations were defined as large contigs (>10 predicted genes and >10 kb) identified as most likely originating from bacterial or archaeal viruses. These 6,322 contigs remained and were then clustered into populations if they shared more than 80% of their genes at >95% nucleotide identity. This resulted in 5,477 'populations' from the 6,322 contigs, where as many as 12 contigs were included per population. For each population, the longest contig was chosen as the 'seed' representative sequence. The relative abundance of each population was computed by mapping all quality-controlled reads to the set of 5,477 non-redundant populations (considering only mapping quality scores greater than 1) with Bowtie2 64 and if more than 75% of the reference sequence was covered by virome reads. The relative abundance of a population in a sample was computed as the number of base pairs recruited to the contig normalized to the total number of base pairs available in the virome and the contig length if more than 75% of the reference sequence was covered by virome reads, and set to 0 otherwise (see Brum et al. [2015] 25

Viral host predictions
The longest contig in a population was defined as the seed sequence and considered the best estimate of that population's origin. These seed sequences were used to assess taxonomic affiliation of each viral population. Cases where >50% of the genes were affiliated to a specific reference genome from RefSeq Virus (based on a BLASTp comparison with thresholds of 50 for bit score and 10 −5 for e-value) with an identity percentage of at least 75% (at the protein sequence level) were considered as confident affiliations to the corresponding reference virus. The viral population host group was then estimated based on these confident affiliations (see Supplementary Table 13 for host affiliation of viral population contigs associated to carbon export).

Viral protein clusters
Viral protein clusters (PCs) correspond to ORFs initially mapped to existing clusters (POV, GOS and phage genomes). The remaining, unmapped ORFs were self-clustered, using cd-hit as described in Brum et al. [2015] 25 . Only PCs with more than two ORFs were considered bona fide and were used for subsequent analyses. To compute PC relative abundance for statistical analyses, reads were mapped back to predicted ORFs in the contigs dataset using Mosaik as described in Brum et al. [2015] 25 . Read counts to PCs were normalized by sequencing depth of each virome. Importantly, we restricted our analyses to 4,294 PCs associated to the 277 viral population contigs significantly associated to carbon export in 37 samples from 22 stations.

Sparse Partial Least Squares analysis
In order to directly associate eukaryotic lineages to carbon export and other environmental traits (Fig. 1b), we used sparse Partial Least Square (sPLS 65 as implemented in the R package mixOmics 29 . We applied the sPLS in regression mode, which will model a causal relationship between the lineages and the environmental traits, i.e. PLS will predict environmental traits (e.g. carbon export) from lineage abundances. This approach enabled us to identify high correlations (see Supplementary Table 1) between certain lineages and carbon export but without taking into account the global structure of the planktonic community.

Co-occurrence network model analysis
Weighted correlation network analysis (WGCNA) was performed to delineate feature (lineages, viral populations, PCs or functions) subnetworks based on their relative abundance 66,67 . A signed adjacency measure for each pair of features was calculated by raising the absolute value of their Pearson correlation coefficient to the power of a parameter p. The default value p=6 was used for each global network, except for the Prokaryotic functional network where p had to be lowered to 4 in order to optimize the scale-free topology network fit. Indeed, this power allows the weighted correlation network to show a scale free topology where key nodes are highly connected with others. The obtained adjacency matrix was then used to calculate the topological overlap measure (TOM), which for each pair of features, taking into account their weighted pairwise correlation (direct relationships) and their weighted correlations with other features in the network (indirect relationships). For identifying subnetworks a hierarchical clustering was performed using a distance based on the TOM measure. This resulted in the definition of several subnetworks, each represented by its first principal component.
These characteristic components play a key role in weighted correlation network analysis. On the one hand, the closeness of each feature to its cluster, referred to as the subnetwork membership, is measured by correlating its relative abundance with the first principal component of the subnetwork. On the other hand, association between the subnetworks and a given trait is measured by the pairwise Pearson correlation coefficients between the considered environmental trait and their respective principal components. A similar protocol has been performed on the eukaryotic relative abundance matrix, the prokaryotic relative abundance matrix, the prokaryotic functions relative abundance matrix and the viral population and PC relative abundance matrices. All procedures were applied on Hellingertransformed log-scaled abundances. Noteworthy, the protocol is not sensitive to copy number variation as observed across different eukaryotic species, because the association between two species relies on a correlation score between relative abundance measurements. Computations were carried out using the R package WGCNA 33 .
Given the nature of the eukaryotic dataset (three distinct size fractions), the sampling process may lead to the loss of size fractions. In particular, samples #1, #3, #17, #37, #39,  #43, #48, #53, #54, #55, #66 are eventually biases by such a loss (Supplementary Table 12). A complementary WGCNA analysis was performed with addition of these samples to evaluate the robustness of our protocol to missing size fractions. The composition of the eukaryotic subnetwork built with an extended dataset (i.e., 67 samples from 37 stations for which size fractions were missing in 11 samples) was compared to the subnetwork as presented above (i.e., 56 samples from 33 stations). Both subnetworks shown an overlap of 75% of lineage, whereas four of the top five VIP lineages with the extended dataset (see Extended data Fig. 5 for details) can be found in the top six VIP lineages of the above subnetwork (Supplementary Table 2), emphasizing highly similar results and a small sensitivity to size fraction loss.

Extraction of subnetworks related to carbon export
For each subnetwork (called modules within WGCNA) extracted from each global network, pairwise Pearson correlation coefficients between the subnetwork principal components and the carbon export estimation was computed, as well as corresponding p-values corrected for multiple testing using the Benjamini & Hochberg FDR procedure. The subnetworks showing the highest correlation scores are of interest and were investigated. One subnetwork (49 nodes) was significant within the eukaryotic network; one subnetwork (109 nodes) was significant for the prokaryotic network; one subnetwork (277 nodes) was significant within the virus network; two subnetworks (441 and 220 nodes) were significant within the prokaryotic functional network, and two subnetworks (1,879 and 2,147 nodes) were significant within the viral PCs network.

Partial Least Squares regression
In addition to the network analyses, we asked whether the identified subnetworks can be used as predictors for the carbon export estimations. To answer this question, we used Partial least squares (PLS) regression, which is a dimensionality-reduction method that aims at determining predictor combinations with maximum covariance with the response variable. The identified combinations, called latent variables, are used to predict the response variable. The predictive power of the model is assessed by correlating the predicted vector with the measured values. The significance of the prediction power was evaluated by permuting the data 10,000 times. For each permutation, a PLS model was built to predict the randomized response variable and a Pearson correlation was calculated between the permuted response variable and in Leave-One-Out Cross-Validation (LOOCV) predicted values. The 10,000 random correlations are compared to the performance of the PLS model that were used to predict the true response variable. In addition, the predictors were ranked according to their value importance in projection (VIP) 68 . The VIP measure of a predictor estimates its contribution in the PLS regression. The predictors having high VIP values are assumed important for the PLS prediction of the response variable. The VIP values of the prokaryotic functional subnetworks are provided in Supplementary Tables 5, 6. For the sake of illustration, only lineages or functions with VIP > 1 68 are discussed and pictured in Figure  2 and 4. Our computations were carried out using the R package pls 69 . All programs are available under GPL Licence.

Subnetwork representations
Nodes of the subnetworks represent either lineages (eukaryotic, prokaryotic or viral) or functions (prokaryotic or viral). Subnetworks related to the carbon export have been represented in two distinct formats. Scatter plots represent each nodes based on their Pearson correlation to the carbon export and their respective node centrality within the subnetwork. The latter has been recomputed using significant Spearman correlations above 0.3 (>0.9 for viral PCs) as edges, this is done for visualization purposes since WGCNA subnetworks (based on the Topology Overlap Measure (TOM) between nodes) are hyperconnected. Size representation of nodes are proportional to the VIP score after PLS. The hiveplots depict the same subnetworks by focusing on two main features: x-axis and y-axis depict nodes of subnetworks ranked by their VIP scores and Pearson correlation to the carbon export, respectively. Guidi  Overview of analytical methods used in the manuscript. a, Depiction of a standard pairwise analysis that considers a sequence relative abundance matrix for s samples (s × OTUs (Operational Taxonomic Units)) and its corresponding environmental matrix (s × p (parameters)). sPLS results emphasize OTU(s) that are the most correlated to environmental parameters. b, Depiction of a graph-based approach. Using only a relative abundance matrix (s × OTUs), WGCNA builds a graph where nodes are OTUs and edges represent significant Europe PMC Funders Author Manuscripts co-occurrence. Co-occurrence scores between nodes are weights allocated to corresponding edges. These weights are magnified by a power-law function until the graph becomes scalefree. The graph is then decomposed within subnetworks (groups of OTUs) that are analyzed separately. One subnetwork (group of OTUs) is considered of interest when its topology is related to the trait of interest; in the current case carbon export. For each subnetwork (for instance the subnetwork related to carbon export), each OTU is spread within a feature space that plots each OTU based on its membership to the subnetwork (x-axis) and its correlation to the environmental trait of interest (i.e., carbon export). A good regression of all OTUs emphasizes the putative relation of the subnetwork topology and the carbon export trait (i.e. the more a given OTU defines the subnetwork topology, the more it is correlated to carbon export). c, Depiction of the machine learning (PLS) approach that was applied following subnetwork identification and selection. Greater VIP scores (i.e. larger circles) emphasized most important OTUs. VIP refers to Variable Importance in Projection and reflects the relative predictive power of a given OTU. OTUs with VIP score greater than 1 are considered as important in the predictive model and their selection do not alter the overall predictive power. Lineage ecological subnetworks associated to environmental parameters and their structures correlating to carbon export. a,b,c, Global ecological networks were built using the WGCNA methodology (see methods) and correlated to classical oceanographic parameters as well as carbon export (estimated at 150 m from particles size distribution and abundance). Each domain-specific global network is decomposed into smaller coherent subnetworks (depicted by distinct colours on the y-axis) and their eigen vector is correlated to all environmental parameters. Similar to a correlation at the network scale, this approach directly links subnetworks to environmental parameters (i.e. the more the taxa contribute to the subnetwork structure, the more their abundance are correlated to the parameter). a, A single eukaryotic subnetwork (n=58, N=1′870) is strongly associated to carbon export (Pearson cor. 0.81, p = 5e −15  WGCNA and PLS regression analyses for the full Eukaryotic dataset. a, A single eukaryotic subnetwork (n=58, is strongly associated to carbon export (Pearson cor. 0.79, p = 3e −14 ). b, The eukaryotic subnetwork structure correlates to carbon export (Pearson cor. = 0.94, p = 4e −27 ). c, The eukaryotic subnetwork predicts carbon export with a R 2 of 0.76. d, Lineages with the highest VIP score (dots size is proportional to the VIP score in the scatter plot) in the PLS are depicted as red dots corresponding to two rhizaria (Collodaria), one copepod (Euchaeta), and three dinophyceae (Noctiluca scintillans, Gonyaulax polygramma and Gonyaulax sp. (clade 4)).  Table 9) are depicted by white dots. Stations at which eukaryotic samples are available are colored in red (Supplementary Tables 10 and 12). b, Eukaryotic lineages associated to carbon export as revealed by standard methods for regression-based modeling (sPLS analysis). Correlations between lineages and environmental parameters are depicted as a clustered heatmap and lineages with a correlation to carbon export higher than 0.2 are highlighted (detailed results in Supplementary Table 1).

Figure 2. Ecological networks reveal key lineages associated with carbon export at 150 m at global scale
The relative abundances of taxa in selected subnetworks were used to estimate carbon export and to identify key lineages associated with the process. a, The selected eukaryotic subnetwork (n=49, see Supplementary