An overview of our study design is provided in Figure 1. We analyzed 48 Prevotella isolates (Table S1), all taken from different sites from 35 subjects: 14 P. intermedia isolates, derived from sites with periodontal disease (Pi-disease group); 6 P. intermedia isolates from healthy sites (Pi-health); 17 P. nigrescens isolates from diseased sites (Pn-disease); and 11 P. nigrescens isolates from healthy sites (Pn-health). In some cases, the same subjects contributed samples of healthy and diseased areas.

Genome sequencing data were generated from four pooled samples, as described below. The general genome features of P. intermedia and P. nigrescens samples from diseased and healthy sites are listed in Table 1. About 80% of reads of clinical samples of P. intermedia could be mapped to the P. intermedia-17 (GenBank number: NC_017860.1/NC_017861.1) reference genome, whereas the remainder had best matches in various species, such as other Prevotella species, Bacteroides species, P. gingivalis, T. forsythia, Capnocytophaga species, Eubacterium species, Flavobacterium species, Alistipes species, and Fusobacterium species. Major differences between Pi-disease and Pi-health sample groups were that partial reads in Pi-disease can be mapped to Acidithiobacillus ferrivorans, Anaerofustis stercorihominis, Mycoplasma fermentans, Ruminococcaceae species and Treponema sp., while some reads in Pi-health mapped to Enterobacter cloacae and Xanthomonas fragariae.

Regarding the P. nigrescens samples, ~20% of the reads mapped to species other than the reference strain P. nigrescens ATCC33563 (GenBank number: NZ_AFPX00000000.1), including Prevotella species, Bacteroides species, P. gingivalis, T. forsythia, Alistipes species, Odoribacter splanchnicus, Neisseria species, Parabacteroides distasonis, Flavobacterium species, Tannerella species, Treponema species, Streptococcus agalactiae, and Faecalibacterium prausnitzii. Major differences between Pn-disease and Pn-health samples were that some reads in the Pn-disease group can be mapped to Bergeyella zoohelcum, Chryseobacterium species, Clostridiales species, while some reads in healthy group mapped to Acinetobacter baumannii and Saccharomyces cerevisiae.

Open reading frames (ORFs) likely to encode proteins were predicted according to NCBI Prokaryotic Genomes Annotation Pipeline (Tatusova et al., 2016). Then all predicted proteins were analyzed for sequence similarities against the Uniprot database (The UniProt Consortium, 2015). The total predicted ORFs in each group were: Pi-disease: 8020; Pi-health: 6255; Pn-disease: 6943; and Pn-health: 7710. There were 813 orthologous genes shared by the Pi-disease, Pi-health, and reference sequences (Figure 2A). In total, 3210 genes were shared between the Pi-disease and Pi-health groups, whereas Pi-disease samples had more group-specific genes than Pi-health samples (1987 vs. 1049). Totally 794 orthologous genes were shared by the Pn-disease, Pn-health, and reference sequences. Both diseased and healthy groups of P. nigrescens had almost the same numbers of group-specific genes (2344 vs. 2263), and shared 3247 genes between Pn-disease and Pn-health groups (Figure 2B).

Based on the lengths of matched segments in sample sequences, the percent identity of matched segments to reference, the full length of the matched contig, the degree of coverage of matched segments to reference, and the degree of coverage of matched segments to contigs, we investigated inter-sample diversity by performing a principal component analysis (PCA) of Pi sample groups (Figure 2C) and Pn sample groups (Figure 2D). Between Pi-disease and Pi-health samples, most contigs could not be differentiated (Figure 2C), indicating that overall there were no obvious differences in genome sequences between the groups that matched the reference. However, two additional small clusters of samples were clearly distinct (Figure 2C), consisting of contigs from each group. The contigs in the two distinct clusters were identified in Table S2, which were all belonged to the predicted insert sequences of Pi-disease or Pi-health samples. We also identified the annotated genes in these distinct contigs, most of which were conjugal transfer proteins (e.g., TraG, TraA, TraB, TraD etc. for Pi-disease group; TraE, TraK, TraM, and TraN for Pi-health group), phage associated proteins (both groups), recombinase (both groups), peptidase (Pi-health group), DNA-binding protein (Pi-disease group), and tetracycline resistance genes (e.g., tet Q and tet R; for Pi-disease group). Regarding the Pn group, most contigs from Pn-health samples clustered with those from Pn-disease; however, a considerable proportion of contigs in the Pn-disease group were clearly distinct from those in the Pn-health group (Figure 2D).

