For the 454 shotgun runs, the first sequencing run generated 110,684 reads (~50 Mb of sequence information) and the other generated 82,726 reads (~35 Mb of sequence information). Both of these shotgun runs generated ~10× shotgun reads and the average read length is about 400 bp. For the 454 paired-end runs, the first sequencing run generated 116,008 reads (~38 Mb of sequence information) and the other generated 121,685 reads (~46 Mb of sequence information). Both of the two paired end runs generated ~10× paired end reads and the average read length is about 330 bp. The raw shotgun reads and paired-end reads were assembled into 123 contigs and 3 scaffolds. The N50 contig size is 145,039 bp and the largest scaffold contains 68 contigs contributing to a size of 5,911,648 bp, which shows that our raw assembly is highly contiguous. Over 99% of the total reads were assembled, resulting in coverage of around 23× for Klebsiella oxytoca HKOPL1. After PCR gap filling and Sanger sequencing of the PCR product, the complete genome of Klebsiella oxytoca HKOPL1 was achieved.

We hereby report that the genome of Klebsiella oxytoca HKOPL1 is a circular genome of 5.9 Mb with 55.92% GC content (Table 1). This data is now available in the NCBI GenBank (accession: CP004887, BioProject ID: PRJNA194061).

A total of 5,772 CDS sequences were predicted from the complete K. oxytoca HKOPL1 genome. About 95.8% (5,528 genes) of the total predicted genes were annotated against the NCBI non-redundant database (Table 1). Thirty-five tRNAs and sixty-two rRNAs were identified in the chromosome sequence. Most RNA elements were clustered to form 6 large RNA islands (Figure 1).

After comparative genome analysis, 1,391 genes were identified as core genes from 10 complete Klebsiella strains. Subsequently, a BMCMC (Bayesian Markov Chain Monte Carlo) phylogenomic tree was constructed based on these identified core genes (Figure 2). Our phylogenomic analysis highlighted that the Klebsiella oxytoca HKOPL1 strain is clustered with K. oxytoca strains 1686 and E718, and it is phylogenetically distant from the K. pneumoniae strains (Figure 2). It is also phylogenetically closer to K. oxytoca KCTC 1686 than K. oxytoca E718 (Figure 2).

K. oxytoca has a different classification from K. pneumoniae because it can grow on melezitose, but not on 3-hydroxybutyrate 10 ; K. oxytoca is also phylogenetically distant from K. pneumonia based on 1,391 core genes, indicating that K. oxytoca may have many other unique features. After comparative analysis between the three K. oxytoca strains and one representative K. pneumoniae strain (K. pneumoniae KCTC 2242), we put together a venn diagram showing 3,134 common genes shared among the four strains (36.7% out of 8,554 genes found in all four Klebsiella strains, see Figure 3). K. pneumoniae KCTC 2242 has more strain-specific genes and shares less common genes (<100) with the other K. oxytoca strains, while the three K. oxytoca strains share 1,343 common genes absent in the genome of K. pneumoniae KCTC 2242 (Figure 3).

The numbers of predicted genes in the chromosomes of all Klebsiella strains range from 4,923 to 5,772 (see Additional file 1); the percentage of paired orthologous genes in all Klebsiella strains shared with K. oxytoca HKOPL1 ranges from 37.85% to 86.22% (see Additional file 1). Additional file 1 also indicates that the percentage of paired orthologous genes in K. oxytoca strains 1686 and E718 (>80%) is much higher than that in both K. pneumonia and K. variicola strains (~40%). Additional file 2 suggests that all paired orthologous genes in Klebsiella strains are shared with K. oxytoca HKOPL1, and that their PAI-like regions and prophage-like regions also share similar positioning.

Both phylogenomic analysis and the percentage of common genes in all Klebsiella strains show that K. oxytoca HKOPL1 is closer to K. oxytoca strains 1686 and E718, rather than K. pneumoniae and K. variicola strains.

