INTRODUCTION

As part of ongoing protist genome projects, we have generated large-scale genomic data from several recently discovered free-living heterotrophic flagellates, including the cryptist Palipitomonas bilix (Yabuki et al., 2010; 2014), the undescribed mantamonad strain SRT-306 (Glücksman et al., 2011; Brown et al., 2018), the phylogenetically deep-branching discobid Ophirina amphinema (Yabuki et al., 2018), and the marine malawimonad strain Imasa heleensis (Heiss et al., 2020). All of these protists have been recognized in the last decade as representing deeply branching lineages amongst all eukaryotes (Yabuki et al., 2010; Glücksman et al., 2011; Brown et al., 2018; Yabuki et al., 2018), and their genomic and cell structural data potentially have significant implications for better understanding early eukaryotic evolution, including the nature of the last eukaryotic common ancestor (LECA). Moreover, these microbes are phagotrophic nanoflagellates, an ecological group that plays an important role as grazers of bacteria (Collier and Rest, 2019). As each protist species investigated here is an obligate heterotroph, bacteria, presumably coisolated together with the protist species, are consumed by the protist as a necessary food source. However, the identity and community composition of the cocultured bacteria have been unknown.

Despite efforts to reduce bacterial load in DNA extraction, such as by physical separation using polycarbonate membrane filters before harvesting, the bacteria could not be completely removed, and their genomic sequences became a major source of contamination in the resulting data, often accounting for the majority of the generated reads. Therefore, the bacterial sequences had to be identified and separated from the eukaryotic reads during the process of protist genome assembly and annotation.

Studies in other systems have shown that there is value in nontarget sequencing data. For example, genomic and transcriptomic datasets generated from animal tissues can contain sequences originating from parasites and/or symbionts of those organisms. Classifying and identifying these sequences can shed light on the distribution of particular pathogens in animal populations (Lopes et al., 2017; Galen et al., 2020). In eukaryote genome projects like ours, these bacterial data, which could be valuable to the microbial research community, are often discarded without being assembled, annotated, and archived at a public data repository. However, these data may contain new genomes and/or species, and depending on sampling, culturing, and DNA extraction conditions, they could reveal population dynamics of the sampled habitat or the culture system. Here, we analyzed such “contaminating” bacterial data: from these four protist cultures, we assembled, annotated, and identified 28 scaffolds representing complete or near-complete bacterial genomes. Our analyses show that the majority of these bacterial scaffolds represent species that are not currently represented in GenBank's bacterial genome databases.

MATERIALS AND METHODS

Culture conditions: The culture strain of Palpitomonas bilix investigated in this study was established from a driftwood sample from Marcharchar Island, Palau, collected in June 2006 (Yabuki et al., 2010). The malawimonad Imasa heleensis and the discobid Ophirina amphinema were collected in September 2013 and March 2016, respectively, from a shallow lagoon on Nusa Lavata in the Hele island chain, in the Western Province of the Solomon Islands (Yabuki et al., 2018; Heiss et al. 2020). Mantamonad strain SRT-306 came from the scraped surface sample of a barracuda that was caught in a lagoon in Iriomote-jima, Okinawa Prefecture, Japan, in September 2013. All of the cultures grew well and were maintained in Erd-Schreiber medium (ESM: UTEX) or ESM fortified with 2.5% Cerophyl medium (ATCC 802) at 23° C, with the exception of SRT-306, which was maintained at 16° C.

DNA extraction and sequencing: For I. heleensis and O. amphinema, DNA was extracted on multiple occasions between January 2018 and May 2019. Cell cultures were scraped using a sterile cell scraper to lift adherent cells, and protist cells were collected on 0.6 µm and 0.4 µm Millipore polycarbonate filters, respectively, under partial vacuum, to reduce bacterial load. Cells were incubated on these filters for 3 hours at 56° C in lysis buffer from the Qiagen MagAttract HMW DNA Kit (Qiagen, Hilden, Germany). DNA was extracted from the lysates using this same kit according to the manufacturer's instructions. An aliquot of DNA (1 µg) from each culture was sent to the New York Genome Center for Illumina Paired End 2x150 bp sequencing on the Illumina HiSeqX platform (Illumina, San Diego, CA). Additional DNA (∼1–4 µg) was prepared for sequencing on the Oxford Nanopore platform using either the SQK-LSK108 or SKK-LSK109 Genomic DNA by Ligation kit (Oxford Nanopore Technologies, Oxford, UK) according to manufacturer's instructions. Libraries were sequenced on the MinION platform using FLO-MIN106 SpotON R9 Flow Cells.

