The domestic moist STPs were purchased from vendors in the Delhi and Uttar Pradesh state of India. The 16S rDNA sequencing techniques were used to analyze the bacterial diversity of 11 moist STPs from 4 different categories including Khaini (K1, K2, and K3), Moist-snuff (MS1, MS2, and MS3), Qiwam (Q1, Q2, and Q3), and Snus (S1 and S2) samples. The STPs were stored at −20°C to inhibit the further growth of microorganisms. The STPs were open in the sterilized conditions and microbial metagenomic DNA was isolated with a Power-soil DNA isolation kit as per the protocol provided by the manufacturer (Qiagen, Bangalore, India). The purity and quantity of isolated metagenomic DNA from the different STPs was checked by NanoDrop (Thermo, Bangalore, India). Further, the quality of metagenomic DNA was confirmed agarose gel (1%) before amplification by PCR. The metagenomic DNA concentration of all tested SLT products was found to be >30 ng/μl.

The PCR was performed to amplify the 16S rDNA V₃–V₄ region of metagenomic DNA isolated from the different moist STPs. The extracted metagenomic DNA (40 ng) was amplified with a pair of universal primers (10 pM of each) (FP: 5′-AGAGTTTGATGMTGGCTCAG-3′ and RP: 5′-TTACCGCGGCMGCSGGCAC-3′) as described earlier (Kroes et al., 1999). Along with primers and metagenomic DNA, a master mix was added containing dNTPs (0.5 mM), MgCl₂ (3.2 mM), high-fidelity DNA polymerase and PCR enzyme buffer. The PCR amplification was performed with the following condition: 95°C for 3 min chase by 25 cycles at 95°C/15 s, 60°C/15 s and 72°C for 120 s and with a final elongation at 72°C for 10 min. The amplified 16S PCR amplicons were purified and subjected to agarose gel (2%) and NanoDrop for quality check.

The Ampure beads (Beckman Coulter Inc., Indianapolis, IN, United States) were used to purify the amplicons of each sample by eliminating the unused primers and to prepare the sequencing libraries an additional 8 cycles of PCR was executed via Illumina barcoded adapters (Supplementary Table 1). Further, prepared libraries were purified by Ampure beads and quantified with the help of a QuDye-dsDNA HS assay kit (Thermo Fisher, Bangalore, India). The Illumina Miseq with a 2 × 300PE v3 sequencing kit was used for sequencing at Biokart India Pvt. Ltd., Bengaluru, India. The raw data was submitted to NCBI Short Read Archive (SRA) under BioProject accession number PRJNA767533.

The raw data binary base call (BCL) file received from the sequencer was de-multiplexed into FASTQ format. Quality control checks of raw sequenced data were performed by FastQC (Version 0.11.9) and MultiQC (Version 1.10.1) tools. Next, removal of contaminant adapters and trimming of low-quality reads were done by TrimGalore (Version 0.6.6)¹. The QC passed samples were again analyzed by Quantitative Insights Into Microbial Ecology (QIIME version 1.9.0) workflow including fusion of paired-end reads, chimera elimination, OTU clustering and taxonomy assignment (Caporaso et al., 2010). The QIIME workflow facilitates precise exploration at the genus level. The Kraken2 with database NCBI was used as a reference for OTU picking (Wood et al., 2019). Further, data was filtered on a web-based platform MicrobiomeAnalyst to eliminate the low quality or uninformative features using minimum count = 4, 20% prevalence and low variance filter (Inter-Quartile range – 10%) (Chong et al., 2020). A total of 148 low abundance features were separated based on prevalence and a total of 28 low variance features were removed based on inter-quartile range. The 247 number features remain after the data filtering step. All analyses like rarefaction curve, α-diversity, principal coordinate analysis (PCoA), core microbiome, cluster study, random forest, Linear discriminant analysis Effect Size (LEfSe) and Sparse Correlations for Compositional data (SparCC) were performed by MicrobiomeAnalyst². Subsequently, the metabolic pathway analysis was executed using the PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) algorithm (Langille et al., 2013). The functional genes associated with bacteriome of moist STPs were derived from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000). Next, Taxon Set Enrichment Analysis (TSEA) module of MicrobiomeAnalyst was performed to identify the biologically or ecologically meaningful patterns of STP associated bacteriome by analyzing them with context to pre-defined taxon set. To check the reproducibility of sequencing results, duplicates of Snus samples (S1_dup and S2_dup) were taken.

