Microbiome succession with increasing age in three oral sites

Characteristics of recruited participants

We collected samples from the gingival crevicular fluid (GCF), saliva (SAL) and tongue back (TB) sites of healthy people of different ages. The subjects were divided into five age groups: Group A (age 11-15), Group B (age 18-20), Group C (age 28-32), Group D (age 38-45) and Group E (age 50–65) (Supplementary Table 1), representing adolescents, youth, the middle-aged and the elderly. All samples were analyzed by 16S rRNA sequencing; detailed information is introduced in the Material and Methods section. After quality filtering, more than 15.05 million clean reads were harvested corresponding to a mean of 84,090 effective tags and 372 OTUs per sample (Supplementary Tables 2 and 3). The rarefaction curve indicated that the OTUs of all samples tended to saturate as the sequence number increased (Supplementary Figure 1), indicating that the OTUs in the data cover most of the reads in the samples.

Lower within-sample and higher between-sample microbial diversities in aging mouth

To compare the differences in oral microbiota among different age groups, the α-diversity of the oral microbiota relative to age and oral site was first explored. Microbial OTUs within each oral site were analyzed by ANOVA, and the results showed that the OTUs within GCF or TB had marked reductions with age (Supplementary Figure 2, p=0.0029 and 0.0001, respectively). Meanwhile, the bacterial OTUs in GCF decreased gradually with age, whereas there was no difference between groups B-E for TB; however, the OTU was significantly larger in group A than in the older groups. The Faith_PD diagrams also demonstrated the distinction of microbial within-sample diversity between different age groups within each site (Figure 1A). The variation trends of Faith_PD in GCF, SAL and TB were similar to that of OTUs. Additionally, an intersite comparison was performed on the bacterial α-diversity within each age group, and the cross-site pattern of Faith_PD revealed differentiation in α-diversity between various oral sites. The α-diversity differed between at least two sites except group A (Figure 1B), further demonstrating the distinction of bacteria in different oral sites and age-related variation.

Subsequently, the β-diversity (between-sample) of the oral microbiota in response to aging was investigated. The principal coordinate analysis (PCoA) by the Bray_Curtis distance illustrated a remarkable shift in microbial distribution between age groups (Figure 1C), and the p-values obtained by ADONIS (permutational MANOVA) analysis of the bacterial variations in GCF, SAL and TB were <0.001, 0.003 and <0.001, respectively. The comparison of the first principal components by Bray_Curtis showed that microbial β-diversity significantly increased with age (Figure 1D), except for that of SAL and TB, which declined in the >50 years group. Furthermore, the principal coordinate analysis (PCoA) by the Bray_Curtis distance illustrated the cross-site variation in microbial distribution, and the p-values obtained by ADONIS (permutational MANOVA) analysis were <0.001 for Group A-D and 0.036 for Group E (Supplementary Figure 3). The first principal components analysis by Bray_Curtis also showed the β-diversity differed at least between two sites (Supplementary Figure 4).

To verify the changes in between-sample diversities, the software STEM (Short Time-series Expression Miner) was used to cluster bacteria into modules and calculate their p-values. The results showed that 3, 1 and 1 microbial clusters had a p-value of <0.001 (non-random) in GCF, SAL and TB, respectively (Figure 1E). Analyzing the tendency of variation in each module, the results revealed a rising trend with various patterns. Meanwhile, the TB and SAL modules increased sharply at first but then slowly while a module at the GCF site decreased in the last age group (Figure 1E and Supplementary Figure 5, Supplementary Table 4). Checking the tendency of variation and bacterial composition of each module, Eubacterium sulci, Stomatobaculum longum, Alloprevotella rava, Porphyromonas gingivalis, Abiotrophia defectiva, Porphyromonas endodontalis, Parvimonas micra, and Dialister pneumosintes increased their abundance with age in GCF; Dialister pneumosintes, Tannerella forsythia, Treponema medium, and Slackia exigua increased their abundance with age in saliva; and only Alloprevotella rava was found with enhanced abundance in TB. In summary, the above results indicate that α-diversity has a tendency to decrease with age, and a cross-site distinction existed. The β-diversity has an increasing trend in three sites, and some bacteria contributed to the increase in β-diversity.

