Background. Recent evidence suggests that multigenic and complex environmentally modulated diseases result from early life exposure to toxicants at least partly via gut microbial influences. Environmental toxicants, polybrominated diphenyl ethers (PBDEs), and polychlorinated biphenyls (PCBs) are breast milk-enriched persistent organic pollutants (POPs) and thus remain a continuing threat to human health despite being banned from production. Recent findings focused on the liver developmental reprogramming capabilities from neonatal BPA exposure; however, little is known on how PBDEs and PCBs regulate the liver transcriptome with respect to the gut microbiome.

Objectives. We investigated whether the gut microbiome can be persistently reprogrammed with the liver following neonatal exposure to POPs, and whether microbial biomarkers associated with disease-prone changes in the hepatic epigenetic and transcriptomic landscape in adulthood.

Methods. C57BL/6 male and female mouse pups were orally administered vehicle, bisphenol A (BPA), BDE-99 (a breast milk-enriched PBDE congener), or the Fox River PCB mixture (an environmentally relevant PCB mixture), between postnatal day (PND) 2 to 4, once daily for three consecutive days. Tissues were collected at PND5 and PND60 for 16S rDNA sequencing and targeted metabolomics.

Results. Neonatal exposure to BDE-99, followed by BPA and PCB, produced the greatest persistent changes in the adult hepatic transcriptome, including an inflammation and cancer-prone transcriptomic signature. BDE-99 exposure resulted in a persistent increase in Akkermansia muciniphila throughout the intestinal sections and feces. We observed persistent increases in acetate and succinate, metabolites A. muciniphila is able to produce. Correspondingly, liver H3K4me1 and H3K27 acetylation were enriched around the loci encoding liver cancer-related genes following neonatal BDE-99 exposure.

Conclusion. Similar to BPA, early life exposure to BDE-99 also produced a cancer-prone hepatic transcriptomic signature corresponding to an increase in permissive epigenetic signatures around cancer-related genes in adulthood. This positively associates with BDE-99 mediated increase in A. muciniphila and its metabolites which are established epigenetic modifiers.

Chemicals. 2,2’,4,4’,5-pentabromodiphenyl-ether (BDE-99) (CAS No. 60348-60-9) and bisphenol A (BPA) were purchased from AccuStandard, Inc. (New Haven, CT). Fox River PCB mixture was prepared from technical PCB mixtures as we described previously (Lim et al. 2020). Corn oil was purchased from Sigma-Aldrich (St. Louis, MO). All other chemicals, unless otherwise noted, were purchased from Sigma-Aldrich St. Louis, MO). 

Animals