The relative abundance percentage of each OTUs of taxonomic classification was calculated and plotted using Origin software. For data set phyla, we removed OTUs < 10 reads and for the data set genus, we removed OTUs < 100 reads. Analysis of Variance (ANOVA) was used to determine differences between products. A p-value < 0.05 was considered statistically considerable.

After all pre-processing steps, an OTU table was generated and the number of amplicons obtained for each moist STPs varies between 10,601 and 131,082 (Supplementary Table 2 and Supplementary Figure 1). Most of the moist STPs showed an increase in observed bacterial species richness and diversity (Supplementary Figure 2). The numbers of OTUs were found to increase as the sequencing depth increases and saturation of species richness was observed in the rarefaction curves for all the tested moist STPs (Supplementary Figure 2). However, the numbers of observed bacterial species were less in K3 and S2 as compared to other moist STPs. Interestingly, the STPs having the highest moisture content such as Moist-snuff and Qiwam exhibited higher overall species diversity than Khaini and Snus which have comparatively low moisture content (Supplementary Figure 2 and Supplementary Table 3).

Further, OTUs richness was estimated by determining α-diversity (within-sample diversity) indices (Chao1, Fisher, Shannon, and Simpson) of STP-associated bacteriome. The Moist-snuff and Qiwan showed increased diversity indices as compared to Khaini and Snus suggesting that the Moist-snuff and Qiwam products have higher expected species richness of the bacteriome (Supplementary Figure 3 and Table 1). However, there was no statistically significant change observed among the four groups of moist STPs (Supplementary Figure 3). All the moist STPs displayed Good’s estimator values > 99% suggesting that the majority of bacterial species in the sample have been detected (Table 1).

The β–diversity indicates differences in the bacterial community profile between the samples. To calculate β–diversity among the moist STPs, the Bray–Curtis dissimilarity metric was determined from the OTU abundance and exploited in Principal Component Analysis (Goodrich et al., 2014). The Permutational Multivariate Analysis of Variance (PERMANOVA) algorithm on Bray–Curtis dissimilarity was applied to construct Principal Coordinate Analysis (PCoA) plots (Kelly et al., 2015). PERMANOVA analysis of Bray–Curtis dissimilarities revealed that the bacteriome of each group was highly dissimilar (PERMANOVA; F = 2.9354, R² = 0.49456, p < 0.004) (Supplementary Figure 4). The 3D-PCoA plot displayed that the Moist-snuff and Qiwam samples showed close association to each other and therefore, restrain more related bacteriome profiles (Figure 1A). Next, one sample from Khaini (K1) and one sample from Snus (S1) was also found to be close to each other and clustered together with Moist-snuff and Qiwam products and the sample K3 and S2 clustered separately from the other samples (Figure 1A). The interactive PCoA 3D plots at the level of the genus were constructed and dissimilar samples K3 and S2 showed a higher abundance of genera Actinobacteria and Prevotella, respectively (Figures 1B,C).

Further, a ‘Random Forest’ algorithm was applied to validate the similarity and dissimilarity in bacteriome among the four groups of moist STPs (Figure 1D). The decision tree constructed from the random forest classification recognized distinctive bacteriome in moist STPs. In the error plots identified from random forest analysis, the red line represents overall genera present in moist STPs while distinct genera present in Khaini, Moist-snuff, Qiwam, and Snus were indicated by yellow, green, blue, and magenta-line, respectively (Figure 1D). Among all three samples of Khaini, one sample contained unique genera and two samples demonstrated overlapping genera with Qiwam and Snus; whereas Moist-snuff and Qiwam products showed resemblance with each other. The Snus products exhibited a unique genera profile (Figure 1D).