For P. bilix and mantamonad strain SRT-306, DNA was extracted using the Qiagen DNeasy Blood & Tissue Kit, following the manufacturer's instructions. Cells were collected onto 0.8 µm filters, as specified above, and washed three times each with 5–10 ml artificial seawater to reduce the bacterial load prior to DNA extraction. The plates on which SRT-306 were grown were scraped to lift adherent cells; no such procedure was required for the free-swimming P. bilix. DNA samples were sent to Cornell Sequencing Core and New York Genome Center for Illumina Nextera library preparations and Paired End 2x150 bp sequencing on the on the Illumina HiSeq2500 platform.

Base calling and assembly: Raw MinION FastQ files were base-called using Guppy v1.1 (Oxford Nanopore Technologies). Hybrid assemblies using both MinION and Illumina data from I. heleensis and O. amphinema were assembled using MaSuRCA v3.2.6 (Zimin et al., 2013). Illumina reads generated from cultures of P. bilix and SRT-306 were assembled using ALLPATHS-LG (Gnerre et al., 2011; Ribeiro et al., 2012).

Identification of bacterial genomes from metagenomic assemblies: The resulting assembled scaffolds and contigs were individually analyzed for completeness using Benchmarking Universal Single-Copy Orthologs (BUSCO) v3 (Simão et al., 2015) with the “Bacteria odb9” dataset. The presence of 5S, 16S, and 23S rRNA genes was searched for using rnammer v1.2 (Lagesen et al., 2007). The scaffolds both (a) with a BUSCO completeness ≥80% and (b) that contained at least one complete rRNA operon were retained as “complete” bacterial scaffolds. The 80% completeness threshold was chosen because there was a sharp cutoff in completeness for scaffolds below 80%, most of which were less than 60% complete and/or did not contain complete rRNA operons. In an attempt to circularize the complete scaffolds, two independent approaches were taken: MUMmer v3.0 (Kurtz et al., 2004) was used to look for overlaps at the ends of each scaffold, and Circlator (Hunt et al., 2015) was used to attempt circularization of each scaffold assembly.

TABLE 1.

Features of the four metagenomic assemblies.

Annotation and classification of bacterial genomes: Bacterial scaffolds were annotated for gene content using the rapid prokaryotic genome annotation tool, Prokka (Seemann, 2014). The scaffolds were classified into taxonomic groups using the Genome Taxonomy Database Tool Kit (GTDB-TK) version 0.3.1 (Chaumeil et al., 2019). Scaffolds were compared to each other using the online ANI calculator tool from the EZ BioCloud database (Yoon et al., 2017). The scaffolds are available on NCBI under the BioProject accession code PRJNA619388.

Abundance estimates: Relative abundance estimates were calculated for each complete genome by aligning the Illumina Paired End reads to the pooled set of complete bacterial genomes isolated from the same culture using the Burrows-Wheeler Aligner (BWA) (Li and Durbin, 2009; Durbin, 2010), retaining only the uniquely mapping read pairs. Read pairs mapping to each scaffold were calculated using SAMtools version 1.9 (Li et al., 2009; Li, 2011). Final calculations and plots were made in R version 3.3.3 and ggplot2 version 3.3.2 (Wickham, 2016).

RESULTS

Identification of bacterial genomes from metagenomic assemblies: We screened mixed assemblies from four mono-isolate cultures of marine protists that were sequenced for the purpose of assembling high-quality reference genomes for the protist species (table 1). We identified 28 scaffolds that represent complete or near-complete bacterial genomes: four from SRT-306, three from P. bilix, six from O. amphinema, and 15 from I. heleensis (fig. 1, table 2). We attempted to circularize all genomes by looking for overlaps at the ends of the scaffolds by aligning each scaffold to itself and plotting the overlaps using MUMmer (Kurtz et al., 2004), and by running each scaffold through the Circlator assembly circularization pipeline (Hunt et al., 2015). MUMmer found no overlapping regions in any scaffold, and Circlator did not circularize any of the 28 scaffolds (data not shown). All but one of the 28 scaffolds contained at least one complete prokaryotic rRNA operon (i.e., with full-length 5S, 16S, and 23S rRNA; table 2). For scaffold CP051235, only one incomplete rRNA operon was annotated, which contained full-length 16S and 23S rRNA genes separated by two tRNA loci, but no downstream 5S rRNA. The 16S and 23S rRNA genes are located midscaffold on CP051235, and genes are annotated both upstream of the 16S and downstream of the 23S loci; thus, the lack of an annotated 5S rRNA is not due to an assembly truncation.

