Gene clustering and identification of Red Sea-associated ortholog groups.

The SAGs described here come from an undersampled region of the ocean (the Red Sea) and in part from undersampled clades of marine bacteria (SAR11 subclades Ib, Id, and II) and therefore provide the opportunity to identify OGs specific for these clades or possibly endemic to this ocean region. To investigate these patterns, we combined the Red Sea SAGs with available cultured genomes and SAGs (separately for Prochlorococcus and SAR11), clustered genes into OGs using a Markov clustering algorithm (OrthoMCL; see Materials and Methods), and identified OGs found only in the Red Sea SAGs (RS-OGs) and OGs found only in certain clades.

We identified 878 SAR11 RS-OGs and 96 Prochlorococcus RS-OGs, that is, OGs not found (in this analysis) in genomes from outside the Red Sea (File S1). These totals represent 16.7% of all (19.1% of noncore) SAR11 OGs and 0.9% of all (1.0% of noncore) Prochlorococcus OGs. Many of the RS-OGs were found only in a single clade: 96 in Prochlorococcus ecotype HLII, 484 in SAR11 subclade Ia, 101 in SAR11 subclade Ib, 101 in SAR11 subclade Id, and 132 in SAR11 subclade II. The numerous clade-specific OGs present targets for understanding ecotype-specific physiology.

The first pattern of note was that there were more RS-OGs in the SAR11 SAGs than in the Prochlorococcus SAGs. This reflects the large contribution of our SAR11 SAGs to the sequenced SAR11 pangenome: the number of SAR11 Red Sea SAGs (21) was nearly as high as the number of SAR11 reference genomes (26). In contrast, the number of Prochlorococcus Red Sea SAGs (5) was only 3% of the number of Prochlorococcus reference genomes (140). Emphasizing the effect of the genome reference database on estimates of OG endemicity, after new Prochlorococcus genomes (9, 28) were added to the clustering, the number of RS-OGs dropped from 1,192 to 96 (Fig. S4). Another explanation for the greater number of new SAR11 OGs is that the SAR11 SAGs span previously unsampled or undersampled clades: 334 of the 878 Red Sea-associated SAR11 OGs were found in only one of subclade Ib, Id, or II. Furthermore, SAR11 is a broader phylogenetic group, based on 16S rRNA diversity, than Prochlorococcus (31); therefore, its pangenome may be expected to be larger. In summary, we suspect that the large number of new SAR11 OGs (878) in general more likely reflects the current dearth of sequence data for SAR11 rather than a significant degree of endemism due to isolation and/or selection.

The second pattern we examined was inspired by our question about possible endemic gene content in the Red Sea: based on the geographic isolation of the Red Sea and its unique combination of physicochemical conditions (simultaneously high irradiance, high salinity, high temperature, and low nutrients), do genomes isolated from the Red Sea exhibit endemic OG content encoding adaptive functions for this environment? The answer that emerged to this question is that there were some indications of possible endemic adaptations to the Red Sea, but there were no new pathways identifiable among the RS-OGs, most of the RS-OGs with annotated functions were found in only one or two SAGs, and the majority of RS-OGs encoded hypothetical proteins with no assigned function.

The majority of RS-OGs were hypothetical proteins: 82% (723 of 878) for SAR11 and 91% (87 of 96) for Prochlorococcus. It was difficult to infer possible adaptive functions for OGs with no predicted functions; however, these OGs may be referenced later when new approaches for annotating conserved hypotheticals are developed. The remaining nonhypothetical OGs (155 SAR11, 9 Prochlorococcus), i.e., those with predicted functions, are listed in Table S2. While we could not detect a widespread signature of adaptation to the Red Sea environment, i.e., RS-OGs with annotated functions represented across multiple SAGs, below we highlight a few sparsely represented RS-OGs that may have adaptive functionality in the Red Sea environment, some with biotechnological potential.