Bacterial populations recognized in moist STPs were first analyzed at the phyla level (Figure 2A). There were 4 major phyla Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes (range 89–98%) were observed in all moist STPs. The other notable phyla were Acidobacteria, Chloroflexi, Cyanobacteria, Fusobacteria, Gemmatimonadetes, Spirochaetes, and Verrucomicrobia. The relative abundance of phylum Proteobacteria was found to be majorly present in all STPs except S2 (Figure 2A). The highest proportion of Proteobacteria was observed in the Khaini group (K1–91%, K2–57%, and K3–50%). Moist-snuff (MS1–41%, MS2–49%, and MS3–43%), Qiwam (Q-41% and Q2–50%), and Snus (S1–39%) also showed an increased level of Proteobacteria phyla compared to other phyla. In addition, phyla Firmicutes was the predominant phyla in S2 (60%) with the lowest presence in K3 (4%) (Figure 2A). The third most prevalent phylum identified in all STPs was Bacteroidetes (Figure 2A). However, Khaini had a lower abundance (range 2–14%) of Bacteroidetes as compared to Moist-snuff, Qiwam, and Snus (range 15–23%). The phyla Actinobacteria was noticeably present in all moist STPs ranging from 1 to 10% (Figure 2A).

A total of 549 genera were identified in all moist STPs and OTUs > 100 reads showed 61 genera were in relative abundance (Figure 2B). The genus Acinetobacter was abundantly observed in Khaini (90% of K3, 34% of K2, and 28% of K1). Another genus Prevotella was found to be significantly high in one of the Snus (36% of S2). However, a moderate level of Prevotella was observed in Moist-snuff (MS1-15%, 10% of MS1 and MS2), Qiwam (14% of Q1, 12% of Q3, and 10% of Q2) and Snus (9% of S1). The other important genera observed in STPs were Faecalibacterium (17% of S2), Bacillus (12% of Q2, 11% of K1 and MS2, and 9% of MS1), Lactobacillus (10% of K2), and Ruminococcus (7% of S2) (Figure 2B).

Despite inter-product variability, there was a core bacteriome identified in the moist STPs that remain unchanged in their composition across different groups of moist STPs. Core bacteriome investigation was executed at the genus level according to sample prevalence ≥ 20% and relative abundance ≥ of 0.2% (Figure 3 and Supplementary Table 4). The 22 core bacterial genera were identified and the prevalence of Acinetobacter, Bacillus, and Prevotella genera was observed in moist STPs (Figure 3). Additionally, Acetobacter, Lactobacillus, Paracoccus, Flavobacterium, and Bacteroides were found to be the dominant core bacteria in moist STPs (Figure 3).

The best correlation among samples of moist STPs at the OTU level was determined using the Bray–Curtis index (Supplementary Figure 5). The dendrogram showed similarities between the products and 11 moist STPs clustered into three groups. Group-I consists of two Khaini products K1 and K2 whereas K3 was clustered with Q2 and MS2 into a subgroup of Group-II. The high moisture containing products such as Moist-snuff and Qiwam clustered together into subgroups of Group-II (Supplementary Table 3 and Supplementary Figure 5). Therefore, the high moisture content may lead to similar bacterial diversity in moist STPs. Group-III comprises Snus samples S2 and S_2 dup whereas S1 clustered with their respective sample in Group-II which confirmed that the sequencing method produced reproducible data (Supplementary Figure 5).

Additionally, a hierarchical clustering heat map was generated for improved visualization of the distinct bacterial genera abundance across different moist STPs (Supplementary Figure 6). For Khaini samples the dominant genera included Acinetobacter, Staphylococcus, Panacibacter, Citrobacter, Pediococcus, Lactococcus, Weissella, Lactobacillus, and Leuconostoc. The Moist-snuff samples showed the abundance of genera such as Mycoplasma, Pantoea, Filifactor, Trueperella, Lysinibacillus, Marinobacter, Frankia, Campylobacter, Olsenella, Aeromonas, Dolosigranulum, Agromyces, Collinsella, and Fuerstia (Supplementary Figure 6). The Qiwam products demonstrated the dominance of Chelativorans, Proteus, Capnocytophaga, Microvirga, Methylobacterium, Turicibacter, Parvibaculum, and Simkania genera. The Snus products illustrated the occurrence of genera Tabrizicola, Sulfuritortus, Fimbrimonas, Anaeromyxobacter, Caloramator, Immudisolibacter, Mannheimia, Dermabacter, Ruminococcus, and Lachnoclostridium (Supplementary Figure 6). The genera having clinical relevance like Pseudomonas, Haemophilus, Actinomyces, Neisseria, Streptococcus, Campylobacter, Corynebacterium, Porphyromonas, and Fusobacterium were abundant in Moist-snuff products while Prevotella, Faecalibacterium and Clostridium were high in Snus product S2. The genus Capnocytophaga and Bacillus was elevated in Q1 and Q2 product, respectively (Supplementary Figure 6).

