Obesity influences composition of salivary and fecal microbiota and impacts the interactions between bacterial taxa

Abstract

Obesity is an increasing global health concern and is associated with a broad range of morbidities. The gut microbiota are increasingly recognized as important contributors to obesity and cardiometabolic health. This study aimed to characterize oral and gut microbial communities, and evaluate host: microbiota interactions between clinical obesity classifications. We performed 16S rRNA sequencing on fecal and salivary samples, global metabolomics profiling on plasma and stool samples, and dietary profiling in 135 healthy individuals. We grouped individuals by obesity status, based on body mass index (BMI), including lean (BMI 18–124.9), overweight (BMI 25–29.9), or obese (BMI ≥30). We analyzed differences in microbiome composition, community inter-relationships, and predicted microbial function by obesity status. We found that salivary bacterial communities of lean and obese individuals were compositionally and phylogenetically distinct. An increase in obesity status was positively associated with strong correlations between bacterial taxa, particularly with bacterial groups implicated in metabolic disorders including Fretibacterium, and Tannerella. Consumption of sweeteners, especially xylitol, significantly influenced compositional and phylogenetic diversities of salivary and fecal bacterial communities. In addition, obesity groups exhibited differences in predicted bacterial metabolic activity, which was correlated with host’s metabolite concentrations. Overall, obesity was associated with distinct changes in bacterial community dynamics, particularly in saliva. Consideration of microbiome community structure and inclusion of salivary samples may improve our ability to understand pathways linking microbiota to obesity and cardiometabolic disease.

INTRODUCTION

Obesity is a growing worldwide epidemic and is linked to a range of health issues including hypertension, type 2 diabetes, asthma, coronary heart disease, Alzheimer’s disease, and cancer (Alford et al., 2018; Avgerinos et al., 2019; Seganfredo et al., 2017; Thompson et al., 2007; Wahba & Mak, 2007). Known risk factors include imbalances between calorie intake and expenditure, genetics, stress, and disruptions in the endocrine system (Han & Lean, 2016; Seganfredo et al., 2017); however much remains unknown. Better characterization of mechanisms predisposing to obesity could enable novel prevention and treatment strategies.

The composition of an individual’s microbiota is increasingly being recognized as a contributor to obesity risk (Benahmed et al., 2021; Crovesy et al., 2020; Stanislawski et al., 2019). Microbiota can influence the host’s metabolic phenotype both by directly affecting energy and nutrient availability (Jumbo-Lucioni et al., 2010; Kaoutari et al., 2013; LeBlanc et al., 2013; Shortt et al., 2018; Turnbaugh et al., 2006), and through modulation of signaling pathways (Bindels et al., 2013; Davison et al., 2017; Fellows et al., 2018; Kimura et al., 2013; Mohammadkhah et al., 2018; Overby & Ferguson, 2021; Tilg & Moschen, 2016; Ye et al., 2017). Previous studies suggested that the fecal symbiotic bacterial community of obese individuals is less diverse than that of lean individuals (Stanislawski et al., 2019; Tilg et al., 2009). In addition, the abundance of several bacterial taxa including Lactobacillus, Pervotella, Alistipes, Akkermansia, and others vary with obesity status (Benahmed et al., 2021; Crovesy et al., 2020). Salivary microbiota of lean and obese individuals also differ in diversity and composition (Araujo et al., 2020; Benahmed et al., 2021; Andrade et al., 2020; Raju et al., 2019; Si et al., 2017). The abundance of several salivary bacterial taxa including Campylobacter, Aggregatibacter, and Veillonella was reported to be positively associated with obesity (Balakrishnan et al., 2021; Schacher et al., 2007; Szafrański et al., 2015). Higher abundances of Bacteroidetes, Spirochaetes, and Firmicutes were observed in lean individuals (Benahmed et al., 2021; Janem et al., 2017; Sohail et al., 2019). However, data are contradictory, even for rather abundant bacteria taxa. For example, the abundance of intestinal Lactobacillus was reported to be both positively and negatively associated with obesity (Azad et al., 2018; Crovesy et al., 2020; Khalili et al., 2019; Million et al., 2012). While some variability might be explained by differences in diet, geographical location, or analytical methods, these discrepancies may be due in part to complex interactions between microbial community members, where the metabolic activity of individual bacterial taxa can vary based on the activity of other microbes in the community (Estrela et al., 2016; Morris, 2015; Thommes et al., 2019; Zomorrodi & Segrè, 2017). Consideration of interactions between members of microbiota might be essential to improve the identification of bacterial mechanisms underlying obesity.