Among Prochlorococcus SAGs, none of the 9 nonhypothetical RS-OGs (Table S2) were found in more than one SAG. One OG (proch20425) found in SCGC AAA795-M23 encodes UvrABC system protein B, responsible for repair of DNA damage. We could posit that this enzyme is found preferentially in the Red Sea because of the year-round high irradiance, which increases the rate of DNA damage in cells.

Among SAR11 SAGs, there were 21 nonhypothetical RS-OGs found in two or more SAGs and another 134 found in only one SAG (Table S2). These OGs show links to high light adaptation, motility, and nitrogen and phosphorus assimilation. One OG (pelag14710, found in one SAG) encodes a photolyase enzyme that repairs DNA damage caused by exposure to UV light. Pyrophosphatase (pelag15064; found in one SAG) is involved in the hydrolysis of inorganic pyrophosphate into two orthophosphates and may have a role in phosphorus utilization. Allantoinase (pelag15247) and urease accessory protein UreF (pelag14490) are each found in one SAR11 SAG. These enzymes involved in phosphorus and nitrogen metabolism may provide an adaptive advantage in the Red Sea, which exhibits colimitation to both elements and may be relatively more nitrogen limited (12, 15). Several of the SAR11 RS-OGs encode enzymes with biotechnological relevance. DNA polymerase I (pelag12679, pelag14776, and pelag14807) from this higher temperature environment could have heat-resistant properties, for example, marginal thermostability conferred by amino acid substitutions (32).

After the major analyses had been completed for this study, two SAR11 genomes (21) and one Prochlorococcus genome (20) derived from cultivated strains were sequenced, and four Prochlorococcus genomes were assembled from metagenomes (33). Of the SAR11 genomes, one was assigned to subclade Ia and the other to subclade Ib (21). Of note, the subclade Ia genome (RS39) contained several OGs also found among the Red Sea-associated SAR11 OGs: 3-oxoacyl-acyl-carrier-protein synthase, ABC branched amino acid transporter, arylsulfotransferase, formate dehydrogenases, glycosyl transferases, methyltransferases, sialic acid synthase, sucrose synthase, sulfotransferases, and a type II restriction-modification system. Several of these functions may play roles in one-carbon and sugar metabolism by SAR11 in the Red Sea (21). The Prochlorococcus genome was assigned to the HLII ecotype and notably contained a pathway for biosynthesis of the osmolyte (compatible solute) glucosylglycerol (20). This pathway represents a possible adaptation to the higher salinity of the Red Sea. However, the three genes in this pathway were not found among the Red-Sea-associated Prochlorococcus OGs, and they were not found elsewhere among the retained or removed contigs from the Red Sea SAGs (BLASTN).

Expression of ortholog groups in the Red Sea water column.

To further test the idea that there could be OGs of ecological importance endemic to the Red Sea, we analyzed metatranscriptomes from the Red Sea. Any RS-OGs with functional roles would be expected to be expressed in the Red Sea water column. We collected seawater and filtered the prokaryotic fraction from a station in the central Red Sea over a broad temporal and depth range: samples were collected at four depths and 13 time points over a 48-h period. We extracted and sequenced RNA from these samples and mapped the reads to the Red Sea SAGs.

We found that around two-thirds of RS-OGs were expressed in one or more sample: 64% for SAR11 (Fig. 3b) and 66% for Prochlorococcus (Fig. 3d). This was more than the fraction of non-RS-OGs expressed: 32% for SAR11 (Fig. 3a) and 20% for Prochlorococcus (Fig. 3c). We were curious if the high fraction of non-RS-OGs that were unexpressed was due to many of these OGs being singletons (OGs having only one member). To the contrary, heatmaps of OG size versus number of metatranscriptomes in which the OG was found (Fig. 3, inset) do not show a high density of singleton OGs having no expression in non-RS-OGs; rather, the trend toward singletons is more common in RS-OGs.

Open in a separate window

FIG 3