To identify potential interactions within the moist STPs bacterial communities, a co-occurrence network at the level of genus was constructed using a compositional robust method SparCC correlation coefficient that formulate a strong assumption of a sparse correlation network (Friedman and Alm, 2012). Overall, the prevalence of 247 genera were found to be considerably different between the groups of moist STPs and to aid interpretation, nodes were colored according to their phylum (Supplementary Figure 7A). Altogether, 245 positive and 238 negative considerable correlations (coefficient correlation > 0.3 and p-value < 0.05) were observed between 247 genera (Supplementary Table 5). The co-occurrence pattern of relevant genera, significantly associated with pre-cancer lesion or oral cancer, in moist STPs were examined in detail (Halboub et al., 2020; Sarkar et al., 2021; Srivastava et al., 2021). The genus Prevotella was found to be correlated positively with Clostridium (SparCC = 0.9719, p = 0.0099), Bifidobacterium (SparCC = 0.9332, p = 0.0099), Treponema (SparCC = 0.9332, p = 0.0099), Mycobacterium (SparCC = 0.8094, p = 0.0099), Rhodococcus (SparCC = 0.7948, p = 0.0099), and Lactobacillus (SparCC = 0.7294, p = 0.0396) while it was negatively correlated with Corynebacterium (SparCC = −0.9622, p = 0.0099) and Capnocytophaga (SparCC = −0.6792, p = 0.0297) (Figure 4 and Supplementary Table 5). Streptococcus was positively correlated with Veillonella (SparCC = 0.9237, p = 0.0099) and Haemophilus (SparCC = 0.8739, p = 0.0099), whereas Fusobacterium showed positive association with Capnocytophaga (SparCC = 0.7415, p = 0.0297) and Lautropia (SparCC = 0.9473, p = 0.0099). Another important genera Pseudomonas displayed positive concurrence with Treponema (SparCC = 0.9149, p = 0.0198), Mycobacterium (SparCC = 0.9027, p = 0.0099), Haemophilus (SparCC = 0.8982, p = 0.0297), Prevotella (SparCC = 0.8000, p = 0.0099), Acinetobacter (SparCC = 0.8247, p = 0.0099), Bifidobacterium (SparCC = 0.8186, p = 0.0099), Streptomyces (SparCC = 0.8592, p = 0.0099), and Clostridium (SparCC = 0.8345, p = 0.0099) (Figure 4 and Supplementary Table 5). The correlation network plot is interactive and the genera Prevotella showed the highest abundance in Snus products compared to Khaini, Moist-snuff, and Qiwam (Supplementary Figure 7B).

The robust biomarker of moist STPs was identified using a non-parametric statistical method LEfSe (Segata et al., 2011). LEfSe method discovers features with considerable differential abundance across the moist STPs and established the biomarker bacteriome at the genus level. Fifteen significant taxa were recognized as per the cutoffs values: FDR-adjusted p-value < 0.1 and linear discriminant analysis (LDA) > 2.0 (Figure 5A). The LDA score for Acinetobacter and Lactococcus were highest in Khaini products whereas Streptococcus and Neisseria had high LDA scores in Moist-snuff samples (Figure 5A). Further, the LDA score of Sphingobacterium, Janthinobacterium, and Pseudomonas was dominant in Qiwam, whilst that of Prevotella, Faecalibacterium, Oscillibacter, Ruminococcus, Blautia, Clostridium, Megasphaera, and Clostridioides was highest in Snus products (Figure 5A).

