Abstract

   Bisphenol A (BPA), a high production volume chemical widely used in
   consumer products, is an endocrine active compound associated with
   complex epigenetic responses in animal models and humans. Developmental
   BPA exposure in mice previously revealed widespread changes in the
   mouse liver methylome. Here, we undertake the first epigenome-wide
   analysis of the effect of BPA concentration on human fetal liver DNA
   methylation. Enzymatic enrichment of genomic DNA for high CG density
   and methylation followed by next-generation sequencing yielded data for
   positional methylation across the genome. Comparing three groups of
   BPA-exposed subjects (n = 18; 6 per group), high (35.44–96.76 ng/g),
   low (3.50 to 5.79 ng/g), and non-detect (<0.83 ng/g), revealed regions
   of altered methylation. Similar numbers of regions of altered
   methylations were detected in pairwise comparisons; however, their
   genomic locations were distinct between the non-detect and low or high
   BPA groups. In general, BPA levels were positively associated with
   methylation in CpG islands and negatively associated with methylation
   in CpG shores, shelves, and repetitive regions. DNA methylation at the
   SNORD imprinted cluster (15q11q13) illustrated both linear and
   non-monotonic associations with BPA levels. Integrated methylation and
   RNA-sequencing gene expression analysis revealed differential
   regulation of transcription at low BPA levels, as well as expression
   changes in RNA for ligand-binding proteins as BPA levels increase. BPA
   levels in human fetal liver tissue are associated with complex linear
   and non-monotonic as well as sequence-dependent alterations in DNA
   methylation. Future longitudinal studies are needed to link these
   changes with altered health risks.

   Keywords: bisphenol A, DNA methylation, environmental epigenetics

Introduction

   Environmental exposures during fetal growth can influence later-in-life
   health risks, including metabolic and phenotypic outcomes. The
   developmental origins of health and disease hypothesis posits that
   chemical and/or nutritional factors during early life result in lasting
   effects on disease risk, even in the absence of chronic exposure
   [[37]1, [38]2]. Growing evidence supports epigenetic inheritance of
   chromatin marks, such as DNA methylation, as a mechanistic link between
   fetal exposure and later susceptibility to disease [[39]3, [40]4].
   Adult phenotypic variation deriving from uneven re-establishment of DNA
   methylation during blastocyst formation is seen across the animal
   kingdom and can be correlated with environment in the form of stress,
   chemical exposures, nutrition, maternal behavior, and stochastic
   effects [[41]5]. Next-generation sequencing enables us to move beyond
   candidate-gene-based approaches by expanding coverage and sensitivity
   to detect previously unknown labile regions of the genome that are
   responsive to early life toxicant exposure. Here, we evaluate bisphenol
   A (BPA) levels and the DNA methylome in human fetal liver, as an
   extension of our previously published controlled BPA exposure in mice
   [[42]6], to uncover BPA’s association with the developing human
   epigenome ([43]Fig. 1).

Figure 1.

   Figure 1
   [44]Open in a new tab

   schematic overview of experimental design from BPA quantification in
   human fetal liver tissue through next-generation sequencing and
   analysis to identify RAMs. ND, non-detect

   BPA entered commercial production in the 1950s, and by 2008, it had
   reached global production of 11.5 billion pounds [[45]7], making it one
   of the highest volume production chemicals in the world. It is used in
   the manufacture of clear plastics, bottles, can liners, eyeglass
   lenses, cell phones, thermal receipt paper, and adhesives, among many
   other products. The safety of BPA has been the subject of long-standing
   scientific and public debate [[46]8]. A large proportion of the
   population has detectible levels of BPA in urine [[47]9]. Consumption
   of canned soup results in >1000% increase in urine BPA [[48]10] and
   dermal exposure to receipt paper can increase urine concentration from
   a baseline of 1.8 to 5.8 μg/l [[49]11]. BPA levels in urine have been
   correlated with obesity in Caucasian children [[50]12]. Pregnant women
   in southeast Michigan were determined to have between 0.5 μg/l and
   22.3 μg/l BPA in circulating blood [[51]13], and our recent study of
   human fetal livers indicated that most BPA is in the unconjugated form
   not readily eliminated from the body [[52]14]. These findings
   strengthen the case that early life development is a particularly
   crucial window to evaluate BPA’s effects on the epigenome.

   Our studies and others have found that developmental exposure to BPA is
   associated with epigenetic changes in specific tissues at specific
   genes [[53]15–20] and results in genome-wide changes in liver in
   rodents [[54]6, [55]21]. In a cross-sectional study of Egyptian girls,
   we identified BPA-associated DNA methylation alterations in saliva.
   Until now, there has been no genome-wide studies of altered epigenetic
   changes associated with fetal BPA exposure in human tissue. By
   combining enzymatic methods to enrich genomic DNA for high methylation
   and high GC content prior to next-generation sequencing, we identified
   regions of altered methylation (RAMs) between groups of fetal human
   livers, stratified according to quantified BPA levels. We validated
   selected top-hit genes to confirm DNA methylation differences and used
   RNA sequencing to show differential regulation of transcription at low
   levels of BPA, as well as expression changes in ligand-binding proteins
   as BPA levels increase. Thus, BPA levels in human fetal liver tissue
   were associated with complex and sequence-dependent alterations in DNA
   methylation, as well as by differentially regulating transcription.

