Temporal Variability of Oral Microbiota over 10 Months and the Implications for Future Epidemiologic Studies

Tooth loss and periodontal disease have been found to be associated with mortality (1, 2) and cancer (3, 4), which bolster the hypothesis that oral microbiota is associated with human health. Oral microbial differences have been observed in individuals with adverse health conditions such as pancreatic cancer (5–7), but the majority of studies have been conducted cross-sectionally, which precludes understanding the temporality of the effects, including an oral sample from a single timepoint.

Use of an oral sample from a single timepoint for epidemiologic studies of the microbiome will be effective if the microbial characteristics are fairly stable over time. However, if the oral microbiota is highly variable, collection of multiple oral samples would improve the power to detect associations with specific outcomes. Some previous studies considering the temporal variation in the oral microbiome within healthy individuals have concluded that oral microbiota is fairly stable over time (8–20) while others have suggested instability in oral microbiota (20–23). However, there have been limitations to many of the published studies including nonstandard assessments of stability, limited sample sizes, few samples collected per individual, and testing of only short time intervals.

Given these limitations, we conducted a study of 40 individuals who provided oral samples for a total of six visits, approximately every two months, using two collection methods. We used these samples to determine the temporal variability in oral microbiota using a variety of conceptual models. The results of this study have important implications for future oral microbiota studies and can be used to help determine the most appropriate sample size and number of sample timepoints for collection from each individual.

We recruited a convenience sample cohort of 40 adults working at the National Cancer Institute in the fall of 2014 by e-mail and word of mouth. Potential participants who were not planning on working at the National Cancer Institute for at least one year were advised not to participate. Participants received no monetary compensation and provided written informed consent. This study was approved by the Special Studies Institutional Review Board of the National Cancer Institute.

Oral samples were collected from participants approximately every two months for a total of six visits. In the 12 hours prior to the scheduled visit, participants were asked to refrain from dental hygiene procedures (including tooth brushing, flossing, and use of mouthwash or other dental rinse products), eating, drinking (other than water), chewing gum, consuming throat lozenges or candies, smoking cigarettes, or chewing tobacco.

At the baseline visit, participants provided two consecutive oral samples, both using the OMNIgene ORAL collection kit (DNA Genotek; Ottawa, Ontario, Canada). At all other visits, the participants provided a single oral sample using the OMNIgene ORAL kit. Participants spit into the kit until the saliva reached the fill line, the kit was closed to release the fixative, a cap was placed on the tube to replace the funnel, and then the tube was shaken. These samples were stored for up to 16 days at room temperature and transferred to the laboratory where they were shaken and then incubated at 50°C for 1 hour in a water bath, aliquoted, and frozen at −80°C. At the second and sixth visit, the participants provided an additional oral sample using Scope mouthwash after completing the OMNIgene ORAL collection. Using the protocol for the NHANES oral wash collection (24), the participant was instructed to swish 10 mL of Scope mouthwash from a sterile cup for 5 seconds and then gargle for 5 seconds and repeat for a total of 30 seconds. After the 30 seconds, the mouthwash was spit back into the sterile cup. The Scope mouthwash samples were stored for up to 2 days at room temperature and then transferred to the laboratory, aliquoted, and frozen at −80°C.

At each sample collection, participants filled out a short electronic questionnaire which included questions regarding adherence to the preparation procedures, oral health practices, and tobacco, alcohol, and antibiotic use. We also obtained demographic information and self-reported height and weight.

