Surface (SUR) seawater samples (5 m deep, not pre-filtered, Sample ID: TARA_G000000266) were collected for single cell analysis from the North Indian Ocean during the circumnavigation expedition Tara Oceans, station TARA_039 ¹ (Figure 1). The environmental setting of the station was inferred from its physical report and can be found at Pangaea website in the following link: http://store.pangaea.de/Projects/TARA-OCEANS/Station_Reports/TARA_039_oceanographic_context_report.pdf.

Replicated 1 mL aliquots of seawater were cryopreserved with 6% glycine betaine (Sigma) and stored at -80°C. Samples used in this study were collected from TARA_039 (SUR; 5 m deep) located in the North Indian Ocean (18.59–66.62) on March 18th, 2010 (Supplementary Table 1). SAGs generation and preliminary 16S rRNA screening with primers 27F and 907R (Supplementary Table 2) were performed at the Bigelow Laboratory Single Cell Genomics ². More information can be found in Swan et al. (2011). Partial 16S rRNA gene sequences (640–850 bp) were compared to previously deposited sequences using the RDP Naive Bayesian rRNA Classifier Version 2.4 tool, and 98 SAGs with >99% similarity to the reference strain K. algicida OT-1 (AY195836) were selected for further analyses. These represented 84% of successfully amplified SAGs (total of 117 SAGs). They were named CFSAG39SUR referring to their taxonomic assignment (Cytophaga-Flavobacteria), retrieval method (SAG) and sampling station (TARA_039, SUR). Other SAGs used in this study (AAA285-F05 and AAA242-P21) were collected at 3,000 m at 4,800 m, respectively, in the North Pacific Sub-tropical Gyre. They belong to the Hawaii Ocean Time-Series (HOT, Figure 1A) and were generated following the same procedures.

Amplification of the 16S rRNA gene, ITS and 23S rRNA gene was done by a single Polymerase Chain Reaction (PCR) amplification with the 358F and CF434R primers (Supplementary Table 2). Amplification was performed in a Biorad Thermocycler by an initial denaturation at 94°C for 5 min, followed by 40 cycles of 94°C for 1 min, 55°C for 1 min and 72°C for 2 min, and a final extension at 72°C for 10 min. Each amplification reaction contained: 3 to 20 ng of template DNA, dNTPs (200 μM each), MgCl₂ (2 mM), primers (0.5 μM each), Taq DNA polymerase (1.25 U), the PCR buffer supplied by the manufacturer (Invitrogen, Paisley, United Kingdom) and MilliQ water up to the final volume of 25.

A set of newly designed primers was used for rpoB amplification (Supplementary Table 2). Each PCR reaction followed an initial denaturation at 94°C for 5 min, 40 cycles of 94°C for 1 min, 40°C for 45 s and 72°C for 1 min, and a final extension at 72°C for 5 min. Each amplification reaction contained: 1.5 to 10 ng of template DNA, dNTPs (0.2 mM each), MgCl₂ (1.5 mM), primers (0.25 μM each), KAPA2G Robust HotStart DNA Polymerase (1 U), KAPA2G Buffer A (KAPA BIOSYSTEMS, Wilmington, MA, United States) and MilliQ water up to the final volume of 25 μl. For proteorhodopsin (PR) amplification, the primers used and the PCR protocol were those described by Yoshizawa et al. (2012) (Supplementary Table 2) using the KAPA2G Robust HotStart DNA Polymerase and KAPA2G Buffer A from KAPA BIOSYSTEMS. Each PCR reaction followed an initial denaturation at 94°C for 5 min, 35 cycles of 94°C for 1 min, 44°C for 45 s and 72°C for 1 min, and a final extension at 72°C for 5 min. Each amplification reaction contained: 3 to 20 ng of template DNA, dNTPs (0.2 mM each), MgCl₂ (3 mM), primers (0.5 μM each), Taq DNA polymerase (1.25 U), the PCR buffer supplied by the manufacturer (Invitrogen, Paisley, United Kingdom) and MilliQ water up to 25 μl.