We hypothesized that the presence of obesity, in the absence of known disease, would associate with differences in microbiome composition and function. We further hypothesized that community structure and bacterial inter-relationships would differ by obesity status. We evaluated the differences in compositional and phylogenetic diversity of salivary and fecal microbiota between obesity groups in a well-characterized sample of healthy individuals. We examined interactions between bacterial taxa based on the obesity status of the host, and showed that predicted bacterial metabolic activity varies between obesity groups and is correlated with intestinal and circulating metabolite concentrations.

MATERIALS AND METHODS Study population

We analyzed data from the ABO Glycoproteomics in Platelets and Endothelial Cells (ABO) Study (n = 135) as described previously (Bagheri et al., 2021; Ferguson et al., 2018; Tang et al., 2019). Demographic information is provided in Supplement Table 0. Briefly, healthy non-pregnant and non-lactating women and men were recruited to a cross-sectional study. Individuals were non-smokers, with no medication or supplement use (apart from oral contraceptives), and no clinical disease. Participants completed dietary profiling (validated 3-day food records, and DHQ II food frequency questionnaires [FFQ]), and provided stool, saliva, and blood samples. Height and weight were measured at the study visit. Individuals were classified based on body mass index (BMI, weight (kg)/height (m)-squared), including lean (BMI 18–24.9; fecal samples n = 76, saliva samples n = 49), overweight (BMI 25–29.9; fecal samples n = 34, saliva samples n = 19), or obese (BMI ≥30, fecal samples n = 25, saliva samples n = 16), to explore differences in composition and function of microbiota by obesity. All participants provided written informed consent. The study was approved by the Institutional Review Boards of the University of Pennsylvania and Vanderbilt University.

Sample profiling

As we have previously described, 16S rDNA sequencing of the bacterial V4 fragment was performed on Illumina MiSeq platform using 135 fecal and 85 saliva samples to identify bacterial community composition (Tang et al., 2019). Global metabolomics profiling of fecal and plasma samples, from a subset of individuals (n = 75) was performed at Metabolon (Metabolon Inc., Morrisville, NC, United States), as previously described (Tang et al., 2019).

Pre-analysis processing Sequences alignment and normalization

Pre-analysis processing of 16SrRNA reads was performed with R v4.0.2 (Team R. Core, 2019). Demultiplexed sequences were filtered, forward and reverse reads were merged, and resulted sequences were assigned to amplicon sequence variants (ASVs), with the default settings of DADA2 pipeline v1.18.0 (Callahan et al., 2016). Chimeric sequences were also removed with the dada2 package v1.18.0 (Callahan et al., 2016). Sequence variants were assigned taxonomy with dada2 and SILVA v138.1 database (Callahan et al., 2016; Quast et al., 2012). ASVs counts were normalized with cumulative sum scaling method implemented in the metagenomeSeq v1.32.0 package (Paulson et al., 2013). In the salivary samples, we identified 1932 ASVs that belonged to 12 phyla, 19 classes, 44 orders, 70 families, 134 genera, and 229 bacterial species. In our fecal samples, we identified 5000 ASVs that belonged to 16 phyla, 26 classes, 55 orders, 86 families, 270 genera, and 338 bacterial species.

Alpha diversity

Normalized ASVs counts were used to calculate species richness, Shannon, and Gini–Simpson alpha diversity indices with the vegan v2.5.7 package (Oksanen et al., 2009).

Beta diversity

Bray-Curtis distances were calculated with vegan v2.5.7 (Oksanen et al., 2009). An unrooted neighbor-joining tree was computed with the ape package v5.5 (Paradis et al., 2004). The tree was optimized based on a generalized time-reversible model implemented in the phangorn v2.5.5 package (Schliep 2019; Schliep, 2011). Lastly, weighted and unweighted Unifrac distances between each sample were calculated with the phyloseq v1.30.0 package (McMurdie & Holmes, 2013).

Functional potential

Functional potential of the bacterial communities was predicted with PICRUSt2 according to the default pipeline (Douglas et al., 2020). Predictions were made for Enzyme Commission numbers (EC), Kyoto Encyclopedia of Genes and Genomes orthologs (KO), and MetaCyc pathways (Bairoch, 2000; Douglas et al., 2020; Kanehisa, 2000; Karp et al., 2002). In accordance with PICRUSt2 authors’ recommendations, the resulting data were transformed with the centered-log ratio transformation implemented in the ALDEx2 v1.24.0 package (Gloor, 2015).

Statistical analysis

