Periodontitis in elderly patients with type 2 diabetes mellitus: impact on gut microbiota and systemic inflammation

Gut dysbiosis in patients with T2DM and periodontitis

We analyzed fecal samples from 78 participants aged ≥ 65 years and evaluated the influence of periodontitis on gut microbiota of elderly patients with T2DM. The clinical characteristics and pathological indexes of participants are summarized in Table 1. There were no significant differences among groups in terms of age, sex, or BMI. All 78 DNA samples were analyzed for 16S rRNA gene amplification and products were sequenced. OTU accumulation boxplots were plotted to confirm that the sample size was adequate for comparative analysis (Figure 1A). The α-diversity indexes (Shannon and Simpson) did not significantly differ among the three groups (Supplementary Table 1). The environmental factors collected for analysis of β-diversity included basic clinical characteristics in this study. We focused on the major digestive tract complications that might exhibit associations with gut microbiota, including dry mouth and gastrointestinal (GI) symptoms (e.g., abdominal pain, abdominal bloating or distension, nausea/vomiting, gastroesophageal reflux, dysphagia, diarrhea, and constipation). Metadata were collected regarding 15 potential confounders, such as age, sex, BMI, duration of T2DM, hypertension, and GI symptoms, as well as use of medications (e.g., acarbose, metformin, aspirin, statins, calcium channel blockers, angiotensin receptor blockers, and pioglitazone). Eleven of the 17 variables exhibited collinearity (variance inflation factor > 2) and were therefore excluded (Supplementary Table 2). The number of remaining natural teeth, CPI score, and five other variables were used as environmental variables for PERMANOVA assessment. After adjustment for confounders, periodontal condition explained 6.6% of the variation in Bray–Curtis dissimilarity, 5.9% of the variation in unweighted UniFrac distance, and 2.3% of the variation in weighted UniFrac distance (Figure 1B, P < 0.05, PERMANOVA, Supplementary Table 3).

The dominant phyla among participants were Firmicutes, Bacteroidetes, and Proteobacteria (Figure 1C). BugBase algorithm-based prediction suggested that the T2DM_P group exhibited preferential enrichment of gram-positive taxa, compared with the healthy volunteers (Figure 1D). Predicted microbiota functions were determined by KEGG pathway analysis (Figure 1E, Supplementary Table 4). Pathways associated with cytotoxicity that were predicted to be upregulated in the T2DM_P group included the HlyD family secretion protein, UDP-N-acetylglucosamine acyltransferase, and chorismate mutase (Kruskal-Wallis test P < 0.05). Pathways predicted to be downregulated in the T2DM_P group included carbamate kinase, threonine dehydratase, NADH-quinone oxidoreductase subunit E, and acetolactate synthase I/II/III large subunit (Kruskal-Wallis test P < 0.05).

Disturbance in gut flora interaction driven by periodontitis

To infer potential interactions among gut bacterial community members, we constructed co-occurrence networks of bacterial genera, based on correlations of their relative abundances (Figure 2). Overall, the co-occurrence network of the T2DM_P group exhibited greater complexity with distinct components and topographies, compared with other groups (Figure 2A). The co-occurrence network of the T2DM_P group was characterized by two unique modules separate from the core module (one comprised Prevotella, Bacteroides, and Alloprevotella, and the other comprised Streptococcus and Veillonella). Specifically, the co-occurrence networks of the T2DM_P and T2DM_NP groups contained 72 and 82 edges, respectively; the co-occurrence network of the control group contained 104 edges (Figure 2B). Despite the presence of some overlapping edges, many edges were specific to each group. Notably, correlations between Romboutsia and Turicibacter, and between Streptococcus and Veillonella, were present in all three groups. Connections containing Butyricimonas, Odoribacter, Oscillibacter, and Eubacterium showed high frequencies of co-occurrence in samples from patients with T2DM, suggesting they were universally associated with T2DM. In addition, the T2DM_P group had more co-exclusion relationships in the bacterial community; in particular, antagonistic relationships between Alloprevotella and Methanobrevibacter, as well as between Bilophila and Actinoplanes, were only observed in the T2DM_P group. We then computed different bacterial network characteristics (i.e., degree, closeness, and betweenness) of shared nodes among bacterial networks to approximate the structural importance of individual bacterial taxa within each respective network (Figure 2C). In terms of the closeness centrality of shared nodes among bacterial networks, several genera appeared to have important effects in the T2DM_P network, including Alistipes, Bilophila, Butyricimonas, Faecalibacterium,Parabacteroides, Streptococcus, Veillonella, Odoribacter, and Pseudomonas.

