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

ABSTRACT

Changes in the oral microbiome may contribute to oral pathologies, especially in patients undergoing cancer therapy. Interactions between oral microbiome and oral mucosa may exacerbate inflammation. We determined whether probiotic-controlled plaque formation could impact proximal oral mucosa gene expression profiles in healthy volunteers. A 3-weeks balanced sample collection design from healthy volunteers (HVs) was implemented. At Week-1 plaques samples and labial mucosa brush biopsies were obtained from HVs in the morning (N = 4) and/or in the afternoon (N = 4), and groups were flipped at Week-3. A fruit yogurt and tea diet were given 2-4hrs before sample collection. mRNA gene expression analysis was completed using RNA-Seq and DESeq2. Bacterial taxa relative abundance was determined by 16S HOMINGS. Bacterial diversity changes and metabolic pathway enrichment were determined using PRIMERv7 and LEfSe programs. Alpha- and beta-diversities did not differ morning (AM) vs. afternoon (PM). The most affected KEGG pathway was Toll-like receptor signaling in oral mucosa. Eighteen human genes and nine bacterial genes were differentially expressed in plaque samples. Increased activity for ‘caries-free’ health-associated calcifying Corynebacterium matruchotii and reduced activity for Aggregatibacter aphrophilus, an opportunistic pathogen, were observed. Microbial diversity was not altered after 8 hours plaque formation in healthy individuals as opposed to gene expression.

KEYWORDS:

• Host–microbe interaction
• dental plaque
• oral mucosa
• oral microbiome
• gene expression

Introduction

The oral microbiome consists of over 800 bacterial species collectively, with approximately 150 to 250 species present in oral cavity per individual [1–3]. Oral cavity habitats include saliva, tongue, buccal mucosa, teeth surfaces, gums, palate and dental plaque presenting core bacterial profiles with modest distinctions at each location [4]. Dental plaque biofilms form in a sequential order with 80% of early colonizers represented by Streptococci able to bind directly to tooth substrates [5]. Secondary colonizers include Corynebacterium species bridging Streptococci with other species Fusobacterium, Capnocytophaga and Leptotrichia spp [6]. Plaque biofilm–host interactions are complex among subgingival bacterial species thriving on inflammation [5]. Subgingival species such as P. gingivalis and T. denticola can evade host immune responses through Toll-like receptor (TLR) pathways and damage oral tissues [7]. These species can also down-regulate epithelial cell interleukin expression to evade phagocytic neutrophils and Th1 cells [5]. Oral microbiome may thus undergo dysbiosis disrupting healthy interactions, thereby leading to lasting local or systemic deleterious effects [8–10].

Dysbiosis may promote periodontal disease, oral ulcers, oral candidiasis, dental caries and altered wound healing [11]. Furthermore, poor oral hygiene may result in plaque accumulation and support growth of opportunistic bacteria. In cancer patients, dysbiosis can lead to increased bleeding and infections [12]. Previous studies have reported increased proportions of oral bacterial species associated with periodontal disease in patients undergoing chemotherapy [13–18]. These species include Actinobacillus actinomycetemcomitans, Porphyromonas gingivalis and Fusobacterium nucleatum. Proper oral hygiene by brushing teeth using a fluoride toothpaste and avoiding spicy, acidic, or sugary foods can reduce the risk for oral pathologies in cancer patients [19].

In addition, oral toxicities associated with cancer treatment may be dose-limiting and may impact a patient’s quality of life and overall survival [12]. Cancer treatment-associated oral toxicities affect approximately 40% of patients receiving chemotherapy and up to 80% of receiving conditioning therapy prior to hematopoietic cell transplantation (HCT), thereby leading to dental alterations, neurological disorders, salivary alterations, dysgeusia, infections, increased bleeding, osteonecrosis and cancer-therapy induced oral mucositis (CTOM) [20–23]. Moreover, the severity of CTOM may also be impacted by a dysbiotic oral microbiome and poor oral hygiene resulting in delayed wound healing [24,25]. Studies have tested probiotics to mitigate CTOM, and one meta-analysis was able to determine that probiotics can mitigate moderate and severe oral mucositis [26,27–29].

