Human oral mucosa and oral microbiome interactions following supragingival plaque reconstitution in healthy volunteers: a diet-controlled balanced design proof-of-concept model to investigate oral pathologies

Figures & data

Figure 1. (a) Experimental balanced design for determination of RNA expression profiles of dental plaque microbiome and oral mucosa, and dental plaque microbiome abundance data, in healthy volunteers. (b) Overall analytical design.

Legend. (a) Supra- and superficial sub-gingival plaque (SSP) was collected by scraping a dental scaler across all teeth of healthy volunteers (N = 8). Oral mucosal epithelial cells from the malar oral mucosa were exfoliated from the same subjects using Oral CDx brush. Only one brush biopsy was performed at any time point followed by at least a 3-week interval before performing a second brush biopsy to minimize the risks associated with healing. No food other than yogurt and tea diet was consumed for at least 10 h before start and during the day of the experiment. The proposed sample collection strategy is designed to enable comparisons of dental plaque microbial and oral mucosa gene expression at two daily time points (i.e., T = 0 h and T = 8 h for dental plaque and oral mucosa epithelium). In addition, baseline dental plaque microbial gene expression (T = 0 hr) can be compared between Week 1 and Week 3. The design implies that Group 2 would not have been subjected to an AM brush biopsy procedure at Week 1. However, the design may not be repeated with the same cohort by reversing the roles of Group 1 and Group 2 for paired comparisons, to avoid repeated mucosal injury. A total of 20 RNA samples for RNA-Seq analysis can be prepared in Week 1 (12 dental plaque and 8 oral mucosa) and 20 more during Week 3. Patients were divided into Group 1 (N = 4) and Group 2 (N = 4) to allow a balanced design of sample collection at Week 1 and Week 3, allowing appropriate time for healing process of the oral mucosa and to limit experimental bias associated with routine hygiene and diet. Plaque samples were also used for 16S metagenomic sequencing. (b) Legend. Supra- and superficial sub-gingival plaque (SSP) and buccal mucosa brush biopsies (BMB) were collected from eight healthy volunteers at T = 0 h (AM) and T =©8 h (PM). HOMINGS v3-v4 16S rRNA gene sequencing (V3-V4 region) was used to identify genus and species level taxa through ProbeSeq program. Alpha- and beta-diversities AM vs. PM and AMavg vs. PM based on SSP relative abundance data were calculated using PRIMERv7. Additionally, RNA was extracted from SSP and BMB samples and sequenced using RNASeq Ion Proton™ System. Resulting FASTQ files were trimmed using Trim Galore v0.6.7. SSP trimmed fastq files were then subjected to KneadData for removal of contaminating human RNA reads. Discarded SSP human reads and BMB sample reads were then separately aligned using Spliced Transcripts Alignment to a Reference (STAR) algorithm with NCBI’s RefSeq GRCh38.p13 as reference. Aligned reads were run through FeatureCounts software to attain counts of coding sequences using annotated GRCh38.p13 file. SSP meta-transcriptomic FASTQ files were trimmed using Trim Galore and analyzed using HMP Unified Metabolic Analysis Network (HumanN2). Additionally, SSP trimmed FASTQ files were analyzed©using Simple Annotation of Meta-transcriptomes by Sequence Analysis 2 (SAMSA2) pipeline. Downstream analysis was completed using various libraries in R v4.3.1.

Table 1. Healthy volunteers’ cohort demographics and multiple sample collection.

Criteria	Combined^a	AM BMB^b	PM BMB^c	AM SSP^d	PM SSP^e
Healthy Volunteer (M/F)f	8 (4/4)	8	8	16	8
Week 1 - not processedg		2			1
Week 3 - not processedh				4	1
Weeks 1 & 3 - processed		6 (2/4)	8 (4/4)	12 (6/6)	6 (3/3)
Agei:
Median	31	35	31	31	41.5
Mean	39.75	37.33	39.75	37.25	41.67
Standard deviation	16.61	13.49	16.61	15.21	18.40
Range	21–68	21–54	21–68	21–68	21–68