Identification of key gut microbiota taxa associated with periodontitis by random forest and LEfSe

To identify key gut microbiota markers that distinguished participants with different periodontal conditions, we performed fivefold cross-validation analysis using a random forest model (Figure 3A). We created a ROC curve to measure the accuracies at which the relative abundances of key gut microbiota related to periodontitis were able to classify two groups of samples (Figure 3B, 3C). The top 34 OTU markers identified by random forest analysis were able to distinguish the T2DM_P group from healthy volunteers with an area under the ROC curve of 0.82 (95% confidence interval: 74.0%–96.0%); they could distinguish the T2DM_P group from the T2DM_NP group with an area under the ROC of 0.64 (95% confidence interval: 52.2%–87.7%). In addition, we performed LEfSe analysis to confirm the composition of significant differences between groups at various taxonomic levels, especially the Faecalibacterium genus, for which significance was confirmed by the Wilcoxon rank-sum test (Figure 3D, 3E). However, because of high inter-individual variabilities, only seven taxa exhibited significant differences in relative abundance between groups after FDR adjustment (Kruskal-Wallis test, adjusted for age, sex, and use of acarbose; FDR-corrected P < 0.05, Supplementary Table 5). Overall, we found that the abundance of the genus Prevotella was significantly elevated in the T2DM_P group, while the abundance of the genus Faecalibacterium was significantly depleted. At the species level, significantly altered abundances of Prevotella copri and Faecalibacterium prausnitzii (the only species thus far defined in the Faecalibacterium genus) were observed in the T2DM_P group.

Clinical correlates of key gut microbiota markers associated with periodontitis

Clinical factors associated with T2DM might shape the gut microbiota through long-term effects on the ecological balance between the host and the normal microbiota, including GI symptoms and a variety of medications that are used to treat T2DM. We performed forward selection analysis and redundancy analysis to evaluate the relationships between these features and the relative abundances of gut microbiota at the genus level (Figure 4A). The results indicated that the number of teeth, CPI, age, sex, and acarbose use exhibited greater power as driving factors; metformin use, statin use, pioglitazone use, and abdominal discomfort exhibited weaker impact. Only the number of teeth and CPI had significant effects (P = 0.003 and P = 0.02, respectively, Supplementary Table 6). Subsequently, we explored the correlation between the key gut microbiota markers identified above and clinical characteristics, including various disease indexes and use of medications (Figure 4B, Supplementary Tables 7, 8). The duration of diabetes was positively correlated with the relative abundances of Sutterella and Dialister and negatively correlated with the relative abundance of Blautia. The CPI score was positively correlated with the relative abundances of Haemophilus, Veillonella,Streptococcus, and Aggregatibacter; the number of teeth was positively correlated with the relative abundances of Faecalibacterium, Oxalobacter, and Eisenbergiella and negatively correlated with the relative abundances of Prevotella, Aggregatibacter, and Streptococcus. Dry mouth, another oral complication of T2DM linked to periodontitis, was positively correlated with the relative abundance of Dialister and negatively correlated with the relative abundance of Coprobacillus. Furthermore, we found that GI symptoms was positively correlated with the relative abundances of Flavonifractor, Parabacteroides,Sutterella, and Bacteroides. Notably, the relative abundances of Flavonifractor and Bacteroides were also positively correlated with acarbose use. In addition, the relative abundance of Eisenbergiella was significantly negatively correlated with hypertension and calcium channel blocker use.

