Plankton networks driving carbon export in the oligotrophic ocean

Abstract

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 uncharacterized. 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 and 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 a few bacterial and viral genes can predict a significant fraction of the variability in carbon export in these regions.

Extended Data Figure 1:

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 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 scale-free. 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.

Extended Data Figure 2:

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⁻¹⁵). b, A single prokaryotic subnetwork (n=109, N=1′527) is moderately associated to carbon export (Pearson cor. 0.32, p = 9e⁻⁰³). c, A single viral subnetwork (n=277, N=5′476) is strongly associated to carbon export (Pearson cor. 0.93, p = 2e⁻¹⁵). d,e,f, The WGCNA approach directly links subnetworks to environmental parameters, i.e. the more the features contribute to the subnetwork structure (topology), the more their abundance are correlated to the parameter. This measure allows to identify subnetworks for which the overall structure, summarized as the eigen vector of the subnetwork, is related to the carbon export. d, The eukaryotic subnetwork structure correlates to carbon export (Pearson cor. = 0.87, p = 5e⁻¹⁶). e, The prokaryotic subnetwork structure correlates to carbon export (Pearson cor. = 0.47, p = 5e⁻⁰⁶). f, The viral population subnetwork structure correlates to carbon export (Pearson cor. = 0.88, p = 6e⁻⁹³). g,h,i, Lineage subnetworks predict carbon export. PLS regression was used to predict carbon export using lineage abundances in selected subnetworks. LOOCV was performed and VIP scores computed for each lineage. g, The eukaryotic subnetwork predicts carbon export with a R² of 0.69. h, The prokaryotic subnetwork predicts carbon export with a R² of 0.60. i, The viral population subnetwork predicts carbon export with a R² of 0.89. j, k, l, Synechococcus (rather than Prochlorococcus) absolute cell counts correlate well to carbon export. j, Prochlorococcus cell counts estimated by flow cytometry do not correlate to carbon export (mean carbon flux at 150m, Pearson cor. = −0.13, p = 0.27). k, Synechococcus cell counts estimated by flow cytometry correlate significantly to carbon export (Pearson cor. = 0.64, p = 4.0e⁻¹⁰). l, Synechococcus / Prochlorococcus cell counts ratio correlates significantly to carbon export (Pearson cor. = 0.54, p = 4.0e⁻⁰⁷).

Extended Data Figure 3:

Prokaryotic function subnetworks associated to environmental parameters and their structure correlate to carbon export. a,b,c Global ecological networks were built for the prokaryotic functions using the WGCNA methodology (see methods) and correlated to classical oceanographic parameters as well as carbon export. a, Two bacterial functional subnetworks (n=441 and n=220 N=37′832) are associated to carbon export (Pearson cor. 0.54, p = 1e⁻⁰⁷ and 0.42, p = 1e⁻⁰⁴). b, The WGCNA approach directly links subnetworks to environmental parameters, i.e. the more the features contribute to the subnetwork structure (topology), the more their abundance are correlated to the parameter. This measure allows to identify subnetworks for which the overall structure, summarized as the eigen vector of the subnetwork, is related to the carbon export. The bacterial function subnetwork structures correlate to carbon export (FNET1 Pearson cor. = 0.68, p = 3e⁻⁶¹, and FNET2 Pearson cor. = 0.47, p = 6e⁻¹³). c, Two functional subnetworks (light and dark green, FNET1 (n=220) and FNET2 (n=441), respectively) are significantly associated with carbon export (FNET1: Pearson cor. 0.42, p = 4e⁻⁰⁹ and FNET2: 0.54, p = 7e⁻⁰⁶). The highest VIP score functions from top to bottom correspond to red dots from right to left. d, PLS regression was used to predict carbon export using abundances of functions (OGs) in selected subnetworks. LOOCV was performed and VIP scores computed for each function. Light green subnetwork (FNET1) functions predict carbon export with a R² of 0.41. Dark green subnetwork (FNET2) functions predict carbon export with a R² of 0.48. e, Cumulative abundance of genus-level taxonomic annotations of genes encoding functions from FNET1 and FNET2 subnetworks and Bacterial function subnetworks predict carbon export. Genes contributing to the relative abundance of FNET1 and FNET2 subnetwork functions were taxonomically annotated by homology searches against a non-redundant gene reference database using a last common ancestor (LCA) approach (see methods).