We have previously shown that oral microbiome profile changes in HCT patients over a one-year time-period involved an increase in Gammaproteobacteria abundance in patients with more severe CTOM [30]. Haverman et al. showed that yeast Candida spp. and Porphyromonas gingivalis inhibit wound healing [31]. However, ulcers in oral mucositis rarely occur within mucosal surfaces of keratinized gingiva surrounding the teeth [22]. Therefore, one may assume distal effects through the release of proinflammatory molecules from the bacteria present in dental plaque [31,32]. Resistance to CTOM may also be explained by an effective fibrinolysis process occurring at the tooth-giving interface despite neutrophil infiltration which correlates with CTOM severity [32–35].

Thus far, there is no study that simultaneously investigated host–microbe interactions in various oral pathologies by analyzing the interaction between mucosal epithelium located at the vicinity of dental plaque and dental plaque bacteria using transcriptomics and metagenomics approaches. Routine daily tooth brushing will cause regular reconstitution of dental plaque. We hypothesized that a plaque-inducing diet would remove some variability during plaque biofilm formation in healthy volunteers which would allow to better characterize changes in the gene expression profile. Using a cohort of healthy volunteers, we hypothesized that diet-controlled and experimentally induced plaque formation can identify gene expression signatures specific to oral mucosa and dental plaque.

Methods

Study population and controlled diet

The study was approved by the Institutional Review Board of Carolinas Medical Center-Atrium Health, Charlotte, NC (IRB file #01–16-09E). Study participants signed an informed consent form. Subjects were healthy volunteers (HVs; N = 8, 4 females and 4 males) [mean age (SD) = 39.75 (16.61)]. Determination of healthy dental and periodontal status included the following criteria: have at least 24 teeth, have less than 4 mm pocket depths and have less than 10% Bleeding on Probing (BOP). At least 10 h prior to probiotic diet breakfast, participants were to have a meal (dinner) and perform fluoride-based tooth brushing and flossing and record times of the meal, brushing and flossing. In this study, 2–4 h before sample collection, subjects were offered two starch-free and blended particle-free Greek-style yogurt (Chobani, Norwich, NY) with fruit jam on the bottom (strawberry or blueberry), containing only natural non-GMO ingredients and without artificial flavors or preservatives. The yogurt contains the probiotic bacterial species Streptococcus thermophilus, Lactobacillus bulgaricus, Lactobacillus acidophilus, Lactobacillus bifidus, Lactobacillus casei and Lactobacillus rhamnosus and the sources of sugar (e.g., cane sugar from the fruit jam) necessary to boost plaque formation [36,37]. To minimally control for bacteria attaching to oral soft tissues such as the buccal mucosa, unsweet oolong, black, or green teas were offered as a beverage (Ten Ren tea Co. Ltd., Taiwan). These teas have been previously shown to minimally affect bacterial attachment to immortalized human gingival fibroblast-1 (HGF-1) cells [38].

Dental plaque sample collection

Supra- and superficial sub-gingival plaque (SSP) samples were collected by scraping a dental scaler from the buccal side of the teeth at T = 0 h (morning; AM) and T = 8 h (afternoon; PM), 2-4hrs after a diet of fruit yogurt and tea was given to boost plaque formation. SSP samples were obtained by scraping a dental scaler across supra- and superficial sub-gingival buccal surfaces. Four double-sided sterile scalers were used, namely, one scaler end per one mouth-quadrant. Plaque from all four scaler heads were pooled into one microcentrifuge tube containing 1 mL cold PBS-NF by vigorously rolling the scaler against the inside surface of the sterile tube for 10 seconds. Samples were transported on ice, homogenized by pipetting, and a 60 µL aliquot was removed and frozen at −80°C for subsequent bacterial DNA isolation.

Bacterial genomic DNA was isolated using QIAamp DNA Mini Kit (Qiagen, Venlo, the Netherlands) as previously described [39,40] following manufacturer’s instructions. Isolated DNA was frozen at −20°C for further microbiome profiling analysis. Microbiome profiling was performed using Human Oral Microbiome Identification Next Generation Sequencing (HOMINGS) after amplification of the 16S rRNA gene (V3-V4 region), using a modified MiSeq NGS method (Illumina, Inc., San Diego, CA) [41]. Sequences were barcoded and saved electronically. Oral taxa identification and relative abundance were determined using ProbeSeq program to identify sequences matching up to 638 species probes and 129 genus probes [30,42,43].

Buccal mucosa sample collection

Buccal mucosa brush biopsies (BMB) were collected from HVs at T = 0 h (AM) or T = 8 h (PM) using a sterile OralCDx brush (CDx Diagnostics, Suffern, NY) following manufacturer’s instructions. Briefly, after the mouth was rinsed with 5 mL tap water, the buccal mucosa was dried with a sterile gauze. An OralCDx brush was rotated 10 times firmly against the buccal surface, then immediately submerged into a microcentrifuge tube containing 250 µL of cold nuclease-free PBS (PBS-NF) (ThermoFisher Scientific, Waltham, MA). The handle of the Oral CDx brush was clipped with sterile, RNase-treated scissors prior to closing the tube lid and samples were transported for immediate processing.

Total RNA isolation

RNA was extracted from BMB and SSP samples using miRNeasy Mini Kit (Qiagen, Venlo, the Netherlands) following the manufacturer’s protocol with minor modifications. Briefly, cell pellets were pipet-resuspended in 250 µL of nuclease-free PBS and lysed with 700 µL of Trizol LS (ThermoFisher Scientific, Waltham, MA) and 140 µL of chloroform. After centrifugation, the aqueous phase was removed and further processed using Qiagen miRNAeasy Mini Kit procedure. Total RNA yield, composition and quality were assessed using NanoDrop (ThermoFisher Scientific, Waltham, MA) and Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA).