Features and gene content of each bacterial scaffold: We used a rapid annotation tool, Prokka (Seemann, 2014), to annotate each bacterial scaffold for gene content and other basic genomic features (table 3). The 28 bacterial scaffolds ranged in length from 3.1 to 6.0 Mb, contained 2873 to 5666 annotated genes, and had an average of 44.5 tRNAs, congruent with what is typical of bacterial genomes (Land et al., 2015). GC content of these genomes ranged from 37.4 to 65.8%, within the range noted for bacterial genomes (Almpanis et al., 2018).

Taxonomic classification of the 28 scaffolds: We used the GTDB-TK tool to classify each genome taxonomically (Chaumeil et al., 2019). From the 28 scaffolds we found three phyla, five classes, eight orders, 10 families, and 21 genera represented (table 3). Eight of the scaffolds were classified to the species level. Two of these genomes represent the same taxon, Alcanivorax sp. DSM 26293 (GenBank Assembly Accession GCF_900107995.1), which is present in data from the cultures of O. amphinema and I. heleensis. These two scaffolds are close to identical at the sequence level, with an average nucleotide identity (ANI) of 99.99%. For the remaining 20 scaffolds, 19 were classified to the genus level and one (CP051235) was classified only to the family level. As such, these 20 scaffolds could represent either novel or known but unsequenced bacterial species. GTDB-TK establishes taxonomy using a combination of the placement of the query into a reference tree, the queries' Relative Evolutionary Distance (RED), and its ANI percent values compared to reference genomes (Matsen et al., 2010; Parks et al., 2018; Chaumeil et al., 2019). Queries placed on terminal branches of the reference tree (i.e., are sister to individual known species or strains) are placed in the genus of their sister taxon. A query is classified to species level if it is placed in an existing genus and its ANI value is within the circumscription radius (typically, 95%) of the closest species; otherwise the query is classified as a novel species in that genus (Chaumeil et al., 2019). For example, scaffold CP051248 had an ANI value of 99.94% and was classified as Pseudooceanicola atlanticus. In contrast, scaffold CP051241 was placed under the genus Marinobacter. Its closest ANI value, to the species Marinobacter similis, was 89.97%, suggesting that it is a novel species of Marinobacter. For queries that are placed on internal branches, taxonomy is instead determined by placement and/or RED (Parks et al., 2018; Chaumeil et al., 2019). Only one of our queries, scaffold CP051235 (which was the scaffold lacking the 5S rRNA gene), was not placed in an existing genus. This scaffold was placed in the family Hyphomicrobiaceae, and potentially represents a new genus therein.

FIG. 1.

Light microscopy images showing the four protist cells mentioned in this study. Arrow indicates protist cell, bac = bacteria. Scale bar = 5 µm. A. Palpitomonas bilix, B. mantamonad strain SRT-306; C. Ophirina amphinema; and D. Imasa heleensis.

TABLE 2.

Features of each bacterial genome. Genomes are organized by protist culture of origin. Asterisks indicate different assemblies of the same bacterial strain.

FIG. 2.

Count of scaffolds in each metagenome that are ≥80% complete or more, as determined by BUSCO.

Relative abundance of bacterial species in each culture: We next evaluated the relative abundance of the 28 bacterial scaffolds by calculating the number of Illumina read pairs that aligned uniquely to each genome from the four sequencing libraries (fig. 2). Two of the metagenomes, from O. amphinema and P. bilix cultures, were each dominated by a single bacterial taxon, Marinobacter and Nitratireductor respectively, that accounted for ∼90% of the bacterial data from each library. The metagenomic data from the I. heleensis and SRT-306 cultures were not dominated by any one species (fig. 2).

DISCUSSION

In this study, we investigated bacterial data that are often discarded in the course of sequencing genomes from xenic protist cultures. From four metagenomic protist genome assemblies, we identified 28 assembly scaffolds that represent complete or near-complete bacterial chromosomes.

We were able to classify eight of the bacterial scaffolds to species level, 19 to genus level, and one to family level. Only two of these 28 bacterial scaffolds represented the same species, Alcanivorax sp. DSM 26293, which was present in the O. amphinema and I. heleensis cell culture metagenomes, making a total of 27 unique bacterial scaffolds from our dataset. Ophirina amphinema and I. heleensis were originally isolated from the same lagoon, so this Alcanivorax bacterium is possibly from this original source.

TABLE 3.

Taxonomic classification of bacterial genomes. Results given by GTDB-TK database. “NA” indicates scaffolds for which no closest hit was given. Asterisks (*) indicate placeholder family and genus names provided by GTDB-TK when no family and/or genus name is present in NCBI database.