Expression of SAG OGs in Red Sea metatranscriptomes. The 52 metatranscriptomes span a broad range of the water column at a station in the central Red Sea: four depths and 13 time points over a 48-h period (every 4h). Histograms show the number of metatranscriptomes found in SAR11 non-RS-OGs (a), SAR11 RS-OGs (b), Prochlorococcus non-RS-OGs (c), and Prochlorococcus RS-OGs (d). Heatmaps (inset) show the density of OGs based on OG size (number of total copies across the SAGs) and the number of metatranscriptomes an OG is found in. RS-OGs were more likely than other OGs to be expressed in one or more samples, and non-RS-OGs that were expressed were more likely to be expressed in a large number of samples.

Of OGs expressed in at least one sample, non-RS-OGs (Fig. 3a and ​andc)c) tended to be expressed in more samples than RS-OGs (Fig. 3b and ​andd).d). This is consistent with many of the non-RS-OGs being core genes, many of which are housekeeping genes that are often constitutively expressed. Overall, the expression patterns indicate that the majority of RS-OGs are transcribed to mRNA, consistent with the synthesis of functional gene products.

Distribution of ortholog groups across the global ocean.

The analysis to this point has focused on the distribution of OGs among cultured and single-cell genomes and their expression in the Red Sea water column. A set of OGs has been found that is exclusive to Red Sea genomes (to date), and a majority of them are expressed in the water column. However, we cannot rule out the possibility that these OGs appear endemic only because more genomes are not available from around the World Ocean. If we extended our search to global marine metagenomes, instead of just genomes, would we in fact find these putative endemic OGs in other seas?

To investigate the possibility that, contrary to our original hypothesis, there may be few OGs truly endemic to the Red Sea microbial community, we analyzed metagenomes collected from across the global ocean by the Tara Oceans expedition. We searched for SAR11 and Prochlorococcus OGs in 139 prokaryote-fraction metagenomes from the Tara Oceans expedition (34), which come from several depths in the water column: surface, mixed layer, deep chlorophyll maximum, and mesopelagic zone. We queried the data set to determine what fraction of all OGs and what fraction of RS-OGs could be found outside the Red Sea. If RS-OGs represent gene content endemic to the Red Sea, we would expect to find them absent from metagenomes from other regions. Our approach was complementary to a recent study that analyzed the global metapangenome of Prochlorococcus in the Tara metagenomes, showing the distributions of gene clusters (OGs) with strain-level resolution across the Tara samples (35). In the work here, we employed rarefaction and ordination techniques, with a particular focus on RS-OGs.

The presence or absence of SAR11 and Prochlorococcus orthologs in Tara Oceans prokaryote-fraction metagenomes (Files S7 and S8) was plotted as rarefaction curves (Fig. 4). Tara Oceans metagenomes were added randomly one by one, and the fraction of SAR11 and Prochlorococcus OGs found was tallied and plotted. The rarefaction curves show the averages ± standard deviations from 1,000 permutations. They also show the best-case (and worst-case) scenarios, that is, the fraction of OGs found if each new metagenome adds the most (or fewest) new OGs. Between 70 and 85% of OGs could be found in one or more Tara Oceans metagenomes (Fig. 4), and in the best-case scenarios it took at most ten metagenomes to find 90% of these OGs (Table S3). The percentage of OGs not found (15 to 30%) was independent of whether they were Red Sea associated or not. This result, combined with the rarefaction analysis, suggests these OGs are unlikely to be found in the Tara samples with deeper sequencing. It is possible that some OGs are rare and/or divergent enough to be undetectable with the current methodological approach.

Open in a separate window

FIG 4