Unique Pi genes are listed in Tables S3, S4. The Pi-disease group had more unique genes involved in DNA modification, DNA recombination, DNA repair, and signal transduction, and more virulence-related factors (Figure 3A). Moreover, there were differences in exopolysaccharide biosynthesis and membrane proteins between the Pi-disease and Pi-health groups. Unique Pn genes are listed in Tables S5, S6. The dominant genes differentiating the groups were related to DNA modification, phage proteins, DNA recombination, and DNA repair (Figure 3B). The Pn-disease group had many more specific genes that functioned in DNA modification, and more specific phage-related proteins were found in the Pn-health samples, but both groups had almost the same number of unique genes related to virulence. There were also differences in exopolysaccharide biosynthesis and membrane proteins between the groups.

To investigate functional differences in orthologous genes between the Pi groups, we detected single-nucleotide variants (SNVs) in coding sequences (CDSs), and screened those whose functions were predicted to be highly affected (Figure 3C). The most affected regions were those encoding CDSs involved in DNA binding, DNA modification, and recombination, including HsdM, HsdS, XerC, XerD, and Type I restriction modification DNA-specific domain protein. The dominant regions that were affected by SNVs in the Pn-health (vs. Pn-disease) group were those encoding CDSs involved in ATPase, DNA modification (HsdS, HsdR), DNA recombination (XerC, XerD), conjugate transposon (TraD, TraN, TraI), and virulence protein (Figure 3D).

We also identified unique sequences in each sample, and tried to find any discrepant virulence gene candidates using BLASTX searches against the virulence factor database (VFDB; last accessed Nov 25, 2016; Chen et al., 2012) under a threshold e-value (≤ 1e-5). The numbers of unique sequences and matched virulence factors in each group are listed in Table S7. Most virulence factors specific to disease groups were capsule and lipopolysaccharide (LPS) synthesis-related factors. Moreover, disease groups had more unique sequences related to secretion systems, where they and the diverse proteins they secreted represent critical determinants of competitive fitness and pathogenic potential. The most abundant specific gene cluster in the Pi-disease group was the Type VI secretion system (T6SS). Pn-disease samples had specific sequences mapped to T3SS. Both the Pn-disease and Pn-healthy groups had T4SS effectors, but with different sequences. BspA and HlyB existed in all groups (health as well as disease), but with unique sequences in each sample group.

Furthermore, disease groups had more specific virulence factors annotated as proteinases and toxins. Only Pi-disease and Pn-disease samples possessed specific sequences that could match proteasome-associated ATPase (mpa gene).

We performed a comparison of metabolic networks between combined disease (Pi, Pn- disease) and combined health samples (Pi, Pn- health). We first determined the functional gene orthologs (KOs) in metabolic and regulatory pathways that shared by Pi-disease and Pn-disease (disease-share), and by Pi-health and Pn-health (health-share). Then we identified the differential KOs involved in networks between disease-share and health-share. We also managed to define special metabolic and regulatory characteristics of each sample group. (Figures 4, 5). The metabolic pathway of Pi-disease is highly consistent with that of Pn-disease, except for a few differences that mainly involved in energy metabolism, carbohydrate metabolism, and amino acid metabolism. In regulatory pathways, the major differences were seen in signal transduction (Figures 4A, 5, Figure S1) Pi-health and Pn-health shared all the KOs so they had exactly the same metabolic and regulatory pathways (Figure 4B, Figure S1B). Disease-share had specific KOs involved in the generation of 3-hydroxy-3-methglutaryl-CoA through valine, leucine and isoleucine degradation (e.g., E.2.3.3.10), and amino acid metabolism (e.g., murI), while health-share had specific KOs involved in pathways including Sulfur metabolism (e.g., cysI), glycolysis/gluconeogenesis (e.g., adhP), amino acid metabolism (e.g., CNDP2), amino sugar and nucleotide sugar metabolism (e.g., neuA, neuB, GME) (Figures 4, 5A). In addition, there was a distinct difference in signal transduction between disease-share and health-share. Both had unique two-component systems: e.g., algR and arcB in diseased samples and evgS, rcsC, and barA in healthy samples (Figure 5B, Figure S1).