Continued

A literature search was done of the top taxonomic hits for each genome to ascertain the typical environmental ranges of the source organisms. Salegentibacter is a genus found in hypersaline lakes, on the surfaces of marine fauna, and in marine sediments (McCammon and Bowman, 2000; Bowman, 2016). Hyphomonas (Abraham, 2020), Epibacterium (Wirth and Whitman, 2020), Maritimibacter (Lee et al., 2007), Aestuariibacter (Yi et al., 2004), Thalassobius (Arahal et al., 2005), Muricauda (Bruns et al., 2001; Bruns and Berthe-Corti, 2015), Nitratireductor (Singh et al. 2012), Balnenola (Urios et al., 2006), Oceanicola (Cho and Giovannoni, 2004), Pseudoceanicola (Lai et al., 2015), Planctomycetaceae (Ward, 2015), and Balneolaceae (Xia et al., 2016) represent groups of widely dispersed marine bacteria isolated from seawater and/or marine sediments or surfaces. Alcanivorax (Yakimov et al., 1998; Golyshin et al., 2015), Marinobacter (Gauthier et al., 1992; Bowman and McMeekin, 2015), Hyphomicrobiaceae (Kesy et al. 2019), and Filomicrobium (Schlesner, 1987; 2015) are marine groups previously isolated from seawater and seawater sediments enriched with crude oil and/or hydrocarbons. Nitratireductor is a nitrate-reducing genus found in various marine habitats (Labbe et al., 2004), members of which have previously been isolated from diatom cultures (Jang et al., 2011). The family Cryomorphaceae is present within a wide range of nonextreme ecosystems, both marine and terrestrial (Bowman, 2015). Lastly, Pararhizobium is a nitrogen-fixing genus that associates with plant roots, which has a worldwide distribution (Mousavi et al., 2015).

In summary, all but one of the described groups (Pararhizobium) are characteristically found in marine habitats, consistent with the fact that the cultures are from marine environments. These cocultured bacteria therefore are likely from the source samples from which each protist was isolated. New microbial genomic data generated from little-sampled habitats has the potential to be highly valuable to the research community, as studies have shown that marine and other infrequently sampled niches may contain more novel microbial diversity than habitats that are much closer at hand, like the human gut or local soils (Bech et al., 2020; Tessler et al., 2017). However, we cannot rule out that they are environmental contaminants, considering that most of the taxa have global marine distributions, the seawater culture medium would select for marine species, and our lab in Manhattan is surrounded by seawater and bacteria are commonly dispersed through the air and on surfaces (Mayol, 2017).

FIG. 3.

Relative abundance of each bacterial scaffold in Illumina sequencing libraries, represented as proportion of length normalized read pairs for the portion of each sequencing library mapping to the bacterial scaffolds. The two Alcanivorax sp. DSM 26293 scaffolds are marked with asterisks (*).

When we calculated the relative abundance of each bacterial scaffold across the four cultures, we found that the species are unevenly represented in the sequencing libraries, with both Ophirina amphinema and Palpitomonas bilix having one dominant bacterial taxon, Marinobacter and Nitratireductor respectively, while the novel malawimonad and mantamonad strains each had a more even distribution of bacterial taxa (fig. 2). This distribution may represent population patterns in the source cultures possibly caused by varying bacterial diversity in the original samples or by differing grazing preferences of each isolated protist in culture. However, each culture was preprocessed using a polycarbonate membrane filter prior to DNA extraction, which may have influenced the distribution of each species, by retaining both larger bacterial cells and bacterial flocs. Therefore, the bacterial abundances presented in figure 2 are not necessarily representative of the bacterial community of each culture.

Importantly, 20 of the 28 scaffolds we identified represent taxa that are not yet present in bacterial genome databases. As mentioned above, these are data that we usually discard from our protist genomics projects, but we find that they can also be a source of novel bacterial genomes and species. We note as well that we obtained these unidentified species from long-term cultures grown under constant environmental conditions and containing relatively few different organisms. This suggests that establishing axenic cultures of at least some of these novel bacteria might be practicable, thereby allowing for their formal description. However, the existence of axenic cultures may not remain the only requirement for formal description. A recent consensus statement (Murray, 2020) suggests mechanisms for the establishment of new prokaryote taxa based on genomic data. Should this consensus be formally accepted, all of our new taxa should be candidates for formal description. In either case, we find that the prokaryotic genomes generated during eukaryote sequencing projects need not be disposable byproducts. Instead, we suggest they are a key source for potentially important new discoveries.