The Oral Microbiome in Children with Black Stained Tooth

Black stain (BS) is a characteristic extrinsic discoloration, which occurs along the third cervical line of the buccal and/or lingual surfaces of teeth, particularly in the primary dentition of humans. BS is produced by oral bacteria and byproducts of saliva, but there is a controversy about related bacteria. The aim of this study was to identify the oral microbiome in tooth BS using pyrosequencing. It was hypothesized that the oral microbiome of BS in children might be related to black-pigment producing bacteria. Supragingival dental plaque was obtained from six children (mean 8.1 years) with BS and four children (mean 8.3 years) without BS. The bacterial metagenome was obtained by pyrosequencing. The BS group contained 348 operative taxonomic units (OTUs), whereas the control group had 293 OTUs. Microbial abundance and diversity were significantly higher in the BS group (p < 0.05). In the heatmap, the correlation between samples was the same as the BS scale. At the genus level, six genera—Abiotrophia, Eikenella, Granulicatella, Neisseria, Porphyromonas and Streptococcus—were significantly different between the two groups (p < 0.05). We suggested that compositional changes in the oral microbiome are essential, and several species in the genus Neisseria, Porphyromonas and Streptococcus may be major contributors for BS formation. Although the number of subjects was relatively limited, our study is the first species-level analysis of pyrosequencing data in BS formation.


Introduction
Black stain (BS) is a characteristic, extrinsic dental discoloration that occurs on the cervical third of the tooth and follows the contour of the gingival margin, particularly in primary dentition [1]. BS is thought to be a special form of dental plaque related to calcification because it contains an insoluble ferric salt, which is likely ferric sulfide, and has high calcium and phosphate contents [2]. The ferric sulfide may be formed through a reaction between hydrogen sulfide produced by bacterial action and iron in saliva or gingival exudates.
BS occurs at any age, but there is a higher prevalence in childhood that decreases in adolescence and adulthood [3]. Studies have shown no difference in occurrence between the sexes [1]. Researchers have been interested in BS because of its correlation with the low caries experience [2,4]. However, a recent study of more than 3000 people found no association between dental caries and BS [5].
Many studies have examined BS-related bacteria, but the causative bacteria is still not clear [1,3]. Early ultrastructural examinations of BS demonstrated that the black material is a ferric salt, probably ferric sulfide, formed by the reaction between sulfide produced by bacterial action and iron in the saliva or gingival exudate [2,3]. The most commonly cultivated bacteria from BS teeth were Actinomycetes with the detection of rarely pigmented Gram-negative rods [6]. Black-pigmented bacteria were found primarily in extrinsic tooth stains [6,7]. Among black-pigmented bacteria, Prevotella intermedia and Prevotella nigrescens are the ones most frequently found in the oral cavity and associated with BS [3]. Soukos et al. [6] reported that Porphyromonas gingivalis, P. intermedia and P. nigrescens were related to BS using DNA-DNA hybridization for 40 taxa, whereas Saba et al. [7] suggested that P. gingivalis and Prevotella melaninogenica were absent and Actinomycetes sp. and Aggregatibacter actinomycetemcomitans were present in BS by means of polymerase chain reaction (PCR) and electrophoresis gel on the agarose medium for 200 people. This conflict was partly caused by limitations in the analysis methods. As a result, more advanced analysis methods are needed.
Recent advances in sequencing technologies, such as 454 pyrosequencing, have provided a deep understanding of the metagenomic diversity of the human oral microbiome [8]. Pyrosequencing can identify bacterial sequences and their amounts to give both qualitative and quantitative information of the oral microbiome, including uncultivable microorganisms [8]. For example, high-throughput pyrosequencing has revealed that dental plaques pooled from adults contain approximately 10,000 microbial phylotypes, with the ultimate diversity of oral microbiomes being estimated to contain approximately 25,000 phylotypes [9].
The aim of this study was to identify the oral microbiome in BS using pyrosequencing and to test our hypothesis that the oral microbiome of BS in children might be related to black-pigmented bacteria.

Sampling
This study was approved by the Institutional Review Board of the School of Dentistry, Kyung Hee University (KHD-IRB-1509-5). Ten children were recruited as study subjects. Six of the children had black stained teeth while four children maintained good oral hygiene and had no black stained teeth. Oral examinations were performed on all the children, and decayed-missing-filled teeth (DMFT) indices were estimated. The BS scale was estimated using the Shourie method [10]. The criteria for classifying BS were as follows: 1; no line, 2; incomplete coalescence of pigmented spots and 3; continuous line formed by pigmented spots. After classification, subjects with parental consent were instructed not to brush their teeth for 24 h before sampling. Sampling was performed in the morning from the subjects, who were instructed to skip breakfast. Sterilized dental explorers were used to collect supragingival plaque on the buccal and lingual surfaces of all teeth. The samples were collected in 1 mL of sterile phosphate-buffered saline and stored at −80 • C after liquid nitrogen treatment.

