Metabolic impact of persistent organic pollutants on gut microbiota

The physiologic and community structure effect of POPs on cecal bacterial mixtures in vitro

Cecal bacterial physiologic status after short-term exposure of POPs in vitro was assessed using a flow cytometry approach (Figure 2). A significant dose-dependent decrease in bacteria with high nucleic acid (HNA) content was observed following exposure to all four POPs (Figure 2(a)). Short-term incubation of POPs resulted in no significant changes in SybrGreen, propidium iodide (Pi), and bis-(1,3-dibutylbarbituric acid) trimethine oxonol (DiBAC₄) stained cells but significantly decreased carboxyfluorescein diacetate (CFDA)-stained bacteria in a dose-dependent manner (Figure 2(a-b)), suggesting a rapid decrease in microbial metabolic activity without significant cell damage. Moreover, changes in the cecal community structure were observed with TCDD exposure via Generalized UniFrac analysis (Supplementary Figure 1a). Short-term TCDD exposure had significantly increased numbers of bacteria from the genus Lactobacillus, Roseburia, and Oscillibacter but decrease in the genus Bacteroides (Supplementary Figure 1b).

The metabolic effect of POPs on cecal bacterial mixtures in vitro

To further explore the influence of POPs on the cecal bacteria, hydrophilic metabolite and lipid profiling were performed using ¹H NMR-based metabolomics (Figure 3). Principal component analysis (PCA) score plots showed distinct clustering of hydrophilic metabolite and lipid profiling obtained from cecal bacteria with two doses of POPs compared to vehicle (Figure 3(a) and 3(c)). PCB 153 exposure (6 µM) resulted in significantly lower levels of nucleic acids, acetate, formate, methionine, and aspartate, but significantly higher levels of glucose, lipids, lactate, butyrate, succinate, alanine, lysine, tyrosine, phenylalanine, and branched-chain amino acids (BCAAs) in cecal bacteria (Figure 3(b)). In contrast, TCDF, TCDD, and PCB 126 did not exhibit a pronounced effect on microbial hydrophilic metabolites (Figure 3(b)). Significant changes in fatty acid metabolism were observed with PCB 153 exposure (6 µM) compared to other POPs (Figure 3(d)). PCB 153 exposure (6 µM) resulted in significantly lower unsaturated fatty acids (UFA) and monounsaturated fatty acids (MUFA) but higher saturated fatty acids (SFA) and the ratio of SFA to UFA; whereas exposure to other POPs resulted in significantly lower SFA but higher UFA and MUFA (Figure 3(d)).

UPLC-MS analysis was performed to identify altered microbial lipids following POP exposure (Figure 4 and Supplementary Table 1–2). Significant increases in bacterial membrane lipids including fatty acid, fatty acid esters of hydroxy fatty acid (FAHFA), phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylinositol (PI), phosphatidylglycerol (PG), sphingomyelin (SM), diacylg lycerol (DAG), triacylglycerol (TAG), lysophosphati dylcholine (LPC), and lysophosphatidylethanolamine (LPE) were observed following two doses (0.6 µM and 6 µM) of PCB 153 and a higher dose of TCDD (0.6 µM) and PCB 126 (6 µM) exposure (Figure 4(a-b) and Supplementary Table 1–2). Together, these results indicate that short-term incubation of POPs directly and rapidly affect cecal bacterial global metabolism that is associated with significant decreases in microbial metabolic activity.

The transcriptional effect of POPs on cecal bacterial mixtures in vitro

