This study was carried out according to the National Guide for the Care and Use of Laboratory Animals. All experimental procedures were approved by the Ethics Committee of Nanjing Stomatological Hospital, Medical School of Nanjing University [IRB Approval Number: 2018NL-008 (KS)] and the Animal Ethical and Welfare Committee of Nanjing University (IACUC-D2102043). Five-week-old Wistar rats, weighing 145–155 g were provided by Vital River Laboratory Animal Technology Co., Ltd (Beijing, China). After arrival at our facilities, each rat was housed in its own cage in a specific pathogen-free facility to avoid microbiota transfer as previously described (Laukens et al., 2016). Rats were housed in cages under controlled ambient temperature (22 ± 3°C) and humidity with a strict 12-h light/12-h dark cycle. Rats were provided with tap water and a standard laboratory chow diet. After one week of acclimatization, the rats were divided into two groups with equivalent average body weights and standard deviations. The rats in the control group (n = 10, male 5 and female 5) were given intraoral delivery of room air. The rats in the active smoking group (n = 9, male 5 and female 4) were given intraoral delivery of cigarette smoke. The rats in the control group and active smoking group were separated in different rooms. Both groups of rats were kept in the animal retainers twice a day (20 min per time in the morning and afternoon) to receive intraoral delivery of room air or cigarette smoke. The total duration of the experiment was 3 months. 3R4F Research Cigarettes (Kentucky Tobacco Research Institute, Lexington, KY, USA) were used in this experiment. Each cigarette contained the following ingredients: tar 9.4 mg, nicotine 0.73 mg, carbon monoxide 12.0 mg.

We developed an intraoral smoking exposure apparatus for rodents that simulates human active smoking. Briefly, the intraoral smoking exposure apparatus included animal retainers, intraoral smoking pipes, cigarette smoke delivery tubes, peristaltic pumps (LongerPumper, Boonton, NJ, USA), and a hood over the smoking rats to evacuate the extra smoke from the environment. During the treatment process, the body of the rats was retained in a plastic hollow cylinder. The delivery tube was connected to the filter tip of a lit cigarette at one end and connected to the intraoral smoking pipe at the other end. The intraoral smoking pipe was sterile and replaced between rats. Then, a peristaltic pump at a speed of 30 rpm was used to automatically draw, deliver and spray cigarette smoke into the oral cavity of rats. Approximately one week after the first smoking treatment, the rats became addicted to active smoking, and it could be observed that the rats smoked with the oral cavity by regular cheek blowing. A video as supplemental information shows that a rat was actively smoking.

After 3 months, microbial samples were collected from the oral cavity by swabbing over the buccal mucosa of rats. The animals were then deeply anesthetized via intraperitoneal injection of sodium pentobarbital (2%, 0.2 ml/100g). The abdomen was sterilized with 75% ethanol, and a midline abdominal incision was made. Blood was collected via abdominal main vein puncture and centrifuged for further serum biochemical analyses. Then, the rats were sacrificed. Subsequently, the stomach, small intestine, cecum and colon were isolated. The gastrointestinal contents were rapidly collected through a sterile incision in the stomach (gastric antrum), small intestine (middle segment of the jejunum), cecum (body of cecum), and colon (middle segment) while avoiding junction sites of these regions (Li et al., 2017). The serum, swab, and content samples were rapidly collected, flash-frozen with liquid nitrogen, and immediately stored at -80°C until analysis.

Serum level of cotinine, the primary metabolite of nicotine, was assayed by using liquid chromatography/mass spectrometry. Serum corticosterone level was measured by rat enzyme immunoassay assay kits (Cayman Chemical, Ann Arbor, MI, USA) according to the manufacturer’s instructions. Serum epinephrine level was measured by rat enzyme immunoassay assay kits (Biovision, Milpitas, CA, USA) according to the manufacturer’s instructions. Serum gastrin 1 level was measured by rat enzyme immunoassay assay kits (Abcam, Cambridge, MA, USA) according to the manufacturer’s instructions. Serum glucose level was determined by rat glucose colorimetric assay kits (Cayman Chemical, Ann Arbor, MI, USA) according to the manufacturer’s instructions. Serum GHbA1c level was measured by rat enzyme-linked immune sorbent assay kits (Cusabio, Houston, TX, USA) according to the manufacturer’s instructions. Serum insulin, leptin, and unacylated ghrelin levels were measured by rat enzyme immunoassay assay kits (Cayman Chemical, Ann Arbor, MI, USA) according to the manufacturer’s instructions. Serum adiponectin level was measured by rat enzyme immunoassay assay kits (Abcam, Cambridge, MA, USA) according to the manufacturer’s instructions.

