DOI: 10.5061/dryad.gf1vhhmpz
title: Neonatal exposure to BPA, BDE-99, and PCB produces persistent changes in hepatic transcriptome associated with gut dysbiosis in adult mouse livers

description: 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.
description: 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).  <strong>Animals</strong>
 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: <br> 16S rDNA sequencing and data analysis<br>
 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).<br> 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| &gt; 0.8 and false discovery rate (FDR)-adjusted
 p-value &lt; 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).
 <strong>RNA Isolation</strong>
 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) &lt;0.05 in the chemical-exposed groups compared with the
 vehicle-exposed control group. <br> 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 &gt;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. <br> 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 &gt; 1.3 or &lt; 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/). 
description: 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:<br> Chemical exposure +
 biocompartment + age + sex + seven digit numbers + lane information +
 strand. <br> 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.