Statistical analysis and data visualization was done with R v3.6.1 (Team R. Core, 2019). Beta diversity distances between obesity groups were compared with pairwise permutational multivariate analysis of variance, based on the vegan package v2.5.7 (Oksanen et al., 2009). The difference in alpha diversity measurements was evaluated with Wilcoxon signed-rank test, implemented in the rstatix v0.7.0 package (Kassambara, 2021). In order to evaluate if the obesity groups can be classified based on differential abundance (taxonomic units with the differential abundance of less than 20 amplicons in the whole data set, were filtered out, in order to avoid constant variables across the groups) of bacterial taxa and inferred functional abundances (based on EC, KO, and MetaCyc classification), we used linear discriminant analysis, implemented in the MASS package v7.3-51.4 (Venables & Ripley, 2002). In addition, we repeated a linear discriminant analysis using only the 15 most abundant bacterial taxa, in order to evaluate if the dominant bacterial taxa were sufficient for discrimination of the communities, with the obesity status. The results were visualized by plotting the first and second linear discriminants, with the ggplot2 v3.2.1 and the ggpubr v0.4.0 packages (Kassambara & Kassambara, 2020; Wickham, 2016). The difference in differential abundances of bacterial taxa and predicted ECs, KOs, and MetaCyc pathways, between obesity groups was evaluated with a pairwise t-test function, implemented in R v3.6.1 (Team R. Core, 2019). The correlations between differential abundances of bacterial taxa were calculated with Spearman’s rank correlation test, included in the Hmisc v4.5.0 package (Harrell & Harrell, 2019). Resulted correlation matrices were used to construct network plots, using the corrr v0.4.3 package (Kuhn & Jackson, 2020). In addition, the absolute values of correlation coefficients were compared between obesity groups with a pairwise Wilcoxon signed-rank test, implemented in the rstatix v.7.0 package (Kassambara, 2021). The influence of 133 recently consumed (from 3-day food records) and 185 habitually consumed (from FFQ) nutrients on beta diversity distances was evaluated with permutational multivariate analysis of variance using a quadratic model (Oksanen et al., 2009). The quadratic model was used as most living organisms, including bacteria have an optimal range of environmental conditions rather than a linear relationship (Bombin & Reed, 2016; Kindt & Coe, 2005; Leboffe & Pierce,). The difference in nutritional profiles between obesity groups was evaluated with the adonis function on Euclidean, Bray-Curtis, and non-binary Jaccard distances (Oksanen et al., 2009).

For enrichment analysis, we calculated the mean abundance of each KEGG ortholog for obesity groups and used them as input for MicrobiomeAnalyst (2021-07-01) shotgun data profiling tool, with the default settings (Dhariwal et al., 2017). False discovery rate (FDR) p-values were adjusted using the Benjamini–Hochberg correction, implemented in rstatix v0.7.0 package (Kassambara, 2021). We note that usage of any particular FDR threshold is ambiguous and often varies between microbiome studies; weaker correlations that fail to hold up to p adjustment methods often have biological relevance. Premature rejection of associations falling below conservative p-value thresholds may lead to loss of biologically meaningful data. (Althouse & Soman, 2017; Bombin et al., 2020; Bruce-Keller et al., 2015; Jehrke et al., 2018; Pawitan et al., 2005; Wu et al., 2015). For this reason, statistical results below 0.05 p-value threshold were considered to be significant. However, taking into account the difference in opinions and for the readers’ convenience, we report both unadjusted and FDR-adjusted p-values in supplementary data.

RESULTS Lean, overweight, and obese individuals can be separated into distinct groups based on their oral and intestinal microbiota

Evaluating beta diversity distances, we observed that salivary microbiota communities of obese and lean individuals were significantly different as measured with Bray-Curtis and Weighted Unifrac distances (Supplement Table 1). Based on linear discriminant analysis (non-overlapping confidence ellipses), obesity classes were separated by the differential abundances of bacterial ASVs (Figure 1a). Obesity groups were also clearly characterized based on the differential abundance of microbial species, genera, families, and orders but weaker based on classes and phyla (Supplement Figure 1).

FIGURE 1. Obesity groups can be discriminated by the abundance of salivary or fecal microbiota. Linear discriminant analysis of (a) ASVs identified in salivary samples (b) ASVs identified in fecal samples. ASVs with abundance of less than 20 sequences were filtered out. Obesity groups are represented by color, lean group by red, overweight group by green, and obese group by blue. Confidence ellipses are shaded. Normal data ellipses are unfilled and leveled to include 50% of the samples

In fecal samples, we did not observe a significant difference in beta diversity distances between any of the obesity groups (Supplement Table 1). However, based on a linear discriminant analysis, obesity groups could be classified based on the differential abundance of bacterial ASVs (Figure 1b). Obesity groups were also clearly characterized based on the differential abundance of bacterial species, genera, families, and orders but weaker at class and phylum ranks (Supplemental Figure 2). We did not observe any significant differences in alpha diversity indices between obesity groups in saliva or feces (Supplement Table 2).

Obesity status influences the differential abundance of individual bacterial taxa In saliva