Results

Analysis Pipeline and Quality Control for Identifying Differential
Methylation

   We used the MethylPlex-next-generation sequencing (M-NGS) platform to
   evaluate genome-wide DNA methylation patterns associated with various
   levels of BPA quantified in human fetal liver samples. This methodology
   requires minimal DNA input (∼50 ng) and enriches methylated DNA using a
   cocktail of methylation-dependent restriction enzymes prior to
   next-generation sequencing ([56]Fig. 1). We confirmed that MethylPlex
   library reads were enriched in genomic regions containing higher
   numbers of genes and CpG islands (CGIs). To estimate the false
   discovery rate (FDR) of the data analysis pipeline (see Methods for
   details), we employed a sex-based analysis comparing methylation
   profiles between female and male subjects. Under the conservative
   assumption that all autosomal regions that pass our filters for
   significance are false positives, the FDR was estimated to be 11.4%
   ([57]Supplementary Fig. S1A); however, the actual FDR may be lower, to
   the extent that true autosomal differences in methylation exist between
   sexes. The difference in mapped reads on chromosomes X and Y was
   clearly distinguishable between male and female subjects with minimal
   background noise observed on chromosome Y from female subjects,
   confirming the noted sex of each sample ([58]Supplementary Fig. S1B).

BPA-Associated RAMs

   We identified BPA-associated RAMs using a moving window approach, the
   edgeR Bioconductor package, and post-processing filtering steps, across
   three BPA categories (non-detect vs. low, non-detect vs. high, and low
   vs. high), and conducted a refined downstream analysis. As part of data
   exploration, overall across autosomal chromosomes, we observed a
   greater number of hypomethylated RAMs with increasing BPA levels, when
   non-detect subjects were compared with either low or high BPA subjects
   ([59]Fig. 2A[60] and B). In contrast, when comparing low BPA to high
   BPA subjects, approximately similar numbers of genomic regions were
   identified as hyper- and hypo-methylated with increasing BPA levels
   ([61]Fig. 2C). Chromosome level window counts were not tested for
   significance.

Figure 2.

   Figure 2
   [62]Open in a new tab

   the chromosomal distribution of hyper-methylated (in red) and
   hypomethylated (in blue) sites as counted by 100 bp sliding windows per
   1 million bp, normalized by chromosome length. (A) Number of windows in
   the non-detect vs. low groups are normalized to chromosome length and
   overall reflect the GC density of each chromosome. The majority of
   windows were hypo-methylated except for chromosome X. (B) Number of
   windows in the non-detect vs. high BPA groups, similar to (A). (C)
   Number of windows in low vs. high BPA groups shows a balanced number of
   hyper- to hypo-methylated regions. Overall data were not tested for
   significance

   Comparing low BPA subjects with high BPA subjects resulted in the
   largest number of 100 bp windows with RAMs (11 194). When subjects with
   non-detectable levels of BPA were compared with those with low and high
   levels of BPA, similar numbers of 100 bp windows with RAMs were
   identified (6286 and 7337, respectively), yet these regions were highly
   mutually exclusive with only 634 windows overlapping ([63]Fig. 3A). A
   majority of RAMs (19 522 out of a total 24 817 windows) were distinct
   from one another, suggesting locus-specific and non-monotonic effects
   of DNA methylation patterns. To identify genic regions associated with
   RAMs, each candidate 100 bp window was mapped to the nearest gene. In
   total, 10 005 RAMs were identified within 5 kb of a TSS ([64]Fig. 3B).
   A total of 296 RAMs were shared across all three comparisons. In
   contrast to the superset of all RAMs, the TSS-associated RAMs showed
   greater overlap.