Rarefaction analysis showing the proportion of all OGs (a and c) and RS-OGs of SAR11 (b and d) and Prochlorococcus observed in Tara Oceans metagenome samples. Curves show the cumulative number of OGs observed in Tara Oceans samples (E value of <1e−5) as more samples are added. Yellow lines show the averages ± standard deviations from 1,000 permutations of randomly added samples. Blue lines show the best-case scenario (each sample added is that with the most number of new OGs observed) and worst-case scenario (each sample added is that with the fewest number of new OGs observed). Red lines show the means from 1,000 permutations of randomly added samples but with Red Sea samples (031_SRF_0.22-1.6, 032_DCM_0.22-1.6, 032_SRF_0.22-1.6, 033_SRF_0.22-1.6, 034_DCM_0.22-1.6, and 034_SRF_0.22-1.6) added last. As more Tara metagenome samples are added to the analysis, the number of new OGs identified approaches a plateau where new samples do not reveal many new OGs. The same is true with RS-OGs, even when samples from the Red Sea are added last, with the exception of 5 Prochlorococcus OGs (proch20367, proch20368, proch20390, proch20423, and proch20438).

Across the 139 Tara Oceans prokaryote-fraction metagenomes, we found 84.9% (4,475/5,272) of all SAR11 OGs in one or more metagenomes (leaving 15.1% not found) (Fig. 4a) and 72.2% (7,537/10,439) of all Prochlorococcus OGs in one or more metagenomes (leaving 27.8% not found) (Fig. 4c). In the best-case scenarios, it took only 5 metagenomes to find 90% of the found SAR11 OGs and 50 metagenomes to find 99%; it took only 10 metagenomes to find 90% of the found Prochlorococcus OGs and 60 metagenomes to find 99% (Table S3). The fractions of OGs found were similar for RS-OGs, where 81.2% (713/878) of SAR11 OGs were found (leaving 18.8% not found) (Fig. 4b) and 69.8% (67/96) of Prochlorococcus OGs were found (leaving 30.2% not found) (Fig. 4d). That is, RS-OGs were about as likely to be found across the World Ocean as non-RS-OGs. For both SAR11 (Fig. S5a) and Prochlorococcus (Fig. S5b), considering the number of Tara metagenomes in which each OG was found, RS-OGs were less likely to be found in a large fraction of metagenomes than all OGs. This is not surprising: the set of non-RS-OGs contains all of the core OGs, which would be expected to be found in most if not all samples.

To evaluate whether Tara Red Sea metagenomes contained any RS-OGs not already found in the non-Red Sea metagenomes, we tested scenarios where the Red Sea metagenomes were added last in the rarefaction analysis. There was no change in the mean curve of cumulative SAR11 OGs found when the six Tara Red Sea metagenomes were added last (Fig. 4b): all of the SAR11 RS-OGs could be found without examining the Red Sea metagenomes. In contrast, there were five Prochlorococcus RS-OGs that were added to the cumulative total when the Tara Red Sea metagenomes were added last (Fig. 4d). These five OGs, all with unknown function, represent a small fraction of the total Prochlorococcus pangenome (10,439 OGs total). Given the available genomes, this study may have uncovered a small set of OGs (Table S2) that possibly reflect gene content endemic to or generally associated with Red Sea environmental conditions, and this marks an area for further research. In light of this metagenomic analysis, however, it appears that the putative RS-OGs provide a relatively minor contribution to the whole and that these new SAR11 and Prochlorococcus genomes from the Red Sea generally reflect global pangenomes.

Finally, we were curious if OG composition as a whole could show the Red Sea metagenomes to be different from the other metagenomes, despite the lack of evidence of endemic OGs. More generally, could the relative abundance of OGs across Tara be used to distinguish populations of Prochlorococcus and SAR11?