Library preparation & RNA sequencing

RNA was converted into cDNA using Ion Total RNA-Seq Kit for AB Library Builder system (ThermoFisher Scientific, Waltham, MA). The library quality was determined on an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA), and RNA concentrations were determined from the analysis profiles and Qubit fluorometer (ThermoFisher Scientific, Waltham, MA). Barcoded libraries were pooled on an equimolar basis, loaded onto Plv2 chips and sequenced on the Ion Proton^TM System (Thermo Fisher Scientific, Waltham, MA).

Bacterial plaque alpha- and beta-diversities

Raw ProbeSeq data was transformed into relative abundance (RA) data. Shannon and Simpson indices were generated using PRIMERv7 (PRIMER-E Ltd., Ivybridge, UK). Mann–Whitney U-tests were used to determine the significance of alpha diversity differences AM vs. PM (α = 0.05) using XLSTAT_v2016 (Addinsoft (2020) XLSTATv2016 Statistical and Data Analysis Solution. New York. (https://www.xlstat.com)). Permutational multivariate analysis of variance (PERMANOVA) was performed AM vs. PM using PRIMERv7 (PRIMER-E Ltd., Ivybridge, UK) as previously implemented [39,42,43]. Species and genera RA data were squared root-transformed and converted into a Bray-Curtis similarity matrix. Metrics were checked to ensure the transformation would not influence the results using the ‘ecotraj’ package in Rv4.3.1 [44]. 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 a fixed factor, and the significance level was set at α = 0.05.

Buccal mucosa brush biopsy RNA-seq analysis

RNA-Seq fastq files (fastqs) were checked for quality control using FASTQC and trimmed using Trim Galore (TrimGalore). QIAseq miRNA NGS 3’ adapter (5’-3’ AACTGTAGGCACCATCAAT) was trimmed from all samples. BMB sequences were aligned to Human RefSeq GRCh38.p13 and a STAR index was created using human genome for alignment [45,46]. The alignment was completed for all brush biopsies’ data using the splice-aware aligner. FeatureCounts program was used to count genes from the annotated files using GRCh38.p13 genomic annotation file to count multi-mapping coding sequence features with a minimum mapping quality score of 10 [47].

The DESeq2 Rv4.3.1 program was used with AnnotationDbi, org.Hs.eg.db, Gage, GageData and Pathview programs to analyze the differential mRNA expression AM vs. PM [48–56]. The significance of associated pathways was determined using the Generally Applicable Gene-set Enrichment (GAGE) algorithm with R Pathview package where gene sets are binned, and differential expression tests were performed based on one-on-one sample comparisons AM vs. PM using log fold change. Two-sample t-tests were used to determine if gene sets were differentially expressed relative to background genes. A global p-value was derived based on the negative log sum of p-values from each one-on-one comparison [54].

STRING online tool was used to visualize significant protein–protein interactions with a confidence score of 0.150 removing disconnected nodes [57]. All nodes were checked using Human Protein Atlas database to ensure gene expression was represented in the proximal digestive tract (i.e., the oral mucosa, salivary gland, esophagus and tongue) [58].

Dental plaque RNA-seq analysis

SSP RNA-Seq data were processed using Human Microbiome Project Unified Metabolic Analysis Network 2.0 (HUManN2) to explore gene families and UniProt Reference Cluster (UNIREF) pathways present [58]. Species identified using HUManN2 were used with RefSeq and annotation files using NCBI [45]. Corresponding fasta files were merged and an index was built. Bowtie2 was used to align the files to the large index [59]. SSP RNA-Seq samples data were merged by time and normalized using counts per million (CPM).

Linear Discriminant Analysis (LDA) Effect Size (LEfSe) online program was used to determine significant pathway abundances using time as the input ‘class’, and sample as the input ‘subject’ [60] with the ‘one-against all’ strategy and α = 0.05. Using the KneadData program, human reads were removed from trimmed plaque fastqs using the GRCh38.p13 reference (KneadData). Removed reads were processed using the same approach as for buccal mucosa brush biopsy samples by running fastqs through GAGE algorithm and R Pathview package to determine the significance of associated KEGG pathways. The Simple Annotation of Metatranscriptomes by Sequence Analysis 2.0 (SAMSA2) pipeline was then utilized. SAMSA2 employs DIAMOND for annotation based on NCBI RefSeq database along with python and R for analysis of aggregation files on functions and organisms present in SSP data remaining reads [45,61,62].

Results

Diet-controlled balanced design for oral mucosa and dental plaque sample collection

Due to the excessive invasiveness of collecting two oral mucosa brush biopsies on the same day from one individual, we used a paired balanced sample collection design of two groups of four HVs, i.e., two males and two females each: Group 1 [AM week 1] and [PM week 3], Group 2 [PM week 1 and AM week 3] (Figure 1a). Table 1 summarizes the demographics, sample collection and sequencing tallies of our HV cohort (N = 8).

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.

Bioinformatics analysis

An overview of the bioinformatics analysis strategy is presented in Figure 1b.

Figure 1. (Continued).

Alpha- and beta-diversities from HOMINGS data

ProbeSeq of 16S sequence data from SSP samples identified 304 species and 49 genera of the possible 638 species and 129 genus probes. An integer adjusted average count of species and genera detected per sample are presented in Supplemental Table S1. AM vs. PM differences in alpha-diversity were not significant (Table 2a). None of the beta-diversity comparisons, AM vs. PM or AM_avg vs. PM, were significant (p > 0.05) (Table 2b). A retrospective power analysis with MicroPower R package for PERMANOVAs [63] based on bacterial abundance showed that we had 90% power to detect an effect size of 0.022 AM vs. PM (AM samples: n = 12; PM samples: n = 6) and 0.015 AM_avg vs. PM (n = 6 AM_avg and n = 6 PM samples).

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.

Buccal mucosa mRNA gene expression

We identified 47 differentially expressed genes when comparing buccal mucosa mRNA gene expression of AM vs. PM samples. A protein–protein interaction network was generated with STRING confidence score of 0.150, returning 27 interacting proteins (Supplemental Figure S1a). The most significant gene ontology (GO) term was ‘Protein insertion into mitochondrial inner membrane’ at a false discovery rate of 0.0183. Significantly under-expressed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways from DESeq2 included olfactory transduction (hsa04740), oxidative phosphorylation (hsa00190), glycerophospholipid metabolism (hsa00564) and ribosome (hsa03010) (Supplemental Figure S2a-d). The only significantly increased KEGG pathway for the PM group was Toll-like receptor-related chemokine signaling (hsa04062) (Supplemental Figure S2e). PathView was able to show STAT1 as downregulated in the Toll-like receptor signaling pathway (Supplemental Figure S2e). Furthermore, 76 Gene Ontology (GO) biological processes (BP) were found to be significant for the lower expressed genes. Significant GO BP for significantly overexpressed genes included DNA replication (GO:0006260), response to virus (GO:0009615), chemokine-mediated signaling pathway (GO:0070098) and DNA-dependent DNA replication (GO:0006261). All significant KEGG pathways and GO biological processes result from DESeq2 analysis of buccal mucosa sequencing data are shown in Data File 1.

Dental plaque human mRNA gene expression

Using DESeq2, human reads removed from dental plaque sequencing data corresponded to 18 differentially expressed genes comparing AM samples to PM samples. A protein–protein interaction network showing 12 of the 18 significant genes at a confidence score of 0.150 is shown in Supplemental Figure S1b. Further analysis revealed 27 significant KEGG pathways including one for under-expressed genes (olfactory transduction) and 26 for over-expressed genes including: protein processing in endoplasmic reticulum (hsa04141), ribosome (hsa03010), phagosome (hsa04145), and endocytosis (hsa04144). The top five significant pathways are presented in Supplemental Figure S3a-e. Additionally, 44 GO biological processes were determined to be significant for under-expressed genes and 419 for over-expressed genes. Top GO terms for under-expressed genes included homophilic cell adhesion (GO:0007156), negative regulation of retinoic acid receptor signaling (GO:0048387), synaptic transmission (GO:0019226), and cell–cell adhesion (GO:0016337) while GO terms for over-expressed genes included protein catabolic process (GO:0030163), innate immune response (GO:0045087), interaction with host (GO:0051701), symbiosis, encompassing mutualism through parasitism (GO:0044403) and interspecies interaction between organisms (GO:0044419). All significant KEGG pathways and GO biological process terms are listed in Data File 2.

Meta-transcriptomic bacterial identification and gene expression

SAMSA2 plaque RNA-Seq analysis identified 7,613 bacterial taxa including 679 genera and 6,934 species. Organism-specific SAMSA2 analysis identified two significant species, Bacillus aryabhattai and Endozoicomonas numazuensis (p = 0.03241) when comparing AM vs. PM samples. A heatmap and principal component analysis (PCA) plot of organisms AM vs. PM samples are presented in Figures 2(a,b). PCA of organism aggregation determined that 24% variance of the first principal component reflected AM vs. PM differences (Figure 2b). A combined plot of the top 30 most abundant taxa is presented in Figure 3a.

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.

Using SAMSA2 with bacterial functional aggregation, the heatmap indicated downregulation of many bacterial genes in PM group with high expression of gallidermin family proteins and PTS glucose transporter subunit IIBC among three PM samples. Additionally, T9SS C-terminal target domain-containing protein and ornithine decarboxylase SpeF were upregulated in four PM samples (Figure 2c). PCA of bacterial functional aggregation in SAMSA2 showed 37% of variance in the first principal component represented differences AM vs. PM (Figure 2d). A heatmap and PCA plot of functional plaque components comparing AM to PM samples are presented in Figure 2(c,d). We identified 18,963 bacterial functional proteins including hypothetical proteins, OmpA family protein, and porin OmpA having the highest base mean and being the most abundant when comparing AM plaque samples to PM plaque samples (Figure 3b and Data File 3). Functional SAMSA2 analysis identified nine significant differential bacterial genes comparing AM vs. PM samples (p < 0.05) including dinucleotide-utilizing protein, lysozyme inhibitor, lysyl aminopeptidase, DUF4743 domain-containing protein, sodium:proton exchanger, sugar tyrosine-protein kinase, flagellar basal body rod modification protein, lauroyl acyltransferase and multispecies: HD domain containing protein. A combined functional plot of the most abundant functional proteins related to these 30 taxa is presented in Figure 3b. The organism and functional DESeq2 results can be found in Data File 3.

Figure 3. (Continued).

Figure 2. (Continued).

Figure 2. (Continued).

Meta-transcriptomic metabolics

Using HumanN2 analysis on meta-transcriptomic dental plaque LEfSe data comparing AM samples to PM samples, no MetaCyc pathways were enriched in the AM group. However, 12 significant MetaCyc pathways were enriched in PM group (Table 3). The 12 significant bacterial metabolic pathways in the PM group included super pathway of pyridoxal 5 phosphate biosynthesis and salvage (PWY-0845) and arginine biosynthesis III via N-acetyl-L-citrulline (PWY-5154) by uncultured species, super-pathway of L-serine and glycine biosynthesis I (SER-GLYSYN-PWY) and 6-hydroxymethyl-dihydropterin diphosphate biosynthesis I (PWY-6147) by Corynebacterium matruchotii, and S-adenosyl-L-methionine cycle I (PWY-6151) by Neisseria elongata [64]. Table 3 lists enriched bacterial metabolic pathways that may also exist in humans, as determined by HumanN2 analysis. HumanN2 also identified unintegrated Haemophilus parainfluenzae which had the largest sum across all samples, followed by unintegrated Fusobacterium nucleatum and Rothia dentocariosa.

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’.

Discussion

We described an approach to investigate gene expression signatures relevant to interactions between dental plaque bacteria or their products, with oral mucosa epithelial cells in healthy individuals. We analyzed the data for 1) buccal mucosa gene expression, 2) dental plaque human gene expression, 3) dental plaque bacterial gene expression, and 4) dental plaque bacterial abundance upon 8-h plaque reconstitution.