Genomic DNA Extraction
The following molecular and sequencing procedures were performed in an outside outsourcing research laboratory. Genomic DNA (gDNA) was isolated from the dental plaque samples using a DNA extraction kit for bacteria (iNtRON Biotechnology, Seongnam, Korea) following the manufacturer's protocol. The collected samples were treated with lysozyme at 37 • C for 15 min and lysed with a buffer containing proteinase K and RNase A at 65 • C for 15 min. The lysates were mixed with a binding buffer and then passed through resin columns for purification. The harvested g DNA from the samples was quantified using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Waltham, MA, USA).

PCR Amplification
PCR was performed to amplify the V1 to V3 region of the 16S rRNA gene from the samples' gDNA. 9F (5 -AGAGTTTGATCMTGGCTCAG-3 ) and 541R (5 -ATTACCGCGGCTGCTGG-3 ) primers were used in the amplification. The amplifications were carried out under the following conditions: initial denaturation at 94 • C for 5 min, followed by 30 cycles of denaturation at 94 • C for 30 s, primer annealing at 55 • C for 30 s and extension at 72 • C for 30 s with a final extension phase at 72 • C for 5 min. The PCR products were visualized under a Gel Doc system (BioRad, Hercules, CA, USA) after electrophoresis on a 2% agarose gel. The amplified products were cleaned with the QIAquick PCR purification kit (Qiagen, Valencia, CA, USA).

DNA Pyrosequencing
Equal concentrations of the PCR products were pooled together and cleaned with Agencourt Amperes magnetic purification beads (Danvers, MA, USA). The quality and size of the pooled products were assessed on an Agilent Bioanalyzer 2100 (Palo Alto, CA, USA) using a DNA 7500 chip. The pooled products were used for emulsion PCR and then loaded onto Picotiter plates for high-throughput pyrosequencing. The sequencing was performed on the GS FLX Plus system (Roche, Branford, CT, USA) according to the manufacturer's instructions at ChunLab, Inc. (Seoul, Korea).

Pyrosequencing Data Analysis
Obtained raw data were categorized into samples according to the barcode information. The barcode, linker and primers were trimmed off from the reads and reads showing ambiguous nucleotides (2≤), a low-quality score (<25) or a short length (<300 bp) were removed. Potential chimera sequences detected by the Bellerophon method were filtered. The final reads were compared with sequence and taxonomic information acquired from the EzTaxon-e database. This database has curated 16S rRNA gene sequences of species and phylotypes of either cultured or uncultured entries in the GenBank database with taxonomy information. The read sequences were queried with BLASTN, and a maximum of five best-hit sequences were compared by pairwise alignment to determine species or phylotypes. Operational taxonomic units (OTUs) were obtained at a 97% similarity level by the CL community program (ChunLab, Inc., Seoul, Korea) after comparing with the EZ taxon database. Alpha diversity was determined using the "phyloseq" R package. A heatmap analysis was performed based on microbial abundance in the "heatplus" R package. Principal component analysis (PCA) was performed using "devtools" and "ggbiplot" R packages.

Statistical Analysis
In clinical findings, dmft indices between the BS and control groups were statistically analyzed by Mann-Whitney test (p = 0.05), and the association between BS formation and caries development were statistically analyzed by Chi-squared test (p = 0.05). After the pyrosequencing, the diversity of the oral microbiome of the subjects was statistically analyzed by the t-test (p = 0.05). In the heatmap, the Bray-Curtis method was used for calculating cluster dissimilarity in the clustering analysis (p = 0.05). The abundance of the predominant bacterial groups was analyzed by t-test (p = 0.05). Jonckheere's trend test was used for the trend analysis (p = 0.05).

Clinical Findings
Children with BS (BS group, n = 6) whose ages ranged from 4 to 11 (mean of 8.1 years) were included in this study. The control (non-BS) group (n = 4) had a mean age of 8.3 (range 6-9) years. The BS group exhibited a mean BS score of 1.5 (BS scale 1 = 3 and BS scale 2 = 3), zero DMFT index and a 3.8 average dmft index. By contrast, four control patients showed zero BS scale, zero DMFT index and a 1.75 average dmft index. The dmft indices between the BS and control groups were not significantly different (p > 0.05). No significant association between BS formation and caries development was observed (p > 0.05).