We used the tables of OG counts in the 63 Tara surface (SRF) prokaryote-fraction metagenomes to do PCA ordination on the Tara metagenomes (Fig. 5; top OGs driving separation among the metagenomes are provided in Table S4). SAR11 OG composition (Fig. 5a) was not obviously structured by temperature differences in the temperate and tropical ranges, although Red Sea samples clustered together, and polar samples were separate from the others. Prochlorococcus OG composition (Fig. 5b), however, was structured by temperature differences in the temperate and tropical ranges. The four Red Sea samples were split, with two samples clustering with the warm samples and two samples with the cooler samples. These Red Sea samples are positioned where they would be expected based on temperature: the two southern samples (latitudes 18.4°N and 22.0°N) were warmer (temperature, 27.6°C and 27.3°C) and clustered with other warm/tropical samples (Fig. 5b, left side of PC1); the two northern samples (latitudes 23.36°N and 27.16°N) were cooler (temperature, 25.8°C and 25.1°C) and clustered closer to the cool/temperate samples (Fig. 5b, right side of PC1). Note these temperatures are lower than average for Red Sea surface waters because the Tara Red Sea samples were collected in winter (January); in contrast, the Red Sea samples in the World Ocean Database (described above) were collected in spring (April). Given that temperature tolerances generally lack known genetic markers (36), these data suggest an area for future investigation.

Open in a separate window

FIG 5

Principal component analysis of Tara Oceans surface samples by the abundance of SAR11 (a) and Prochlorococcus (b) OGs. The ordination shows the similarity of Tara Oceans samples to each other along the first two principal components. Samples are colored by Tara temperature categories: polar samples (<10°C) are dark blue, temperate samples (10 to 20°C) are light blue, tropical samples (>20°C) are orange, and Red Sea tropical samples are orange with black edges. Red Sea samples and Tara samples generally show more separation based on temperature when ordinated by Prochlorococcus OG composition than by SAR11 OG composition.

In summary, the analysis of Prochlorococcus and SAR11 OGs in Tara Oceans metagenomes shows that (i) most Red Sea-associated OGs are actually widely distributed across the World Ocean, not endemic to the Red Sea, and (ii) OG distribution patterns as a whole, taking relative abundance into account, place the Red Sea on a continuum with other seas, with patterns explained by environmental factors, including temperature. Supporting this idea, differences in the relative abundance of OGs, with physicochemical properties covarying with OG functions, have been observed among the North Pacific, Sargasso Sea, Mediterranean Sea, and Red Sea in previous comparative metagenomics studies (11, 12). Despite the Red Sea existing at the periphery of multiple physicochemical parameters in the World Ocean, its distinctiveness may best be revealed by the relative abundance of OGs rather than in the wholesale presence or absence of OGs. In addition to this general pattern, this effort also identified a small set of putative and nonhypothetical proteins that warrant further ecological and biotechnological study.

Nucleic acid extraction and amplification.

Single bacterioplankton cells in the preserved samples were flow sorted, whole-genome amplified (MDA, or multiple displacement amplification), and PCR screened at the Bigelow Laboratory Single Cell Genomics Center (SCGC; Boothbay Harbor, ME, USA), by following previously described protocols (37), with SYTO-13 nucleic acid stain used to stain cells for flow sorting. SAG identification was carried out with SCGC protocol S-102 for bacteria using 16S rRNA primers 27F and 907R (38, 39). Totals of 21 and 5 cells were identified from 16S PCR screening and subjected to a second round of MDA before sequencing.

The RNA extraction protocol for metatranscriptomics was adapted from references 40,–42. After expelling RNAlater from the Sterivex filter, 2ml lysozyme solution (1mg/ml in lysis buffer, containing 40mM EDTA, 50mM Tris, pH 8.3, 0.73M sucrose) was added and then filter incubated at 37°C with rotation for 45min. Proteinase K solution (50 μl at 20mg/ml; Qiagen/5PRIME) and SDS solution (100 μl at 20%) were added, and then the filter was incubated at 55°C with rotation for 2h. Lysate was expelled to a separate tube; meanwhile, 1ml lysis buffer was added to the filter to wash at 55°C for 15min. The two lysates were pooled, to which was added 1.5ml absolute ethanol. RNA was then extracted from this solution using the RNeasy Protect bacterial midi kit (Qiagen). RNA was eluted with two volumes of RNase-free water. RNA sample was concentrated using a speed vacuum, from 250 μl to 60 μl. To this volume we added DNase (1 μl Ambion TURBO DNA free, 6 μl 10× buffer, 60 μl RNA) and incubated at 37°C for 30min. This solution was purified using the RNeasy MinElute cleanup kit (Qiagen) and eluted with RNase-free water. The final yield was 1 to 2ng total RNA. Total RNA was amplified using the C&E Version ExpressArt bacterial mRNA amplification nano kit, which preferentially amplifies mRNA [independent of poly(A) tail] and selects against rRNAs. A single round of amplification was performed on 2 to 4ng of total RNA, which yielded about 10μg final amplified RNA.