The pathways we identified as enriched for human genes of dental plaque, namely, hsa04141, hsa03010, and hsa04145, likely stem from immune and epithelial cells. Endoplasmic reticulum stress has been associated with amelogenesis imperfecta, an inherited disease with characteristics such as reduced enamel and unusually small or discolored teeth that are prone to rapid breakage and wear [65,66]. The ribosome-related pathway hsa03010 has also been implicated for its role in diseases due to its ability to activate p53 pathways that can result in cell cycle arrest and apoptosis [67].

Furthermore, the phagosome related pathway hsa04145 plays a role in inflammation and degradation of biological materials and involves OmpA family of proteins [68]. Reduced levels of these proteins indicate a successful bacterial breakdown [69,70]. In our study, OmpA family proteins had the highest abundance (base mean = 8449.74), and lysozyme inhibitors were the most significant bacterial function using SAMSA2 (p = 1.77×10⁻¹¹) in the AM vs. PM comparison. In cases of excessive presence of reactive oxygen species, the phagosome may lead to cytotoxic effects on periodontal tissues [71]. In addition, species associated with periodontal disease progression and infected root canals, such as Treponema denticola, can influence the release of pro-inflammatory cytokines inducing host inflammation [71–73]. In plaque samples, we were able to show tight clustering of T. denticola,T, vincentii, T. medium, and T. maltophilum having higher abundance in PM vs. AM group. This could explain the enriched phagosome pathway hsa04145 in the PM group.

