Oral Pathobiont-Induced Changes in Gut Microbiota Aggravate the Pathology of Nonalcoholic Fatty Liver Disease in Mice

This study was approved by the Institutional Animal Care and Use Committee of Niigata University (permit number; SA00328). All experiments were performed in accordance with the Regulations and Guidelines on Scientific and Ethical Care and Use of Laboratory Animals of the Science Council of Japan, enforced on June 1, 2006.

All authors had access to the study data and had reviewed and approved the final manuscript.

P. gingivalis strain W83 and P. intermedia ATCC25611 maintained in our laboratory were cultured in modified Gifu anaerobic medium broth (Nissui, Tokyo, Japan). Veillonella rogosae JCM 15642^T and Actinomyces naeslundii ATCC19039 obtained from Dr. Mashima at Aichi-Gakuin University, Nagoya Japan and maintained in our laboratory, respectively were cultured in brain–heart infusion broth. P. gingivalis and P. intermedia were used as periodontopathic bacteria whereas A. naeslundii and V. rogosae were used as commensal controls.

Six-week-old male C57BL/6N mice were obtained from Japan SLC (Shizuoka, Japan). After acclimatization under specific pathogen-free conditions and feeding regular chow and sterile water for 1 week, the mice were orally administered with any one of the following: vehicle (PBS with 2% carboxymethyl cellulose; Sigma-Aldrich, St. Louis, MO), or a total of 1 × 10⁹ colony-forming units of each bacterial species suspended in vehicle through a feeding needle five times a week for 3 weeks. The number of administered bacteria was determined by considering the body weight and the number of bacteria in the saliva of periodontitis patients (26–28). At 1 week after the commencement of infection, the diet was changed to CDAHFD60 (#A06071302, Research Diets Inc., New Brunswick, NJ), except for the negative control group that was only administered the vehicle only until the end of the experiment. To analyze the effect of bacteria alone, an experimental group without diet change was also set (Figure 1).

Part of the left and medial lobes were fixed in neutral-buffered formalin. After deparaffinization and rehydration, paraffin-embedded sections (5 μm in thickness) were stained with hematoxylin and eosin (H&E) or were subjected to Masson’s trichrome staining to visualize collagen fibrils.

Liver triglycerides were determined using a Triglyceride Quantification Colorimetric kit (BioVision Inc., Milpitas, CA, USA). Triglyceride levels were expressed as the concentration of triglycerides divided by the total protein concentration. The hepatic hydroxyproline level was photometrically measured using a commercially available kit (Quickzyme Bioscience, Leiden, Netherlands), according to the manufacturer’s instructions.

Feces were collected and stored at −80°C until ready for use. Fecal DNA extraction was performed as described previously (29, 30). DNA was extracted from fecal samples by using lysozyme (Wako Pure Chemical Industries, Osaka, Japan), achromopeptidase (Wako Pure Chemical Industries), and proteinase K (Merck & Co., Inc., Kenilworth, NJ).

Bacterial DNA from fecal samples was amplified by PCR, as described previously (25). Primers (515F and 806R) with the adaptor sequencing for the Illumina MiSeq platform (Illumina, San Diego, CA, USA) were used. PCR amplification using Ex Taq Hot Start Version (Takara Bio, Shiga, Japan) was performed for 25 cycles. Amplicons were purified with AMPure XP (Beckman Coulter, Brea, CA, USA) and sequenced using the Illumina MiSeq platform.

Metagenome shotgun libraries (insert size of 500 bp) were prepared using the TruSeq Nano DNA kit (Illumina) and sequenced on the Illumina NovaSeq platform. After quality filtering, reads mapped to the reference human genome (HG19), and the phiX bacteriophage genome were removed. For each individual, filter-passed NovaSeq reads were assembled using MEGAHIT (v1.2.4). Contigs with a length <500 bp were removed. Gene prediction was performed using the Prodigal software (31). The nucleotide sequences of the predicted genes were clustered using CD-HIT-est v.4.6, with ≥95% sequence identity and ≥90% coverage. The clustered gene sequences were translated into proteins. Functional annotation of the predicted proteins against the KEGG database (release 63) was performed using DIAMOND v0.8, with an E-value of 1e−5. A total of 1 million high-quality reads were mapped onto the nonredundant gene set using bowtie2 to quantify the functional composition of each sample.

Serum levels of ALT and AST were analyzed using commercially available kits (BioVision, Inc., Milpitas, CA, USA), according to the manufacturer’s instructions.

Serum endotoxin levels were determined using a Limulus Amoebocyte Lysate Test (Toxicolor™ LS50M, Seikagaku Co., Tokyo, Japan), according to the manufacturer’s instructions. Serum samples were 1:4 diluted for the assay. Optical densities were measured using an enzyme-linked immunosorbent assay (ELISA) plate reader (Emax Plus, Molecular Devices, San Jose, CA, USA) at 545 nm.

Total RNA from the tissue samples was extracted using the TRIzol reagent (Molecular Research Center) 24 h after the final bacterial or sham administration, and quantified using a NanoDrop 2000 (Thermo Scientific, Wilmington, DE, USA). The total RNA was labeled and then hybridized to an Agilent SurePrint G3 Mouse Gene Expression 8 x 60K mRNA microarray chip (Agilent Technologies). All microarray experiments were conducted in Macrogen, Japan (Kyoto, Japan).

Microarray results were extracted using the Agilent Feature Extraction software v11.0 (Agilent Technologies). Hierarchical cluster analysis was performed using complete linkage and Euclidean distance as a measure of similarity. Gene enrichment and functional annotation analysis for the significant probe list was performed using the GO resource (www.geneontology.org/↗). All data analyses and visualization of differentially expressed genes were conducted using R 3.3.2 (www.r-project.org↗).

Hierarchical cluster analysis was carried out using Ward’s linkage (ward.D2) with Euclidean distance (32), and transcriptomic data were clustered into 21 groups. Among the 21 clusters, we selected seven characteristic clusters with up- or downregulated genes in Pg mice. Enrichment analyses were performed to obtain more information about the biological functions and pathways significantly enriched in up- and downregulated genes by focusing on the GO term biological process and KEGG pathways using the DAVID enrichment analysis system (33). P-values were corrected for multiple testing using the Benjamini–Hochberg method implemented in DAVID (34).

cDNA was synthesized using the Transcriptor Universal cDNA Master (Roche Molecular Systems, Pleasanton, CA, USA). Primers and probes for real-time PCR were purchased from Life Technologies (Waltham, MA, USA). Reactions were carried out in a LightCycler 96 System (Roche) using TaqMan Gene Expression Assays (Life Technologies) as described previously (24). The LightCycler 96 software (Roche) was used to analyze the standards and perform quantification. The relative quantity of each mRNA was normalized to that of glyceraldehyde-3-phosphate dehydrogenase mRNA.

Serum samples diluted to one-sixth of their original concentration in 100 mmol/L potassium phosphate buffer (in deuterium oxide containing 1 mmol/L sodium 2,2-dimethyl-2-silapentane-5-sulfonate, pH = 7) were quantified using an Nuclear Magnetic Resonance (NMR) spectrometer (Bruker AVANCE II 700, Bruker Biospin, Rheinstetten, Germany) as described previously (25, 35). Intact liver samples were placed in zirconia 4 mm diameter zirconia rotors and analyzed by ¹H high-resolution magic angle spinning (hr-MAS) NMR spectroscopy at 500.132 MHz, with a spin rate of 4000 Hz (36). To annotate the signals detected in the ¹H NMR spectra, two-dimensional J-resolved (J-res) NMR measurements (Bruker standard pulse program “jresgpprqf”) and hetero-nuclear single quantum coherence (HSQC) measurements (Bruker standard pulse program “hsqcetgpsisp2.2” were performed as described previously (37–39). The detected signals were annotated using the SpinCouple program (40) (http://dmar.riken.jp/spincouple/↗), and InterAnalysis program (41) (http://dmar.riken.jp/interanalysis/↗) based on HSQC and J-res cross peaks.

Tissue sections above mentioned were deparaffinized with xylene, and rehydrated. This was followed by heat-induced antigen retrieval in BD Retrievagen A (pH 6.0, BD Biosciences, San Diego CA, USA). After blocking in fetal bovine serum (FBS), tissues were incubated with fluorescence-labeled anti-E-cadherin antibody (Alexa FluorTM 594 anti-mouse/human CD324, Biolegend, San Diego, CA, USA). After mounting with VECTASHIELD HardSet Mounting Medium with DAPI (Vector Laboratories, ABurlingame CA, USA), the stained sections were visualized by fluorescence microscopy (Biozero BZ-X710; Keyence Corporation, Osaka, Japan).

HepG2 cells were obtained from the American Type Culture Collection (Manassas, VA) and maintained in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 1% penicillin/streptomycin and 10% FBS in an atmosphere of 5% CO₂ at 37°C. The cells were seeded in a 48-well plate at 3 × 10⁵ cells/well and cultivated for 24 h. Thereafter, the medium was changed to DMEM without FBS and the cells were cultured overnight. The cells were treated with a free fatty acid solution (final concentrations: 0.5 mM oleic acid and 0.25 mM palmitic acid) and fat-free bovine serum albumin for 24 h. The cells were then stimulated with LPS (1 μg/mL P. gingivalis or P. intermedia LPS, and 1 ng/mL E. coli LPS) for 4 h. Total RNA was extracted from the cells as described above and used for quantitative PCR.

Taxonomic assignments and estimation of relative abundance from sequencing data were performed using the analysis pipeline of the QIIME2 version 2020.6.0 (42). Amplicon sequence variants (ASVs) were inferred from the denoised reads using DADA2 (43) implemented in QIIME2. The ASV taxonomy was assigned based on a comparison with the SILVA version 138 (44). β-Diversity was calculated using weighted UniFrac distances based on the operational taxonomic unit distribution across samples and visualized by principal coordinate analysis (PCoA).

From the results of metagenomic analysis at TP2 and TP3, KEGG orthologies (KOs) with significant differences determined by the t-tests, which were carried out by using R, between two groups, except for RC mice, were extracted under the condition of P < 0.01. The KOs were mapped to the reference pathway in the KEGG database. Significant pathways were enriched with Fisher’s exact probability tests using the number of KOs mapped in each pathway. The p-values were adjusted by using the Benjamini–Hochberg method, which are used as q-values in the Figure 6. The quantified metabolome data were normalized using an autoscaling method and statistically analyzed using PCA.

Statistical analyses were performed using the GraphPad Prism version 9 (GraphPad Software, Inc., La Jolla, CA, USA) and R (version 4.0.4.). Randomization or blinding was not performed in the present study. All data are expressed as the mean ± standard error of the mean. Statistical analyses were performed using one-way analysis of variance with Tukey’s correction. Analysis of similarity was performed to identify differences in bacterial community compositions and PERMANOVA (PERmutational Multivariant Analysis Of Variance) was used for comparison of microbes between groups. Statistical significance was set at P <0.05.

Oral administration of bacteria did not affect body weight until commencement of the diet change, and no difference was observed among the vehicle-, Actinomyces naeslundii-, Veillonella rogosae-, P. gingivalis- and P. intermedia-administered mice (hereafter referred to as Sham, An, Vr, Pg, and Pi mice, respectively). After changing the diet from regular chow (RC) to choline-deficient, l-amino acid-defined, high-fat diet with 60 kcal% fat and 0.1% methionine (CDAHFD60) (CDAHFD60; HFD), a gain in body weight was suppressed in all experimental groups in comparison with RC-fed mice. In contrast, the liver-to-body weight ratio was significantly increased in Sham, An, and Pg mice compared with that in RC-fed mice (Figure 2A). HFD feeding promoted hepatic steatosis in all experimental groups, and the degree of steatosis was greater with P. intermedia administration, and it was further aggravated by P. gingivalis administration (Figure 2B). Similarly, the degree of fibrosis was progressively exacerbated with increasing bacterial burden (sham < P. intermedia < P. gingivalis) (Figure 2C). In contrast, no additional histological changes were observed with A. naeslundii and V. rogosae administration compared with those in Sham mice. Bacterial administration induced minimal histological changes in the liver of RC-fed mice except for Pg mice, in which slight steatosis was observed (Supplementary Figure 1). Therefore, further analyses were focused on Pg and Pi mice.

Under these conditions, the content of hepatic hydroxyproline increased with the increasing bacterial burden (Figure 2D). Although triglyceride content and aspartate transaminase (AST) and alanine transaminase (ALT) activities were significantly higher in HFD-fed groups, bacterial administration had no effect on them (Supplementary Figure 2).

Dysregulation of gut barrier function and subsequent endotoxemia are major contributors to NAFLD. Therefore, we analyzed whether the barrier function was altered in Pg mice, and if the barrier function was disorganized, whether endotoxemia was induced. As shown in Figure 3A, the expression of Tjp1, encoding tight junction protein, tended to be lower in Pg mice compared with other groups. Moreover, a decrease in the expression of E-cadherin was observed in the colon of Pg mice (Figure 3B). Additionally, the serum endotoxin level was increased with the increasing bacterial burden. The level was significantly higher in Pi mice compared with that in RC-fed and Sham mice, and it was further elevated in Pg mice (Figure 3C).

The global differences among experimental groups were evident from baseline (time point 1; TP1) to TP2, after 1 week of bacterial administration, as well as after additional HFD-feeding (TP3), according to the analyses of α- (Figure 4A) and β-diversities (Figure 4B).

At TP2, the proportion of the phylum Firmicutes was significantly lower in Pg mice than that in Sham mice (Figure 5A). However, the proportion of phylum Bacteridota (previously known as Bacteroidetes) did not differ among the experimental groups. At the genus level, the proportion of Eubacterium fissicatena group was significantly higher in Pg mice than in Sham mice. In contrast, the proportion of Lactobacillus was significantly lower and tended to be lower in Pg mice compared to RC-fed and Sham mice, respectively (Figure 5B).

HFD feeding had a strong effect on the gut microbiota composition. Despite the considerable effect of diet, bacterial administration induced a substantial change in the gut microbiota composition. Although no difference was observed in the phylum Firmicutes, HFD feeding significantly decreased the proportion of the phylum Bacteridota (Figure 5C). At the genus level, significantly elevated proportions of Atopostipes and Colidextribacter were observed in Pg mice compared with those in RC and Sham mice. The proportions of Lactobacillus and Muribaculaceae were significantly lower in HFD-fed mice and there was no effect of bacterial administration (Figure 5D). Although HFD feeding increased the abundance of many genera, there were some notable patterns of relative abundance among the different experimental groups. While P. intermedia administration further increased the relative abundance of some genera, it suppressed the effect of HFD feeding in some genera (Supplementary Figure 3).

Thus, hepatic pathology observed in Pg mice appeared to be directly related to gut dysbiosis. Consistent with these findings, the expression profile of genes in the intestine showed distinct patterns when compared between any two experimental groups (Supplementary Figure 4), suggesting that diet as well as the administered bacteria have a characteristic effect on gene expression, possibly through changes in the gut microbiota composition.

Enrichment analysis was applied to the relative abundance of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways via metagenomic analysis to study the function of gut microbiota in different groups at different time points. At TP2, in Pg and Pi mice, the enriched KEGG pathways overlapped with elevated amino acid metabolism, particularly, with respect to phenylalanine, tyrosine, and tryptophan biosynthesis. In Pi mice, the lipopolysaccharide biosynthesis pathway was also notably overrepresented. A total of eight enriched KEGG pathways, including the NOD-like receptor signaling pathway, were observed in Pg mice relative to the pathways observed in Pi mice (Figure 6A). Although the effect of bacterial administration on the gut metabolic pathways was diminished after HFD feeding (the number of pathways enriched at TP2 decreased compared to that at TP2), the difference between Pg mice and Pi mice was even more pronounced. Similar to TP2, amino acid metabolic pathways were the characteristic pathways when the two groups were compared. In addition, the pathways for fatty acid degradation and the two-component system were enriched in Pg mice relative to those in Pi mice (Figure 6B).

We further analyzed liver tissue and serum samples to assess whether the murine metabolome was perturbed by bacterial administration, in addition to being affected by diet. We performed an untargeted analysis of all the data acquired to examine a wider pool of metabolites. Principal component analysis (PCA) was conducted using nuclear magnetic resonance (NMR)-derived data to obtain an overview of the differences among the groups. This analysis not only differentiated HFD-fed mice from RC-fed mice, but also highlighted the effect of P. gingivalis administration among HFD-fed mice in the liver tissue and serum metabolites (Supplementary Figure 5A and Figure 7A), which was consistent with the histological findings. Furthermore, the machine-learning model allowed a discrete classification within the four groups in both tissue and serum metabolites; maltose, choline, a fatty acid (FA.CH3.n.3), choline metabolite (CHONCH3), and methionine were among the characteristic metabolites in the tissue (Supplementary Figure 5B), whereas metabolites such as ROI.38 [which has been suggested to be a sugar-phosphate (sugar-P)], tyrosine, and choline were revealed to be important in the serum (Figure 7B). The tissue profile of short-chain fatty acids (SCFAs) showed no difference among the groups (Supplementary Figure 6A). In contrast, the levels of acetate and citrate were significantly higher in HFD-fed mice compared with those in RC-fed mice and tended to increase with an increasing bacterial burden, with the highest level in the serum of Pg mice. However, there was no difference in the level of the other SCFAs (Supplementary Figure 7A). With respect to the amino acid levels in the tissue, the level of lysine was significantly decreased in HFD-fed mice, with a further decrease in Pg mice compared with Sham mice. A similar trend was observed for threonine (Supplementary Figure 6B). The level of ChoNCH3, a choline metabolite, showed a pattern similar to that of choline (Supplementary Figure 5C). The levels of lipid metabolites were increased in HFD-fed mice; however, bacterial administration had a minimal effect on these levels (Supplementary Figure 5C). In the sera, HFD feeding significantly elevated the serum levels of the annotated signals, except for isoleucine, glutamine, and phenylalanine (Supplementary Figure 7B). Notably, tyrosine levels were significantly elevated in Pg mice compared with those in Sham mice. Similarly, sugar-P levels were significantly higher in Pg mice compared with those in the other groups. The choline levels were drastically lower in Pg mice compared with those in the other groups (Figure 7C). The levels of other annotated molecules were also elevated in HFD-fed mice, except those of allantoin, which increased with increasing nutrition and bacterial burden, and there was a significant difference in these levels between RC-fed and Pg mice (Supplementary Figure 7C).

Although the effect of HFD feeding on the liver was robust, it was evident that the administration of both bacteria induced additional and substantial changes in the expression profile of genes. Moreover, the expression profile of genes in the liver was clearly distinct between Pi and Pg mice (Supplementary Figure 8). Therefore, we evaluated the expression of genes in the liver using functional enrichment analysis.

Based on hierarchical clustering, we extracted seven clusters in which genes in clusters 1–3 were downregulated, whereas those in clusters 4–6 were upregulated in Pg mice compared with the other groups. Cluster 7 included genes, the expression levels of which were higher in Pg mice than those in RC-fed mice, but lower compared with those in Sham and Pi mice (Supplementary Figure 9). The genes in these clusters are listed in Supplementary Table 1.

In cluster 1, significantly downregulated genes were mostly annotated as those involved in the biosynthesis and metabolic processes of lipids, organic acids, oxoacids, steroids, and fatty acids. No gene ontology (GO) terms were enriched in clusters 2 and 3. In cluster 4, genes involved in the cell cycle, cell death, nuclear division, DNA replication, responses to oxidative and ER stress, and regulation of intrinsic apoptosis were upregulated. In cluster 5, genes involved in the inflammatory response, such as in the regulation of leukocyte migration, response to lipopolysaccharide (LPS), cellular response to inflammatory cytokines, and nuclear factor-κB (NF-κB) pathways, were significantly enriched. In cluster 7, genes involved in the biological processes for circadian regulation of gene expression were significantly enriched (). 1

The significantly enriched KEGG pathways in clusters 1, 2, and 4 are shown in Supplementary Table 3. These included various pathways for steroid hormones, retinol, primary bile acids, arachidonic acid, amino sugars, and nucleotide sugars. Some metabolic pathways were consistent with the enriched Gene Ontology (GO) terms in cluster 1. The top 15 significant GO terms and the complete list of enriched KEGG pathways in each cluster are shown in Figure 8.

Progression from simple steatosis to steatohepatitis, fibrosis, and hepatocellular carcinoma is believed to be mediated by multiple parallel factors, including inflammation, ER stress, lipotoxicity, and altered circadian rhythms. The GO terms and KEGG pathways implicated in these events were found to be significantly enriched in the clusters (Supplementary Tables 2, 3). To confirm the data of DNA microarray analysis, quantitative real-time PCR was conducted for representative genes involved in the pathological mechanisms of NAFLD. Genes involved in the inflammatory responses were upregulated by HFD feeding, irrespective of the bacterial administration. The expression of Tsc22d3, which mediates the anti-inflammatory response, was only elevated in Pg mice, reflecting a weaker elevation of inflammatory genes in this group (Supplementary Figure 10A).

Fibrosis- and ER stress-related genes were enriched and classified in cluster 4. Whereas Col1a1 and Timp1 were upregulated by HFD feeding, Ctgf was upregulated specifically by P. gingivalis administration (Supplementary Figure 10B). ER stress-related genes Chop/Ddit3 and Ddit4 were also elevated by P. gingivalis. Consistent with these findings, Fgf21 and Trib3, both of which are induced by ER stress and have been reported to be elevated in NAFLD, were significantly upregulated in Pg mice, but not in Pi mice (Supplementary Figure 10C).

The GO terms associated with the cell cycle process, which are potentially implicated in carcinogenesis and end-stage NAFLD, were also enriched in cluster 4. The expression of Hnf6 and Hhex was upregulated in Pi mice, but not in Pg mice (Supplementary Figure 10D). The expression of Hnf6 in hepatocellular carcinoma (HCC) cells is negatively associated with their malignancy. Delivery of Hhex into a hepatoma cell line has been reported to decrease the expression of several proto-oncogenes and to increases the expression of some tumor suppressor genes (45).

In cluster 7, GO terms associated with controlling the circadian rhythm were enriched. Fluctuations in the circadian rhythm affect metabolism and alter the expression of liver clock genes in NAFLD pathology. In this connection, the clock gene, Per1, was significantly elevated in Pg mice. Pg and Pi mice demonstrated contrasting expression patterns of Bmal1 and Dbp. (Supplementary Figure 10E).

These distinctive expression profiles of hepatic genes in RC-fed, Pi, and Pg mice may not be due to the direct effect of administered bacteria, but rather due to differences in the gut microbiota compositions triggered bacterial administration. In support of this assumption, the gene expression of ER stress-related genes in HepG2 cells stimulated with LPS from P. intermedia, P. gingivalis, and E. coli was different from that in the liver of each experimental group (Supplementary Figure 11).

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/↗, GSE136937↗.