A total of 4,337 (~75.1%) out of 5,772 predicted genes of K. oxytoca HKOPL1 were identified in the NCBI COG database. Comparative study of single letter COG analysis for the three closely related K. oxytoca strains is shown in Figure 4, and the result for K. oxytoca HKOPL1 shows that the 9 COG classes ([G] (10.15%), [K] (8.58%), [S] (8.12%), [E] (8.65%), [R] (9.2%), [P] (6.13%), [C] (6.04%), [M] (4.61%) and [J] (4.2%)) have over 4% of the total annotated COG genes that comprise over half of the total annotated COG genes. The corresponding functions of these 9 COG classes are: carbohydrate transport and metabolism [G], transcription [K], function unknown [S], amino acid transport and metabolism [E], general function prediction only [R], inorganic ion transport and metabolism [P], energy production and conversion [C], cell wall/membrane/envelope biogenesis [M] and translation, ribosomal structure and biogenesis [J]. After the comparative COG study, the percentages of all COG classes were found to be similar within the three K. oxytoca strains (Figure 4).

Only 211 (~3.7%) out of 5,772 predicted genes of K. oxytoca HKOPL1 were identified in the VFDB database. Comparative study of virulence genes between K. oxytoca HKOPL1 and other Klebsiella strains shows that all Klebsiella strains have 3.5%-4.07% of virulence genes in their genomic sequences (Table 2); the location of all potential virulence genes in the K. oxytoca HKOPL1 genome is shown in Figure 1. Thirty-five Pathogenicity island-like regions were identified in the chromosome and their features are shown in Additional file 3.

K. oxytoca was reported as a causative organism of antibiotic-associated Hemorrhagic Colitis in human patients 13 . However, their pathogenicity in giant pandas is still unknown. In the K. oxytoca HKOPL1 genome, 35 pathogenicity island-like regions were identified, of which 14, PAI-2, PAI-5, PAI-7, PAI-8, PAI-9, PAI-11, PAI-14, PAI-18, PAI-20, PAI-21, PAI-22, PAI-24, PAI-27 and PAI-33, make up over 50% of PAI-virulence genes out of the total number of ORFs (see Figure 1 and Additional file 3).

A comparative CDS study (see Additional file 2) indicates that PAI-5 exists in all Klebsiella strains, while PAI-2, PAI-8, PAI-9, PAI-11, PAI-14, PAI-18, PAI-20, PAI-21, PAI-24 and PAI-33 only exist in K. oxytoca strains. PAI-7, PAI-22 and PAI-27 exist in all K. oxytoca strains, and some of these genes are partially shared among other Klebsiella strains. Although virulence genes and pathogenicity island-like regions were observed in the K. oxytoca HKOPL1 genome, over-interpretation of its pathogenicity should be avoided (especially as the strain was isolated from a healthy giant panda) and further targeted fragment knockout experiments are required for confirmation.

Only 23 (~0.4%) out of 5,772 predicted genes of K. oxytoca HKOPL1 were identified in the ARDB database and were matched to their corresponding resistant antibiotics. A comparative study of K. oxytoca HKOPL1 and two other closely related K. oxytoca strains shows that all three strains potentially have antibiotic resistance against gentamicin_b, acriflavin, chloramphenicol etc. plus multidrug resistance; for the full list, please refer to Additional file 4.

The analysis of potential antibiotic resistance shows that the three K. oxytoca strains have similar antibiotic resistance profiles with the exceptions of resistance to gentamincin_b, kasugamycin, trimethoprim, and dibekacin. Dibekacin resistance was found only in K. oxytoca E718, while resistance to kasugamycin was not found in K. oxytoca KCTC 1686. This comparative study provides detailed antibiotic resistance profiles for K. oxytoca strains that will be useful for future gene/gene cluster knockout experiments.

K. oxytoca was reported as a bacterial species resistant to several antibiotics (Högenauer 2006). Potential antibiotic resistant genes were identified in all three K. oxytoca complete genomes. The number of potential antibiotic genes and the antibiotic categories are similar in all three K. oxytoca strains. For K. oxytoca HKOPL1, 12 potential antibiotic resistant genes may come from plasmid and one potential antibiotic gene may come from prophage (data not shown); this information is important for studying antibiotic resistant genes existing in the gut microbiota of giant pandas. Future experiments are required to confirm all potential antibiotic resistance in K. oxytoca HKOPL1.