Nucleic acid sequencing.

For single-cell genome sequencing, genomic library preparation with Illumina TruSeq and sequencing with Illumina GAIIx and Illumina HiSeq 2000 was done at the KAUST Bioscience Core Laboratory, generating paired 105-bp reads. The assembled contigs (assembly methods are described below) are available from NCBI with accession numbers PRJEB9287 (BioProject) and SAMEA3368552 to SAMEA3368577 (BioSample) and from Integrated Microbial Genomes (43) with accession numbers 2630968236, 2630968238 to 2630968254, 2630968277 to 2630968281, and 2630968285 to 2630968287.

For metatranscriptomics, sequence data were processed as described in reference 20. Amplified RNA was used to construct sequencing libraries using the TruSeq stranded RNA LT sample preparation kit (Illumina) according to the manufacturer’s protocol. Libraries were paired-end sequenced with the Illumina HiSeq 2000 platform (2×100bp). A link to the raw FASTQ sequences is provided at the end of Materials and Methods. Low-quality reads and sequencing adapters were removed using Trimmomatic v0.32 (44). Sequence reads shorter than 50bp were discarded. Bowtie 2 v2.2.4 (45) was used to identify and remove PhiX contamination sequences. The remaining sequences were error corrected using the BayesHammer algorithm (46) implemented in SPAdes v3.5.0 (47), followed by removal of putative rRNA gene transcripts with SortMeRNA v2.0 (48).

Genome assembly and annotation.

De novo assemblies were generated using CLC Genomics Workbench 4.9. The genomes were assembled independently, and unless otherwise specified, the following applies to all of the SAGs. The reads were first imported and quality trimmed with a limit of 0.01. They were then assembled using CLC’s de novo assembler with a word size (k-mer) of 64 and with the minimum/maximum of the insert size set to 100/1,000bp. Only those contigs greater than 200bp in length were included in downstream analyses. The reads were mapped to the consensus sequence of the assembled contigs using CLC’s default parameters but with the length fraction set to 1.0 and the similarity set to 0.95.

Assembled SAG contigs were ordered and oriented relative to SAR11 HTCC1062 (NC_007205.1) or Prochlorococcus MIT 9202 (NZ_DS999537) using ABACAS 1.3.1 (49). The ordered sequences were then imported into GAP4 (50), and additional joins were made between overlapping contigs if conserved synteny supported the arrangement. To identify and remove possible contaminating sequences from the assemblies, each contig was retained only if it met one or both of the following criteria: (i) the contig was binned into a bin annotated as SAR11 or Prochlorococcus using Metawatt 3.5 (22), using the medium bin level, with a minimum bin size of 50 kbp and minimum contig size of 500bp; (ii) the contig had a top-10 BLASTN hit against the GenBank nucleotide database, with an E value of <1e−5 to SAR11 or Prochlorococcus.

Prediction of gene open reading frames (ORFs) and functional annotation of SAGs were performed by the RAST web service (51) with FIGfam, release 59.

OG clustering.