In human buccal mucosa, we determined that toll-like receptor chemokine signaling pathway hsa04062 was upregulated. Interestingly, the toll-like receptor chemokine signaling pathway is induced when oral mucositis occurs after bacteria colonize in ulcers [74]. Toll-like receptors have been previously implicated for their involvement in periodontitis since bacteria invading oral tissues can cause the release of pro-inflammatory cytokines leading to host tissue degradation [75]. Indeed, this pathway is likely upregulated due to the induction of plaque formation in the PM group after a diet that encourages the constitution of plaque.

We observed that the olfactory transduction KEGG pathway hsa04740 was suppressed in both buccal mucosa and dental plaque in the AM vs. PM comparison. A recent study by Orecchioni et al. showed that bacterial lipopolysaccharides can boost RNA expression of olfactory receptors in bone marrow-derived macrophages [76]. It was also previously suggested that the microbiome, diet, oxidative stress, and inflammation may have the ability to influence olfactory receptors or their ligands leading to disease [77–79].

The bacterial pathways identified using HumanN2 included PWY-0845 known to involve non-pathogenic Escherichia coli K-12 substr. MG1655 (p = 0.012118) [64]. This metabolic pathway can activate UDP-glucose-P27 signaling and upregulate signal transducer and activator of transcription 1 (STAT1) [80]. Persistent Candida infections of the mucosal membranes have been proposed to cause mucosal ulcers by way of a gain-of-function mutation in STAT1 [81]. STAT1 was found downregulated in Toll-like receptor signaling pathway of the PM group. Studies have shown that green teas can downregulate the expression of STAT1 and that STAT1 may be a promising target for anti-inflammatory treatments [82,83].