Extended Data Figure 4:

Viral protein cluster networks reveal potential marker genes for carbon export prediction at global scale. a, A viral protein cluster (PC) network was built using abundances of PCs predicted from viral population contigs associated to carbon export (Fig. 2c) using the WGCNA methodology (see methods) and correlated to classical oceanographic parameters. Two viral PC subnetworks (n=1′879 and n=2′147, N=4′678, light and dark orange, VNET1 and VNET2, left and right panel respectively) are strongly associated to carbon export (VNET1: Pearson cor. 0.75, p = 3e⁻⁰⁷ and VNET2: 0.91, p = 3e⁻¹⁴). b, The viral PC subnetwork structures correlate to carbon export (VNET1 Pearson cor. = 0.91, p < 1e⁻²⁰⁰, and VNET2 Pearson cor. = 0.96, p < 1e⁻²⁰⁰). c, Size of dots is proportional to the VIP score computed for the PLS regression. d, Viral PC subnetworks predict carbon export. PLS regression was used to predict carbon export using abundances of viral protein clusters (PCs) in selected subnetworks. LOOCV was performed and VIP scores computed for each PC. Light orange subnetwork (VNET1, left panel) PCs predict carbon export with a R² of 0.55. Dark orange subnetwork (VNET2, right panel) PCs predict carbon export with a R² of 0.89.

Extended Data Figure 5:

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⁻¹⁴). b, The eukaryotic subnetwork structure correlates to carbon export (Pearson cor. = 0.94, p = 4e⁻²⁷). c, The eukaryotic subnetwork predicts carbon export with a R² 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)).

Figure 1. Global view of carbon fluxes along the Tara Oceans circumnavigation route and associated eukaryotic lineages

a, Carbon flux in mg.m⁻².d⁻¹ and carbon export at 150 m estimated from particles size distribution and abundance measured with the Underwater Vision Profiler 5 (UVP5). Stations at which environmental data are available (Supplementary 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 Table 2) can predict carbon export with high accuracy (PLS regression, LOOCV, R²=0.69, see Extended data Fig. 2g). 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 three Rhizaria (Collodaria, Collozoum inerme and Sticholonche sp.), one copepod (Oithona sp.), one siphonophore (Lilyopsis), three Dinophyceae and one ciliate (Spirotontonia turbinata). b, The selected prokaryotic subnetwork (n=109, see Supplementary Table 3) can predict carbon export with good accuracy (PLS regression, LOOCV, R²=0.60, see Extended data Fig. 2h). c, The selected viral population subnetwork (n=277, see Supplementary Table 4) can predict carbon export with high accuracy (PLS regression, LOOCV, R²=0.89, see Extended data Fig. 2i). Two viral populations with a high VIP score (red dots) are predicted as Synechococcus phages (see Supplementary Table 4).

Figure 3. Integrated plankton community network built from eukaryotic, prokaryotic and viral subnetworks related to carbon export at 150 m

Major lineages were selected within the three subnetworks (VIP > 1) (Supplementary Tables 2, 3 and 4). Co-occurrences between all lineages of interest were extracted, if present, from a previously established global co-occurrence network (see methods). Only lineages discussed within the study are pinpointed. The resulting graph is composed of 329 nodes, 467 edges, with a diameter of 7, and average weighted degree of 4.6.

Figure 4. Key bacterial functional categories associated with carbon export at 150 m at global scale

A bacterial functional network was built based on Orthologous Group/Gene (OG) relative abundances using the WGCNA methodology (see Methods) and correlated to classical oceanographic parameters. Two functional subnetworks (FNET1 (n=220) and FNET2 (n=441), respectively, Extended data Fig. 3a) are significantly associated with carbon export (FNET1: Pearson cor. 0.42, p = 4e⁻⁰⁹ and FNET2: 0.54, p = 7e⁻⁰⁶, see Extended data Fig. 3b). Higher functional categories are depicted for functions with a VIP score >1 (PLS regression, LOOCV, FNET1 R²=0.41 and FNET2 R²=0.48, see Extended data Fig. 3d) in both subnetworks.