Incense Burning is Associated with Human Oral Microbiota Composition

Incense burning is common worldwide and produces environmental toxicants that may influence health; however, biologic effects have been little studied. In 303 Emirati adults, we tested the hypothesis that incense use is linked to compositional changes in the oral microbiota that can be potentially significant for health. The oral microbiota was assessed by amplification of the bacterial 16S rRNA gene from mouthwash samples. Frequency of incense use was ascertained through a questionnaire and examined in relation to overall oral microbiota composition (PERMANOVA analysis), and to specific taxon abundances, by negative binomial generalized linear models. We found that exposure to incense burning was associated with higher microbial diversity (p < 0.013) and overall microbial compositional changes (PERMANOVA, p = 0.003). Our study also revealed that incense use was associated with significant changes in bacterial abundances (i.e. depletion of the dominant taxon Streptococcus), even in occasional users (once/week or less) implying that incense use impacts the oral microbiota even at low exposure levels. In summary, this first study suggests that incense burning alters the oral microbiota, potentially serving as an early biomarker of incense-related toxicities and related health consequences. Although a common indoor air pollutant, guidelines for control of incense use have yet to be developed.


Results
Demographics relative to incense use. Our analysis included 303 subjects who provided a valid informed consent, a completed questionnaire on incense use, and a mouthwash sample from the United Arab Emirates Healthy Future Study-Pilot (UAEHFS-pilot) ( Supplementary Fig. S1). The study participants included 6.6% never, 24.1% occasional (one time or less/week), 33.7% frequent (2 to 5 times/week) and 35.6% daily (5-7 times/week) incense users (Table 1). Incense use was reported significantly more often by females then males (p = 0.006). No significant differences in incense use were found with respect to age, education or tobacco smoking habits (p > 0.05).
Richness, diversity and overall oral microbial community structure. After filtering of poor quality sequences and exclusion of low count OTUs, we obtained a dataset of 13,470,071 high quality 16S rRNA sequences representing 13 phyla, 20 classes, 26 orders, 41 families, 57 genera and 1102 OTUs. Measures of richness (observed) and diversity (Shannon entropy) revealed that diversity of the oral microbiota was significantly increased in daily incense users, compared to never (p = 0.013) and occasional (p = 0.011) users (Fig. 1). Overall oral microbiota structure and composition was also significantly different among the incense frequency use groups, largely driven by differentials between never and daily users of incense (p = 0.003, Fig. 1C), as revealed by Permutational multivariate analysis of variance (PERMANOVA) of UniFrac weighted distances adjusting for age, gender, smoking habits and assay batch. Interestingly, gender (whether modeled as a single factor or via an interaction with incense use) was not significantly associated with microbiota composition (Supplementary Table S1).
Differences of taxa relative abundances as exposure to incense increases. We investigated whether increased incense exposure was associated with differences in the relative abundance of specific bacterial taxa in the oral microbial community using negative binomial generalized linear models, as implemented by the DESeq function in the DESeq2 R package 20 (Fig. 2, Table 2, Supplementary Tables S2 and S3) adjusting for age, gender, smoking habits and assay batch. Trends of depletion and enrichment, with respect to incense use, tended to follow taxonomic affiliations within the phyla, particularly for Firmicutes (within the class Clostridia, 15 out of 16 taxa (94%) were enriched, and within Bacilli, 8 out of 9 (88%) were depleted) and Proteobacteria (all taxa but one (92%) were depleted) ( www.nature.com/scientificreports www.nature.com/scientificreports/ Greater incense use was associated with depletion of the Firmicutes class Bacilli, largely due to depletion of the high-abundance genus Streptococcus, while Firmicutes order Clostridiales exhibited increased abundance, related partially to a number of minor genera (Schwartzia, Mogibacterium, Parvimonas, Peptostreptococcus and Selenomonas) ( Table 2, Fig. 2, Supplementary Table S2). Greater incense use was associated with increased abundance of Bacteroidetes order Bacteroidales and its genus [Prevotella], while the Bacteroidales genus Paludibacter and order Flavobacteriales genus Capnocytophaga exhibited decreased abundance. Among Proteobacteria, genera Aggregatibacter and Cardiobacterium exhibited depletion; only the genus Bifidobacterium in Actinobacteria was identified as significantly enriched as exposure to incense increased. The same trend analysis across incense groups was independently performed for nonsmokers (n = 190) and smokers (n = 83); the direction of association between incense use and oral bacterial abundance was largely replicated in both of these groups (e.g., Streptococcus, Clostridiales, [Prevotella]), although some inconsistent associations between these two groups were observed (e.g., Mogibacterium, Peptostreptococcus and Cardiobacterium) (Supplementary Table S3). A variance partitioning analysis 21 attributed 0.7% of the variance to incense use, and 0.6% to smoking, with only 0.09% shared by both variables. These two independent variance components were highly statistically significant (p ≤ 0.005). We also found that even occasional incense use (once a week or less) tended to be associated with taxon abundance differentials (never vs. occasional, Supplementary Table S4).
Correlation networks amongst selected taxa. We built a correlation network with all incense users by means of the nonparametric Spearman correlations coefficient with the thirteen genera that were identified as being significantly associated with incense exposure in the trend analysis (Table 2, Fig. 3). A correlation network allows to detect whether changes in taxa abundances are occurring in the same direction or not. Two subcomponents in the correlation network were observed. The first one comprised most of the depleted taxa, Aggregatibacter, Cardiobacterium, Capnocytophaga, Paludibacter, all of which are positively correlated. The second and larger one comprised six of the eight enriched taxa and Streptococcus (depleted with incense use). The latter subcomponent comprised all six Firmicutes identified as associated with incense exposure and the Bacteroidetes [Prevotella]. Streptococcus was negatively associated with most of the other members of this network subcomponent. No significant correlations involving the taxa Bifidobacterium and Desulfobulbus were observed.

Discussion
This study explored for the first time the potential effect that household incense use has on the oral microbiota. We found that incense use was associated with increased diversity and changes in overall structure and composition of the oral microbial community. As the frequency of incense use increased, several taxa tended toward lower abundance, including Streptococcus, the dominant genus of the oral microbiota, and Paludibacter, www.nature.com/scientificreports www.nature.com/scientificreports/ Capnocytophaga, Aggregatibacter and Cardiobacterium, while others tended toward increased abundance, including Bacteroidales, and its genus [Prevotella], and the Actinobacteria genus Bifidobacterium. Overall, bacterial abundances were altered with increasing levels of incense use, but notably abundance differentials were observed even in occasional users (once a week or less), suggesting that incense use impacts the oral microbiota even at relatively low levels of exposure.
Incense burning has a long history in many cultures and is commonly used for religious ceremonies, aromatizing and other reasons. In the Arabian gulf, and the UAE in particular, incense is frequently burned in households, but also in mosques and shopping malls, being one of the most common sources of indoor smoke exposure 4,22 . Incense burning produces fine and ultrafine airborne particulates in large quantities when compared to other indoor air pollutant sources (i.e. cooking) and emits gaseous pollutants, including carbon monoxide (CO), sulfur dioxide (SO 2 ), nitric oxide (NO) and volatile organic compounds (VOCs) in relatively high concentrations 1,8 , as well as other aromatic, irritant and toxic compounds 10 .
Incense smoke presents many similarities in composition to tobacco smoke. Hence, exposure to incense could potentially have similar effects on the oral environment and on the function and composition of saliva, as observed in tobacco smokers 23,24 . Thus, we could expect that incense use, as with tobacco use, could result in the reduction of saliva by the depletion of oxygen 25 , increased acidity 26,27 , decreased activity and abundances of salivary proteins 28,29 , and lower abundance of Immunoglobulin A (IgA) molecules 30 . We recently reported 19 that Only taxa that were identified in the trend analysis by DESeq. 2 with statistical significance (q < 0.1, i.e., after FDR correction) are reported. Red and green nodes indicate enrichment and depletion, respectively, for that particular taxon as exposure to incense increases. The first inner ring represents a heatmap showing mean abundance of the genera adjusted for size factors in DESeq2 (>0.1%; 0.1-0.99%; 1-5%; 5-10%; >10%). Bars on the external ring represent log2 fold changes of taxa mean counts per unit of increased incense exposure at the genus level. Red bars indicate enrichment of a genus while green bars indicate depletion, as detected by the DESeq2 analysis. Triangles underneath the log2 fold change bars indicate statistical significance with q < 0.1. The cladogram was built using Graphlan 53  www.nature.com/scientificreports www.nature.com/scientificreports/ cigarette use in this population is strongly related to Proteobacteria and Fusobacteria depletion and enrichment of several taxa within the phyla Synergistetes and Actinobacteria, among others 19 . Despite incense and tobacco smoke potentially having similar effects on the oral environment, our results indicate that the effect of incense exposure on microbiota composition differ from those of tobacco exposure. First, the impacted bacteria are largely dissimilar from those associated with tobacco use and second, the microbial differentials observed with incense use are similar whether an individual used tobacco products or not. In fact, our variance partitioning showed that the effects of incense use and smoking on the oral microbiota are largely independent. This could in part be due  Table 2. Selected differentially abundant taxa of the oral microbiome as exposure to incense increases. a Only those taxa that had a significantly differential abundance with q < 0.1 and a Cook's distance < 10 for the trend analysis are shown. b Mean values refer to mean normalized counts of taxa according to incense burning group. c FDR adjusted p value implemented independently at each level (i.e. phylum, class …). d Occasional users report to burn incense in the household one time or less a week. e Frequent users report to burn incense in the household 2-5 times a week. f Daily users report to burn incense in the household 5-7 times a week. *Trend analysis corresponds to the log2 fold change per unit of change of the continuous-valued incense variable. www.nature.com/scientificreports www.nature.com/scientificreports/ to toxicants present in incense use that are absent in tobacco smoke and that might have differential cytotoxic or mutagenic effects on the bacterial community 10 , potentially affecting the oral environment in distinctive ways.
One of the most striking results of our study was the observed depletion of the most abundant taxa of the oral microbiota, the Streptococci. These are, for the most part, commensal bacteria with high ability to adhere to hard and mucosal oral surfaces, as most possess a diverse array of adhesins 31 . If incense use decreases the activity and abundance of salivary proteins 28,29 , this may compromise the ability of Streptococci to efficiently adhere and colonize oral surfaces and instead, facilitate their agglutination and removal by swallowing 31 . In addition to enabling bacterial adherence, proteins present in the saliva are also an important source of nutrients. When salivary proteins become a limiting factor, highly proteolytic bacteria will outcompete low proteolytic bacteria 32 . This is in agreement with the patterns observed ( Table 2, Supplementary Table S2), as depleted Streptococci have a low proteolytic potential 33 , while Bifidobacterium (Actinobacteria), Parvimonas and Peptostreptococcus (Firmicutes; Clostridia), which comprise highly proteolytic species [34][35][36] were enriched.
Bacteria are organisms that interact within and between species forming complex communities that influence the niche they inhabit and that can respond to external stimuli as well 37 . The interactions between members of a bacterial community are for the most part non-random and thus, shape the structure and composition of that particular community 38 . While synergistic interactions between members of the community will stimulate the growth of one another, antagonist effects will lead to the depletion of one or more of the organisms 32 . For example, a synergistic interaction is known for Prevotella, as it interacts with Peptostreptococcus by providing amino acids to the medium which are used by Peptostreptococcus promoting its growth 39 . Similarly, Schwartzia is a bacterium limited by being able to metabolize only succinate 40 , but its association with a succinate producer such as Selenomonas will ensure that the growth of one will enable the growth of the other 40 . Both of the latter associations (Prevotella-Peptostreptococcus and Schwartzia-Selenomomas) were observed as positive correlations in the network analysis (Fig. 3). Nevertheless, further research is needed, as positive correlations could also mean that both taxa respond in a similar way to a perturbation or that they consume similar resources but are in different niches of the oral cavity.
This study is the first to our knowledge that investigates the association between use of incense and the oral microbiome. We recognize that the amount of variance explained by incense was not large, however, this likely reflects the large individual differences in overall microbiome composition among the participants in our study. The use of the structured questionnaire was important to reveal the extent to which different frequencies of incense burning in the household where associated to oral microbiome dysbiosis. Our study had some limitations such as (1) the low number of nonusers of incense, since household incense use is traditional in UAE households, (2) the absence of data related to dental history and dental hygiene habits, an important factor influencing the presence of particular non-pathogenic and pathogenic bacteria, (3) the lack of information regarding the type of incense used and (4) the low number of participants that smoke but did not burn incense (only two), preventing the comparison between exclusively incense users and exclusively smokers. Also, we were unable to track the progression of microbiota changes with continued incense use, as only one oral sample was collected per individual. Finally, although 16S rRNA sequencing provides information on the bacterial composition and abundances of the oral taxa, information on functional potential of these microbes can only indirectly be inferred from these data. The latter will be assessed in future studies.
In summary, our results demonstrate that indoor incense use is associated with oral microbiota structure and composition, even with relatively infrequent use (as observed in occasional users). Furthermore, this association is largely independent of tobacco-related microbial perturbations (as no significant differences were observed between nonsmoker and smoker incense users). As the oral microbiome serves important functions in health 12,13 , the observed impact on oral microbiota may serve as an early biomarker of incense-related toxicities and related health consequences. Although an important indoor air pollutant, guidelines for control of incense use have yet to be developed.

Materials and Methods
Study subjects. This study used mouthwash samples and metadata from 303 subjects recruited by the UAEHFS-pilot study undertaken between the months of December 2014 and April 2015 41 . Study subjects were eligible to participate if they were Emirati nationals aged 18 and above. Participants completed physical and clinical exams, including measurements of body composition and blood pressure, provided blood, urine and mouthwash samples and completed a self-administered questionnaire including information on socio-demographic factors, lifestyle and medical history 41 . From 517 consented study subjects, 303 subjects completed baseline questionnaires including incense data and provided mouthwash samples (Supplementary Fig. S1). This study was approved by the Institutional Review Boards (IRB) of Sheikh Khalifa Medical City (SKMC), Zayed Military Hospital (ZMH), New York University Abu Dhabi (NYUAD) and NYU Langone Health, New York. All individuals participating in the study read and signed an informed consent. All experiments were performed in accordance with relevant guidelines and regulations.
Measurements. Incense burning classification of study subjects. Incense burning habits were ascertained through a self-reported structured questionnaire and were classified as follow; Never users, occasional users (usually one or less times a week), frequent users (2 to 5 times a week) and daily users (5 to 7 times a week). (Supplementary Metadata).

Definition of smoker and nonsmoker.
Smoker and nonsmoker status was defined as in 19 . Briefly, study subjects provided detailed information on smoking habits in a structured questionnaire. In addition, study subjects provided urine samples that were used to test for presence of cotinine, a nicotine metabolite, using the COT rapid www.nature.com/scientificreports www.nature.com/scientificreports/ test cassette (International Biomedical Supplies), with a cut off concentration of 200 ng/ml. Smokers were defined as all study subjects that self-reported as smoker independently of cotinine results. Nonsmokers were those study subjects that self-reported as nonsmokers and were further validated by a cotinine negative result.
Mouthwash sample collection. As reported previously 19 study subjects were provided with 10 ml of pharmaceutical grade normal saline (0.9%) solution and asked to vigorously swish for 30 seconds and spit it out into a collection tube. Samples were initially stored at 4 °C. Once in the lab, samples were vortexed for 20 seconds, pipetted up and down and aliquoted into 1 ml cryotubes to be stored at −80 °C until further processing. DNA extractions were performed on two blank saline samples alongside the mouthwash samples to confirm that the saline solution used for collection of the mouthwash samples contained no detectable levels of DNA. DNA concentrations were quantified using the high sensitivity Qubit assay. Only mouthwash samples yielded measurable amounts of DNA 19 .
Microbiome assay. Briefly, 2 ml were used per mouthwash sample for DNA extraction. Cell pellets were collected by spinning samples at 6000 g for 3 min and then at 10000 g for 10 min. DNA was extracted using the Mo BioPowerSoil PowerLyzer kit following manufacturer's instructions (Mo Bio Laboratory Inc, California, USA). Genomic DNA was visualized on a gel and quantified using the Qubit HS kit (Thermo Fisher Scientific). Amplification of DNA from the V4 region of 16S rRNA gene (515F -5′GTGCCAGCMGCCGCGGTAA3′ -806R -5′GGACTACHVGGGTWTCTAAT3′) was performed using specifically designed primers with Illumina adaptor sequences and a 12 bp index (reverse primer only) added for posterior multiplexing. The FastStart enzyme (Roche, IN) was used to amplify the 16S rRNA gene for 27 cycles using approximately 12.5 ng of DNA per sample. PCR products were visualized in an agarose gel, purified using Agencourt AMPure beads (Beckman Coulter Life Sciences, IN) and quantified using the Qubit BR kit (Thermo Fisher Scientific). Samples were then pooled for sequencing on an Illumina MiSeq platform in two batches 42 .
Quality control. As previously reported, samples were sequenced in two batches 42 . In addition, three samples were triplicated in each batch to assess sample quality control. A blank was used in each batch for both DNA extractions and PCR amplifications. DNA concentrations were measured each time to verify that blanks were devoid of significant amounts of DNA. Only one of them was sequenced with each sequencing batch. Quality control samples showed good consistency, with the coefficient of variability ranging from 1.65-2.32% for the Shannon entropy and 1.02-7.21% for specific phyla relative abundances 19 . Statistical analysis. Sequence processing and taxonomic assignment. The QIIME 1 bioinformatic pipeline was used to process the sequence data. Sequences were de-multiplexed and trimmed using the split_libraries_ fastq.py QIIME script 43 . Poor quality sequences were excluded from further analyses (minimum average base score quality per read was 20, minimum read length 200 bp and no mismatches in adaptor or barcode sequences was permitted). The pick_de_novo_otus.py workflow as implemented in QIIME 43 to cluster sequences into operational taxonomical units (OTU) using a 97% pairwise-identity cutoff. UCLUST 44 , PyNAST 45 and the Greengenes database were further implemented and used respectively to obtain taxonomical assignment of the sequences. Removal of chimeric sequences was accomplished using ChimeraSlayer as implemented in the QIIME workflow 46 . Low count OTUs were filtered from the analyses if they were singletons and absent in more than 10% of the samples.
Alpha diversity. Richness and Shannon entropy were estimated to assess the within-sample diversity. Both indexes were estimated for 200 iterations of rarefied OTU datasets (16738 sequences per sample) followed by computing the average for each sample using the vegan library in R 47 . Multiple linear regression adjusting for age, gender, smoking status and batch was implemented to compare alpha-diversity scores between different incense burning groups (never, occasionally, frequently and daily).
Microbial composition according to incense burning. We examined the relationship between oral microbiome composition and incense burning by conducting a permutational multivariate analysis of variance (PERMANOVA) of weighted (taxa relative abundances) and unweighted (binary input, absence/presence) UniFrac distance matrices 47 . The UniFrac distance metric incorporates phylogenetic distances between community organisms into the calculations of the dissimilarity matrices and were computed using the UniFrac function in the Phyloseq library in R 48,49 . The latter was visualized through Principal Coordinate Analyses (PCoA) after correction for negative eigenvalues. All PERMANOVA analyses were adjusted for age, gender, smoking status and sequencing batch, and were performed using the Adonis function in the vegan R library 47 .
Identification of differences in taxa abundances as exposure to incense burning increases. We studied differences of taxa relative abundances associated to increasing incense exposure using a model based on negative binomial distribution as implemented by the DESeq function in the DESeq2 R package 20 . The trend analysis was executed by using the continuous-valued incense variable. Since the average use of incense by the participants was not reported in the questionnaire, value 1 was attributed to "Never", 2 to "Occasionally", 3 to "Frequently, and 4 to "Daily" users. Log 2 fold change results reported for the trend analysis correspond to the log 2 fold change per unit of change of the continuous-valued incense variable. Comparisons between incense frequency categories were done by including in the model the contrast argument and specifying the comparison of interest. A p-value < 0.05 was considered of nominal statistical significance, and an FDR-adjusted 50 p-value < 0.10 (named hereafter q-value) was considered significant after adjustment for multiple comparisons. Analyses were done at all taxonomical levels and models were adjusted for age, gender, smoking status and sequencing batch. All analyses