We identified the KOs that merely exist/absent in Pi-disease or Pn-disease samples as well. Compared to other three sample groups, Pi-disease samples had gained specific genes involved in carbohydrate and lipid metabolism (e.g., lacZ, crtB, idi, and glmM), metabolism of cofactors and vitamin (e.g., arfC, UBIG, pncA, and lplA) and amino acid metabolism (e.g., speB, dapE, and E6.4.1.4B) (Figures 4, 5), while Pn-disease samples had extra genes with specific pathways including lipopolysaccharide biosynthesis (e.g., gmhB), pentose and glucuronate intercoversions (e.g., DCXR), and ubiquinone biosynthesis (e.g., menF). In regulatory pathways, Pi-disease samples had more unique genes involved in two-component system (e.g., phoN, ampC, and cikA) and membrane transport (e.g., ABCB10). Pn-disease had more unique genes involved in mismatch repair (e.g., MSH6) and two-component system (comA) (Figure 5). On the other hand, Pi-disease samples have also lost several genes, which involved in metabolism of cofactors and vitamin (e.g., K00077 apbA, K01719 hemD, K01495 folE, and K00652 bioF), fatty acid metabolism (e.g., K02372 fabZ, and K09478 ACADSB), and nucleotide metabolism (e.g., K01588 purE and K01493 comEB). There was only one gene ortholog that absent in Pn-disease samples only, which is involved in lipopolysaccharide metabolism (K02535 LPXC).

To understand the functions of different genes among the disease-derived (combined Pi,Pn- disease) and healthy-derived (combined Pi, Pn- health) isolates, protein functions were examined for clusters that were unique to isolates from diseased and healthy sites, based on the clusters of orthologous groups (COG) database (Tatusov et al., 2001) (Figure 6A). Diseased samples were enriched in 17 COGs, mostly related to gene regulation, signal transduction, biosynthesis, and function of cell wall structure, cell growth processes, and transport or metabolism functions, whereas healthy samples were enriched in only 1 COG related to energy production and conversion. The specific genes are listed in Tables S8, S9 and Figure 6. Diseased isolates had more unique factors related to cell wall and outer membrane structure and function, such as teichuronic acid biosynthesis glycosyltransferase (TuaC), D-alanyl-lipoteichoic acid biosynthesis protein (DltB), and the outer membrane cobalamin receptor protein. Others included the lipid A export ATP-binding/permease protein (MsbA), acetyltransferase (EpsM), lipid-A-disaccharide synthase, and dTDP-4-dehydrorhamnose 3,5-epimerase (dTDP_sugar_isom), which are involved in the pathway of LPS biosynthesis (Figure 6B).

Regarding microbial growth, disease-derived samples had more unique factors that influence cell division, cell metabolism, cell morphology, cell motility, and exopolysaccharide biosynthesis (Figure 6C). Diseased samples of both species contained an unusual YpdA/YpdB two-component regulatory system involved in a carbon control network whereas healthy isolates did not. In particular, isolates from diseased groups had a unique anaerobic ribonucleoside-triphosphate reductase activating protein, which is activated under anaerobic conditions, whereas healthy groups had a unique fumarate hydratase (class I, aerobic) that functions as an aerobic enzyme, involved in the citric acid cycle. Moreover, isolates from both diseased and healthy sites possessed different DNA methylases (Tables S8, S9). Diseased samples had more unique gene clusters that could aid resistance to stress factors, such as arsenate reductase and arsenic-resistance protein, which can enhance resistance to arsenate; heat-shock protein, which may help cells recover from stress; Na⁺/drug antiporter; and cold-shock DEAD box protein. Isolates from both types of isolates had different FeS cluster assembly proteins including SufB, which helps facilitate iron uptake from extracellular iron chelators under iron limitation.