In saliva, we observed that abundances of Campylobacterota, Firmicutes, and Spirochaetota were significantly different between obesity groups at the phylum rank. Obesity groups were significantly different in the differential abundances of 5 bacterial classes, 10 orders, 17 families, 33 genera, 52 species, and 409 individual ASVs (Supplement Table 3A). Across all taxonomic ranks, obese and lean individuals had the highest number of taxa that were significantly different in their differential abundances (Supplement Table 3A). We evaluated which of the 15 most abundant bacterial taxa were the most influential for defining each of the obesity groups with a linear discriminant analysis. At the genera taxonomic rank, Campylobacter, Veillonella, Aggregatibacter, and Prevotella defined the obese group (Figure 2). Although lean and overweight groups were not distinct from each other, Actinomyces and Haemophilus were characteristic for the overweight group (Figure 2). Overall, we note that across all taxonomic ranks the 15 most abundant bacteria taxa contribute only modestly to discrimination of obesity groups (Supplemental Figure 3).

FIGURE 2. Obese and lean groups can be characterized by the abundance of dominant bacteria genera. Linear discriminant analysis of the 15 most abundant bacterial genera identified in (a) Salivary samples (b) Fecal Samples. Obesity groups are represented by color, lean group by red, overweight group by green, and obese group by blue. The higher abundance of bacterial genera in the obesity groups is indicated by the direction of the vector rays. The intensity of vector rays’ color corresponds to the strength of the impact. Confidence ellipses are shaded. Normal data ellipses are unfilled and leveled to include 50% of the samples

In feces

In feces, at the phylum rank, only differential abundance of Fusobacteriota was significantly different between overweight and lean groups. Obesity groups were significantly different in the differential abundances of 2 bacterial classes, 8 orders, 10 families, 35 genera, 45 species, and 690 individual ASVs (Supplement Table 3B). The highest number of significant differences between groups varied with taxonomic rank but was always between lean and one of the overweight/obese groups. Linear discriminant analysis indicated that at the genus taxonomic rank Agathobacter and Parabacteroides were influential in discriminating obese from lean groups (Figure 2). Although lean and overweight groups were not clearly separated, lean group was primarily characterized by Blautia and Ruminococcus (Figure 2). Similar to what we observed in salivary samples, the most abundant fecal bacteria taxa were not the most influential variables for discriminating samples based on obesity status (Supplement Table 4).

The number of strong correlations between bacterial taxa varies by obesity status

We hypothesized that microbial community inter-relationships, as evidenced by correlations between taxa, would differ by obesity status. We assessed the number of strong correlations (> = |0.7|) between differential abundances of microbial taxa in saliva and stool samples by obesity group and found evidence for increasing inter-dependence in the setting of obesity (Figure 3). Among microbiota genera in saliva, there were 67 strong correlations in the obese group, 32 in the overweight, and only five strong correlations in the lean group. The absolute means of correlation coefficients were significantly different between all groups, and this observed pattern remained across all taxonomic ranks (Supplement Table 4). We observed a similar pattern in fecal samples, with 52 strong correlations between microbiota genera in the obese group, 20 in the overweight group, and only 8 in the lean group. The absolute values of the correlation coefficients, for differential abundances of the bacterial taxa were significantly different between all obesity groups. Obese individuals had more strong correlations between bacterial taxa than lean individuals across all phylogenetic ranks except phylum, at which no group had strong inter-bacterial correlations. (Supplement Table 4).

FIGURE 3. Number of strong connections between bacterial genera increases with the obesity status. Spearman’s rank correlation network between (a) Salivary bacterial genera of lean individuals; (b) Salivary bacterial genera of overweight individuals; (c) Salivary bacterial genera of obese individuals; (d) Fecal bacterial genera of lean individuals; (e) Fecal bacterial genera of overweight individuals; (f) Fecal bacterial genera of obese individuals. For (a–c) included genera had a minimum abundance of 30 sequences and for (d–f) minimum abundance of 20 sequences

Nutritional factors influencing bacterial communities

We examined recent (3-day food records) and habitual (food frequency questionnaire) dietary consumption by obesity status, and found that overall nutritional profiles were not significantly different between the obesity groups. We then examined the relationships between dietary variables and the overall bacterial community in all individuals, to identify influential nutrients from recent and habitual consumption. We applied Bray-Curtis, weighted Unifrac, and unweighted Unifrac distances, and assessed both linear and quadratic relationships. For recently-consumed nutrient, xylitol and pectins had significant linear relationships across all three methods, while inositol, glucose and omega-3 polyunsaturated fatty acids approached significance for quadratic relationships across all three methods (Supplement Table 5). For habitually-consumed nutrients, no nutrients displayed consistent linear relationships across all methods, while for quadratic relationships, sorbitol and pinitol, as well as dairy cheese and yogurt were consistently associated (Supplement Table 6). In the fecal bacterial community, recently-consumed pectins, folate, and fiber had consistent significant linear relationships, while oxalic acid, formononetin, biochanin A, and the ratio of polyunsaturated to saturated fat had consistent quadratic relationships (Supplement Table 5). For habitually-consumed foods, there were consistent linear relationships with cheese and vegetables, in addition to vegetable-derived nutrients (beta carotene, oxalic acid, Vitamin K). Significant quadratic relationships were observed for grains and processed meats, in addition to xylitol, caffeine, sodium, and potassium (Supplement Table 6).