Sample collection and processing tallies: At Week 1, four SSP samples from Group 1 and four from Group 2, collected in the morning, were sequenced successfully. Four BMB samples collected in the morning for Group 1 and four collected in the afternoon for Group 2 were also processed successfully. Three of four possible SSP samples that were collected in the afternoon.

from Group 2 were successfully sequenced. At Week 3, four AM SSP samples were not processed successfully, precluding a powered paired Week 1 vs. Week 3 comparison. Meanwhile, four BMB samples collected in the morning for Group 2 and four collected in the afternoon for Group 1 were sequenced. Additionally, three SSP samples of four collected in the afternoon from Group 1 were sequenced. Sequencing data were used for AM vs. PM comparisons only.

BMB: buccal mucosa brush biopsy; SSP: supra- and superficial sub-gingival plaque.

M: male; F: female.

aCombined number of healthy volunteer participants.

bCounts of T = 0 h (AM) BMB samples.

cCounts of T = 8hrs (PM) BMB samples.

dCounts of T = 0 h (AM) SSP samples. Four healthy volunteers (Group 2) had two samples each of AM SSP from weeks 1 and 3 (data were averaged or not). Four healthy volunteers (Group 1) had only AM samples processed from Week 1.

eCounts of T = 8hrs (PM) SSP samples.

fHealthy volunteers (left panel) and anticipated AM and PM samples to be collected and processed for sequencing (right panel).

gWeek 1 samples not processed due to technical error.

hWeek 3 samples not processed due to technical error.

iMedian, mean, standard deviation and range of healthy volunteers’ ages.

Figure 1. (Continued).

Table 2. Alpha- and beta-diversity analysis.

Comparison	Monte-Carlo p-value
a. Alpha-diversity of T = 0 h (AM) plaque vs. T = 8 h SSP samples and the average of AM samples (AMavg) vs. PM samples (α = 0.05)
AM vs. PM	0.653
AMavg vs. PM	0.195
b. Beta-diversityT = 0 h (AM) SSP vs. T = 8hrs SSP samples and the average of AM samples (AMavg) vs. PM samples (α = 0.05)
AM vs. PM	0.8694
AMavg vs. PM	0.7794

(a) Shannon and Simpson indices were generated using PRIMERv7 (PRIMER-E Ltd., Ivybridge, UK). T = 0 h (AM) healthy volunteer (HV) SSP ProbeSeq data (n = 12) were compared to T = 8 h HV SSP ProbeSeq data (n = 6). Mann–Whitney U-tests were used to determine significance of differences at α = 0.05 using XLSTATv2016. Alpha-diversity of averaged HV AM samples (AMavg; n = 6 datapoints) compared to PM samples (n = 6) was also determined.

(b) Significance of beta-diversity differences was determined by PERMANOVA in the comparison AM group (n = 12) vs. PM group (n = 6) using PRIMERv7 (PRIMER-E Ltd., Ivybridge, UK). Species and genera relative abundance data were square root transformed and converted into a Bray-Curtis similarity matrix. PERMANOVA was completed using a mixed model with unrestricted permutation of raw data, 9,999 permutations and type III partial sum of squares. ‘Time’ (AM or PM) was used as the fixed factor in this design and significance level was set at α = 0.05.Beta-diversity of averaged HV AM samples (AMavg; n = 6 data points) compared to PM samples (n = 6) was also determined.

Figure 2. Heatmap and PCA plots comparing AM to PM SSP data as determined by DESeq2 analysis using SAMSA2 pipeline. a. organism Heatmap b. organism PCA c. function Heatmap d. function PCA.