We also investigated the influences of periodontitis on metabolic and inflammatory profiles. The abundances of periodontitis-associated gut microbiota taxa (e.g., Prevotella, Faecalibacterium, Haemophilus, Veillonella,Streptococcus, Aggregatibacter, Oxalobacter, and Eisenbergiella) also exhibited significant correlations with the blood levels of pro-inflammatory cytokines (e.g., PGE₂, IFN-γ, IL-17, IL-6, and TNF-α) and plasma metabolic parameters (e.g., HbA1c, fasting blood glucose [FBG], total cholesterol [TCHO], high-density lipoprotein cholesterol [HDL-c], and low-density lipoprotein cholesterol [LDL-c]). Specifically, the T2DM_P group exhibited a significantly elevated Hemoglobin A1c (HbA1c) level, compared with the T2DM_NP group (P = 0.009, Figure 5A); there were no significant differences in other serum metabolic parameters between the T2DM_P and T2DM_NP groups (Figure 5B). Adverse systemic inflammation was observed in the T2DM_P group, compared with the T2DM_NP group. Several inflammatory factors were significantly elevated in the T2DM_P group, including serum levels of tumor necrosis factor α (TNF-α) and Prostaglandin E2 (PGE₂) (P = 0.030 and P = 0.036, respectively; Figure 5C and 5D). In addition, bone loss was more serious in the T2DM_P group, as indicated by higher levels of bone alkaline phosphatase and osteocalcin (Figure 5E). Circos plot analysis revealed correlations between clinical characteristics and peripheral risk factors. Correlations were observed between CPI and peripheral concentrations of both PGE₂ and HbA1c; a significant correlation was also observed between IL-22 and duration of diabetes. Figure 5G shows a heatmap of the correlation coefficients between circulating metabolic markers and key gut microbiota taxa, as well as between circulating inflammatory markers and key gut microbiota taxa (Supplementary Table 9, 10). The relative abundance of Faecalibacterium was significantly negatively correlated with serum levels of PGE₂ and interferon-γ (IFN-γ), whereas the relative abundance of Prevotella was significantly positively correlated with plasma low-density lipoprotein level. Serum levels of PGE₂ and IL-6 were negatively correlated with relative abundance of Eisenbergiella and positively correlated with the relative abundance of Aggregatibacter. Additionally, a positive association was observed between the relative abundance of Veillonella and IL-17, while a negative association was observed between the relative abundance of Dorea and IFN-γ; these findings suggested probable contributions to systemic inflammatory responses in periodontitis.

Study design and participants

Ethics approval for the study was obtained from the ethics board of the First Affiliated Hospital, School of Medicine, Zhejiang University (IRB no. 2014-114). Written informed consent and questionnaire responses were obtained from all participants in this study. Patients with T2DM, with or without periodontitis, were enrolled in the study by their physicians. Patients provided fecal and blood samples if they fulfilled the following criteria: (1) age ≥ 65 years with body mass index (BMI) of 17–38; (2) no antibiotic use in the 12 weeks prior to enrollment. Patients were included in the T2DM_P group if they exhibited periodontitis (e.g., chronic periodontitis, aggressive periodontitis, gingivitis, periodontal abscess, or combined periodontic-endodontic lesions). Patients were included in the T2DM_NP group if they did not exhibit periodontitis. The exclusion criteria for this study were: (1) presence of acute illness (e.g., upper tract respiratory infections, pneumonitis, or enterocolitis); (2) presence of any other organ-specific diseases (e.g., intestinal diseases, pancreatic diseases, tumor recurrence, heart failure, renal failure, stroke, and/or peripheral artery disease); (3) presence of severe T2DM complications (e.g., stroke, retinopathy, and nephropathy); and (4) consumption of alcohol, tobacco, Chinese herbal medicine, and/or recreational drugs. We reviewed the medical records of 49 patients aged ≥ 65 years with T2DM (28 with periodontitis and 21 without periodontitis) and the records of 29 age- and sex-matched healthy volunteers from the same community. All healthy volunteers were within the normal ranges upon physical examination; had no history of smoking or alcohol abuse; had no airway infection, runny nose, sputum production, or other systemic disease; and had not received antibiotics, probiotics, or prebiotics in the month before enrollment. At each clinic visit, participants described their use of medications, including acarbose, metformin, aspirin, statins, calcium channel blockers, angiotensin receptor blockers, and pioglitazone.