We investigated the bacterial functional response to POPs by analyzing the metatranscriptome from cecal bacteria with TCDF (6 µM), TCDD (0.6 µM), and PCBs (6 µM) exposure compared to vehicle using RNA-seq. PCA showed distinct separations between cecal bacteria with vehicle and POP exposure with the greatest differences occurring between PCB 153 and other POP (p-value = 0.002 of PCA models from each POP compared to vehicle) (Figure 5(a)). The number of differentially expressed gene orthologs differed among samples exposed to different POPs (Figure 5(b)). The number of upregulated gene orthologs varied from 60 to 212 for TCDD (0.6 µM) and PCB 153 exposure (6 µM), while the number of downregulated gene orthologs varied from 19 to 84 for PCB 126 (6 µM) and TCDD exposure (0.6 µM) (Figure 5(b)). Among the differentially expressed gene orthologs, most were found to be specifically changed only by one POP (Figure 5(c)). We identified five gene orthologs that were significantly regulated in response to all four POPs corresponding to phosphate butyryltransferase (ptb), small acid-soluble spore protein D (sspD), putative pyrimidine permease RutG (rutG), 2-oxoisovalerate dehydrogenase E1 component (bkdA), and transposase IS5 family (Figure 5(d)), which are associated with butyrate metabolism,²¹ BCAAs degradation,²² and pyrimidine degradation²³ as well as DNA cleavage and modification,²⁴ respectively.

Analysis of the KEGG module and pathway enrichment confirmed and extended these trends (Table 1). Pathways including tricarboxylic acid (TCA) cycle, glycolysis, amino acid biosynthesis, nucleotide metabolism, and fermentation pathways were integrated and visualized using metatranscriptomics data combined with metabolomics data (Figure 5(e-f) and Supplementary Fig. 2–4). Although each POP induced specific transcriptomic responses, some general trends were shared, such as transcript levels related to the TCA cycle, carbon metabolism, carbon fixation pathways, and protein export (Table 1). For example, the TCA cycle and glycolysis pathways were significant enriched by exposure to most POPs, especially PCB 153 (6 µM), supported by significantly higher levels of gene orthologs and intermediates in pathways (Figure 5(e-f)). Of particular note, PCB 153 exposure (6 µM) resulted in distinct changes in bacterial transcript levels compared to other POPs, such as fatty acid biosynthesis, lipopolysaccharide (LPS) biosynthesis, vitamin B6 metabolism, mismatch repair, and base excision repair (Table 1).

To further explore the transcriptomic responses of specific bacteria following POP exposure, five representative bacteria including Bacteroides, Clostridium, Lactobacillus, Bifidobacterium, and Fusobacterium were mapped to their reference sequences (Supplementary Table 3). The differentially expressed gene orthologs induced by four POPs compared to a vehicle from these five genera are listed in Supplementary Tables 4–8. In total, we identified 125 gene orthologs differentially expressed after TCDF exposure (6 µM), 96 gene orthologs differentially expressed after TCDD exposure (0.6 µM), 98 gene orthologs differentially expressed after PCB 126 exposure (6 µM), and 234 gene orthologs differentially expressed after TCDF exposure (6 µM) in Bacteroides (Supplementary Table 4). Bifidobacterium and Fusobacterium have less differentially expressed gene orthologs by POPs exposure, which could be due to the reduced mRNA that mapped to their reference sequences (Supplementary Tables 3 and 7–8). Consistent with isolated cecal bacteria results, PCB 153 exposure resulted in a greater impact on these five individual bacteria than other POPs at the transcriptional levels (Supplementary Tables 3–8).

The physiologic and metabolic effect of POPs on individual bacteria species in vitro

Having defined the transcriptomic responses of specific bacteria with POPs exposure, we sought to explore the physiologic and metabolic effect of POPs on individual bacteria. Gram-negative and Gram-positive bacteria are have different cell wall structures that play a variety of protective and adaptive roles.^25,26 In order to explore the specific effect of POPs on Gram-negative and Gram-positive bacteria, we combined growth rate measurements (OD₆₀₀) and flow cytometry with UPLC-MS-based metabolomics to investigate the responses of POPs on six individual bacteria species. Six bacterial species were used including three Gram-positive bacteria (Lactobacillus paracasei, Bifidobacterium longum, and Clostridium ramosum) and three Gram-negative bacteria (Fusobacteria nucleatum, Bacteriodes fragilis (638 R), and Bacteriodes fragilis (ATCC®25282^TM)), which are found in the gut of healthy individuals and play a key role in host metabolism.^27,28