DNA extraction was performed as described previously (Clos-Garcia et al., 2019). Bacterial genomic DNA was extracted from the buccal swab and gastrointestinal content samples using a DNA extraction kit (Omega Bio-tek, Norcross, GA, USA) according to the manufacturer’s protocols. Agarose gel electrophoresis was used to assess the DNA integrity and size. DNA concentration and purity were evaluated using a Qubit 2.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The V4 region of the bacterial 16S ribosomal RNA gene was amplified by PCR using the primers 515 F (GTG CCA GCM GCC GCG GTA A) and 806 R (GGA CTA CHV GGG TWT CTA AT) (Kaakoush et al., 2017).

The PCR products were separated by gel electrophoresis and purified using an AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) according to the manufacturer’s instructions. Bacterial 16S rDNA libraries were constructed according to the 16S Metagenomic Sequencing Library Preparation guide from Illumina (Forest City, CA, USA). Only libraries without primer dimers and contaminant bands tested with an Agilent High Sensitivity DNA Kit were used for sequencing via Illumina MiSeq. The equimolar purified products were pooled and paired-end sequenced (2 × 250) on an Illumina MiSeq platform according to standard protocols at BGI-Shenzhen (Shenzhen, China). Prior to submission for sequencing, libraries were quality checked. As an added quality control measure, the software package QIIME pipeline was used to filter out and discard poor-quality sequence reads using the default settings (Liu et al., 2018).

Operational Taxonomic Units (OTUs) based on the 16S rRNA gene sequences were determined using QIIME’s UCLUST algorithm, and taxonomy assignment was performed using the RDP classifier (Bokulich et al., 2018). The relative abundance of each OTU was determined as a proportion of the sum of sequences for each sample.

Alpha diversity analyses estimated from observed species (Sobs), Chao 1, abundance-based coverage estimate (ACE), Shannon, and Simpson indices-based measurements were used to evaluate microbial richness, evenness and community diversity based on the OTUs. The alpha diversity metric was determined using mothur software v.1.31.2, and rarefaction curves were visualized with R program v.3.1.1 (Hall et al., 2017; Lopez-Garcia et al., 2018).

Beta diversity analyses were used to determine the degree of dissimilarity between pairs of bacterial communities using QIIME v.1.80. The unweighted, weighted UniFrac and Bray-Curtis distance matrix, which measures pairwise taxonomic dissimilarity between microbial populations, was analyzed with an unsupervised clustering algorithm (Wu et al., 2018). The beta diversity metric was visualized with a matrix heatmap constructed using the aheatmap function from the NMF package in R program v.3.1.1 (Gaujoux and Seoighe, 2010).

Taxonomic assignments of the microbiome populations were made according to the OTU annotation within QIIME v.1.80, presenting the taxonomic profiling of the samples at each taxonomic level (Liu et al., 2018).

Heatmaps of bacterial abundance showed the bacterial distribution in five regions along the digestive tract. Log10-transformed relative abundances of bacteria were visualized using the gplots package in R program v.3.1.1 (Li et al., 2020b).

The linear discriminant analysis (LDA) effect size (LEfSe) method was used for statistical analysis at different taxonomic levels (Segata et al., 2011). Cladogram constructed using the LEfSe method to indicate the phylogenetic distribution of active bacteria that were remarkably enriched. LDA scores showed significant bacterial differences within groups at different taxonomic levels.

The common and unique OTUs in different groups were determined and compared using the VennDiagram package in R program v.3.1.1 to generate the Venn diagram (Lam et al., 2016).