Polymerase Chain Reaction products were purified and sequenced by Genoscreen (Lille, France) with OneShot Sanger sequencing and aligned against NCBI’s nBLAST and xBLAST databases for identification. Reference sequences used in this study were retrieved from Integrated Microbial Genomes and Metagenomes (IMG) 4 Data Management or NCBI database. Using Geneious Pro 4.8.5 software all sequences were aligned (using default Geneious progressive pairwise aligner) trimmed and manually checked. The software was also used to obtain the complete 16S rRNA gene assembling the amplicon sequences from MDA screening and the phylogenetic markers amplification. Alignments were processed into phylogenetic trees with Mega 5.2.2 Software choosing Maximum Likelihood statistical methods along with 1000 bootstrap replications. The best-fit substitution models for Maximum Likelihood phylogenies of the studied markers were: GTR with Gamma distribution (16S rRNA gene), Kimura-2 with Gamma distribution (23S rRNA gene) and WAG with Gamma distribution (PR and RpoB). Representatives from the Proteobacteria and Cyanobacteria phyla were used as outgroups in all phylogenies. Phylogeny reconstruction of functional genes used the closest sequences to the SAGs’ genes obtained from the Ocean Microbiome Reference Gene Catalog (Sunagawa et al., 2015).

Nucleotide-Nucleotide BLAST 2.2.28+ was used to recruit metagenomic reads from several metagenomic samples similar to all the available sequenced Kordia genomes. Considering the aim of this part of the study, which was to determine the presence and distribution of the novel Kordia sp. CFSAG39SUR in different oceanic regions, depths and size fractions, a database containing all 4 Kordia genomes was generated. This ensured a competitive recruitment between the genomes for each metagenomic sample analyzed. The quality and genome completeness of the four reference genomes in the database was measured with software checkM (Parks et al., 2015) and fetchMG (Sunagawa et al., 2013) (Supplementary Table 3). The available genomic sequences of Kordia AAA285-F05 did not code for any of the marker genes or COGs used by the software, resulting in a genome completeness estimation of 0% despite being 283.58 Kb long. The other three reference genomes (K. algicida, K. jejudonensis, and K. zhangzhouensis) were estimated to be complete and free of contamination. Recruitment regarded only one read per target gene (-max_target_seqs 1), an identity percentage higher than 70% (-perc_identity 70), as well as an e-value lower than 0.000001 (-e-value 0.000001). All other nBLAST settings were set on default. In order to avoid random alignments, a filtering process was applied using R software (R Core Team, 2013) (i) excluding alignments with coverage lower than 90% of the sequence length, (ii) removing duplicated reads within each genome and (iii) masking reads belonging to the ribosomal operon. This operon does not follow the species delimitating recruitment patterns as the rest of the genome does and its presence would overestimate read recruitment (Caro-Quintero and Konstantinidis, 2012). All contigs from each genome were concatenated to a single sequence and the file containing the four sequences (one for each genome) was used to generate the database. This step was done using NCBI’s makeblastdb application. As query we selected those metagenomic samples which contained the most metagenomic reads 16S rRNA and metagenomic Illumina tags (miTags) (Logares et al., 2013) annotated as Kordia (data not shown) using SILVA database Release 115 (Quast et al., 2013). These were stations TARA_039 (0.2–20 μm) and TARA_085 (0.2–3 μm) from Tara Oceans (Supplementary Table 1 and Figure 1A) (Sunagawa et al., 2015). In addition, we used 4 metagenomes from the deep waters of Malaspina 2010 circumnavegation expedition (Acinas et al., in preparation), two of them from the Brazil Basin in the Atlantic Ocean and two of them from Circumpolar Deep Waters (Supplementary Table 1 and Figure 1A). Due to the flexibility of Bacteroidetes’ lifestyle, the metagenomic selection was done from a limited pool of deep ocean metagenomes covering the free-living fraction of prokaryotes (0.2–0.8 μm) and the fraction of prokaryotes attached to particles or aggregates (0.8–20 μm) of a same seawater sample. The coverage of different deep ocean biomes and the highest abundance of Kordia related metagenomic read recruitments was also taken into account for the metagenomes’ selection for the study. All queries consisted on merged, paired-end reads with different sequencing depths. Merging was done using Mothur v.1.33.3 (Schloss et al., 2009). Recruited data normalization was done: (i) by reference genome size (resulting in recruited reads/genomic bp), then (ii) upscaling this ratio to 1 Mb (as done in Swan et al., 2011), and finally (iii) by metagenomic sequencing depth. As sequencing depths varied between metagenomes, normalization was performed to the smallest metagenomic sequencing depth of the study.