We combined flow cytometry with optical density (OD₆₀₀) over 18 h to monitor the impact of POPs on the growth and physiologic status of six bacteria species (Figure 6 and Supplementary Figure 5–6). No significant effects were observed in the growth rate and physiologic status in most bacteria species including Bacteriodes fragilis (ATCC®25282^TM), Bacteriodes fragilis (638 R), Fusobacteria nucleatum, and Clostridium ramosum following exposure to two doses of four POPs for 18 h (Supplementary Figure 5–6), indicating no significant cell damage following POP exposure to those bacteria species. The significant decreases in the growth of Bifidobacterium longum were observed with a higher dose of TCDF (6 µM) and TCDD (0.6 µM) and two doses of PCB 126 and PCB 153 exposure, while subtle decreases were observed with a lower dose of TCDF (0.6 µM) and TCDD (0.06 µM) exposure (Figure 6(a)). Consistent with the OD₆₀₀ results, two doses of PCB 153 (0.6 µM and 6 µM) and a higher dose of PCB 126 (6 µM) exposure significantly increased the proportions of damaged cells (Pi and DiBAC₄-stained cells) and decreased CFDA-stained bacteria (Figure 6(c)), suggesting dramatic changes to the physiological status of Bifidobacterium longum after 18 h POPs exposure. POPs exposure did not result in dramatic changes in the growth of Lactobacillus paracasei but a subtle increase in the growth of Lactobacillus paracasei, which is indicated by the increased growth rate and a significant decrease in Pi-stained bacteria with POPs exposure (Figure 6(b) and 6(d)).

Having determined the significant changes in lipid profiling in cecal bacteria with POP exposure, we sought to investigate the influence of POPs on the lipid metabolism of these six bacteria species (Figure 7 and Supplementary Figure 7–18). To visualize the patterns of lipid metabolism of six bacteria species following POP exposure, the log₂ fold change values of lipid profiles obtained from UPLC-MS were assembled using self-organizing maps (SOMs) and projected onto suprahexagonal landscapes (Figure 7 and Supplementary Figure 7–18). These maps display the “metabolic fingerprint” for six bacteria species following POP exposure and provide a means to visually interpret complex lipid changes. The six bacteria species showed extreme variation in lipid composition in response to TCDD, TCDF, or PCBs exposure in both positive and negative modes, with increases of most lipids in Lactobacillus paracasei with POPs exposure and decreases of most lipids in Clostridium ramosum with POPs exposure (Figure 7, Supplementary Figs. 8–9, and Supplementary Figs. 14–15). Notably, two strains of Bacteriodes fragilis also showed different responses to POP exposure, with Bacteriodes fragilis (ATCC®25282^TM) being more sensitive (Figure 7, Supplementary Figs. 11–12, and Supplementary Figs. 17–18). Moreover, most individual bacteria species, such as Bifidobacterium longum, Clostridium ramosum, Fusobacteria nucleatum, and Bacteriodes fragilis (ATCC®25282^TM) are more sensitive to (non-coplanar) PCB-153 than (coplanar) PCB-126 (Figure 7). It is interesting to note that limited significant effects in hydrophilic metabolites were observed following exposure to two doses of POPs in Lactobacillus paracasei and Bifidobacterium longum (Supplementary Figs. 19–20 and Supplementary Tables 9–10).

Chemicals

POPs including TCDD, TCDF, PCB-126, and PCB-153 (Figure 1(c)) were purchased from Cambridge Isotope Laboratories, Inc. (Tewksbury, MA) and suspended in dimethyl sulfoxide (Sigma-Aldrich, St. Louis, MO). The fluorescent dyes including SybrGreen, propidium iodide (Pi), carboxyfluorescein diacetate (CFDA), and bis-(1,3-dibutylbarbituric acid) trimethine oxonol (DiBAC₄) were ordered from Sigma-Aldrich (St. Louis, MO) and Invitrogen (Carlsbad, CA).

In vitro bacterial culture and incubation

