Microbiome analysis

The oral microbiota were examined by parallel pyrosequencing of the 16 S rRNA gene on a GS FLX. In a bar-coding approach, 12 distinct pools of samples were sequenced and analyzed - deciduous, early mixed, late mixed, and permanent dentition at healthy, treated or caries state. Of all reads, 130572 reads (78.3%) passed quality filtering and 126174 reads remained after a cut-off of 5 reads per sequence was applied. On average, 10515 reads per sample pool were analyzed. A total of 1045 unique sequences were identified. The sequences represented 8 phyla and 113 higher taxa (genus or more inclusive taxa when sequences could not be confidently classified to the genus level) (Additional file 1). Four phyla - Firmicutes, Bacteriodetes, Proteobacteria and Actinobacteria - predominated in all groups (Figure ​(Figure1),1), while three phyla - Fusobacteria, candidate division TM7 incertae sedis and Spirochaetes - were found in relatively low proportions. Representatives of the eighth phylum, Cyanobacteria, possibly reflecting plant chloroplasts, were found only sporadically, with the highest number of reads in the pool of late mixed dentition with caries (Additional file 1).

Open in a separate window

Figure 1

Relative abundance of the main bacterial phyla (99.9% of the reads) identified in saliva of children at various developmental stages of their dentition. Average and standard deviations from three pools of samples - healthy, treated and caries - per each developmental stage. * - Adult saliva data are added for comparison from Keijser et al. 2008 [1].

The deciduous dentition harboured a relatively high proportion of Proteobacteria compared to Bacteroidetes. In all other groups, the Bacteroidetes were found to be more abundant than Proteobacteria. In comparison to adult saliva [1], child saliva appears to have a higher proportion of Firmicutes and Actinobacteria and a lower proportion of Bacteroidetes, Fusobacteria, TM7 (Figure ​(Figure1)1) and Spirochaetes at all stages of the dentition.

Representatives of the Veillonellaceae family (unclassified Veillonellaceae, Veillonella and Selenomonas) and the genus Prevotella increased with increasing age (Figure ​(Figure2),2), while Carnobacteriaceae (mainly representatives of genus Granulicatella) showed the opposite trend.

Open in a separate window

Figure 2

Relative abundance of the predominant (92-95% of the reads) bacterial groups (family or genus if only one genus comprised the family, or unclassifiable bacteria within one phylum) in saliva samples of children at various developmental stages of the dentition. Average and standard deviations from three pools of samples - healthy, treated and caries - per each developmental stage. NA - not assigned.

The higher proportion of Proteobacteria in the deciduous dentition than at any other stage was mainly due to a higher proportion of Gammaproteobacteria at this stage of dentition. The order Pseudomonadales with Pseudomonadaceae (genus Pseudomonas) and Moraxellaceae families (genus Acinetobacter, Moraxella and Enhydrobacter) contributed the most to this difference. At the deciduous dentition stage, about 8% of all reads in the healthy dentition and 4% of the reads in the group with untreated caries belonged to genus Acinetobacter while sequences belonging to this genus were nearly absent from all other stages (Additional file 1, Figure ​Figure2).2). Other groups that also showed a higher abundance in the deciduous dentition included the Enterobacteriaceae family and the genus Aggregatibacter within the Pasteurellaceae family.

No significant differences were observed with respect to species richness between the different sub-groups. Rarefaction curves however (Additional file 2) did show a decreased estimate in species richness in children with caries compared to healthy or treated group.

Individual phylogenic profiles of saliva

