This study was conducted on the campus of Bengbu Medical University, with the research protocol approved by the University’s Ethics Committee. All participants provided their written informed consent before their participation. The sample size was estimated based on data from the Fourth National Oral Health Epidemiology Survey of China, using the formula N = (Z_α/2/δ)² * P(1 − P) to ensure sufficient statistical power for estimating the prevalence of the target diseases. For the dental caries survey, the parameters were set as α = 0.05, Z_α/2 = 1.96, δ = 0.05, and p = 44.42% (derived from the prevalence among 15-year-olds) (Cheng et al., 2020), resulting in a minimum required sample size of 378 participants. For the dental fluorosis survey, P was set to 13.4% (Zhou et al., 2018), while other parameters remained unchanged, yielding a minimum sample size of 178 participants. A total of 957 college students were recruited. The study population is the same as that described in our previously published article (Liu et al., 2023), though the current research addresses a distinct scientific objective. A final analytical sample of 299 participants was selected for microbiome analysis. This cohort comprised 100 healthy controls (HC), 99 with dental fluorosis but no caries (DF), and 100 with dental caries (DC). Based on the frequency of sweets consumption, each group was stratified into high-frequency (>1 time/day) and low-frequency (≤1 time/day) subgroups, yielding the following groups: HC-L, HC-H, DF-L, DF-H, DC-L, and DC-H. All participants met the following inclusion criteria: no history of recurrent oral ulcers, pulp or periodontal diseases, or other oral pathologies; absence of systemic or genetic disorders; and no use of antibiotics or exposure to tobacco products within the past 6 months.

The oral health examination was conducted in accordance with the World Health Organization (WHO) Oral Health Surveys Methods (Ben Yahya, 2024). Standardized instruments, including plane mirrors, tweezers and CPI probes, were applied to assess dental caries (documented by ICDAS index) and dental fluorosis (classified by Dean’s index). To ensure data robustness, all examinations were carried out by dentists from the Department of Stomatology, First Affiliated Hospital of Bengbu Medical University. Prior to the survey, all examiners received uniform training. Calibration for diagnosing caries and fluorosis showed Kappa values exceeding 0.8, indicating strong inter-examiner reliability. A questionnaire was then administered to the participants to gather information on their parental education level, use of fluoride toothpaste, toothbrushing and flossing frequency, and sweet food intake frequency.

Saliva samples were collected from all participants before the clinical assessment. Participants were instructed to refrain from eating or drinking for 2 h prior to sampling and to rinse their mouths with water before collection. Unstimulated whole saliva, free of sputum, was collected by placing a 5 mL centrifuge tube under the lower lip until at least 3 mL of saliva was obtained. One milliliter of each sample was reserved for microbiome analysis. All samples were stored at −80 °C within 4 h of collection for subsequent processing.

Genomic DNA was extracted using the CTAB or SDS method. After assessing purity and concentration by agarose gel electrophoresis, appropriate aliquots were transferred into centrifuge tubes and diluted with sterile water to a final concentration of 1 ng/μL for use as PCR templates. The 16S rRNA gene harbors nine hypervariable regions (V1–V9). As reported by Digvijay Verma and colleagues, in NGS-based oral microbiome studies, the V1–V2 and V3–V4 regions are commonly targeted for amplification owing to their pronounced sequence variability and rich taxonomic information. These regions provide robust resolution for distinguishing bacterial groups and have been shown to be highly representative and reliable in revealing the phylum-level architecture of the oral microbiome (Verma et al., 2018). Guided by this evidence, we selected the V3–V4 region as the target for amplification in this study. The V3–V4 region was amplified by PCR using individual-specific barcoded primers, 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′), with primer design and selection guided by reference (He et al., 2022).

Libraries were constructed with the TruSeq® DNA PCR-Free Sample Preparation Kit (Illumina, USA) and paired-end sequencing (2 × 250 bp) was performed on the Illumina NovaSeq platform. Initial processing removed barcodes and primers with Perl v5.0, after which paired-end reads were merged into complete sequences using FLASH v1.2.7¹ (Magoč and Salzberg, 2011). Quality control was conducted with fastp,² producing high-quality clean reads, which were further processed with VSEARCH³ (Yu et al., 2022) to eliminate chimeras, resulting in valid reads for subsequent analyses. Using QIIME2, denoising was performed via the DADA2 module, with sequences of abundance <5 excluded, thereby producing ASVs and their corresponding feature tables. Taxonomic classification of ASVs was achieved via the classify-sklearn module in QIIME2 by aligning sequences against reference databases (Davis et al., 2019). Taxonomic analysis was conducted mainly with the eHOMD database (Escapa et al., 2018) and the NCBI Taxonomy database.⁴