Figure 3.

   [65]Figure 3
   [66]Open in a new tab

   overlap of 100 bp windows with RAMs in non-detect vs. low, non-detect
   vs. high, and low vs. high BPA groups (A) Genome-wide and (B) within
   5 kb of TSSs of genes

   We also examined the DNA methylation landscape of RAMs around CGIs, CGI
   shores (0–2 kb from a CGI), and CGI shelves (2–4 kb from a CGI) among
   the three BPA groups (non-detect, low, and high). We observed a trend
   of increased methylation in CGIs and decreased methylation in CGI
   shores and shelves with increasing BPA levels ([67]Fig. 4). These
   trends did not reach statistical significance.

Figure 4.

   [68]Figure 4
   [69]Open in a new tab

   characterization of genome-wide DNA methylation patterns in CGIs, CGI
   shores, and CGI shelves. Increasing BPA levels were positively
   associated with DNA methylation in CGIs and negatively associated with
   DNA methylation in CGI shores and shelves

   In addition, we examined the genomic distribution of RAMs relative to
   various genomic regions: promoters, or exons, introns, and intergenic
   regions, or repeat elements. The overview of the statistics of genomic
   regions overlapping with RAMs is available in [70]Supplementary Table
   S3. In all three comparison pairs, there were more hypermethylated than
   hypomethylated regions in the higher dose group in transcriptional
   start regions (TSRs), exons, and promoters ([71]Fig. 5A). A similar
   pattern is observed in genomic areas including gene locus, defined as
   the entirety of the gene transcript, as well as for introns alone
   ([72]Fig. 5B). However, the reverse pattern is observed in intergenic
   regions; there were more intergenic hypomethylated regions than
   hypermethylated in higher dose groups ([73]Fig. 5B). In the non-detect
   vs. low comparison, around 40.4 of hypo- and 43.2% of hypermethylated
   regions overlapped with repetitive regions. In both non-detect vs. high
   and low vs. high comparisons, we observed a striking increase in the
   percentage of RAMs in repetitive regions among hypomethylated regions
   in the higher dose group (73.5 and 72.2%, respectively), while the
   percentage of hypermethylated RAMs in repetitive regions remained
   similar, at 45.2 and 49.6%, respectively ([74]Fig. 5B).

Figure 5.

   Figure 5
   [75]Open in a new tab

   genomic distribution of RAMs. TSRs are defined genomic regions where
   experimentally verified transcriptional initiation has taken place.
   Exons are regions of coding sequence. Promoters are defined as −500 to
   +100 of the TSR. Gene locus encompasses the entirety of the gene,
   including introns. Intergenic regions are defined as regions outside
   genic and regulatory regions. Repeats are matched against the UCSC
   Genome Browser RepeatMasker track. (A) A greater percentage of RAMs
   overlap with exons than TSRs or promoters and have >2-fold higher RAMs
   overlapping in the relatively higher exposure group, for all
   comparisons. (B) Among gene locus, intergenic regions, introns, and
   repeats, only hypermethylated RAMs overlapping with repetitive regions
   are less in the relatively higher exposure groups

Validation of RAMs

   To investigate the correlation of BPA level with imprinted gene DNA
   methylation, we closely examined the snoRNA cluster surrounding
   SNORD116, which is a maternally imprinted locus implicated in
   Prader–Willi syndrome. In MethylPlex sequencing, SNORD116 was
   hypomethylated in the non-detect group and hypermethylated in both the
   low and high BPA groups ([76]Fig. 6). Validation via pyrosequencing of
   four CpG sites in the original samples showed linear increased
   methylation between non-detect and low sample groups (P < 0.1) but also
   detected a non-monotonic response with the average of the low BPA
   showing higher methylation than the high BPA group (P < 0.05).

Figure 6.

   [77]Figure 6
   [78]Open in a new tab

   differential methylation at the SNORD gene cluster on chr15. (A) Each
   row in the upper panel represents the average of the six samples in the
   group. ND, non-detected BPA; low, low BPA level; high, high BPA level.
   (B) Four CpG sites within the SNORD cluster were averaged for
   validation using pyrosequencing. Exact locations are listed in
   [79]Supplementary Material and [80]Table S2

   A top hit intergenic locus, located on chromosome
   7:152888851–152888950, was selected for quantitative DNA methylation
   validation via EpiTYPER analysis. MethylPlex sequencing resulted in a P
   value < 1.5 × 10^−^10 (FDR < 4 × 10^−^5) between the low BPA and high
   BPA groups with the high BPA group displaying less methylation. This
   pattern was validated by quantitative bisulfite sequencing verifying
   the non- monotonic DNA methylation ([81]Supplementary Fig. S2).