Periodontal evaluation was performed promptly by trained dentists and the severity of periodontitis was assessed by simple World Health Organization Community Periodontal Index (CPI) codes [30]. A detailed medical and dental history was taken for each participant, followed by an oral examination that included a periodontal examination. The CPI was scored on a scale of 0 to 4, as follows: 0 = healthy periodontal tissue (no bleeding, calculus, or pocket depth ≥ 3.5 mm); 1 = bleeding in periodontal tissue (bleeding on probing, but no calculus or pocket depth ≥ 3.5 mm); 2 = supragingival or subgingival calculus; 3 = 3.5–5.5-mm-deep periodontal pockets; and 4 = ≥ 5.5-mm-deep periodontal pockets.

The Diabetes Bowel Symptom Questionnaire was used to assess the severity of GI symptoms (Supplementary Table 11). Details of the questionnaire have been described previously [31]. The questionnaire consisted of parts A and B, which contained detailed information regarding seven symptoms/symptom complexes occurring within the prior 3 months: abdominal pain; abdominal bloating or distension; nausea/vomiting; gastroesophageal reflux; dysphagia; diarrhea; and constipation. Thirty-two questions were designed to examine the presence and severity of abdominal pain/discomfort (Part A) and upper GI symptoms (Part B). Higher scores indicated higher frequency or severity of symptoms.

Laboratory evaluations of plasma biochemical profile and inflammatory markers

Venous blood samples were taken for measurement of standard hematological and clinical chemistry variables. Plasma biochemical profiles were assessed in all participants using standard in-hospital methods, including HbA1c, FBG, TG, TCHO, HDL-c, LDL-c, VLDL-c, BALP, and OC. The plasma concentrations of TNF-α, IL-6, IFN-γ, IL-17, IL-22, IL-25, and PGE₂ were measured by enzyme-linked immunosorbent assay kits (4A Biotech Co., Beijing, China). All assays were conducted in duplicate and in accordance with the manufacturer’s instructions. Biochemical indexes were measured by the hospital central diagnostic laboratory.

16S rRNA gene tag sequencing

Fecal samples (250 mg, wet weight) were immediately stored at −80° C until DNA extraction. The protocols for fecal bacterial DNA extraction, V3–V4 amplification, and sequencing were performed as described in our previous study [32]. Briefly, total fecal bacterial DNA was isolated using the Qiagen Mini Kit (Qiagen, Hilden, Germany), and molecular sizes were estimated by means of agarose gel electrophoresis. V3 and V4 regions were amplified by polymerase chain reaction using the following primer pair: 338F, 5’-barcode-ACTCCTACGGGAGGCAGCA-3’ and 806R, 5’-GGACTACHVGGGTWTCTAAT-3’. Polymerase chain reaction amplification was then performed as described in our previous study [33]. DNA libraries were constructed using kits provided by Illumina Inc. (San Diego, CA, USA), in accordance with the manufacturer’s instructions; DNA sequencing was performed using the Illumina MiSeq 2000 platform. Sequence data have been deposited under NCBI BioProject accession number PRJNA655208.

