Genomic sequences used in this study were downloaded from the NCBI FTP site (ftp://ftp.ncbi.nlm.nih.gov/genomes/all). The versions that were downloaded are also available online at ftp://www.homd.org/publication_data/20160425. A summary of all the meta information for each genome is available in the Excel file PG_Genome_Summary.xlsx in the above FTP folder. This file lists all the detail information that are provided by NCBI, such as methods for sequencing, assembling and annotation, as well as various IDs for the same genome including GenBank Accession, GenBank Assembly Accession, Refseq Accession, Refseq Assembly Accession. Table 1 lists the basic information and sources of the sequence data of the 19 P. gingivalis genomes analyzed in this report.

To facilitate further comparison and future studies, all the data and results generated in this study, including the original downloaded sequences, annotations, the comparative results presented in this paper, as well as additional complete results that were not mentioned or discussed, are available for download from this FTP data repository site: ftp://www.homd.org/publication_data/20160425.

The first P. gingivalis genome released was that of the strain W83 in 2003 and the latest one was released in February 2016. Of the 19 genomes, eight were assembled into a single contig and were considered complete and finished genomes; the remaining were released as various numbers of sequence contigs assembled from whole genome shotgun (WGS) sequence reads. The sequence of JCVI SC001 appears to have a 1-contig circular sequence under the Genbank Accession number CM001843, however it is a pseudo-contig generated by ordering the 284 unassembled contigs (accession number APMB01000000) based on the homologous matches to the genome of TDC60 (McLean et al., 2013) and joining the ordered contigs with 282 100-N spacer sequences (total N length is 28,200 bps). Thus, it is not considered a complete or finished genome. Examining the sequences for the presence of Ns reveals the “completeness” of the genomes. Table 2 shows the reported length, non-N length, total number of Ns and the distribution of the N fragments in the genomic sequences. Overall strain A7A1-28 is the smallest of the completed P. gingivalis genomes with a size of 2,249,024 bps. HG66 has the largest size of all the sequenced P. gingivalis genomes at 2,441,680 bps after removing the 100 Ns placed at the end of the sequence. The placement of the 100 Ns at the end of the sequence was due to the unsuccessful attempt to circularize the sequence with the minimus2 software used by the PacBio sequencer at default settings (personal communication). For this reason the HG66 genome should not be considered complete. Almost all the unfinished draft genomes consist of various numbers of Ns ranging from 698 Ns in SDJ2 to 7,200 Ns in F0569 (Table 2). It is likely that some of these published contigs were assembled based on a reference genome and the Ns had been filled in the gaps. Hence the true order of genes identified by the annotation process may not be correct.

Table 3 gives a numeric summary of the genome annotation results by the NCBI Prokaryotic Genome Annotation Pipeline (released 2013, http://www.ncbi.nlm.nih.gov/genome/annotation_prok/). The NCBI pipeline is capable of identifying more than just the protein-coding genes, rRNAs and tRNAs, including several interesting types of genes such as binding sites, repeat sequences, pseudo-genes, and several types of non-coding RNAs (ncRNAs). However, since the NCBI pipeline is quite new, more features are still being added and since some of the annotations of these P. gingivalis genomes were done prior to 2013, the annotation results may not be comprehensive until the annotation is updated again based on the latest NCBI pipeline.

In addition to the NCBI annotations, RAST (Rapid Annotations using Subsystems Technology) is also a popular pipeline for annotating microbial genomes (Aziz et al., 2008). All the 19 P. gingivalis genomes, as well as the chosen outgroup P. asaccharolytica DSM20707 were subjected to the RAST pipeline and the results were compared with those done by the NCBI pipeline. As shown in Table 4, both the RAST and NCBI pipelines identified almost the same number of rRNA and tRNA genes. However, the numbers of protein-coding genes varied quite significantly between the two pipelines. Although most of the genes were commonly identified, up to hundreds of protein-coding sequences can be missed by either system. Moreover, 86% (6,422 of 7,382 for all the 19 genomes) of these uniquely identified genes code for hypothetical proteins and 80% are shorter than 100 amino acids in length (in fact, only 94 have lengths ≥ 500 amino acids), thus the impact due to the annotation discrepancy may not be as significant especially when drawing conclusions in genome-wide systematic analysis or metabolic pathway capability.

A list of the 960 (7,382–6,422) non-hypothetical proteins is provided at the link (ftp://www.homd.org/publication_data/20160425/2_Summary_of_Genome_Annotations/Non-overlap_Non-hypothetical_protein_identified_by_NCBI_or_RAST.fasta).

The 16S rRNA sequences have been used to infer the evolutionary relatedness of the prokaryotes due to its slow rate of evolution (Woese et al., 1990). However, multiple rRNA genes including 16S rRNAs are common in prokaryotic genomes (Klappenbach et al., 2000) and the genomic copy number of 16S rRNA varies greatly among species from 1 to 15 (Vetrovsky and Baldrian, 2013). The number of rRNA genes was reported to correlate with the rate at which phylogenetically diverse bacteria respond to resource availability (Klappenbach et al., 2000). As shown in Table 4, all of the eight genomes which had been assembled to a single contig contain four copies of 5S, 16S, and 23S rRNA genes respectively, thus it is reasonable to believe that all P. gingivalis genomes have four copies of the rRNA operons. The lower number of rRNA genes in the unfinished genomes is likely due to the incompleteness of the sequences and is also likely due to the fact that genomes sequenced by short reads sequencing platforms such as those of the Illumina sequencers cannot be easily assembled across the repeated regions such as the highly conserved rRNA operons.

The 16S rRNA sequences of all the 19 genomes annotated by NCBI were extracted and aligned for the construction of a phylogenetic tree. Based on the annotation, there are a total of 24 unique 16S rRNA gene sequences identified from the 19 genomes (Table 5, first column), plus the sequence of a close species P. asaccharolytica strain DSM 20707 (Accession Number CP002689), making it a total of 25 unique sequences in the study. However, many of the sequence differences are due to different annotated lengths. After aligning all the 24 (25 if including the outgroup sequence) unique sequences and trimming off the leading and trailing sequences not present in all copies (trimmed aligned length = 1,425 bps), the aligned portion of several sequences are identical and the number of unique sequences was reduced to 20 (second column of Table 5). Strains 381, A7A1-28, ATCC 33277, and W83 all have four copies of identical sequences in their genomes (Table 5, last column). The aligned regions of all the 16S rRNA gene sequences from ATCC 33277, 381, and HG66 have identical sequences, indicating the close evolutional distance of these strains. Strain A7436 shared three of its four copies of 16S rRNA sequences identically with those of W83. Together with the single copy from W50, they formed an identical group of sequences. W50 has been known to be a close strain of W83, thus the identical sequences between these two are not surprising. The explanation of identical copies of the 16S rRNA sequence in the genome is apparently due to the gene duplication event and the fact that several strains shared identical duplicated sequences suggested that the duplication event occurred after the speciation. Strains A7436, AJW4, and HG66 had three strain-specific identical sequences with the 4th copy different from the other three. Overall, all the P. gingivalis 16S rRNA gene sequences were extremely similar and often have only a single number of nucleotide mismatches between any two strains (if not identical). Altogether only 16 loci on the gene had nucleotide variations (some could have variations in more than one strain), with the exception of one copy in TDC60, which had a series of A → C or G → C transversions between position 50 and 130 and two single nucleotide insertions at position 174 and 233. It is thus the most divergent sequence of among all 19 P. gingivalis genomes. These aligned and trimmed sequences, including the outgroup sequence from P. asaccharolytica strain DSM 20707, were used to construct a phylogenetic tree based on Kimura's nucleotide substitution model and the result is shown in Figure 1. The phylogenetic tree depicts a likely evolutionary path for these different P. gingivalis strains. The strains 381, ATCC 33277, and HG66 appeared to be closer to the potential common ancestor, based on the tree topology inferred with a close species as the outgroup sequence. The other strains gradually diversified into deeper branching nodes with two of the sequences from strains F0566 and TDC60 as the most deeply branched and mutated from the common ancestor, which was inferred by using the sequence of a neighboring species (PaDSM20707) as an outgroup.

The phylogenetic relationship inferred based on the 16S rRNA gene sequences reported above can only represent the evolution of this particular gene, hence a gene tree. A more comprehensive way of studying the evolutionary relatedness of different genomes is to use as much genomic information as possible in the analysis (i.e., phylogenomics). A popular approach is to use the core proteins for the construction of a tree that may be closer to a true species tree, if such a tree exists, or if there is no true species tree, may reflect more on the relatedness of these strains at the genomic level. The concept of “core” proteins is ideally defined as proteins that are present and required by all the genomes in study, however the identification of such a group of proteins, namely orthologs, is not straightforward and the results vary depending on the criteria used. It is challenging, if not impossible, to identify all the orthologous proteins among a group of genomes. In general, genomes of closer species or strains share more orthologs; however any percent protein sequence identity chosen as the cutoff to test whether a group of homologous proteins are truly orthologs (or paralogs) can always include some false positive and negative orthologs. Nevertheless, one can still hypothesize that a more reliable evolutionary relationship of a group of genomes can be obtained if the use of higher or lower percent identity constrains does not affect the overall tree topology.

To test this hypothesis, the “core” proteins were first identified among all P. gingivalis genomes under different cutoffs. Based on the NCBI annotation, a total of 39,926 protein sequences were identified. However, for some unknown reasons, some of the annotated protein lengths were as short as one or two amino acids. For example, proteins with Genbank IDs GAP82138.1, GAP81676.1, and GAP81848.1 in strain Ando were identified with only 1, 2, and 2 amino acids in length respectively. These clearly were annotation errors caused by the computational bugs in the annotation pipeline. In this analysis, only protein sequences with a minimal length of 50 amino acids were used (a total of 37,667 proteins) for identifying core proteins. They were subject to the “blastclust” program (http://www.ncbi.nlm.nih.gov/Web/Newsltr/Spring04/blastlab.html) to identify clusters of proteins that share a certain degree of sequence homology and with specified alignment length coverage. In this analysis, if a protein is present (i.e., meets the % identity and alignment cutoffs specified) in all 19 genomes (or 20 genomes if PaDSM20707 was used as the outgroup in some results) it is considered as a core/shared protein, and if a protein is only present in a single genome it is considered as a strain-specific unique protein.

Figure 2 shows the potential numbers of both core and unique proteins in the 19 genomes analyzed with “blastclust” by varying two parameters: Sequence percent identity cutoffs (from 95 to 10%) and percent alignment length (90 and 50%). Figure 2A shows that regardless of sequence identity cutoffs; the number of core proteins stays relatively constant around 1,000 with 90% as the alignment length cutoff. The number of core protein groups increased gradually from 1,037 at 95% identity, maximized at 1,045 at 60%, then decreased to 910 at 10%. The reason for the increase from 95 to 60% was due to more core protein groups clustered together at lower % identity. The decrease after 60% identity was due to the fact that different protein groups identified with identity ≥ 60% began to merge into fewer groups. The 1,037 shared proteins were detected under most stringent conditions thus it is reasonable to state that at least 1,037 core proteins were detected based on 19 strains. This number is expectedly smaller than the 1,476 detected in the core genome based on eight P. gingivalis strains (Brunner et al., 2010). It should also be noted that the core genes/proteins are not the same as the “essential” genes, for which only ca. 400 were experimentally detected previously (Klein et al., 2012).

For the purpose of identifying a core/shared set of proteins for constructing a phylogenomic tree, the 1,045 core proteins identified at 60% sequence identity and 90% alignment length cutoffs were used for sequence alignment and tree building. This set of sequences is available for download in the data repository FTP site mentioned in the Material and Methods. In addition, as expected when the percent alignment length was decreased from 90 to 50%, more proteins were identified as core proteins, e.g., from 1,289 at 95% identity cutoff to 1,301 at 60%, due to the fact that more proteins share the same percent identity over shorter sequence length.

Figure 2B shows the number of protein groups that are shared by 2 to 18 genomes (partially shared proteins). The number decreases with the lower % identity because similar protein groups that were identified as separate groups merged into a single (but larger) group due to the more relaxed (lower) % identity (e.g., from 1,927 at 95% to 1,651 at 10%). However, contrary to the core proteins above, which require proteins present in all 19 genomes, when the percent alignment decreased, fewer partially-shared proteins were identified. This is as expected because when the percent alignment cutoff was lowered, a protein group which consists of only members from for example 18 genomes, at higher cutoff, now may find a member in the 19th genome thus disqualifying it as the 18-genome partially-shared group.

Figure 2C shows the number of strain-specific proteins that are present in only one genome with a single copy. Similar to the partially-shared groups, as the % identity decreases the number of unique proteins becomes smaller because more proteins from different genomes were lumped together as a homologous group under a lower % identity, resulting in the loss of the “uniqueness”. In addition, for example, at 60% sequence identity and 90% alignment cutoffs, there were 2,289 proteins identified as present in a single genome, but the number was reduced to 1,044 at 50% alignment cutoff–1,245 proteins lost their uniqueness due to the presence of more “similar” proteins found in other genomes.

For the unique proteins identified, it would be interesting to observe their distribution in the 19 genomes and the result may help understand which genomes possess more or fewer unique proteins. Figure 3 shows the distribution of the 1,044 unique proteins identified with the 50% alignment cutoff (Figure 3A) and 2,289 with 90% cutoff (Figure 3B). Regardless of the sequence identity and percent alignment cutoff, the results show that some strains possess significantly more unique proteins than others. Four strains, F0566, F0568, F0569, and JCVI SC001 have a significantly higher number of unique proteins under all identification conditions - as high as 96–249 unique proteins (for percent identity 10–95% at 50% alignment length) in the case of F0566 (Figure 3A). On the other hand, strains 381, ATCC 33277, A7436, and W83 are the four strains with the lowest number of unique proteins, only 10–15 unique proteins (10–95% identity; 50% alignment) in the case of W83 (Figure 3A). Interestingly, strain W50, the closest strain to W83, encodes more unique proteins (30–46) even though it is an unfinished draft genome. Apparently the incompletes of the draft genomes are not the cause for the difference in the number of unique proteins (JCVI SC001 and all the F strains are draft genomes). This further suggests that the gaps in the draft genomes may contain only repeated sequences that either do not encode for proteins, or encode repeated proteins that do not contribute much to the genome's uniqueness. However, we cannot rule out the possibility that these regions could still have some unique proteins which would drive up the number already observed.

Another noteworthy observation is that there is a consistent gap between the data points 95% and 90% identity when searching for unique proteins in all strains under both alignment conditions (Figure 3). This suggests that the more proteins identified as unique at 95% became “similar” at 90%. Hence 90% sequence identity may be an ideal cutoff for differentiating homologs and unique proteins, at least at the strain level.

Table 6 lists the percentage of proteins that were annotated in NCBI as the hypothetical proteins (functionally unknown) and the percentage of the unique proteins that were identified with 80% as the sequence identity and 50% alignment length. The total percent hypothetical proteins range from 26% (W50) to as high as 46% (F0566 and F0568), whereas the majority of the unique proteins are hypothetical, from 68% (TDC60) to 100% (W83). Thus, until more functions of the hypothetical proteins are understood, it will still be challenging to understand what each genome's overall “specialty” functions conferred by the unique genes. To give a glimpse of what each genome's most unique functions are, based on the currently available information, Table 7 lists the functional annotations of the non-hypothetical unique proteins for each of the 19 genomes (with default BLASTP parameter, i.e., expected e value ≤ 10) (Altschul et al., 1997). All of them are among the proteins identified above under the most stringent parameters in terms of uniqueness–50% sequence identity and 50% alignment length. Strain JCVI SC001, an environment isolate from a hospital sink drain, has the most diverse functions encoded by these unique proteins. The unique toxin-antitoxin system detected in strain F0569 (a clinical isolate from subgingival plaque biofilm) is also of interest. The toxin-antitoxin system genes when carried on a plasmid is often referred to as the post-segregational killing (PSK) system (Gerdes, 2000) while when carried in the chromosome such a system is involved in stress response and “programmed cell death” (Hayes, 2003). Although indigenous plasmids have never been detected in P. gingivalis this does not rule out finding one in the future. Thus, which type of function that this toxin-antitoxin system is involved in can only be speculated at this point of time. Whether these annotations translate to unique functions of the genome, require further investigation to ensure there are no other non-homologous proteins that play similar functions. All of the unique proteins identified are available by strain in the FTP data repository.

Once a group of putative core proteins is identified, they can be concatenated and aligned together and used for compiling a phylogenomic tree to infer a possible evolutionary relationship at a level closer to the species than just any single gene. In this analysis, the 1,045 proteins shared by all 19 genomes at 60% sequence identity and 90% alignment cutoffs (Figure 1A), were first aligned individually with the “mafft” software (Katoh and Standley, 2013). Each of the 1,045 protein sets contained exactly 19 aligned sequences, one from each of the 19 genomes. The aligned proteins were concatenated in the same protein order. This generated a set of 19 mega protein sequences with each consisting of 1,045 concatenated aligned sequences. The poorly aligned sequence regions, including leading and trailing unaligned portions of the sequences, as well as low-confidence parts of the alignment, such as positions that contain many gaps, were removed with the “Gblock” tool (V 0.91) (Talavera and Castresana, 2007). After the Gblock screening, a final set of 19 aligned protein sequences, each with a length of 395,174 amino acids were used for constructing an unrooted tree. However, among the 395,174 aligned amino acids, only 17,389 positions had at least two different amino acids across proteins of all 19 P. gingivalis genomes, the remaining 377,785 were all the same amino acids across all genomes. Thus, only those 17,389 informative or effective positions contributed to the pairwise distances calculated among all genomes. Figure 4A is the result of the unrooted tree compiled based on the 1,045 shared proteins processed as described above. The overall topology is quite different from that of the 16S rRNA tree (Figure 1) with the exception of two very closely related groups of strains, one consists of strains 381, ATCC 33277, and HG66 and another A7436, W50, and W83. This is not surprising because both groups have members with identical 16S rRNA sequences hence their shared protein sequences are closer to each other in the group than other genomes.

To test whether including proteins from the outgroup species will result in a tree more similar to that of 16S rRNA, i.e., a tree that is rooted at a potential common ancestor for these strains, ortholog candidates were first identified from the genome of P. asaccharolytica DSM 20707, of which the 16S rRNA sequence was also used for the 16S rRNA tree. At 90% alignment length cutoff, the number of homologous proteins in P. asaccharolytica decreases as the percent sequence identity cutoff increases. The numbers of protein homologous to any of the 1,045 core proteins used for the unrooted tree above are 436, 271, 146, 36, 7, 1, and 0 respectively for percent identity cutoffs 50, 60, 70, 80, 85, 90, and 95%. Figures 4B,C are the two rooted phylogenetic trees constructed based on the 436 (50% identity) and 36 (80% identity) proteins shared between P. asaccharolytica DSM 20707 and all 19 P. gingivalis strains. After Gblocks screening, the length of the aligned sequences were 12,646 (80% identity) and 177,272 (50% identity) amino acids respectively and the number of effective amino acids positions are 154 and 4,771 respectively. In general, the branch lengths increased with more effective amino acids positions which resulted in greater distances. Again, the only consistent close clusters were the two grouped with identical 16S rRNA, i.e., the group of 381, ATCC 32277, and HG66, and of A7436, W50, and W83.

Figure 4D is the rooted tree constructed using the software PhyloPhlAn (Segata et al., 2013) version 0.99 (8 May 2013). All 41,625 proteins annotated for the 20 genomes were subject to PhyloPhlAn with the default parameters that excluded proteins shorter than 30 amino acids in length. PhyloPhlAn finds among the input protein matches to a pre-set of the 400 most conserved proteins for extracting the phylogenetic signals. A total of 264 query proteins were matched to the 400 preset core but only 225 were present in all 20 genomes. These proteins were then aligned individually and subsampled based on a sophisticated procedure provided by PhyloPhlAn, which emphasizes regions both universally conserved and phylogenetically discriminating. The final aligned, subsampled, and concatenated sequences had a length of 3,082 aligned amino acids with 840 effective positions. The PhyloPhlAn tree is shown in Figure 4D. Similar to the two rooted trees (Figures 4B,C) and the 16S rRNA tree (Figure 1) the PhyloPhlAn tree also placed the three strains ATCC 33277, 381, and HG66 closest (but much closer) to the root and the remaining strains in a more linearly nested topology.

In summary, the only consensus based on interpretation of the three rooted protein trees and the 16S rRNA tree is that the group ATCC 33277, 381, and HG66 is less evolved and closest to the common ancestor of this species (inferred based on the distance to the root). Strains W83, W50, and A7436 consistently formed a close group regardless of how the trees were built, but their exact phylogenetic position is inconclusive based on these analyses. Strains F0185, F0568, F0570, W4087, and MP4-504 are also found to be in the proximity of each other, although not as close as the two groups mentioned above. In general more effective/informative aligned amino acid positions resulted in longer branches and pairwise distances. To this end, the unrooted tree (Figure 4A) has the best resolution to reveal the similarity/differences among these strains, in the most genome-wide manner. Until a group of true orthologous proteins are identified (together with the outgroup) a true phylogenetic tree that infers the evolutionary path for this species will not be accessible.

Moving up the scale for comparison, one possible way is the whole genome nucleotide alignment with a commonly used software NUCleotide MUMmer (NUCmer), which identified nucleotide MUMs–minimal unique matches between two genomic sequences (Delcher et al., 2002). Figure 5 shows some of the pairwise alignment results of the 19 P. gingivalis genomes. Figure 5A is the nucleotide alignment between strains 381 and ATCC 33277 and the almost perfect diagonal high similarity (red) match line indicates highly similar sequences, with only two visible exceptions – one inversion and one insertion (to 381)/deletion (to ATCC 33277). Interestingly the inverted sequence almost matches the inserted sequence; apparently the inverted sequence was duplicated in the 381 genome and inserted somewhere else in the genome, where the ATCC 33277 genome shows no counterpart. The high DNA sequence similarity between 381 and ATCC 33277 is also supported by the identical 16S rRNA gene sequence and copy numbers (Figure 1) as well as the protein-based phylogenetic relationships (Figure 4), even though their genomes are not far from identical. The phenomenon that a fairly large chunk of genomic sequence was duplicated and inserted elsewhere in the genome is only observed in strain 381, as evidenced by the NUCmer self-alignment of its genome (data available from the FTP site), but very similar to the alignment between 381 and ATCC 33277). No duplication event was observed in the self-alignment of the other 18 genomes.

Strain HG66 is the genome that is closest to 381 and ATCC 33277 based on 16S rRNA genes and protein sequences, on the other hand it shows the disconnected high similarity match lines, which indicates more large-scale genomic arrangement between the two close strains – between 381 and HG66 (Figure 5B) and between ATCC 33277 and HG66 (Figure 5C).

The second closest groups of strains are A7436, W50, and W83 and their nucleotide sequences are also highly similar based on the NUMMER plots (Figures 5D,E). However, the contigs of the unfinished draft genome of W50 were rearranged by NUCmer in the order based on the similarity to the W83 sequence. Whether there is a large scale genomic rearrangement between W83 and W50 cannot be known until the genome of W50 is completed. Strain A7436, a finished genome, shows only one inversion of the genome when compared to that of W83 (Figure 5E). The fact that A7436 is not as close to W50 and W83 as the distance between HG66 and 381 (or ATCC 33277) based on 16S rRNA and protein phylogeny (Figures 1, 4), suggests that the genomes of the group of HG66, 381, and ATCC 33277 have higher genomic sequence rearrangement activity than the A7436-W50-W83 group. The next genome which is closest to the A7436-W50-W83 group is strain AJW4, with several visible (larger fragments) of insertions/deletions and inversions when compared to A7436 (arrows heads in Figure 5F). This relationship is also consistent with the 16S rRNA gene tree (Figure 1).

Another interesting observation is the alignment between JCVI SC001 and TDC60. These two strains are not among the closest groups based on the 16S rRNA and protein sequences (Figures 1, 4). The NUCmer plot between these two genomes appears to be a straight diagonal red line (Figure 5G), similar to that between 381 and ATCC 33277. However, since the genomic sequence of JCVI SC001 was not really completed and closed to a circular chromosomal format, the 284 de novo assembled contigs were mapped to the genome of TDC60 and the gaps were filled with Ns to form a single pseudo-contig (Genbank Accession CM001843) (McLean et al., 2013). Thus, the contig order in the published single contig genomic sequence of JCVI SC001 may not be correct and the sequence similarity between JCVI SC001 and TDC60 may not be as “straight” as indicated in the NUCmer plot. In fact when the plot was filtered to show only the region with percent identity ≥ 99%, the red line became fragmented with large gaps (Figure 5H), indicating that a large portion of the genomic sequences between these two strains are under 99% similarity.

The complete pair-wise NUCmer plots of the 19 P. gingivalis genomes can be viewed in a specifically designed interactive webpage at http://bioinformatics.forsyth.org/publication/20160425. The web page provides interactive tools to choose any two P. gingivalis genomes for the NUCmer results, as well as the possibility of viewing the alignment at various percent sequence identity cutoffs.

The NUCmer plots above are limited to viewing comparisons only between two given genomes. To view and compare nucleotide difference/similarity for all genomes on the same plot, the overall oligonucleotide composition and frequency can be measured along the entire genome and the results can be plotted out and visually compared to each other. This analysis started by collecting all the possible 20-mer sequences in all 20 genomes and then count for each 20-mer how many genomes have each particular sequence. The number of genomes (genomic frequency) for each 20-mer thus ranges from 1 (unique) to 20 (universal). The frequencies can be calculated and plotted along the entire genome by taking every 20-mer from the beginning to the end of the genome. Figure 6A depicts the results of the 20-mer oligonucleotide frequencies among all 20 genomes (including the out-group P. asaccharolytica DSM 20707). If a region of a genome is shared by all other 20 genomes, it is colored black; and if a region is unique to the genome itself, it is colored bright yellow. In other words, a black region means that all the possible 20-mer sequences appeared in all tested genomes, whereas the brightest yellow regions have unique 20-mer sequences that are only found in one genome. For easy comparison, the order or the genomes shown in Figures 6A,B were arranged according to that of the 16S rRNA tree (Figure 1) with a dendrogram reflecting similar tree topology. As expected, the two closest strains ATCC 33277 and 381 share almost identical 20-mer frequency patterns, with the exception of a small insertion at nucleotide position ca. 1,400,000, which is also detected by the NUCmer plot in Figure 5A. The genome of strain 381 is ca. 24 Kbps longer than that of ATCC 33277 due to this insert and the length difference is illustrated in Figure 6A because the length of the bars were based on actual genome sizes. Another interesting example observation is that even though strain JCVI SC001 is closest to SDJ2 due to their identical 16S rRNA sequences (Figure 1), their oligonucleotide frequency patterns are quite different, with each showing unique regions (brighter colors) at different places. This can be due to two possibilities: (1) the artificial order of the unfinished sequence contigs (in the plots, contig order was the same as that in the downloaded sequences); and (2) the bona fide differences in sequence. This is also true for the other three genomes A7436, W83, and W50, which share identical 16S rRNA sequences but exhibit distinct frequency patterns.

When the frequencies were plotted only in the open reading frames (ORFs), most of the ORFs appeared dark in color and are shared by the majority of the P. gingivalis genomes (Figure 6B). Those ORFs with colors close to yellow should account for the differences of the number of unique proteins previously identified (Figure 3). These differences in nucleotide sequences are essentially reflected in the genes and then translated to proteins, and ultimately accounting for differences in biological functions.

Finally the out-group P. asaccharolytica DSM 20707 shows mostly yellow colors in the plot, which is as expected and means that P. asaccharolytica does not share much of the 20-mer oligonucleotide sequences with P. gingivalis. Interestingly, by lowering the oligomer size to 14 bases (14-mer), most of the DSM 20707 genome appears black, meaning that these two different species share most of the 14-mer sequences (data not shown). When the plot was generated with 15-mer sequences, the P. asaccharolytica DSM 20707 genome started to show patches of yellow area (data not shown), indicating that some unique 15-mer sequences are present between these two species. With 15-mer all P. gingivalis genomes are black in the plot (data not shown), meaning 15-mer is too short and have not enough resolution power to differentiate unique regions among P. gingivalis genomes. Hence whether the oligomer frequency analysis can detect unique/shared regions in a group of genomes, depends on the size of the oligomer. The choice of 20-mer was able to identify unique regions among the strains of P. gingivalis, as shown in Figures 6A,B, yet is too sensitive for a different species.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer DY and handling Editor declared their shared affiliation and the handling Editor states that the process nevertheless met the standards of a fair and objective review.