The alteration of oral microbiota with the increase in age

To clarify the effects of age on oral microbiota, we then compared the structure of oral microbiome at the genus and species levels. At the genus level, Streptococcus, Neisseria, Haemophilus, Prevotella and Fusobacterium were dominant in all three sites (Figure 2A). Comparing the abundance of bacterial genera between the age groups, 88, 44 and 47 genera were found to have marked variation with age in the GCF, SAL and TB sites, respectively (Supplementary Figure 6). Meanwhile, 39, 6 and 7 genera enhanced or reduced their contents continuously with age in the GCF, SAL and TB sites, respectively (Figure 2B and Supplementary Figure 7, Supplementary Table 5); i.e., their trends were not changed in the age groups. Among all the genera shifting up or down gradually in all three sites, Delftia, Escherichia/Shigella, and Herbaspirillum were reduced in both the GCF and SAL sites while the Rothia content increased in the two sites. Similarly, Bradyrhizobium, Brevundimonas, Pelomonas and Serratia decreased their abundance with age in both the GCF and TB sites (Figure 2B and Supplementary Table 5). There was no overlap between the SAL and TB sites. At the species level, Haemophilus parainfluenza, Veillonella dispar, Porphyromonas catoniae and Prevotella jejuni were dominant, with Haemophilus parainfluenza having the highest abundance in all three sites (Supplementary Figure 8). Among the species of bacteria included in all age groups, 75, 24 and 49 species were found to have significant age-related variation in the GCF, SAL and TB sites, respectively (Supplementary Figure 9). Meanwhile, 14, 2 and 4 species varied continuously with upward or downward trends in the GCF, SAL and TB sites, respectively (Supplementary Figure 10 and Supplementary Table 5). In GCF, Acinetobacter spp., Actinomyces spp., Anaeroglobus geminatus, Bifidobacterium longum subsp. infantis, Propionibacterium acnes, Pseudomonas beteli, Sphingomonas xenophagum, Treponema maltophilum, Oribacterium sinus, Prevotella nanceiensis, Solobacterium moorei, and Stomatobaculum longum changed their contents with age continuously; in SAL, Faecalibacterium prausnitzii and Fusobacterium mortiferum had a trend of variation; and in TB, Eubacterium sulci, Feacalibacterium prausnitzii, Porphyromonas gingivalis and Tannerella forsythia changed their contents gradually (Supplementary Table 5). Among these species with trends, only Feacalibacterium prausnitzii decreased in abundance with age in more than one site (SAL and TB) (Supplementary Table 5). The phyla to which the different genera belong were displayed by the Sankey diagram, and the results demonstrated that the different genera harboring GCF and SAL with prevalence all belong to the phyla Firmicutes and Proteobacteria (Figure 2C); meanwhile, other genera in TB belonging to Bacteroidetes were also found. The above results suggest that the microbial composition changes with age, and the patterns of variation in different oral sites were divergent. Afterward, we performed an intersite comparison of the microbial community within each age group at the genus and species levels. At the genus level, 12, 56, 28, 18 and 9 genera were differentially distributed across the sites for the age groups A to E, respectively (Figure 3A and Supplementary Figure 11, Supplementary Table 6); the number of differential genera demonstrated a trend of a bell-shaped curve (increasing first and then decreasing) with age (Figure 3A and Supplementary Figure 11). The cross-site differential genera included Prevotella, Rothia, Solobacterium, etc., but only Solobacterium was shared by all age groups (Supplementary Figure 12, Supplementary Table 6). At the species level, 10, 72, 43, 17 and 25 species showed significantly varied intersite distribution with age, respectively (Supplementary Figures 13, 14, Supplementary Table 6). The number of cross-site differential species also formed a bell-shaped curve in aging mouth. The intersite differential species included Prevotella spp., Actinomyces spp., and Solobacterium spp., etc., but only Actinomyces odontolyticus and Solobacterium moorei were shared by all age groups (Supplementary Figure 12 and Supplementary Table 6). The phyla to which the intersite differential genera belong were displayed by the Sankey diagram, and the results demonstrated that the different genera in Group A all belonged to the phyla Actinobacteria, Firmicutes and Proteobacteria (Figure 3B), whereas more differential genera belonged to the phyla Bacteroidetes and Teneriquetsas as the host age increased. The results above indicate the separation in the patterns of bacterial composition in various sites of aging mouth, and the number of intersite differential bacteria followed bell-shaped trend in response to aging.