The oral microbiota were analyzed at an individual level using a taxonomic microarray based on discriminating regions of the small subunit (16S) ribosomal RNA gene. Of the 350 microarray probes, 156 gave a positive signal (signal/background ratio at least 3) in any of the samples analyzed, with, on average, 77 probes (range 48 - 93, s.d. 10.4) per individual sample (Additional file 3). All samples were positive for probes targeting the phyla Actinobacteria (genus Actinomyces and Rothia), Bacteroidetes (genus Prevotella and order Bacteroidales), Firmicutes (genus Streptococcus, Veillonella, Lactobacillus, Granulicatella) and Fusobacteria (genus Fusobacterium). All but one sample were positive for probes targeting the phylum Proteobacteria. Betaproteobacteria was the most prevalent class (genus Neisseria - in 97% of the samples), followed by Gammaproteobacteria (genus Haemophilus) and Epsilonproteobacteria (genus Campylobacter). Probes targeting candidate division TM7 were positive in 40% of the samples; the probes targeting the domain Archaea (Additional file 4: probe o1311, designed to target Crenarchaeota, Class Thermoprotei) were positive in 44% of the samples. The occurrence of Pseudomonas in the group of children with deciduous teeth was confined to one individual, and did not reflect a common feature of the oral microbiota of this group.

To assess which probe-signals in the salivary profiles associate with any of the groups (the stage of the dentition and the oral health status), we performed the Significance Analysis of Microarrays (SAM analysis) [8]. No significant associations of microarray probe signals were found with the stage of dentition (deciduous, early or late mixed and permanent dentition). The health status (healthy, treated or carious) did show significant association with the signal of three probes targeting the Porphyromonas group (Additional file 4: probes o1503, o1504, o1506; the closest match with known species - P. catoniae, P. gingivalis, P. gulae, P. macacae), and one probe targeting Neisseria flavescens (o1486) (Figure ​(Figure3).3). The healthy group had higher signal values of these probes than the group with active caries.

Open in a separate window

Figure 3

Results of the Significance Analysis of Microarrays (SAM) on individual sample profiles obtained by phylogenetic microarray. Individual saliva samples are shown in the columns, the probes - in the rows. Sample groups by health status and by dentition are shown in different colours. High probe signal abundance is depicted red, low - dark blue. Only probes which had significant association with the oral health status are shown. No association with the stage of dentition and individual probes was found.

Dimensional reduction of the salivary profile data by principal component analysis (PCA) explained 45% of the total variance among the individual samples by the first three components (Figure 4A, B). Principal components PC2 and PC3 (14% and 11.4% of variation, respectively) partly discriminated the treated and carious samples from the samples that belonged to the healthy group. The three probes that contributed most to this difference were the same probes identified by the SAM analysis - probes (probe o1486, p = 0.04; probe, o1503 p = 0.04; probe o1506, p = 0.015; ANOVA, Games-Howell post-hoc test) targeting N. flavescens and Porphyromonas.

Open in a separate window

Figure 4

Results of the Principal Component Analysis (PCA) on individual sample profiles obtained by phylogenetic microarray. A) the plot of the PCA axis 1 (accounting for 21% of intersample variation) and the axis 2 (14% of intersample variation) B) the plot of the PCA axis 2 and the axis 3 (11.4% of intersample variation). Green dots - samples from the caries-free group, blue dots - samples from the treated group, red dots - samples from children with dental caries. Samples identified with the same number belonged to the siblings from one family. Samples without numbers belonged to unrelated individuals.