Diversity of the Oral Microbiome
A total of 10 samples were sequenced, and 60,371 reads in total were obtained. Of all the reads, 39,013 reads (64.6%) passed quality filtering (low-quality sequence reads and a minimum length of 300 bp) and remained. Both the BS and control groups contained 28,142 and 10,871 reads with an average length of 472 and 465 bp, respectively (Supplementary Table S1). The microbial abundance based on the read number was significantly higher in the BS group than the control group ( Figure 1A, p < 0.05). OTUs were identified by pairwise alignment to the EzTaxon database at a 3% distance. Totals of 15 phyla, 28 classes, 45 orders, 80 families, 145 genera and 423 OTUs (species) were identified in these samples. The BS group contained 348 OTUs in total (without duplication), whereas the control group had 293 OTUs. The alpha diversity of the microbiome was estimated by Simpson and Shannon indices, indicating that the diversity in the BS group was significantly higher than in the control group ( Figure 1B,C, p < 0.05).

Diversity of the Oral Microbiome
A total of 10 samples were sequenced, and 60,371 reads in total were obtained. Of all the reads, 39,013 reads (64.6%) passed quality filtering (low-quality sequence reads and a minimum length of 300 bp) and remained. Both the BS and control groups contained 28,142 and 10,871 reads with an average length of 472 and 465 bp, respectively (Supplementary Table S1). The microbial abundance based on the read number was significantly higher in the BS group than the control group ( Figure  1A, p < 0.05). OTUs were identified by pairwise alignment to the EzTaxon database at a 3% distance. Totals of 15 phyla, 28 classes, 45 orders, 80 families, 145 genera and 423 OTUs (species) were identified in these samples. The BS group contained 348 OTUs in total (without duplication), whereas the control group had 293 OTUs. The alpha diversity of the microbiome was estimated by Simpson and Shannon indices, indicating that the diversity in the BS group was significantly higher than in the control group ( Figure 1B,C, p < 0.05). Shannon index (C) were compared between the black stain and control groups. Significance was determined by t-test (p = 0.05).

Different Genera and Species between BS and Control Group
Heatmap analysis was employed to see overall patterns of microbial abundance at different levels of taxonomy. At all levels from species to order, the control group was isolated from the BS group (Figure 2), suggesting that the microbial composition was quite different between the two groups. The only exception was the sample "C12", which formed its own lineage within the BS group. Since the diversity pattern of C12 was much different from the other samples in the control group, we excluded the sample from further statistical analysis. To test the hypothesis that changes in microbial diversity affected BS formation, the abundance (read number) of each OTU was compared between the groups using the Mann-Whitney U test. At the genus level, six genera-Abiotrophia, Eikenella, Granulicatella, Neisseria, Porphyromonas and Streptococcus-were significantly different between the two groups (p < 0.05). Neisseria, Porphyromonas and Streptococcus contained 11, 14 and 31 species within them, respectively ( Figure 3). All six genera were more abundant in the BS group than the control group. In addition, the genus Selenomonas exhibited the opposite trend and was more abundant in the control group (p = 0.053). At the species level, 19 species were identified as significantly different between the two groups (p < 0.05, Table 1). The number of each sample represented the identified OTUs. Except for Selenomonas noxia, Corynebacterium_uc and JQ454562_s (Prevotella), all the species were higher in the BS group than the control group. Twelve out of 19 species were included in the genus showing significance. Shannon index (C) were compared between the black stain and control groups. Significance was determined by t-test (p = 0.05).

Different Genera and Species between BS and Control Group
Heatmap analysis was employed to see overall patterns of microbial abundance at different levels of taxonomy. At all levels from species to order, the control group was isolated from the BS group (Figure 2), suggesting that the microbial composition was quite different between the two groups. The only exception was the sample "C12", which formed its own lineage within the BS group. Since the diversity pattern of C12 was much different from the other samples in the control group, we excluded the sample from further statistical analysis. To test the hypothesis that changes in microbial diversity affected BS formation, the abundance (read number) of each OTU was compared between the groups using the Mann-Whitney U test. At the genus level, six genera-Abiotrophia, Eikenella, Granulicatella, Neisseria, Porphyromonas and Streptococcus-were significantly different between the two groups (p < 0.05). Neisseria, Porphyromonas and Streptococcus contained 11, 14 and 31 species within them, respectively ( Figure 3). All six genera were more abundant in the BS group than the control group. In addition, the genus Selenomonas exhibited the opposite trend and was more abundant in the control group (p = 0.053). At the species level, 19 species were identified as significantly different between the two groups (p < 0.05, Table 1). The number of each sample represented the identified OTUs. Except for Selenomonas noxia, Corynebacterium_uc and JQ454562_s (Prevotella), all the species were higher in the BS group than the control group. Twelve out of 19 species were included in the genus showing significance.