Variation of the “core microbiome” in GCF, SAL and TB of aging mouth

Next, the “core microbiome” of the five age groups within each site, i.e., bacteria shared by all the age groups, was obtained by Venn analysis at the genus and species levels. At the genus level, 93, 90 and 85 genera were included in the “core microbiome” of GCF, SAL and TB, respectively (Figure 4A). At the species level, the number changed to 93, 119 and 96, respectively (Figure 4B). In all three sites, the number of genus and species contained in the “core microbiome” all decreased gradually with age (Figure 4A, 4B). The cross-site “core microbiome,” i.e., the bacteria shared by the GCF, SAL and TB within each age group, was then obtained by Venn analysis. At the genus level, the cross-site “core microbiome” consisted of 300, 162, 148, 108 and 82 genera for Groups A-E, respectively (Supplementary Figure 15). At the species level, the number was 237, 162, 139, 110 and 97, respectively (Supplementary Figure 15). The number of genera and species contained in the cross-site “core microbiome” decreased with age, and 73 genera and 85 species were shared by all age groups (Figure 4C). The results above suggest that the intra- and intersite “core microbiome” reduced their bacterial diversity significantly in aging mouth, which is consistent with the results of α-diversity.

Correlations of the age- and site-related bacteria

To study the correlation between the age-differentiated bacteria within each site, |Spearman correlation| ≥0.7 and q value≤0.01 were employed to analyze the prevalent taxa (≥0.02%) at the genus and species levels. The results showed that close correlations existed among the different bacteria in the age groups (Figure 5 and Supplementary Figure 16), and clear variations in the pattern of bacterial correlations were observed in each oral site, indicating the shift in bacterial correlation in response to aging and oral site variation (Figure 5 and Supplementary Figure 16).

Using differential bacteria to distinguish age groups

To explore the feasibility of using oral bacteria to assess mouth age or maturity, the random forest classifier model was employed for distinguishing different age groups. The data preparation of the random forest model was performed with 5-fold cross-validation, i.e., the original data was randomly divided into 5 sub-groups in each age group firstly; then, using a random sub-group as the test set, the remaining 4 sub-groups were used as training groups. The differential bacteria species with trends of GCF were first evaluated by the random forest model (Supplementary Figure 17), and the results showed that application of trending differential bacteria as a biomarker to distinguish Group A, B, C, D, and E could reach an accuracy of 0.998, 0.809, 0.668, 0.675 and 0.956, respectively (Figure 6A). The SAL and TB sites had many fewer different bacteria species than GCF, so the accuracy rates were relatively lower; the value was 0.804, 0.677, 0.835, 0.642 and 0.620 for Group A-E for SAL, respectively; and 0.831, 0.597, 0.515, 0.663 and 0.631 for Group A-E for TB, respectively (Supplementary Figures 18, 19). Meanwhile, Actinomyces odontolyticus, Oribacterium sinus and Solobacterium moorei were used for linear fitting of the age groups B, C, and D, and much smaller p-values were obtained (Figure 6B).

DNA extraction, library construction and sequencing

Extraction of DNA

Total DNA was extracted with CTAB/SDS method. 1% agarose gels was used for checking the DNA concentration and purity. DNA was diluted to 1ng/μl using sterile water based on the obtained concentration.

Library construction and sequencing

Genome DNA from all the samples was used as amplification templates. PCR primers were from the V3-V4 region of 16S rDNA, forward primer, 5'-ACTCCTACGGGAGGCAGCA-3'; and reverse primer, 5'-GGACTACHVGGGTWTCTAAT-3'. All PCR reactions were carried out in 30μL total volume with 15μL of Phusion® High-Fidelity PCR Master Mix (New England Biolabs), 0.2μM of forward and reverse primers, and about 10ng templates DNA. Thermal cycling included an initial denaturation step at 95°C for 5 min, followed by 30 cycles of denaturation at 95°C for 30s, annealing at 50°C for 30s, and elongation at 72°C for 40s. Finally 72°C for 7 min.