In our study group, 35 children had one or more siblings (brothers or sisters) participating in this study, representing 16 families. To quantify the similarities between salivary microbiota profiles of siblings and unrelated individuals we calculated the sample distance matrix (Additional file 5) and visualized it using hierarchical clustering (Figure ​(Figure5).5). If two sample profiles were identical, they would receive similarity value of 1. The average similarity between the respective family members was 0.87 (range 0.64 - 0.96; s.d. 0.09). The similarity between saliva profiles of unrelated children was slightly lower than the similarity between siblings (0.85; range 0.51 - 0.98; s.d. 0.08), but this difference was not statistically significant (Independent samples T-test, p = 0.275). For a few siblings, a higher similarity score was obtained, as for example members of family number 9 that showed nearly identical oral microbial communities (Figure ​(Figure44).

Open in a separate window

Figure 5

Visualization of similarities between salivary profiles by hierarchical clustering analysis (average linkage method, Pearson correlation). Samples are depicted in the columns. Siblings, who participated in the study, are marked by the same colour, each colour representing different family. Unrelated individuals have no colour code.

Sample Collection

Unstimulated saliva, as the representation of an average sample of the whole oral ecosystem [3,4], was collected at the day of the dental check-up, at home, before breakfast and toothbrushing, by drooling into a DNA-free, sterile vial for 5 minutes. The parents were instructed to store the saliva refrigerated (4-7°C) and take the vial to the dental clinic, i.e., within 6 h after collection.

DNA Extraction

A 0.1-ml quantity of saliva sample was transferred to a sterile screw-cap Eppendorf tube with 0.25 ml of lysis buffer (AGOWA mag Mini DNA Isolation Kit, AGOWA, Berlin, Germany). Then 0.3 g zirconium beads (diameter, 0.1 mm; Biospec Products, Bartlesville, OK, USA) and 0.2 ml phenol were added to each sample. The samples were homogenized with a Mini-beadbeater (Biospec Products) for 2 min. DNA was extracted with the AGOWA mag Mini DNA Isolation Kit and quantified (Nanodrop ND-1000; NanoDrop Technologies, Montchanin, DE, USA).

PCR Amplification, Sample Pooling and Pyrosequencing

PCR amplicon libraries of the small subunit ribosomal RNA gene V5-V6 hypervariable region were generated for the individual samples. PCR was performed using the forward primer 785F (GGA TTA GAT ACC CBR GTA GTC) and the reverse primer 1061R (TCA CGR CAC GAG CTG ACG AC). The primers included the 454 Life Sciences Adapter A (forward primer) and B (reverse primer) fused to the 5' end of the 16 S rRNA bacterial primer sequence and a unique trinucleotide sample identification key per each sample group. This resulted in 12 distinctly labelled pools of samples - deciduous, early mixed, late mixed, and permanent dentition at healthy, treated or caries state. The amplification mix contained 2 units of Pfu Ultra II Fusion HS DNA polymerase and 1× PfuUltra II reaction buffer (Stratagene), 200 μM dNTP PurePeak DNA polymerase Mix (Pierce Nucleic Acid Technologies, Milwaukee, WI), and 0.2 μM of each primer. After denaturation (94°C; 2 min), 30 cycles were performed that consisted of denaturation (94°C; 30 sec), annealing (50°C; 40 sec), and extension (72°C; 80 sec). DNA was isolated by means of the MinElute kit (Qiagen, Hilden, Germany). The quality and the size of the amplicons were analyzed on the Agilent 2100 Bioanalyser with the DNA 1000 Chip kit (Agilent Technologies, Santa Clara, CA, USA) and quantified using Nanodrop ND-1000 spectrophotometer. The amplicon libraries were pooled in equimolar amounts and sequenced unidirectionally in the reverse direction (B-adaptor) by means of the Genome Sequencer FLX (GS-FLX) system (Roche, Basel, Switzerland).

Processing of the Pyrosequencing Data

GS-FLX sequencing data were processed as previously described [35]. In brief, we trimmed sequences by removing primer sequences and low-quality data (sequences that did not have an exact match to the reverse primer, that had an ambiguous base call (N) in the sequence, or that were shorter than 50 nt after trimming). We then used the GAST algorithm [36] to calculate the percent difference between each unique sequence and its closest match in a database of 69816 unique eubacterial and 2779 unique archaeal V5-V6 sequences, representing 323499 SSU rRNA sequences from the SILVA database [37]. Taxa were assigned to each full-length reference sequence using several sources including Entrez Genome entries, cultured strain identities, SILVA, and the Ribosomal Database Project Classifier [38]. In cases where reads were equidistant to multiple V5-V6 reference sequences, and/or where identical V5-V6 sequences were derived from longer sequences mapping to different taxa, reads were assigned to the lowest common taxon of at least two-thirds of the sequences. Only sequences that were found at least 5 times were included in the analyses. This strict and conservative approach was chosen to preclude inclusion of sequences from potential contamination or sequencing artefacts.

Probe Design for the 16S rDNA Microarray

For the design of the taxonomic microarray literature inventory was made on the description of the normal oral microbiota [2,3,39-41]. The literature based inventory was compared to the 454 pyrosequencing data that were obtained in this study. While the majority of taxa observed by 454 sequencing were already included in the literature inventory, a number of taxonomic groups were found to be missing. These were included in the list. Based on the list, probes were designed using ARB [42]. The probes were 20-22 nucleotides in size, with a predicted melting temperature of 60°C ± 5, and a GC level between 40 and 60%. The full list of oligonucleotides is provided in Additional file 4.

Microarray Printing

Microarrays for this study were produced in house, using the ArrayIt Nanoprint60 Microarray Robot. The 5'-amino-modifier-C6-linked oligonucleotides were diluted to a final concentration of 25 μM in a 50 mM Phosphate buffer (pH 7), and printed onto CodeLink Activated Slides. Following incubation at 80% humidity under ambient conditions, slides were blocked in a buffer of 0.1 M Tris Base and 50 mM ethanolamine (pH 9) for 45 minutes at 50°C while shaking in accordance with the manufacturer's recommendations. After that, the slides were rinsed twice in MilliQ purified water (MilliQ, Millipore), washed in 4 × SSC (600 mM sodium chloride/60 mM sodium citrate), 0.1% sodium dodecyl sulfate (SDS) for 30 minutes at 50°C, and washed in pre-warmed (50°C) MilliQ water. Following that, the slides were washed twice in MilliQ water at a room temperature, and dried under a stream of nitrogen. Slides were stored under nitrogen until use.

DNA Labelling and Hybridization

For taxonomic microarray analysis, 16 S rDNA was amplified by PCR as described above, using forward primers 8F (AGA GTT TGA TCH TGG YTC AG) and 8F-bif (TGG CTC AGG ATG AAC GCT G) and reverse primer 1061R (TCA CGR CAC GAG CTG ACG AC). Following PCR amplification, DNA was amplified by random Klenow amplification with the BioPrime DNA Amplification kit (Invitrogen) according to the manufacturer's recommendations. Klenow-amplified DNA was passed through an Illustra AutoSeq G-50 column (GE Healthcare G50) for purification and concentrated in a speedvac. Following, the DNA was labelled by Terminal Transferase coupling of Cy3-dUTP (Promega). Following Illustra AutoSeq G-50 column purification, and vacuum concentration, the labelled DNA was dissolved in 40 μl Easyhyb hybridization buffer (Roche), and denatured for 2 min at 95°C. Printed slides were pre-hybridized in 0.45 μm-filtered pre-hybridization buffer [1% BSA, 5 × SSC, and 0.1% SDS] at 42°C for 45 min with rotation, then washed twice with MilliQ purified water, dried with nitrogen, and pre-warmed at 42°C. The microarrays were placed in the ProPlate multi-array system (Grace Bio-Labs). The hybridization mixture was then pipetted in the individual wells, and incubated in a hybridization chamber for four hours. Following hybridization, slides were then thoroughly washed sequentially in 1 × SSC, 0.2% SDS for 10 sec at 37°C, 0.5 × SSC for 10 sec at 37°C, and twice in 0.2 × SSC for 10 min at room temperature. Slides were dried with nitrogen and scanned using a Scanarray Express 680013 Microarray Analysis System (Perkin Elmers Life Analytical Sciences Inc.). Images were obtained and quantified with ImaGene 4.2 software (Biodiscovery).

Microarray validation

To validate microarray performance, a number of tests were performed. These included: 1) Replicated analysis of samples to test robustness of microarray data; 2) Spiking DNA of underrepresented species to a complex mixture followed by microarray analysis; 3) performing a direct comparison between qPCR and microarray data, and 4) performing a direct comparison between sequencing and microarray data.