Next, to identify bacterial genera that differentiate between phenotypes, a random forest algorithm was applied to bacteriome data (Knights et al., 2011). The significant genera identified with random forest showed a pattern of changes across different groups of moist STPs (Figure 5B). At the genus level, the random forest model brought up Sphingobacterium, Ruminococcus, Blautia, Lactococcus, Janthinobacterium, Streptococcus, Moraxella, Achromobacter, Veillonella, Haemophilus, Faecalibacterium, Clostridioides, Acidiphilium, Prevotella, and Megasphaera (Figure 5B). The genera Sphingobacterium, Ruminococcus, and Blautia were the most decisive discriminated genera because of their higher predictive values.

The 16s rDNA sequencing data were used to infer the metabolic potential of bacteriome with PICRUSt which is derived from a phylogenetic distance or sequence similarity of identified microbes with the related microorganism whose whole genome has been sequenced and based on Greengenes annotated OTUs (Langille et al., 2013). The MicrobiomeAnalyst utilizes the PICRUSt algorithm and metagenome contributions were computed for all moist STPs based on Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology (KO term) (Kanehisa et al., 2013; Dhariwal et al., 2017). The result containing KO abundance level was generated and in sum 3695 KO terms were observed in the imputed metagenome of moist STPs (Supplementary Table 6). Several KEGG metabolic pathways were identified in moist STPs and their relative abundance in different moist STPs was monitored by MicrobiomeAnalyst (Supplementary Figures 8A,B).

Microbial reduction of nitrates to nitrite that leads to formation of TSNAs involves the nitrogen metabolism pathway.

The two important pathways involved in the extracellular accumulation of nitrite are (i) dissimilatory nitrate reduction nar operon that includes regulators (narXL), transporters (narK) and nitrate reductases (narZHJI); (ii) periplasmic nitrate reductase nap operon (González et al., 2006). The dissimilatory nitrate reduction pathway genes (narK, narZ, narJ, narI, and narH) were abundant in the Q3 product and Q1, Q2, and MS2 products showed a significant number of imputed genes of the nitrate reduction pathway (Figure 6 and Supplementary Table 6). The periplasmic nitrate reductases gene napA was prevalent in one Snus product (S2) and noticeably present in Q3, Q1, and MS2 products. The assimilatory nitrate reductases genes nasA and nasB were also predicted by PICRUSt and abundantly monitored in Qiwam (Q1, Q2, and Q3), Moist-snuff (MS2) and Snus (S2) products. The predicted genes of nitrite reductases including nirA, nirB, nirD, nirK, and nrfA were found to be abundant in all three Qiwam products and one Moist-snuff (MS2) product, whereas nirB was also prevalent in one Snus product (S2) (Figure 6 and Supplementary Table 6). Another important step in nitrogen metabolism is denitrification in which nitrogenous compounds (nitrate, nitrite, and ammonia) were converted to nitrogen gas (N₂). All moist STPs contained imputed genes related to denitrification nosZ (K00376), norB (K04561), and norC (K02305) but their abundance was very low (Supplementary Table 4). Further, nitrogen fixation related genes nifD1 (K02586), nifH (K02588), and nifK (K02591) were present in all moist STPs and the abundance of nifD1 and nifK were high in Q3, Q1, MS2, Q2, and K1 whereas nifH was prevalent in Q3, K2, K3, Q1, MS2, Q2, and K1 (Figure 6 and Supplementary Table 6).