Then the PCR products were purified with GeneJET Gel Extraction Kit (Thermo Scientific) and qualified by electrophoresis on 2% agarose gel, samples with single amplification product were chosen for further experiments. The library was sequenced on an Illumina Hiseq 2500 platform at Novogene Company (Beijing, China).

Data processing

The quality of the Raw sequence was filtered and analyzed using the next-generation microbiome bioinformatics platform (QIIME2 version 2018.6 pipeline) for data preprocessing. A software package included with Usearch was used to identify the exact sequence variants (ESV). α- and β-diversity analyses were performed using the phyloseq software in R package. α-diversity was calculated by Faith's Phylogenetic Diversity, Shannon index, and observed OTUs. Principal coordinate analysis (PCoA) was analyzed based on the Bray_Curtis distance, Bray-Curtis variability is an indicator used to measure differences in taxonomic composition in ecology. PCoA is a method to explore and to visualize similarities or dissimilarities of data. It sorts the data through a series of eigenvalues and eigenvectors, selects the main eigenvalues in the first few places, finds the most important coordinates in the distance matrix, and perform a rotation on the data matrix. PCoA does not change the mutual positional relationship between the sample, but only changes the coordinate system. PCA is based on the sample similarity matrix (such as Euclidean distance) to find principal components, while PCoA is based on dissimilar distance matrix (other distances than Euclidean distance, including binary_jaccard, bray_curtis, unweighted_unifrac, and weighted_unifrac distance) to find principal coordinates. Statistical differences in β-diversity were analyzed using permutation-arranged multivariate analysis (PERMANOVA, R function Adonis (vegan, 999 permutations)). Usearch was used to cluster the Effective Tags of all samples, and the sequences cluster into 372 operational taxon units (OTUs) with 97% identity, and then perform representative sequences of OTUs according to the database Ribosomal Database Project (The species annotation is RDP version 11.5). At the genus and species levels, the differences between age groups were analyzed. The “pheatmap” for the heat map, and the average relative abundance in the group were displayed with a histogram (R package: “ggplot2”). Correlation analysis was performed at the genus and species levels, all the differential genera and species identified in the variation analysis were pooled for spearman correlation analysis of their correlation with all the bacteria and clinical parameters after 0.002% and 30% screening (R package: “Hmisc”). The genera and Species with padj<=0.05, |cor|>=0.5 were shown with a network map (software cytoscape, R package “pheatmap”). To further explore the relationship between oral microbiota and age, we used a random forest algorithm to predict the classification of samples. Random forest regression was performed with 1,000 regression trees based on 5-fold cross-validation, 80% of the samples were randomly selected for model training and the remaining 20% were used for validation. The predicted result is displayed by the roc curve (R package: “proc”). In the random forest prediction classification algorithm, the contribution of different species can be known. The top species were selected and the random forest algorithm prediction classification was performed.

Ethics statement and samples collection

All saliva samples used in this study were collected from the Stomatological Hospital of Shandong University (Jinan, Shandong). Please refer to Supplementary Table 1 for more information of the sampling patients. This study was reviewed and approved by the Ethics Committee of the Stomatological Hospital of Shandong University. Informed consent was signed by all volunteers prior to inclusion in the study. All methods are implemented in accordance with the relevant guidelines and regulations.

Saliva specimens were collected 2 hours after a meal: mouth was rinsed with sterile double distilled water for 3 times firstly before collecting saliva, then the saliva samples were collected with 1.5μL centrifuge tubes. To collect a gingival crevicular fluid sample, 2 paper points (Henry Schein, Almere, Netherlands) were inserted to the bottom of the gingival pocket for 10 s afer supragingival plaque removal, and the paper points after sampling were immersed in sterile phosphate buffer saline. For tongue back sample, a cotton swab was use to collect fluid at the soft tissue of the tongue back, and the cotton swab after sampling were immersed in sterile phosphate buffer saline. All samples were immediately transferred to the laboratory within 20 minutes and stored at -80 °C for later usage.

Availability of data

The sequencing data will uploaded to the NCBI website with a BioProject ID PRJNA number and available at https://www.ncbi.nlm.nih.gov/sra/PRJNA.