About 28% (1,615 out of 5,772) of the predicted genes have been identified in the ACLAME database. For K. oxytoca HKOPL1, 20 genes potentially come from viruses, 1,398 genes potentially come from plasmids and 197 genes potentially come from prophages. The total number of horizontal gene transferring (HGT) elements in K. oxytoca HKOPL1 is less than the other two closely related strains (Table 3). Twenty prophage associated genes were further curated to show their functions as well as their sources. Results indicate that all 20 prophage associated genes are from different phages infecting a wide range of bacteria genus including Aeromonas, Enterobacteria, Xanthomonas, Staphylococcus, Geobacillus, Enterococcus, Prochlorococcus, Vibrio, Mannheimia, Klebsiella and Salmonella (see Additional file 5). Some of the phage genes, such as contig00001_orf00136/contig00001_orf00138 (NrdD anaerobic NTP reductase large/small subunit) and contig00001_orf05069/contig00001_orf05192/contig00001_orf05193/contig00001_orf05194/contig00001_orf05202 (Gp22, Gp18, Gp18, Gp20 and Gp62), are clustered together. Such evidence suggests that these clustered genes were probably obtained from ancient phage integration events even though these regions lack complete phage structure (see Additional file 5).

HGT events were very common in bacteria species 14 and prophages, viruses and plasmids were supposed to be major vectors in HGT events 14 . For K. oxytoca HKOPL1, 28% of its genes were identified as potential HGT genes. The number of potential HGT genes in the other two K. oxytoca strains was similar. Most of the potential HGT genes of K. oxytoca HKOPL1 might have come from plasmids (~80%).

A region ranging from 4,612,823bp to 4,641,509bp was identified as an incomplete prophage region (see Additional file 2). The length of this region was 28.6kb with 39 CDS and it was similar to Enterobacterial phage mEp390. The CDS in this region was classified into different categories and illustrated in Additional file 6.

Horizontal gene transfer analysis suggests that HGT events are common in all K. oxytoca strains and various plasmids can be potential players in such events. In addition, 20 prophage associated genes were identified in K. oxytoca HKOPL1 from various phages infecting a wide range of different bacteria species. Some of them were clustered together, suggesting that they were obtained from ancient phage integration events. An incomplete prophage integrated region (~28.6 kb) was also identified in K. oxytoca HKOPL1 which might have been obtained from a previous prophage integration event.

Functional pathways of K. oxytoca HKOPL1 were annotated using KEGG Automatic Annotation Server. Among all these pathways, the starch and sucrose metabolism pathway (ko00500) and the butanoate metabolism pathway (ko00650) were present in the annotation result.

For the starch and sucrose metabolism pathway (ko00500), 32 genes were identified and two of them are related to cellulose degradation. They are Endoglucanase (K01179, EC:3.2.1.4) and Beta-glucosidase (K05349, EC:3.2.1.21). The former enzyme is able to transform cellulose to cellobiose, while the latter enzyme is able to transform cellobiose to β-D-Glucose. Further comparative analysis also shows that endoglucanase and beta-glucosidase are core genes existing in all 10 Klebsiella strains. This evidence reveals that Klebsiella may be able to degrade/utilize cellulose as an energy source and might have the ability to survive cellulose-rich environments such as inside the gut of giant pandas.