For DNA extraction, samples were grouped by collection method (OMNIgene ORAL or Scope mouthwash) and then batched into 23 samples for extraction, keeping all samples from each individual, of the same source material, within one extraction batch to limit technical variation due to batch effects. Each batch also included three quality control samples: (i) either an oral artificial community or a chemostat community (25); (ii) a blank; and (iii) an extraction duplicate of a randomly selected sample. For the OMNIgene ORAL samples, the vials were thawed at 4°C and then incubated for one hour at 50°C with vortexing every 20 minutes. An aliquot of 1,000 μL was then incubated for an additional 15 minutes at 75°C before DNA extraction. For the Scope mouthwash samples, the vials were thawed at 4°C and then 1,000 μL was transferred into a Pathogen Lysis Tube-L (Qiagen) and pelleted. After removing the supernatant, 894 μL of Buffer ATL (Qiagen) and 6 μL of Reagent DX solution (Qiagen) were added to each tube. Cells were lysed by vortexing the bead and buffer mixture for 20 minutes. For the artificial community and chemostat samples, similar procedures to the Scope mouthwash were used, except there was no pelleting step. After pretreatment, DNA was extracted from all samples using the DSP DNA Virus Pathogen kit (Qiagen) on a QIAsymphony instrument (Qiagen) following a modified Complex800_V6_DSP protocol.

For each sample, 20 ng of DNA, as measured by Quant-iT PicoGreen dsDNA (Thermo Fisher Scientific) quantitation, was split into 10 ng (1 ng/μL) aliquots for two separate PCR reactions. PCR was performed in 25 μL reaction volumes consisting of: 10 ng (10 μL) of DNA, 10 μL of 2.5× HotMasterMix (5 PRIME), 3 μL of MBG water, and 2 μL of the 5 μmol/L 16S rRNA V4 (515F/806R) barcoded primer mix that was composed of equimolar forward and reverse primer pairs. The primers were ordered from Integrated DNA Technologies. A total of 192 unique barcodes were used, resulting in 192 unique forward/reverse primer pair mixes that could be used for amplification. The 515F forward PCR primer sequence was AATGATACGGCGACCACCGAGATCTACAC TATGGTAATT GT GTGCCAGCMGCCGCGGTAA, consisting of the 5′ Illumina adapter, forward primer pad, forward primer linker, and forward primer. The 806R reverse PCR primer sequence was CAAGCAGAAGACGGCATACGAGAT XXXXXXXXXXXX AGTCAGTCAG CC GGACTACHVGGGTWTCTAAT, consisting of the reverse complement of the 3′ Illumina adapter, Golay barcode (defined by X, 12 bp barcode identifier generated specifically for this primer set), reverse primer pad, reverse primer linker, and reverse primer (26). Thermal cycling was performed on a GeneAmp PCR 9700 system (Thermo Fisher Scientific) where PCR conditions consisted of: 94°C hold for 3 minutes, denature at 94°C for 45 seconds, anneal at 50°C for 1 minute, extend at 72°C for 1 minute 30 seconds for 25 cycles, followed by a 72°C hold for 10 minutes.

Sample PCR replicates were then pooled and purified using a 1:1 AMPure XP (Beckman Coulter Genomics) ratio, performing the final elution in 30 μL of Buffer EB (Qiagen). Amplified sample libraries were quantified using Quant-iT PicoGreen dsDNA Reagent and up to 96, with unique barcoded adapters, were combined in equal amounts (100 ng each) and pools normalized to 10 nmol/L with Buffer EB for pooled sequencing.

For DNA sequencing, the Illumina MiSeq v2, 500 cycle kit was used (Illumina) following the manufacturer's protocol with modifications (27). Pooled libraries were diluted to 5 pmol/L in a serial dilution, and 25% denatured 5 pmol/L PhiX was spiked-in and added to the “Load Sample” well. A total of 3.4 μL of Index Sequencing Primer at 100 mmol/L, 3.4 μL of Read 1 Sequencing primer at 100 mmol/L, and 3.4 μL of Read 2 Sequencing Primer at 100 μmol/L was added to wells 13, 12, and 14 of the MiSeq sequencing cartridge. 2 × 250 bp paired end sequencing was performed on the MiSeq, with up to 96 samples per run.

From the Illumina MiSeq, the sequence files were demultiplexed using the bcl-to-fastq converter from CASAVA. The DADA2 pipeline 1.2.1 was used to generate sequence variant tables and phylogenetic trees based on paired-end sequence reads (28). For quality filtering, the sequences were trimmed at 10 and 220 base pairs for forward reads and 220 base pairs for reverse reads. The reads were then merged using the default mergePairs function in DADA2. After merging and error correction, sequence variants (i.e., 100% OTUs) were identified. After removal of chimeras, using removeBimeraDenovo, a total of 90% of the sequence reads were retained. The detected sequence variants were then aligned with the SILVA v123 database for taxonomy assignment (29). Sequences which aligned to mitochondria were removed. From 505 samples, 1,531 sequence variants were identified in total that included 16 phyla, 25 classes, 50 orders, 86 families, and 225 genera. An average of 84,787 sequences were generated per sample.