Quality-filtered sequences were clustered into unique sequences and sorted in order of decreasing abundance to identify representative sequences using UPARSE [34]. Sequences were grouped into operational taxonomic units (OTUs) using the clustering program VSEARCH [35]. The phylogenetic affiliation of each 16S rRNA gene sequence was analyzed using RDP Classifier (http://rdp.cme.msu.edu/) against the Silva (SSU132) 16S rRNA database, with a confidence threshold of 70%. 16S rDNA V3-V4 region amplicon sequencing generated 5,154,465 high-quality reads, with an average of 69,915 reads (range: 52,536–89,360) per sample.

Bioinformatics and statistical analysis

OTUs were grouped at different levels of taxonomy classification (from phylum to species) following the Bayesian approach with a 97% cutoff value. OTU-based α-diversity was estimated by calculating the following indexes: Chao1, observed OTUs, Simpson, and Shannon. Differences in these indexes among the three groups were analyzed by using Kruskal–Wallis tests. β-diversity was estimated by computing the Bray–Curtis dissimilarity, weighted UniFrac distance, and unweighted UniFrac distance; it was visualized using principal coordinate analysis, and the results were plotted using the WGCNA, stats, and ggplot2 packages in R software (Version 3.4.4, http://www.R-project.org/). To probe the microbial metabolism and predict metagenome functional content from marker genes, PICRUSt analysis was performed to explore differences in KEGG pathways among groups [36].

The variables were tested for collinearity with the CPI using variance inflation factor in the HH package in R. To ensure that the choice of the metric did not affect the results, the distances were calculated using three metrics: unweighted UniFrac distance, weighted UniFrac distance, and Bray–Curtis dissimilarity. Constrained ordination redundancy analysis was implemented in the vegan package in R, in accordance with previously described methods [37]. The significance of environmental variables was checked by forward selection using the ordistep function in the vegan package, with 99,999 permutations. A rarefied OTU table was used for all three metrics. UniFrac distances were calculated in QIIME 1.9.1 and Bray–Curtis dissimilarity was calculated in the vegan package in R. Permutational MANOVA (PERMANOVA) was performed on the UniFrac distances and Bray–Curtis dissimilarity using the adonis function in the vegan package in R, with 99,999 permutations. To test for confounding, adjusted PERMANOVA was conducted with T2DM_P and seven covariates in the model; the marginal effects were tested in accordance with previously described methods [38].

adonis(dis~P_status + Age + Sex + Abdominal_discomfort + Acarbose + Pioglitazone + Metformin + Statin, group, permutations = 99999)

Differences were tested at species, genus, family, order, class, and phylum levels. Taxa present in <10% of samples were removed; thus, 83 species, 113 genera, 45 families, 30 orders, 20 classes, and 12 phyla were included. Taxon abundances were compared between T2DM_P and controls, as well as between T2DM_P and T2DM_NP. Linear discriminant analysis effect size (LEfSe, http://huttenhower.sph.harvard.edu/lefse/) was used to characterize microorganismal features that distinguished gut microbiota specific to the T2DM_P group at multiple taxonomic levels, grade biomarkers according to statistical significance, and visualize the results in bar charts. Community analysis and differential abundance of OTUs were performed using the Statistical Analysis of Metagenomic Profiles software. Both analyses incorporated false-discovery rate (FDR) correction for multiple testing (FDR<0.05). To test for potential confounding, significant taxa were retested with adjustment for covariates, using a generalized linear model with negative binomial distribution and controlling for zero-inflation as appropriate, in the glmmADMB package in R:

Taxon~P-status Age+ Sex + Acarbose use

As previously described [39], identification of signature bacteria was performed by fivefold cross-validation analysis on a random forest model (R 3.4.1, randomForest 4.6–12 package) and receiver operating characteristic (ROC) curves were constructed using the pROC package.

The co-occurrence networks of gut microbiota were visualized using the igraph package in R. Closeness and eigenvector of the nodes were calculated to measure node centralities in each network, as described previously [18]. Genera with relative abundance > 0.05% were subjected to Spearman correlation analyses of their occurrence patterns. The bacterial correlations in samples from each group were computed based on the respective abundances of individual genera. Only strong correlations (P < 0.05 after FDR correction, correlation coefficient > 0.6) were visualized in this analysis. The microbiota networks of these three groups were compared; the numbers of edges were counted and centralities of nodes were determined in each co-bacterial occurrence network to quantify these differences.

Categorical variables were analyzed using the Chi-square or Fisher’s exact test. Continuous variables were compared using the Wilcoxon rank-sum test. Categorical variables were compared using Fisher’s exact test. The correlations of signature bacteria (identified by the random forest algorithm) with various clinical parameters (e.g., combined disease, metabolic profile, and inflammatory profile) were assessed using Spearman correlation analysis in R. Correlations were identified by Spearman’s rank correlation coefficient (significance thresholds were P < 0.05). Statistical analyses were performed using IBM SPSS Statistics for Windows V.20.0 (IBM Corp., Armonk, NY, USA). GraphPad Prism V.6.0 (GraphPad Software Inc., San Diego, CA, USA) and R were used for various analyses and preparation of graphs.