LPSs are important factors involved in the virulence of Prevotella species, and we found that both diseased and healthy samples had several different factors directly or indirectly related to LPSs. Thus, we carried out a comprehensive analysis of LPS biosynthesis proteins in type culture strains from the NCBI database and clinical strains from each group of this study by reconstructing a neighbor-joining phylogenetic tree with MEGA 6 (Tamura et al., 2013) (Figure S2). LPS biosynthesis proteins from P. intermedia and P. nigrescens isolates clustered together. Notably, four disease-derived isolates clustered within a highly distinct and uniform clade.

Each group had many mobile genetic elements. We counted the number of conjugate transposons (CTns), transposases (Tns), and integrases, and divided by the number of isolates contained in each group to calculate the mobile genetic elements (MGE)/isolates ratio (Table 2). The ratio of CTns to isolates in Pi-health was higher than that in Pi-disease group, while the ratio in Pn-disease and Pn-health groups were similar. The ratios of Tns and integrase to isolates were higher in healthy groups (both in Pi-health and Pn-health) than that in disease groups. We also found that each group contained several prophages or degraded prophages, with Pn-disease and Pi-health samples having more than Pi-disease and Pn-health samples.

The numbers of spacers in Pi-disease, Pi-healthy, Pn-disease, and Pn-healthy samples were 577, 183, 980, and 522, respectively. In Pi-disease and Pn-disease samples, 4.7 and 4.8% spacers showed significant nucleotide similarity to the phage database, slightly higher than for Pn-healthy (4.4%) and Pi-healthy (2.7%) samples. More spacers matched the plasmid database in Pi-disease (3.3%) and Pi-healthy (3.3%) samples than in Pn-disease (2.7%) and Pn-healthy (2.9%) samples (Figure S3).

Patients with chronic periodontitis were recruited from the Department of Periodontology, Peking University School and Hospital of Stomatology. Enrollment occurred between December 2013 and May 2015. The inclusion criteria were: good general health and not pregnant, aged 18 to 65 (male or female), had not taken any antibiotic or antimycotic compounds in the past 3 months or undergone any periodontal therapy within 1 year, no aggressive form of periodontal disease, and periodontal disease sites at least 5 mm probing depth and attachment loss >0 mm. Subjects were excluded if they were pregnant, had any systemic condition which could affect the progress of periodontal disease. Periodontally healthy sites require a probing depth not >4 mm and no attachment loss.

The study protocol was approved by the institutional review board of Peking University School and Hospital of Stomatology (Beijing, China) (approval number: PKUSSIRB-2012063).

For each subject, subgingival plaque samples were collected from gingival pockets with a sterile Gracey curette from the bottom of the pocket to the coronal portion of the pocket. Samples were placed in individual tubes containing 500 μL Ringer's solution and then inoculated onto trypticase soy agar supplemented with 5% sheep blood, 1 mg/mL yeast extract, 1 μg/mL menadione, and 5 μg/mL hemin. Plates were grown anaerobically (90% N2, 5% H2, 5% CO2) at 37°C for 5 days, and black-pigment clones were selected and subcultured for purity. P. intermedia and P. nigrescens isolates were identified through full-length 16S rDNA sequencing (Sangon Biotech, Shanghai, China). Briefly, the full-length 16S rDNA region was amplified with a universal primer set (27forward, 5′-AGAGTTTGATCCTGGCTCAG-3′; 1492reverse, 5′-GGTTACCTTGTTACGACTT-3′). The PCR condition is as follows: 5 min initial denaturation at 95°C; 35 cycles of denaturation at 95°C (30 s), annealing at 60°C (30 s), elongation at 72°C (90 s); and final extension at 72°C for 10 min. The resulting 16SrRNA sequences were searched against the National Center of Biotechnology Information (NCBI) database, and assigned to species based on the best match at greater than 97% nucleotide identity over at least 95% length of the query.