Cecal microbiota incubation studies were modified from previously described protocols.³³ Briefly, cecal contents were isolated from 7-week-old C57BL/6 J male wild type mice and diluted with brain heart infusion broth (1 g in 10 ml). The cecal suspensions were incubated with two doses of TCDD (high: 0.6 µM and low: 0.06 µM), TCDF (high: 6 µM and low: 0.6 µM), and PCBs (high: 6 µM and low: 0.6 µM) at 37°C for 4 h. The high dose (6 µM) of TCDF and PCBs is equivalent to the dose of TCDF that we have previously published in vivo studies.^13,20 After incubation, bacterial samples were washed and stained for flow cytometry analysis or kept at −80°C for other analyses.

Six individual bacteria species including Bifidobacterium longum (ATCC®15707^TM), Lactobacillus paracasei (ATCC®25303^TM), Fusobacteria nucleatum subsp. nucleatum (ATCC®25586^TM), Clostridium ramosum (ATCC®25582^TM), Bacteriodes fragilis (ATCC®25282^TM), and Bacteriodes fragilis (638 R) were cultured and treated with two doses of TCDD (high: 0.6 µM and low: 0.06 µM), TCDF (high: 6 µM and low: 0.6 µM), and PCBs (high: 6 µM and low: 0.6 µM) at 37°C for 18 h (for more detailed methods about bacteria culture, see Supplemental Material). The effects of POPs on the growth rate of different bacteria in vitro were measured by optical density (OD₆₀₀) using a Multiskan Sky Microplate spectrophotometer (Thermo Fisher Scientific, Waltham, MA). Flow cytometry and metabolite analysis were performed after 18 h of POPs incubation. All experiments were performed in a monitored anaerobic chamber (Coy Laboratory Products, 95% N₂, 5% H₂) and repeated at least three times.

Flow cytometry analysis

After incubation, the bacterial mixtures and individual bacteria species were washed twice and resuspended with 1X reduced PBS (1 mg/ml L-cysteine). Bacterial physiology was assessed by four fluorescent dyes: SybrGreen that stains all nucleic acids regardless of cellular damage, Pi that stains only dead or damaged cells, CFDA that stains for enzymatic/metabolic activity, and DiBAC₄ that stains damaged bacteria.^33,49 All cytometric analyses were performed using a BD Accuri™ C6 plus flow cytometer (BD Biosciences, San Jose, CA) and data were analyzed with FlowJo V10 software (FlowJo, LLC). The percentages of Pi, CFDA, DiBAC₄, low nucleic acid (LNA), and high nucleic acid (HNA) stained cells were calculated relative to the total bacterial counts obtained by SybrGreen staining.

¹H NMR-based metabolomics experiments

The hydrophilic metabolites and lipids from 1 ml of bacteria were extracted twice with 1 ml pre-cooled methanol/H₂O (2:1, v/v) or chloroform/methanol (2:1, v/v), followed by three freeze-thaw cycles. After evaporation, the extracts were reconstituted in 500 µl of 0.1 M PBS containing 100% D₂O and 0.005% (v/v) 3-(trimethysilyl) [2,2,3,3–²H₄] propionate (TSP) (hydrophilic metabolites) or deuterated chloroform containing 0.03% (v/v) tetramethylsilane (TMS) (lipids) and analyzed with a Bruker Avance NEO 600 MHz spectrometer equipped with an inverse cryogenic probe (Bruker Biospin, Germany) at 298 K. A typical 1D NMR spectrum named NOESYPR1D was acquired for both hydrophilic and lipid extracts. Metabolite assignments were carried out on the basis of a set of 2D NMR spectra and published results.^19,49

