We evaluated individuals from local communities in the United Arab Emirates who smoked cigarettes for over 5 years and non-smokers that were otherwise healthy, excluding those with chronic respiratory illness, reported antibiotic use, and prescribed probiotic use for the past 3 months.We have also assessed the behavioral activation for depression scale (BADS), as previously described (17). Briefly, this scale contains 9 items with scores ranging from 0 to 54, with high total scores representing higher activation while a high score on the avoidance subscale means more avoidance. We also collected buccal swabs from all participants using Isohelix DNA/RNA Buccal Swabs (Isohelix Ltd. Harrietsham, United Kingdom) following the manufacturer's instruction (Isohelix Ltd.). Participants did not eat, drink, smoke, clean their teeth, or chew gum for the past 1–3 h prior to collecting the samples. Furthermore, participants were asked to rinse their mouths with water just before sample collection. Samples where then transferrdinto a sterile container, stored immediately in liquid nitrogen, and then transferred to a – 80°C freezer. All participants in the study signed written informed consent, and the Research Ethics Committee at the University of Sharjah approved the study (REC-18-02-18-01).

DNA was extracted using the Qiagen MagAttract PowerSoil DNA KF Kit (Formerly MO Bio PowerSoil DNA Kit) using a KingFisher robot. DNA quality was evaluated visually via gel electrophoresis and quantified using a Qubit 3.0 fluorometer (Thermo-Fischer, Waltham, MA, USA). Libraries were prepared with the Illumina Nextera library preparation kit using an in-house protocol (Illumina, San Diego, CA, USA). Paired-end sequencing using NextSeq 500 followed by shotgun metagenomic with the Sunbeam pipeline as described before (10). Raw sequencing reads were submitted to ENA and are available under accession number PRJEB53577.

All analyses were performed in R (v4.1.0). If necessary, p-values were corrected for multiple testing using Benjamini-Hochberg correction and are presented as q-values (18). For non-parametric testing of two groups, Wilcoxon tests were used and for contingency tables, X²-tests were applied with 1,000 replicates to compute p-values using Monte Carlo simulation. For analysis and visualization, the following additional R packages were applied: tidyverse (v1.3.1) for data handling; ggplot2 (v3.3.6), ggpubr (v0.4.0), and patchwork (v1.1.1) for visual presentation of the results. Additional R packages and further information about data analysis are given below in the subsections.

Already processed sequencing data from Pwas imported as phyloseq object [phyloseq v1.36.0 (19)] in R v4.1.0 and operational taxonomic units occurring in less the 20 individuals with just 1 contig were removed (2,292 OTUs left after filtering) (10). For further analysis, data was summarized at genus level using the phyloseq::tax_glom function.

Alpha diversity on genus level was estimated on non-normalized counts using DivNet v0.3.7 (sample-wise Shannon diversity; OTU6740 selected as base taxon) to account for observed and unobserved taxa. Differences in alpha diversity between groups were calculated using the betta_random function (breakaway v4.7.2) (20, 21) with BADS group (activation or avoidance) as corresponding fixed effects and gender as random effect to adjust for gender differences. To assess beta diversity on genus level, data was centered log-ratio (clr) transformed and distances were calculated using Euclidean distance (Aitchison distance). Permutational multivariate analysis of variance using distance matrices (PERMANOVA) was used to analyze differences in beta diversity (adonis function, as implemented in the vegan package v2.5-7, with 99,999 permutations and BADS group as corresponding fixed effect).

Differentially abundant taxa between BADS groups (model: H ~ BADS_group) were identified using Analysis of Compositions of Microbiomes with Bias Correction (ANCOM-BC v1.6.0; zero_cut = 0.8, tol = 10x10⁻⁶, max_iter = 10,000, struct_zero = TRUE, Benjamini-Hochberg correction to correct for multiple testing) (22). Beta coefficients (effect sizes) and standard errors were obtained from the ANCOM-BC log-linear model and were used for visualization. To verify ANCOM-BC results and account for potential gender effects, a count regression for correlated observations with the Beta-binomial approach implemented in corncob (v0.2.0) was applied with gender in the NULL model (H₀ ~ gender) to account for gender effects (23). The full model comprised gender and BADS group to test for differential abundant taxa (H₁ ~ BADS_group + gender). Resulting p-values were corrected for multiple testing using Benjamini-Hochberg correction.