Predicted proteins from SAGs were clustered with proteins from published cultured and SAG genomes (see File S1 in the supplemental material) into OGs using OrthoMCL 2.0 (52). OrthoMCL configuration settings were a percentMatchCutoffof50 and evalueExponentCutoff of −5. This yielded 5,272 SAR11 OGs and 10,439 Prochlorococcus OGs. After OrthoMCL clustering, OGs were assigned as core and noncore based on copy number in the non-Red Sea, cultured (non-SAG) genomes: core OGs are those found at least once in each of the non-Red Sea, cultured genomes, and noncore OGs are those not found in at least one of the non-Red Sea, cultured genomes. Among SAR11, there were 683 core OGs and 4,589 noncore OGs. Among Prochlorococcus, there were 1,152 core OGs and 9,287 noncore OGs. Protein sequence identifiers and FASTA sequences for each OG have been archived at https://zenodo.org (https://doi.org/10.5281/zenodo.2634561).

Estimation of genome completeness.

Completeness of SAGs was assessed using two methods. First, completeness was assessed using single-copy core OGs, i.e., those OGs found once and only once in each complete genome based on the OrthoMCL clusters (analyzed separately for SAR11 and Prochlorococcus). Completeness was calculated as the number of core orthologs present in each SAG out of 649 SAR11 or 1,144 Prochlorococcus single-copy core OGs. Second, genome completeness of the SAGs was assessed using CheckM 1.0.13 (53) using lineage-specific workflow (lineage_wf) with database file checkm_data_2015_01_16.tar.gz downloaded from https://data.ace.uq.edu.au/public/CheckM_databases; CheckM was also used to estimate genome redundancy (called contamination in CheckM). For comparison, CheckM completeness and redundancy were calculated for the reference genomes used in this study (Table S1).

Genome taxonomy and phylogenetics.

A total of 89 SAR11 and 96 Prochlorococcus shared single-copy orthologous genes were identified using the GET_HOMOLOGUES software (54). Amino acid sequences translated from gene sequences were aligned using the MAFFT software (55). These alignments were concatenated, sites with gaps were deleted, and the concatenated data were partitioned using the PartitionFinder software (56) to account for variations of evolutionary processes among gene families. With the Bayesian information criterion (BIC) statistic, a 16-partition framework was chosen to optimally describe the variability, in which the LG rate matrix with gamma distribution of rate variation (LG+G) was selected for 15 partitions and the VT rate matrix with gamma distribution of rate variation (VT+G) was selected for the remaining partition. This partition model was used in the maximum-likelihood phylogenomic construction using the RAxML software (57).

Ordination of SAGs and genomes using k-mer composition and ortholog composition.

