PCR amplification

We amplified the eluted bacterial DNA using a previously described set of degenerated primers [27] which cover the V1-region of the 16S rRNA gene (F: 5′-AGAGTTTGATCMTGGCTCAG-3′ and R: 5′-TTACTCACCCGTICGCCRCT-3′), generating ~113 bp amplicons, compatible with the sequencing platform. For a better representation of the bacterial DNA, we performed a total of six simultaneous PCR amplifications for each subject. PCR mixes consisted of 2 μM of each primer, 5 μl of the 2X high-fidelity Taq DNA polymerase master-mix (Roche, Mannheim, Germany) and 2 μl of eluted DNA in a final volume of 10 μl. The cycling conditions included an initial denaturation step (94°C for 2 mins) followed by 35 cycles at 94°C for 15 secs and 60°C for 30 secs and a final extension step at 72°C for 5 mins.

Sequencing

We evaluated 2 μl of the amplification products in 8% silver-stained polyacrilamyde gels [28] and quantified the amplicons of each individual using a Qubit high-sensitivity double stranded DNA kit (Qiagen, Crawley, UK). The amplicons were purified using AMPure XP beads (Beckman Coulter, Indianapolis, USA), end-repaired and Ion Torrent adaptors with 10 bp barcodes were ligated to their ends using proprietary kits (Life Technologies, Carlsbad, USA). For emulsion PCR we pooled together equimolar amounts of amplicons from each sample, as indicated by a Bioanalyzer High Sensitivity chip (Agilent Technologies, Santa Clara, USA). All samples were sequenced on the Ion Torrent PGM platform (Life Technologies).

Sequence filtering

Sequences processed by the latest version of the Ion Torrent server (version 3.6.2) were used as input into the Qiime (Quantitative insights into microbial ecology) software package (Version 1.6.0) [29]. We first removed sequences that had an average quality score < 20, then identified barcodes used for subject-assignment accepting no more than 2 mismatches and discarded sequences smaller than 30 nt or greater than 200 nt. PCR primers identified at the start or at the end of the reads, allowing a maximum of 1 mismatch, were trimmed. Sequences with no identifiable primers were discarded. After primer trimming we again removed all sequences under 30 nt. We used the remaining sequences for downstream analysis.

Sequence clustering and OTU filtering

After primer- and quality-filtering, we clustered the remaining sequences with 97% identity using UPARSE (implemented in USEARCH v7) [30], and picked the less dissimilar sequence of each cluster as its’ representative. Chimeric sequences (and clusters) were identified using UCHIME [31] and the Broad Institute’s chimera slayer database (version microbiomeutil-r20110519) and excluded. We then used the Qiime interface to perform BLAST [32] analysis (using default parameters, with a maximum e-value of 0.001) against the formatted SILVA 111 SSU database [33] for taxonomic rank assignment of each sequence and, subsequently, to each operational taxonomic unit (OTU). In order to decrease PCR- and sequencing-derived errors, OTUs that were not present in at least 25% of all samples, had less than three reads, or had no taxonomic classification, were discarded from further analysis.

Data normalization

Alpha diversity

To investigate species diversity, richness and evenness, we used the same sequence sampling depth for all samples by rarefying the OTU table. Rarefaction curves, the Shannon-Weaver diversity index [35] and equitability for species evenness were calculated using Qiime. Observed species richness and richness estimators Chao1 [36] and ACE [37] were calculated using the R (Open sourced statistical software, Vienna, Austria - Version 3.0.1) Phyloseq package (Version 1.4.5) [38].

Beta diversity

The Euclidean distance matrix, the coordinates for the Principal Coordinate Analysis (PCoA), and the intra- and inter-dissimilarity indexes were calculated using Qiime.

Sequence filtering

After size, quality and primer trimming, a total of 3,843,174 sequences remained, with an average of 174,689 sequences/sample and a mean sequence length of 72 ± 4.9 nt (standard deviation) after primer trimming and barcode removal.

Sequence clustering and OTU filtering

When all individuals were considered, a total of 2,931 clusters were obtained. Twenty-one (0.7%) of these were removed after being identified as chimeras by UCHIME, and 112 (3.8%) had no assigned taxonomy. After filtering to remove OTU clusters with less than 3 sequences and not present in at least 25% of all samples, 1,849 OTUs remained.

Beta diversity