Yielding an average of 84,624 reads per sample across 299 samples. After quality control, an average of 68,515 high-quality reads per sample were retained. We further evaluated the impact of rarefaction thresholds on the representation of community structure and confirmed that a depth of 36,365 reads per sample was sufficient to capture the majority of microbial diversity in most samples. All downstream analyses were conducted using this rarefied dataset to ensure robustness and reliability of the results.

During the experimental process, extracted deionized water was used as a negative control to identify potential contaminating sequences derived from reagents or the environment. Following sequencing, any ASVs detected in the negative controls were removed from the actual samples to eliminate low-frequency signals that may represent technical noise or incidental contamination.

Analyses of microbial communities were conducted using the QIIME2 platform. α-diversity was estimated using the Chao1 index, Observed ASVs and Shannon index (Li et al., 2013), and statistical differences among groups were tested by Wilcoxon test. Microbial community composition (β-diversity) was computed from Bray–Curtis distance matrices, with principal coordinate analysis (PCoA) performed and visualized in R. β-diversity was compared using Permutational Multivariate Analysis of Variance (PERMANOVA), implemented with the vegan package (version 2.5–7) in R. Differential taxa were identified using a combination of LEfSe analysis (LDA score >3.0) and statistical testing (Wilcoxon test) across phylum, genus, and species levels, with significance set at false discovery rate (FDR) < 0.05 and fold-change (FC) > 2 or < 0.05. Venn diagram was used to compare the differential microorganisms at the genus and species level among different groups of samples. Receiver operating characteristic (ROC) curves were employed to assess the efficacy of microbial markers in distinguishing among different groups.

An LDA score threshold of 3.0 (corresponding to a 1,000-fold effect size) was applied, consistent with common practice in microbiome research to balance sensitivity and specificity when identifying biomarkers with large effect sizes and clear biological relevance. To evaluate LDA assumptions, the analysis was performed using normalized relative abundance data, and the intrinsic log-transformation in LEfSe helped approximate normality. Although microbial data rarely satisfy all theoretical assumptions of LDA, we cross-validated the LEfSe results using non-parametric tests (Wilcoxon test with FDR correction) and observed strong concordance between the methods. This agreement reinforces the robustness of the identified biomarkers and indicates that our findings are not dependent on the assumptions of any single method (see Figure 1).

Based on the microbiome analysis results, ASV-level analysis was performed on the α-diversity indices (such as Shannon index, Chao1 index, observed_ASVs) for high and low sweet consumption frequency groups in the healthy population (HC) (Figure 2A). The results revealed that the p > 0.05, suggesting no significant differences in α-diversity metrics such as richness and evenness between the high and low sweet consumption frequency groups in the healthy population. Based on Bray–Curtis distance, PCoA was used to reduce the complex microbiome structure to a two-dimensional scatter plot. The PERMANOVA method was used to calculate the statistical indicators R² and p-value to verify whether the observed separation trend was statistically significant. It was observed that the HC group did not exhibit a clear separation trend on PC1 (explaining 19.76% of the variance) and PC2 (explaining 7.67% of the variance) (R² = 0.013, p = 0.14), suggesting that sweet consumption frequency does not significantly change the overall composition of the oral microbiome (Figure 2B). We further selected high-abundance bacteria specific to the healthy population and compared them at the phylum, genus, species levels (Top 5—Top 20 abundance) (Figure 2C). The results showed that at the phylum level, there were no significant differences between the two groups in the top five dominant phyla. At the genus and species levels, the proportion of Leptotrichia (a genus of Fusobacteria) and Lautropia mirabilis was higher in the low sweet consumption group (HC-L) than in the high sweet consumption group (HC-H), suggesting they may be the dominant bacteria in the low sweet group. The LEfSe algorithm was used to further confirm the differential microbiomes between the two groups, and the results showed that Desulfobulbus sp._HMT_041 (Desulfobulbus HMT_041) had a higher LDA score in the low sweet consumption group (HC-L), suggesting its significant representativeness in this group.

