Is oral microbiome of children able to maintain resistance and functional stability in response to short-term interference of ingesta?

A series of studies had focused on the ecological stability of human microbiome (Lozupone et al., 2012; Faith et al., 2013; Moya and Ferrer, 2016). Despite the continuous perturbation and the highly personalized composition within the human microbiome (Human Microbiome Project, 2012), healthy adults stably maintain their microbial communities in terms of space and time (Faith et al., 2013; Moya and Ferrer, 2016; Oh et al., 2016). This stability is proved to be critical for the well-being of human body (Lozupone et al., 2012). On the contrary, major shifts in microbial community composition are often related to diseases (Lynch and Pedersen, 2016). Of the many habitats throughout the human body, the bacterial community of the oral cavity is highly complex, constituting the second most diverse microbiome in the human body (Human Microbiome Project, 2012). These bacteria have a wide variety of functions, and many of them are important in maintaining oral health (Gao et al., 2018). Moreover, it is proved that some of the oral bacteria could colonize and persist to live in the gut, suggesting that oral bacteria might contribute to the emergence of chronic inflammatory diseases in gut (Lira-Junior and Bostrom, 2018). Oral microbiome dysbiosis is associated with oral diseases such as caries and periodontal diseases. Furthermore, a number of systemic diseases is associated with the dysbiosis of oral microbiome, including diabetes, rheumatoid arthritis and Alzheimer’s disease (Gao et al., 2018). Because oral cavity is the gateway to the human body for both food and air intake with frequent disturbances undergone, it is important to figure out whether oral microbiome can maintain stability in daily life within a short period. Considering the rapid and evenly distributed effects of beverages and the same standard schedule in kindergarten, we took three kinds of popular beverages as representatives of ingesta (orange juice [Group J], sugar-free tea [Group T], sugar-free liquid yoghurt [Group Y], representing the effects of sugar, tea polyphenols, and exogenous bacteria; and water [Group W] as control) to design a study to investigate the change of the oral microbial communities of children and their functional variations, using 16S rRNA gene, metagenomic sequencing and the Biolog technology, to further explore the characteristics of the human oral microbiome in maintaining stability in response to beverage intake (Fig. 1A). It was hypothesized that the species composition of oral microbiome could resist the interference of short-term beverage intake with a quick recovery to maintain functional stability. We first focused on the resistance of oral microbiome to beverage interference. Initially, we investigated the changes of oral microbiota in response to several factors using permutational multivariate analysis of variance (PERMANOVA) (Fig. S1A). According to the analysis, beverage intake was closely associated with the constitution of salivary microbial communities second to individual specificity. We also carried out a single nucleotide polymorphism (SNP)-based phylogenetic analysis to assess differences among the samples at the strain level. The sequences representing the top-9 abundant phylotypes from the same individual shared the most SNPs, so most samples from the same person were adjacent to one another in the phylogenetic analysis (Fig. S1B), indicating that an individual’s microbiome profile could be relatively stable over time though short-term beverage consumption altered the salivary microbial communities. This stable state was found in a number of sites in the healthy population (Human Microbiome Project, 2012; Lloyd-Price et al., 2017).

during which ~80 mL of the corresponding beverage was given to each child at four time points (09:30, 10:00, 11:00, and 14:00). They were instructed to only drink water during the 5 days except the 3 beverages. Whole stimulated saliva (~2 mL) was generated through passive tooth-brushing by trained dentists from Peking University School and Hospital of Stomatology at 14:30 each day (Fig. 1). The samples were collected in sterile, ice-chilled tubes and placed in liquid nitrogen for temporary storage. One milliliter of saliva was then aspirated and centrifuged at 10,000 × g for 10 min at 4°C. The pellets were stored at -80°C until DNA extraction.

DNA extraction and 16S rRNA V3-V4 region sequencing
DNA was extracted using the TIANamp Bacteria DNA Kit DP302 (Tiangen, Beijing, China) according to the manufacturer's protocol. DNA purity was evaluated based on the A260/A280 ratio using a NanoDrop 8000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, NC, USA). DNA integrity was verified by agarose gel electrophoresis. A negative control containing only buffer was included during DNA extraction and quantification. DNA was stored at -80°C in Tris-EDTA buffer before use. To enable the amplification of the V3-V4 region of the 16S rRNA gene and the addition of barcode sequences, we designed unique fusion primers based on the universal primer set, 338F (5′-GTACTCCTACGGGAGGCAGCA-3′) and 806R (5′-GTGGACTACHVGGGTWTCTAAT-3′), along with barcode sequences. The PCR mixtures contained 1 μL of template DNA (30 ng), 2 μL of each forward and reverse primer (10 µM), 4 μL of deoxynucleotides (2.5 mM), 5 μL of 10× Pyrobest buffer, 0.3 μL of Pyrobest DNA polymerase (2.5 U/μL, TaKaRa Code: DR005A), and 35.7 μL of ddH 2 O in a 50 μL reaction volume. Thermal cycling consisted of an initial denaturation step at 95°C for 5 min, followed by 25 cycles of denaturation at 94°C for 30 s, annealing at 57°C for 30 s, and extension at 72°C for 40 s, with a final extension step at 72°C for 4 min. PCR amplicons were subsequently purified using the PCR Product Rapid Recovery and Purification Kit (ZhiAng Biotechnology, Changchun, China). Quality was assessed by agarose gel electrophoresis after the first two steps.
An amplicon library for high-throughput sequencing on the Illumina MiSeq platform was constructed and subsequently quantified (KAPA Library Quantification Kit KK4824) according to the manufacturer's instructions.