Network analysis was conducted to identify correlations between core genera. The inter-generic associations were identified using Spearman’s rho correlation coefficient. Network analyses were visualized in network plots generated by Cytoscape software v.3.5.1 (Xi et al., 2019).

The biological functions of the microbial community were predicted and annotated with KEGG pathways using the PICRUSt algorithm (Li et al., 2017). The OTU contribution to the abundances of functional KEGG pathways was calculated by summing the abundance contribution of each OTU for each KEGG orthology that belonged to the same pathway and then normalizing by the number of samples. Differentially enriched KEGG pathways or modules were identified according to their reporter score (Chen et al., 2017b).

Most biochemical data and the abundances of the bacteria are presented as the median and interquartile range. Non-parametric statistical methods were used throughout the study. Statistical comparison of two groups was performed with the nonparametric Mann-Whitney U-test as described previously (Chu et al., 2017; Clos-Garcia et al., 2019). In each case, p-values ≤ 0.05 were considered statistically significant. The analyses were performed with the Statistical Package for the Social Sciences, version 18.0 (SPSS, Chicago, USA).

Figure 1 shows a schematic diagram of the control and rat model of active smoking using in the present study. The serum cotinine concentration was 2.813 ± 0.944 ng/ml in rats in the active smoking group, while that in the control group rats was zero (p < 0.001) ( Table 1 ). There was no statistically significant differences in serum gastrin 1, corticosterone and epinephrine levels between the control group and active smoking group ( Table 1 ). Our results showed that rats in the active smoking group gained less weight compared with rats in the control group (p < 0.01) ( Supplementary Figure S1 ).

Biochemical analyses indicated that the serum glucose level in the active smoking group was significantly higher than that in the control group (p < 0.01). Moreover, the serum insulin and leptin levels in the active smoking group were remarkably lower than those in the control group (p < 0.01 and p < 0.05, respectively). However, the serum glycosylated hemoglobin A1c (GHbA1c), adiponectin and ghrelin levels were not significantly altered in the active smoking group compared with the control group (p > 0.05) ( Table 1 ).

Except for two colonic content samples that were unavailable in a rat in the control group and a rat in the active smoking group, all 93 specimens met the requirements for library establishment. In total, 4,245,900 sequence reads with a mean length of 252 bp were obtained from all the samples (the oral swabs and gastric, small intestinal, cecal, and colonic contents) in the control and active smoking groups. Each sample was covered by an average of 45,655 reads, arrange from 36,789 to 46,773. Figure 2 shows a schematic of five regions along the digestive tract and the microbiome constituents in each region in the control and active smoking groups.

Regardless of the smoking status, bacterial species richness and diversity were significantly different between most regions along the digestive tract. The cecal microbiota had the greatest diversity (596 unique OTUs per sample on average), followed by the colonic microbiota (582 unique OTUs per sample on average), and the oral microbiota had the lowest diversity (102 unique OTUs on average). The gastric and small intestinal microbiota shared the same diversity (both with 225 unique OTUs on average). Our analyses indicated that alpha diversity was generally increased along the digestive tract. The stomach and small intestine harbored a more diverse microbial community than the oral cavity, while the cecum and colon possessed a more diverse microbial community than the stomach and small intestine ( Figure 3A ).

The alpha diversity in each region did not remarkably differ between the control group and active smoking group (p > 0.05) ( Figure 3A ). Based on analyses of the unweighted UniFrac distance metrics, the bacterial communities in the gastric, small intestinal, cecal or colonic contents in the control and active smoking groups were markedly different in beta diversity (p < 0.001, p < 0.001, p < 0.05, and p < 0.001, respectively) ( Figure 3B ). Based on analyses of the weighted UniFrac distance metrics, the bacterial communities in the gastric or cecal contents in the control and active smoking groups were significantly different in beta diversity (p < 0.001 and p < 0.001, respectively) ( Figure 3B ). Based on the Bray-Curtis distance metrics, the bacterial communities in the gastric, cecal or colonic contents in the control and active smoking groups were also clearly distinct by beta diversity (p < 0.001, p < 0.001, and p < 0.05, respectively) ( Figure 3B ).