Trend Analysis among the Three BS Scale Groups
In the heatmap at the species level, clusters based on microbial abundance were correlated to BS scale severity (Figure 2). For example, the samples "B1", "B2" and "B9", which had been determined as having "BS scale 2" severity according to Shori's method, were clustered together in similarity (Figure 2). A similar correlation was observed at all taxonomical levels, suggesting that BS scale severity may be correlated with the bacterial composition ( Figure 2). Therefore, a dose-dependent relationship between the BS scale and abundance was examined using a nonparametric Jonckheere's trend test in the "clinfun" R package. A total of 47 species showed a trend among the BS scale clusters: 39 species exhibited an "increasing" trend when the BS scale increased (p < 0.05), whereas eight species were found to exhibit a "decreasing" tendency (p < 0.05). Principal component analysis (PCA) for these selected species revealed that species showing an "increasing" trend contributed to an increase in BS scale severity (Figure 4). This microbial community may play an essential role in the formation of BS. The most abundant genera were Neisseria, Porphyromonas, Streptococcus and Capnocytophaga. In contrast, the "decreasing" population, which showed abundance in the control group, consisted of Corynebacterium, Prevotella, Selenomonas, Capnocytophaga and Leptotrichia. Four genera were found in both groups: Capnocytophaga, Corynebacterium, Leptotrichia and Prevotella.

Trend Analysis among the Three BS Scale Groups
In the heatmap at the species level, clusters based on microbial abundance were correlated to BS scale severity (Figure 2). For example, the samples "B1", "B2" and "B9", which had been determined as having "BS scale 2" severity according to Shori's method, were clustered together in similarity (Figure 2). A similar correlation was observed at all taxonomical levels, suggesting that BS scale severity may be correlated with the bacterial composition ( Figure 2). Therefore, a dose-dependent relationship between the BS scale and abundance was examined using a nonparametric Jonckheere's trend test in the "clinfun" R package. A total of 47 species showed a trend among the BS scale clusters: 39 species exhibited an "increasing" trend when the BS scale increased (p < 0.05), whereas eight species were found to exhibit a "decreasing" tendency (p < 0.05). Principal component analysis (PCA) for these selected species revealed that species showing an "increasing" trend contributed to an increase in BS scale severity (Figure 4). This microbial community may play an essential role in the formation of BS. The most abundant genera were Neisseria, Porphyromonas, Streptococcus and Capnocytophaga. In contrast, the "decreasing" population, which showed abundance in the control group, consisted of Corynebacterium, Prevotella, Selenomonas, Capnocytophaga and Leptotrichia. Four genera were found in both groups: Capnocytophaga, Corynebacterium, Leptotrichia and Prevotella.