Using oral wash and broncheoalveolar lavage samples, a previous study found significant differences in bacterial communities when comparing smokers and non-smokers, suggesting that tobacco may significantly affect the oral microbiome [39]. In our study we first sought to evaluate if the 16S rRNA data generated would be capable of clustering individuals according to their use (or not) of alcohol/tobacco. We opted for the binary Euclidean distance metric due to the size of the 16S rRNA hypervariable region. Figure 1A shows the results of the PCoA plot using the distance matrix obtained by the binary Euclidean metric, demonstrating that 16S rRNA profiles clustered all individuals from the SD group together, as well as all samples from the S group. One of the control samples clustered closer to smokers/drinkers, whereas the other eight remaining control samples formed an independent group. Also, one of the smokers/drinkers samples clustered closer to only smokers, whilst the other six formed an independent group. The cladogram in Figure 1B indicates that all samples of the control group (except sample C01 which clustered closer to smokers/drinkers in the PCoA plot) belonged to the same main branch and that samples from smokers/drinkers were closer to only smokers than to controls. Also, all members of the only smokers-group clustered in the same tree branch. This consistent clustering suggests that alcohol and principally tobacco might influence the growth of certain bacterial species and may impact the normal abundance of bacteria commonly found in the oral cavity’s mucosa. It also suggests that smokers/drinkers have bacterial communities more similar to that of only smokers, rather than to controls.

Species richness

In order to measure species richness in the three groups we used the number of distinct OTUs observed in each sample when the same sequence sampling depth was applied at two different OTU-abundance thresholds (Figure 3). Saturation was not achieved when all OTUs represented by at least three reads were considered. This is likely to be derived from the combined effect of a still low coverage of very rare OTUs (those with frequencies below 0.03%, i.e. less than three reads/10,040 reads), together with PCR-amplification and/or sequencing artifacts. On the other hand, using this cutoff of at least three reads, which might underestimate very rare bacteria and precludes the representation of microorganisms non-amplifiable by this set of 'universal primers’, as well as bacteria that cannot be discriminated on the basis of sequences derived from the variable region 1 (V1) of the 16S rRNA gene, we clearly see that the microbiome of all groups exceeds 300 distinct OTUs (Figure 3A), indicating a rich community. Indeed, bacterial communities of the oral cavity are one of the most complex across all body sites, second only to the colon [41]. When we considered a more stringent OTU abundance threshold (Figure 3B), a saturation trend for OTUs with at least 100 reads (1% of the sequences used for each subject) could be seen and indicates the presence of at least 200 OTUs with frequencies above this threshold (Figure 3B). When considering OTUs with at least three reads, we observed that the average species richness of controls was about 35% and 17% larger than that of only smokers and smokers/drinkers, respectively (Figure 3A-B). When the same sequencing depth was used to evaluate species richness between pairs of the three groups, using three distinct metrics, we observed a significant decrease in species richness for only smokers when compared to controls. A significant reduction in OTU richness between controls and smokers/drinkers was also found at the sampling depth of 142,161 reads. The reduction in richness (in both SD and S when compared to C) was observed when using the Observed Richness (only smokers p-value= 0.0005; smokers/drinkers p-value= 0.0446) and also for richness estimators Chao1 (only smokers p-value= 0.0013; smokers/drinkers p-value= 0.0465) and ACE (only smokers p-value= 0.0009; smokers/drinkers p-value= 0.0474) (Figure 4). Surprisingly, we found no statistical difference in species richness between only smokers and smokers/drinkers using any of these metrics (p > 0.05).

Species evenness

Species diversity

Also, we evaluated species diversity using the Shannon-Weaver index and found no significant differences in diversity between each of the possible group pairs, but did see an increase in controls (mean C: 5.19, SD: 4.76 and S: 4.59). This suggests that the key difference in alpha diversity between controls, smokers/drinkers and only smokers lies solely in species richness. A previous study found a more diverse microbiome in the oropharynx of smokers than non-smokers [43], however it was restricted to tobacco-induced effects and no alcohol-drinking habits were reported. Another study evaluated the effects of chronic alcohol feeding on the intestinal microbiome of mice [44]. In this case, a reduction in species richness was observed along three time-points in alcohol-fed versus pellet-fed animals, and this decrease showed a negative correlation between species richness and the cumulative time mice were alcohol-fed. Another recent work [45] also found a reduction in species richness, this time in the oral biofilms covering the teeth of rats that received 30-days of a 20% ethanol-containing diet. Besides technical limitations due to the approach used for bacteria detection and quantification in the biofilms of rat teeth (restricted to probes for only 33 bacteria of the human oral microbiome) these authors were able to demonstrate a significant association between alcohol consumption and lower levels of bacterial richness. A reduction in microbial diversity has been associated with other diseases such as irritable bowel syndrome [46], eczema [47], Crohn’s disease [48], obesity [49] and recurrent Clostridium difficile- associated diarrhea [50], evidencing the crucial role that commensal bacteria play on keeping pathogenic microbes at bay and maintaining general homeostasis.