For the butanoate metabolism pathway (ko00650), 33 genes were identified and some of them were related to 2,3-Butanediol production. 2,3-Butanediol is widely used in plastics, solvents and antifreeze preparation. The 2,3-Butanediol production pathway is as follows: 1) Pyruvate dehydrogenase E1 component subunit alpha/beta (K00161 and K00162, EC:1.2.4.1) and Pyruvate dehydrogenase E1 component (K00163, EC:1.2.4.1) convert pyruvate to 2-(α-Hydroxyethyl)-ThPP, 2) Acetolactate synthase I/II/III large/small subunit (K01652 and K01653, EC:2.2.1.6) and acetolactate synthase II small subunit (K11258, EC:2.2.1.6) convert 2-(α-Hydroxyethyl)-ThPP to 2-Acetolactate, 3) Acetolactate decarboxylase (K01575, EC:4.1.1.5) converts 2-Acetolactate to (R)-2-Acetoin, and 4) (R,R)-butanediol dehydrogenase / diacetyl reductase (K03366, EC:1.1.1.4) convert (R)-2-Acetoin to 2,3-Butanediol. Further comparative analysis shows that all genes in the 2,3-Butanediol production pathway are core genes existing in all 10 Klebsiella strains. This evidence reveals that Klebsiella may be a bacterial genus that has industrial potential for 2,3-Butanediol production (e.g. a 2,3-Butanediol producing strain K. pneumonia KCTC 2242 7 ).

Pathway analysis identified a cellulose degradation pathway in K. oxytoca HKOPL1. Further comparative study shows that this cellulose degradation pathway is present in all 10 Klebsiella strains. This evidence reveals that Klebsiella may have the ability to degrade/utilize cellulose as an energy source and survive cellulose-rich environments, and it also provides insight into speculation that the gastrointestinal bacteria of giant pandas help to transform cellulose to β-D-Glucose, which is absorbed by the gastrointestinal cells as an energy source. Since giant pandas cannot degrade bamboo cellulose using its own digestive enzymes 2 , K. oxytoca HKOPL1 may play a key role in supplying energy to giant pandas. This is the first report with evidence of cellulose metabolism in a complete bacteria genome isolated from giant panda feces. In addition, a 2,3-Butanediol production pathway was identified in K. oxytoca HKOPL1 and found in all 10 Klebsiella strains as well. This discovery is of particular industrial interest because it means the Klebsiella genus could potentially produce 2,3-Butanediol in a biologically sustainable way.

To provide functional evidence of K. oxytoca HKOPL1’s ability to degrade cellulose, carboxymethyl cellulose (CMC) was used as a sole carbon source. After 72 hours of anaerobic incubation at 37°C, the appearance of a clear zone around the colony was observed (see Additional file 7). This result infers the presence of cellulolytic activity of isolated K. oxytoca HKOPL1 from the gut of giant pandas.