Analysis of inferred metabolic pathways reveals enrichment in 2-oxocarboxylic acid metabolism in lean individuals in oral and intestinal microbiota

We hypothesized that the functional activity of microbiota, as predicted using PICRUSt2, would differ by obesity status. We assessed differences in inferred function between obesity groups, and found that obesity served as a good classifier for enzyme counts (ECs), KEGG orthologs (KOs), and MetaCyc pathways abundances in saliva (Figure 4). There were 969 significant differences in ECs, 3915 in KOs and 177 significant differences in the abundance of MetaCyc pathways across all groups (Supplement Table 7). In all cases, lean and obese individuals had the highest number of differences. 2-oxocarboxylic acid metabolism, terpenoid-quinone biosynthesis, and D-glutamine and D-glutamate metabolism KEGG pathways were enriched in lean individuals but not in the obese group (Supplement Table 8). The obese group was uniquely enriched in fluorobenzoate, sulfur, and several amino acid metabolic pathways.

Similarly, obesity groups could be characterized based on the abundance of MetaCyc pathways, KOs, and ECs in fecal samples (Figure 4). We observed 128 significant differences between the obesity groups in ECs, 391 in KOs, and 19 in MetaCyc pathways (Supplement Table 7), spread across lean, overweight, and obese groups. The lean group was uniquely enriched in 2-oxocarboxylic acid metabolism, D-glutamine and D-glutamate metabolism, and pentose and glucuronate interconversions, when compared with obese group. The obese group was enriched in C5-branched dibasic acid, lipoic acid, and one-carbon KEGG metabolic pathways (Supplement Table 8).

FIGURE 4. Obesity groups can be discriminated by metabolic potential predicted by PICRUSt2. Linear discriminant analysis of relative abundances of (a) ECs inferred from saliva samples (b) KOs inferred from saliva samples, (c) MetaCyc pathways inferred from saliva samples, (d) ECs inferred from fecal samples, (e) KOs inferred from fecal samples, (f) MetaCyc pathways inferred from fecal samples. Obesity groups are represented by color, lean group by red, overweight group by green, and obese group by blue. Confidence ellipses are shaded. Normal data ellipses are unfilled and leveled to include 50% of the samples

Abundance of inferred bacterial metabolic enzymes/pathways influences the host’s metabolites’ concentrations

We were interested in whether predicted functional activity would associate with measured metabolic activity, as assessed by metabolomic profiling of plasma and stool. We observed high numbers of correlations with predicted saliva microbial activity across all three databases (EC: 78,635 with plasma, 82,722 with stool; KO: 249,473 plasma, 263,616 stool; MetaCyc: 15,633 plasma, 17,915 stool). The highest number of correlations was observed with valerate and isoeugenol sulfate in plasma samples and with inosine in stool samples (Supplement Table 9). We similarly observed high numbers of correlations between predicted stool microbial activity and metabolites (EC: 92,852 with plasma, 109,830 with stool; KO: 299,557 plasma, 332,789 stool; MetaCyc: 18,179 plasma, 17,728 stool). The highest number of correlations was observed with 1-palmitoyl-GPE and CMPF in plasma samples and steviol in stool samples (Supplement Table 9).

DISCUSSION

Obesity has been linked to alterations in microbiota, however, the relative importance of gut and oral microbiota is unclear. We aimed to identify microbial signatures of obesity using both stool and salivary samples in healthy individuals classified as normal weight, overweight or obese based on their BMI. We observed that obesity status was associated with differences in bacterial community composition and shifts in inter-microbial relations that were especially evident in the salivary bacterial community. Although salivary and fecal microbiota were largely impacted by different nutrients, dietary sweeteners were associated with both composition and phylogenetic diversity of both the oral and gut bacterial communities. In addition, samples from obese and lean individuals were enriched in several unique metabolic pathways, inferred activity of which was correlated with plasma and stool metabolite concentrations.

Obesity influences microbial community composition, especially in saliva