Legend. (a and b) Heatmaps and principal component analysis (PCA) of organism and (c and d) function analysis results using DESeq2 in the SAMSA2 pipeline. Heatmaps were plotted using RColorBrewer, ggplot2 and pheatmap libraries in R. In each plot, ‘control’ refers to AM (T = 0 h) supra- and superficial sub-gingival plaque (SSP) samples, while ‘experimental’ refers to PM (T = 8 h) samples. Heatmaps (a and c) are shown clustered using complete Euclidean clustering. The color scheme goes from blue (low abundance) to red (high abundance). (b and d) PCA plots are shown with AM (control) samples in red and PM (experimental) samples in blue.

Figure 2. (Continued).

Figure 3. Combined histograms showing top 30 most abundant taxa and functions using the SAMSA2 pipeline on dental plaque bacterial RNA-Seq data a. organism level b. function level.

Legend. Stacked graphs showing the top (a) organisms (bacterial taxa) and (b) functions within all dental plaque bacterial RNA-Seq data comparing the AM (T = 0 h; left) group to the PM (T = 8 h; right) group. For (a) and (b) the top graph shows the relative activity of the total samples, and the bottom plot shows the total read composition per sample.

Figure 3. (Continued).

Figure 2. (Continued).

Figure 2. (Continued).

Table 3. Significant bacterial gene expression-derived pathway enrichment determined by LEfSe analysis from AM to PM.

Pathway Name^a	PWY-ID^b	Species^c	Abundance^d	logLDA score^e	p-value^f	PMIDs^g
Super pathway of pyridoxal 5-phosphate biosynthesis and salvage	PWY-0845	Unclassified	NA	1.26	0.01212	31586394; 32286295
Super pathway of L-serine and glycine biosynthesis I	SER-GLYSYN-PWY	Corynebacterium matruchotii***	Increase; PM>AM, 4.30FC	1.14	0.03385	31570555; 31586394
L-arginine biosynthesis III via N-acetyl L-citrulline	PWY-5154	Unclassified	NA	1.43	0.03385	21247782; 31182718; 31586394; 32409559
S-adenosyl L-methionine cycle I	PWY-6151	Neisseria elongata	Increase; PM>AM, 1.11FC	1.24	0.03959	29102113; 31586394
6-hydroxymethyl dihydropterin diphosphate biosynthesis I	PWY-6147	Corynebacterium matruchotii***	Increase; PM>AM, 4.30FC	1.60	0.03959	31586394; 33173145
Super pathway of glycolysis pyruvatede hydrogenase TCA and glyoxylate bypass	GLYCOLYSIS-TCA-GLYOX-BYPASS	Unclassified	NA	1.22	0.03959	31586394
TCA cycle V 2-oxoglutarate ferredoxin oxidoreductase	PWY-6969	Unclassified	NA	1.85	0.03959	31586394
Gondoate biosynthesis anaerobic	PWY-7663	Aggregatibacter aphrophilus**	Decrease; AM>PM, 3.30FC	1.21	0.03959	21910872; 22223671; 22369781; 30899518; 31037175; 31586394; 33654282;
5-amino imidazole ribonucleotide biosynthesis I	PWY-6121	Lautropia mirabilis*	Increase; PM>AM, 2.94FC	1.93	0.03959	20622961; 32286295; 31586394
UNINTEGRATED	UNINTEGRATED	Aggregatibacter aphrophilus**	Decrease; AM>PM, 3.30FC	2.82	0.03959	21910872, 22223671; 30899518; 31037175; 31586394
Pyridoxal 5-phosphate biosynthesis I	PYRIDOXSYN-PWY	Unclassified	NA	1.10	0.04274	31586394
Adenine and adenosine salvage III	PWY-6609	Unclassified	NA	1.56	0.04296	31586394

Pathway analysis annotation results:.