Using QIIME 1.9.1, relative abundances, alpha and beta diversity metrics were calculated (30). The relative abundance tables ranging from the phylum to genus level were generated for each sample without rarefaction. Alpha diversity metrics including Chao1 (31), observed species, Shannon index (32), and PD whole tree (33) were calculated on the basis of the average of 10 subsamples with rarefaction at 30,000 reads. These alpha diversity metrics measure richness and/or evenness of the bacteria in the samples. The beta diversity matrices, Bray-Curtis (34), unweighted UniFrac, and weighted UniFrac (35), in addition to the principal coordinates from those matrices were generated on the basis of rarefaction at 30,000 reads. Beta diversity measures the microbial composition of the samples and utilizes relative abundance data in the Bray–Curtis and weighted UniFrac matrices and presence/absence data for unweighted UniFrac, while phylogenetic tree data is included for both UniFrac matrices. Of the blank samples with any amplification, the median number of reads was 18 and all blank samples were removed after rarefaction at 30,000 reads. A principal coordinates analysis (PCoA) plot using the first two principal coordinates from each beta diversity matrix was generated to visually inspect clustering of the oral artificial community, the chemostat community, and the individual participant samples. The intraclass correlation coefficients (ICC) for the extraction duplicates were all greater than or equal to 0.94 for the alpha diversity metrics, 0.93 for the first two principal coordinates of the beta diversity matrices, and 0.96 for the relative abundances of the top four phyla. For the artificial community and chemostat samples, the coefficient of variation was less than or equal to 8% and 6%, respectively, for all alpha diversity metrics.

To quantify temporal variability, we calculated ICCs for the alpha diversity metrics and the square root of the relative abundance of the most abundant taxa at the phylum and genus level. For each microbiome metric, we estimated the between subject variance |\sigma _b^2$| and the within subject variance |\sigma _w^2$| using a linear mixed effect model from the lme4 package and estimated the ICC as |\sigma _b^2/( {\sigma _b^2 + \sigma _w^2} )$|⁠. The ICCs were calculated for the OMNIgene ORAL samples from the 6 visits (7 samples for each subject including 2 at baseline) and Scope mouthwash from 2 visits (2 samples for each subject). We calculated 95% confidence intervals (CIs) for each ICC based on 100 random bootstrap samplings. To diminish a potential inflation of the ICC for the repeated measurement of the OMNIgene ORAL samples at the first visit, we included a fixed effect for visit in the mixed model. We estimated the impact of adjusting the model individually for use of antibiotics, having a dental visit in the past 2 months, and not adhering to the preparation procedures prior to the visit. For the beta diversity metrics, a distance-based ICC was calculated on the basis of the three distance matrices. The 95% CIs for each distance-based ICC was also calculated on the basis of 100 random bootstrap samplings.

As another measure of temporal variability, we used a Dirichlet multinomial mixture model (DMM; ref. 36) to cluster the OMNIgene ORAL samples into community types by modeling the abundances of all taxa at each taxonomic level. The number of microbial community types was determined by minimizing the Bayesian Information Criterion (BIC). The OMNIgene ORAL samples were assigned to their community type based on the maximum posterior probability. We then characterized the temporal dynamics of microbiome community types using a Markov model (22). The transition probabilities of the Markov model were estimated from the maximum posterior probability using the microbial data from all subjects and a high probability of remaining in the same cluster would indicate temporal stability.

We also investigated the impact of the estimated ICC from the OMNIgene ORAL and Scope mouthwash samples for total sample sizes (N) required to detect differences between cases and an equal number of controls. We calculated the required sample size to detect a specific OR comparing the top quartile to the bottom quartile for a given microbial measurement from the estimated ICCs for the OMNIgene ORAL and Scope mouthwash samples assuming 80% power, a disease prevalence of 1%, and 1, 2, or 3 samples per individual. The required sample sizes were calculated for a range of P value thresholds to account for multiple testing.