We also found C. matruchotii increased in abundance in PM dental plaque samples. C. matruchotii has been previously shown to be associated with oral health and dental plaque homeostasis [84,85]. Specifically, C. matruchotii has higher abundance in the plaque than in other oral niches and may counteract the negative effects of early plaque colonizers, such as Streptococcus species, since it is seen as a marker of a caries-free state [84]. A study by [6] observed Haemophilus and Aggregatibacter species only appeared near Corynebacterium species when Corynebacterium species were bound to Streptococcal early colonizers of dental plaque [6]. We identified H. parainfluenzae using RefSeq and HumanN2 with both showing the largest sum across all samples. H. parainfluenzae is normally prevalent in dental plaque and upper respiratory tract, exhibiting health associated or pathogenicity properties depending on the strains involved [86–88].

Our results are likely related to experimentally induced plaque formation through yogurt and green tea diet stimulation in healthy volunteers. Overall, plaque reconstitution had a large impact on mRNA expression of dental plaque microbiome. The pathways/proteins identified from human genes/dental plaque and their involvement in inflammation provide insight into the importance of good oral hygiene to reduce the deleterious effects of plaque formation on human health.

Limitations

The healthy volunteer cohort in this pilot study was small. We also encountered technical issues for this longitudinal multiple sampling design involving RNA analysis. Extending this study to a larger sample size-independent repeat cohort would increase accuracy. In addition, although a healthy core microbiome is relatively stable, sampling procedures being 3 weeks apart may have introduced unaccounted variability in the microbiome data.