Genomic DNA of each isolate was extracted using the Qiagen DNA Purification Kit following the protocol. After checking quality, genomic DNA was pooled according to the study design, and then a sequencing library was built and sequenced on an Illumina HiSeq 2500 system (Biotechnology Corporation, Shanghai, China). Approximately 23–40 million paired-end reads (2 × 125 nt multiplex) were obtained, representing >100-fold genome coverage. Genome mapping to reference genomes was performed using the Burrows-Wheeler Aligner software (ver. 0.7.12). Contigs were assembled de novo using SPAdes-3.5.0 (Bankevich et al., 2012), and filtered by deleting contigs with coverage ≤ 10% and length < 500 bp. Then we performed a cluster analysis of remained contigs using CD-HIT version 4.6 (Fu et al., 2012). Contigs with similarity ≥90% were compared, and only the longest one reserved. To understand the relationship between the contigs and the reference genome, we aligned the contigs to the reference genomes (P. intermedia samples to P. intermedia 17 and P. nigrescens samples to P. nigrescens ATCC 33563), and the lengths of matched segments in sample sequences, the percent identity of matched segments to reference, the full length of the matched contig, the degree of coverage of matched segments to reference, and the degree of coverage of matched segments to contigs were recorded for principal component analysis.

The draft genome was annotated using the NCBI Prokaryotic Genomes Annotation Pipeline (Tatusova et al., 2016). The remaining un-mapped reads of four samples were searched by BLASTN program against NCBI non-redundant nucleotide database with an E-value threshold of 1E-5.

To determine the unique genes in each sample group, we first performed alignments using the pan-genomes analysis pipeline (PGAP) among those four sample groups. The chosen standard was GF, score: 40, identity: 0.6, and coverage: 0.8. We further carried out a group-to-group BLASTN search. For example, the genes obtained from Pi-disease sample group were aligned with the genome sequences in Pi-health sample by BLASTN search, and the results were further filtered by deleting genes with identity > 30% and coverage > 30%, then we obtained a unique genes list for Pi-disease sample group.

To identify unique sequences in each sample, the sequences in each sample group were compared with all the other sample groups through MUMmer 3.23(Kurtz et al., 2004) (for example, Pi-disease vs. Pi-health; Pi-disease vs. Pn-disease; and Pi-disease vs. Pn-health). The sequences (≥500 bp) that existed in Pi-disease but absent in any other groups were considered as label sequences of Pi-disease group. Label sequences in each group were further annotated by BLASTX searching against the virulence factors database (VFDB, last accessed: November 25, 2016; Chen et al., 2012) with a word size of 10 under the thresholds of the highest bit score and e-value (≤ 1e-05).

Clustered regularly interspaced short palindromic repeat (CRISPR) regions were identified from contigs using the CRISPRs Finder online (Grissa et al., 2007). The sequences of CRISPR direct repeats and spacers in the four sequenced samples were extracted. To characterize a target of each spacer sequence, the spacer list was subjected to a BLASTN search with a word size of seven against the plasmid and the viral database, as described previously (Pride et al., 2012; Rho et al., 2012; Watanabe et al., 2013; Zhou et al., 2015; Wang et al., 2016). Hits were considered significant under the threshold of the e-value (≤ 1e-02).

The metabolic and regulatory pathways were reconstructed using 146 KEGG pathways and 22 KEGG regulatory pathways. An overview of the complete metabolism in biological systems was provided through the web-based interactive tool “Pathways Explorer” (iPath) (Yamada et al., 2011).

The assembled contigs were compared to reference genomes (P. intermedia 17 and P. nigrescens ATCC 33563) through MUMmer 3.23 (Kurtz et al., 2004), and sequences under the thresholds of length ≥500 bp and sequence identity < 90% were predicted as insert sequences. The predicted insert sequences were further searched using BlastX analysis in MvirDB database (MvirDB; last accessed April 21, 2012; Zhou et al., 2007) to investigate the possible mobile genetic elements (e.g., conjugate transposons, transposases, and integrases) (E value ≤ 1e-7). Prophages or degraded prophages were predicted using PHAge Search Tool (PHAST), which is a web server that available at (http://phast.wishartlab.com) (Zhou et al., 2011).

The raw reads have been deposited in the Sequence Read Archive under the following accession numbers: Pi-disease SRR5830884; Pi-health SRR5830885; Pn-disease SRR5830886; Pn-health SRR5830889.

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.