RNA-Seq for Differential Expression

   Using RNA-seq, we identified seven, one, and three differentially
   expressed genes in fetal liver for non-detect vs. low, non-detect vs.
   high, and low vs. high, respectively, using an FDR < 0.10 and 2-fold
   cutoffs. Although these numbers were small, several high ranked genes
   were noted to be involved in estrogen processes, e.g. AHR [P = 0.012
   (non-detect vs. low) and P = 0.0046 (non-detect vs. high)], PGR (P
   value = 0.020), GPER (P value = 0.025), and PRL (P = 0.046). Several
   genes had known BPA-related interactions, e.g. genes in the Wnt
   signaling pathway were significantly enriched (WNT2, WNT5a, WNT7A, and
   WNT11). Genes related to pregnancy or placental development were also
   found (CGA, CSH1, CSH2, INSL4, FLT1, ADAM12, PSG4, and PSG5). The
   remaining genes came from a variety of categories and include SERPINE1,
   INDO, IGFBP1, AREG, LOXL1, LTBP2, FGF9, FGF23, EGFL6, LRP2, JAG1,
   NOTCH3, THBS2, and APOD.

Pathway Enrichment Analysis

   We performed pathway enrichment analysis with hyper and hypo-methylated
   RAMs using ChIP-Enrich ([82]http://chip-enrich.med.umich.edu), which
   assigns peaks to genes based on a chosen method (we used nearest
   transcription start site, TSS) and tests peaks from ChIP-seq
   experiments for enrichment of biological pathways, GO terms, and other
   types of gene sets. Pathway analysis found a strong enrichment of genes
   involved in metabolism in the low versus high comparison, including
   regulation of nucleotide metabolic process (q value = 0.023), which was
   enriched among hypermethylated regions, and regulation of metabolic
   process (q value = 0.017), which was enriched among the hypomethylated
   regions. Genes predicted to be targeted by the aryl-hydrocarbon nuclear
   receptor translocator were also significantly enriched with
   hypermethylated regions (q value = 0.0098). Among the genes regulated
   by this transcription factor is the Wilms tumor gene (WT1), an
   imprinted locus. The most consistent enrichment observed across
   comparisons and both hypo- and hyper-methylation was sets of genes
   down-regulated in response to estrogen. These gene sets were derived
   from GEO experiments [83]GSE11324 (3, 6, and 12 h of estrogen
   treatment), [84]GSE11791, and [85]GSE10879, and all were performed
   using MCF7 cells ([86]Supplementary Table S4). In particular, breast
   carcinoma amplified sequence 4 and breast cancer anti-estrogen
   resistance 3, both involved in the growth and maintenance of breast
   cancers. In addition, an enrichment analysis based on transcriptomic
   data from RNA-seq was performed using a logistic regression-based
   pathway enrichment analysis tool (LRpath) available at
   [87]http://lrpath.ncibi.org. When enriched concepts among
   hypomethylated regions and over-expressed genes were compared, the
   metabolic and catabolic processes as well as cell cycle were
   significantly enriched in both analyses (q < 0.05).

Integrative Analysis of DNA Methylation and Gene Expression

   From the RNA-seq data, we identified 273 (non-detect vs. low), 430
   (non-detect vs. high), and 484 (low vs. high) candidate genes, whose
   expression levels are inversely correlated with methylation levels
   ([88]Supplementary File S1). When the top enriched biological functions
   were compared across the three exposure comparison pairs, there were
   overlapping categories of enrichment shared between the non-detect vs.
   low and non-detect vs. high groups: across all three comparisons,
   calcium ion binding genes were enriched (P < 5 × 10^−^4). Genes in this
   category included several protocadherins, which are expressed in
   nervous tissue and includes protocadherin19 which, when mutated,
   results in an epileptic phenotype in humans [[89]22]. The shared
   categories between the non-detect vs. high and the low vs. high
   comparisons included ligand-binding functions, specifically calcium and
   sulfur compound binding functions (P < 4 × 10^−^4). Several genes in
   the sulfur compound binding function have also been shown to interact
   with BPA according to the Comparative Toxicogenomics Database [[90]23],
   including lactotransferrin, chemokine ligand 10, follistatin, glypican
   4, r-spondin 1, and acyl-coa binding domain containing 5 [[91]24].
   Non-detect vs. low or high compared with the low vs. high groups shared
   only the calcium binding function.

Discussion

   Adaptations of next-generation sequencing to assess DNA methylation and
   chromatin states are still rapidly evolving, with a variety of
   currently available approaches. In this study, genomic DNA from
   BPA-characterized human fetal liver tissue was enriched for methylated
   regions and high CG density, giving a proxy for differential DNA
   methylation via regional read-counts. We previously used this approach
   to show that the mouse liver epigenome responds to BPA exposure at
   comparable levels found in humans [[92]6]. Our bioinformatics analysis
   pipeline has been used in mouse [[93]6] and human [[94]25] studies and
   is validated here by both RNA-seq and quantitative bisulfite sequencing
   via EpiTyper and pyrosequencing methylation analysis in selected loci.
   In contrast to the mouse, we find large shifts in RAMs at CGIs. We
   identified an asymmetric difference in hypo vs. hypermethylated loci in
   the two exposure groups compared with the non-detect group, while
   finding similar numbers when comparing low versus high BPA groups. This
   suggests that BPA exposure biases genes towards hypomethylation in
   general, as we have previously seen in saliva DNA from a cohort of
   prepubescent girls [[95]26] and in mouse tail DNA [[96]16]. However,
   the total number of regions with altered methylation was largest
   between low and high exposure groups, suggesting that there is a
   dose-dependent effect. Though our analysis found RAMs between groups to
   be largely distinct, enrichment analyses pointed to overlapping sets of
   affected pathways; future studies should be conducted with more
   granular exposure groups and larger sample sizes to confirm whether the
   number of genes shifting in methylation concomitant with BPA exposure
   is nonlinear with dose.

   Growing interest in DNA methylation centers on CGI shores and shelves;
   e.g. in cancers, CGIs have been observed to become hypermethylated,
   while the surrounding regions tend to become hypomethylated, with an
   overall loss of signal distinguishing CGIs from their shores and
   shelves [[97]27]. We add to this data by showing increasing CGI
   methylation with BPA exposure and decreasing methylation in shores and
   shelves with BPA. These general results cannot indicate whether a
   specific locus will have altered expression but would be consistent
   with reduced expression. Overall, we observed hypomethylation averaged
   across genic regions; such patterns are also seen in aging cells and
   tissues [[98]28]. Similarly, long-term epigenetic drift can be impacted
   by early life exposures and also acts in a locus specific manner
   [[99]29].

   Consistent with the CGI results, RAM location analysis showed higher
   methylation with higher BPA exposure in the areas most relevant to gene
   expression: the promoters and TSRs. By raw read count, the highest
   enrichment by fold-increase of RAMs between all intergroup comparisons
   (when compared with the genome) occurs in promoter regions; however,
   the global response of BPA is hypomethylation. The RAM analysis narrows
   down the location of where the bulk of this hypomethylation is
   occurring. We observed a decrease in repeat region reads in the high
   BPA group compared with the non-detect group. Hypomethylation in LINE1
   repeat regions with occupational BPA exposure has been documented
   previously [[100]30]. There was only a small decrease in comparing the
   low BPA to the non-detect group, but the low vs. high group mirrored
   the non-detect vs. high group. Given the growing concern over
   transposon methylation and the importance of silencing [[101]31], BPA's
   activity in these genomic regions must be carefully examined in depth.

   Imprinted genes are an important bellwether of environmentally induced
   epigenetic change and the SNORD cluster identified here reflects BPA's
   complex effects. Mouse studies have demonstrated the effect of BPA on
   imprinted gene methylation through the estrogen receptor signaling
   pathway [[102]32] and on imprinted gene expression in embryos and
   placenta [[103]33]. Even low levels of neonatal BPA exposure
   (2.4 μg/pup) in rats cause persistent hypomethylation in imprinted
   regions in adult male spermatozoa [[104]34]. Here, despite the SNORD
   cluster's relatively small size and contiguous maternal imprint, SNRPN
   and adjacent imprinted genes exhibited both linear and non-monotonic
   responses to BPA exposure. The difference between the MethylPlex and
   pyrosequencing may be explained by their different regions of
   assessment: MethylPlex assessed the average of the region, while
   pyrosequencing probed four specific CpG sites. Interestingly, recent
   examination of genome-wide DNA methylation patterns in paternal sperm
   of autism risk children identified differentially methylated genes
   involved in developmental processes, including many genes in the SNORD
   family [[105]35].

   Pathway analysis was utilized to help understand the likely biological
   impact of the RAMs, supporting the known estrogenic role of BPA by
   correlating RAMs with estrogen stimulus response genes. In all
   comparisons of altered methylation, estrogen responding gene categories
   were statistically significantly enriched, suggesting that this class
   of genes is both up- and down-regulated by BPA exposure. Genes targeted
   by aryl-hydrocarbon nuclear receptor translocator—AHR complex were
   hypermethylated in low vs. high BPA, which is unsurprising given that
   BPA is an aryl-hydrocarbon receptor ligand. The RNA-seq GO results
   taken alone are less clear; however, since methylation does not
   correspond perfectly to gene expression, integrating the two datasets
   provides greater insight. Another potential impact may be seen in the
   detection of pregnancy and placenta-related genes in the RNA-seq GO
   results, suggesting that differential timing of fetal development may
   account for the enrichment of these categories across exposure groups.
   Thus, by integrating the results of methylation data with the RNA-seq
   data and analyzing GO term enrichment, we found several categories of
   biological function overlapping between the pairwise comparisons of the
   non-detect vs. low and non-detect vs. high, especially those relating
   to gene transcription, suggesting that BPA exposure at both levels has
   wide-ranging effects on gene expression. In contrast, the genes
   enriched in the low vs. high group are involved primarily in ligand
   binding. The only exception to this was in the calcium binding
   category. Importantly, given the sample source, our low and high groups
   both represent physiologically relevant concentrations; thus, these
   results promote the need to examine low-dose BPA exposures in animal
   studies where the dose ranges may exhibit large gaps, potentially
   failing to detect discrete but relevant changes.

   Ultimately, these results begin to reveal the complexity of BPA's
   response in human fetal liver at both the DNA methylation and the
   transcription level. While the global response to BPA exposure is
   hypomethylation in general, especially in repeats, the CGI
   hypermethylation effects appear to be both dose and locus dependent.
   Thus, BPA levels in human fetal liver tissue are associated with
   complex linear and non-monotonic as well as sequence-dependent
   alterations in DNA methylation. Furthermore, although the RNA-seq
   results represent a functional outcome specific to developmental stage
   (e.g. gestational age) and tissue, the DNA methylation changes may
   persist across developmental stages and tissues, with only a small
   subset of the RAMs having a functional consequence for any one
   time-tissue combination. These limitations can be addressed by future
   pregnancy studies with larger sample sizes and/or longitudinal studies
   to link changes in DNA methylation or gene expression with altered
   health risks and to evaluate target tissue DNA methylation profiles
   with surrogate DNA, such as saliva or blood, to best inform studies of
   epigenetic epidemiology, when target tissue DNA is not always readily
   or ethically available.

Materials and Methods

Human Fetal Liver Tissue Samples

   Human fetal liver samples (n = 50), ranging from gestational days 70 to
   120, were procured from the NIH-funded University of Washington Birth
   Defects Research Laboratory Fetal Biobank (2R24 HD000836-47). As
   previously described [[106]14, [107]15, [108]36], these healthy tissue
   specimens were collected from voluntary pregnancy terminations after
   surgery and proper consent from donors and were flash frozen and stored
   in polycarbonate-free tubing at −80°C until processed for BPA analysis
   and RNA/DNA extraction ([109]Fig. 1, top box). No identifying clinical
   data were available on subjects except for gestational age and
   occasionally sex and race. Thus, samples met the criteria for IRB
   exemption for human subjects research (UM IRB Exemption: HUM00024929).
   As described previously, sex was determined from subjects with missing
   data (N = 10) using nested polymerase chain reaction (PCR) assays
   specific for the Y-chromosome SRY and the X-chromosome ATL1 genes
   [[110]14].

   Total (free plus conjugated) BPA concentrations were measured via
   high-performance liquid chromatography coupled with an API 2000
   electrospray triple-quadrupole mass spectrometer (ESI-MS/MS) by the
   Kannan Laboratory at the Wadsworth Center (New York State Department of
   Health) in 0.5 g of n = 50 human fetal liver tissue samples and ranged
   from below the limit of quantification of 0.1 up to 96.8 ng/g [[111]12,
   [112]14]. From the wide range of total BPA concentrations quantified in
   the fetal liver samples, we trichotomized 18 samples with gestational
   ages ranging from 80 to 115 into non-detect (total BPA concentrations
   ranging from non-detect to 0.83 ng/g; mean gestational age: 96.2 days),
   low- (3.5 to 5.79 ng/g; mean gestational age: 105.8 days), and high-BPA
   (35.44 to 96.76 ng/g; mean gestational age: 103.3 days) exposure groups
   (n = 6 samples per group). The n = 18 samples for this analysis were
   selected from the available set (n = 50) based on high DNA and RNA
   quality, as well as age and gender match per exposure group. Thus, no
   significant exposure group differences by gestational age were present.

M-NGS Library Generation

   A MethylPlex library synthesis and GC-enrichment kit (Rubicon Genomics
   Inc., Ann Arbor, MI) was obtained, and the experiments were carried out
   according to the manufacturer’s protocol ([113]Fig. 1, second box).
   Methylplex enzymatic enrichment in combination with next-generation
   sequencing (M-NGS) is capable of identifying RAMs. The method used here
   was previously described for both human and mouse in Kim et al.
   [[114]6, [115]25]. Briefly, a patented cocktail of
   methylation-sensitive restriction enzymes were used with 50 ng of DNA
   prior to ligation to universal PCR primer and amplification to create a
   Methylplex library. A second step enzymatic treatment depleted most
   non-GC-rich DNA sequences, after which, DNA was re-amplified in a
   second round of PCR. The adaptor sequences were removed and the product
   purified prior to incorporation into the Illumina Hi-Seq sequencing
   platform at the end-repair step of the sample preparation kit according
   to manufacturer's instructions (Illumina Inc., San Diego, CA). After
   ligation to Illumina adaptors, the product was run out on a 2% agarose
   gel with the DNA excised and extracted at the 400 bp position using
   Qiagen gel extraction kit (Qiagen Inc., Valencia, CA).

M-NGS Sequencing and Alignment

   The purified MethylPlex library was analyzed by Bioanalyzer (Agilent
   Technologies, San Diego, CA) prior to flow cell generation, where 10 nM
   of library was used to prepare flowcells with approximately 30 000
   clusters per lane. Sequencing was performed by the University of
   Michigan DNA Sequencing Core on an Illumina HiSeq, and the raw
   sequencing image data using 100 cycles of single ends were analyzed by
   the Illumina analysis pipeline ([116]Fig. 1, middle box). An average of
   104 million reads per sample (ranging from 78 to 125 million reads)
   were obtained, where approximately 83.9% of these were mapped uniquely
   to the human hg19 reference genome using the Burrows-Wheeler Aligner
   tool (BWA version 0.5.9-r16, parameters: aln -q 6 -t 6)
   ([117]Supplementary Material and [118]Table S1).

Differential Methylation Analysis

   We applied our tiered-based profiling pipeline developed for mouse BPA
   exposure studies to identify human locus-specific RAMs ([119]Fig. 1,
   bottom two boxes), which is described in detail in Kim et al. [[120]6].
   Briefly, we scanned the entire genome using a window size of 100 bp
   with a 50 bp moving shift size, which accounts for over 61 million
   windows for each sample and obtained the number of mapped reads per
   100 bp window per sample. Genomic regions were then filtered to those
   with at least three samples having a read count greater than 20. The
   resulting 1.34 million regions were then tested for differential
   methylation using the edgeR Bioconductor package, which we used to test
   for differences between each pair of exposure groups [[121]37]. The
   edgeR analysis was run using the glmFit function with tagwise
   dispersion estimation, which uses a negative binomial model with
   moderated dispersion estimates, and identified the regions with
   differential methylation in three different comparisons; the
   methylation levels from the non-detect group (n = 6) against the low
   exposure group (n = 6), non-detect group against high exposure group
   (n = 6), and low exposure group against high exposure group. To
   minimize the sample-specific methylation variation affecting the
   results, and because edgeR is sensitive to outliers, we further
   analyzed the data with filtered RAMs that (i) exhibited methylation
   change in at least half of the samples per exposure group and (ii)
   displayed differential methylation either in at least one out of two
   flanking windows (shift window size is 50 bp) or two 100-bp windows
   within a 500 bp stretch.

   We used the X and Y chromosomes to confirm the sex of each sample and
   to estimate the overall FDR resulting from using the filtering steps
   described above. Using the methylation reads mapped to chromosome X and
   Y, the sex of each human fetal sample (8 males vs. 10 females) was
   re-confirmed, and the underlying methylation difference among male and
   female subjects was examined and easily visually distinguished. The
   above edgeR analysis and filtering steps were used to test for male
   versus female methylation differences. Eight random sets of six males
   and six females were generated and used for the analysis. Using the
   conservative assumption that all autosomal methylation differences
   detected were false positives, we calculated the average percent of
   significant sites that were autosomal as 11.4% ([122]Supplementary Fig.
   S1). Thus, we estimate that our FDR is no greater than 11.4% and less
   to the extent that true autosomal sex-specific differences exist. Three
   separate BED files containing 6286 unique genomic coordinates from
   non-detect vs. low, 7337 unique genomic coordinates from non-detect vs.
   high, and 11 194 unique genomic coordinates from low vs. high analysis
   were uploaded to Genomatix genome analyzer software (v2. Release 4.7)
   ([123]https://www.genomatix.de) and mapped to the human genome (NCBI
   build 37). The classification of regions identifying the overlap with
   exons, introns, promoters, and intergenic regions, as well as
   transcription start regions (TSRs) was performed using the RegionMiner
   workflow. Genomatix defines repeat regions using in-house libraries.

Gene Set Enrichment Testing

   Genome-wide region enrichment of GO terms was performed using the
   chipenrich Bioconductor package ([124]http://chip-enrich.med.umich.edu)
   [[125]38] with all genomic-regions that passed the filter described
   above. The nearest TSS locus definition and human reference genome
   assembly (hg19) were used. The results were visualized using the Reduce
   and Visualize Gene Ontology (REViGO) web application
   ([126]http://revigo.irb.hr), which removed redundant GO terms and
   linked highly similar GO terms with the similarity cutoff value of 0.7
   using the Homo sapiens database [[127]39].

Quantitative Methylation Validation

   Top candidate regions were selected based on various factors, including
   P values, the number of samples with RAMs, the number of reads, and the
   methylation status of adjacent regions. Of the two candidate regions
   selected for validation, the first was in the region located 1.5 kb
   upstream of SNORD116-1 in the SNORD116/SNURF-SNRPN cluster
   (chr15:25050592-25050745), and the other was intergenic on chromosome 7
   (chr7:152888651-152889150, hg38) proximal to an AluSg4 transposon
   exclusive to Hominidae between the genes ACTR3B and DPP6. Genomic DNA
   (500 ng) from the fetal liver tissue used for M-NGS was bisulfite
   treated using the EpiTect bisulfite kit (Qiagen Inc., Valencia, CA) to
   allow for the conversion of unmethylated cytosines to uracil (read as
   thymine during PCR amplification), whereas the methylated cytosines
   remain unconverted [[128]40]. Bisulfite-converted DNA (2 µl) was then
   amplified using Bio-Rad (Model #C1000) thermal cyclers (see
   [129]Supplementary Table S2 for primer information and PCR conditions)
   using primers targeting the SNORD gene cluster and the intergenic
   region on chromosome 7 (chr7:152888651-152889150). PCR products were
   run on a 1.5% agarose gel to ensure PCR quality, correct product
   length, and lack of contamination. Following manufacturers’ protocols,
   amplified SNORD products were analyzed for quantitative DNA methylation
   levels via the PyroMark Q96 MD Pyrosequencing system (Qiagen, Valencia,
   CA), and amplified chromosome 7 products were analyzed via the EpiTYPER
   platform (Sequenom, San Diego, CA). For each primer set, the
   methylation percentage across individual CpG sites was visualized along
   with the amplicon average for each sample. The differences in mean
   amplicon methylation levels in each paired group were tested using a
   two-tailed t-test.

RNA Sequencing and Analysis

   RNA-sequencing was performed on 12 human fetal liver samples (N = 4 per
   group), a subset of the 18 samples profiled with M-NGS, for
   differential expression analysis using 150 bp single end reads on the
   Illumina HiSeq analyzed by the Illumina analysis pipeline. Over 438
   million reads across 12 samples were obtained, where the average number
   of reads was 36.5 million per sample. The FastQC tool was used to
   perform quality control checks on raw data. Because of poor quality
   (<28) in the second half of the reads, we trimmed reads to 70 bases in
   length and aligned the reads to the hg19 human reference genome with
   TopHat, using parameters that accept only the best alignment for reads
   that align in more than one location and increased searching time to
   improve sensitivity. After using SamTools to remove duplicate reads to
   eliminate the effects of PCR duplicates, read counts were generated by
   CuffDiff software. Differential expression between the three pairwise
   comparisons (non-detect vs. low, non-detect vs. high, and low vs. high)
   was tested with edgeR Bioconductor package [[130]37], following the
   developers’ protocol (calcNormFactors, estimateCommonDisp,
   estimateTagwiseDisp, exactTest, and topTags), using q values (Benjamini
   and Hochberg FDRs) to adjust for multiple testing. Gene set enrichment
   testing for the RNA-seq data was performed using LRpath
   ([131]http://lrpath.ncibi.org) [[132]41] against GO and KEGG. Gene set
   enrichment for the combined Methylplex and RNA-seq datasets were
   analyzed with GOrilla ([133]http://cbl-gorilla.cs.technion.ac.il)
   [[134]42] when compared with the complete complement of human protein
   coding genes ([135]http://genenames.org). GO testing was performed
   against molecular function at a P value threshold of >10^−^3 and
   resulting categories are listed in [136]Supplementary File S1.

Supplementary Material

   Supplementary Data

   Supplementary Data
   [137]Click here for additional data file.^ (1.3MB, zip)

Acknowledgements