Comparative genome analysis and identification of competitive and cooperative interactions in a polymicrobial disease.

Polymicrobial diseases are caused by combinations of multiple bacteria, which can lead to not only mild but also life-threatening illnesses. Periodontitis represents a polymicrobial disease; Porphyromonas gingivalis, Treponema denticola and Tannerella forsythia, called 'the red complex', have been recognized as the causative agents of periodontitis. Although molecular interactions among the three species could be responsible for progression of periodontitis, the relevant genetic mechanisms are unknown. In this study, we uncovered novel interactions in comparative genome analysis among the red complex species. Clustered regularly interspaced short palindromic repeats (CRISPRs) of T. forsythia might attack the restriction modification system of P. gingivalis, and possibly work as a defense system against DNA invasion from P. gingivalis. On the other hand, gene deficiencies were mutually compensated in metabolic pathways when the genes of all the three species were taken into account, suggesting that there are cooperative relationships among the three species. This notion was supported by the observation that each of the three species had its own virulence factors, which might facilitate persistence and manifestations of virulence of the three species. Here, we propose new mechanisms of bacterial symbiosis in periodontitis; these mechanisms consist of competitive and cooperative interactions. Our results might shed light on the pathogenesis of periodontitis and of other polymicrobial diseases.

Introduction

Infectious diseases are the second leading cause of death in the world (http://www.who.int/whr/2000/en/). A polymicrobial disease is defined as an infectious disease caused by the mixed infection of several microbes. It is difficult to elucidate onset/progress mechanisms of a polymicrobial disease because this nosological entity is not eligible for the basic concept established by Robert Koch. Periodontitis (Ingham and Sisson, 1984) and respiratory illness (Hament ) are representative examples of polymicrobial diseases.

Periodontitis, which is known as a main cause of tooth loss (Papapanou ; Irfan ), leads to the destruction of periodontal ligaments and adjacent supportive alveolar bone. According to a report by the World Health Organization, 15–30% of 35- to 44-year-old human beings are afflicted by periodontitis (http://www.who.int/en/), and the prevalence of periodontitis is the highest of all the infectious diseases in the world. More than 300 existing bacterial species are implicated in the initiation and progression of periodontitis, especially the following three bacterial species: Porphyromonas gingivalis, Treponema denticola and Tannerella forsythia, which are classified as ‘the red complex'. They are detected together at a high rate in deep periodontal pockets (Socransky ; Darveau, 2010).

In bacterial communities, one bacterial species tends to obtain different functions from the other species (Scott ). The red complex species might interact with one another to adapt to the niche and thus may cause periodontitis in a cooperative manner. Many studies revealed the virulence factors of T. denticola and P. gingivalis, whereas the reports of T. forsythia virulence factors are rarer (Ishihara, 2010; Grenier and La, 2011; Dashper ; Bostanci and Belibasakis, 2012). Some reports are suggestive of synergistic effects of virulence factors between T. denticola and P. gingivalis; T. denticola does not form biofilms by itself in vitro (Vesey and Kuramitsu, 2004), but coculture of P. gingivalis and T. denticola results in biofilm formation and coaggregation (Kuramitsu ). On the other hand, T. forsythia grows slowly and requires special precise conditions for its growth (Honma ; Sakakibara ). The information about T. forsythia is limited, and there is no genome-wide analysis to identify any interactions among the red complex species.

Major research projects on the diagnosis of polymicrobial disease are mainly based on the metagenomic analyses, which can identify a wide range of bacteria from various environments (Brogden ). Nonetheless, to the best of our knowledge, no data in the literature describe simultaneous comparative genome analysis and focuses on the specific bacterial species that are strongly related to polymicrobial diseases. In this study, we performed comparative genome analyses among the red complex species and demonstrated reciprocal relationships among the three species: competitive and cooperative relationships. Our results are suggestive of significance of competitive and cooperative interactions in the periodontitis as a novel mechanism of bacterial symbiosis.

Materials and methods

Detailed information is described in Supplementary Materials and methods.

Bacterial strains and culture conditions

We used 18 T. forsythia strains that were isolated at Tokyo Medical and Dental University, Tokyo, Japan.

Determination of T. forsythia complete genome sequences

Complete genome sequences of two T. forsythia strains (KS16 and 3313) were determined using a combination of the GS Junior Titanium instrument (454 Life Sciences a Roche Company, Branford, CT, USA) for shotgun/paired-end sequencing, Genome Analyzer IIx (GAIIx; Illumina, San Diego, CA, USA) for 101-bp paired-end sequencing and traditional Sanger sequencing (ABI 3730; Applied Biosystems, Foster City, CA, USA). Then, these reads were assembled in Newbler v.2.6 to obtain scaffolds.

Determination of T. forsythia draft genome sequences

Draft genome sequences of 16 T. forsythia strains (sit3, 15, 20, 2612, MH6, L7, 273, KM4, 3322, 291, 3114, TR1, 2444, 1224, 2442 and K8) were determined using GAIIx paired-end sequencing. Then, the draft genomes were annotated using Rapid Annotation using Subsystem Technology server v.4.0 (http://rast.nmpdr.org/; Aziz ).

Characterization of T. forsythia species by genome information

We characterized 19 T. forsythia strains by estimating the pan-genome size, by dotplot analysis and via construction of phylogenetic tree, as described below in the methods for the three species of the red complex. The findings of these experiments are described in Supplementary Results.

Estimation of pan-genome size

The gene sequences were clustered among the strains of each species by PGAP v.1.11 (Zhao ). Following the clustering of the genes, we estimated the number of new genes when sequencing more genomes and the pan-genome size for each species. Regression analysis for new genes and the pan-genome was conducted as described previously (Tettelin ).

Multiple dotplot analysis

Multiple dotplot analysis was conducted as described previously (Watanabe ).

Construction of phylogenetic tree and single-nucleotide polymorphism analysis

In each species in the red complex, the core genes that were single copy in all the strains were tested to see whether rearrangement in a gene was significant, by using the Phi test (Huson and Bryant, 2006). Then, the concatenated amino-acid sequences from all the single-copy core genes were used for construction of phylogenetic trees. For all the gene clusters in each species where a gene was located in more than two genomes, a dN/dS ratio was calculated by means of codeml in the PAML v.4.6 (Yang, 1997).

CRISPR analysis

Clustered regularly interspaced short palindromic repeat (CRISPR) regions and CRISPR-associated (cas) genes were identified, and a target sequence of each spacer was characterized as described previously (Watanabe ).

Metabolism analysis

In each species, the core genes were mapped on metabolic pathways using KAAS (http://www.genome.jp/kaas-bin/; Moriya ), and were visualized using the iPath v.2 (Yamada ).

Virulence gene

To test for the presence of known virulence genes in all the genomes of strains in each species, BLASTP was used for all the genes of 38 strains (i.e., the sum of the three species strains) against the virulence genes described in the previous reports (Holt ; Sharma, 2010; Dashper ). In addition, we tried to find any novel virulence gene candidates using BLASTP against the virulence factor database (last accessed 25 Novemer 2013; Chen ) and the microbial database of protein toxins, virulence factors and antibiotic resistance genes for biodefense applications (MvirDB; last accessed 21 March 2007; Zhou ).

Results and discussion

Comparative genomics of the red complex species

The red complex species are reported to persist together in the periodontal pocket (Socransky ; Paster ; Kumar ; Griffen ; Casarin ; Ge ), and other oral bacteria are not always detected (Supplementary Table S1). Nevertheless, there are no reports in the literature about the precise relationship among the red complex species. Comparative genome analysis among the red complex species may be a significant clue not only to the characteristic features of each species but also to the interactions among them. For the comparative genome analyses, we applied the same annotation criteria (see Materials and methods) to all the strains of the red complex species including the publicly available genome sequences (Supplementary Tables S2 and S3).

Gene content

An average genome size of T. forsythia, T. denticola and P. gingivalis was 3.38, 2.83 and 2.34 Mbp, respectively; the genome size of T. forsythia was the largest in the red complex species (Supplementary Table S3). In the complete genome sequences, the average number of coding DNA sequences (CDSs) of T. forsythia, T. denticola and P. gingivalis were 3005, 2693 and 2175, respectively. Within each species, all the strains had the same numbers of rRNA and tRNA genes; the number of rRNA operons of T. forsythia, T. denticola and P. gingivalis were 2, 2 and 4, respectively, and the number of tRNA genes of T. forsythia, T. denticola and P. gingivalis were 46, 45 and 54, respectively. The core genes (shared within all the strains in the species) and the species-specific genes of T. forsythia, T. denticola and P. gingivalis were 1733/1130, 1966/1899 and 1610/1005, respectively (Supplementary Figure S1 and Supplementary Tables S4–S6). The number of core genes accounted for 60% (1733/2912), 73% (1966/2693) and 74% (1610/2175) of average number of CDSs in 19 T. forsythia strains, 5 T. denticola strains and 3 P. gingivalis strains, respectively. Forty-five genes were shared among the red complex species, including housekeeping genes and molecular chaperones such as ClpB (Hsp104) and DnaK (Hsp70; Supplementary Table S7). A previous report indicated that these chaperones regulate overproduction of reactive oxygen species that is associated with the inflammatory response (Steeves ). When comparing the genes within each pair of species, the number of gene homologs was 548, 10 and 12 between T. forsythia and P. gingivalis, between T. forsythia and T. denticola, and between T. denticola and P. gingivalis, respectively (Supplementary Figure S1). The homologs shared between T. forsythia and P. gingivalis were 19% (548/2912) of all T. forsythia genes and 25% (548/2175) of all P. gingivalis genes.

Next, we estimated pan-genome size based on the 19 T. forsythia genomes, 14 T. denticola genomes and 5 P. gingivalis genomes to determine the genetic landscape (Supplementary Figure S2). First, we estimated the number of new genes that would be discoveredby using power-law regression analysis, n=κ·N−α (Tettelin ). We then estimated the number of pan-genomes that would be accumulated by using the same analysis, n=ɛ·Nγ (Tettelin ). The values were α=0.494 and γ=0.254 in T. forsythia, α=0.733 and γ=0.256 in T. denticola and α=0.834 and γ=0.246 in P. gingivalis, respectively. These values in red complex species were indicative of an open pan-genome because the exponent α was <1.0 and γ was >0. Recent studies reported that some bacterial species also show an open pan-genome, including Propionibacterium acnes and Corynebacterium pseudotuberculosis (Soares ; Tomida ). Although the pan-genomes of the red complex species tend to be open, these results suggested that almost 20% genes were homologs shared between T. forsythia and P. gingivalis. These results may be the consequence of the fact that P. gingivalis and T. forsythia belong to the same phylum Bacteroidetes.

Genomic collinearity

Genomic rearrangements have a significant effect on bacterial evolution, including the genome reductive processes and the creation of new DNA regions (Sun ). We observed only one or two continuous lines in the dot plots of T. forsythia and T. denticola. On the contrary, the dot plots of P. gingivalis contained many parallel and numerous discontinuous lines in its plot (Figure 1). Similar to the previous report (Naito ; Watanabe ), we identified frequent genome rearrangements in P. gingivalis, which contained inversions, translocations, simple insertions, deletions and replacements. Simple insertions, deletions and replacements were observed two times more frequently than inversions or translocations.

Figure 1

Dotplot analysis of T. forsythia, T. denticola and P. gingivalis. A dot is plotted when the 20-bp sequences exhibited similarity between two genomes. The dots are colored according to the similarity score, with gradation from 70% (blue) to 100% (red).

Because mobile genetic elements (MGEs) are intimately involved in the extensive genome rearrangement of P. gingivalis (Watanabe ), we next counted the number of MGEs: insertion sequences, conjugative transposons and miniature inverted-repeat transposable elements (Supplementary Figure S3). The number of each MGE per CDS was the lowest in T. denticola and highest in P. gingivalis among the red complex species. The ratio of insertion sequences/CDS of T. denticola was lower than that of the other two species, with statistical significance (P<0.01 for T. denticola and P. gingivalis; P<0.05 for T. denticola and T. forsythia; Student's t-test).

These results suggest that the number of MGEs is high when extensive genome rearrangements are frequently observed. On the other hand, the number of MGEs is low when genome rearrangements are rare. The red complex species might be classified into two types, based on their genome structure: one possible type, including T. forsythia and T. denticola, has maintained stable genome structure, and the other type (represented only by P. gingivalis) had undergone frequent changes in the genome structure as a result of complicated genomic rearrangements mediated by transposition of MGEs, demonstrated in the previous reports (Naito ; Watanabe ).

The ratio of nonsynonymous substitutions to synonymous substitutions

To assess the dynamics of molecular evolution and to detect proteins undergoing adaptation to environmental changes in each species (Kimura, 1983; Gillespie, 1991; Ohta and Ina, 1995), we estimated the ratio of the non-synonymous substitution rates (dN) to the synonymous substitution rates (dS; the dN/dS ratio) for the core genes of each species. In T. forsythia, T. denticola and P. gingivalis, 5.4% (54/995), 0.7% (10/1 320) and 2.6% (11/425) of the core genes exhibited dN/dS>1, respectively, indicating positive selection (Figure 2 and Supplementary Table S8). The 54 genes in T. forsythia contained 26 hypothetical genes; among these hypothetical genes, no protein motif was found in 16 genes, whereas the other 10 genes encoded several functional domains including those in peptidase and outer membrane proteins. In T. denticola, all 10 genes with dN/dS>1 were hypothetical genes, including four genes that contained the domains such as the RHS repeat-associated core domains and cysteine protease domains. In P. gingivalis, 11 genes with dN/dS>1 included seven hypothetical genes, the nitroimidazole resistance gene and the SOS mutagenesis/repair gene. Among the hypothetical genes, six genes encoded several domains including TlpA-like family related to the cytochrome maturation and phosphopantetheinyl transferase superfamily related to the biosynthesis of butyric acid.

Figure 2

The dN/dS ratio of T. forsythia, T. denticola and P. gingivalis. The dN/dS ratios of all the gene clusters in this analysis are shown in an ascending order for (b) 19 T. forsythia, (b) 14 T. denticola and (c) 5 P. gingivalis strains, respectively. Red squares indicate the genes that correspond to the condition dN/dN>1.

As a conclusion from the above three analyses, we propose different evolutionary strategies among the red complex species; intraspecies genetic diversity is low in T. forsythia reflecting stable genome structure, and amino-acid sequence alterations may be a major driving force of diversification in T. forsythia. T. denticola shows the most stable genome structure and the lowest selective pressure among the three species. With respect to P. gingivalis, the main evolutionary strategy may be based on frequent genome rearrangement by transposing MGEs rather than amino-acid sequence variation.

Competitive interactions between T. forsythia and the other red complex species via CRISPR/Cas systems

Previous reports showed that the CRISPR/Cas system prevents acquisition of conjugative plasmids (Marraffini and Sontheimer, 2008), integrative conjugative elements (Lopez-Sanchez ) and environmental DNA via natural transformation (Bikard ; Zhang ). Similarly, the CRISPR/Cas system might affect the molecular interactions and the survival strategies among the red complex species.

The red complex species had different sets of CRISPR-associated genes (cas genes) and CRISPR repeat sequences (Supplementary Table S9). When we investigated the complete genome sequences of three T. forsythia strains, according to the most recent classification by Makarova , all three strains carried a type I-B CRISPR/Cas locus with six cas genes (cas6, cas7, cas3, cas4, cas1 and cas2; Supplementary Figure S4). Strain KS16 had another CRISPR/Cas locus, which was a type II CRISPR/Cas system with three cas genes (cas2, cas1 and cas9). We next identified the CRISPR/Cas loci in the 19 T. forsythia genomes including draft genome sequences. The CRISPR/Cas loci were detected in all 19 strains. Five strains (KS16, 2612, 291, 3114 and 3322) had both type I-B and II CRISPR/Cas systems, and strain TR1 had a type I-B and an unclassified type of CRISPR/Cas locus, whereas the other 13 strains possessed only a type I-B CRISPR/Cas system (Supplementary Table S9). In the type I-B CRISPR/Cas locus, a gene deletion was observed in six strains without differences in clinical parameters (sit3, KS16, 3322, 2612, 3114 and 2444; Supplementary Figure S4). Furthermore, a 3.5-kbp insertion was found between the csm4 and csm5 genes in 3313. The insertion region contained hypothetical genes and phage-encoded genes, and hereafter we call this region a ‘phage-like structure'. Such structures were found in five strains (3313, 2444, 1224, 2442 and TR1) and were not found in any regions other than this CRISPR region. In the phylogenetic trees (see Supplementary Materials and methods and Supplementary Results for the construction procedure), strain 3313 was phylogenetically distant from the four strains (2444, 1224, 2442 and TR1; Supplementary Figure S5), suggesting that strain 3313 obtained the phage-like structure via a different mechanism compared with the above-mentioned four strains.

On the other hand, the cas gene sets of T. denticola were highly conserved among the six draft genomes. The CRISPR/Cas locus in T. denticola was classified as a type II-A CRISPR/Cas system because of the genes csn2, cas2, cas1 and cas9 (Supplementary Figure S6). Strain H1-T had only one repeat sequence, whereas we identified no repeat sequences in strain H22. The structure of the cas gene in P. gingivalis genomes was also conserved at two CRISPR loci, including type I-C and III-B (Watanabe ).

Polymorphisms of the repeat sequences were observed among the strains in each species (Supplementary Figures S7 and S8 and Supplementary Tables S10 and S11), in line with a previous report (Watanabe ). We identified a few mutations at the end and in the middle of the repeat sequences of T. forsythia and T. denticola, respectively. Three types of T. forsythia CRISPR repeats had different secondary structures and were grouped into different categories based on superclasses, classification of the phylum, cas gene types and sequence families (Supplementary Figure S9). Protospacer-adjacent motifs were not clearly detected in all CRISPR types of T. forsythia and T. denticola, which may be because of the lesser amount of the spacers exhibiting high nucleotide similarity to the sequences in the databases (Supplementary Figure S10).

The total number of spacers and unique spacers in each species (T. forsythia, T. denticola and P. gingivalis) were 1631/1149, 78/66 and 267/244, respectively (Supplementary Table S12). Among these spacers, none was common for the three species. Dendrogram-based analysis of spacers showed that strains 3322 and 3313 were clustered as the same group, and this finding was consistent with the structure of phylogenetic trees of the core genes (Supplementary Figure S11). In P. gingivalis, CRISPR typing was shown to be capable of clustering the strains beyond the fimbrial gene types (Watanabe ); in this study, we demonstrated that the CRISPR typing was applicable to T. forsythia to show strain clusters that can be formed in the core gene-based phylogenetic tree.

A previous study reported that the CRISPR spacers in P. gingivalis possibly inhibit both genomic rearrangements and intercellular recombination among P. gingivalis strains (Watanabe ). Nonetheless, there are no studies focused on CRISPR targets in T. forsythia and T. denticola. To identify such targets of the CRISPR/Cas system in T. forsythia and T. denticola, we analyzed the spacers in these species by using nucleotide similarity searches against the seven databases (see Materials and methods). In T. forsythia and T. denticola, 106/1631 (6.5%) and 7/78 (9.0%) non-redundant spacers exhibited high nucleotide similarity to the sequences in the databases, respectively (Figure 3, Table 1, Supplementary Figure S12 and Supplementary Table S13).

Figure 3

Similarity searches of CRISPR spacers in T. forsythia. In the T. forsythia CRISPR spacers, the presence of the targets is shown as a pie chart, where the targets of the 106 spacers are shown as other charts. The targets of the 25 spacers against the bacteria other than the three red complex species are indicated by the species names.

Table 1

Sequences in the databases exhibiting high nucleotide similarity to unique spacers

Spacer ID					Species	Gene name
T. forsythia
546					T. forsythia ATCC 43037	Putative membrane protein
744					T. forsythia ATCC 43037	Hypothetical protein
750					T. forsythia ATCC 43037	Putative membrane protein
975					T. forsythia ATCC 43037	Phenylalanine and histidine ammonia lyase
1090					T. forsythia ATCC 43037	Hypothetical protein
1256					T. forsythia ATCC 43037	Recombination protein O
1416	1558	1683	1806		T. forsythia ATCC 43037	Hypothetical protein
1434	1436	1448			T. forsythia ATCC 43037	Bacterial group 2 Ig-like protein
1435	1437	1438	1440	1442	T. forsythia ATCC 43037	Bacterial group 2 Ig-like protein
1439					T. forsythia ATCC 43037	Bacterial group 2 Ig-like protein
1820					T. forsythia ATCC 43037	Putative nuclease sbccd, subunit C
1859					T. forsythia ATCC 43037	Hypothetical protein
157	1032				P. gingivalis TDC60	IS (12 bp terminal inverted repeat)
545					P. gingivalis TDC60	Conjugative transposon protein traq
578	756	1258			P. gingivalis TDC60	Hypothetical protein
765					P. gingivalis TDC60	50S ribosomal protein L10
769					P. gingivalis ATCC 33277	30S ribosomal protein S16
1165					P. gingivalis ATCC 33277	Hypothetical protein
1171					P. gingivalis ATCC 33277	Hemagglutinin protein haga
1182					P. gingivalis ATCC 33277	DNA-binding protein HU
1275					P. gingivalis ATCC 33277	Intergenic region
1389	1531	1656	1779		P. gingivalis ATCC 33277	Decarboxylating precorrin-6Y C5,15-methyltransferase
265					T. denticola ATCC 35405	Hypothetical protein
267	268				T. denticola ATCC 35405	Hypothetical protein
547					P. denticola F0289	Hypothetical protein
199					Actinomyces sp. oral taxon 178 str. F0338 contig00144
118	1146	1294			Bacteroidetes bacterium oral taxon 274 F0058
795	867				Canis familiaris chromosome 9
1309	1451	1576			Human chromosome 14 DNA sequence BAC R-662J14 of library RPCI-11
1310	1452	1577			Homo sapiens chromosome
1917					Acetobacterium woodii DSM 1030
T. denticola
229					T. denticola ATCC 35405	Hypothetical protein
251					T. denticola ATCC 35405	Hypothetical protein
265	266	267	268		T. denticola ATCC 35405	Hypothetical protein

Abbreviations: Ig, immunoglobulin; IS, insertion sequences.

In T. denticola, 6/7 spacers showed significant nucleotide similarity to hypothetical genes in the genome of T. denticola (Table 1, Supplementary Figure S12 and Supplementary Table S13). Four out of those six spacers were located only in the ATCC 35404 genome, and their targets were the genes in the ATCC 35405 genome that perfectly matched the nucleotide sequence. Meanwhile, the remaining two of the six spacers were present only in ATCC 35405 and targeted the genes on the same chromosome, with one single-nucleotide polymorphism present in their protospacer regions. This phenomenon can be called ‘self-targeting', in the sense that the CRISPR spacer targets the same chromosome endogenously. One study reported that 100 of 23 550 spacers (0.4%) are self-targeting spacers; these spacers are widely distributed over diverse phylogenetic lineages (Stern ). The self-targeting spacers may have influenced genome-scale evolution by participating in regulation of gene expression (Aklujkar and Lovley, 2010) or by being lethal for the strain harboring such spacers (Westra ). Generally, a strict nucleotide match (no more than 1-bp mismatch) between a CRISPR spacer and its target site is required for CRISPR interference (Sorek ), although it was recently reported that a maximal 5-bp mismatch is tolerated in type I if these nucleotides are outside the seed region in a protospacer (Semenova ). Because strain ATCC 35405 is culturable even harboring the self-targeting spacer, it is possible that ATCC 35405 may evade lethality via unknown mechanisms that do not involve the mismatch being in the target site and may use the CRISPR/Cas system for its evolution.

In T. forsythia, 55/106 nonredundant spacers showed high nucleotide similarity to the annotated sequences in the databases (Figure 3, Table 1 and Supplementary Table S13). Of these 55 spacers, 22 exhibited significant nucleotide similarity to T. forsythia's own genome, which contained outer membrane protein genes and bacterial group 2 immunoglobulin-like protein genes. Meanwhile, 25/55 spacers showed similarity to other bacterial genomes; 16/25 and 3/25 spacers exhibited significant nucleotide similarity to genomes of P. gingivalis and T. denticola, respectively. These results suggest that the CRISPR/Cas system of T. forsythia might target the genomes of other red complex species. In particular, the spacers of strains 2444, 1224, 2442 and TR1 showed significant nucleotide similarity to the methyltransferase gene of P. gingivalis. According to an investigation of the adjacent regions of the methyltransferase gene on P. gingivalis genomes, this gene was not likely to be associated with any transposable elements (data not shown). In addition, it is notable that the spacers of strains 3313, 3322 and 15 exhibited significant nucleotide similarity to several MGEs (conjugative transposons and insertion sequences) in P. gingivalis. In P. gingivalis, DNA transfer is reported via conjugation or natural transformation (Naito ; Tribble ). Competitive interactions between T. forsythia and the other two species in the red complex are possible; P. gingivalis delivers its DNA into the cells of other species to assert its dominance in the niche, whereas T. forsythia attacks the methyltransferase gene of P. gingivalis by delivering the spacer and related Cas proteins of T. forsythia CRISPR/Cas system into P. gingivalis cells, to impede their persistence.

Bacterial interactions may be important as a competitive interaction among the species that causes growth inhibition (Rao ), metabolic inactivation (Stringfellow and Aitken, 1995) and downregulation of genes (Xie ). We demonstrated a possibility of competitive interactions among the red complex species via the CRISPR/Cas system. The CRISPR/Cas system in T. forsythia may work as both a weapon toward the other members in the community and a barrier against them. Further research will be needed to prove our hypothesis, for example, by characterizing targets of the CRISPR spacers that are not characterized in this study or by demonstrating DNA/protein transfer among the red complex species.

On the other hand, we also demonstrated that evolutionary strategies were different among the three species living in the same niche; this observation is consistent with that in other environments (Venter ; Eckburg ) because environmental conditions may affect the evolutionary strategy of bacteria (Roszak and Colwell, 1987). The periodontal pocket is an environment that can be changing continuously; the environment is altered by medical treatments with the removal of subgingival calculus and biofilm to reduce pocket depth (Hellström ), although the red complex species tend to form a biofilm to protect themselves from deleterious conditions. The different evolutionary strategies of the red complex species may be essential for adaptation to this variable environment. One study suggested that P. gingivalis produces a number of surface molecules and structures to interact with T. denticola and T. forsythia (Holt and Ebersole, 2005). We propose different roles for P. gingivalis and its two partner species to persist in the periodontal pocket; P. gingivalis might establish an ecological niche in the pocket, followed by participation of T. denticola and T. forsythia carrying biofilm-related genes to form a community in which the three species are the main members. It is still an open question as to why they dominate other kinds of oral bacteria.

Cooperative interactions among the red complex species

Several relations between the metabolism and niche habitat of organisms have been shown by both experimental and in silico analyses (Kastenmüller ). Metabolic pathways are important for both viability of bacterial cells and expression of virulence, which is achieved by production of virulence substances via multiple metabolic reactions (Panos and Ajl, 1963). We mapped the core genes in each species in the red complex onto the metabolic pathways (Supplementary Table S14, see Materials and methods), and obtained the complementation ratio of the number of mapped genes to the total number of genes involved in the respective pathways.

All 29 main pathways exhibited an abundance pattern of the mean complementation ratio when the number of species that were included in the mapping analysis was increased (Supplementary Figure S13 and Supplementary Table S15). Notably, the complementation ratio of the fatty acid biosynthesis pathway was the highest among the 29 pathways for each number of the species (one species: 68.84±6.99% two species: 76.81±2.51% and three species: 78.26%, respectively). Fatty acids are known to stimulate bacterial growth (Grenier, 1992a) and to be associated with inflammatory conditions (via increased concentration), such as periodontitis (Kurita-Ochiai ; Niederman ) and anaerobic infections (Ladas ; Mills ). We examined the complementation patterns in the fatty acid biosynthesis pathway, and found that only the butyric acid metabolism pathway was the case where all essential genes were complemented in the three species (Figure 4 and Supplementary Figure S13). Moreover, the fold increase was the highest in the butyric acid pathway among the 29 main pathways described above. Our results are suggestive of a possibility that the biosynthesis of fatty acids is an essential metabolic pathway for the red complex species, and the fatty acid produced via cooperative interactions among the red complex species might stimulate growth of these species and development of periodontitis. We also noted that the pathways except for the butyric acid biosynthesis were incomplete even if the red complex species coexisted.

Figure 4

Metabolic analyses of T. forsythia, T. denticola and P. gingivalis. (a) The stand-alone pathways in the small pathway groups are colored if they are mapped by the core genes of either one, two or three species, as follows: gray (a known metabolic pathway), green (T. forsythia), blue (T. denticola), red (P. gingivalis), purple (T. forsythia and T. denticola), yellow (T. denticola and P. gingivalis), pink (P. gingivalis and T. forsythia), light blue (T. forsythia, T. denticola and P. gingivalis), respectively. (b) A fatty acid synthesis pathway group, which includes a butyric acid pathway, is shown by enlarging the corresponding region in (a). (c) In the butyric acid pathway, the gene required for each change of substance is shown, with the colors described in (a).

Butyric acid is reported to be responsible for progression of periodontitis because this substance affects the host immune system in the periodontal pocket (Kurita-Ochiai ; Kurita-Ochiai and Ochiai, 2010). Previous reports have shown that there were several pathways for the synthesis of the butyric acids. Metabolic pathways vary among Fusobacterium species with distinct pathways (Gharbia and Shah, 1991) and metabolic shift from butyrate formation to lactate and acetate formation in the fermentation is found to be associated with changes in the activities of several key enzymes (Zhu and Yang, 2004). However, the pathway for the synthesis of the butyric acids among the red complex species might not apply to any of them. Furthermore, the fatty acid cross-feeding has previously been demonstrated in P. gingivalis and T. denticola (Grenier, 1992b), which is consistent with the downregulation of three genes (fabG, fabF and acpP) encoding enzymes participating in the initial stage of fatty acid synthesis (Tan ). These results suggest that the presence of T. denticola helps P. gingivalis reduces energy-consuming processes (Tan ). Complementation of essential genes in the butyric acid metabolism pathway may be important for the three species to facilitate their persistence in the niche. Because of the hypothetical protein transfer among the three species as described the above section, proteins encoded by those essential genes may be secreted to extracellular space and may participate in sequential reactions with the intermediates. Alternatively, it is also possible that the intermediates are transferred from the extracellular space into bacterial cells, and vice versa, and undergo reactions inside a bacterial cell possessing a relevant protein. The intermediates may be transferred via endocytosis or exocytosis; an endocytosis-like process was reported in Gemmata obscuriglobus (Fuerst and Sagulenko, 2010), whereas an outer membrane vesicle can be an example of exocytic transfer (Haurat ).

Furthermore, one study showed that cell-to-cell bacterial communication affected their growth (Toyofuku ). Despite these scenarios, there is still a possibility that one chromosome harbors all essential genes in a metabolic pathway, which can be shown via functional assignment of hypothetical genes. To find them, we searched the genes that are not associated with metabolic pathways for the presence of any protein motifs from the pathways, in NCBI's Conserved Domain Database (CDD; last accessed 21 March 2013; Marchler-Bauer ). There were no genes with significant probability of participating in metabolic pathways. Our findings are based on the analyses focused on the three red complex species; a biofilm in a periodontal pocket harbors hundreds of other bacterial species, which can compensate the metabolic pathways described above. Nevertheless, because the species other than the red complex are not always detected in periodontal pockets (Supplementary Table 1), they may not be capable of compensation of the metabolic pathways.

No study has shown that mixed cultures of the red complex species produce higher amounts of butyric acid compared with sole cultures. However, several studies have shown that the red complex species showed a clear increase in abundance and prevalence with increase in the severity of periodontal diseases (Loozen ) and the concentration of butyric acid in periodontal pockets was correlated with the severity of the diseases (Margolis ; Botta ). Furthermore, one study has shown that the isobutyric acid was found to promote the growth of P. gingivalis and T. denticola when both were inoculated together (Grenier, 1992a). The possible biological significance of our findings is that the mutual compensation of the butyric acid metabolism genes could facilitate persistence and virulence of the three red complex species.

In addition to the analyses of metabolic pathways, we examined the presence of virulence genes in the genome of each species. Most core genes of T. forsythia species were found to be housekeeping genes, and more than 90% of known virulence genes were included in the core genes (Supplementary Information and Supplementary Tables S4 and S16). In addition, by searching for amino-acid sequence similarities between the core genes and the virulence factors in the database, we identified four novel virulence gene candidates, which encode L-asparaginase, a hypothetical protein, an ATP-binding protein, and a translocase (Supplementary Table S17). We searched the genes of T. denticola and P. gingivalis strains to find such candidates too; there was no candidate in both species. Then we searched for the known virulence genes in each strain genome in the two species. In T. denticola and P. gingivalis, all the known virulence genes were conserved among all the strains; however, one strain of T. forsythia species lacked one virulence gene, which was the forsythia cell-detaching factor (Supplementary figure S14). It is notable that all the virulence genes were species specific except for a novel dipeptidyl peptidase dpp, whose mechanism of virulence in T. forsythia and T. denticola remains unknown (Ohara-Nemoto ), and that there was little redundancy in function of the virulence genes among the species (Table 2). Aside from compensation of metabolic pathways, it is likely that each red complex species may provide its virulence factors to the other species in which such factors are absent: this may be an example of synergistic interactions among the three species.

Table 2

Function of virulence genes of the red complex

Species	Virulence gene	Virulece function
T. forsythia	Exo alpha sialidase, trypsin-like proteases	Degradation of host proteins
	PrtH protease, forsythia cell-detaching factor protein	Inducing detachment of adherent cells and interleukin-8 production in human fibroblasts
	Surface layer protein A	Promoting epithelial cell adherence, invasion and biofilm formation
	BNR Asp-box repeat protein	Bacterial adherence and invasion of epithelial cells
	Surface antigen BspA	Binding to fibrinogen and fibronectin and the induction of proinflammatory cytokine expression in host cells by activating Toll-like receptor 2
	Methylglyoxalsynthase	Enolization of DHAP to produce methylglyoxal
Novel virulence gene candidates of T. forsythia	L-asparaginase I	Inhibiting T-cell responses of host
	Hypothetical protein	Protein of unknown function*
	ATP-binding protein	ABC transporters involved in nutrient or metal ion uptake may result indirectly from attachment of the bacteria to host cells
	Protein translocase subunit secA 1	Secretion of factors associated with pathogenesis or stress response
T. denticola	Leucine-rich repeat protein	Binding to T. forsythia
	Dentilisin	Disrupting or modulating intercellular host signaling pathways and degradation of intercellular adhesion proteins
	Major outer sheath protein Msp	Mediating colonization of host tissue
	Lipoprotein OppA, FhbB	Epithelial cell binding and invasion, subversion of the complement cascade or tissue invasion
	Chemotaxis CheR, CheB	Damaging host tissue
P. gingivalis	Novel dipeptidylpeptidase dpp	Abscess formation and lethality
	Hemagglutinin HagA	Periodontal tissue/cell invasion
	Hemagglutinin HagB	Stimulate strong IgG and Th immune responses; induce immune protection
	Gingipains RgpA, RgpB, Kgp	Tissue destruction and alter cytokine/chemokine and Igs bioactivity (i.e., IL-12,TNFa, C3 and C5, IgG/A)
	Fimbriae FimA	Bacterial colonization, induce host IgA, IgG and Th1 immune responses
	Capsular polysaccharide	Stimulating immune protection in mouse oral challenge model (via IgG activity)

Abbreviations: ABC, ATP-binding cassette; DHAP, dihydroxyacetone phosphate; IgG, immunoglobulin G; IL, interleukin; Th, T helper; TNF, tumor necrosis factor.

*Structurally this family is part of the six hairpin glycosidase superfamily, suggesting a glycosyl hydrolase, which is one of the virulence gene (Niu ).

Bacterial interactions may enhance pathogenicity in natural host systems (Raymond ), improve sharing of a nutritious substance (Mashburn ), and promote metabolism (Shoaie ). Moreover, one study showed that the rats infected with the red complex species exhibit a significant increase in the level of serum immunoglobulin G (specific antibodies) and suggested that the interaction of the red complex species might promote periodontal tissue destruction (Kesavalu ). Nonetheless, genetic mechanisms of the actual molecular interactions among the red complex species are unknown at present. In this study, we demonstrated complementary relationships between the red complex species both in the metabolic pathways and in the virulence genes. Such relationships may have an important role not only in persistence in the deep periodontal pocket, where nutrient flexibility is limited, but also in the manifestation of virulence toward the host, as a progression mechanism of periodontitis. We showed the diversity of the symbiotic strategies among the red complex species and an interesting role of competition among the three species, although they persist together. During the initial stages of periodontitis, we believe that the interaction between the red complex species is important. Despite our findings, the precise mechanism of onset is yet to be elucidated; thus, further studies are needed for the validation of the interactions using well-established animal experimental periodontitis models and applying meta-analysis, which includes metatranscriptomic and metagenomic approaches.

In conclusion, this is the first piece of genome-wide evidence that the red complex species coexist in an oral environment based on the cooperative and competitive interactions among bacterial species. The interactions may allow bacteria to use different evolutionary strategies to survive in the variable environment. Our results will shed light on the mechanisms of interaction among not only the red complex species in periodontitis but also the bacteria involved in other polymicrobial diseases. Because this study is focused exclusively on the red complex species, behavior of the orange complex, which is a bacterial group believed to persist at high rates in the niches of the red complex (Socransky ), will be examined in future studies.