Analysis of the NCBI genome database

On 3 April 2017, we accessed all data through 2016 in the NCBI genome database (https://www.ncbi.nlm.nih.gov/genome/browse). For analysis purposes, we considered the "FCB Bacteroidetes/Chlorobi," "Terrabacteria Firmicutes," "Terrabacteria Actinobacteria," and "All Proteobacteria" subgroups to represent the phyla Bacteroidetes, Firmicutes, Actinobacteria, and Proteobacteria, respectively. The levels "Complete" and "Chromosome" were considered finished or complete genomes, and "Scaffolds" and "Contigs" were considered draft or incomplete genomes. Analysis was performed using R (v3.3.2) [48] and in-house parsing scripts in Python (v2.7.10). We excluded all entries marked “Candidatus” and entries that did not include full species identification specifying both the genus and the species.

Analysis of Bacteroidetes genomes

Bacteroidetes genomes were categorized by class and sequencing status. The sequences of all copies of genes encoding 16S rRNA were extracted from the complete genomes. Alignment was performed using the MAFFT (v7.222) plug-in [49] in Geneious (v10.2.3) [50] set to the following parameters: automatic detection for the algorithm; a scoring matrix of 200 PAM/k = 2; a gap opening penalty of 1.53; and an offset value of 0.123. A consensus sequence was then created for each species using the default Geneious settings. Finally, a phylogenetic tree was constructed using the PhyML plug-in [51]. The tree was based on an HKY85 substitution model [52] and features 100-bootstrap branch support; optimized topology, branch lengths, rates, and nearest-neighbor interchange (NNI); and a subtree pruning and regrafting (SPR) topology search. The tree was simplified to the class level according to the NCBI taxonomy. The isolation sources of all complete genomes were identified via the NCBI BioProject and BioSample databases. The sources were classified into three categories: environmental (soil, fresh or marine water, sludge, mud, plants and algae samples), animal (insects, mollusks, fish, birds, cattle, and domestic animals), and human (isolated from various body sites of healthy or sick individuals).

Complete Bacteroidia genomes and genomic repeats

For complete genomes, we used Repeatoire to identify genomic repetitions that were at least 95% similar to each other and longer than 500 base pairs (bp) [54], visualizing them with Circos (v0.69) [55]. In each genome, the repeat’s initial location was fixed as the starting point for all of the links to the other positions of that repeat, and we color-coded the copy numbers of each repetition.

Simulated read assembly

For each genome, artificial reads were produced using ART software (v2.5.8) [56] set for paired-end reads of 250 nucleotides (nt) each, 500 nt insert size, and a coverage of 40; the built-in MiSeq simulation profile (v1) was used. For the P. gingivalis genomes, eleven de novo assemblers were tested using the default parameters: A5-miseq (v20160825) [57]; CodonCode Aligner (CCA v7.0.1); CLC Genomics Workbench (CLC-GW v8.5.1); fermi (v1.1) [58]; Geneious (v10.2.3) [50]; Minia (v2.0.3) [59]; MIRA (v4.0.2) [60]; PERGA (v0.5.03.02) [61]; SOAPdenovo (v2.04) [62]; SPAdes (v3.10.1) [63]; and Velvet (v1.2.10) [64]. For CCA and Geneious, reads were preassembled with PEAR (v0.9.10) [65]. Where required, KmerGenie (v1.7016) was used to determine the k-mer length parameter [66].

After the initial test step, we selected three software packages for assembly of all 24 genomes. The first was A5-miseq, in which we used the default parameters. For Geneious, reads were preassembled with PEAR and then assembled using the default Medium Sensitivity/Fast setting and generating consensus with a 50% strict threshold for calling bases. Finally, SPAdes was used with the default parameters but with the “careful” option and k-mer lengths of 21 to 127.

For each genome and assembly tool, we rejected contigs of less than 1 Kbp because they are not informative. We assessed the assemblies with QUAST (v4.5) [67]. The fidelity of each assembly was evaluated by mapping each set of contigs to the reference genome using Geneious mapper (Medium Sensitivity/Fast option); the unmapped contigs were rejected, and QUAST was then used again. The correlation between the number of repeats (copy number > 3) and each assembly was tested using the “ggplot2” and “nortest” packages in R. We identified the P. gingivalis reference genome segments that correspond to assembly gaps and classified these into five categories: genomic islands, rrn operons, coding sequences (CDS) with repeated domains, intergenic sequences, and insertion sequences or miniature inverted-repeat transposable elements (IS/MITE).

Bacterial strain cultures and DNA extraction

We purchased the ATCC 33277 and W83 (also known as BAA-308) P. gingivalis strains from the American Type Culture Collection-LGC Standards (Manassas, VA, USA) in September 2006; a low (< 20) passage number was used. P. gingivalis TDC60 (also known as JCM 19600) was purchased from the Japan Collection of Microorganisms (Riken BioResource Center, Koyadai, Japan) in November 2015, and a low (< 10) passage number was used. All strains were cultured on Columbia European Pharmacopoeia agar plates (Conda, Madrid, Spain) supplemented with 5% (v/v) defibrinated horse blood (Eurobio, Courtaboeuf, France), 5 g/L yeast extract (Conda), 25 mg/L hemin (Sigma-Aldrich, Saint-Quentin Fallavier, France), and 10 mg/L menadione (Sigma-Aldrich). The cultures were incubated in an anaerobic chamber in a Whitley DG500 Workstation (Don Whitley Scientific, Shipley, UK) for 5 days at 37 °C in an atmosphere composed of 80% N₂, 10% H₂, and 10% CO₂.

For DNA extraction, each strain was cultured for 48 h at 37 °C under the same atmospheric conditions described above. This was done in 50 mL BHI broth (BioMérieux, Marcy l'Etoile, France) enriched with 5 g/L yeast extract (Conda), 25 mg/L hemin (Sigma-Aldrich), and 10 mg menadione (Sigma-Aldrich). After harvesting, the cells were washed twice in Dulbecco's PBS (Dominique Dutscher, Brumath, France). A QIAamp DNA Mini Kit (QIAGEN, Courtaboeuf, France) was used for cell lysis and protein denaturation. The following steps were performed using standard methods: DNA precipitation with 5 mol/L NaCl (Sigma-Aldrich) and 0.7 volumes of cold isopropanol (VWR Chemicals, Fontenay-sous-Bois, France), two washes in 70° ethanol (WWR), and resuspension in sterile Milli-Q ultrapure water (Merck, Darmstadt, Germany).

DNA library preparation and PacBio SMRT sequencing

Barcoded DNA library preparation and single molecule real-time (SMRT) sequencing were performed using the Genome et Transcriptome (GeT) GénoToul platform (Toulouse, France) according to the manufacturers’ instructions. Quality control was performed at each step. DNA was purified using AMPure PB beads (Pacific Biosciences, Menlo Park, CA, USA), and the mass of dsDNA was verified using a Qubit fluorometer (Thermo Fisher Scientific, Villebon sur Yvette, France). The purity of the DNA was determined based on absorbance ratio using a Nanodrop spectrophotometer (Thermo Fisher Scientific). Sizing measurements were performed using a Fragment Analyzer (Advanced Analytical Technologies, Evry, France). In brief, each sample was diluted to 10 μg/mL and sheared on a Megaruptor (Diagenode, Seraing, Belgium). Using 5 μM of various barcoded adapters, a SMRTbell Barcoded Adapter Prep Kit (Pac Bio) was used to repair and ligate 150 ng of DNA fragments. After end-repair and ligation using a SMRTbell DNA Damage Repair Kit, the samples were pooled. To remove unligated DNA fragments, the library was treated with an exonuclease cocktail consisting of 1.81 U/μL Exo III and 0.18 U/μL Exo VII (PacBio).

Library selection in the 6-50 Kbp range was performed using BluePippin 0.75% agarose cassettes (Sage Science, Beverly, MA, USA). Primers were annealed to the size-selected SMRTbell with the full-length libraries. The primer-template complex was then bound to the P5 enzyme using a 10:1 ratio of polymerase:SMRTbell for 4 h at 30 °C. The magnetic bead loading step was conducted at 4 °C for 1 h. The complexes were placed into the PacBio RS II sequencer, which was configured to run continuously for 6 h at a sequencing concentration of 80 pM. The sequencing results were validated using the NG6 integrated next-generation sequencing storage and processing environment [68].

Genome assembly and finishing strategy

We mapped the contigs obtained from the in silico reads to each reference genome. We identified and retrieved the sequences not covered by the contigs, which corresponded to the assembly gaps. Because these gaps are repeated regions, we clustered the sequences at 99% nucleotide identity and generated a consensus sequence. Thus, each consensus represents one type of repeated region. We used canu (v1.3) to correct, trim, and assemble the raw reads [69]. We mapped the corrected long reads to each consensus sequence that was previously identified; for each, we selected only reads overhanging at least 500 nt at both ends. The selected reads were de novo assembled using Geneious at 100% nucleotide identity. They were then used to scaffold the contigs generated by canu and to finish the assembly and reconstruct the genome organization.

Genome annotation and biocuration

Genomes were annotated using Prokka (v1.12-beta) [71], Genix online [72], and RASTtk (v1.3.0) [73]. The annotations for the gene starts/ends, gene names, gene product descriptions, gene status (gene, pseudogene by stop in-frame or frame-shift), EC numbers, and functional descriptions were all manually biocurated [74]. For this purpose, we performed NCBI BLAST searches [75] against non-redundant databases as well as domain searches using the Conserved Domains search tool [76]. The results were then compared to the precomputed annotations from MicroScope [77].

The presence of signal peptides was analyzed using the SignalP 4.1 server [78] and SOSUIsignal [79]. Protein subcellular localization predictions were made using PSORTb (v3.0) [80] and CELLO (v2.5) [81]. The outer membrane proteins predicted by these tools were then confirmed using BOMP [82], and LipoP (v1.0) [83] and DOLOP [84] were used to predict the lipoproteins.

Insertion sequence transposases were renamed as per the ISfinder [85] nomenclature. Clustered Regularly Interspaced Short Palindromic Repeats (CRISPRs) were identified using CRISPRfinder [86], CRISPRDetect [87], CRISPI [88] webtools and the CRT plug-in [89] with Geneious. We predicted genomic islands using IslandViewer (v4) [90], which has precomputed predictions from SIGI-HMM [91], IslandPath [92], and IslandPick [93]. Island prediction was also performed using EGID [94], which implements AlienHunter [95], SIGI-HMM, IslandPath, INDeGenIUS [96], and PAI-IDA [97].

All genomic sequences and the reads used to produce them were deposited in GenBank and in the NCBI BioProject database with links to the BioProject accession number PRJNA393092.

Comparative genomics

Genome-level comparisons were made with progressiveMauve [98] using the default settings. The same program was used to identify locally collinear blocks (LCBs) and single-nucleotide polymorphisms (SNPs).

For each strain, we compared the available annotation in the NCBI database to our manually biocurated annotations feature-by-feature: rRNA; tRNA; tmRNA; ncRNA; regulatory; repeat region; and coding DNA sequence (CDS). The CDS were separated into pseudogenes and “true” coding regions. For the latter category, the CDS/pseudogenes present in both annotation versions were termed “common,” and those that were present only in our new version were called “new.” We noted all instances of changes from a CDS to a pseudogene and vice versa, fusion (more than one interval becoming one interval) and separation (one interval becoming more than one interval), and changes in the coding strand. Finally, we also analyzed the changes in start/stop codons.

We concluded our study with a pangenomic analysis of the P. gingivalis strains ATCC 33277, TDC60 and W83. The CDS and pseudogenes were classified into six groups: multiple-copy genes (nucleoid-associated proteins, tra genes, xer genes, and transposases in ISPg); those that appeared in only one copy in all of the strains and were highly conserved (> 97% nucleotide identity); those present in only one copy in all of the strains but with sequence divergence (< 97% nucleotide identity); those present in all three strains but present in more than one copy in at least one strain and/or having pseudogenes; those present only in two strains and either copied or with pseudogenes; and those present in only one strain.

The Bacteroidetes phylum is represented by few genomes, mostly incomplete and with variable numbers of contigs

Bacterial genomes represent approximately 85% of the available genomes in the NCBI genome database. Despite the fact that Proteobacteria represents a relatively minor constituent of human microbiota, it is the most sequenced bacterial phylum, even in comparison to the two most abundant phyla, Bacteroidetes and Firmicutes (Additional file 2: Figure S1a). Furthermore, the accumulation of incomplete draft bacterial genomes is notable, comprising 85-95% of all entries depending on the phylum (Additional file 2: Figure S1b).

Within the phylum Bacteroidetes, the classes Bacteroidia and Flavobacteriia constitute 90% of the listed genomes. Flavobacteriia is mostly associated with environmentally obtained samples, particularly with aquatic species such as water and fish pathogens, and rarely (18% of complete genomes) with humans. In contrast, Bacteroidia isolation sources are human in more than 75% of cases (Additional file 3: Figure S2a). Although generally considered part of the normal microbiota, Bacteroidia are pathobionts that can become pathogens upon dysbiosis. Prior to the introduction of massively parallel sequencing, only four complete Bacteroidia genomes were published. Since 2007, this number has grown exponentially, mainly due to the large number of draft genomes that have accumulated to such an extent that today 10 times as many draft genomes as complete genomes are available (Additional file 3: Figure S2b).

In addition to the preponderance of draft genomes, the Bacteroidia draft genomes are extremely fragmented; half of them contain more than 75 contigs. The number of contigs per draft ranges widely from 2 in Prevotella oryzae DSM17970 to 4357 in Bacteroides acidifaciens 1a3B. Observing the number of contigs per draft in bacterial genera with at least five draft genomes, we noted that 50% of all entries are Bacteroides and that Bacteroides, Parabacteroides, and Prevotella all have at least one draft genome featuring more than 500 contigs. Based on their interquartile ranges (IQRs), the genera with less variation were found to be Alistipes and Tannerella; coincidentally, these were the genera with the smallest number of draft genomes. At 93 contigs per draft, Bacteroides and Tannerella have the highest median number of draft genomes (Additional file 4: Figure S3a). Even when only species having at least two complete genomes (B. dorei, B. fragilis, B. ovatus, B. thetaiotaomicron, P. gingivalis, and T. forsythia) are considered, this variability is still quite large. B. fragilis represents almost 20% of all Bacteroidia drafts and is the most variable in terms of the number of contigs per draft. B. dorei and T. forsythia, although the species with the smallest number of drafts, are the least variable. Finally, P. gingivalis has an intermediate profile (Additional file 4: Figure S3b).

Draft variability cannot be explained by sequencing and assembly methods

Several hypotheses could explain the observed variability in the features of the draft genomes. It could be linked to differences in the assembly strategies (algorithm types or specific assembly tools) and/or the massively parallel technologies used. To test these hypotheses, we collected all of the available information on sequencing and assembly methodology for all of the draft genomes of the six species studied here (Additional file 5: Table S2). Of the 166 drafts, 144 were sequenced using Illumina technology, and only six assemblers were used for de novo assembly. Five assemblers used the De Bruijn graphs (Allpaths, Velvet, ABySS, CLC Genomics Workbench), covering 40% of the drafts. The final MaSuRCA assembler, which is based upon the de Bruijn and Overlap-Layout-Consensus (OLC) approaches, was used for the remaining 60% of the drafts. However, the overrepresentation of MaSuRCA is due to a single sequencing project of different strains of B. fragilis at the Institute for Genome Sciences (University of Maryland, Baltimore, MD, USA). Once again, the results (whether per assembler or per bacterial species) were very diverse. In B. fragilis projects using Illumina for which an assembler was specified (99 of 107 drafts), MaSuRCA generated 31 to 2566 contigs (n = 69), Allpaths produced 5-14 contigs (n = 10), SPAdes generated 73-343 (n = 8), ABySS produced 150-1290 (n = 6), CLC Genomics Workbench produced 5-156 (n = 5), and Velvet, which was only used once, generated 52 contigs. The 20 P. gingivalis draft genomes were sequenced using Illumina technology and then assembled 17 times with Velvet (generating 22-192 contigs) and once each with SPAdes (92 contigs), Celera (104 contigs), and SOAPdenovo (117 contigs). For all other species, the number of drafts with known assemblers is too small (< 10) to draw conclusions. At this point, it therefore seems that technological differences in sequencing and assembly cannot explain the extensive differences in the number of contigs per draft.

The hypothesis that intra-species diversity might be responsible for this variability remains to be explored. Because this hypothesis involves comparative genomics, only complete genomes could be used to explore its validity.

Complete Porphyromonas gingivalis genomes are the most diverse and repeated genomes within Bacteroidia

By calculating the average nucleotide identity (ANI), we separated two groups of species. The first group contained the Bacteroides genus, whereas the second group included P. gingivalis and T. forsythia (Fig. 1a). ANI values offer a robust and sensitive way to measure the evolutionary relationship of bacterial strains [99]. With the exception of the P. gingivalis group (11 of 21 pairs), the ANI values of all the species are high and uniform (Fig. 1b). These values suggest that only minor genomic differences exist within species and refute the hypothesis that global intra-species diversity could explain the variability in the number of contigs per draft genome. However, because small variations such as repeated genomic regions cannot be detected by the ANI method and since repeats have been declared to be the main cause of assembly gaps, we decided to evaluate the possible role of repeats in the fragmentation of draft genomes.

We determined the number and locations of each repeat type. Because Illumina is the industry standard and its MiSeq platform generates 2 x 250 nt paired reads, we set the low threshold to 500 bp and looked for the repeats. The histogram in Additional file 6: Figure S4 shows the mean number of repeats in each species (2 to 10 copies at 95% identity). Notably, although B. fragilis displays the greatest variation in number of contigs per draft, it has an intermediate number of repeats, and these have low copy numbers and occur a maximum of only six times. P. gingivalis can once again be distinguished as having the lowest total number of repeats but the highest number of copies for each repeat type, with a total of approximately 40 different repeats copied more than 10 times (Additional file 6: Figure S4).

We used a Circos figure to visualize the genomic repeats with more than 3 copies in each complete genome. The plot illustrates the distribution of the repeated loci and their copy numbers in the form of a heat map, the color of which changes from blue to red as the count increases (Fig. 2). The plots show that all six available strains of B. fragilis possess few types of repeats and that repeats are not frequent. In contrast, all of the other Bacteroides species genomes (B. dorei, B. ovatus, and B. thetaiotaomicron) have more repeats that occur more often and display greater variability within the strains. Finally, we observed that the seven P. gingivalis strains had the highest genomic complexity and diversity with respect to repeat frequency.

The intra-species diversity of P. gingivalis and its large number of highly frequent repeats make this species an interesting model for analysis of the impact of genomic repeats on bacterial genome assembly, especially in Bacteroidetes. The P. gingivalis strains appear as three branches on the DNA-DNA distance tree. Two of these branches correspond to the previously mentioned closely related strains: ATCC 33277 with 381, and W83 with A7436 (Fig. 3). The ATCC 33277/381 branch contains the genomes with the most repeats, with some loci repeated more than 25 times, followed by the W83/A7634 branch. The remaining branch contains TDC60, AJW4, and A7A1-28, which display the lowest numbers and frequencies of repetition. With the exception of W83 (from Bonn, Germany) and TDC60 (from Tokyo, Japan), all of the cited strains were originally isolated in the USA. Despite this common origin, there is no apparent link between isolation populations and genomic repetition frequencies. To test whether the variation in genomic repeat counts affects assembly completion and since the real sequencing reads from the initial sequencing projects were not available, we produced in silico reads based on the complete published genome of the seven P. gingivalis strains and in silico simulated paired-end reads.

Simulated sequencing reveals a correlation between contig and repeat counts

To generate the artificial reads, we used ART software. This set of tools mimics the real sequencing process and permits simulation of the data that are produced by massively parallel methods; the simulated data display each technology’s inherent empirical errors and profile qualities. For the seven complete P. gingivalis genomes, 11 assemblers were tested. The assembler spectrum was restricted to software that can treat Illumina reads but was otherwise chosen to be as wide as possible. Three of the assemblers are part of commercial software suites (Geneious, CLC Genomics Workbench, and CodonCode Aligner); the others are freely available for academic purposes. Since de Bruijn graphs are inescapable (they are used by A5-miseq, CLC Genomics Workbench, Geneious, Minia, MIRA, SOAPdenovo2, SPAdes, and Velvet), we made an effort to also include different assembly methods, including Overlap-Layout-Consensus (CodonCode Aligner), string graph (fermi), and greedy algorithm (PERGA).

We assessed each assembly’s fragmentation and compared it to that of the others using QUAST’s N50 parameter. For all P. gingivalis strains, Geneious produces the longest contigs and the fewest per assembly. SPAdes and A5miseq are just behind, with very similar performances. The other assemblers generate more fragmented assemblies, and some (Velvet and PERGA) do not even seem suitable for this bacterial group (Fig. 4a). Closer examination of the results of the three assemblers producing the highest N50 (Fig. 4b, upper panel) shows that the number of contigs obtained is below the median (90.5) of the published draft genomes assembled with real reads (Additional file 4: Figure S3b).

N50 can be misleading, however, because it provides no information on assembly accuracy; it can therefore give the impression that one tool is more algorithmically efficient than another even when its results do not agree with the biological sequence. Since the simulated reads were created from the complete NCBI reference genomes, mapping the resulting contigs against the genomes should be biologically consistent. The unmapped contigs would therefore be rejected as believed to be incorrect, resulting in the graph shown in the lower panel of Fig. 4b. The plot demonstrates that despite their high N50 values, Geneious contigs, especially the longer ones, are not consistent with the real sequence. The plot line will therefore fall short of the expected genome size, resulting in drafts that are not at all complete. In contrast, A5miseq and SPAdes yield similar plots since they accurately reproduce the biological sequences, despite the fact that their resulting contigs are shorter. Evaluating assembly completeness solely based on the N50 metric is a bad idea, and it is preferable to obtain a draft genome that is correct even if it is more fragmented. Consequently, A5miseq and SPAdes produce the best results for the P. gingivalis genomes (modeled with MiSeq reads), and we confirmed this to be true for the other Bacteroidia species as well (Additional file 7: Figure S5a).

A direct positive linear correlation can be observed in the number of P. gingivalis genomic repeats (at least three copied loci) and the number of biologically sound contigs produced by A5miseq or SPAdes from the in silico simulated reads (Fig. 4c). This correlation is also observed for the other complete Bacteroidia genomes studied here (Additional file 7: Figure S5b). In the case of P. gingivalis, it is true even when all genomic repeats, including duplications, are included (Additional file 7: Figure S5c).

We can classify the functions of the elements that are annotated in the breaking points of the A5miseq and SPAdes P. gingivalis assemblies. These gaps coincide with the repeated regions in the genome, and two thirds of them correspond to copied insertion sequences (IS) or miniature inverted-repeat transposable elements (MITEs) with a maximum length of 1.1 Kb. In order of their occurrence, the other four categories are: intergenic unannotated regions; repeats found in coding regions (such as gingipains) that have several paralogs and internal repeated motifs; rrn ribosomal RNA operons (present in four copies in P. gingivalis); and, the least frequent, genomic islands (Fig. 4d).

This study shows that to correctly assemble a genome from Illumina MiSeq reads of a few hundred nucleotides, it is helpful and even indispensable to have at least a good estimate of the quantity and frequency of genomic repeats. This knowledge permits the calculation of contig counts with little or no misassembly expected. It could also help in the choice of the best finishing strategy for the assembly. However, as seen here, this knowledge is often impossible to obtain in advance, since even among bacterial species the strains are very diverse.

Choosing Porphyromonas gingivalis strains for resequencing and reassembly

Taking into account our initial results and observations, we chose to resequence three P. gingivalis strains (ATCC 33277, TDC60, and W83) using long-read technology from PacBio. We chose these strains because they are commercially available, allowing other researchers to reproduce our results, and because of their geographic diversity (they appear in North America, Europe, and Asia). Moreover, they were the first P. gingivalis strains to be sequenced and are considered reference strains. In addition, ATCC 33277 and W83 are frequently used to study mutants, in functional analysis, in adherence/invasion tests involving different types of human cells, and for other purposes.

For each strain, the mean coverage was approximately 30X (29.7 to 30.4), and the median corrected/trimmed read length obtained from canu was 6.3 Kb. The reads were assembled with canu, producing 18, 28, and 11 contigs for ATCC 33277, TDC60, and W83, respectively, and the finishing strategy described in the Methods section enabled assembly completion. All of the resulting complete assemblies differ from the initially published sequences. To validate the results, we performed PCR to confirm the organization of the large genomic multicopy regions. These included the rrn operons (Additional file 8: Figure S6a) and, for ATCC 33277, the duplication of the CTnPg1 genomic island and the orientation of the two copies (Additional file 8: Figure S6b-c). All of the RNA ribosomal operons from the ATCC 33277 and W83 strains were validated and found to correspond to the published structures. However, we corrected 3 of the 4 rrn in the TDC60 published genome, and this result confirms the observations made by Naito et al. in 2008 (see their Additional file 4: Figure S3) [100]. The most surprising result was the reorganization of 2 rrn loci in W83 compared to the other two strains (Additional file 8: Figure S6a). To determine whether this reorganization occurs only in this strain, 22 additional strains were subjected to PCR verification; we found that it is in fact restricted to W83 (Additional file 8: Figure S6d). The general arrangement of rrn operons in P. gingivalis chromosomes therefore seems remarkably stable despite the varied positions of these operons relative to the origin of replication (oriC) due to the high genomic mosaicism of the species. Homologous recombination at the rrn extremities seems extremely rare and probably accidental since of approximately 30 strains only W83 is affected.

Differences from published genome sequences

With all three circular genomes confirmed, we proceeded to identify the ways in which these genomes differed from the previously published sequences for these strains. We used progressiveMauve, which identifies locally collinear blocks (LCBs) to compare the published genomes with our constructions. Insertions and deletions in the de novo assembly compared to the published sequence were noted (Additional file 9: Figure S7). The differences are reported in detail below.

In ATCC 33277, CTnPg1 is found to be completely duplicated, with both copies oriented in the same direction (Additional file 9: Figure S7). The other differences in this strain are three deletions and one insertion (Additional file 10: Figure S8a). Since the ANI values showed a significant similarity of the ATCC 33277 and 381 strains and the 381 strain also has two complete CTnPg1 copies, we compared our ATCC 33277 genome reconstruction with the 381 genome [101]. We observed two collinear genomes with a length difference of only 1 kb and approximately 100 SNPs. The dissimilarities, which are minor, occur mainly in repeats and in the CRISPR1 region. As the complete 381 genome was assembled from 454 reads using Velvet and Newbler, there could be assembly errors. By mapping the Illumina HiSeq reads of the recent 381 sequencing [102] to our de novo ATCC 33277 construction, we positioned 99.5% of the reads without any gaps, suggesting that these American strains could be variants of the same strain.

TDC60 is the most changed, since, as discussed, 3 of the 4 rrn were incorrectly assembled during the initial construction. This means that in our new assembly, large sections of the genome are translocated (LCB2 and LCB5) and inverted (LCB3) (Additional file 9: Figure S7). The other differences are one insertion into a BrickBuilt 7 MITE [103] and four deletions (Additional file 10: Figure S8b).

Finally, the new W83 strain has a central inversion (LCB2) (Additional file 9: Figure S7). The new sequence consists only of a 523-bp insertion, which corresponds to eight additional direct repeats and the same number of new spacers in the CRISPR2 region (Additional file 10: Figure S8c). The genome length differences and numbers of SNPs in the three strains are shown in detail in Additional file 9: Figure S7.

Annotation and manual biocuration of three selected P. gingivalis strains

We automatically de novo annotated the three Porphyromonas gingivalis strains using three pipelines, after which they were manually biocurated. This allowed us to standardize the structural (syntactic) annotation, which consists mostly of CDS start/stop codon positions, and the functional annotation, using a single ontology for gene names and functional descriptions. The number and nature of all 12 rRNA genes (in four operons), 53 tRNA genes, and 1 tmRNA gene are identical to those in the published genomes, although the positions of some relative to oriC vary due to the reassembly. Seven of the riboswitches (5 cobalamin and 2 thiamine pyrophosphate riboswitches) were already annotated in all of the strains, and we added 1 S-adenosyl methionine (SAM)-II long-loop riboswitch per strain. It is noteworthy that all of these riboswitches were described by Hovik et al. [104]. For the small non-coding RNA (ncRNA), only rnpB (bacterial RNAse P) is present in the previous annotation. We positioned and annotated additional ncRNAs: one each ctRNA (antisense RNA), Bacteroidales-1 RNA, and bacterial signal recognition particle (SRP RNA) for each strain. We also added various group II catalytic introns: 5 in ATCC 33277; 4 in TDC60; and 3 in W83. For all riboswitches and ncRNA, the synteny is conserved. Two group II catalytic introns are present in haeR and in punA. The variation in numbers in this group is explained by the paralogy in traE (a gene containing this ncRNA is present once in W83 and twice in ATCC 33277 and TDC60) and by an additional ncRNA in the ATCC 33277 traG gene that is absent from the other two strains.

For ATCC 33277, TDC60, and W83, we added 213, 199, and 188 peptide signals and 23, 19, and 21 mobile elements (genomic islands, transposons, and conjugative transposons), respectively, to the original annotations. We further completed the annotation of intergenic sequences with repeated elements: BrickBuilt/MITEs (complete or partial), dispersed genomic repeats (> 500 bp with at least 95% identity), and CRISPR regions. We did not annotate any short repeats (3 to 31 bp) or sequence-tagged sites.

Finally, the longest and probably the most important biocuration work involved the CDS and pseudogenes, which we were able to dramatically improve. To facilitate traceability in publications and databases, we kept the initial locus_tags when the DNA sequence was unchanged or very similar (only having a new start and/or stop position). When the differences were larger (e.g., changes in the coding strand, ORF fusion, or a new CDS), we created new locus_tags marked PGN_n, PGTDC60_n, and PG_n for ATCC 33277, TDC60, and W83, respectively. In the previous versions, the pseudogenes had coding sequences and gene annotations, with pseudogene status annotated in the gene qualifier “note.” These could be classified into the categories “frameshifted,” “incomplete,” and “internal stop.” Following the NCBI Prokaryotic Genome Annotation guidelines (https://www.ncbi.nlm.nih.gov/genbank/genomesubmit_annotation/), we annotated the pseudogenes with the gene feature only (no CDS) and added an asterisk to the gene name field so that pseudogenes could be easily distinguished from normal coding regions. The gene_desc qualifier field describes the function of the pseudogenized gene, and the reason for the pseudogenization (“fragment” or “frameshift”) is listed in the note qualifier. Where they are relevant, gene fusions and coding strand changes are also mentioned in the “note” field.

To compare the original and curated annotations, we began by summarizing the common CDS and genes in the two versions as well as the number of curated start and stop positions and eliminated ELFs (Table 3). Additionally, for ATCC 33277, we added 55 new coding sequences; these usually corresponded to the re-establishment of a complete CTnPg1 copy or the reconstruction of complete ISPg CDS copies. For the TDC60 strain, we added 2 new pseudogenes (ISPg) and 6 new CDS (coding for 2 ISPg, an NAP, a peptidase, a rhodanase, and a PF07877 domain-containing protein). Finally, in the W83 strain, we de novo annotated 6 pseudogenes (5 ISpg and a PF07877 domain-containing protein) and 12 CDS (these encode 7 ISPg, 2 integrases, an NAP, a transcriptional regulator HxlR, and a virulence-related RhuM protein).

To compare the functional annotations in the two versions, we grouped the common features into five categories: nucleic acids and protein metabolism (replication, transcription, translation, histones, proteases, etc.); metabolism and transport of ions and other macromolecules (transporters, porines, ion channels, ferritin, kinases, hydrolases, etc.); mobile elements (conjugation, competence, phages, etc.); proteins without assigned biological functions but with a predicted subcellular localization or possessing an identified functional domain (lipoproteins, inner/outer membrane, or proteins which contain GLPGLI or zinc finger domains, etc.); and proteins that are conserved in all available P. gingivalis genomes but to which we could not assign a function (DUFs, FIGs, and some COGs, all marked as hypothetical).

Our biocuration and re-annotation work on the coding sequences and pseudogenes shared by the two versions resulted mainly in the de-anonymization of hypothetical CDS/proteins of unknown function (Fig. 5). These constituted approximately 22% of the original annotation and now represent only approximately 9%. The newly described functions were placed into four other categories as follows: 6% were added to macromolecule/ion transport and metabolism, an additional 3% were placed in nucleic acids/protein metabolism, mobile elements increased by 2%, and 2% were added to the list of CDS with a predicted conserved domain and/or subcellular localization. Finally, we noted that for the new pseudogenes and coding sequences as well as all of the other changes (features, fusions, separations, etc.), the new annotations are particularly enriched in mobile elements (Additional file 11: Table S4).

Pangenome analysis

All gene annotations (both true CDS and pseudogenes) from the three strains were classified into the five previously cited categories (Fig. 6). As shown in Additional file 12: Figure S9, with the exception of genes related to mobile elements, the absolute number of genes in each category is uniform. These transposases, integrases, phages, and conjugation genes are overrepresented in ATCC 33277, and underrepresented in TDC60. Even after manual biocuration, 16.8% of all genes remain poorly characterized; 7.7% have only subcellular localization predictions or functional domains, and 9.1% of the proteins conserved in P. gingivalis still have unknown functions.

• The P. gingivalis ATCC 33277, TDC60 and W83 pangenome can be divided into four categories

The first class is genes that are present in at least 2 copies in all of the studied strains; these represent 10.1% of the total number of genes in ATCC 33277 and 7.5% of the total genes in the other two strains. Some, such as DNA helicase, transcriptional regulators of the helix-turn-helix (HTH) 17 family, and nucleoid-associated proteins (NAP), are involved in nucleic acid metabolism. The others are associated with horizontal gene transfer and include integrases, tetracycline-resistance elements, tra genes of conjugative transposons, and transposases encoded by several types of insertion sequences (Additional file 13: Figure S10a). Transposases represent more than a third of the multicopy genes and are highly pseudogenized, particularly in W83 (Additional file 13: Figure S10b). Their distribution in the strains is variable; ISPg8 (previously known as ISPg1) occurs the most frequently, followed by IS195 (formerly ISPg3). It is notable that W83 is the only strain rich in ISPg4. In ATCC 33277, we annotated a pseudogenized ISPg5 despite this gene’s being described as absent by Califano et al. [109] (Additional file 13: Figure S10c).

The next category is the “dispensable/accessory genome,” which consists of genes that are present in only two of the three strains (70 genes in ATCC 33277, 82 in TDC60, and 49 in W83). Genes with unknown function are more abundant, with 40% (38.6 to 44.9%) poorly characterized genes in the accessory genome vs. 17% in the core (Fig. 5). Many of these genes are in chromosomal regions that have been annotated as genomic islands.

The third gene class (“strain-specific or unique genes”) consists of genes that are present in only one of the three strains. This class includes 3.3% of the genes in ATCC 33277, 2.9% of the genes in TDC60, and 4.7% of the genes in W83. Within this class, genes related to mobile elements per se are rare (< 5%), and genes of unknown function represent approximately 45% in ATCC 33277 and approximately 35% in TDC60 and W83. The singularity of these “unique” genes is relative since BLAST analysis of Porphyromonas genomes yields homologs in at least one other genome. This is often the case for ATCC 33277 genes in the 381 or HG66 strains and for W83 genes in the A7A1-28, A7536, and AJW4 strains. Some of the unique genes present in the three strains have hits in Tannerella forsythia genomes. Many of the unique genes are clustered in the genome.

The final gene class can be referred to as the “core” genome because it contains the genes that are present in all three strains. This class includes more than 80% of all genes (82.8% in ATCC 33277, 85.0% in TDC60, and 84.7% in W83) and can be further subdivided into a constant core, a variable core, a paralogy core, and a pseudogenized core (Additional file 14: Figure S11a). Due to paralogy, the total number of genes in the core genome varies, but it starts at 1522 genes.

The first subgroup of genes in the core genome is the “core genome sensus stricto” or the “constant core genome,” which contains genes with only one ortholog in each strain and more than 97% nucleotide identity. The constant core represents the majority of the core genes (n = 1424, 93% of the core); we classified these genes into the functional categories described above (Additional file 14: Figure S11b). There are only four genes related to horizontal gene transfer, and all four (pir, yhaI, vrgG, and p12) correspond to phage-related proteins. Since no phage infecting P. gingivalis has been described [110], it would be interesting to study the significance of this finding. We also observed four pseudogenes that are present in the three strains: COG4335, a hypothetical protein with a conserved domain associated with DNA alkylation activity [111]; pspB, which encodes alpha-ribazole-5’-phosphate phosphatase, an enzyme involved in coenzyme B12 biosynthesis; pyrD, coding for dihydroorotate dehydrogenase, which is involved in the de novo biosynthesis of pyrimidine; and finally tetR, which encodes a transcriptional regulator.

The second subgroup is the “variable core”, consisting of genes that also only have one ortholog but display less than 97% nucleotide identity (ranging from 96.9 to 66.3%). This subgroup represents only 3% of the core or 47 variable genes (Additional file 14: Figure S11c).

The third subgroup is the “core with paralogy”; it contains genes with only one ortholog in a strain but with duplications in at least one of the strains. This subgroup represents only 2.5% of the core genes (n = 36). The genes in this subgroup are related to DNA transcription, modification and repair (alkD, betI, btr, and ytxK), and transport (irtA, ndvA, and ydfJ). There are also four genes of known function (ctpA protease, epsJ glycosyltranferase, era GTPase, and the rhuM virulence protein) and four hypothetical proteins (Additional file 14: Figure S11d).

The fourth and final subgroup is the “core with pseudogenization”, which includes genes with only one ortholog in each strain but also at least one pseudogene in another strain. This group contains 36 genes, represents 2.5% of the core genes, and includes genes with unknown function as well as genes involved in adaptation and/or pathogenicity, for example, genes involved in transport, nucleic acid metabolism and cellular appendages (Additional file 15: Table S5).