Our analyses indicated that the community structure of the stomach was analogous to that of the small intestine, while the cecum and colon shared a similar microbial community ( Figure 4 and Supplementary Figures S2, S3 ). Moreover, the community structure of the oral cavity was closer to that of the stomach and small intestine in the unweighted, weighted UniFrac or Bray-Curtis distance ( Figure 4 and Supplementary Figures S2, S3 ). A matrix heatmap shows distinctive clustering of the cecal and colonic bacterial communities by active smoking status based on the unweighted UniFrac distance analyses ( Figure 4 ). A matrix heatmap shows distinctive clustering of the oral cavity, gastric and small intesinal bacterial communities by active smoking status based on the Bray-Curtis distance analyses ( Supplementary Figure S2 ). A matrix heatmap shows relative weak clustering of the digestive tract bacterial communities by active smoking status based on the weighted UniFrac distance analyses ( Supplementary Figure S3 ).

Matrix heatmaps based on the unweighted UniFrac analyses show distinctive clustering in the control group and active smoking group in the oral cavity, cecal and colonic bacterial communities ( Figure 5 ). Matrix heatmaps based on the Bray-Curtis and weighted UniFrac analyses show distinctive clustering in the control group and active smoking group in the gastric, cecal and colonic bacterial communities ( Figure 5 ).

From the phylum to the species level, the taxonomic analysis revealed the enriched taxa in the control group or active smoking group ( Supplementary Tables S1-S5 ). Not surprisingly, both the relative abundance heatmap and taxonomic profiling bars showed that the overall microbial composition significantly varied along the digestive tract ( Figure 6 ). Despite this, Firmicutes was the most abundant phylum in the gastrointestinal regions.

At the phylum level, rats in the active smoking group possessed significantly lower relative abundances of Bacteroidetes (p < 0.01) and higher relative abundances of Cyanobacteria (p < 0.05) in the oral cavity. When comparing the composition of the gastric and small intestinal microbiota between the groups, there were no differences in abundance at the phylum level. Rats in the active smoking group bore significantly lower relative abundances of Firmicutes (p < 0.05) in the cecum. Rats in the active smoking group possessed significantly higher relative abundances of Tenericutes and TM7 (p < 0.05 and p < 0.05, respectively) in the colon ( Supplementary Tables S1–S5 ).

Overall, at the genus level, the results suggested that active smoking significantly altered the enrichment of low-abundance genera ( Supplementary Tables S1–S5 ). Genera with relative abundance higher than 5% seemed to not be directly affected by active smoking. Interestingly, our results showed that the proportions of the class Clostridia, the order Clostridiales and the order Turicibacterales, the family Clostridiaceae and the family Turicibacteraceae, the genus Clostridium and the genus Turicibacter, and the species Clostridium perfringens were reduced in the active smoking group compared with the control group ( Supplementary Tables S1–S5 ).

In the class level, the relative abundance of Clostridia was lowered in the oral, small intestinal and cecal samples from the active smoking group (p < 0.05, p < 0.05 and p = 0.05, respectively) ( Supplementary Tables S1, S3, S4 ). In the order level, the relative abundance of Clostridiales was decreased in the oral, small intestinal and cecal samples from the active smoking group (p < 0.05, p < 0.05 and p = 0.05, respectively) ( Supplementary Tables S1, S3 and S4 ). At the family level, the relative abundance of Clostridiaceae was decreased in the cecal samples from the active smoking group (p < 0.05) ( Supplementary Table S4 ). At the genus level, the relative abundance of Clostridium was reduced in the gastric and cecal samples from the active smoking group (p < 0.05 and p < 0.01, respectively) ( Figure 7 and Supplementary Tables S2, S4 ). At the species level, the relative abundance of Clostridium perfringens was lower in the gastric, cecal and colonic samples from the active smoking group (p < 0.05, p < 0.05, and p < 0.01, respectively) ( Supplementary Tables S2, S4, S5 ).