Accession numbers for PCR products of Kordia CFSAG39SUR: 16S rRNA gene (MF187452), ITS (MF187454), 23S rRNA gene (MF187453), rpoB (MF187456), and proteorhodopsin (MF187455). Accession numbers for partial 16S rRNA gene sequences are: MF187458 for AAA285-F05 and MF187457 for AAA242-P21.

The genomes used as reference for the FRA are the following: Bacteroidetes bacterium SCGC AAA285-F05 (JGI GoldStamp Id Go0060979), K. algicida strain OT-1 (NCBI RefSeq NZ_ABIB00000000.1), K. jejudonensis strain SSK3-3 (NCBI RefSeq NZ_LBMG00000000.1), K. zhangzouensis strain MCCC 1A00726 (NCBI RefSeq NZ_LBMH00000000.1). The Tara Oceans metagenomes used as query are the following (INSDC run accession numbers):TARA_039_DCM_0.22-1.6 (ERR599145), TARA_039_MES_0.22-1.6 (ERR599037| ERR599172), TARA_085_SRF_0.22-3 (ERR599090| ERR599176), TARA_085_DCM_0.22-3 (ERR599104| ERR599121), TARA_085_MES_0.22-3 (ERR599008| ERR599125). The Expedición Malaspina metagenomes are the following (JGI Genome Portal, genome.jgi.doe.gov/): MP1493 (/DeeseametaMP1493_FD), MP1494 (/DeeseametaMP1494_FD), MP0326 (/DeeseametaMP0326_FD), MP0327 (/DeeseametaMP0327_FD) with the previous agreement of the PI of the JGI CSP 602 grant.

On February 2010, 1 month prior to sampling, a moderate bloom of phytoplankton occurred in the Oman Gulf due to a coastal upwelling under strong wind conditions. It propagated to the central basin of the Arabian Sea reaching only up to station TARA_039 (see Figure 3 in Roullier et al., 2014), which was the last mesotrophic station occupied by the Tara Oceans expedition before entering the warm oligotrophic waters of the central basin. Previous stations occupied during this leg in the North Indian Ocean (TARA_037 and TARA_038) were in the cooler waters of the upwelling, whereas station TARA_039 was at a mesoscale frontal region with the oligotrophic zone (Figure 1B). Particles produced during the phytoplankton bloom decreased in size as the cruise reached the central basin (see Figure 10 in Roullier et al., 2014). Diatom contribution to total phytoplankton decreased from 20 to 5% from stations close to the Gulf of Oman down to station TARA_40, outside the upwelling waters. Sunagawa et al. (2015) found that the bacterial assemblage at station TARA_039 was dominated by Proteobacteria in both surface and DCM (Deep Chlorophyll Maximum) layers, followed by Cyanobacteria and other, less abundant groups including Euryarchaeota, Actinobacteria, Deferribacteres, and Bacteroidetes (Sunagawa et al., 2015). In the mesopelagic (MES) zone there was a decrease in the abundance of Proteobacteria, a small increase of Deferribacteres and Actinobacteria and a remarkable increase in Thaumarchaeota and unclassified Bacteria. There was an increase in Bacteroidetes abundance toward the end of the Tara Oceans’ sampling leg in the North Indian Ocean (Sunagawa et al., 2015).