In agreement with published research, we observed that oral bacterial community composition was distinct between lean and obese individuals (Araujo et al., 2020; Andrade et al., 2020; Raju et al., 2019; Si et al., 2017). In our work, we also observed that the difference in salivary bacterial composition between obese and lean individuals extends to phylogenetic diversity measurements. Consistent with previous research, we also observed some differences in gut bacterial communities between obese and lean groups, however in our work, the differences were not supported by Bray-Curtis or weighted Unifrac distances (Del Chierico et al., 2018; Palmas et al., 2021). There were no significant differences in overall dietary consumption, as assessed using both 3-day food records and food frequency questionnaires, suggesting that differences in microbiota were not attributable to differences in diet between lean, overweight or obese individuals. Our results suggest that at the level of the whole community, salivary microbiota composition better reflects the difference in obesity status than fecal microbiota.

With the analysis restricted to the dominant bacterial taxa, we observed a strong influence of Campylobacter, Aggregatibacter, Veillonella, and Prevotella on characterizing the obese group in salivary samples. Interestingly, all of these bacterial genera have been shown to be correlated not only with obesity but also with oral diseases, especially periodontitis (Balakrishnan et al., 2021; Durbán et al., 2013; Maciel et al., 2016; Mashima et al.,; Schacher et al., 2007; Szafrański et al., 2015). Considering the whole bacterial community (abundance >20 reads), we observed that some of the bacteria taxa with lower differential abundance had a stronger effect on the differentiation of the obese group than dominant bacteria, including Shuttleworthia at the genus rank and Mycoplasmataceae at the family rank that were also significantly more abundant in the obese group. Previous studies identified a correlation between Mycoplasmataceae and obesity (Huang et al., 2015; Kim et al., 2021). Although to the best of our knowledge, no previous works associated Shuttleworthia with obesity in humans, it was associated with obesity and elevated weight in model organisms (Henning et al., 2018; Lee et al., 2017; Xie et al., 2016). In addition, similar to what we observed with the dominant bacteria taxa, Shuttleworthia and Mycoplasmataceae are associated with periodontitis (Krishnan et al., 2017; Toyama et al., 2021).

In the fecal samples, the dominant bacterial genera that characterized the obese group were Agathobacter and Parabacteroides. Agathobacter and Parabacteroides were shown to be associated with metabolic disorders in humans and a murine model (Del Chierico et al., 2017; Liu et al., 2016; Salah et al., 2019; Schroeder et al., 2020). Similar to what we observed in the saliva samples, several less abundant bacterial taxa that were previously associated with obesity, including Mitsuokella and Neisseria, at the genus rank and Fusobacteriaceae and Gemellaceae, at the family rank, produced more impact on the separation of obese and lean categories than dominant bacterial taxa (Moreno-Indias et al., 2016; Palmas et al., 2021; Peters et al., 2018; Uberos et al., 2010; Zhang et al., 2021). Proportionally to all identified taxa, more organisms were significantly different in abundance between lean and obese groups in saliva samples, when compared with fecal samples, which might suggest that sampling oral microbiota may be more informative in identifying microbial biomarkers of obesity. Given the relative ease of collection of saliva as compared with stool, this could facilitate increased accessibility for research into the microbial contributors to obesity and cardiometabolic disease; however, this remains to be confirmed in independent studies.

Number of strong correlations between bacterial taxa increases with the obesity status

In saliva samples, bacterial taxa exhibited the highest inter-microbial connectivity (strong correlations ≥0.7) in obese individuals. In the obese group, the highest connectivity was observed for Fretibacterium (8 connections), F0058 (7 connections), Mycoplasma (7 connections), and Tannerella (7 connections). Several of these genera, including Fretibacterium, F0058, and Tannerella were shown to be correlated with metabolic disorders (Belstrøm, 2020; Haffajee & Socransky, 2009; Janem et al., 2017; Silva-Boghossian et al., 2018; Thomas et al., 2021). In addition, all of the most connected bacterial taxa were associated with periodontitis (Krishnan et al., 2017; Kwek et al., 1990; Nóvoa et al., 2020; Silva-Boghossian et al., 2018). In the lean group, the most connected bacteria exhibited less strong connections than in the obese group and were Atopobium (3 connections), Megasphaera (2 connections), and Prevotella 7 (2 connections). The abundance of Atopobium was shown to be reduced in obese individuals (Nardelli et al., 2020). Previous research indicated that the abundance of Megasphaera might increase after anti-obesity treatments (Federico et al., 2016; Kang et al., 2019). Prevotella was shown to be associated with plant rich diet and increase in abundance after antidiabetic treatment, however, the genus is very diverse (Ding et al., 2019; Jang et al., 2017; Precup & Vodnar, 2019).