PWY-0845 is a pathway wild-type Escherichia coli K-12 substr. MG1655 is known to possess, which involves the metabolism of glycogen and amino acids. Glycogen metabolism can regulate inflammatory responses in patients with sepsis and macrophages can activate glycogen metabolism regulating the macrophage-mediated inflammatory response. SER-GLYSYN-PWY involves the production of L-serine and glycine. E. coli K-12 substr. MB1655 and humans are known to possess this pathway. L-serine supplementation has been shown to reduce inflammation in the lungs and increase the survival rate of mice infected with Pasteurella multocida. Taxa known to possess PWY-5154 include Cytophaga hutchinsonii, Prevotella ruminicola, Tannerella forsythia and Bacteroides thetaiotaomicron. L-arginine biosynthesis contributes to complex pathogen–host interactions. Clinical trials have shown supplementation of L-arginine in diet reduces inflammation. Many oral hygiene companies have started using arginine in toothpaste as it may enrich alkali-producing bacteria, preventing tooth decaying pathogens compared to sodium-fluoride toothpaste. PWY-6147 is induced by the 6-hydroxymethyl-dihydropterin diphosphate biosynthesis III (PWY-7359) pathway. This pathway has been shown to be enriched in a group of overweight individuals characterized as metabolically unhealthy. PWY-7663 involves the production of a long-chain fatty acid, gondoate by Aggregatibacter aphrophilus. It is suggested that long chain fatty acids produced by bacteria may alter cell membranes causing modified cell signaling and gene expression. Long chain fatty acids may be associated with persistent infection of opportunistic pathogens in respiratory samples in a cohort of children with cystic fibrosis. A. aphrophilus has been previously described as an opportunistic pathogen. PWY-6121, a pathway of Lauratropia mirabilis, represents a key intermediate step in the biosynthesis of purine nucleotides and thiamine. This pathway has been shown to be a differential pathway between healthy subjects and subjects with inflammatory bowel disease and a differential pathway between subjects with inflammatory bowel disease and subjects with colorectal cancer in a study investigating human intestinal bacteria. It has been previously shown that thiamine deficiency in gastric cancers can lead to Wernicke’s encephalopathy, a neurologic syndrome characterized by paralysis of the extraocular muscles, impaired balance/coordination, and acute confusion. Highly conserved pathways include GLYCOLYSIS-TCA-GLYOX-BYPASS, PWY-6969, PYRIDOXSYN-PWY, PWY-6609, and PWY-6151. GLYCOLYSIS-TCA-GLYOX-BYPASS pathway is a super pathway essential for energy metabolism. PWY-6969 involves aerobic respiration to generate energy and is the foundation for biosynthesis. PYRIDOXSYN-PWY is key for amino acid and carbohydrate metabolism. PWY-6609 is essential for the salvage and recycling of adenine to produce inosine monophosphate resulting in the synthesis of nucleotides. PWY-6151 is crucial in regulating S-adenosyl-I-homocysteine in bacteria. This pathway is known to control pathogenicity and can influence the rate of E. coli bacterial growth in bird models. gNCBI PubMed identification numbers (PMIDs) related to pathways annotations are shown.

aMetaCyc pathway abundance name.

bMetaCyc pathway identifier.

cTaxa contributing to the pathway. Mann-Whitney U-test significance levels of SAMSA2 taxa abundance are denoted with asterisks (*): p < 0.05 as *; p < 0.01 as **; p < 0.001 as ***.

dSAMSA2 abundance change from AM to PM. Abundance refers to taxa reads as a percentage of total reads. NA is not applicable. FC is fold change. PM>AM denotes mean PM abundance was greater than the mean AM abundance by x FC and vice versa.

eThe log Linear Discriminant Analysis (LDA) score.

fp-value of the ‘one-against-all’ strategy for LDA.

Note: Unclassified reads per HumanN2 refer to ‘coding reads that are highly diverged from the corresponding sequence in the pangenome database, spanned multiple genes, and/or are from apparent non-coding regions’.

Data availability statement

All data and accompanying files are available via the Translational Research Lab Github repository (https://github.com/mbeckm01/RNASeq).