K. oxytoca HKOPL1 was isolated from a healthy adult giant panda living in the zoological institution of Hong Kong. All samples were obtained with permission and under memorandum through the Ocean Park Corporation, Hong Kong. First, fecal swabs were taken from inside fresh fecal masses using sterile plain swabs (Copan®#155C) and then streaked on different primary culture media. Afterwards, the agar plates were incubated aerobically at 37°C for 18–24 hours. Initial colonies were separated from the rest according to the colonial morphology on the primary plates, and then Gram staining was performed on these colonies. Subcultures of each colony from each agar medium were performed in the appropriate media (e.g. brain heart infusion agar plate - Oxoid® #CM1136) to get the pure culture. This was followed by initial rapid identification procedures (catalase, oxidase, indole tests etc.) and biochemical identification using Mini API® and API® systems (bioMerieux ®sa France).

An isolated single bacteria colony was grown in LB medium overnight to reach the exponential growth phase. Bacteria cells were then harvested by centrifugation before performing the bacteria genomic DNA extraction according to JGI’s “Bacterial genomic DNA isolation using CTAB” protocol 15 .

To confirm the purity of genomic DNA, the 16S rDNA specific region of K. oxytoca HKOPL1 was amplified and cloned. Then, 20 positive clones were submitted for Sanger sequencing. BLASTN analysis revealed that the K. oxytoca HKOPL1 16S rDNA sequences highly correlate with the current Klebsiella genus database. Evaluation of the quality of genomic DNA was done by using the Quant-iT™ PicoGreen dsDNA kit (Invitrogen).

A whole genome shotgun library of K. oxytoca HKOPL1 was generated with 0.5 μg of the genomic DNA; the shotgun sequencing was done with 454 GS Junior General Library Preparation Kit (Roche). In addition, an 8kb insert paired-end library was generated with 15 μg of the genomic DNA; the paired-end sequencing was done with 454 GS Junior Paired-end Library Preparation Kit (Roche). Paired-end reads were used to orientate the contigs into scaffolds.

The DNA libraries were amplified by emPCR and sequenced with FLX Titanium Sequencing Chemistry (Roche). Two shotgun runs and two paired-end runs were performed for each library. After sequencing, the raw data was assembled by Newbler 2.7 (Roche) with default parameters. Primer pairs were designed to amplify the gaps between contigs, and the resulting PCR products were directly sequenced by Sanger sequencer ABI 3130.

To predict the genes in the complete genome of K. oxytoca HKOPL1, Glimmer 3.02 16 , a prokaryotic gene prediction software, was used. First, all predicted CDS (Coding DNA Sequences) were translated into protein sequences by in-house Perl scripts. These protein sequences were then aligned against the NCBI non-redundant database (January 2013) with BLASTP 17 . Filtering of protein sequences with >90% alignment length and >40% identity was applied, and the description of the best hit was assigned to each corresponding predicted gene. Intergenic regions were annotated by RepeatMasker 18 with default parameters.

In order to identify the evolutionary position of K. oxytoca HKOPL1, phylogenomic analysis was applied to this strain together with other publicly available Klebsiella strains. Nine complete genomes of Klebsiella strains were downloaded from NCBI GenBank (see Additional file 8 for accession numbers). The orthologous genes were identified with BLAT (BLAST-like alignment tool, 19 ) under default parameters, by aligning the predicted genes of K. oxytoca HKOPL1 against all annotated genes of the nine Klebsiella strains.

The single copy genes with >90% of alignment length of the total Klebsiella strains were considered core genes. All the core genes were then aligned by MUSCLE (multiple sequence comparison by log-expectation) under default parameters 20 and concatenated together. Finally, the concatenated aligned genes were submitted to MrBayes 21 with the GTR + G + I substitution model for BMCMC (Bayesian Markov Chain Monte Carlo) phylogenomic tree construction. The chain length was set to 10,000,000 (1 sample per 1,000 generations). The first 2,000 samples were discarded as burn in after scrutinizing the trace files of two independent runs with Tracer v1.4 22 .

For comparative studies, the orthologous genes that were found to be shared between K. oxytoca HKOPL1 and other Klebsiella strains were plotted alongside all the predicted genes in K. oxytoca HKOPL1 in order to compare pathogenicity island-like regions and prophage-like regions in K. oxytoca HKOPL1. For comparative study within K. oxytoca species only, a Venn diagram was constructed based on the common genes found in the three K. oxytoca strains (K. oxytoca HKOPL1, K. oxytoca KCTC 1686 and K. oxytoca E718). Circular figures were generated with CIRCOS software 23 .

The COG (Clusters of Orthologous Groups of proteins) database 24 is a collection of functionally orthologous proteins from both prokaryotic and eukaryotic genome sequences. This database is often applied to functional gene classification of complete bacterial genome sequences.

In this study, the predicted ORF (Open Reading Frame) sequences of K. oxytoca HKOPL1 were translated into protein sequences by in-house Perl scripts. BLASTP 17 was subsequently applied to align these translated protein sequences against the COG database. Protein sequences with >90% alignment length and >20% identity were filtered, with the description of the best hit assigned to the corresponding predicted gene. Then, all predicted genes were classified into COG classes. The K. oxytoca HKOPL1 COG annotation was also compared to the other two K. oxytoca strains.

The Virulence Factor Database (VFDB) 25 is a collection of virulence factors (VFs) of various medically significant bacterial pathogens. In silico identification of virulence factors is helpful prior to any form of virulence-based animal experiment or gene knock-out experiment.

As in COG analysis, the predicted ORF sequences of K. oxytoca HKOPL1 were translated into protein sequences by in-house Perl scripts. BLASTP 17 was subsequently applied to align these protein sequences against the VFDB database. Protein sequences with >90% alignment length and >40% identity were filtered, and the description of the best hit was assigned to the corresponding predicted gene. K. oxytoca HKOPL1 virulence factors were compared to other publicly available Klebsiella strains.

Pathogenicity islands (PAIs) are mobile microbial genetic elements essential for the process of disease development; they are important research targets of microbial pathogenesis. They have been observed in evolutionary distant bacteria and may have been horizontally transferred among microbes. The pathogenicity island database (PAIDB) 26 is a comprehensive database of all reported pathogenicity islands (PAIs) and predicted/potential PAI regions. In silico identification of PAIs is useful prior to subsequent confirmation by gene knock-out experiments. In this study, all predicted genes of K. oxytoca HKOPL1 were preprocessed by in-house Perl scripts before submission to PAIDB for the identification of pathogenicity island-like regions.

Bacterial antibiotic resistance is often caused by bacteria acquiring various drug resistant genes. The Antibiotic Resistance Genes Database (ARDB) 27 is a summary of most antibiotic resistance genes from publicly available information. Each gene (including its resistance type) is annotated with its resistance profile and mechanism of action.

To analyze the drug resistant genes of K. oxytoca HKOPL1, in-house Perl scripts were used to translate predicted ORF sequences into protein sequences. These protein sequences were then submitted to BLASTP 17 to be aligned against the ARDB database. Protein sequences with >90% alignment length and >40% identity were filtered, and the description of the best hit is assigned to the corresponding predicted gene. Lastly, in-house Perl scripts were used again to summarize all annotated genes for their corresponding antibiotics. K. oxytoca HKOPL1’s drug resistant genes were compared to the other two K. oxytoca strains.

The major source of prokaryotic horizontal gene transferring (HGT) elements includes genes from phage genomes, plasmids, viruses and transposons. Classification of these genes can be found in the ACLAME (A CLAssification of Mobile genetic Elements) database 28 , which helps scientists to discover potential HGT genes from various sources. Psi-BLAST and HMM (Hidden Markov Model) methods were applied in a similarity search against the public database to classify HGT genes into different sources.

As before, analysis starts with translating predicted ORF sequences of K. oxytoca HKOPL1 into protein sequences by in-house Perl scripts. BLASTP 17 was then applied to align all these protein sequences against the ACLAME database. Protein sequences with >90% alignment length and >40% identity were filtered, and the description of the best hit is assigned to each corresponding predicted gene before matching all annotated genes to their corresponding horizontal transferring vectors. K. oxytoca HKOPL1’s potential horizontal transferring genes were again compared with the other two K. oxytoca strains.

In relation to prophage-like regions, phages are viruses which infect bacteria and they can be classified into two categories: lytic and temperate. Integrated phages are technically termed as prophages. Temperate phages are able to integrate into the plasmid or chromosome of their bacterial host. The PHAge Search Tool (PHAST) is a web server for identification of prophage existence (both complete and incomplete) 29 . The complete genome sequence of K. oxytoca HKOPL1 was submitted to PHAST for prophage-like region identification. Afterwards, the features of the prophage-like region were summarized.

Kyoto Encyclopedia of Genes and Genomes (KEGG) is a high level pathway analysis database which collects metabolic pathway information and integrate genomic, chemical and systemic functional information 30 . KEGG gene catalogs from genomes are associated with higher systemic cell functions. To analyze the metabolic pathway of K. oxytoca HKOPL1, the predicted ORF sequences were first translated into protein sequences by in-house Perl scripts. Then, these protein sequences were submitted to the KEGG database for automatic pathway annotation (http://www.genome.jp/kaas-bin/kaas_main). Afterwards, all annotated pathways were manually downloaded and curated by in-house Perl scripts.

Five microliters of pure culture grown overnight (OD600 ~ 1) was spotted on CMC (Carboxymethylcellulose) agar plates, which were prepared according to the procedure reported by Kasana 31 . Plates were incubated anaerobically at 37°C for 72 hours, followed by flooding with Gram’s Iodine solution for 3 minutes. The appearance of a clear zone around the colony is recognized for the bacteria’s ability to degrade CMC.

The data set supporting the results of this article is available in the NCBI BioProject repository (accession: CP004887, BioProject ID: PRJNA194061), unique persistent identifier and hyperlink to dataset in http://www.ncbi.nlm.nih.gov/bioproject/194061.