Oral microbiome characteristics in children with and without early childhood caries

2.1 Cohort enrollment and volunteer clinical information

Forty children aged 3–6 years (20 children with ECC and 20 healthy children) were recruited for this study after obtaining permission and informed consent from their guardians. The severity of primary caries of ECC subjects was 5–6 according to the International Dental Caries Detection and Assessment System (ICDAS-II) and the numbers of decayed, missing, and filled teeth (DMFT) of the subjects were determined by an experienced dental surgeon [13]. Children with mixed dentition were excluded from the study group. The inclusion criteria were as follows: (I) no antibiotics used within 1 month; (II) no fluoride toothpaste used and fluoride treatment within 1 month; (III) the children had no other diseases except for caries; and (IV) DMFT ≥6 in children with ECC, and DMFT = 0 in the healthy controls [14]. The clinical data of the volunteers are shown in Table 1.

2.2 Sample collection

The subjects were told not to eat or drink 2 h before sample collection. The plaques were collected from the gap of the first and second primary molars of healthy controls using dental excavators (Shanghai Kangqiao Dental Instruments Factory, China) by a dental surgeon and were assigned to the HH cohort. For each ECC patient, plaques were collected mainly from two sites. One site was the gap of the first and second primary molars (without caries), and the other was from the cavitated site of carious lesions (the CH and CC cohorts, respectively). The samples were transferred into a labeled sterile centrifuge tube (Maisinuo, China) containing 1 mL sterile phosphate-buffered saline (PBS), then transported to a laboratory refrigerator at −80°C for freezing within 2 h since sampling. The samples were kept on ice during the transfer process to avoid contamination.

2.3 DNA extraction PCR and Illumina sequencing

The genomic DNA of dental plaque samples was extracted using the Fast DNA^TM Spin Kit for Soil (16560200, MP, CA, USA) following the manufacturer’s instructions. Subsequently, the V3–V4 region (forward primer, 5’-CCTACGGGNGGCWGCAG-3’; reverse primer, 5’-GGACTACHVGGGTATCTAAT-3’) of the 16S rDNA gene fragments were amplified [15, 16]. The 50 μL reaction system contained 25 μL Premix Taq^TM plus dye (TaKaRa, China), 0.5 μL forward primer, 0.5 μL reverse primer, 2 μL DNA template, and 22 μL ddH₂O. The amplification conditions were as follows: 95 ℃ for 5 min; 40 cycles of 95 ℃ for 30 s, 52 ℃ for 30 s, 72 ℃ for 30 s, and maintained at 12 ℃ for 10 min.

The PCR products were purified after amplification using the QIAquick PCR Purification Kit (28106, Qiagen, Germantown, USA), and paired-end sequencing was performed on the Illumina Miseq platform [17].

After sequencing, Quantitative Insights into Microbial Ecology 2 (QIIME2, Boulder, CO, USA) was used to analyze the high-quality sequences.

2.4 Microbial diversity assessment

Microbial diversity, including alpha-diversity and beta-diversity, was analyzed using QIIME2. Alpha diversity was evaluated using Shannon and Simpson indices [18]. Beta-diversity was determined based on Bray-Curtis dissimilarity and unweighted-UniFrac dissimilarity [19]. Microbial community clustering was arrayed using the principal coordinate analysis (PCoA), followed by visualization with the R (version 3.6.0, Ross Ihaka & Robert Gentleman, Oakland, New Zealand) “ggplot” package.

2.5 Linear discriminant analysis effect size (LEfSe) analysis