All ¹H NMR spectra were phase- and baseline-corrected and referenced to TSP (δ = 0.0) for hydrophilic metabolites or TMS (δ = 0.0) for lipids using TopSpin 3.6 (Bruker Biospin). The spectra were integrated into 0.004 ppm-width buckets using the AMIX 3.8 (Bruker Biospin). Principal component analysis (PCA) was performed using the SIMCA 13.0.3 (Umetrics, Umea, Sweden). Heatmaps were plotted using RStudio (pheatmap). The relative levels for representative fatty acids including the molar percentages of unsaturated fatty acids (UFA%), saturated fatty acids (SFA%), polyunsaturated fatty acids (PUFA%), monounsaturated fatty acids (MUFA%), and SFA-to-UFA ratio (SFA/UFA) were also calculated as previously reported.⁵⁰ These calculations were based on the spectral integral areas (Supplementary Fig. 21) for – CH= CH – (from UFA, δ 5.38), – CH₃ (from all fatty acids, δ 0.84), = CH–CH₂–CH = (from PUFA, δ 2.76) taking into consideration of proton numbers; the signal area for SFA and MUFA was calculated by subtracting that for UFA from 1 and PUFA from UFA.

LC-MS based metabolomics experiments

Bacteria species samples (0.6 ml) were extracted with 1 ml of pre-cool chloroform/methanol (2:1, v/v). After homogenization, the samples were added 250 µl of HPLC water and vortexed, followed by centrifugation (22,000 rpm, 4°C for 10 min). The top phase was collected for hydrophilic metabolites and the upper phase was collected for lipid analysis. Two phases were dried in a vacuum and reconstituted in 60 μl of 3% methanol containing 1 µM chlorpropamide (hydrophilic metabolites) or 100 μl of isopropanol/acetonitrile/H₂O (50:25:25, v/v/v) containing 1 µM 15:0–18:1-d7-PC (lipids). After centrifugation, supernatants were transferred to autosampler vials for LC-MS analysis.

Hydrophilic metabolite analysis was performed with a Dionex Ultimate 3000 quaternary HPLC system connected to Exactive^TM Plus Orbitrap mass spectrometer (Thermo Fisher Scientific, Waltham, MA) with a Hydro-RP C18 column (2.1 × 100 mm × 2.5 µm particle size; Phenomenex, Torrance, CA). Samples for lipid profiling were separated by reverse-phase HPLC using a Vanquish UHPLC system (Thermo Fisher Scientific, Waltham, MA) with a Waters (Milford, MA) CSH C18 column (1.0 × 150 mm × 1.7 um particle size). The eluate was delivered into an Orbitrap Fusion Lumos Tribrid™ mass spectrometer using a H-ESI™ ion source (all Thermo Fisher Scientific) (for more detailed methods about LC-MS, see Supplemental Material). LC-MS data were analyzed by the open-source software MS-DIAL.⁵¹ MetaMapp network analysis for hydrophilic metabolite was performed to visualize and integrate biochemical and chemical similarities as previously described.⁵² The colors of each node such as red and blue represent up- or down-regulation compared to the vehicle. The size and shape of nodes represent fold-change and class of metabolites. The log2 fold changes of lipids from each POPs exposure compared to the vehicle group were trained a self-organizing map (SOM) and projects onto a two-dimensional suprahexagon using the Suprahex package for R/Bioconductor.⁵³ These suprahexagons arrange related metabolites within nodes or small hexagons, that are further arranged based on vector weight. The map preserves the input information and provides the underlying cluster details such as the dimensionality, distribution, distance, clusters, and identity of metabolites.

Metatranscriptomic analysis

The cecal bacterial metatranscriptomic analysis was done as previously described²⁰ (Figure 1(b)). Briefly, total RNA was isolated from 1 ml of cecal bacterial mixtures by 1 ml of Trizol and purified with RNeasy® Mini Kit twice (QIAGEN, Hilden, Germany). The total RNA was measured and checked by an Agilent Bioanalyzer. The 16S and 23S rRNA fractions were removed from total RNA using the RiboMinus^TM Bacteria Transcriptome Isolation kit (Invitrogen, Carlsbad, CA) and purified with RNeasy® MinElute® Cleanup Kit (QIAGEN, Hilden, Germany). The deletion of 16S and 23S rRNA was checked again with the Agilent Bioanalyzer. The rRNA depleted bacterial RNA samples were submitted to the Pennsylvania State University Genomic Sequencing core for sequencing. TruSeq Stranded mRNA kit was used to make a uniquely indexed library from each sample and then sequenced on Illumina NextSeq High-Output 150 nt single read sequencing run.