Next, we compared the differences in salivary microbiome composition between the high sweet consumption frequency (DC-H) and low sweet consumption frequency (DC-L) groups in the caries group (DC). The results indicated that no significant differences were observed between DC-L and DC-H in both α-diversity (ASV level) and β-diversity analysis (R² = 0.009, p = 0.55), implying that variations in sweet consumption frequency do not significantly alter the internal diversity or overall structure of the oral microbiome in the caries group (Figures 3A,B). Further comparison of the high-abundance microbiota (Top 5—Top 20, covering phylum, genus, species levels) between the two groups (Figure 3C) showed that there were no significant differences in high-abundance bacteria at the phylum level. At the genus and species levels, Porphyromonas was significantly enriched in the DC-L group, while Streptococcus parasanguinis_clade_411 was notably enriched in the DC-H group. These results indicate that in the caries group, although differences in sweet consumption frequency did not change the overall microbiome structure, the differences at the specific genus/species levels might reflect the dominant microbial traits under varying sweet exposure conditions. To further identify the specific differential microbiomes based on varying sweet consumption frequencies in the caries group, we applied the LEfSe method (LDA score >3.0, FDR < 0.05) to compare the microbiome composition between the high sweet consumption frequency group (DC-H) and the low sweet consumption frequency group (DC-L). The results identified 7 significant differential taxonomic units: In the DC-H group, the species enriched included Prevotella salivae (s_salivae), Veillonella atypica (s_atypica), Streptococcus parasanguinis clade 411 (s_parasanguinis_clade_411), Rothia dentocariosa (s_dentocariosa), and Veillonella dispar (s_dispar). These bacteria are closely related to sugar metabolism, acid production, and the progression of caries. The species enriched in the DC-L group included Haemophilus sp. HMT 908 (s_sp_HMT_908) and Porphyromonas gingivalis (g_Porphyromonas; s_gingivalis). They are typically related to the maintenance of oral health or periodontal-associated microbiota. In summary, the results indicate that in the caries group, high sweet consumption frequency significantly changes the salivary microbiome structure, promoting the enrichment of acid-producing and acid-tolerant bacteria; while low sweet consumption frequency helps maintain microbiota traits associated with periodontal health.

In the fluorosis group without caries (DF group), we compared the salivary microbiome composition between the high sweet consumption frequency (DF-H) and low sweet consumption frequency (DF-L) groups. α-diversity analysis showed that at the ASV level, there were no significant differences in Shannon, Simpson, and Chao1 indices between DF-L and DF-H (Figure 4A), suggesting that differences in sweet consumption frequency did not alter the richness or evenness of the community. β-diversity analysis based on Bray–Curtis distance showed separation of DF-L and DF-H samples in coordinate space (Figure 4B), and PERMANOVA testing suggested that the difference was statistically significant (R² = 0.02, p = 0.02), indicating detectable differences in the overall community structure between the two groups. Community abundance comparisons showed (Figure 4C): At the phylum level, DF-L was enriched in Proteobacteria and Bacteroidetes, while DF-H was enriched in Firmicutes and Actinobacteria. At the genus level, Streptococcus and Rothia were significantly more abundant in DF-H, while Pseudomonas was enriched in DF-L. At the species level, Rothia mucilaginosa and Streptococcus salivarius were enriched in DF-H, while Pseudomonas fluorescens was enriched in DF-L. These results suggest that the DF group exhibits notable differences in dominant microbiota profiles under varying sweet consumption frequencies. LEfSe analysis (LDA score > 3.0, FDR < 0.05) further confirmed this trend (Figure 4D): The DF-H group was enriched in Streptococcaceae (S. salivarius, S. parasanguinis clade 411) and Rothia (R. mucilaginosa), while the DF-L group was enriched in Pseudomonas (P. fluorescens). In conclusion, in the fluorosis group without caries, high sweet consumption frequency did not lower overall diversity but significantly altered the community structure, with the enrichment of potential cariogenic or acid-producing bacteria like Streptococcus and Rothia. Low sweet consumption frequency was dominated by microbiota like Pseudomonas fluorescens, which may play a role in maintaining oral homeostasis. Therefore, even without caries, individuals with fluorosis may be at risk of their oral microbiome shifting toward a cariogenic direction under a high-sugar diet.

