The following all encompassing search term was entered into PubMed and the NCBI Sequence Read Archive (SRA) in Oct 2023 to generate an unbiased representation of studies studying the relationship between bacteria community and apical periodontitis.

“‘apical periodontitis”[All Fields] AND “microbiome”[All Fields] OR “apical periodontitis”[All Fields] AND “microbiota”[All Fields] OR “apical periodontitis”[All Fields] AND “bacteria community”[All Fields] “

Among these, we identified 7 studies as sequencing-based research. We further examined the sampling and sequencing methods, selected only paired-end sequencing data sets, and removed samples from the maxillofacial region or samples that had not been externally sterilized during the sampling. Of these 7 studies, 4 datasets(SRA ID: SRP075560, SRP121389, ERP108053, SRP361111) from Siqueira et al. (2016); Bouillaguet et al. (2018); Qian et al. (2019); Zhang et al. (2022) encompass a collection of 224 samples that include 110 dentin and periapical microbiota from teeth with primary apical periodontitis, 95 dentin and periapical microbiota samples with secondary apical periodontitis, and 19 periapical microbiota samples from healthy molar. we filtered out 104 samples from the dentin, thereby ensuring that the characteristics we ultimately observed were associated exclusively with the disease process, rather than with the tissue type. All SAP samples were identified by imaging evidence, and none of the selected samples had been on antibiotics for two weeks.

120 final samples were derived via V3-V4 16S rRNA amplicon sequencing on Illumina-platform. Raw reads were downloaded for 120 samples from the NCBI Sequence Read Archive (SRA), then were filtered on quality in qiime2 2023.9 (Bolyen et al., 2019). Pair reads were aligned and denoised using deblur following parameters –p-trim-length 200/–p-min-reads 10. A total of 6,636 ASVs were observed across 110 samples, those ASV numbers varies from data sources ( Supplementary Figure S1 ). The feature table was annotated using naive-bayes classifer trained on eHOMD v15.23 from (Escapa et al., 2018) release in qiime2. The ASV feature table was converted to a biom file and processed with PIRCRUSt (Langille et al., 2013) for ko and ec recognition. Tree based on feature table was generated using qiime2 fasttree pipeline. We further refined the dataset by filtering out near-zero variance ASVs using the ‘nearZeroVar’ function from the ‘caret’ package and ASVs not presenting in at least 3 samples with at least a total of 10 reads with ‘Confidence.Filter’ function from ‘MicrobeR’ package, resulting in 2063 ASVs. Additionally, we filtered out datasets whose depth is less than 1000 reads. This process resulted in a final feature table containing 1740 ASVs. In order to demonstrate our data processing, we also set up 2 group with data filter1 presenting in at least 1 sample with at least a total of 10 reads, and data filter2 presenting in at least 3 samples with at least a total of 500 reads.