Amplification of the three phylogenetic markers of the ribosomal operon was successful in 78 out of 98 Kordia SAGs, and their alignment indicated a 100% identity among them (Table 1). The resulting sequences contained full 16S rRNA gene (1,495 bp), ITS-1 (536 bp) and partial 23S rRNA gene (415 bp) (Table 1). The alignment of the resulting full 16S rRNA gene sequence against NCBI database assigned the SAGs to the genus Kordia, family Flavobacteriaceae. Its best hit with a nucleotide identity of 99% and sequence coverage of 97% was an uncultured Bacteroidetes sampled in the Puerto Rico Trench (clone PRTBB8540, acc. HM798967) at 6,000 m depth. The closest sequenced cultured relatives were K. algicida strain OT-1 with a 96.8% of nucleotide identity and 99% coverage (complete sequence, acc. NR027568) and K. jejudonensis strain SSK3-3 with the same identity percentage but 96% coverage (partial sequence, NR_126287). A 100% nucleotide identity was found with partial 16S rRNA gene sequences from 17 SAGs collected at 3,000 m (AAA285) and one SAG from 4,800 m depth in the North Pacific Sub-tropical Gyre (AAA242-P21), which belong to the Hawaii Ocean Time-series (HOT) (Stepanauskas, unpublished results). The phylogenetic reconstruction based on the representative sequences of the 16S rRNA genes from the 78 SAGs (Figure 2) strongly indicated that they belonged to the genus Kordia (phylum Bacteroidetes, family Flavobacteriaceae). High bootstrap values supported the cluster formed by the 78 CFSAG39 SAGs with the 17 AAA285-SAGs, the one AAA242-P21 SAG and the uncultured Bacteroidetes clone PRTBB8540. This cluster’s sequences were all retrieved from deep waters except for the SAGs reported here. The resulting sequence from the 78 identical partial 23S rRNA genes, when aligned against the NCBI database, showed a lower identity against K. algicida OT-1 (90.5%). Nevertheless, it was still its closest culture hit and the phylogenic tree supported the assignment of the SAGs to the genus Kordia ITS with high bootstrap values (Supplementary Figure 1).

The 78 SAGs’ ITSs were identical and coded for Ala-tRNA and Ile-tRNA. These tRNAs showed only one polymorphism each when compared to those of K. algicida. The nucleotidic change was located at the next-to-last base, not affecting the structure of the functional molecule. The box A conserved region present at the end of the ITS was identical to that of K. algicida although the former was located slightly upstream (17 bp) (Supplementary Figure 2).

Amplification of the beta subunit of RNA polymerase (rpoB) was successful in 18 SAGs, obtaining ∼1,400 bp amplicons with 100% of both amino acid and nucleotide identity among them. Their closest hit in NCBI database was the rpoB gene of K. algicida OT-1 with a 98% amino acid identity. High bootstrap values supported the location of the SAGs’ rpoB sequence in the phylogenetic tree, within the Flavobacteriaceae cluster and closest to K. algicida OT-1 (Supplementary Figure 3).

An amplification summary can be found in Table 1.

The PR gene was amplified in 34 out of the 98 Kordia SAGs. The resulting 167 amino acid sequences turned out to be 100% identical to one another. The sequence was closest to the Flavobacteriaceae Flagellimonas sp. DIK-ALG-169’s proteorhodopsin (89% identity, 100% coverage) (AHN13811) and Kordia sp. PC-4 (Yoshizawa et al., 2012) (88% identity, 91% coverage) in the NCBI database. When comparing sequences in the Ocean Microbiome Reference Gene Catalog (OMRGC), the closest hit shared 84.2% amino acid similarity (unclassified Flavobacteria, NOG136807 ID.v1.022398661). Phylogenetic reconstruction placed the SAGs’ proteorhodopsin in a cluster of Flavobacteriaceae PR genes (Figure 3). We tried to amplify the proteorhodopsin gene from the 17 SAGs retrieved from 3,000 m at the North Pacific Gyre (AAA285-F05) but the results were negative.

