Probability estimations of genomic features

Genomic signatures have been shown to display a species-specific pattern [13–16] and have been applied to bin sequences from metagenomic datasets. The most widely used genomic signature is tetranucleotide frequencies, which has been used in a number of metagenomic studies [4–8]. To distinguish whether two sequences are sampled from the same species based on their tetranucleotide frequencies, we downloaded 3,181 bacterial and archaeal genomes from the IMG website (http://img.jgi.doe.gov/), simulated genomic sequence fragments, and calculated Euclidean distance between the extracted tetranucleotide frequencies of the two sequences. The lengths of the simulated sequences were randomly chosen between 1,000 bps to 1,000,000 bps. The simulation was performed one million times for intra-genome (sequences sampled from the same genome) and one million times for inter-genome (sequences sampled from different genomes) comparisons. The histogram of intra-genome and inter-genome simulations, as shown in Additional file 1: Figure S1(A) and S1(B), revealed a large difference, in which intra-genome distances were grouped below 0.2 while inter-genome distances were more evenly distributed between 0.02 and 0.1. We estimated both the mean and variance values separately for both intra-genome and inter-genome Euclidean distances. The mean and standard deviation values for intra- and inter-genome distances were 0.015, 0.010 and 0.068, 0.034, respectively. We must note that both distributions of intra- and inter-genome distances were rejected by the Shapiro-Wilk test to be normally distributed (W value is 0.81 and 0.93 for intra- and inter-genome distances, respectively); however, because the difference between the histograms of intra- and inter-genome distances is very clear, we still applied their mean and standard deviation values to test whether the measured Euclidean distances are between sequences sampled from the same genome or between different genomes. We discuss this circumstance and potential improvement in the Discussion section.

in which $\mathit{G}\left({\mathit{S}}_2\right)$ represents the genome that sequence S₂ belongs to, $\mathit{dist}\left({\mathit{S}}_1,{\mathit{S}}_2\right)$ is the distance between the extracted tetranucleotide frequencies of S₁ and S₂ using the Euclidean distance function, N_intra and N_inter are the Gaussian distributions with estimated intra- and inter-genome mean distance values (μ _intra and μ _inter ) and distance variance values $\left({\mathit{\sigma }}_{\mathit{intra}}^2\mathrm{and}{\mathit{\sigma }}_{\mathit{inter}}^2\right)$, respectively (note that P _dist (S₁ ∈ G(S₂)) = P _dist (S₂ ∈ G(S₁)) since dist(S₁, S₂) = dist(S₂, S₁). The distribution P _dist for Euclidean distance is shown in Additional file 1: Figure S1(D). One can easily observe from the figure that the lower the distance, the more probable two sequences are sampled from the same genome.

where cov(S) indicates the coverage of sequence S, and Poisson (cov(S₁) | cov(S₂)) is the Poisson probability density function given mean value λ = cov (S ₂ ).

Expectation-maximization algorithm

The EM algorithm calculates the probability that a given scaffold belongs to any genome at the same time. The maximum number of iterations of the EM algorithm was determined by running MaxBin on simulated datasets with different maximum iteration numbers and measuring the performances of the binning results by precision and sensitivity (see below). As shown in (Additional file 1: Figure S2), the precision and sensitivity of the two simulated datasets were very stable for all maximum iteration number settings. We have set the maximum number to perform EM algorithm to 50 in case some larger metagenomes need more iterations to achieve reasonable binning results.

After the EM algorithm is finished, the scaffolds are assigned to the bin with the highest probability as long as the probability values surpass the minimum probability threshold, which is set to 80%. Sequences that do not meet the threshold are discarded as ‘unclassified.’

Before applying the EM algorithm on any metagenomic datasets, sequences shorter than minimum length threshold are removed since these sequences are likely to produce skewed tetranucleotide frequencies, which confuse the binning algorithm and erroneously classify shorter sequences into wrong bins. To find the best minimum length cutoff threshold, we performed binning using different length cutoff settings on simulated datasets. The result demonstrated that MaxBin achieves the best sensitivity and comparable precision with 1,000 bps cutoff setting [see Additional file 1: Figure S3]. Therefore, we set the minimum length threshold to 1,000 bps to achieve best performances.

Initialization of the algorithm

A common practice for an expectation-maximization algorithm is to randomize all parameters. This initial condition permits the possibility of converging parameters into local maxima. At the same time, the number of bins, which is one of the most crucial parameters, is usually unknown. We employed the single-copy marker genes to estimate the number of bins and initialize all required parameters. Genes from the scaffolds were predicted using FragGeneScan [19], and HMMER3 [20] was used (with --cut_tc option) to scan the predicted genes for 107 single-copy marker genes that are conserved in 95% of all sequenced bacteria, which has been used to determine the genome completeness of bins [7, 21]. After filtering out all mapped genes that do not meet the coverage threshold, which is 40%, the median number of scaffolds containing each of the marker genes is identified, considering that some marker genes may be fragmented into several pieces that may distort the estimation of the number of bins. The shortest marker gene that corresponds to the median number of bins is selected, and the tetranucleotide frequencies and read coverages of the scaffolds harboring this shortest marker gene are extracted to generate the initialization parameters for the algorithm. The reason we use the shortest marker gene for initialization is that shorter genes are less likely to be split between two scaffolds.

Recursive classification of bins

Despite careful selection of initialization conditions, the EM algorithm sometimes may still group scaffolds from several composite genomes into one bin. To alleviate this problem, all bins are recursively checked for the median number of marker genes. If the median number of marker genes of any bin is at least 2, the bin will be treated as a dataset waiting to be binned, and the whole EM algorithm will be applied to split the bin. All bins (including those created by reapplying EM algorithm on bins) will be checked for the number of marker genes until no bins can be split further.

Estimation of genome completeness

After the EM process, MaxBin scans the 107 marker genes in all bins for genome completeness, which is measured as the fraction of unique marker genes versus all marker genes. Not all bacteria are equipped with all 107 marker genes, such as the reconstructed genomes from the uncultivated TM7 lineage, in which representatives of this candidate phylum have only 100 out of all 107 marker genes [7]. However, the metric of 107 marker genes is maintained as the standard in the absence of specific phylum level counts of conserved marker genes.

Simulation of the test metagenomes

Simulated metagenomes with 10 species and 100 species were generated by MetaSim [22] and assembled by Velvet assembler 1.2.07 [23] with K = 55 and coverage cutoff set to 1. The detailed genome simulation settings can be found in Additional file 1: Table S1 and S2. For 10-species metagenomes, we sampled 5 million and 20 million paired-end reads referred to as the 20X and 80X datasets, respectively. For 100-species metagenomes, the settings used by [24] were mimicked, and 100 million paired-end reads were sampled to create three datasets: simLC+, simMC+, and simHC+. The 80-bps error model was downloaded from the MetaSim website (http://ab.inf.uni-tuebingen.de/software/metasim/errormodel-80bp.mconf/view) and used in simulating all metagenomes.

Performance evaluation of the binning results of simulated datasets

in which R _ij indicate the number of sequences (in terms of base pairs) that belong to genome j appears in bin i. If M > N, the majority of sequences in each bin likely belong to a single genome and the precision will be high; however, the sensitivity will be low since some genomes are represented by more than one bin. On the other hand if M < N, the sensitivity will tend to be high while precision is likely to be low. We also note that scaffolds lower than 1,000 bps, which is our minimum length cutoff for MaxBin algorithm, are not included in the unclassified sequences since these sequences cannot be binned and will be discarded in the first step of the binning algorithm.

Besides precision and sensitivity, we also evaluated the amount of correctly binned and mis-assigned sequences for individual species in the simulated datasets. We evaluated the number of sequences of each species that appear in each bin, and assigned the bins to the species with the highest amount of sequences in the bins. If two or more bins are assigned to the same species, only the bin with the greatest amount of sequences belonging to that species is kept, and only sequences belonging to that species in that bin are regarded as correctly binned; all sequences appear in other bins and sequences that do not belong to the assigned species in any bin are regarded as mis-assigned. Unclassified sequences are defined as sequences that pass the minimum length threshold but are not classified into any bins. Sequences lower than the minimum length threshold are not considered here since these sequences are not treated as part of the binning result.

Test environment

MaxBin was tested on a Linux operation system with 128G memory space and 16 AMD Opteron™ CPU cores at 2.2 GHz. With the exception of HMMER3, which MaxBin utilized to extract marker gene information, MaxBin itself and other affiliated software are not multi-threaded. The running time of MaxBin for all simulated and real datasets is reported in Additional file 1: Table S8. Note that MaxBin did not consume more than 1GB of memory space for all our test datasets except when we map reads against scaffolds to obtain scaffold coverage information, suggesting that the MaxBin algorithm could be executed on a personal computer.

Sequencing the enriched cellulolytic compost samples

Green waste compost samples were obtained from the City of Berkeley, CA, on 30 June 2012. Replicate 50 mL cultures (37A and 37B) containing MES-buffered M9TE minimal media [25] were established with microcrystalline cellulose (500 mg, 1% v/v) (Sigma-Aldrich) as the carbon substrate. These cultures were incubated on rotary shakers for two weeks at 37°C and 200 rpm. Five milliliters of this culture (10%) was transferred to a second replicate set of microcrystalline cellulose-containing cultures and incubated for an additional two weeks. After the second passage, DNA was extracted from culture biomass using previously described methods [26]. DNA fragments for Illumina sequencing were created using the Joint Genome Institute standard library generation protocols for Illumina HiSeq 2000 platforms. Illumina sequencing was performed on a HiSeq 2000 system. The DNA fragments were assembled using the Joint Genome Institute assembly pipeline. Raw reads were trimmed using a minimum quality cutoff of Q10. Trimmed, paired-end Illumina reads were assembled using SOAPdenovo v1.05 (http://soap.genomics.org.cn/soapdenovo.html) with default settings (-d 1 and -R) at different Kmer sizes (85, 89, 93, 97, 101 and 105 respectively). Contigs generated by each assembly (a total of six contig sets from the six Kmer sizes), were merged using in-house Perl scripts as following. Contigs were first de-replicated and sorted into two pools based on length. Contigs <1,800 bps were assembled using Newbler (Life Technologies, Carlsbad, CA, USA) to generate larger contigs (-tr, -rip, -mi 98, -ml 80). All assembled contigs >1,800 bps, as well as the contigs generated from the Newbler assembly, were combined and merged using minimus2 (-D MINID = 98 -D OVERLAP = 80) (AMOS: http://sourceforge.net/projects/amos). The average fold coverage (or read depth) of each scaffold was estimated by mapping all Illumina reads back to the final assembly using BWA (version 1.2.2) [27].

Extraction of glycoside hydrolase genes from the Sorangium sp. bin

Glycoside hydrolase (GH) genes of the Sorangium sp. binned genomes and the two Sorangium cellulosum strains were extracted using dbCAN [28], which was based on the protein families generated by CAZy database (http://www.cazy.org/) [29]. The extracted GH families were then grouped based on [30], in which the GH families were classified into cellulases (GH5, 6, 7, 9, 44, 45, 48), endohemicellulases (GH8, 10, 11, 12, 26, 28, 53), cell wall elongation enzymes (GH16, 17, 74, 81), de-branching enzymes (GH51, 54, 62, 67, 78), and oligosaccharide-degrading enzymes (GH1, 2, 3, 29, 35, 38, 39, 42, 43, 52). A new protein class, ‘lignin-degradation enzymes,’ was defined by including the protein families from AA1 to AA8 as suggested by [31].

Construction of phylogenetic trees

The 16S ribosomal RNA genes were extracted and aligned using MUSCLE [32]. The alignments were refined using Gblocks [33] and loaded into MEGA5 [34] to construct the maximum-likelihood tree using default settings (Tamura-Nei model, uniform rates, and complete deletion) with 1,000 bootstraps. The marker gene trees were built using the 35 marker genes previously reported [35]. The 35 marker genes were extracted from the downloaded or binned genomes, translated to amino acid sequences, and aligned by MUSCLE separately. The alignments were then concatenated and refined using Gblocks. MEGA5 was again used to build the maximum-likelihood marker gene tree using default settings (JTT model, uniform rate, and complete deletion) with 1,000 replicates.

Identifying clusters of orthologous groups families

Multiple alignments of all clusters of orthologous groups (COGs) were downloaded from eggNOG website (http://eggnog.embl.de/) [36] and converted to hidden Markov models using HMMER3 [20]. The functional categories of all COGs were also downloaded from eggNOG website for counting the numbers of genes for each functional category.

Pathway mapping

Proteins extracted from Sorangium sp. were searched for their KEGG (Kyoto Encyclopedia of Genes and Genomes) Orthology (KO) numbers using KEGG2 KAAS web service [37]. The resulting KO numbers were inputted into the Search&Color Pathway web service (http://www.genome.jp/kegg/tool/map_pathway2.html) available on the KEGG2 website [38] to identify the associated pathways.

Data access

The MaxBin program is available at https://sourceforge.net/projects/maxbin/. Metagenomic data, including raw sequencing reads and assembled sequences, for the enriched cellulolytic compost consortia are available at JGI IMG/M website (https://img.jgi.doe.gov/cgi-bin/m/main.cgi) under JGI taxon id 3300000869 (37A) and 3300001258 (37B). The MetaSim setting files, assembled scaffolds, and coverage files for replicating the simulation results, the binning results of HMP datasets that were mentioned in the text, and the binning results of the enriched cellulolytic compost metagenomes can be downloaded from the MaxBin download page (http://downloads.jbei.org/data/MaxBin.html).

Testing MaxBin on simulated metagenomes

MaxBin has been designed as an automated metagenomic binning software, which allows binning of assembled metagenomic scaffolds after the assembly of metagenomic sequencing reads with minimal human intervention. MaxBin was initially tested by binning several simulated metagenomic datasets produced by MetaSim [22] to evaluate its effectiveness. MaxBin was applied to two simulated metagenomes containing 10 species with different overall sequencing coverage (20X versus 80X average coverage). The species used in this simulation and their relative abundance ratios, defined by the actual coverage levels divided by summed coverage levels of all genomes, were summarized in Additional file 1: Table S1. MaxBin was first interrogated for its ability to classify sequences correctly into corresponding bins. MaxBin successfully binned the 80X metagenome into 10 bins, in which the majority of sequences were correctly classified (Figure 2(A)). High abundance genomes were binned almost perfectly; most of the erroneously binned sequences occurred in low abundance genomic bins with similar coverage levels. The precision was estimated to be 96.9%, as shown in Table 1. These results demonstrated the ability of MaxBin to estimate correctly the number of bins as well as utilize tetranucleotide frequencies and scaffold coverage information to bin most of the sequences accurately.

The metagenome with 20X average coverage was binned into three bins, each consisting of sequences from the most abundant three genomes (Figure 2(B)). Due to the lower sequencing coverages of the seven low abundance genomes, the majority of assembled scaffolds that belong to these low abundance genomes did not pass the minimum length threshold and cannot be binned by MaxBin. Therefore, only three bins were produced from the 20X metagenome, each belonging to the three genomes with highest abundances. Most of the scaffolds belonging to the seven low abundance genomes were discarded because these scaffolds did not pass the length threshold. For most binning methods, 1,000 bps is the minimum length to bin scaffolds successfully (see Discussion for details). Despite these limitations, MaxBin still binned the three most abundant genomes with precision measured at 97.01%. We note that for both 80X and 20X samples, the unclassified sequences were both lower than 1 Mbps (Table 1; sequences below minimum length threshold were not included), suggesting that MaxBin was able to classify most of the sequences while committing very few mis-assigned errors.

The binning capabilities of MaxBin were further tested on more complex metagenomic datasets. Genome coverage settings of three previously published simulated metagenomic datasets were generated (simLC, simMC, and simHC) [24], and three datasets were produced with similar settings: simLC+, simMC+, and simHC + [see Additional file 1: Table S2]. MaxBin isolated 57 bins from simLC+, 11 bins from simMC+, and 78 bins from simHC+. The precision, sensitivity, and the amount of unclassified sequences are also summarized in Table 1. Overall, sequences from high abundance genomes were binned more accurately and more comprehensively than low abundance genomes as demonstrated for the simLC + and simMC + datasets (Figure 3(A)-(B)). The first few genomes with high abundance levels were binned very accurately; the majority of mis-assigned scaffolds occurred in low-abundance bins. Investigation of the assembled sequences of simLC + and simMC + demonstrated that the correctly assigned scaffolds were significantly longer than mis-assigned ones [see Additional file 1: Figure S4]. This phenomenon accounts for the high accuracy of binning of the high abundance genomes, which generally consisted of longer scaffolds that can be binned more accurately. This effect can also be observed in the simMC + dataset, in which only a small portion of scaffolds from low abundance genomes passed the MaxBin minimum length threshold, and therefore only the high abundance genomes were binned. Due to the poor assembly quality of simMC + (the N50 for simMC + is only 383 while the N50 for simLC + and simHC + are 1293 and 17169, respectively), only 11 bins were generated from the simMC + dataset compared to 57 and 78 bins from the simLC + and simHC + datasets.For simHC + dataset, which had an evenly distributed species abundance levels, MaxBin yielded 78 bins (Figure 3(C)) with 88.93% precision. In the simHC + dataset, there were fewer mis-assigned scaffolds than in simLC + (22.24 Mbps in simHC+; 53.24 Mbps in simLC+). The appearance of unclassified sequences in some bins, which was measured at 62.5 Mbps, was probably due to the similarity of abundance levels between different species, which resulted in a lowered sensitivity of 73.46%.

Binning of datasets from the Human Microbiome Project

The Human Microbiome Project (HMP) was designed to document the microbial populations that occupy habitats in or on the human body. A total of 749 samples were generated for selected corporeal habitats, including the gut, the mouth, the vagina, and the skin [39]. MaxBin was applied to identify genomic bins from three selected HMP datasets with very different sample sizes: the tongue dorsum sample SRS013705, the subgingival plaque sample SRS014477, and the stool sample SRS018656. The reads and assembled sequences were downloaded from HMPDACC website (http://hmpdacc.org/). These three samples were chosen to compare the performance of MaxBin on real metagenomic datasets with very different amounts of raw sequence reads (12 GB for SRS013705, 1.4 GB for SRS014477, and 6.4 GB for SRS018656). MaxBin yielded 31 bins for SRS013705, 4 bins for SRS014477, and 10 bins for SRS018656, consistent with the total base pairs of scaffolds that passed the minimum length threshold (83 MB for SRS013705, 14 MB for SRS014477, and 45 MB for SRS018656). The number of unclassified sequences for SRS013705, SRS014477, and SRS018656 was 11.18, 0.37, and 2.17 Mbps. The taxonomy of the bins was analyzed using MEGAN4 [40] and compared to the results of the HMP community profiles (also downloaded from HMPDACC website) generated by fragment recruitment of the metagenomic sequencing data to related sequenced isolates [41]. The binning results and the closest species are listed in Additional file 1: Table S3. In general, MaxBin successfully generated bins for the high abundance populations in each sample. For instance, three of the four most abundant predicted species in SRS013705 were successfully binned, including populations related to Prevotella melaninogenica, Streptococcus salivarius, and Veillonella dispar. Similarly, in SRS014477, species among the four recovered genomes included Treponema denticola, Treponema vincentii, Corynebacterium matruchotii, and Bacteroidetes oral taxon 274, which were identified as the most abundant species in this sample. In the SRS018656 dataset, three of the four most abundant species were isolated as well, including Dialister invisus, Bacteroides cellulosilyticus, and Ruminococcus sp., [see Additional file 1: Table S3]. The MaxBin threshold on the minimum scaffold lengths limited the number of isolated bins that could be extracted from these datasets, including several strains of Rothua dentocariosa in SRS014477 and Fusobacterium sp. in SRS013705. Nevertheless MaxBin was still capable of isolating most of the high abundance genomes from real metagenomic samples with variable sequencing depths and different species composition complexity.

We compared our binning results to an approach using emergent self-organizing maps (ESOM), [2, 8, 9]. An ESOM graph was generated for the tongue dorsum sample SRS013705 based on tetranucleotide frequencies [2], and the sequences on the graph were colored based on the MaxBin binning results [see Additional file 1: Figure S5]. Since ESOM relied on the contour boundaries to distinguish binned genomes, we focused on observing whether the bins generated by MaxBin were consistent with the original ESOM boundaries. Indeed, comparison between ESOM [see Additional file 1: Figure S5(A)] and MaxBin binning results [see Additional file 1: Figure S5(B)] demonstrated that most MaxBin bins (shown as different colors) overlaid with the ESOM boundaries, confirming the ability of MaxBin to separate genomes according to their genomic signatures. In addition we also observed that closely related populations, such as the three bins representing three Prevotella species and two bins belonging to the Leptotrichia species, were not successfully resolved using the ESOM approach, suggesting that MaxBin could further distinguish genomes with similar tetranucleotide frequencies based on scaffold coverage levels and single-copy marker genes.

Recovering genomes from cellulolytic consortia using MaxBin

We found that even though the two replicates were enriched from the same compost inoculum, the species composition and the relative abundance ratios of these species diverged (Figure 4). In the 37A dataset, the most abundant bin was classified as Sorangium sp. followed by bins classified as Niastella sp. and Opitutus sp. In 37B, the most abundant species was classified as Niastella sp., which was a nearly identical bin to the one found in 37A, followed by Teredinibacter sp. and Sphingomonas sp. The Sorangium sp. bin, which occupied nearly half of the population in 37A (48.1%), was only found at 2.5% abundance in 37B. Furthermore, the second and third most abundant bins in 37B (Teredinibacter sp. and Sphingomonas sp.) were not observed in 37A. Note that the second most abundant species in 37B, Teredinibacter sp., is distantly related to Teredinibacter turnerae (with amino acid identity at 67.4%), an endosymbiotic cellulolytic gammaproteobacteria isolated from the gill tissue of a shipworm, Lyrodus pedicellatus[44].

The binning results obtained for MaxBin were compared with two other binning methods: ESOM (tetranucleotide frequencies) and differential coverage binning [7]. The ESOM graphs demonstrated that the MaxBin binning results fit extremely well with the ESOM boundaries for both replicates [see Additional file 1: Figure S6]. Differential coverage binning was performed on the Sorangium sp., Niastella sp. and Opitutus sp. bins since initial inspection of the metagenomic datasets revealed that the assembled sequences for these populations were nearly identical in 37A and 37B.

The sequences that were found only in MaxBin-derived bins compared to the differential coverage binning method were collected and analyzed by MEGAN to identify the taxonomy of these sequences. MEGAN analysis of the scaffolds not classified by differential coverage binning demonstrated that most of the scaffolds missed by the differential coverage binning approach were classified into the Proteobacteria lineage, as indicated in Additional file 1: Figure S8. Since the bins closest in abundance were affiliated with Bacteroidetes and Verrucomicrobia, the sequences affiliated with the Proteobacteria lineage shown in Additional file 1: Figure S8 likely belong to the Sorangium bin. Therefore, MaxBin collected a more complete Sorangium sp. genome compared to the differential coverage binning approach. We manually inspected the five scaffolds grouped into the Bacteroidetes lineage and found that the coverage levels of these scaffolds (974.7, 1990.5, 1019.9, 1027.4, and 1709.5) were much higher than predicted coverage of second abundant species, Niastella sp. (270.36). Since the coverage levels of these five scaffolds were much closer to that of Sorangium sp., their assignment may be consistent with belonging to Sorangium sp. bin. An expanded genome size obtained by MaxBin compared to differential binning was also observed for the Niastella and the Opitutus spp. genomes. In both cases, the majority of the additional sequences recovered by MaxBin were affiliated with the predicted lineage [see Additional file 1: Figure S8].

Sorangium sp. bin represents the genome of an unusual myxobacterium

A complete 16S rRNA gene (90% identical to S. cellulosum) was recovered from the bin containing the Sorangium sp. genome and a phylogenetic tree was constructed to classify the bin (Figure 5(A)). Analysis of the phylogenetic tree demonstrated that the novel myxobacterial population was a Deltaproteobacterium in the Myxococcales order affiliated with the suborder Sorangiineae, but was distinct from the family Polyangiaceae, which contains the validated species Sorangium cellulosum, Byssovorax cruenta and Chondromyces apiculatus[46]. This new family in the Sorangiineae has no cultivated members and consists of 16S rRNA clones representing uncultivated species. The two 16S rRNA clones in this family that are most similar (99% identity) to that of Sorangium sp. were recovered separately from earthworm guts and large-discharge carbonate springs [Genbank: HM459718 and KC358117]. The phylogenetic classification of this bin was confirmed by construction of a concatenated gene tree from the genomic bin with 35 single-copy marker genes, which confirmed that it was distantly related to Sorangium cellulosum (Figure 5(B)). Surprisingly, the MaxBin binning results, supported by complementary binning by ESOM and differential coverage binning methods, demonstrated that the Sorangium sp. genome was approximately 5 MB, while the two sequenced strains of Sorangium cellulosum have genomes of 13.0 MB (strain So ce56) and 14.7 MB (strain So0157-2). Genomes of 11 myxobacterial genomes were compared, and 193 genes were identified as universally shared. For those 193 genes, 158 genes were found to be present in Sorangium sp., suggesting that despite its significantly smaller size, this genome still contains most of the common genes found in myxobacteria.

We compared the gene content between the three Sorangium species to further understand their differences. COG families of all extracted genes from the three Sorangium species were identified and compared. The identified COGs were classified into different functional categories for the three species, as depicted in Figure 6 (the actual numbers can be referred to in Additional file 1: Table S6). Since the genome sizes of the two Sorangium cellulosum genomes were about 2.5 times larger than the recovered Sorangium sp. genome, gene numbers for individual functional categories were expected to be approximately 2.5 times larger for S. cellulosum isolates. Indeed, we observed that the number of genes in COG categories for the two Sorangium cellulosum species were, on average, 2.35 times more than those in Sorangium sp. (7,535 and 7,255 genes in all COG categories for the two Sorangium cellulosum strains compared to 3,155 genes for Sorangium sp.). Among the COG categories, gene numbers differed most for COG group Q (Secondary metabolites biosynthesis, transport, and catabolism); the two Sorangium cellulosum isolate genomes averaged 4.17 times more assigned genes than the recovered Sorangium sp. genome. Production of large numbers of secondary metabolites is a characteristic of myxobacterial metabolism, and Sorangium cellulosum isolates in particular produce antifungal and antibacterial compounds [47]. The only other myxobacterial isolate with a similarly sized genome to the uncultivated Sorangium sp. genome is Anaeromyxobacter dehalogens (5.01 MB) [48]. A. dehalogens is characterized by an extensive genomic repertoire for anaerobic respiration, including growth with nitrate, halogenated organics and metals as terminal electron acceptors. The uncultivated Sorangium sp. genome lacked genes for denitrification, sulfate reduction, and metal reduction, indicating that it may primarily have an aerobic lifestyle [see Additional file 1: Figure S9].

In contrast to the general survey provided by COG family comparisons, extraction of glycoside hydrolase (GH) and auxiliary activity (AA) enzymes using the CAZy database demonstrated that the number of genes relevant to biomass deconstruction was, on average, doubled in the uncultivated Sorangium sp (Figure 7). In particular, genes that cluster with glycoside hydrolase families specifically involved in cellulose hydrolysis (GH5, 6, 7, 9, 44, 45, 48) were twice as abundant (29 versus 13 and 16) in the uncultivated Sorangium sp. genome compared to the S. cellulosum isolate genomes (the actual number of all GH genes are listed in Additional file 1: Table S7).