Analysis of oral microbiota diversity and composition revealed no significant differences in Shannon, Observed, or Chao1 indices between the low-score (HC-L, DC-L, DF-L) or high-score (HC-H, DC-H, DF-H) groups (p > 0.05, Wilcoxon test). Principal coordinates analysis based on Bray–Curtis distances showed no clear separation in microbial community structure among groups. Furthermore, the relative abundances of major genera and species at different taxonomic levels did not differ significantly between groups (Supplementary Figure S1). To explore the oral microbiome response patterns under the strong influence of high-frequency sweet food, we integrated the differentially expressed microorganisms at the genus and species levels (selection criteria: FC > 2 or < 0.5, FDR < 0.05) into a Venn diagram for comparison (Figures 5A,B). The analysis identified a total of seven intersecting sets. Notably, at both the genus and species levels, a unique core bacterium, Ralstonia, was observed to show differential abundance across all three group comparisons, with abundance ordered as DC-H (caries) > HC-H (healthy) > DF-H (fluorosis). This result indicates that Ralstonia may have a close association with the occurrence of caries in a high-sugar environment. At the species level, we further identified 8 differential bacteria between DC-H and DF-H (e.g., s_Parvimonas micra, s_Streptococcus mutans, s_Prevotella melaninogenica) and 7 DC-specific bacteria (e.g., s_Lachnoanaerobaculum orale, s_Veillonella dispar, s_Prevotella oris, s_Prevotella denticola). These results indicate that in a high-sugar environment, caries-associated dominant microbiota are enriched, forming a typical caries-dominated core microbiota.

In contrast, under low-frequency sweet food conditions, we also combined the differential microbiota at the genus and species levels into a Venn diagram (Figures 6A,B). The analysis revealed that the three groups showed significant differences in their microbiome characteristics. First, the core bacteria showing differences in all three group comparisons included Stenotrophomonas and Fastidiosipila. The abundance of Stenotrophomonas was in the order of HC-L > DF-L > DC-L, indicating it may be a dominant bacterium in health, helping to maintain oral ecological stability and suppress the expansion of cariogenic microbiota in a low-sugar environment. Its reduction could be associated with microbiota imbalance in the caries group. The abundance of Fastidiosipila was in the order of DF-L ≫ HC-L > DC-L, showing significant enrichment in the fluorosis group. This result suggests that it may represent a fluoride-adapted microbiota, with its dominance reflecting the selective pressure of fluoride on the oral microbiome and potentially granting the fluorosis group a distinct ecological niche. Furthermore, pairwise comparisons further revealed community differences: 17 differential bacteria between DC-L and DF-L, suggesting that even in low-sugar conditions, distinct microbial differentiation exists between caries and fluorosis. These microbial groups may be key microorganisms responsible for the phenotypic differences between the two diseases, reflecting the fundamental differences in ecological adaptation between caries and fluoride exposure.13 differential bacteria between HC-L and DF-L, indicating that the fluorosis group is different not only from the caries group but also from healthy individuals. This suggests that the microbiota in the fluorosis group is not simply “in between health and caries,” but rather forms an independent microbial ecological pattern.

To evaluate the biological relevance of between-group differences and assess the reliability of our statistical conclusions, we quantified the magnitude of these differences by calculating Hedges’ g effect size along with its 95% confidence interval. The results indicated only minimal effects, with all g values below 0.2 and confidence intervals encompassing zero. Consistently, the β-diversity analyses for HC-L vs. HC-H and DC-L vs. DC-H groups revealed exceptionally low R² values (< 2%). Those may indicate that the effect of this factor on α and β-diversity is weak, and its potential ecological role needs to be further evaluated by combining functional data, network analysis, or multi-model inference. This further supports the statistical conclusion that no substantial differences exist between these group pairs.

Although a statistically significant difference in β-diversity was observed between DF-L and DF-H groups based on sweet food consumption frequency (R² = 0.02, p = 0.02), several considerations suggest limited biological importance. As emphasized by Nakagawa and Cuthill (2007), biological interpretation should prioritize practical effect magnitude over statistical significance alone, since minimal effects may achieve significance in large samples yet lack ecological relevance. Additionally, methodological work by Caporaso and colleagues highlights that microbial community data typically exhibit high baseline variation (high within-group heterogeneity) (He et al., 2015). Reliable detection (e.g., with >80% statistical power) of effects explaining less than 2% of variance generally requires sample sizes exceeding 200, often reaching 500 or more—far beyond the sample size available for the DF-L and DF-H groups in this study. Therefore, despite the statistical significance, we consider the biological relevance of this difference to be limited, as the grouping factor explains only approximately 2% of community variation and no clear separation is visible in PCoA ordination. These findings suggest that microbial community variation is primarily driven by other, unmeasured factors.