In the order level, the relative abundance of Turicibacterales was decreased in the oral, cecal and colonic samples from the active smoking group (p < 0.05, p < 0.05 and p < 0.05, respectively) ( Figure 7 , Supplementary Tables S1, S4, S5 ). At the family level, the relative abundance of Turicibacteraceae was decreased in the oral, cecal and colonic samples from the active smoking group (p < 0.05, p < 0.05 and p < 0.05, respectively) ( Figure 7 and Supplementary Tables S1, S4, S5 ). At the genus level, the relative abundance of Turicibacter was reduced in the oral, cecal and colonic samples from the active smoking group (p < 0.05, p < 0.05 and p < 0.05, respectively) ( Figure 7 and Supplementary Tables S1, S4, S5 ).

Besides reduced abundance of potentially beneficial genera (Clostridium and Turicibacter), our results also demonstrated a clear enrichment in sulfate-reducing bacteria (Desulfovibrio and Bilophila) and some opportunistic pathogens (Staphylococcus, Jeotgalicoccus, and Odoribacter). At the genus level, the relative abundance of Desulfovibrio was increased in the gastric and cecal samples from the active smoking group (p < 0.01 and p < 0.05, respectively) ( Supplementary Tables S2, S4 ). Also, the relative abundance of Bilophila was increased in the gastric and cecal samples from the active smoking group (p < 0.05 and p < 0.05, respectively) ( Supplementary Tables S2, S4 ). The genus Staphylococcus abundance was increased in the oral, cecal and colonic samples from the active smoking group (p < 0.01, p < 0.01 and p < 0.01, respectively) ( Supplementary Tables S1, S4, S5 ). Also, the genus Jeotgalicoccus abundance was increased in the oral and cecal samples from the active smoking group (p < 0.01 and p < 0.01, respectively) ( Supplementary Tables S1, S4 ). The genus Odoribacter abundance was increased in the cecal and colonic from the active smoking group (p < 0.001 and p < 0.01, respectively) ( Supplementary Tables S4, S5 ).

Surprisingly, the alpha diversity of the gastric microbiome was significantly higher than that of the small intestinal microbiome in the active smoking group (observed species and ACE metric) (p < 0.05), while those of the gastric and small intestinal microbiome in the control group were comparable (p > 0.05) ( Figure 8A ). Intriguingly, in the control group, no significant differences in the relative abundance of the Firmicutes, Bacteroidetes, Proteobacteria, Spirochaetes, and Tenericutes phyla were observed between the gastric contents and the small intestinal contents (p > 0.05). However, in the active smoking group, the relative abundances of these phyla were significantly distinct between the gastric contents and the small intestinal contents (p < 0.05) ( Figure 8B ). At the OTU level, fewer unique OTU-annotated bacterial taxa resided in the oral cavity, small intestine, and colon of the active smoking group compared with the control group. In contrast, more unique OTU-annotated bacterial taxa inhabited the stomach and cecum of the active smoking group compared with the control group ( Supplemental Figure S4 ).

In agreement with the taxonomic profiling and alpha diversity findings, the network analyses revealed prominent shifts in microbial community structure and altered relationships between the gastric and small intestinal flora ( Figure 9 ). It was also noted that robust inter-generic networks were translocated from the small intestine flora in the control group into the stomach flora in the active smoking group. Moreover, in contrast with the control group, the core genera of the stomach were significantly enriched while the core genera of the small intestine were profoundly depleted in the active smoking group. Although Lactobacillus was the keystone species of the stomach and small intestine in the control group and the active smoking group, its role was significantly reduced within the microbial associations in the stomach and conspicuously enhanced within the microbial associations in the small intestine in the active smoking group.

Functional prediction analyses ( Figure 10 ) suggested that microbiota associated with pathways for carbohydrate, lipid, fatty acid, and amino acid metabolism and lipid or amino acid biosynthesis was more pronounced in rats in the control group (p < 0.05). However, microbiota associated with metabolism pathways for some antioxidants and vitamins, such as glutathione, lipoic acid, retinol, taurine, hypotaurine, riboflavin, and cofactors, was enriched in the active smoking group (p < 0.05).