Network reconstruction on clr-normalized genus abundances was performed using the SparCC approach (24), as implemented in the R package NetCoMi (v1.0.2). For this analysis, the 50 genera with the highest variance between groups were selected. Significant edges were selected using Student's t-test and community structures were estimated using greedy optimization of modularity. Hub node detection was performed using a threshold of 0.8, and a quantitative assessment of the network was performed using a permutation approach (100,000 bootstraps) with an adaptive Benjamini–Hochberg correction to adjust p-values for multiple testing.

Previously inferred functional profiles were imported into R as phyloseq object (10). Alpha diversity was estimated on subsampled data (phyloseq::rarefy_even_depth function with seed set to 1 and without replacement; 3,000 counts per sample) using Shannon index (diversity function; vegan package v2.6-2). Beta diversity on the full data was estimated as described above for the sequencing data, and differential abundant functional profiles were identified using ANCOM-BC (described above).

This study assessed buccal swabs from 105 healthy individuals sampled from the local communities in correlation with their cigarette smoking status, among other criteria such as BADS. We further divided participants based on BADS into avoidance and activation groups. We calculated 53 with more avoidance scores and 52 with higher activation, all matched for age, gender, and BMI (Table 1).

For sequencing analysis of the oral microbiome, we aggregated taxa abundances at phylum and genus level and plotted the relative abundances of all phyla and the most abundant genera (abundance >5%, Supplementary Figure 1).

Bacteroidota, Fusobacteriota, and Firmicutes, Proteobacteria, and Actinobacteria were found to be the most abundant phyla, and the most abundant genera were Streptococcus, Haemophilus D, and Veillonella. Next, we estimated the alpha diversity of the oral microbiome between BADS avoidance and activation groups using Shannon's diversity index. The BADS avoidance risk group showed significant higher diversity (betta test with gender as a random variable to correct for gender effects; p < 0.001, estimate = 0.2935 ± 0.0850; Figure 1A). Community composition at genus level (beta diversity estimated using Aitchison distance) was not significantly different between BADS groups (PERMANOVA; p = 0.0664, R² = 0.0120; 99,999 permutations; Figure 1B).

To further explore compositional differences between both groups, we conducted a differential abundance analysis using a compositional data approach with bias-correction (ANCOM-BC) to identify taxa with significantly different abundances between the avoidance and the activation group (Figure 2A).

The phyla Bacteroidota (effect size 0.5047, q = 0.0037), Campylobacterota (effect size 0.4012, q = 0.0276), Firmicutes A (effect size 0.3646, q = 0.0128), Firmicutes I (effect size 0.3581, q = 0.0268), and Fusobacteriota (effect size 0.6055, q = 0.0018) were found to be significantly increased in the avoidance risk group, but Verrucomicrobiota (effect size −0.6544, q = 0.0401), was found to be significantly decreased in the avoidance risk group (Figure 2B), whereas no genera were found to be significantly different between the groups. We applied a second approach using count-regression for correlated observations (corncob) to validate our findings further to determine the differentially abundant taxa between BADS avoidance and activation groups after correction for gender effects (Supplementary Figure 2).

Besides the phyla identified in Figure 2, eight additional phyla were found to be significantly between the two groups (q < 0.01; Supplementary Figure 2A). On genus level, 40 genera showed significant differences between the two groups (q < 0.01; Supplementary Figure 2B). Next, we sought to explore the difference between both conditions using differential network analysis. Network analysis of the top 50 genera displayed the highest variation between the two conditions (Figure 3A).

Network properties are not significantly different between both conditions in terms of degree, betweenness centrality, and closeness centrality (p > 0.05). However, we found that each network has specific hub nodes. Campylobacter B, Centipeda, and Veillonella were identified as hub nodes in the avoidance group, whereas Haemophilus and Streptococcus were identified as hub nodes in the activation group (Figure 3B).

Functional profiling between both groups using shotgun metagenomics shows no significant differences with respect to alpha diversity (Shannon index; Mann-Whitney U test p = 0.1029) or beta diversity (Aitchison distance; PERMANOVA, p = 0.3862, R² = 0.0076) (Figures 4A,B, respectively).

Average functional profile abundance for profiles with an abundance above 1% (19 profiles) is shown in Figure 4C. Lysine degradations, which are not among the top 19 profiles, showed a significantly different abundance between the two conditions (ANCOM-BC, q = 0.0692) (Figure 4D). Further, we did not detect a significant difference in GABA abundance between both groups (Supplementary Figure 3).