SAGs and reference genomes (Table S1) were analyzed using PCA of nucleotide composition and OG composition. Nucleotide composition of the SAGs and reference genomes (SAR11 and Prochlorococcus scaffolds of >200 kbp from Integrated Microbial Genomes; https://img.jgi.doe.gov) was determined as 6-nucleotide words or k-mers (6-mers). k-mer frequencies were calculated using Jellyfish 2.2.5; the main command used was jellyfish count –m 6 –t 8 –s 1M. This resulted in a table of 6-mer frequencies in the SAGs and genomes, one table each for SAR11 and Prochlorococcus. OG composition was derived from tables of OrthoMCL clusters, which, as the SAGs had variable levels of completeness and gene counts (Table 1), were randomly subsampled so that all genomes had the same number of gene counts in the tables; five replicate subsamples produced very similar results (only one SAR11 subsample and one Prochlorococcus subsample are shown). The number of OG counts subsampled was chosen to balance the number of OG counts with the number of genomes retained (less complete SAGs were excluded): the OG composition tables (with counts of 5,272 unique SAR11 OGs and 10,439 unique Prochlorococcus OGs) were subsampled down to 800 gene counts per SAR11 SAG (keeping 12 of 21 SAGs) and 1,400 gene counts per Prochlorococcus genome (keeping 5 of 5 SAGs). Prior to PCA, a pseudocount of 1 was added to k-mer and OG count tables to account for zero values; k-mer counts were then converted to relative abundances for each genome (unnecessary for OG counts because of the subsampling procedure); k-mer relative abundances were then standardized to z scores (not done for OG counts because this reduced the resolving power of PCA). PCA was then performed using the Scikit-Learn function sklearn.decomposition.PCA (58).

Mapping of metatranscriptomic reads to OGs.

The quality-filtered mRNA reads from the 52 samples were mapped against the SAGs using Bowtie 2 (45) with default settings. Each read mapping above the threshold was assigned to exactly one gene in a SAG contig. The resultant read counts were normalized based on the FPKM (fragments per kilobase of gene per million mapped reads) metric. Per-sample FPKM counts for each gene were then summed by OGs, resulting in per-sample FPKM counts for each OG. For downstream analysis, counts were converted to a simple presence/absence measure: if any gene belonging to the OG had one or more mapped transcript, that OG was marked as present in that sample.

Detection and rarefaction analysis of OGs in Tara Oceans metagenomes.

A set of 139 prokaryote-enriched Tara Oceans metagenomic gene files (34) was downloaded from the European Nucleotide Archive (https://www.ebi.ac.uk/ena; accession numbers ERZ096909 to ERZ097150). Each file contains nucleotide sequences for genes predicted on Tara Oceans metagenomic contigs that were assembled from shotgun sequencing reads from individual Tara Oceans samples. The prokaryote fraction was 0.22 to 1.6μm for stations 004 to 052 and 0.22 to 3μm for stations 056 to 152; the environmental features of the samples were indicated as SRF (surface), MIX (mixed layer), DCM (deep chlorophyll maximum), and MES (mesopelagic zone). The metagenomic gene sequences were queried against a database of translated proteins from the SAGs and genomes with DIAMOND 0.8.26 (59) using the program blastx with parameters –p 40 –k 25 –e 1e−3. The top hit (SAG or genome protein sequence) for each Tara gene sequence (E value of <1e−5) was retained. E value cutoffs of 1e−10 and 1e−15 were also tested, which showed the same trends as an E value of <1e−5 but with fewer total OGs identified. Counts of the number of times each protein was a top hit were then summed across each OG. This resulted in a table of OGs by samples where each OG was either present (at least one constituent protein was a top hit at least once) or absent in each sample. These presence/absence tables (one for SAR11, one for Prochlorococcus) were used to generate rarefaction curves: samples were added one-by-one randomly (1,000 permutations), and the cumulative number of OGs found was recorded.

Ordination of Tara Oceans metagenomes by OG composition.

OG counts (total, not presence/absence) in Tara Oceans surface (SRF) sample metagenomes were used for ordination using PCA. Prior to PCA, a pseudocount of 1 was added to OG count tables to account for zero values; counts were then converted to relative abundances for each metagenome; OGs with an average relative abundance across all metagenomes of less than 0.0001 (0.01%) were removed; relative abundances were then standardized to z scores. PCA was then performed using the Scikit-Learn function sklearn.decomposition.PCA (58).

World ocean temperature and salinity data.

Surface temperature and salinity data (WOD13_ALL_SUR_OBS) from the World Ocean Database 2013 (https://www.nodc.noaa.gov/OC5/WOD13/) were downloaded from the Research Data Archive at the National Center for Atmospheric Research (https://rda.ucar.edu/datasets/ds285.0/).

We thank Haiwei Luo for assistance building genome trees, Mamoon Rashid for consultation about decontamination methods, Qiyun Zhu for assistance with genome analysis, and Ramunas Stepanauskas and Nicole Poulton for assistance with the single-cell genomics protocol.

This paper was prepared by L.R.T. under award NA06OAR4320264 06111039 to the Northern Gulf Institute by NOAA's Office of Oceanic and Atmospheric Research, U.S. Department of Commerce. This work was supported by the NOAA Atlantic Oceanographic and Meteorological Laboratory 'Omics Program.