The imputed metagenome of moist STPs through PICRUSt assists in the prediction of genes involved in antibiotic drug resistance. The highest number of antibiotic resistance genes was observed in all three products of Qiwam, one Moist-snuff product (MS2) and one Snus product (S2) as compared to Khaini products (K1, K2, and K3), MS1, MS3, and S1 products (Figure 7 and Supplementary Table 6). The multiple antibiotic resistance protein marc (K05595), penicillin-binding protein 1A mrcA (K05366) and penicillin-binding protein-4 (serine-type D-Ala-D-Ala carboxypeptidase/endopeptidase) dacB (K07259) was the two most abundant antibiotic resistance genes (Supplementary Table 6). The Q3 product showed genes having higher predicted prevalence such as macrolide-specific efflux system (macA, K13888), macrolide phosphotransferase (mph, K06979), macrolide efflux protein (ykuC, K08217), zinc D-Ala-D-Ala carboxypeptidase (vanY, K07260) and multidrug resistance protein (bcr/tcaB, K07552) (Figure 7 and Supplementary Table 6). The products Q2 and MS2 showed the prevalence of antibiotic resistance genes like multidrug-resistance protein (mdtG, K08161), aminoglycoside 3-N-acetyltransferase (aacC, K00662), penicillin-binding protein 1B (mrcB, K05365), and penicillin-binding protein 5/6 (dacC/dacA, K07258) (Figure 7 and Supplementary Table 6). The sample Q1 displayed the presence of zinc D-Ala-D-Ala dipeptidase (K08641), beta-lactamase induction signal transducer (K082184), and multidrug resistance proteins (emrB, K03446 and norB, K08170). The S2 products displayed antibiotic resistance genes including vancomycin resistance associated response regulator (vraR, K07694), small multi-drug resistance pump (emrE, K03297), macrolide transport system ATP binding/permease protein (macA, K05685), fosmidomycin resistance protein (tsr, K08223), and penicillin-binding protein activator (lpoB, K07337) (Figure 7 and Supplementary Table 6).

Gram-negative bacteria lipopolysaccharide (LPS) is involved in the progression and migration of oral squamous cell carcinoma (OSCC) (He et al., 2015). The LPS inflammatory activity is due to the Lipid-A component that is known to activate pro-inflammatory cytokines (Zhang and Ghosh, 2000; Karpiński, 2019). The ABC transporter complex (lptBFG) responsible for LPS transport from the inner to the outer membrane was found in moist STPs and its abundance was high in Qiwam and MS2 products. In moist STPs several genes related to LPS biosynthesis were observed (Figure 8 and Supplementary Table 6). The abundance of genes Lipid-A-disaccharide synthase (lpxB) and Kdo2-lipid IVA lauroyltransferase/acyltransferase (lpxL) which were involved in LPS synthesis significantly increased in Qiwam (Q1, Q2, and Q3), MS2 and S2 products (Figure 8 and Supplementary Table 6).

Another potent pro-inflammatory molecule is Flagellin which participates in the motility of bacteria. Flagellin coding gene fliC was found in all moist STPs and their abundance was high in S2, Q3, Q2, MS2, and Q1 as compared to K2, K1, S1, MS3, K3, and MS1 (Figure 8 and Supplementary Table 6). However, the genes involved in the synthesis of other pro-inflammatory molecules like peptidoglycan, teichoic acid and lipoteichoic acid (K03739 and K03740) were not observed in any moist STPs (Supplementary Table 6).

Further, an important category of genes is related to bacterial toxins because these toxins can participate in the pathogenesis of inflammation leading to several diseases (toxinoses) and can induce genomic damage that can lead to neoplastic transformation of epithelial cells (La Rosa et al., 2020). The STP-associated bacterial toxins identification based on KEGG orthology displayed that bacterial species having toxins genes were present in moist STPs. The highest sequence hits were observed for alpha-Hemolysin/cyclolysin transporter gene tolC (K12340), followed by putative hemolysin gene (tlyC, K03699), phospholipase-C gene (plcC, K01114), Neuraminidase1 gene (NEU1, K01186), hemolysin III gene (hly III, K11068) and thiol-activated cytolysin gene (slo, K11031) (Figure 8 and Supplementary Table 6). The moist STPs Q1, Q2, MS2, and S2 showed a high abundance of toxin genes whereas Q3, K1, and K2 also contain a significant prevalence of these genes compared to K3, MS1, MS3, and S1 products (Figure 8 and Supplementary Table 6).

Next, we performed a Taxon Set Enrichment Analysis (TSEA) of all significant genera using the MicrobiomeAnalyst tool and examined them across 239 known taxon set linked with host-intrinsic factors such as age and diseases (Chong et al., 2020). The bacteriome of moist STPs showed strong correlations with colorectal cancer, Crohn’s disease and Rheumatoid arthritis (Table 2). Further, we parsed the bacteriome of moist STPs with 53 taxon sets associated with microbiome-intrinsic factors such as microbe motility and shape. The bacteriome of moist STPs displayed a significant correlation with indole producers and mucin degraders (Table 3).