Discussion
Conventional BS therapy has been focused on the mechanical removal of black pigments themselves because the exact mechanism of BS prevalence had yet to be identified [1]. In this study, we used pyrosequencing analysis of supragingival dental plaque of children with and without BS to understand the microbiological aspects of BS, which may lead to better prevention and treatment for BS. The genera Porphyromonas, Streptococcus and Neisseria, were abundant in the microbiome of children with BS, while the genera Corynebacterium, Prevotella and Selenomonas_g1 were abundant in children without BS in this study. This study suggests that dynamics in the oral microbiome may cause the formation of BS in children, and the microbial population in the above genera would play important roles in changing microbial dynamics.
The genus Porphyromonas was significantly abundant in the BS group ( Figure 2). All the species belonging to this genus exhibited a similar pattern in abundance (Figure 3). Four Porphylromonas sp. showed a trend to increase as the BS scales increased (Figure 4). Porphyromonas is an obligate anaerobic bacterium that produces a black pigment called porphyrin [11]. This pigment can protect the bacterial cell from oxidative stress by reducing oxygen. P. gingivalis is one of the most well-known species that produces black pigments [12] and is implicated in the progression of periodontal disease [4]. However, this species was not one of the pathogens observed in this study. Porphyromonas catoniae is the most common species in this genus ( Figure 3A). It is frequently found in the oral microbiota of healthy children [9] and is known to be a pioneer bacterium that is established in the mouths of babies [13]. In addition, it is a non-pigmented species [14]. However, it is possible that this species may play an opportunistic role in the formation of BS because this species exhibits activities of diverse cellular enzymes like trypsin, which is important in bacterial colonization [13]. The other Porphyromonas spp. may contribute to the accumulation of black pigments.
The genus Streptococcus was also observed abundantly in the BS group (Figures 2 and 3B). S. cristatus, S. sanguinis, S. oralis, S. mitis and S. gordonii are known as oral biofilm bacteria that colonize the tooth surface [15,16]. Since these species are normal tooth microbiome and do not produce black pigments, they may contribute to the formation of microbial beds on teeth and interact with other species producing black pigments. Although the genus Streptococcus makes up 10% of the total microbiome, S. mutans was detected in only one subject (BS2, five reads). This rare detection of S. mutans is consistent with a recent report of BS patients where S. mutans was not found. This may be due to the antagonistic effect of the predisposed Streptococcus spp. S. sanguinis and S. gordonii that are known to inhibit the growth of S. mutans by producing hydrogen peroxide [17]. A recent study using pyrosequencing suggests that the genus Streptococcus was most abundant but showed no difference between the groups [18]. Another pyrosequencing study reported that this genus was observed frequently in subjects having both caries and pigments, but not in subjects having only pigments [19]. Interestingly, a recent microbiome study showed that the population of P. catoniae was decreased in children with caries and increased in children with healthy oral status [9]. Since the P. catoniae population was found to be increased in our BS group, it may have a negative influence on caries development ( Figure 3A).
Veillonella dispar and V. rodentium were identified as positive candidates for BS formation (Figure 4). V. dispar produces hydrogen sulfide (H 2 S) from L -cysteine [20]. The genus Actinomyces can produce H 2 S and has been known to be a strong candidate for BS formation. However, Actinomyces gerencseriae was detected in this study, and there is no report about the H 2 S production of this species.
Microbial composition is continuously affected by the environment of the oral cavity, and each population interacts with others in the same community. Therefore, a decrease in population density may result in the breakdown of one community, presenting the opportunity for a new population or community to thrive. We hypothesized that different microbial compositions would appear in different BS scales. The BS scale indicates different environmental factors affecting microbial composition. In that hypothesis, our trend analysis exhibited positive or negative candidates in BS formation. In Figure 4, some species decreased as BS scales increased, and the others were increased in proportion to the BS scale. These species may play core roles in BS formation. However, their exact roles must be identified in further studies.
Recently, two pyrosequencing studies were performed to analyze the microbiota of BS [18,19]. One was performed for 25 subjects [18], while the other study consisted of 111 subjects [19]. Li et al. [18] reported that the genera of Actinomyces, Tannerella, Treponema, Corynebacterium, Cardiobacterium and Hemophilus were more abundant in children with BS. In another study, the genera Leptotrichia and Fusobacterium contributed to BS formation, while the genera Streptococcus and Mogibacterium were important in the formation of both caries and pigments. Interestingly, three studies, including ours, did not reveal a strong consistency between them. Different ethnic backgrounds and dietary styles may provide confounding effects, suggesting that a more careful sampling design is required for this kind of study.
The limitation of this study is the small sample size. Larger sample size is needed for safe experiments and the potential confounding factor that the pyrosequencing may present. Further research will be designed by reflecting this limitation. In addition, it is meaningful to compare the study with existing healthy subjects in place of the small control group. In the study of healthy children, Proteobacteria was abundant in the group without dental caries, and Firmicutes and Bacteroides were abundant in the group with dental caries [21]. A study of healthy elderly people found that Firmicutes and Veillonella were abundant at the phyla level [22,23]. This was markedly different from the BS group in this study, in which Neisseria, Porphyromonas and Streptococcus were significant.

Conclusions
In conclusion, this study investigated the composition of the oral microbiome of supragingival plaque in relation to BS formation. Although the number of subjects was relatively limited, our study is the first species-level analysis of the pyrosequencing data in BS formation. We highlighted that compositional changes in the oral microbiome are essential, and several species in the genus Neisseria, Porphyromonas and Streptococcus may be major contributors to BS formation. The proposed species information will provide a foundation for understanding the complex ecology of BS and for identifying the causal agent for BS formation.