Raw counts were normalized using limma (Ritchie et al., 2015) package and then use sva package (Leek et al., 2012) to reduce batch effect with command (ComBat(dat=voomdata, batch=BatchVariable, mod=NULL, par.prior=TRUE), in which ‘BatchVariable’ refers to data sources. Principal variance components analysis was used to quantified these changes between raw count data and SNM-corrected data using pvca (Bushel, 2024) package in R.

For analysis on a per study basis, samples were rarefied to 5000 depth samples for generating alpha diversity metrics. The diversity and estimateR functions of Vegan (Dixon, 2003) were used to generate Shannon’s diversity index (log base e) and Chao1 estimates respectively and Picante (Kembel et al., 2010) was used to generate Faith’s phylogenetic distance. UniFrac and Jensen-Shannon divergence were calculated using the parallel-enabled distance function of Phyloseq (McMURDIE and Holmes, 2012) on subsampled proportional abundances. Bray-Curtis dissimilarity was also calculated (vegdist, Vegan) on subsampled proportional abundances. The CLR Euclidean distance was calculated by carrying out a centered log2-ratio transformation (Make.CLR, MicrobeR) with count zero multiplicative replacement [zCompositions (Palarea-Albaladejo and Martín-Fernández, 2015)] followed by calculating the Euclidean distance (dist, base). The PhILR Euclidian distance was calculated by first carrying out the phylogenetic isometric log ratio transformation [philr, PhILR (Silverman et al., 2017)] and calculating the distance matrix as before. Principal coordinates analysis was carried out using the pcoa function of APE (Paradis and Schliep, 2019). ADONIS calculations were carried out using adonis2 in Vegan on each distance/dissimilarity metric.

Our study evaluated machine learning models commonly used in genetics, including Random Forest (RF), Support Vector Machine (SVM) with radial and linear kernels, and Logit Boost (LB). RF is known for its robust performance in datasets with a high feature-to-sample ratio and is effective in handling weak predictors and complex interactions. SVM applies kernel functions for robust modeling, even with outliers. LB, part of the boosting family of algorithms, creates strong classifiers by combining multiple weak ones, known for improved accuracy and robustness. Analyses were performed in R Version 4.3.1 using the ‘caret’ package.

Models were built 5 times cross validation by random data splitting, training of the models, making predictions and recording of accuracies after each run using the caret packages for the R statistical environment (Kuhn, 2008). Cross-validation and parallel processing were enabled by the inclusion of a train control parameter. Analysis of variance, at significance alpha value of 0.05 was used to analyze the differences in mean accuracy between the models.

To enhance model performance and interpretability, our approach incorporated a rigorous search for key features within the dataset. This process involved analyzing a wide array of potential predictors, identifying those with the most significant impact on the predictive accuracy of our models. The identification of these key features is critical, especially in complex datasets, as it aids in refining the model and focusing on the most relevant variables. Since our classification target here has three categories, we need to introduce the multiROC package (Pérez-Fernández et al., 2021). The performance of our models was evaluated using two metrics: Area Under the Receiver Operating Characteristic (AUROC), which measures the model’s ability to distinguish between classes, and Area Under the Precision-Recall Curve (AUPR), important for assessing models in imbalanced datasets by evaluating precision and recall. Macro-average ROC/AUC was chosen in this study to describe each model, and it was calculated by averaging all groups results (one vs rest) and linear interpolation was used between points of ROC.

The number of predictor variables was determined by selecting the point of saturation in minimizing error rate and selecting the features based on ranked MeanDecreaseGINI, all other features were excluded from the model. Mtry and Ntree were left as default values (Mtry = sqrt(Nfeatures) and Ntree = 500).

Our exploration of clinical outcome classes revealed significant class imbalances (as shown in Supplementary Figure S1 , Supplementary Table S4 ). Such imbalances can lead to models with poor class-specific performance, as the training process tends to favor patterns associated with larger classes (Velez et al., 2007). To mitigate the effects of this imbalance on our trained model, we employed post hoc sampling approaches (Kuhn et al., 2013). Given the considerably low number of samples in the smallest class, we opted for up-sampling, utilizing methods available in the R environment. The Logit Boost learning method, along with 5-fold cross-validation, was employed to evaluate the effectiveness of our resampling approach.

63 unique studies were retrieved by our search methodology. Among these, 7 studies were identified as sequencing-based and met the eligibility criteria for inclusion in our meta-analysis. Of these 7 studies, 3 lacked public, accessible sequencing data or metadata sufficient for pairing sequencing. As an illustrative example highlighting the difficult acquiring data post-publication, we attempted to contact the corresponding author via email for sequencing data without any response. This left 4 studies for inclusion in our meta-analysis, encompassing a collection of 120 samples (Siqueira et al., 2016; Bouillaguet et al., 2018; Qian et al., 2019; Zhang et al., 2022).

The afflicted samples were limited to those presenting radiographic evidence of periapical lesions, and excluding any with crown damage, severe oral health conditions, or antibiotic treatment within one month prior to extraction. Contamination from saliva was minimized by strictly experimental condition in each study. To establish robust controls and focus our comparison on the apical, we filtered out 105 samples from the dentin, thereby ensuring that the characteristics we ultimately observed were associated exclusively with the disease process, rather than with the tooth structure itself ( Figure 1A ). From the initial 6663 ASVs identified from the sequenced data, we got 1740 ASVs which are present in at least 3 samples with at least a total of 10 reads. We examined the distribution of these remaining ASVs in different study sources ( Supplementary Figure S1 ) and further analyzed their community composition below. We also checked the current data structure and found that there are two datasets with small sample sizes (¡20), and the number between different disease stages is unbalanced. Therefore, we integrate them together for processing instead of applying independent validation on each dataset latter.

Sequencing data obtained from different instruments and studies often exhibit significant batch effects and contamination, especially in microbiota (Salter et al., 2014; Glassing et al., 2016; Poore et al., 2020). We applied PCoA to visualize the hidden batch effect in each study by comparing their community compositions. We first calculated common metrics for beta diversity for PCoA through Bray-Curtis dissimilarity, weighted/unweighted UniFrac, Jensen-Shannon divergence, PhILR Euclidean distance, and CLR Euclidean distance. Due to matrix sparsity, significant distance saturation was observed when all studies were aggregated ( Supplementary Figure S6 ), so only 3 metrics were employed: CLR-Euclidean, unweighted UniFrac, PhILR Euclidean, and associated scree plot. We then employed visualization strategies and statistically tested the effect of data source on community composition. Clear visual clustering independent of study was observed by using principal coordinates analysis of all distance metrics, and ADONIS analysis showed significant difference between data sources ( Supplementary Table S1 ).

Given the clear evidence the existed batch effort in each study, we implemented a pipeline to eliminate observed batch effort in all datasets. Subsequently, we applied supervised normalization (SNM) to reduce noise from data sources variability, while preserving and highlighting the biological variability of interest. Principal Variance Components Analysis (PVCA) demonstrated that SNM-correction mitigated batch effects in main technical variance from 0.504 to 0.013, and enhance the biological signal, specifically the ‘disease type’ from 0.001 to 0.089. Meanwhile, 89.5% of the residual variance in SNM-corrected data suggests the limitations of differential analysis in addressing this dataset. To extract relevant features, we need to employ more sophisticated models, indicating the necessity of applying machine learning approaches ( Figure 1B ; Supplementary Figure S3 ).

SNM-corrected data was applied to generate an unbiased distribution of periapical microbial community ( Figure 2A ), which is dominated by Firmicutes, Bacteroidetes, Proteobacteria, and Actinobacteria (Figure 1C). Within the Firmicutes, orders such as Lactobacillales have been implicated in the fermentation of dietary carbohydrates, a process integral to oral microbial homeostasis. Similarly, the abundance of Bacteroidetes, particularly the Bacteroidia class, aligns with recent insights into their enzymatic capabilities and their potential to degrade polysaccharides, a function essential in the oral ecological balance ( Supplementary Figure S4 ).

Then, we estimate whether the biodiversity of oral microbial environment relates to AP development. We calculated common metrics for alpha diversity, Chao1 richness, Shannon’s diversity, and Faith’s phylogenetic diversity. Taking all studies into consideration, it was observed that the alpha diversity of the periapical microbiota consistently decreased with the progression of apical periodontitis ( Figure 1D ; Supplementary Figure S11 ) Significant differences were found in the ACE and Chao1 richness indices between healthy samples and primary apical periodontitis samples (p-values < 0.001). This suggests that the occurrence of apical periodontitis is associated with a reduction in the alpha diversity of the periapical microbial community. Furthermore, significant differences (p-values < 0.01) in Shannon and Simpson richness indices were observed between PAP and SAP samples. This is likely attributed to the effects of root canal treatment and indicates that in apical periodontitis, a small number of pathogenic bacteria play a dominant role in disease progression.

We next re-tested the relationship between apical periodontitis and beta diversity, i.e. community composition. Multiple distance metrics acquired from statistical testing via ADONIS (analysis of variance using distance matrices) are used for principal coordinates analysis. All distance demonstrated a significant correlation between disease stage and community composition (p ¡0.01, ADONIS, Supplementary Table S1 ), albeit with the variance explained ranging from 0.032 to 0.121 (R square, Supplementary Table S1 ). In summary, both alpha and beta biodiversity is highly correlated to AP development, yet principle coordinate analysis cannot provide simple prediction accordingly.

Next, we attempted to apply machine learning on the corrected ASV table, in order to create a predictive model of AP stages. To provide better interpretability and reduce the complexity of model training, we removed ASVs with near-zero changes among all samples and discarded any ASVs un-annotated at the genus level, which left 1740 ASVs. We selected five machine learning (ML) methods—random forest (RF), support vector machine (SVM) with radial and linear kernels, gradient boosting machine (GBM), and logit boost (LB)—to cover a spectrum of approaches commonly employed in bioinformatics for analyzing complex datasets. These methods are renowned for their capacity to handle high-dimensional data and discern intricate variable interrelationships. By integrating these methods, our aim was to construct a robust, accurate, and interpretable predictive model for AP stages, a standard practice in microbiomics for classification tasks. Datasets were randomly divided into 7:3 for training and test sets, with training details in Methods. Initially, class imbalances in dataset (19:58:43 in health:PAP: SAP) resulted in poor performance with average accuracy of 5 cross-validations under 0.5 in all models. Thus, we upsampled cases from the minority classes with replacement until each class had approximately the same number, and retrained models.

Meanwhile, we employ other 3 feature representation methods on raw 6663 ASVs to enrich the interpretability of the data and elucidate the functioning and potential impacts of microbial communities. The PhilR (Phylogenetic Isometric Log-Ratio) approach transforms species abundance data into balanced evolutionary ratios using phylogenetic trees, providing a more stable and interpretable feature space that captures the true biological signals within the evolutionary context. KEGG Orthology (KO) features, derived from the KEGG database (Kanehisa et al., 2017), annotate the gene content of the microbiome to identify genes associated with known metabolic pathways and biological processes, revealing the functional potential of microbial communities. Enzyme Commission (EC) numbers annotate genes encoding enzymes in the microbiome genomes, reflecting the enzymatic repertoire and metabolic activities present within the community.

After applying up-sampling to all four features (ASV/PhilR/KO/EC), we compared the aforementioned performance metrics. Upon 3 repetitions, methods such as Gradient Boosting Machine (gbm), Random Forest (rf), and Logistic Regression (lb) demonstrated higher average accuracy. Furthermore, we depicted multi-class ROC and PR curves and computed the corresponding AUROC and AUPR values. Notably, gbm and lb showed higher values for these metrics ( Figure 2B ). Likely due to data sparsity and interspecies variation in ASV content, the model trained on ASVs counts had a high AUROC(AUROC = 0.897 in gbm model), indicating a good separation between classes, but shown lower Area Under the Precision-Recall Curve AUPR(AUPR < ¡0.85, Supplementary Figure S10 ), suggesting that the model’s performance on predicting the minority class was not as strong, particularly in cases where the positive class is less prevalent.

Logistic Regression showed to be the best-performing model according to 5-fold cross-validation accuracy, followed by gbm and svm with a radial basis function kernel, respectively ( Figure 2B ). We then employed a Logistic Regression classifier to define reliable biomarkers of AP stage based on a training set consist 70% of whole dataset. This model exhibits the highest performance on the philR features (multiclass AUROC = 0.967, Figure 2C ), indicating that the additional evolutionary node information confers enhanced its predictive power .( Figure 2C ) 10-fold cross validation was applied to determine the optimal number of each features included in the model required to minimize error rates. We noted that even with as few as 14 ASVs, 0.3% classification error rates could be obtained, emphasizing the high predictive power of the top features. ( Supplementary Figure S5 ) To visualize these, a phylogenetic tree of the 14 most informative ASVs was created. Among these ASVs, we scrutinized the abundance of the Lactobacillus genus, a group frequently encountered in food and cause contamination in microbiota analysis, ensuring that the contamination from food-borne bacteria was not shown in our features. To validate the result under varying data cleaning processes, we also compare the number of key ASVs from RF model training on dataset with different filtering criteria ( Supplementary Figure S12 ). 10 key ASVs were commonly recognized from 3 filtered data indicating the robustness of our data processing.

The most predictive ASV, as indicated by the mean decrease in GINI coefficient ( Figure 3A ), belonged to Enterococcus faecalis, which also showed significant differences in raw ASV counts ( Figure 3C ). It is a commonly reported multidrug-resistant pathogen that is prevalent in secondary apical periodontitis (SAP) in previous studies (Peciuliene et al., 2001; Johnson et al., 2006; Wang et al., 2012). However, its role has been controversial in differential 16SrRNA-sequencing-based analysis (Bouillaguet et al., 2018; Qian et al., 2019), and our analysis has now provided further confirmation across datasets. It is recognized that E.faecalis serves as the strongest indicator signal for distinguishing different stages of apical periodontitis, especially in SAP, where it exhibits significantly high abundance ( Figures 3B, C ). Cutibacterium acnes (C.acnes) and Delftia acidovorans(D.acidovorans) also displayed a similar characteristic of a substantial increase in abundance in SAP. ( Figure 3C ) This suggests that treatment failures in root canal therapy may be closely related to the persistent presence of those bacteria, which could be responsible for the development of SAP during the root canal filling process. In recent transcriptome-based differential analyses of primary and secondary AP, C.acnes unexpectedly contributed the largest number of differential genes at the transcriptomic level (Pinheiro et al., 2022). This specie was also highlighted as an active member of the persistent community in a RT-PCR based study (Nardello et al., 2020), and popular in root canal after chemomechanical procedures (Nardello et al., 2022). Remarkably, our batch-corrected data endorses D. acidovorans as a crucial contributor to the progression of SAP, as proposed in Anderson et al. (2012), while some studies have recognized its protective role in root caries (Abram et al., 2022).

Prevotellaceae family, with the highest proportion in ASVs features, was once again highlighted in PhILR features. ( Figure 3A ; Supplementary Figure S5 ) Previous studies have also emphasized the role of Prevotellaceae in periodontal inflammation (Könönen et al., 2022). Within this family, Prevotella nigrescens and Prevotella oris, initially isolated from periodontal pockets and considered important pathogenic bacteria in periodontal inflammation, were found to have a significant impact on the progression of apical periodontitis in our analysis. Notably, Alloprevotella tannerae from the Prevotellaceae family, previously reported to induce dentin caries, showed an increase in prevalence in primary apical periodontitis but a reduction in SAP. This suggests that these bacteria play a role in the pathogenesis of primary apical periodontitis, but they are effectively cleared by root canal treatment. Similarly, bacteria such as Burkholderia cepacia and Fusobacterium nucleatum, which can be effectively managed by root canal treatment, were enriched in PAP and have been previously detected in periapical abscess samples (Pinheiro et al., 2022). The role of Fusobacterium has been emphasized in primary infections at the transcriptomic level in a recent study, in which the Fusobacterium nucleatum carries over 10 antibiotic resistance genes homologs in multiple cases (Pinheiro et al., 2022).

Utilizing KO and EC features, we were able to perform an enrichment pathway analysis that honed in on particular metabolic pathways of interest. Our findings notably spotlight the phosphotransferase system (PTS), a pathway familiarly documented in biofilm of Streptococcus and Enterococcus species, which is intricately connected to the etiology of dental caries (Suriyanarayanan et al., 2018) ( Figure 3B ; Supplementary Table S2 ). This discovery is in harmony with our cross-validation results concerning Streptococcus and Enterococcus in the ASVs feature. Complementing this, our data derived from EC features also showed a significant enrichment in fructose and mannose metabolism ( Supplementary Table S3 ). These pathway-based analyses collectively underscore the significance of carbohydrate metabolism pathways, highlighting their pivotal role in the pathogenesis of AP and indicating potential biological mechanisms.