Conclusions

In this study, we showed that there were no significant changes in bacterial diversity upon dental plaque reconstitution over an 8 hr timeframe. However, significant bacterial and human mRNA expression changes occurred at oral sites following a plaque inducing diet. Using this methodology will allow for a better understanding of the impact of changes in the oral microbiome on cancer and non-cancer related oral pathologies and contribute to the development of probiotic prophylaxis. Geng et al. [89,90];

Authors’ contributions

JLM and FBM contributed to the conceptualization of this study. JLM, FBM, MTB, DSM and JN participated in protocol development and implementation of the study. NMS implemented transcriptomics RNA-Seq experiments. JLM wrote the manuscript and directed the statistical analyses implemented and verified by MFB. JLM, FBM, and MFB contributed to the overall analysis and/or biological interpretation. All authors participated in the revisions of the manuscript and/or interpretation of the results. All authors gave their final approval and agreed to be accountable for all aspects of the work.

Consent

This study has been approved by the Institutional Review Board of Carolinas Medical Center-Atrium Health, Charlotte, NC. All patients participating in this study have signed an informed consent (IRB file #01–16-09E).

Acknowledgments

We would like to thank CDx Diagnostics (Suffern, NY) for providing brush biopsy test kits for research purposes. We also thank Dr Alex Mira (Department of Health and Genomics, Center for Advanced Research in Public Health, FISABIO Foundation, Valencia, Spain) and Dr Alphonso Benittez-Paez (Centro de Investigación Príncipe Felipe, Valencia, Spain) for providing detailed protocols and critical advice to conduct microbiome transcriptomics of the dental plaque samples. We thank Alexis Kokaras for the next-generation sequencing and ProbeSeq program, Molecular Biology Core lab for RNA sequencing and Kathleen Sullivan for her editorial support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

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

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/20002297.2023.2246279.

Additional information

Funding

Funding was obtained from the Atrium Health Foundation.