The chronic and simultaneous use of alcohol and tobacco appears to significantly decrease species richness in the human oral microbiome, which could lead to dysbiosis, the unbalancing of the relative abundance of individual components of the microbiota compared to their abundance in health [51]. It has been hypothesized that dysbiosis could lead to a variety of diseases [23]. Our data indicates that tobacco, together or not with alcohol, affects the microbiome leading to reduced species richness, through a mechanism that probably includes pro- and anti-inflammatory molecules [52]. To our knowledge, the work presented here is the first to evaluate microbial diversity and composition in biofilms of the oral mucosa in groups of human subjects with regards to the simultaneous and chronic use of alcohol and tobacco and the use of tobacco alone.

Phyla log abundances

The most abundant phyla found in all groups were (in decreasing order of abundance): Firmicutes, Actinobacteria, Proteobacteria, Bacteroidetes, TM7, Fusobacteria and Tenericutes (Figure 5A). At the phylum level, six out of 16 phyla were found to have significant differential log abundances between the three groups (Additional file 1: Table S1). Only smokers had an increase of Actinobacteria and Fusobacteria and a decrease of Firmicutes, Cyanobacteria, TM7 and BD1-5 when compared to the two other groups, whilst smokers/drinkers had a decrease of Fusobacteria. TM7 is a novel candidate bacterial division and previous studies have shown that microbes from this division are commonly found in the human oral flora at relatively low abundance, generally around 1%[53],[54] and up to 8% [55] of the oral bacterial population, and a possible role of this phylum in periodontal disease has been previously suggested [53],[56]. In only smokers, a significant reduction in log abundance of Proteobacteria could be observed (Kruskal-Wallis test, p=0.03), but after correcting for multiple taxa testing, this decrease was no longer significant (p=0.06). A reduction in the relative abundance of this phylum in broncheoalveolar-lavages of smokers when compared to non-smokers has been reported [39], pointing to a role that tobacco may play on reducing representatives of this phylum.

Genera log abundances. At the genus level, 47 out of the 279 (16.8%) identified genera were found to have significant differential log abundances between the three groups (Additional file 2: Table S2). The most abundant genera found in controls were (in decreasing order of abundance): Veillonell a, Actinomyces, Streptococcus, Haemophilus and Neisseria. In smokers/drinkers the most abundant genera were: Veillonell a, Actinomyces, Rothia, Granulicatella and Streptococcus and in only smokers were: Actinomyces, Rothia, Veillonella, Streptococcus and Leucobacter. All these genera have been reported before as common residents of the oral cavity and of oral biofilms [20],[24],[57]. We found a significant decrease in log abundance of Neisseria in both only smokers and smokers/drinkers when compared to controls (Kruskal-Wallis test p= 0.026; Figure 5B), but found no difference in the log abundance of this genus between SD and S. Whereas Neisseria is one of the most dominant genera found in the larynx, including laryngeal tumors [58], other studies have reported that the abundance of this genus is modulated by tobacco as its levels are reduced in the oral wash, oro-/nasopharyngeal swabs and bronchoscopic alveolar lavages of cigarette smokers [39],[43],[59].

Only smokers had significant increases in Prevotella and Capnocytophaga and decreases in Granulicatella, Staphylococcus, Peptostreptococcus and Gemella when compared to the other two groups. A significant decrease in abundance of Peptostreptococcus in smokers has been evidenced before [43], suggesting the susceptibility of this genus to smoke exposure. It is of interest to note that this particular reduction may be significant as several species belonging to this genus have shown to interfere in the growth of pathogenic bacteria in the upper respiratory tract [60]. Another genus that also seems to be modulated by smoking is Gemella, with a previous study also finding a decrease in the abundance of this genus [39]. In our analysis, the genus Porphyromonas, which has been shown to be increased in smokers [39],[43] and to have a role in periodontitis [61],[62], was also found to have higher abundances in only smokers (mean C: 0.67, SD: 0.29 and S: 1.8). However, after multiple taxa correction, the p-value was no longer significant (Kruskal-Wallis test p-value=0.08). Subsequent reports using 16S rRNA sequence profiling of subgingival plaque identified an increase in several disease-associated organisms in smokers, including Parvimonas, Fusobacterium, Campylobacter, Bacteroides, Dialister, and Treponema spp. and a decrease in potential health-promoting taxa from the Veillonella, Neisseria, Streptococcus, and Capnocytophaga genera [63]. Aggregibacter levels were higher in controls, whilst smokers/drinkers had lower log abundances of Fusobacteria compared to C and S (Kruskal-Wallis test; p-values= 0.04417 and 0.02541, respectively).