qPCR for sample selection and metagenome sequencing
The human DNA in samples collected on days 1, 3, and 5 was quantified by real-time read length, 125 bp). Low-quality and host reads were removed, and the remaining high-quality reads were used for further analyses.

Profiling of 16S rRNA gene sequencing data
Using the Quantitative Insights into Microbial Ecology pipeline (Caporaso et al., 2010), the raw sequences were processed to concatenate reads into tags according to their overlaps, after which reads belonging to each sample were separated with barcodes and low-quality reads were removed. The processed tags were clustered, and chimaeras were removed prior to analysis. Operational taxonomic units (OTUs) were defined at 97% sequence similarity to taxa by matching to the Greengenes database (DeSantis et al., 2006).

Taxonomic assignment of metagenomic data
A total of 32,605 sequenced bacterial genomes available in GenBank as of June 2015 were downloaded from the NCBI database (http://www.ncbi.nlm.nih.gov/).
Clade-specific marker genes were identified, and all genomes were clustered into phylotype groups according to gene similarities using the procedure reported by a previous study (Segata et al., 2012).

Mantel test
A Mantel test was carried out to assess the correlation between 2 matrices, the 16S rRNA OTU profile and the metagenome abundance profile. The test was processed by R software (3.3.1, vegan package; permutations = 999, P value = 0.001).

KEGG annotation
The catalogue of non-redundant genes was translated to putative amino-acid sequences and aligned with the proteins/domains in the KEGG database using usearch ublast algorithm 9 (Version v9.2.64, e-value < 1e-5, query coverage > 0.70). Each protein was assigned to one KO based on the highest-scoring annotated hit(s) containing at least one segment pair scoring >60 bits. The Wilcoxon rank-sum test (P value < 0.05) was used to identify significantly changed KOs and KEGG level-3 functional pathways.

PERMANOVA analysis
Based on the profiles of 16S rRNA genes and metagenomic phylotypes, PERMANOVA assessed by 999 permutations for Bray distance was used to estimate the variation effects (R package vegan). In the analysis to compare the effects of beverages, the profile of metagenomic phylotypes was split by beverage, and the PERMANOVA interpretation of randomly-ordered time factors (999 permutations, Bray distance) was calculated for 1000 iterations. The PERMANOVA interpretation value of the time factor in its actual chronological order for each beverage was then compared with the distribution of the randomly-ordered time interpretation values.

SNP-based phylogenetic analysis
We selected the 9 most abundant phylotypes, each accounting for >2% of all known microbes generated from the metagenomic data, and then aligned the reads from all 60 samples to the genomes of strains representing the selected phylotype groups using the Burrows-Wheeler Aligner 0.7.15 (Li and Durbin, 2009). SNPs were called using Samtools 0.1.19. (Li et al., 2009) SNPhylo software (Lee et al., 2014) was used to remove low-quality data and generate a phylogenetic tree with a minor allele frequency of 0.05.

Analysis of alterations in genus counts as indicated by 16S rRNA
Samples with complete time sequences were selected to determine the retained genus counts from day 1 to day 5. Divided according to the four beverages, the mean values of genus counts were calculated and visualized in the form of a Sankey diagram using R software (3.3.1, networkD3 package).

PCA
R software (3.3.1, ade4 package) was used to perform PCA on 16S rRNA genus profiles for the different beverages.

STEM and correlation analysis
The relative abundance data of all the phylotypes were used to cluster time-series dynamic patterns by STEM (Bar-Joseph, 2006). The maximum number of clusters was 8, and the maximum unit change between time points was set to 2. Significance was calculated using the Wilcoxon rank-sum test (P value <0.05) between days 1 and 3, and between days 1 and 5.

Biolog assays and data analysis
Whole stimulated saliva (~1mL) of another 3 children was generated through passive tooth-brushing by trained dentists at 9:00 (Table S4). The samples were collected in sterile, ice-chilled tubes and placed in ice for temporary storage. One milliliter of saliva was then aspirated from each sample and all of them were mixed together. The mixed saliva sample was centrifuged at 500 rpm for 30s to remove food debris, with the supernatants resuspended in 9mL PBS (0.01M, pH=7.2-7.4) and vortexed thoroughly for 60s. Beverage solutions with concentration gradients of 100%, 10%, 5%, and 2% were obtained by diluting the original beverage with PBS (0.01M, pH=7.2-7.4). After measuring ODs at baseline, 99μL saliva and 1μL beverage were inoculated into each well on 96-well Biolog MT2 plates (Biolog Inc.) which contained redox dye and a buffered nutrient medium optimized for a wide variety of bacteria. For the sake of repeatability, the experiment for each concentration was performed for three times. The plates were incubated in a 5% CO 2 incubator at 37 C° for up to 4 days. The OD value at wavelength of 590 nm (OD590) in each well was recorded after 5h, 24h, 48h and 72h using a Biolog micro-station and software kit (Biolog OmniLog version 4.1).
All Biolog assay data were analyzed using SPSS 20.0 software (IBM, Armonk, NY). Differences in OD values among the 4 time points were evaluated using one-way ANOVA test. The threshold for statistical significance was set at P<0.05.    Vertical bars represent microbiome samples in day 1, day 3, and day 5 with metagenomic data; bars indicate relative abundances colored by microbial genera with an abundance higher than 0.01.