In the fecal samples, the most connected bacterial genera identified in the obese group were Christensenellaceae R7 group (8 connections) and Ruminococcaceae UCG-005 (5 connections). Christensenellaceae R7 and Ruminococcaceae UCG-005 were shown to be associated with plasma lipoproteins and triglycerides (Vojinovic et al., 2019). Ruminococcaceae UCG-005 was also shown to be positively correlated with body weight and weight gain in a swine model (Gaukroger et al., 2020; Tang et al., 2020). In addition, several bacterial taxa previously implicated in metabolic disorders, including Actinomyces, Ruminiclostridium, and Lachnospiraceae exhibited strong inter-bacterial correlations in the obese but not in the lean group (Del Chierico et al., 2018; Lee et al., 2019; Liu et al., 2019; Zeng et al., 2016). The most connected genus in lean individuals was Ruminococcaceae NK4A214 (three connections). Previous research identified a negative correlation between Ruminococcaceae NK4A214 and high fat diet and hypertension (Calderón-Pérez et al., 2020; Yang et al., 2020). However, Christensenellaceae R-7 group and Ruminococcaceae UCG-005 were also among a few genera (total three) that had more than one strong correlation in lean individuals.

The impact of the higher degree of microbial interconnectivity observed in obese individuals is unclear but may represent a shift from the relative independence of bacterial taxa to a state more reliant on mutualistic relationships. Obesity is often associated with several physiological and environmental conditions that have the potential to act as stressors for the microbial community, including micronutrient deficiency, increased levels of reactive oxygen species, and an increase in c-reactive protein concentrations and inflammatory response in the host (Du Clos, 2000; McMurray et al., 2016; Via, 2012.; Yanoff et al., 2007). In accordance with the stress gradient hypothesis, several studies demonstrated that the presence of environmental stressors often increases positive facilitation between microbial taxa in the community (Hammarlund & Harcombe, 2019; Hernandez et al., 2021; Li et al., 2013; Lu et al., 2020). In addition, it was demonstrated that nutritional stress could increase the number of connections, in a co-occurrence network of the microbiota members (Ghosh et al., 2014). In agreement with these observations, we found that in the obese individuals, almost all of the strong inter-microbial correlations were positive.

Sweeteners and other nutrients influence compositional and phylogenetic diversity of salivary and fecal bacterial communities

We observed that recently and habitually consumed nutrients influenced bacterial communities across all individuals. For salivary samples, recently consumed nutrients influenced the bacterial community more than habitually consumed nutrients, for both compositional and phylogenetic beta diversity distances. Sugars and sugar alcohols, especially xylitol, mannitol, sorbitol, and pectin were especially influential factors impacting the bacterial community, based on compositional and phylogenetic diversity measurements. Interestingly, all of the listed compounds with the exception of pectin are used as sweeteners (Chattopadhyay et al., 2014; Pasha et al., 2002). Although the effect of sweeteners on gut microbiota was extensively shown in humans and animal models, the studies on oral bacteria community are limited (Gultekin et al., 2020; Söderling & Pienihäkkinen, 2020). Previous work has highlighted changes in the oral microbiota in response to consumption of dietary sweeteners (Štšepetova et al., 2019). To the best of our knowledge, this work is the first report specifically on the correlation between dietary sweeteners and phylogenetic diversity of the human’s salivary bacterial community.

Fecal microbiota community was consistently more influenced by habitual nutrient consumption than recently consumed nutrients, which might suggest a more stable microbial community. Similar to the saliva samples, consumption of xylitol and pectin influenced compositional and phylogenetic diversity of fecal microbiota. Consumption of sweeteners, including xylitol was reported to influence intestinal bacterial community composition (Gultekin et al., 2020). Pectin consumption was also shown to be correlated with compositional changes in the intestinal microbiota (Jiang et al., 2016; Larsen et al., 2019). In our study, compositional and phylogenetic measurements of the fecal microbiota were also consistently influenced by the consumption of vegetables and plant-derived compounds including fiber, oxalic acid, formononetin, and daidzein. Consumption of fiber, formononetin, and daidzein were shown to have microbiota-mediated beneficial effects on host’s metabolic health (Carrera-Quintanar et al., 2018; Makki et al., 2018).

Bacterial communities of obesity groups are associated with enrichment in predicted metabolic pathways, which are correlated with host’s metabolite concentrations

In both saliva and fecal samples, the microbiota of the lean individuals were enriched in 2-oxocarboxylic acid metabolism and D-glutamine and D-glutamate metabolism, based on functional prediction. 2-Oxocarboxylic acid metabolism is involved in ornithine and lysine biosynthesis, supplementation of which were shown to have a potential for improving metabolic health (Kalogeropoulou et al., 2009; Kanehisa, 2019; Park et al., 2012). D-Glutamine concentrations were shown to be decreased in obese individuals and glutamine supplementation may alleviate obesity symptoms (Abboud et al., 2019; Ren et al., 2019). Metabolic pathways enriched in the microbiota of obese individuals included one-carbon metabolism, which was previously shown to contribute to the development of obesity (Arnoriaga-Rodríguez et al., 2021). In addition, steatosis was shown to be associated with one carbon metabolism’s gene expression (Christensen et al., 2010). Enrichment in other pathways such as lipoic acid metabolism and degradation of valine, leucine, and isoleucine might be a response to increase in oxidative stress and branched-chain amino acids concentrations, often associated with obesity (Allam-Ndoul et al., 2015; Furukawa et al., 2017; Rochette et al., 2013).