The obtained sequences were subjected to quality check using FastQC with a minimum Phred score of 30. Low-quality sequences and adapters were removed using Cutadapt⁵⁴ with a quality cutoff of Phred score 30 and adapter sequences recommended for TruSeq Stranded mRNA kits. The quality-controlled trimmed reads were then subjected to ribosomal RNA (rRNA) filtering using SortMeRNA.⁵⁵ Using the Silva database,⁵⁶ all the sequences identified as rRNA were removed to obtain a pool of non-ribosomal rRNA. The non-ribosomal rRNA were then subjected to functional annotation using FMAP.⁵⁷ Briefly, the non-ribosomal rRNA was mapped to the Uniref90 protein database using Diamond.⁵⁸ The uniport ids were mapped to KEGG Orthologs and the abundances of KEGG Orthologs were calculated using FMAP with default parameters. The counts were normalized to the total number of reads in each sample. Exploratory analysis was done using MicrobiomeAnalyst.⁵⁹ PerMANOVA (adonis function in vegan R package) was used to evaluate the effect of treatment on Bray-Curtis distance matrices with 999 random permutations. Each treatment group (TCDD, TCDF, PCB-126, and PCB-153) was individually compared to the vehicle group. Differential abundance analysis was tested using Mann-Whitney-Wilcoxon tests. Multiple test correction was performed by calculating False Discovery Rate-adjusted p-value using the Benjamini-Hochberg method with a threshold of p-value <0.05 and p.adjust <0.5. The significantly different KEGG Orthologs were selected for visualizations and further analysis.

Investigations into the effect of POPs on specific bacteria in the gut microbiome was carried out using a custom pipeline. Sequences pertaining to five representative genera of gut microbiota: Bacteroides, Clostridium, Lactobacillus, Bifidobacterium, and Fusobacterium were downloaded from NCBI Reference sequence databases. Filtered mRNA were mapped to these reference sequences using Diamond⁵⁸ using default parameters. The abundance of each gene was calculated by counting the number of hits in the reference database. The counts were normalized to the total number of reads in each sample. For each bacterium, differentially expressed gene orthologs for each POP were calculated using the student t-test. Multiple test corrections were performed by calculating False Discovery Rate-adjusted p-value using the Benjamini-Hochberg method with a threshold of p-value <0.05 and p.adjust <0.05.

Metatranscriptomic and metabolomics data integration and visualization

Metabolomics data and metatranscriptomics data were visualized using MetaboMaps.⁶⁰ Pathways were condensed into a custom figure. The enzymes were mapped to gene orthologs. Log2Fold changes in metabolite and gene orthologs were plotted on the graph as circle and box heatmaps.

Sanger sequencing and 16S rRNA gene sequencing analysis

Bacteria DNA from bacterial mixtures and individual bacteria species was extracted using the E.Z.N.A. stool DNA kit (Omega Bio-Tek Inc., Norcross, GA). Bacterial species used throughout the study were confirmed by Sanger sequencing as previously described⁶¹ and the result was analyzed with BLAST (Supplementary Table 11). The DNA from cecal bacterial mixtures was analyzed using 16S rRNA gene sequencing as described.⁶² Briefly, the extracted bacterial DNA was amplified using the V4V4 primer set (515 F and 806 R). The product quality was assessed using 1% agarose gel electrophoresis and a DNA 7500 LabChip on the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clare, CA). The sequencing was performed on the Illumina MiSeq platform at Penn State Genomics Core Facility (University Park, PA) and analyzed using the mothur platform.⁶³ The Generalized Unifrac analysis was performed using a phylogenic tree and operational taxonomic unit (OTU) table by R studio.⁶³

Statistics

All data values are expressed as mean ± standard deviation (SD). Graphical illustrations and statistical analyses were performed using GraphPad Prism 6.0 (GraphPad, San Diego, CA). One-way analysis of variance (ANOVA) followed by Dunnett’s multiple comparisons was used.