dataesr/scientific-paragraphs-categorization huggingface dataset:
6. Finally, investigators are encouraged to report all N 2 fixation rate estimates, clearly denoting those below the method LODs, as the latter provide valuable information about environments in which N 2 fixation is not detected. We have created templates for calculations of nitrogen fixation rates and associated LODs (see "spreadsheet": N-Fixation Working Group 2019c; Supporting Information Templates B). We encourage investigators to use these templates both for calculations and to provide metadata and ancillary data pertinent to the rate calculations in peer-reviewed publications, either as supplementary material or through an accession number within the publicly available OCB N 2 fixation rate database (see N-Fixation Working Group 2019a). Standardization of detection limit calculations and the publication of raw data used to calculate N 2 fixation rates will set a precedent of transparency and foster progress towards understanding the patterns and controls of oceanic N 2 fixation.
Acknowledgments This work was supported by a collaborative  NSF   EAGER grant  to P. M. ( OCE-1733610 ), J. G. ( OCE-1736659 ), and A. E. W. ( OCE-1732206 ) and a workshop supported by the  Ocean Carbon & Biogeochemistry (OCB) program. We thank  Brittany Widner  for insights on detection limits, and  Holly Westbrook  for technical assistance. We also acknowledge the  Hawaii Ocean Time-series program  for support and access to the sea aboard HOT 312. Data from this work are archived with BCO-DMO (Moisander et al. 2019).
Therefore, dedicated methods are needed to test the accuracy of NCLDV host prediction with co-occurrence networks and to improve the performance of co-occurrencebased predictions. The Tara Oceans expedition is a global-scale survey on marine ecosystems that expands our knowledge of microbial diversity, organismal interactions, and ecological drivers of community structure  (32). The present study used Tara Oceans metagenomic and metabarcoding datasets to predict virus-host relationships between NCLDVs and eukaryotes by constructing co-occurrence networks using different methods. To quantitatively assess the performance of network-based host prediction, we employed the positive likelihood ratio (LR+) using reference data for known NCLDV-host relationships. We developed a phylogeny-based enrichment analysis approach, Taxon Interaction Mapper (TIM), to enhance the performance in detecting positive signals in the intricate inferred networks.
NCLDV-eukaryote co-occurrence networks From five datasets that corresponded to five size fractions (Fig. S1 ), we generated five co-occurrence networks on a global scale (Fig. 1, S2A ). Altogether, these networks were composed of 20,148 V9 and 5,234 polB OTUs (nodes) and 47,978 polB-V9 associations (edges). Out of these associations, 47,296 had positive weights, and 682 had negative weights (Fig. 2A ). The associations that involved the family Mimiviridae were numerically dominant (n = 36,830) among the different NCLDV families. The second largest family was Phycodnaviridae, with 5,521 edges involving eukaryotes. No other family had more than 2,000 associations with eukaryotes. Marseilleviridae, forming the least associations in the networks, had 132 edges with eukaryotes. Taxonomic annotation of eukaryotic OTUs indicated that Alveolata, Opisthokonta, Rhizaria, and Stramenopiles were the major four eukaryotic groups connected to NCLDVs (with 21,167, 9,179, 6,521, and 5,327 edges, respectively).
Discussion NCLDVs can infect a wide range of eukaryotes, from unicellular to multicellular organisms  (34). However, we are still far from a comprehensive knowledge of their hosts because few have been isolated so far. Therefore, better host prediction algorithms are needed to understand the ecological functions and evolutionary significance of NCLDVs. To make these predictions, we constructed global ocean co-occurrence networks based on the marine metagenome and metabarcoding datasets from 85 stations of the Tara Oceans expedition, which cover all major oceanic provinces across an extensive latitudinal gradient from pole to pole. The edges (associations) between polB and V9 nodes (OTUs) in the networks were generated using FlashWeave. The networks were particularly dense (Fig. 2A, S2A ), thus suggesting that NCLDVs interact with numerous eukaryotes in the ocean.
In a previous study of bacteriophage host prediction, ROC curves were used as an assessment metric to compare different prediction methods  (39). However, the number of known virus-host pairs of NCLDVs is not sufficient to generate a dataset for ROC assessment. Therefore, in this study, we carried out an alternative method, the LR+, to assess the performance. LR+ is calculated with two relative values, sensitivity and specificity (Fig. 1 ). The LR+ of co-occurrence-based host predictions for positive associations was higher than 1 and increased along with increasing cut-off values for the edge weights (Fig. 4A ), which demonstrated that positive prediction results were more likely to be true positives than those based on random prediction. In high-weight regions (> 0. 6 and > 0. 4), LR+ values were larger than 10 and 4, respectively.
e. , currently unknown true interactions), indirect, or false relationships. We also found that true positive predictions only existed in positive weight associations, whereas negative weight associations did not contribute to NCLDV-host detection (Fig. S5B ). This result indicates that the abundance dynamics of NCLDVs and their potential hosts were positively correlated with each other in the analyzed samples, which were collected at a global scale; this might be because NCLDVs detected in the dataset were active viruses that replicate locally in their hosts. Similar results were obtained in other cooccurrence-based host prediction studies  (28, 41). However, several experimental studies showed that the abundance dynamics of NCLDVs and hosts showed a delay in time  (30, 31). It is possible that the global-scale samples did not have sufficiently high resolution to detect negative correlations (or correlations with a time delay) due to lack of time-resolution  (42).
S7 ). This result might be because the true positive and false positive rates defined in this study were too low; therefore, the validation may not be sufficiently sensitive to reflect the change between different abundance trimming strategies. However, it is also possible that low-abundance NCLDV OTUs are indeed network participants, as was demonstrated in a study showing that rare cyanobacterial species might play fundamental roles in blooming  (46). Our result also revealed that FlashWeave-S was better than FlashWeave-HE at predicting NCLDV-host interactions (Fig. S7 ). The difference between FlashWeave-HE and FlashWeave-S is that HE mode can remove structural zeros during network inference. Structural zero is a typical property of heterogeneous datasets, like Tara Oceans datasets, and may lead to false-positive edges  (47).
Metagenomic and metabarcoding data The microbial metagenomic and eukaryotic metabarcoding data used in this study were previously generated from plankton samples collected by the Tara Oceans expedition from 2008 to 2013  (57, 58). Because our research requires paired metagenomic and metabarcoding datasets, we used data derived from the euphotic zone samples, namely those from the surface (SRF) and Deep Chlorophyll Maximum (DCM) layers  (59). Type B DNA polymerase (polB) was used as the marker gene for NCLDVs. A total of 6818 NCLDV polB OTUs were extracted from the metagenomic datasets (i. e. , the second version of the Ocean Microbial Reference Gene Catalog, OM-RGC. v2) using the pplacer phylogenetic placement method (ML tree)  (22, 60, 61). These polB sequences were classified into seven NCLDV families (Mimiviridae, Phycodnaviridae, Marseilleviridae, Ascoviridae, Iridoviridae, Asfarviridae, and Poxviridae) and two other giant virus groups ("Medusaviridae" and "Pithoviridae").
For eukaryotes, we employed used the metabarcoding data for eukaryotes, which targeting the 18S ribosomal RNA gene hypervariable V9 region (V9)  (62). Taxonomic annotation of the eukaryotic metabarcoding data was previously performed by the Tara Oceans consortium using an extensive V9_PR2 reference database  (59) , which was derived from the original Protist Ribosomal Reference (PR2) database  (63). The diversity index of eukaryotic communities was calculated using the package "vegan"  (64). Processed frequency data are available from GenomeNet (ftp://ftp. genome. jp/pub/db/community/tara/Cooccurrence).
Data processing A relative abundance matrix for the NCLDV polB OTUs was extracted from OM-RGC. v2 for the samples derived from the pico-size fractions (0. 22-1. 6 or 0. 22-3. 0 µm). We converted the relative abundances of polB OTUs back to absolute read counts based on gene length and read length (assumed to be 100 nt). This process was required because small decimal numbers cannot be used by FlashWeave and because relative abundance data suffer from apparent correlations, which reduce the specificity of co-occurrence networks in revealing microbial interactions  (44). To build comprehensive interaction networks involving eukaryotes of different sizes, we extracted the V9 read count matrices from the metabarcoding dataset for the following five size fractions: 0. 8-5 μm; 5-20 μm and 3-20 μm (hereafter referred to as "5-20 μm" for simplicity); 20-180 μm; 180-2000 μm; and > 0.
Co-occurrence-based network inference Network inference was performed using FlashWeave [v0. 15. 0 (47) ]. FlashWeave is a fast and compositionally robust tool for ecological network inference based on the local-toglobal learning framework. Meta-variables (such as environmental parameters) can be included in the FlashWeave network to remove potential indirect associations. We used temperature, salinity, nitrate, phosphate, and silicate concentrations as meta-variables in our network inferences to determine their correlations with polB OTUs. FlashWeave provides a heterogeneous mode (FlashWeave-HE), which helps overcome sample heterogeneity. However, FlashWeave-HE may not be appropriate for the Tara Oceans data because it was shown to predict an insufficient number of known planktonic interactions  (47).
Therefore, we mainly used FlashWeave-S with default settings except for the FlashWeave normalization step and comparison between FlashWeave-S and FlashWeave-HE. A threshold to determine the statistical significance was set to alpha < 0. 01. All detected pairwise associations were assigned a value called "weight" that ranged between -1 and +1. Edges with weights > 0 or < 0 were referred to as positive and negative associations, respectively. To compare the performance of FlashWeave-S to other co-occurrence methods, we used FlashWeave-HE, Spearman, and FastSpar  (65). The FlashWeave-HE settings were the same as FlashWeave-S but with a command "heterogeneous". For Spearman, we used stats. spearmanr in package "Scipy"  (66). In FastSpar, we used 50 iterations, 20 excluded iterations, and a threshold of 0.
Network validation We validated the virus-host associations in inferred networks based on a confusion matrix defined by the known NCLDV-host information (Fig. 1 ). Briefly, we manually compiled 69 known virus-host relationships for NCLDVs (Table  S4 ). In the validation process, eukaryotic taxonomic groups were annotated at the level of the "Major lineages" in the extensive PR2 database (updated after publication)  (62). The "Major lineages" were used in the present study because 1) the deficiency of known virus-host relationships limited the use of lower eukaryotic taxonomy ranks, such as genus, for assessment, and 2) these lineages adequately represented marine eukaryotes by covering the full spectrum of cataloged eukaryotic V9 diversity at a comparable phylogenetic depth  (62). Then, we performed BLASTp [2.
10. 1  (67) ] searches from the Tara Oceans PolB sequences against the NCLDV reference database to define groups of metagenomic PolBs with a threshold of 65% sequence identity by retaining only the best hit for each environmental PolB sequence. This threshold was determined because, by using reference PolB sequences and RefSeq protein sequence databases, we found that 60-70% of sequence identity could distinguish whether the NCLDVs infected hosts of the same major lineages; this was mainly tested for Phycodnaviridae because of the lack of host information for closely related viruses in other NCLDV families (Table  S5 ). Then, 65% was chosen because it could provide a better LR+ (as described below) than 60% and 70%. The positive likelihood ratio was used in for assessment to estimate the predictions accuracy.
Phylogeny-guided filtering of host predictions and its assessment We developed Taxon Interaction Mapper (TIM) to improve host predictions by cooccurrence approaches  (33). TIM assumes that evolutionarily related viruses tend to infect evolutionarily related hosts  (17, 70) , and extract the most likely virus-host associations from the co-occurrence networks. TIM requires a phylogenetic tree of viruses (based on marker genes) and a set of connections between viruses and eukaryotes (co-occurrence edges), and then tests whether leaves (i. e. , viral OTUs) under a node of the virus tree is enriched with a specific predicted host group compared with the rest of the tree using Fisher's exact test and Benjamini-Hochberg adjustment (Fig. S9A )  (33). TIM is available from https://github. com/RomainBlancMathieu/TIM.
Virophage-NCLDV associations We inferred the networks between NCLDVs and virophages using. mcp was used as the marker gene for virophages. First, 47 reference MCP amino acid sequences were collected from public databases and used to build an HMM profile. The HMM profile was used to search against the amino acid sequences of OM-RGC v2 using HMMER hmmsearch [version 3. 3. 1] with the threshold of E-value < 1E-90  (75). This threshold was determined by searching reference sequences against the GenomeNet nr-aa database. The search detected 195 Tara Oceans virophage MCP sequences in the OM-RGC database. Together with 47 reference MCPs, a phylogenetic tree of MCP amino acid sequences was built using MAFFT and FastTree. We extracted the abundance profiles for the 195 MCP sequences from the pico-(0.
22-1. 6 or 0. 22-3. 0 µm) and femto-size (< 0. 22 µm) fractions. We used samples from the SRF and DCM depths. PolB and MCP abundance profiles were merged into two matrices corresponding to the two virophage size fractions. Then, network inference was conducted using the FlashWeave default settings after Q1 filtration. In the MCP phylogenetic tree, three virophage clades contributed most of the NCLDV connections. Thus, an NCLDV enrichment analysis for the three clades was carried out using Fisher's exact test, and the p-value was adjusted by the Benjamin-Hochberg method. This approach was the same as TIM, but we did not use the TIM software because the current version of TIM requires inputs of eukaryotic nodes with NCBI taxonomy annotations. We used another approach, HGT, to predict the virophage-NCLDV interactions.
First, we generated an NCLDV genome database, which includes 56 reference NCLDV genomes corresponding to our polB dataset and 2,074 metagenome-assembled genomes from a previous study  (24). A total of 827,548 coding sequences were included in this database. Then, 195 virophage MCPs from the metagenomic data were BLASTp searched against this database using an E-value cut-off of 1E-10 (with a minimum query coverage of 50% and a minimum sequence identity of 50%). If a virophage MCP obtained a hit in the NCLDV genome database with a lower E-value compared with hits in the MCP database (the hit to itself was removed), the hit in the NCLDV genome database was considered an HGT candidate.
Acknowledgments This work was supported by  JSPS/KAKENHI  (Nos. 18H02279  and  19H05667  to H. O. and Nos. 19K15895  and  19H04263  to H. E. ),  Kyoto University Research Coordination Alliance  (funding to H. E. ), and the  Collaborative Research Program of the Institute for Chemical Research, Kyoto University  (Nos. 2019-30  and  2020-27 ). Computational time was provided by the  SuperComputer System, Institute for Chemical Research, Kyoto University. We further thank the  Tara Oceans consortium , and the people and sponsors who supported  Tara Oceans. Tara Oceans  (including both the Tara Oceans and  Tara Oceans Polar Circle  expeditions) would not exist without the leadership of the  Tara Expeditions Foundation  and the continuous support of 23 institutes (https://oceans. taraexpeditions.
, 2006). As head size is expected to directly impact feeding efficiency in gape-limited predators like snakes  (Esquerré & Keogh, 2016; Forsman, 1996; Grundler & Rabosky, 2014) , we also quantified the evolutionary allometry in our dataset. In the second part of this study, we explored the functional implications of the observed shape variation. All considered species successfully capture aquatic prey despite the hydrodynamic constraints they face  (Segall et al. , 2019; Van Wassenbergh et al. , 2010). As these constraints are related to head shape  (Fish, 2004; Godoy-Diana & Thiria, 2018; Koehl, 1996; Polly et al. , 2016) , we expected the observed morphological disparity to have functionally converged (i. e. , have the same hydrodynamic profile) which would indicate a many-to-one-mapping of form to function  (Wainwright et al.
Vogel (1994)  suggested that propulsion-based organisms should be under a selective regime that favors a reduction in acceleration reaction by reducing mass and/or added mass. Both drag and acceleration reaction are linked to the properties of the fluid, the kinematics of the motion, a scaling component, and a shape component. If snakes are under selection and yet display a large head shape disparity, then we can expect a many-to-one mapping of form to function, with several shapes resulting in a reduction in both drag and added mass. To test this hypothesis, we measured the shape-component of both hydrodynamic forces (i. e. , drag and added mass coefficient) of shapes that are representative of the morphological disparity of our dataset. We first compared the head shape of 62 species of snakes that capture elusive aquatic prey under water by scanning the surface of the head of more than 300 specimens from museum collections.
| Geometric morphometrics We created a template consisting of a set of 921 landmarks with 10 anatomical landmarks, 74 semilandmarks on curves corresponding to anatomical features, and 837 surface semilandmarks (Figure  2 ). We manually placed all the landmarks on the template (anatomical, curve, and surface landmarks), and only the anatomical landmarks and curve semilandmarks on all specimens using the Landmark software  (Wiley et al. , 2005). We ensured the reliability and repeatability of the landmark positioning (see Supplementary Material 2). Next, the surface semilandmarks were projected on each specimen, and both curve and surface semilandmarks were relaxed and slid by minimizing the bending energy  (Gunz & Mitteroecker, 2013)  using the "Morpho" package  (Schlager, 2017). We then obtained a consensus shape for each species by performing a generalized Procrustes analysis (GPA) for symmetrical shapes on all the specimens of each species using the function "procSym" of the "Morpho" package (R script available in Supplementary Material 3).
03, p < 10 -12 , R =. 75), we only used the centroid size to test for allometry. We tested for an interaction between size and diet by adding interactions to the model. We assessed the statistical significance of the variables by using 10,000 simulated datasets obtained by permuting the phenotypic data across the tips of the phylogeny. We extracted the shapes associated with allometry (named "smaller" and "larger") by using the function shape. predictor in "geomorph"  (Adams, 2014b). The shapes associated with the different dietary groups (named "bulky" and "elongated") were obtained by performing a GPA on the species belonging to each dietary group. We extracted the resulting consensus along with their centroid sizes. Then, we performed another GPA on the rescaled consensus of the groups to obtain the models in the same morphospace.
We then generated meshes from the different configurations using MeshLab  (Cignoni et al. , 2008 ) and compared them using the function meshDist in "Morpho. " Because the shape variability might be structured by other factors than diet and size, we used an unsupervised pattern recognition based on Gaussian Mixture Modelling (GMM) implemented in the  Pyron & Burbrink, 2014) , head shape, shape of the preferred prey of the 62 selected species (oval: elongated prey, square: bulky prey; see Supplementary Material 1 for references)  The GMM detects whether our dataset can be decomposed in subgroups. As this method is sensitive to the number of variables, we only used the first seven principal components (PC; 90% of the shape variability) as input. This model-based clustering algorithm assumes that the input variables (here the PCs) have a Gaussian distribution.
The function searches for clusters in the dataset, based on the distribution of species in the morphospace by trying to fit several predefined distribution models (for details on models, see  Fraley et al. , 2012). It uses a hierarchical procedure that first considers each species as a single cluster and agglomerates the clusters based on a maximum-likelihood approach. The process stops when all species are gathered into one single cluster. Then, the Bayesian information criterion of each cluster model is calculated to determine which model best fits the repartition of species in the dataset  (Cordeiro-Estrela, Baylac, Denys, & Polop, 2008; Fraley & Raftery, 2003). F I G U R E 1 Phylogenetic relationship (pruned from.
| 3D models To test our hypothesis of many-to-one mapping of form to function, we characterized the hydrodynamic profile of five head models that describe the main axes of variability in our dataset. Thus, we chose to work on the extreme shapes described by the first two PCs, as these components represent 65. 1% of the overall head shape variability, and the mean shape (Figures  3  and 4 ). PC1 represents more than 54. 6% of the variability and separates species with long and thin heads on its negative part from species with bulkier and shorter heads on its positive part (Figure  3 ). PC2 represents 10. 5% of the variability and separates species having a horizontally flattened head from species with a more circular head (Figure  3 ). Aquatic snakes strike at their prey with the mouth open at various angles depending on the species, ranging from 40°  (Alfaro, 2002)  to 80°  (Herrel et al.
| Morphometric analyses The head shape of snakes showed a significant phylogenetic signal (p =. 001, K mult = 0. 37). Prey shape, size, and the interaction between the two factors show a significant impact on head shape (D-PGLS: P prey = 0. 008, P size = 0. 002, P prey*size = 0. 004). Allometry and diet, respectively, represent 7. 6% and 5. 6% of the overall variation in our dataset (D-PGLS R 2 ). In snake species that prefer bulky prey, the region between the eyes and the snout is enlarged compared to elongated prey-eating snakes (Figure  5 ). The upper jaw is slightly enlarged at its rear part for bulky-prey eaters, whereas the lower jaw appears more robust in elongated prey eaters. The rear part of the head is enlarged in elongated prey eaters, especially on the sides, resulting in a more tubular shape, while the bulky-prey eaters show a reduction in the head girth in this region.
The interaction between size and dietary preference highlighted in the linear model suggests that elongated prey-eating species have smaller heads and a shape that is a combination between Figure  5c  and Figure  6c. The Gaussian mixture model applied to 90% of the variability (i. e. , 7 first principal components) returned a unique component suggesting little or no structure described by mixtures of datasets with normal distributions: The species in our dataset form a single cluster, which is visible in the morphospace (Figure  3 ). The variability in head shape is mostly carried by outlier species.