All mice were housed according to the Association for Assessment and Accreditation of Laboratory Animal Care International guidelines (https://aaalac.org/resources/theguide.cfm), and studies were approved by the Institutional Animal Care and Use Committee at the University of Washington. Eight-week-old specific pathogen-free C57BL/6J mice were purchased from the Jackson Laboratory (Bar Harbor, ME) and were acclimated to the animal facility at the University of Washington for at least two breeding generations prior to the experiment. Mice were housed in standard air-filtered cages using autoclaved bedding (autoclaved Enrich-N’Pure) (Andersons, Maumee, OH), and were exposed ad libitum to non-acidified autoclaved water, as well as LabDiet #5021 for breeding pairs, or LabDiet #5010 for post-weaning pups (LabDiet, St. Louis, MO). All chemical solutions were sterilized using the Sterflip Vacuum-Driven Filtration System with a 0.22 μm Millipore Express Plus Membrane (EMD Millipore, Temecula, CA). As shown in the study design diagram (Fig. 1A), litters were randomly assigned to chemical or vehicle exposures between postnatal day (PND) 2 to 4, pups were supralingually exposed to corn oil (vehicle, 5 ml/kg), BPA (250 μg/kg), BDE-99 (57 mg/kg), or the PCB Fox River Mixture (30 mg/kg), once daily for three consecutive days (n= 5/sex/exposure). Mice were euthanized using CO2 and decapitation (PND5), or CO2 followed by exsanguination via cardiac puncture (PND60). Various bio-compartments were collected at PND5 (24 hours after the exposure) and 60 (56 days after the last exposure), including small and large intestine sections (duodenum, jejunum, ileum, and colon), feces, liver, and serum. Samples were stored in a -80˚C freezer until further analysis as described in Fig. 1A and detailed below: 
16S rDNA sequencing and data analysis
Microbial DNA was isolated from duodenum, jejunum, ileum, colon, and feces at PND5 and 60, as well as liver and serum at PND60, using an OMEGA E.Z.N.A. Stool DNA Kit (OMEGA Biotech Inc., Norcross, Georgia). To note, the microbial DNA in liver and serum serves as an indirect biomarker for gut permeability. The concentration of DNA was determined using Qubit (Thermo Fisher Scientific, Waltham, Massachusetts). The integrity and quantity of all DNA samples were confirmed using an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, California). Bacterial 16S rDNA hypervariable region 4 (V4) amplicon sequencing was performed using a HiSeq 2500 second-generation sequencer (250 bp paired-end) (n = 5 per group) using a similar method as we described before (Cheng et al. 2018; Gomez et al. 2021; Lim et al. 2020; Scoville et al. 2019).
Microbial 16S rDNA paired-end sequence reads, in quintuplicate, were merged, de-multiplexed, quality-checked, and chimera-filtered using QIIME 2 (https://qiime2.org/)(Bolyen et al. 2019). The Silva 99 version 132 reference was used (Quast et al. 2013). The operational taxonomic unit (OTU) table annotated with classified bacteria information was read into R for further analysis (R Core Team, 2019). Differential abundance testing was performed using the ANCOM plugin in QIIME 2 (Mandal et al. 2015). Differentially abundant bacteria (relative abundance of OTU) and differentially regulated hepatic genes were correlated (Pearson’s correlation), and significant correlation was defined as |r| > 0.8 and false discovery rate (FDR)-adjusted p-value < 0.1. Correlation results were plotted using ComplexHeatmap (Gu et al. 2016). Functional predictions of the microbial composition were conducted using PICRUSt2 (Douglas et al. 2020).

RNA Isolation

Total RNA was isolated from frozen livers using the RNA-Bee reagent (Tel-Test lnc., Friendswood, TX) according to the manufacturer’s protocol. RNA concentrations were quantified using a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA) at 260 nm. The integrity of total RNA samples was evaluated by formaldehyde-agarose gel electrophoresis with visualization of 18S and 28S rRNA bands under UV light and confirmed by an Agilent 2100 Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA). Samples with RNA integrity (RIN) values above 8.0 were used for RNA-Seq.

Whole transcriptome Sequencing of liver samples. In triplicates, the cDNA library was constructed using a ribosomal depletion method, and reads were sequenced (75 bp paired end) per the Illumina manufacturer’s protocol. FASTQ files were de-multiplexed and concatenated for each sample. Quality control the FASTQ files was performed using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Reads were then mapped to the mouse reference genome (National Center for Biotechnology Information [NCBI GRCm38/mm10]) using HISAT2 version 2.1 (Kim et al. 2015). The sequencing alignment/map (SAM) files were converted to binary alignment/map (BAM) format using SAMtools version 1.8 (Li et al., 2009) and were analyzed by Cufflinks version 2.2.1 to estimate the transcript abundance (Trapnell et al. 2010) using Gencode mouse version 19 (vM19) gene transfer format (GTF). Long non-coding RNA (lncRNA) GTF from Gencode (vM19) was used to estimate the lncRNA transcript abundance. The abundance was expressed as fragments per kilobase of transcript per million mapped reads (FPKM) which was then converted to transcripts per million (TPM). Differential expression analysis was performed using Cuffdiff (Trapnell et al. 2010). The differentially expressed genes were defined as false discovery rate Benjamini-Hotchberg adjusted p value (FDR-BH) <0.05 in the chemical-exposed groups compared with the vehicle-exposed control group. 
RNA-Seq data analysis. All analyses were done in R unless stated otherwise. For RNA-Seq, samples were categorized by age, sex, and exposure. Genes were considered expressed if the average TPM was >1 for at least one group. Protein coding genes (PCGs) were separately categorized based on the Ensembl BioMart gene category. Differentially expressed PCGs were also categorized and matched with genes in categories of interest (xenobiotic biotransformation, epigenetic modifiers, nuclear receptors, and oxidative stress). These specific gene categories were based on literature searches and extracting gene sets from the Gene Ontology consortium and the Kyoto Encyclopedia of Genes and Genomes (KEGG). Venn diagrams were made for differentially regulated PCGs and lncRNAs for each sex and age group (PND5 and 60) using the R package VennDiagram (Chen and Boutros 2011). Hierarchical clustering was generated clustering using R package ComplexHeatmap (Gu et al. 2016). Up- and down-regulations were defined as the absolute value of the fold change greater than 1.5. Lists of up- and down-regulated genes were used as input for gene ontology enrichment using the R package topGO (Alexa and Rahnenfuhrer, 2020) for all groups. The list of genes in the unfiltered expression table was used as the background for gene ontology. Upstream regulators for all differentially regulated genes were determined using IPA. Gene expression values were plotted using bar plots for specific genes of interest using the arithmetic mean TPM in SigmaPlot. 