All statistical analyses were conducted using R (version 3.3.2), except for the Dirichlet multinomial mixture model which was fitted by using MicrobeDMM (https://code.google.com/archive/p/microbedmm/downloads).

As seen in Table 1, the majority of participants were female (62.5%) and non-Hispanic white (67.5%). The average age was 39.2 years old (standard deviation 10.9) and few participants had ever smoked cigarettes (12.5%) but most had drank alcohol (97.5%). All participants completed all six study visits. Prior to each study visit, a small proportion of participants had used antibiotics in the past 30 days (ranging from 5.0% to 10.0%), had a dental visit in the past 2 months (ranging from 20.0% to 50.0%), or did not adhere to the preparation procedures (ranging from 12.5% to 27.5%; Table 2).

Overall, the OMNIgene ORAL samples from all six visits tended to cluster together using the first two principal coordinates from each beta diversity matrix although there was overlap between some participants. The two quality control samples, the oral artificial community and the chemostat community, tended to cluster separately from the study samples and separately from each other (Supplementary Figs. S1–S3).

For the alpha diversity measures, the ICCs were generally high, particularly for Chao1 and observed species with OMNIgene ICCs of 0.73 (95% CI: 0.65–0.81) and 0.74 (95% CI: 0.66–0.82), respectively, with slightly higher ICCs for Scope mouthwash. Compared with the alpha diversity metrics, the ICCs were generally lower for phylum-level relative abundances and Scope mouthwash had a low ICC for Fusobacteria (ICC 0.29; 95% CI: 0.00–0.83). For the beta diversity matrices, the highest ICC was observed for Bray–Curtis of the OMNIgene ORAL samples (ICC 0.70; 95% CI: 0.66–0.74; Fig. 1; Supplementary Table S1). When we included adjustment for use of antibiotics, having a dentist visit in the past 2 months, or not adhering to the requested preparation procedures, the ICCs were not materially affected (Supplementary Fig. S4). At the genus level, the ICCs for the OMNIgene ORAL kit were generally higher than for the Scope mouthwash, but the ICCs were still high for a number of genera in the Scope mouthwash sample such as Veillonella (ICC 0.67; 95% CI: 0.49–0.86) and Megasphaera (ICC 0.77; 95% CI: 0.64–0.89; Supplementary Fig S5; Supplementary Table S2).

For most levels from phylum to genus, the DMM clustered the OMNIgene ORAL samples into three community types (Supplementary Fig. S6). Figure 2 shows the temporal dynamics of phylum-level clusters with transition probabilities labeled on the plot. Cluster C appeared to be the most stable with a probability of 80.7% of remaining in Cluster C. The probability of remaining in Cluster A was only 59.5% and a proportion of Cluster A transitioned to Cluster B with a probability of 34.5%. However, Cluster A was the least common among participants with an average prevalence of 18.1%. Results for the transition probabilities at all taxonomic levels are presented in Supplementary Table S3.

Sample size requirements for future case–control studies considering the association between the microbiota and case status depend on the estimated ICC, the number of specimens collected per participant, the minimum detectable OR, and the P value threshold for statistical significance. From the OMNIgene ORAL ICCs, we estimated that given an ICC for the relative abundance of Actinobacteria was 0.46, with 1 specimen collected per participant, we would require 1,431 cases and 1,431 controls to detect an OR of 1.50 with a P value threshold of 0.05. The required sample sizes decrease for measures with higher ICCs, for larger ORs, and for increased specimens per participant. For example, with 3 specimens per participant and an OR of 3.50, a study would only require 92 cases and 92 controls for the relative abundance of Actinobacteria with a P value threshold of 0.05 (Table 3). Similar sample size requirements are presented for the OMNIgene ORAL samples with other P value thresholds (Supplementary Table S4) and from the Scope mouthwash ICCs (Supplementary Table S5).

In this study, the oral microbiota, measured from saliva and Scope mouthwash, appear to be relatively stable over a period of 10 months. The ICCs for the alpha diversity metrics were high with ICCs greater than 0.70 for both Chao1 and observed species from the OMNIgene ORAL and Scope mouthwash. The ICCs for the relative abundance of specific phyla and beta diversity matrices were generally lower and only the ICC for Bray–Curtis from the OMNIgene ORAL reached 0.70. When we adjusted for factors which may disrupt the oral microbiota, the ICCs were similar. From the DMM, we observed that the probability of remaining in a cluster was high, such that for participants in the clusters with a high prevalence in this study there was a high probability of remaining in those clusters. The DMM analysis provided further evidence for general temporal stability of the oral microbiota. Using the estimated ICCs, for large effect sizes, multiple samples per individual offer modest advanatages over a single sample, but for smaller effect sizes, multiple samples per individual offers larger gains in power.

A number of previous studies have evaluated the temporal variation in oral microbiota using saliva (8–10, 14, 15, 18, 20, 23, 37), tongue swabs (11, 13, 18, 20, 21), oral wash (11, 12, 16), and plaque samples (15, 18–20, 38). Similar to our findings, the majority of previous studies generally concluded that the oral microbiota is relatively stable over time, but these measures of stability varied between studies (8–20). It has been demonstrated previously that antibiotic use disrupts gut microbiota (39–41), but antibiotics appear to only have a modest effect on the oral microbiota (13, 42, 43). In our study, adjusting for the use of antibiotics did not appear to impact the estimates of temporal stability. In addition, adjustment for not adhering to the requested preparation procedures or visiting the dentist in the past month did not substantially improve the estimates of temporal variability.

This study has some limitations. The participants in our study are unlikely to be representative of the general population given the high education status which likely leads to healthier lifestyles. In addition, we did not collect dietary information so we were unable to determine how changes in diet may affect temporal variability. However, our findings are generally in concordance with previous literature that suggests that the oral microbiota may be relatively stable over time. In addition, we only collected samples over the course of approximately 10 months, so our study cannot evaluate longer term stability. Similarly, the sample size calculations were based on ICCs derived from sampling over 10 months and these required sample sizes may differ for samples collected over longer durations.

This study also has numerous strengths. We collected data on potential exposures which may disrupt the oral microbiota to evaluate whether these affected the temporal variability. We also used two different collection methods that may be used in large-scale epidemiologic studies to determine whether temporal variability differs by method of oral sample collection. Finally, we utilized multiple microbial measures and statistical methods used in different studies to evaluate oral microbial stability.

In conclusion, in this study of 40 individuals with oral samples collected every two months over the course of about 10 months, the oral microbiota appears to be relatively stable over time, but some measures, particularly the relative abundance of specific phyla and genera, were less stable. In planning future studies, investigators will need to take into account the temporal stability of specific measures if only collecting one sample from participants and for smaller effect sizes, collecting multiple samples per individual may provide some future gains in statistical power.

No potential conflicts of interest were disclosed.

Conception and design: E. Vogtmann, C.C. Abnet

Development of methodology: E. Vogtmann, C.L. Dagnall, K. Jones, R. Sinha, J. Shi

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): E. Vogtmann, C.L. Dagnall, A. Hutchinson, K. Jones, B.D. Hicks

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E. Vogtmann, X. Hua, L. Zhou, Y. Wan, S. Suman, B. Zhu, B.D. Hicks, R. Sinha, J. Shi, C.C. Abnet

Writing, review, and/or revision of the manuscript: E. Vogtmann, L. Zhou, C.L. Dagnall, B.D. Hicks, R. Sinha, C.C. Abnet

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): C.L. Dagnall, A. Hutchinson

Study supervision: E. Vogtmann, A. Hutchinson, C.C. Abnet

This study was supported by the Intramural Research Program of the National Cancer Institute (to E. Vogtmann, X. Hua, L. Zhou, Y. Wan, S. Suman, B. Zhu, C. Dagnall, A. Hutchinson, K. Jones, B. Hicks, R. Sinha, J. Shi, C. Abnet).

We would like to thank all of our participants.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.