The LEfSe analysis (https://huttenhower.sph.harvard.edu/galaxy/) was used to identify the taxa and functions of the CC, CH, and HH cohorts [20]. The Kruskal-Wallis test was applied, and the p-value was set at 0.05. For taxa differences, the LDA (linear discriminant analysis) score was set at 3.5. As a supervised learning method, LDA analysis can distinguish taxa and functions to the maximum extent [21]. The analysis was performed using the R “mass” package.

2.6 Classification model construction by randomforest

The R “random forest” package was used to classify the different cohorts based on the relative abundance at the genus level, and to identify the potential diagnostic biomarkers [22]. Subsequently, 10-fold cross-validation and 10-time repeats were used to assess the model’s performance. Lastly, the areas under curve index (AUC) and receiver operating characteristic (ROC) curves were analyzed by the “pROC” package, and the model was visualized with the “ggplot” package.

2.7 Statistical Analysis

Clinical indicators were statistically analyzed using the SPSS software (version 26.0; International Business Machines Corporation, Amonk, NY, USA). The non-parametric Mann-Whitney test was used to test the differences in age, height, weight, and other indicators of subjects in each cohort. The R platform (http://www.r-project.org/) was used for statistical analysis of the sequencing data. The difference in the diversity of oral microbiota was compared using the function of the Kruskal-Wallis test and Permutational multivariate analysis of variance (PERMANOVA) of the R “vegan” package [23]. The data are presented as mean ± SD for each cohort. Statistical significance was set at p < 0.05.

3.1 The diversity of oral microbiota in caries and healthy teeth

The V3–V4 regions of the 16S rDNA of the 60 investigated samples (20 from healthy individuals and 40 from ECC children) were sequenced. A total of 837,136 high-quality sequencing reads were obtained, ranging from 2594 to 42,458 sequences, (mean = 13,952 sequences, and median = 10,997 sequences). The Shannon and Simpson indices were used to characterize the microbial richness and uniformity of the CC, CH, and HH cohorts. The results of comparative analyses showed no significant difference among the three groups (Fig. 1A,B; Kruskal-Wallis test, p > 0.05). The beta-diversity values based on Bray-Curtis (Fig. 1C) and Unweighted-UniFrac (Fig. 1D) distance metrics were calculated to reflect the levels of similarities of composition among these cohorts. The results showed that the beta-diversity values of the CC and CH cohorts were similar. The lower value of the HH cohort (Fig. 1C,D; Kruskal-Wallis test, p < 0.001) indicated that the community similarity of the HH cohort was more stable.

To determine the similarity of the microbiota composition (Fig. 1E) and functional composition (Fig. 1F) in healthy and carious sites of each child with ECC, beta-diversity was characterized using the common PCoA based on Bray-Curtis distance metrics. The two ends of each line represent the plaque samples of different health statuses of one child with ECC. The results showed a significant difference in the taxonomical and functional composition (PERMANOVA, p = 0.003 and p = 0.001, respectively) between healthy and carious sites in patients with ECC.

Fig. 1.The diversity and oral microbiota composition in carious and healthy sites of ECC (CC and CH) and healthy controls (HH). (A) Shannon index and (B) Simpson index of microbial diversity revealed no significant difference among the CC, CH and HH cohorts. (C) Bray-Curtis analysis and (D) Unweighted-UniFrac analysis revealed that the beta-diversity of the HH cohorts was lower than that of other cohorts. The taxonomical composition (E) and functional composition (F) in each child (CC and CH cohorts) based on PCoA revealed significant differences between different sites in one child. The dotted line indicates that the oral samples were from the same individual (Bray-Curtis distance, PERMANOVA). (*p < 0.05; **p < 0.01; ***p < 0.001; NS: No Significance; Kruskal-Wallis test and PERMANOVA).

3.2 Shifts of the oral microbiota in children with ECC and healthy individuals

To understand the species distribution of different healthy states, the microbiota was analyzed at the phylum and genus levels. A total of 17 phyla, 26 classes, 51 orders, 104 families, and 273 genera were detected in the oral plaque samples. At the phylum level, Firmicutes, Proteobacteria, and Fusobacteria were the main phyla among the three cohorts (Fig. 2A). Specifically, the relative abundance of Firmicutes in the CC cohort (50.6%) was higher than that in the CH and HH cohorts (34.1% and 34.5%, respectively). In the CH cohort, the relative abundance of Proteobacteria (34.4%) was higher than that in the CC cohort (21.2%) and was similar to that of the HH cohort (31.3%). Fusobacteria was mainly concentrated in the HH cohort, with a relative abundance of 15.3%. In the CC and CH cohorts, the relative abundance of Fusobacteria was 8.2% and 9.9%, respectively.

At the genus level, the top three genera in the three cohorts were Streptococcus, Neisseria, and Leptotrichia (Fig. 2B). The relative abundance of Streptococcus was roughly equal among the three cohorts, while the relative abundance of Neisseria decreased with the degree of ECC. For each cohort, the relative abundance of Neisseria in the CC, CH and HH cohort was 6.1%, 10.6%, and 14.9%, respectively. The relative abundance of Leptotrichia was 10.5% in the HH cohort, which was higher than that in the CC (6.2%) and CH cohorts (6.5%).

Fig. 2.Shifts of the oral microbiota in children with ECC and healthy individuals. The relative abundance of oral microbiota at the phylum (A) and genus (B) levels showed the dominant microbiota among CC, CH and HH cohorts. The top 10 phyla and genera are presented.

3.3 Difference of oral microbial taxonomical and functional composition in caries and healthy children

LEfSe analysis was performed to identify the dominant microbes and their functions. We obtained 29 taxa with LDA > 3.5 (Fig. 3A) and 34 signaling pathways with LDA >2.0 (Fig. 3B). The main microbes were different in all the three cohorts. As shown in Fig. 3A, the main genera of the HH cohort were Neisseria, Leptotrichia, Porphyromonas, and Gemella. The CC cohort included Lactobacillus, Veillonella, and Prevotella 7. While the dominant bacteria in the CH cohort were Actinomyces, Bifidobacterium, and Abiotrophia.

According to the third level of KEGG pathway analysis, the functions were mainly involved in cellular processes, environmental information processing, genetic information processing, human diseases, metabolism, and organismal systems (Fig. 3B). For example, the HH cohort was involved in the biosynthesis of terpenoids and steroids, lipoic acid metabolism, and biotin metabolism, while the CC cohort was associated with histidine metabolism, biosynthesis of amino acids, and fructose and mannose metabolism. In the CH cohort, the main functions included biosynthesis of unsaturated fatty acids and the adipocytokine signaling pathway.

LDA, a supervised discrimination method, was used to assess the microbial composition (Fig. 3C) and functional characteristics (Fig. 3D). The results confirmed that the microbial composition and function varied due to differences in the oral health states.

Fig. 3.Difference of oral microbial taxonomical and functional composition in caries and healthy children. Taxonomical composition (A) and functional composition (B) for identification of the CC, CH and HH cohorts presented the oral microbiota taxa and functions with significant differences. Linear discriminant analysis was performed to maximize the separation of the CC, CH and HH cohorts based on (C) taxonomical composition and (D) functional composition.

3.4 Analysis of oral microbiota competition in different cohorts

Co-occurrence networks help us to better understand the relationship among the oral microbiotal components from an ecological perspective. Hence, the interactions of dominant genera based on the SparCC algorithm were studied. In the genus-genus co-occurrence network, 35 different abundant genera with positive or negative relationships were identified in the CC cohort (Fig. 4). Similarly, 30 and 39 different abundant genera were obtained in the CH and HH cohorts, respectively (Fig. 4). In terms of edges, two related genera were connected, and the number of edges was positively correlated with complexity. The edges of the CC, CH, and HH cohorts were 56, 59, and 63, respectively. Genera of the same color were grouped into the same module. The HH cohort had the largest number of modules and edges, which indicated that the relationship among genera in the healthy oral environment was more complicated.

Streptococcus showed a significant positive correlation with many other genera in the CC cohort, such as Bergeyella, Abiotrophia, and Rosia (Fig. 4A). While in the CH cohort, Streptococcus was positively correlated with Actinomyces and negatively correlated with the Candidy Division SR1 bacterium MGEHA (Fig. 4B). Streptococcus had the most complex network relationship in the HH cohort, but this relationship was mainly negatively correlated (e.g., Campylobacter, Prevotella, and Fusobacterium; Fig. 4C). This demonstrated that diverse microbes dominated in different groups, which caused the diverse functional modules to be enriched in the different cohorts.

Fig. 4.Correlations between enriched microbiota of the CC (A), CH (B) and HH (C) cohorts. The nodes represent different genera, and their sizes represent the relative abundance; same color represents a module. The edges among them indicate correlations; red indicates a negative correlation, and green indicates a positive correlation.

3.5 Identification of the biomarkers for ECC

Considering that each cohort had distinctive microbiota, it was possible to search for biomarkers representing different healthy states. A random forest model was constructed at the genus level. Through 10-fold cross-validation and 10 repetitions, the results showed that the error rate decreased sharply as the number of genera increased. The lowest classification error rate was observed when 10 genera were selected as biomarkers (Fig. 5A). According to the Mean Decrease Accuracy, the ranking index of the importance of characterizing genera was used to determine the top 10 genera (Fig. 5B), which were chosen as biomarkers. Their relative abundances are shown in Fig. 5C. The biomarkers included Firmicutes (Gemella, o-Lactobacillus, Veillonella, Lactobacillus and Abiotrophia), Proteobacteria (f-Neisseria andCardiobacterium), Actinobacteria (Actinobacillus), and Bacteroidetes (Prevotella 7). The accuracy of these biomarkers for model classification was assessed (Fig. 5D), which shown an AUC of 89.8% (95% CI: 73.1%–100%). These results indicated that the classification model had a high predictive ability, and the biomarkers explored could be of great value in preventing and diagnosing caries in children.

Fig. 5.Identification of biomarkers for children with caries and healthy. (A) Correlation between the number of genera and the cross-validation error. The error rate was lowest when the number of genera was 10. (B) Mean Decrease Accuracy indices of important genera presented the contribution to the accuracy of the randomForest model. The lower the value, the smaller the effect on accuracy. (C) Specific relative abundance of biomarkers of each cohort. (D) The ROC of randomForest model constructed by 10 genera. The diagonal line in the graph markers an AUC of 0.5.