The SAGs’ proteorhodopsin alignment with other closely related PR sequences showed that this proteorhodopsin had the same structural features as its relatives. It had the predicted transmembrane helices, key amino acids for functionality and the location of the spectral tuning amino acid, which in this case was methionine (Figure 4), indicating absorbance of light spectrum between 518 and 535 nm (green light).

Metagenomic reads from different stations were recruited through a competitive FRA with the available sequenced genomes of the genus Kordia: the deep-sea Kordia SAG AAA285-F05 (283.58 Kbp), the seawater surface isolated in culture K. algicida OT-1 (5.01 Mbp), the freshwater K. zhangzhouensis MCCC 1A00726 (4.03 Mbp) and K. jejudonensis SSK3-3 (5.3 Mbp), isolated from a region where spring freshwater and seawater meet. We counted reads with 70–100% identity to these reference genomes.

For Kordia AAA285-F05, the number of recruited reads increased with depth in all stations tested. The percentage of reads recruited per genomic Mbp was several orders of magnitude larger in the free-living than in the particle associated size fractions (Figure 5; numeric data in Supplementary Table 4).

Highest relative abundances of reads per Mbp with 95–100% identity for Kordia AAA285-F05’s occurred in the free-living fraction of bathypelagic waters from the Brazil Basin (MP0327, 0.51%) and the Circumpolar Deep Waters of the Pacific Ocean (MP0326, 0.46%). Recruitment percentages for this SAG in these two locations decreased down to 0.04 and 0.01%, respectively. In TARA_039 highest recruitment % per Mbp was 0.0007% at 270 m, while in the Southern Ocean (TARA_085), recruitment increased with depth reaching 0.06% at 790 m. Recruitment was virtually zero for the isolates genomes at this recruitment identity percentage.

Recruitment reads with 70–94.9% identity might be indicative of populations related to, but different from, the reference genome present at the tested metagenomic sample. The recruitment trend observed in the higher identity percentage for Kordia SAG AAA285-F05 was consistent at the lower percentage as well. Recruitment % per genomic Mbp increased at all depths of both TARA_039 and TARA_085, now reaching values of 0.05% at 270 m in TARA_039 and the maximal of 0.38% at 790 m in TARA_085. Similar abundances were obtained for the Brazil Basin (MP0327, 4,001 m) and the Circumpolar Deep Waters (MP1494, 2,150 m) in their free-living fraction (0.3 and 0.27%, respectively). Recruitment decreased in the larger size fraction of the two Malaspina 2010 circumnavigation expedition metagenomes down to 0.04 and 0.01%, respectively. For the other three Kordia genomes competing for recruitment in the same metagenomic samples, recruitment values decreased significantly to values ranging from 0.008% per Mb for K. zhangzhouensis in TARA_085 surface waters to values down to 10^-7% per Mb. The recruitment for these three genomes followed a similar pattern, with highest values in TARA_085 surface decreasing with depth. In their case, bathypelagic samples and North Indian ocean’s DCM (25 m) recruited similar values (∼0.0002%). The lowest values were found at TARA_039 mesopelagic waters (0.00001%).

FRA plots (Figure 6) showed a similar pattern of Kordia AAA285-F05 recruitment for the meso- and bathypelagic metagenomes (TARA_085_MES, MP0327, MP1494) of the free-living prokaryotic fraction. There was a very good coverage of the whole reference sequence, with the exception of some fragments where there was a visible decrease in the recruitment at all identity percentages. The highest read densities were found above 95% identity. In the TARA_039 OMZ recruitment plot, an overall decrease in read density was observed, identities ranged from 85 to 99%. The two genomic regions with very low recruitment were also apparent in this plot.

The three reference Kordia genomes (Supplementary Figure 4) showed read clouds at the lower identity percentage in all cases. There was a good coverage of the length of the genome but not as intense as that observed in Figure 6 for reference genome Kordia AAA285-F05.