Multiple host’s metabolites were significantly correlated with the abundance of KOs involved in enriched pathways. For example, the abundance of KOs, predicted in salivary samples and involved in 2-oxocarboxylic acid metabolism influenced the concentration of 435 plasma and 326 stool metabolites. Alpha-ketobutyrate was shown to be a biomarker of insulin resistance and glucose intolerance and in our study exhibited a negative correlation with more than half of the 2-oxocarboxylic acid metabolism pathway’s KOs, predicted from saliva samples (Gall et al., 2010; Syed Ikmal et al., 2013). In addition, KOs involved in 2-oxocarboxylic acid metabolism were correlated with adenosine and steviol in stool samples, both of which were shown to be beneficial for patients with metabolic disorders (Panagiotou et al., 2018; Pardo et al., 2017).

Our study had considerable strengths, including the availability of salivary and fecal microbial profiling, in addition to metabolic phenotyping. There were also some limitations inherent in all microbiome projects that are based on 16S rRNA sequencing. PCR reaction with degenerate primers is often used for sequence amplification (Kumar et al., 2011). For the taxonomic identification of bacterial samples and their phylogenetic analysis, primers for nine hypervariable 16S rRNA regions V1-V9 are broadly used (Kumar et al., 2011). However, different 16S regions are differentially conserved between bacterial groups and therefore might better suit the identification of particular bacterial taxa, which may lead to discrepancies between taxa identifications and diversity estimations (Chakravorty et al., 2007; Poretsky et al., 2014). In addition, due to the limited amplicon size, 16S microbiome analysis might have a lower confidence in determining deep taxonomic levels such as genus and species (Chakravorty et al., 2007; Poretsky et al., 2014). In addition, different sequence filtering methods, reference database for taxonomic identification, and even normalization methods are all known to cause a degree of bias between studies. Within our work we took several precautions to minimize the limitations described above. In order to identify sequences that belong to different bacterial taxa, instead of using 97% sequence similarity operational taxonomic units identification (OTUs), we used amplicon sequence variants (ASVs), also known as exact sequence variants (ESVs) that are based on exact sequence matches. This methodology allows for more accurate taxonomic identifications even with short sequences (Callahan et al., 2017; Glassman & Martiny, 2018; Prodan et al., 2020). Furthermore, knowing of the limitations in taxonomic identification resolution of 16S rRNA sequencing, we have not limited ourselves to analyzing the differences in microbiota communities at any given taxonomic rank. Most of the analyses were repeated at all taxonomic ranks and main conclusions are based on the patterns that are observed at multiple taxonomic ranks including family and higher. In addition, understanding that different 16S rRNA regions might be more suitable for identification of particular bacterial taxa, we refrained from the direct comparison of fecal and salivary samples, restricting ourselves to the comparisons of changes in bacterial communities based on the obesity status of the host, separately for oral and fecal microbiota.

In addition, results presented in this study are largely based on differential abundances of the identified microbial taxa and therefore might not be interpreted as causative. Therefore, future studies would be necessary to demonstrate the directions of interactions between the host and its oral and intestinal microbiota.

CONCLUSIONS

In this study, we identified differences in salivary and fecal symbiotic bacterial communities based on obesity status, in a population of otherwise healthy individuals. Our results suggest that inter-correlations between bacterial taxa are altered in the setting of obesity and suggest distinct differences in community dynamics at increasing levels of obesity. Consideration of microbial community correlation structure might be more informative than measurement of relative abundances of bacteria taxa or diversity measurements alone. In addition, across multiple comparisons, salivary microbiota provided a more distinct pattern of differentiation between obese and lean individuals, than fecal microbiota. Previous studies have primarily focused on the analysis of gut microbiota in obesity, however, our data suggest that sampling oral microbiota might provide further insights and patterns that are easier to detect across diverse groups.

ACKNOWLEDGMENTS

This work was supported by funding from NIH R01 HL142856, and the Layton Family Fund.

CONFLICT OF INTEREST

The authors declare no conflicts of interest.

AUTHOR CONTRIBUTIONS

Andrei Bombin and Jane F. Ferguson designed the research; Jane F. Ferguson acquired funding and generated the data; Andrei Bombin analyzed data and drafted the manuscript; Shun Yan, Sergei Bombin, and Jonathan D. Mosley advised on analytical methods; Andrei Bombin, Shun Yan, Sergei Bombin, Jonathan D. Mosley, and Jane F. Ferguson reviewed and edited the manuscript.