Replicated Analysis

Technical replicates were included to verify robustness of data. In the examples described below, three replicates were compared. Two were analyzed on slides of the same printing series - sMF02-11, and one was analyzed on an older slide series (sMF02-9). The Pearson correlation (r) between the three microarray analyses was high, and varied between 0.91 and 0.93, and was independent of the slide series. For comparison, correlation between different saliva samples ranged from 0.4 to 0.7.

DNA Spiking

To test the performance of the microarray with the complex oral microbial community, we performed a spiking experiment. For this, DNA was isolated from saliva samples of a healthy adult donor as well as from pure cultures of Lactobacillus casei (ATCC 334). Following, DNA was quantified, and pure culture DNA was added to saliva DNA at levels ranging between 0.01 and 1% of the total community. Then the DNA was labelled, and used to hybridize the microarray. A clear signal was detected in all spiked samples for the L. casei spot on the microarray in a dose responsive manner. At the family level, also a significant dose responsive increase in microarray spot intensity was detected (with lower sensitivity).

Direct Comparison between qPCR and Microarray Data

A third method for validating microarray data was by cross comparing species levels established by quantitative PCR with the results obtained using the microarray. In the current study, this was performed for Porphyromonas catoniae in saliva samples obtained from the children. During the study, 74 saliva samples were analyzed. The q-PCR data of P. catoniae showed that saliva samples contained between 10⁵ and 10⁸ of P. catoniae cells/ml saliva. The fluorescence intensity for the P. catoniae probes on the microarray ranged between 0 and 110. The number of P. catoniae cells correlated significantly with the signal of both Porphyromonas probes - o1503 and o1506. The Spearman's correlations were p < 0.001; r = 0.618 and p < 0.001; r = 0.580, respectively.