Chromatin immunoprecipitation coupled with second-generation sequencing (ChIP-Seq) and data analysis. 
Frozen pooled livers from PND60 males and females neonatally exposed to BDE-99 were subjected to ChIP with validated antibodies for the following three histone marks associated with active gene transcription: histone 3 lysine 4 monomethylation (H3K4me1) (Active Motif 39297) and trimethylation (H3K4me3) (Active Motif 39159), and H3K27 acetylation (H3K27ac) (Active Motif 39133). 75-nt single-end (SE75) sequence reads were generated by Illumina sequencing (using NextSeq 500) and mapped to the mm10 mouse reference genome using the BWA algorithm (Li and Durbin 2009). Reads that passed Illumina’s purity filter, aligned with no more than 2 mismatches, and mapped uniquely to the genome were used for the subsequent analysis. Peak calling, standard normalization, spike-in adjustment, merged region analysis, and genomic annotations were performed by Active Motif (Carlsbad, CA). Genes that had positive epigenetic mark enrichment within +/- 10 kb of the gene loci were considered target genes of that particular epigenetic mark. The union of overlapping intervals (merged region) was used to make Venn diagrams to illustrate common and unique epigenetic marks in control and BDE-99 exposed groups. A 30% change in the average peak value (i.e., BDE-99/ vehicle > 1.3 or < 0.7) was considered as a significant alteration in the epigenetic mark associated with the gene. ChIP-Seq and RNA-Seq datasets were overlaid to examine the association between the alterations in epigenetic marks and changes in target gene expression, and data are visualized as pie charts (Huang et al. 2013) and analyzed using topGO (Alexa A, Rahnenfuhrer J, 2020) and STRING (https://string-db.org/). 

Please refer to the article and Readme linked for full references of the methods.

RNA-seq, ChIP-seq, and 16s rRNA-seq fastq files are in respective zip compressed files. 

The metadata files, in each compressed files, for the RNA-seq, ChIP-seq, and 16s rRNA-seq are in the following respective spreadsheets: RNA-seq_metadata, ChIP-seq_metadta, 16s_metadata.

The following describes the raw data for each experiment:

The first six numbers indicate individual samples for the paired-end RNA sequencing fastq files, i.e. 205643-205690. This is followed by the strand, i.e. 1 for forward and 2 for reverse. 

Fastq files that start with a two-digit number followed by 08XX_0150UWash_ indicate fastq files for ChIP-seq. JL1, JL2, JL3, and JL4 indicate different samples. "male" in the fastq file title represents males and "fem" represents female samples. "cont" indicates vehicle control, and "treat" indicates BDE-99 exposed. Histone marks, i.e. H3K4me1, H3K4me3, H3K27ac, are indicated in the file names.

The paired-end fastq files for 16s rRNA-seq are named as follows:
Chemical exposure + biocompartment + age + sex + seven digit numbers + lane information + strand. 
For example, PCB.60.M.0016393_6_L001_R1_001.fastq represents Day 60 old males exposed to PCB and the fastq file is for the positive strand, and 

M60.jej.CO.0016160_222_L001_R1_001.fastq indicates Day 60 old males exposed to corn oil (vehicle) in the jejunum section. 

Note, duo = duodenum, jej = jejunum, ile = ileum, col = colon, fec = feces, liv = liver, ser = serum. Blank biocompartments in the file names are for feces.