Direct Comparison between Sequencing and Microarray Data

To enable a direct comparison between microarray and sequencing data, both data sets were coupled by using the probe and sequence data. For this, a string search was used to identify matching probe sequences within the 454 pyrosequencing data. For 26 probes a direct match was identified. This relatively low number of matching probe sequences is due to the fact that the pyrosequencing data span the V5-V6 region, while the majority of probes were defined in the V2, V3 and V4 region of the 16 S ribosomal gene. We then calculated the Pearson correlation between the microarray and 454 sequencing data. Good (r = 0.85) to excellent (p > 0.9) correlation was found for probes corresponding to relatively abundant taxa (>0.1% of all reads; signal/background fluorescence value >20). The correlation for less abundant species was weak. No false positive probes were found. Only two probes (o1402 and o1445) were classified as false negatives.

Targeting of specific microorganisms by quantitative PCR

Quantitative PCR was performed on the Applied Biosystems 7500 Fast Real-Time PCR System. The total microbial load was determined using the universal primer set described by Nadkarni et al [43]. The P. catoniae and P. gingivalis primer-probe sets were designed using Primer Express (Applied Biosystems). To quantify P. catoniae, the primer-probe set of P. catoniae-16S-F (CGG TTG CCA TCAG GTA ATG C), P. catoniae-16S-R (CAC CTT CCT CAC GCC TTA CG) and P. catoniae-16S-probe (TCC GTA GAG ACT GCC G), a minor-groove binding probe (MGB) labeled with 6-carboxy-fluorescein (FAM) was used. The P. gingivalis primer-probe set was composed of P. gingivalis-16S-F (GCG CTCA ACG TTC AGC C), P. gingivalis-16S-R (CAC GAA TTC CGC CTG C) and the FAM labeled minor-groove binding (MGB) probe P. gingivalis-16S-probe (CAC TGA ACT CAA GCC CGG CAG TTT CAA). Quantitative PCR was performed using the Diagenode Universal Mastermix, in accordance with the manufacturer's recommendations. Standard controls included a serial dilution series of purified genomic DNA of P. gingivalis ATCC BAA-308 and P. catoniae ATCC 51270.