DNA methylome changes by estradiol benzoate and bisphenol A links early-life environmental exposures to prostate cancer risk.

Developmental exposure to endocrine-disrupting chemicals (EDCs), 17β-estradiol-3-benzoate (EB) and bisphenol A (BPA), increases susceptibility to prostate cancer (PCa) in rodent models. Here, we used the methylated-CpG island recovery assay (MIRA)-assisted genomic tiling and CpG island arrays to identify treatment-associated methylome changes in the postnatal day (PND)90 dorsal prostate tissues of Sprague-Dawley rats neonatally (PND1, 3, and 5) treated with 25 µg/pup or 2,500 µg EB/kg body weight (BW) or 0.1 µg BPA/pup or 10 µg BPA/kg BW. We identified 111 EB-associated and 86 BPA-associated genes, with 20 in common, that have significant differentially methylated regions. Pathway analysis revealed cancer as the top common disease pathway. Bisulfite sequencing validated the differential methylation patterns observed by array analysis in 15 identified candidate genes. The methylation status of 7 (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, Tpd52, Creb3l4) of these 15 genes exhibited an inverse correlation with gene expression in tissue samples. Cell-based assays, using 5-aza-cytidine-treated normal (NbE-1) and cancerous (AIT) rat prostate cells, added evidence of DNA methylation-mediated gene expression of 6 genes (exception: Paqr4). Functional connectivity of these genes was linked to embryonic stem cell pluripotency. Furthermore, clustering analyses using the dataset from The Cancer Genome Atlas revealed that expression of this set of 7 genes was associated with recurrence-free survival of PCa patients. In conclusion, our study reveals that gene-specific promoter methylation changes, resulting from early-life EDC exposure in the rat, may serve as predictive epigenetic biomarkers of PCa recurrence, and raises the possibility that such exposure may impact human disease.

Abbreviations

African-American

Acrosin binding protein 21

Rat prostate cancer cell line

Protein kinase B

5-aza-cytidine

Bisphenol A

BTB domain containing 3

Body weight

Coiled-coil domain containing protein 67

Chondoadherin

Carbohydrate sulfotransferase 14

Cyclic AMP responsive element binding protein 3-like 4

Control

Differentially methylated region

Deoxyribonuclease II β

DNA methyltransferase 3a

DNA methyltransferase 3b

17β-estradiol

17β-estradiol-3-benzoate

Endocrine-disrupting chemical

Estrogen receptor

Extracellular-signal-regulated kinase 1/2

Estrogen related receptor gamma

Estrogen receptor α

Estrogen receptor β

G protein-coupled estrogen receptor 1

High mobility group nucleosome binding domain 5

Hyaluronan binding protein 2

Hippocalcin-like 1

Ingenuity Pathway Analysis

Keratin 83

Methyl-CpG binding domain protein 2

Methyl-CpG binding domain protein 4

Methylated-CpG island recovery assay

Rat prostate epithelial cell line

Nuclear factor I/A

Oxysterol binding protein-like 6

Progestin and adipoQ receptor family member 4

Prostate cancer

Phosphodiesterase type IV variant

Phosphohistidine phosphatase 1

Paired-like homeodomain 3

Prostatic intraepithelial neoplasia

Postnatal day

Protein kinase cAMP-dependent type 1 regulatory subunit α

Recombination signal binding protein for immunoglobulin kappa J region-like

Ring finger protein 186

Ribosomal protein L19

Sprague Dawley

Standard error of mean

Sex determining region Y box 2

Tumor-associated calcium signal transducer 2

T-box 4

The Cancer Genome Atlas

Transforming growth factor β

Transmembrane protein 27

Triple negative breast cancer

Tumor protein D52

Thymidine phosphorylase

Wingless-type MMTV integration site family, member 10B

Introduction

Exposure to environmental agents is a risk factor for multiple diseases. Evidence is mounting in support of the causal link between exposure to xenoestrogens and human diseases, including cancer. Bisphenol A (BPA) is a ubiquitous environmental xenoestrogen widely used during the production of polycarbonate plastics, epoxy resins, carbonless receipt paper, and hundreds of other manufactured products. The fact that urinary BPA is detectable in >95% of the US population indicates that BPA readily leaches into the environment, contaminates our food and water, and enters our bodies. Chronic exposure to BPA in rodents and in humans is associated with prostate cancer (PCa) and its pre-lesions. Moreover, BPA is detectable in umbilical cord blood and sera of newborns, which elevates the concern of early-life BPA exposure and development of later-life disease. As an endocrine disrupting chemical (EDC), BPA binds to estrogen receptors (ERs) ESR1, ESR2, and GPER1, which are expressed in the rodent and human prostate epithelial cells. While affinity for nuclear ERs is low compared to estradiol-17β, BPA has equivalent activational capacities for membrane ERs and rapid actions of low-dose BPA have been documented in prostate cells. Finally, BPA can also signal through non-ER pathways.

Like BPA, 17β-estradiol-3-benzoate (EB) is a ubiquitous environmental agent. It is an estradiol analog commonly used in livestock for inducing weight gain and synchronizing estrus cycles in heifers and cattle. Since EB binds ERs with high affinity, unintended exposures may trigger estrogenic responses that compromise health. In rodents, neonatal EB exposure (30–125 ng/day in mouse and 125–500 ng/day in rats) has been shown to induce the production of morphologically abnormal sperm during adulthood. Similarly, rats neonatally exposed to 25 µg EB developed prostates with a greater risk of malignant changes with aging, including severe prostatic intraepithelial neoplasia (PIN), a precancerous condition. Since production of morphologically abnormal sperm and induction of increased susceptibility to hormonal carcinogenesis are similarly observed following neonatal BPA exposure, these findings strongly implicate early-life as a susceptible window for environmental xenoestrogens to predispose later-life disease risk and EB can serve as a positive estrogen control in studies of weaker estrogenic agents such as BPA. While emerging evidence indicates that epigenetic modifications are involved, the molecular underpinnings of this process remain to be fully clarified.

To initially interrogate whether neonatal exposure to EB and BPA reprogrammed the prostate epigenome, we exposed rats to an environmentally relevant dose of BPA [10 μg/kg body weight (BW)] or EB (2,500 μg/kg BW) on postnatal days 1, 3, and 5. We observed that this transient developmental exposure increased the dorsal and lateral prostate lobe susceptibility to adult-onset E2-induced carcinogenesis, markedly augmenting the incidence and severity of lesions as compared to oil-treated controls. Using methylation-sensitive restriction fingerprinting to identify altered DNA methylation marks, we found that neonatal EB and BPA exposures were associated with the dysregulation of phosphodiesterase type IV variant 4 (Pde4d4), hippocalcin-like 1 (Hpcal1; also known as visinin-like protein-3, Vilip-3), and high mobility group nucleosome binding domain 5 (Hmgn5; formerly known as nucleosome binding protein 1, Nsbp1) through aberrant promoter methylation detected on postnatal day (PND)10, 90, and 200. Together, these findings provided the first evidence that early–life environmental exposure to EDC reprogrammed the prostate epigenome and identified a developmental basis of PCa risk with aging. Moreover, the expression of DNA methylation transferase3a and 3b (Dnmt3a and Dnmt3b) and methyl-CpG binding domain protein2 and 4 (Mbd2 and Mbd4) was upregulated in PND10 and 90 prostate tissues, providing a mechanistic basis for reprogramming of DNA methylation marks upon early-life EDC exposure. Thus, we propose that differential methylation of genes that persist in the PND90 prostate primes the tissue for heightened sensitivity to a secondary exposure to rising estradiol later in life, as occurs in aging males. This is particularly relevant since elevated estrogens have been associated with increased PCa risk in men and are sufficient to transform the human prostate epithelium.

While these findings had a major impact on the field of EDC research and developmental basis of carcinogenesis, a genome-wide search for additional DNA methylation targets of neonatal exposure to xenoestrogens is warranted for an unbiased discovery of other epigenetic marks that may underlie increased PCa risk with aging. To accomplish this goal, we herein used the methylated-CpG island recovery assay (MIRA)-assisted genomic tiling and CpG island array analysis and identified distinct and common EB-/BPA-associated genes in PND90 prostate tissues from rats with neonatal exposure to these xenoestrogens. Among 25 epigenetically regulated candidate genes, the promoter methylation status of 7 genes (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, Tpd52, Creb3l4) was inversely correlated to gene expression. These validated genes have functional connectivity associated with stem cell pluripotency. Of clinical relevance, expression of these genes was found to be associated with recurrence-free survival of 497 patients in The Cancer Genome Atlas (TCGA) PCa cohort, suggesting that they may have utilities for predicting disease progression and patient stratification based on disease aggressiveness.

Results

Methylation array revealed differential promoter methylation of genes associated with neonatal EB or BPA exposure in the PND90 prostate

We performed genome-wide methylation analysis, using MIRA-assisted genomic tiling and CpG island array, in dorsal prostate tissues from PND90 rats neonatally exposed to EB or BPA (Fig. 1). Using a 750 bp sliding window approach and a selection criteria of P < 10−5 and mBar >0.6 or <−0.25 (Fig. 1), we identified a total of 177 differentially methylated regions (DMRs), which were randomly distributed among chromosomes (Supplemental Figure S1). Since this array was designed to target CpG sites only at the gene promoter region, the genes described hereafter refer to those identified with DMRs at the 5′-promoter region. Of the 177 genes with identified DMRs, 111 genes were EB-associated, 86 genes were BPA-associated, and 20 genes were common between the 2 groups (Fig. 1; Supplemental Table S1).

Figure 1.

Schematic diagram of the experimental and analytical procedures. Sprague-Dawley (SD) rats were neonatally [postnatal day (PND) 1, 3, 5] treated with 17β-estradiol-3-benzoate (EB) at 25 μg/pup or 2,500 µg/kg body weight (BW), bisphenol A (BPA) at 0.1 µg/pup or 10 µg/kg BW, or corn oil as a control (Ctrl). DNA extracted from PND90 dorsal prostate tissue was subjected to promoter methylation array analysis, and genes with differentially methylated regions and inverse expression correlation were identified and validated.

Biological significance of differentially methylated regions associated with neonatal EB and BPA exposure in PND90 prostate

To identify biological processes related to the identified exposure-associated DMRs, we performed functional network analysis using Ingenuity Pathway Analysis (IPA). As shown in Table 1, the top 3 networks related to the EB-associated genes were identified as: 1) “tissue morphology, embryonic development, organ development;” 2) “cell-to-cell signaling and interaction, cell-mediated immune response, cellular growth and proliferation;” and 3) “post-translational modification, cellular assembly and organization, cellular function and maintenance.” Upon merging the molecular interactions of these top 3 networks, the EB-associated genes were found to converge at protein kinase B (AKT) and extracellular-signal-regulated kinase (ERK)1/2. While the EB-associated genes had implications in the development and function of several physiological systems, including embryonic, tissue, organ, and reproductive (Table 1), the top related disease was “cancer,” which included PCa.

Table 1.

Top networks and bio-functions of the genes associated with neonatal 17β-estradiol-3-benzoate (EB) or bisphenol A (BPA) exposure in postnatal day (PND)90 prostate.

EB-associated genes	BPA-associated genes
Number of genes = 111	Number of genes = 86
Top networks	Top networks
1. Tissue morphology, Embryonic development, Organ development	1. Cell-to-cell signaling and interaction, Cell-mediated immune response, Cellular growth and proliferation
2. Cell-to-cell signaling and interaction, Cell-mediated immune response, Cellular growth and proliferation	2. Nucleic acid metabolism, Small molecule biochemistry, Molecular transport
3. Post-translational modification, Cellular assembly and organization, Cellular function and maintenance	3. Cellular assembly and organization, Cellular function and maintenance, Cellular compromise
4. Molecular transport, Nucleic acid metabolism, Small molecule biochemistry	4. Cell cycle, Cell-to-cell signaling and interaction, Cellular growth and proliferation
5. Lipid metabolism, small molecule biochemistry, Endocrine system development and function	5. Lipid metabolism, Molecular transport, Small molecule biochemistry
Top bio-functions	Top bio-functions
Diseases and disorders	Diseases and disorders
1. Cancer	1. Antimicrobial response
2. Reproductive system disease	2. Cancer
3. Metabolic disease	3. Connective tissue disorders
4. Auditory disease	4. Developmental disorder
5. Cardiovascular disease	5. Gastrointestinal disease
Physiological system development and function	Physiological system development and function
1. Embryonic development	1. Digestive system development and function
2. Organ development	2. Organ morphology
3. Organismal development	3. Cardiovascular system development and function
4. Reproductive system development and function	4. Organismal development
5. Tissue development	5. Tissue morphology

Data were analyzed using Qiagen's Ingenuity Pathway Analysis (IPA; Qiagen, www.qiagen.com/ingenuity).

Similarly, the top 3 networks related to the 86 BPA-associated genes were identified as: 1) “cell-to-cell signaling and interaction, cell-mediated immune response, cellular growth and proliferation;” 2) “nucleic acid metabolism, small molecule biochemistry, molecular transport;” and 3) “cellular assembly and organization, cellular function and maintenance, cellular compromise” (Table 1). The molecular interactions of these networks converged at transforming growth factor β (TGFβ), and the 86 genes had implications in digestive and cardiovascular system development and function, organ and tissue morphology, and organismal development (Table 1). Like the EB-associated genes, BPA-associated genes were also highly associated with “cancer,” including PCa.

Interestingly, the top 3 networks related to the 20 associated genes common of both EB- and BPA-exposure were: 1) “cancer, organismal injury and abnormalities, renal and urological disease;” 2) “cell-to-cell signaling and interaction, nervous system development and function, cardiovascular disease;” and 3) “cancer, connective tissue disorders, dermatological diseases and conditions” (Supplemental Table S2C). While these 3 networks converged at tumor protein 53 (TP53), overall, the 20 genes were significantly related to digestive, hepatic, and nervous system development and function, and embryonic and organismal development (Supplemental Table S2). They also had implications in “connective tissue disorders,” “developmental disorder,” “hereditary disease,” “inflammatory disease,” and “metabolic disease” (Supplemental Table S2).

Selection of top candidate genes for methylation validation by bisulfite sequencing

To tighten the selection of differentially methylated candidate genes, we initially sorted the 177 identified genes based on P < 10−10 and mBar (EB/BPA – Control) ≥1 .5 or ≤−0.5 (Selection criteria A, Fig. 1). Under this stringent selection criteria, 20 candidate genes were identified (Supplemental Table S3), of which the top 9 most differentially methylated genes, based on differences in mBar values (Fig. 2), across exposure groups (Btbd3, Chst14, Creb3l4, Paqr4, Phpt1, Rbpjl, Sox2, Tacstd2, and Tpd52) were selected for validation. As shown in Figure 2, primers for bisulfite sequencing validation for each gene were designed to amplify the region with the greatest methylation changes between the EB/BPA-treated groups and the control group. In most cases, the interrogated regions were closely aligned to the predicted CpG island based on MethPrimer program.

Figure 2.

Representative results from genome-wide methylation study. A) Predicted CpG islands (light blue shaded areas) in the promoter region of differentially methylated genes (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, Tpd52, and Creb3l4) identified in this study. TSS stands for transcriptional start site whereas ATG stands for translational start codon. Individual CpG sites are represented by red vertical lines. Dark blue horizontal line marks the region selected for BS-sequencing; B) Position of the BS-sequenced region (blue line) relative to the NimbleGen probes covered regions (dirty yellow lines) of each gene; C) Significant methylation level of each gene in vehicle- (control), EB- and BPA-treated groups, which were measured by NimbleGen array probes (mBar values). The height of the mBar represents the probe intensity; red and green bars represent positive and negative methylation value, respectively, relative to their respective input control.

Promoter methylation status was validated by performing bisulfite sequencing analysis on the same PND90 prostate tissue samples used in the methylation array (Table 2; Supplemental Table S3). In accordance with the methylation patterns observed in the array, the DMRs of 7 out of 9 genes were confirmed. The DMRs of Chst14, Paqr4, Phpt1, Rbpjl, and Sox2 were hypomethylated, while those of Creb3l4 and Tpd52 were hypermethylated in the EB/BPA-exposure groups when compared to controls (Table 2).

Table 2.

Validation of promoter methylation status and gene expression in postnatal day (PND)90 prostate.

	Selection criteria#		Validated in EB and/BPA
Gene	A	B	Methylationa	Expressionb	Gene name
Acrbp		√	Hyper		Acrosine binding protein
Btbd3	√				BTB domain containing protein 3
Ccdc67		√			Coiled-coil domain containing protein 67
Chad		√	Hyper		Chondroadherin
Chst14	√		Hypo	↑	Carbohydrate sulfotransferase 14
Creb3l4	√		Hyper	↓	cAMP responsive element binding protein 3-like 4
Dnase2b		√			Deoxyribonuclease II beta
Habp2		√			Hyaluronan binding protein 2/FSAP
Krt83		√			Keratin 83
Nfia		√			Nuclear factor I/A
Osbpl6		√	Hyper		Oxysterol binding protein-like 6
Paqr4	√		Hypo	↑	Progestin and adipoQ receptor family member 4
Phpt1	√		Hypo		Phosphoshistidine phosphatase 1
Pitx3		√	Hypo	↑	Paired-like homeodomain 3
Prkar1a		√			Protein kinase cAMP-dependent type 1 regulatory subunit alpha
Rbpjl	√		Hypo		Recombination signal binding protein for immunoglobulin kappa J region-like
Rnf186		√	Hyper		Ring finger protein 186
Sox2	√		Hypo	↑	Sex determining region Y box 2
Tacstd2	√				Tumor-associated calcium signal transducer 2
Tbx4		√			T-box 4
Tmem27		√	Hyper		Transmembrane protein 27
Tpd52	√		Hyper	↓	Tumor protein D52
Tymp		√	Hyper		Thymidine phosphorylase
Uox		√			Urate oxidase
Wnt10b		√	Hypo	↑	Wingless-type MMTV integration site family, member 10B

Selection criteria used for identifying genes with differential promoter DMR methylation in PND90 prostate of rats neonatally exposed to 17β-estradiol-3-benzoate (EB) or bisphenol A (BPA) vs. control group

Selection criteria A: P < 10−10 and mBar (EB/BPA - Control) ≥1.5 or ≤−0.5; Selection criteria B: P < 10−5 and mBar (EB/BPA - Control) ≥1.2 or ≤−0.5.

Promoter methylation validated by bisulfite sequencing: Hyper=hypermethylated; Hypo=hypomethylated

Gene expression validated by qPCR analysis: ↑ = up-regulated and ↓ = downregulated

Highlighted in light grey, genes with promoter methylation status validated

Highlighted in dark grey, genes with promoter methylation status and gene expression validated to be inversely correlated.

To include more genes in the analysis, we also sorted the 177 exposure-associated genes using wider selection parameters of P < 10−5 and mBar (EB/BPA – Control) ≥1.2 or ≤−0.5 (Selection criteria B, Fig. 1). Using selection criteria B, 91 candidate genes were identified (Supplemental Table S4). Bisulfite sequencing analysis was performed on the top 16 most differentially methylated genes across exposure groups (Acrbp, Ccdc67, Chad, Dnase2b, Habp2, Krt83, Nfia, Osbp16, Pitx3, Prkar1a, Rnf186, Tbx4, Tmem27, Tymp, Uox, and Wnt10b) (Table 2; Supplemental Table S4). Promoter methylation validation of the 16 genes, using PND90 dorsal prostate tissues, confirmed promoter hypomethylation of Pitx3 and Wnt10b, and promoter hypermethylation of Acrbp, Chad, Osbp16, Rnf186, Tmem27, and Tymp (Table 2).

Overall, using 2 selection criteria, there were 25 top differentially methylated candidate genes identified, of which the promoter methylation status of 15 were validated (Chst14, Creb3l4, Paqr4, Phpt1, Rbpjl, Sox2, and Tpd52 from selection criteria A; and Acrbp, Chad, Osbp16, Pitx3, Rnf186, Tmem27, Tymp, and Wnt10b from selection criteria B) (Table 2).

Validation of promoter methylation status and gene expression correlation

To elucidate the biological relevance of the DMRs within these gene promoters, we next performed gene expression analyses on the PND90 prostate tissues. As shown in Figure 3, expression of Pitx3, Wnt10b, Paqr4, Sox2, and Chst14 was significantly upregulated in the EB-exposure group when compared to the control. Of these genes, only Sox2 was significantly upregulated in the BPA-exposure vs. control group. In contrast, expression of Tpd52 was significantly downregulated in both exposure groups, and Creb3l4 showed a trend of downregulation in exposure groups, when compared to the control group (Fig. 3). Since gene expression of Acrbp, Chad, Osbp16, Rnf186, Tmem27, and Tymp in the prostate tissues of the exposure groups did not associate with their validated promoter methylation status, they were not studied further (Table 2; data not shown).

Figure 3.

Effect of neonatal exposure to EB or BPA on promoter methylation and gene expression in PND90 dorsal prostate. Promoter methylation status (left panel: Scatter Plot) and expression (right panel: Bar Graph) of candidate genes, in PND90 prostate tissues from SD rats treated with either corn oil (Ctrl; white), EB (green), or BPA (red), were analyzed using bisulfite sequencing and qPCR, respectively. Each circle in the scatter plot represents mean ±SEM of methylation percentage (averaged from 6 individual samples/animals) at a single CpG site in the gene promoter region. The % methylation of each CpG site in each sample was determined from bisulfite sequencing data derived from 8–12 clones. Gene expression data were expressed as mean ± SEM from 3 individual samples. Statistical significance was determined by one-way ANOVA and Tukey test when compared to Ctrl. *P < 0.05, **P < 0.01, and ***P < 0.001.

To investigate the collective biological relevance of the 7 candidate genes (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, Tpd52, and Creb3l4), we performed functional network analysis using IPA. We found that the candidate genes were associated with 2 networks: 1) “embryonic development, organismal development, gene expression” and 2) “cell cycle, embryonic development, gene expression” (Supplemental Table S5), which converged at Sox2 (Fig. 4). While these 7 candidate genes had implications in embryonic, organ, and organismal development, as well as auditory, vestibular, and nervous system development and function, “cancer” was determined to be the top-related disease associated with the differentially methylated genes with inverse expression correlation (Supplemental Table S5). Complementary reported bio-functions are summarized in Table 3.

Figure 4.

Signaling pathways associated with differentially methylated candidate genes were predicted by Ingenuity Pathway Analysis™. Genes with promoter hypermethylation are shown in red and genes with promoter hypomethylation are shown in green, with color intensity signifying the magnitude of differential methylation. Gray arrows indicate predicted association pathways; purple arrows indicate reported directional pathways.

Table 3.

Reported bio-function(s) and PCa associations of the seven identified candidate genes with DMRs inversely correlated with gene expression as a result of neonatal 17β-estradiol-3-benzoate/bisphenol A (EB/BPA) exposure.

	Promoter methylation		Gene expression
Gene	EBa	BPAb	EBa	BPAb	Reported bio-function(s)	Gene expression in PCa	Other pathological diseases	Reported stem/progenitor cell bio-function
Pitx3	*		*		Knockdown of Pitx3 promoted the loss of lens epithelial cells by inducing early activation of cell cycle inhibitors71	Not reported	Downregulated in Chinese patients with Parkinson's disease72	Expression of Pitx3 initiated differentiation of mid-brain derived neural stem cells73
Wnt10b	*		*		Upregulating Wnt10b expression in endometrial cells Ishikawa increased cell proliferation and reduced apoptosis74Overexpression of Wnt10b in dermal papilla cells increased cell proliferation75Overexpression of Wnt10b induced S-phase cell progression in mouse mammary epithelial and human BCa cells50Overexpression of WNT10B in TNBC increased cell proliferation by upregulating oncogene high mobility group A hook-2.76	Upregulated in prostate tumors when compared to normal prostate tissues57Upregulated in prostate tumors in a stage-dependent manner77	Upregulation is associated with poor survival of patients with osteosarcoma5Upregulated in endometrial cancer tissues when compared to hyperplastic or normal tissue samples74	Treatment of human embryonic stem cell colony culture induces the formation of prostate organoids78Wnt10b conditioned medium promoted proliferation, morphological changes, and chemotaxis of immortalized and primary mouse spermatogonial stem cells79
Paqr4	*		*		Not reported	Not reported	Pancreatic cancer80	Not reported
Sox2	*	*	*	*	Overexpression of Sox2 in DU145 increased cell migration and proliferation, reduced apoptosis,51,52 and promoted tumorigenesis48Conditional Sox2 deletion in pre-existing skin papilloma and skin squamous-cell carcinoma (SCC) reduced tumor growth and cancer cell propagation in serial transplantation81Overexpression of Sox2 in chick embryos inhibited neural differentiation82	Upregulated in PCa83Upregulated in PCa with higher Gleason score48,49	Downregulated in colon cancer84 and gastric cancer7Upregulated in squamous cell lung cancer85, SCC81 nasopharyngeal carcinoma86, and Frequently expressed in basal-like breast cancer tumors when compared to luminal and HER2+ tumors87	Overexpression of Sox2 protein in mouse embryonic stem (ES) cells trigger their differentiation into cells that express markers for various differentiated cell types88Overexpression of Sox2 suppresses neural differentiation through downregulating GFAP expression4.Sox2 expression maintains osteoblast self-renewal89Deletion of Sox2 in neural stem cells (NSC) reduced the size and formation of neurospheres90Inhibiting Sox2 signaling decreased progenitor marker expression and promoted neuronal differentiation82
Chst14	*		*		Loss of Chst14 in neuronal stem cells reduce cell proliferation91	Not reported	Downregulated in colorectal cancer84	Loss of Chst14 decreased neurogenesis and proliferation of neural stem cells91
Tpd52	*		*	*	Overexpression of TPD52 increased LNCaP cell proliferation and migration52Mice inoculated with TPD52 overexpressing fibroblasts have been shown to develop lung metastasis92	Upregulated in high-grade PIN53,93 and PCa52-54Higher expression in Gleason score 4 than 3 prostate tumors,93 in localized PCa with Gleason score 4+4=8 than 3+3=6 localized PCa94, and in hormone naïve and refractory metastatic PCa95	Upregulated in ovarian cancer96Upregulated in testicular germ cell tumors97	Higher expression in the earlier passages of differentiating human dental pulp stem cells98
Creb3l4	*				Overexpression of Creb3l4 in LNCaP cells upregulated genes related to PCa99.	Upregulated in PCa when compared to adjacent normal tissues55,56	Not reported	Loss of Creb3l4 increased apoptosis of meiotic/post-meiotic germ cells in the mouse seminiferous tubules100

Treatment with 2,500 µg EB/kg BW/day on PND1, 3, and 5.

Treatment with 10 µg BPA/kg BW/day on PND1, 3, and 5.

indicates a statistically significant change (P < 0.05) in promoter DMR methylation (validated by bisulfite sequencing) or gene expression (validated by qPCR) between treatment and control groups

Light grey, promoter hypomethylation or decreased gene expression; dark grey, promoter hypermethylation or increased gene expression

Effect of 5-aza-cytidine on the expression of differential methylation candidate genes in rat prostate cells

To determine whether the expression of the 7 identified candidate genes (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, Tpd5, and Creb3l4) were directly regulated by DNA methylation, we compared gene expression in a rat prostate epithelial cell line NbE-1 and the rat PCa cell line AIT, either with or without a 4-day treatment of the DNA methylation inhibitor, 5-aza-cytidine (5-aza) at doses of 0.5 µM and 1 µM. When compared to the untreated control, 5-aza treatment significantly upregulated Pitx3, Wnt10b, Chst14, and Creb3l4 in NbE-1 cells, and Pitx3, Wnt10b, Sox2, Chst14, Tpd52, and Creb3l4 in AIT cells (Fig. 5). Although the same treatment increased the expression of Sox2 and Tpd52 in NbE-1 cells, the change was not significant. Overall, the 5-aza treatment had no effect on the expression of Paqr4 in both NbE-1 and AIT cells.

Figure 5.

Effect of 5-aza-cytidine treatment on gene expression in NbE-1 and AIT cells. Gene expression was analyzed by qPCR in rat (A) normal prostate epithelial NbE-1 cells and (B) prostate cancer AIT cells treated with DMSO (Ctrl), or 0.5 µM or 1 µM 5-aza-2-deoxycytidine (5-aza), a DNA methylation inhibitor, for 8 d. Data (mean ± SEM ) is normalized to corresponding Rpl19 levels, and is expressed as fold change vs. Ctrl. *P < 0.05, and **P < 0.001 by one-way ANOVA and Tukey test when compared to Ctrl.

Association of the identified candidate genes with the recurrence of PCa

To gain clinical significance, RNAseq data from 497 PCa specimens and their associated clinical data (Biotab) were retrieved from the TCGA database. Extensive data transformation and biostatistical analyses revealed that the expression of these 7 genes in human PCa samples correlated with recurrence-free survival, i.e., the cohort of patients was segregated into those with longer (Group 1) and shorter (Group 2) duration to cancer re-appearance after initial treatment. For the initial 5-year survival after diagnosis, Group 2 patients had a shorter time to recurrence (poorer recurrence-free survival) when compared with Group 1 (Fig. 6A, Hazard ratio = 1.786, log rank test P = 0.0411), suggesting Group 2 had significantly more aggressive disease. Consistent with this finding, we found that Group 2 had significantly more cases with higher tumor Gleason score (≥7) when compared with Group 1 (Odds ratio = 0.3423; P < 0.0001) (Fig. 6B). It is well established that men with higher grade tumors have higher risk of dying from PCa.

Figure 6.

Expression of the 7 candidate genes was associated with shorter recurrence-free survival of PCa patients. TCGA data consisting of 497 PCa patients were dichotomized into Groups 1 and 2 by K-means clustering analysis based on the 7 candidate genes. (A) Group 1 patients have longer time to recurrence than Group 2 patients; the two groups differ in recurrence-free survival. (B) Proportion of PCa patients having tumors with high (≥7) and low (<7) Gleason score.

Discussion

In the present study, we conducted a comprehensive genome-wide search for DNA methylation targets of developmental exposure to EB and BPA, with a view to identifying candidate epigenetic biomarkers for predicting PCa risk. We reported the genome-wide promoter methylation changes in PND90 dorsal prostate tissues related to neonatal exposure to EB or BPA. Overall, we identified 111 EB-associated and 86 BPA-associated DMRs, with 20 in common between the 2 groups. Using two selection criteria, 25 candidate genes were selected, of which 15 were validated to have treatment-associated promoter methylation changes using bisulfite sequencing analysis, and 7 identified to exhibit an inverse correlation between promoter methylation status and gene expression. EB treatment upregulated the expression of Pitx3, Wnt10b, Paqr4, Sox2, and Chst14 through promoter hypomethylation, and downregulated Tpd52, and Creb3l4 through promoter hypermethylation. Similarly, BPA treatment upregulated Sox2 and Chst14 through promoter hypomethylation, and downregulated Tpd52 through promoter hypermethylation, but had no significant effects on the other genes. Of the 7 candidate genes identified in the rat, WNT10b, SOX2, TPD52, and CREB3L4 were upregulated in PCa (Table 3), and PAQR4, WNT10b, PITX3, and TPD52 were differentially methylated between PCa and adjacent normal prostate tissues. Clustering analysis, based on the expression of these 7 genes in PCa tissues of 497 PCa patients from the TCGA data sets, further revealed an association with shorter recurrence-free survival. Taken together, DNA methylation-mediated dysregulation of this unique signature of 7 genes established a causal link between early-life EDC exposure and PCa risk, and raises the possibility that such exposure may impact the human disease.

Neonatal development is a susceptible window for the DNA methylation reprogramming effects of EDCs. Since the prostate is estrogen responsive, it is highly vulnerable to insult by estrogenic agents during development. Using a neonatal EB/BPA exposure rat model, we previously reported that treatment-associated dysregulation of Pde4d4, Hmgn5, and Hpcal1 through promoter methylation changes, persisted in the prostate throughout life. While we observed an increase in the incidence and severity of PIN in PND200 rats upon E2-induced carcinogenesis, the current study aimed to interrogate the underlying mechanism for the early-life EDC exposure-associated increase in PCa risk, by performing genome-wide promoter methylation analysis on PND90 prostate tissue. The 111 EB- and 86 BPA-associated genes with identified DMRs were mostly related to cell-to-cell signaling and interaction, cell-mediated immune response, and cellular growth and proliferation, with cancer as the major disease outcome. Specifically, most of these treatment-associated genes were related to 2 oncogenic rapid signaling molecules AKT and ERK1/2. AKT is upregulated in prostate tumors when compared to benign prostatic tissues and is associated with a stage-dependent increase in tumor cell proliferation. Its activation has also been shown to suppress androgen/androgen-receptor induced apoptosis in PC3 cells. Unlike AKT, expression of ERK1/2 is similar between PCa and normal tissues, yet its activation is associated with increased cell proliferation in PCa tissues. Moreover, the concomitant activation of Akt and Erk1/2 promotes PCa cell growth and tumorigenicity in a rodent PCa model. Taken together, these reports suggest a possible role for the identified treatment-associated genes in the malignant transformation of rat prostate epithelial cells.

Interestingly, functional connectivity of the 7 candidate genes was linked to embryonic stem cell pluripotency, further evidenced by convergence at Sox2. SOX2 is a stem cell marker upregulated in PCa. Overexpression of SOX2 in human DU145 PCa cells promotes tumorigenesis, increases cell proliferation and migration, and reduces apoptosis. Of the other identified genes, TPD52, an oncogene, has been reported to be upregulated in high-grade PIN and PCa. CREB3L4 and WNT10B have also been found to be upregulated in prostate tumors. On the contrary, the expression of PITX3 and CHST14 has yet to be studied in PCa tissues. While only 4 of the 7 genes have been functionally characterized in PCa cells and tissues, all genes, except PAQR4, are related to stem cell proliferation, cancer cell proliferation migration, and/or tumorigenesis (reported bio-functions are summarized in Table 3).

We have previously shown that prostate stem/progenitor cells isolated from the human adult prostate are more proliferative and able to retain their “stemness” properties when treated with E2 or BPA. Suppression of a class of non-coding RNAs, the small nucleolar RNA C/D box, associated with distinct histone modifications has been postulated as a mechanism causing a slow-down in biosynthesis and stem cell differentiation, thus extending the proliferation phase of these cells. In concordance, our current study showed that neonatal exposure to xenoestrogens led to the dysregulation of a panel of stem cell function-related genes through DNA methylation in the adult prostate, lending further credence to this hypothesis. In this regard, neonatal exposure to xenoestrogens may allow for an extended proliferation phase and increase the number of stem/progenitor cells in the prostate. As such, these cells are known targets for carcinogen-induced malignant transformation in adult-life, leading to increased cancer susceptibility.

In addition to their reported involvement in PCa, of clinical relevance, we report here that the expression of the novel 7-gene set is associated with shorter recurrence-free survival of PCa patients in the TCGA cohort. These patients were also mostly diagnosed with prostate tumors of high Gleason score (≥7). Of interest, we also found that 4 out of these 7 genes were differentially methylated between PCa tissues and their adjacent normal tissues in the cohort reported by Kim et al. Thus, this gene signature may have diagnostic/prognostic value by serving as biomarkers for predicting disease initiation and progression.

Although neonatal exposure to either EB or BPA increased PCa risk in the rodent model, their DNA methylation mediated effects on the prostate are different. There were more EB-associated than BPA-associated genes identified from our promoter methylation array. Genes with methylation changes validated were also mostly EB-associated. Moreover, the EB-associated genes linked to AKT and ERK1/2 as the major signaling molecules, whereas the BPA-associated genes linked to TGFβ, an oncogenic molecule upregulated in PCa. This indicates that they activate different signaling pathways to increase PCa risk. Furthermore, EB-treatment dysregulated Pitx3, Wnt10b, Paqr4, and Creb3l4 in PND90 prostate tissues through promoter methylation changes, but BPA exposure had no significant effects. This suggests that EB had a stronger effect on the prostate than BPA, which could in part be related to a markedly higher dose of EB utilized (2500 µg/kg BW) when compared to the lower environmentally relevant dose of BPA (10 µg/kg BW). Alternatively, it has been reported that the action of BPA can be independent of the classical estrogen receptors (ERα and ERβ) and mediated via estrogen-related receptor γ (ERRγ) and G protein-coupled receptor 30 (GPER1), to name a few. However, this does not explain why the expression of Acrbp, Chad, Osbpl6, Phpt1, and Tmem27 did not correlate to their confirmed promoter methylation status, nor why the expression of Paqr4 was unaffected by 5-aza treatment in NbE-1 and AIT cells. In human prostaspheres, E2/BPA-mediated small nucleolar RNA box C/D suppression was shown to be associated with the occupancy of histone marks the trimethylated histone 3 at lysine 9 and lysine 27, and the loss of trimethylated histone 3 at lysine 4, rather than DNA methylation. Thus, EB and BPA may dynamically alter the PND90 prostate transcriptome through other epigenetic mechanisms not yet explored.

In summary, we have identified a set of 7 genes (Pitx3, Wnt10b, Paqr4, Sox2, Chst14, and Tpd52, and Creb3l4) with differential methylation at the promoter region, in PND90 dorsal prostate tissue of rats neonatally exposed to EDCs. The methylation pattern of these 7 genes was inversely correlated to gene expression. Cell-based studies using 5-aza-cytidine-treated normal (NbE-1) and cancerous (AIT) prostate epithelial cell lines further confirmed that gene expression was regulated by DNA methylation, in 6 of the 7 genes. To provide clinical relevance, 4 of the 7 genes (PITX3, WNT10B, PAQR4, and TPD52) have been reported to be differentially methylated between PCa and normal adjacent tissues. Moreover, expression of the 7 genes is associated with the recurrence-free survival of PCa patients. Taken together, we conclude that DNA methylation changes in our novel 7-gene signature may be an epigenetic signature associated with increased cancer susceptibility in the adult gland due to early-life exposure.

Materials and methods

Animal housing and treatments

All animal treatments in this study were in accordance with the accepted standards of humane animal care and approved by the Animal Use Committee at the University of Illinois. Pregnant Sprague-Dawley (SD) rats used on postnatal day (PND) 90 were from Zivic-Miller Laboratories (Pittsburgh, PA) and shipped on gestational day 12. Animal husbandry and conditions were controlled to avoid inadvertent BPA and phytoestrogen exposures, as reported previously.

A schematic diagram of the experimental and analytical procedures is outlined in Figure 1. Pregnant dams were monitored and the day of birth was designated as PND0. Anogenital distance was used to segregate male and female pups. Male pups were randomly assigned to one of 3 treatment groups, n = 30/group: (a) controls given tocopherol-stripped corn oil vehicle alone (Ctrl), (b) 25 µg EB/pup or 2,500 µg EB/kg BW (Sigma-Aldrich Chemical Co., E8515), or (c) 0.1 µg BPA/pup or 10 µg BPA/kg BW (Sigma-Aldrich Chemical Co., 239658). To avoid litter effects, pups from each litter were equally divided among treatment groups and toe clipped for permanent identification. Pups were subcutaneously injected in the nape of the neck with the respective steroids on PND1, 3 and 5, weaned on PND21, and siblings were housed until PND90. The dose and route of EB and BPA was based on our published study and utilized to maintain continuity of results. We had previously characterized the BPA pharmacokinetics in PND3 neonates following a single subcutaneous injection of 10 µg/kg BW of the endocrine disruptor. The mean sera unconjugated BPA level reached a Cmax of 1.77 ng/ml in 0.5 h but rapidly declined to 0.7 ng/ml and 0.54 ng/ml at 1 h and 2 h, respectively, post-injection. These findings suggest a rather rapid rate of bio-degradation of BPA in PND3 neonates, making it unlikely that BPA will reach higher levels upon multiple injections (PND1, 3, 5). However, since we did not have actual measurements of circulating BPA in the neonates during this period, this conclusion may simply be a conjecture. Yet, it should be noted that these levels and those reported before in rodent neonates are well within the range reported for human developmental exposures.

Genomic DNA extraction from dorsal prostate

Genomic DNA was extracted from the rat dorsal prostate using DNeasy Blood & Tissue kit (Qiagen, Valencia, CA) and in the presence of RNase A (Qiagen 19101) according to the manufacturer's protocol. All reagents were supplied in the kit. In brief, the tissue was lyzed in Buffer ATL with proteinase K at 56°C until the tissue was completely lyzed. The lysate was incubated with 100 mg/ml RNase A and mixed thoroughly with equal portions of Buffer AL and 100% ethanol. The mixture was loaded on to a DNeasy Mini spin column, and washed with Buffer AW1 and Buffer AW2 consecutively by centrifugation. Genomic DNA was eluted with DEPC water diluted Buffer AE.

Methylated CpG island Recovery Assay (MIRA)-assisted Methylation array analysis

A total of 5, 4, and 5 PND90 prostate genomic DNA samples from EB, BPA, and Ctrl group, respectively, were used to perform the MIRA-assisted genome-wide methylation analysis. In brief, 4.5 µg genomic DNA was sonicated by Bioruptor (Diagenode, Belgium) to obtain 300–500 bp fragments. The methylated DNA fraction was enriched using MethylCollector Ultra kit (ActiveMotif, 55005) via MIRA-assisted methylated DNA enrichment with proper enrichment controls. It should be noted that the MIRA-, also known as MBD2-, assisted procedure enriches CpG-dense sequences that are likely in CpG islands. The methylation-enriched DNA and its corresponding Input were then amplified using GenomePlex WGA kit (Sigma-Aldrich, WGA2), labeled with Cy3/Cy5 dyes, and co-hybridized in Rat ChIP 385K Promoter 2 array set (Roche NimbleGen, Madison, WI) according to the manufacturer's protocol. The arrays have ∼780,000 probes with probe sizes ranging from 50–70 mers in length and median probe spacing of 105 bp, covering ∼16,000 regions primarily in promoter regions of ∼23,000 transcripts from known genes. These regions principally contain CpG-dense DNA sequences. The probes are typically placed at 4,500 bp upstream and 1,125 bp downstream from the transcriptional start site (TSS). The estimated coverage approximates 97Mb and 3.7% coverage of the rat genome. Thus, the design of the experiments was to focus on CpG rich regions within or near known or predicted gene promoters. This design has the advantage of increasing the hit rate and reducing noise, but does have the deficit of biasing toward gene promoters and having low coverage of the overall rat genome and DMRs not assisted with gene promoters.

Methylation array data analysis

Data analysis was performed using R and Bioconductor packages. The log-ratios of Cy5 to Cy3 (M values) were first normalized using GC-loess normalization. The average M value was then calculated for each probe across all samples in the same group. The significantly enriched regions were identified using sliding window analysis for each group separately. In this analysis, a window of 750 bps around each probe was considered at a time and RandomSet statistic was calculated on the average M values of all probes in the window to estimate the enrichment of the probes in the window when compared to the background. Since there is no enrichment of immunoprecipitated DNA among the probes in the window, the RandomSet statistic, which is defined as the average of the average M values in the window, is approximately distributed as the normal random variable. The P-value of each window was obtained based on the RandomSet statistic and adjusted by the false discovery rate. Windows that were significantly enriched (P < 10−5) in any of the groups (EB, BPA, and Ctrl) were mapped to genes and were further filtered, where differentially methylated genes were selected based on the difference of mBar values of treatment group (EB/BPA) when compared to Ctrl and at different cutoff levels of P-value (Fig. 1). Using this approach, a methylation bar plot was generated for each promoter region in the 3 groups (Fig. 2). In a given promoter region, the difference in height of the bars (mBar) between treated and Ctrl groups represents the difference in probe intensity, i.e., the difference in promoter methylation status. When compared with the control, a positive and negative mBar difference in the treated group indicates promoter hyper- and hypo-methylation, respectively. Primers were then designed based on these differentially methylated regions for validation using bisulfite sequencing analysis.

Bisulfite PCR sequencing analysis

Based on the methylation array bar plot, significant differential methylated regions flanking the gene transcription start sites were selected for methylation validation using bisulfite sequencing. Primers used for bisulfite sequencing were designed with MethPrimer. (Supplemental Table S6) Genomic DNA was bisulfite modified using an EZ DNA Methylation kit (Zymo Research, D5001). In brief, 500 ng genomic DNA diluted with M-Dilution buffer was incubated at 37°C for 15 min, mixed with CT reagent and incubated for 16 hours at 50°C. The reaction was then mixed with Binding buffer, loaded to a spin column, washed, incubated with desulphonation buffer at room temperature for 20 min, and eluted with 40 µl elution buffer. Bisulfite PCR was performed using 2 µl bisulfite modified DNA and Platinum Taq DNA Polymerase (Invitrogen, 10966026) in a 25 µl reaction according to the manufacturer's protocol. Following 40 cycles of PCR amplification with the annealing temperature at 55°C, the amplicons were gel-purified and TA-cloned into pGEMT-easy vector (Promega, A1360). Plasmids from a single E. coli colony were amplified using TempliPhi DNA amplification kits (GE Healthcare, 25640010) and sequenced (Macrogen USA, Rockville, MD). The methylation status of each CpG site was analyzed using BiQ Analyzer.

To validate the promoter methylation status of genes identified from the PND90 prostate microarray analysis, the same amount of DNA from each animal of the same treatment group was mixed for bisulfite sequencing with ∼12 clones per group. To further compare methylation pattern of promoter DMR of the 7 selected gene candidates, bisulfite sequencing analysis on 6 samples each from individual animal was performed with ∼8–12 clones selected per sample.

Treatment of NbE-1 and AIT cells

Rat prostatic epithelial cell lines, the immortalized normal prostatic epithelial cell line NbE-1 and the tumorigenic cell line AIT were maintained in DMEM/F12 medium (Invitrogen, 11330057) supplemented with 5% fetal bovine serum (HyClone, 35-010-CV), 1X insulin-transferrin-selenium (Invitrogen, 41400045), 1 mM sodium pyruvate (Invitrogen, 11360070), and 100 mM minimum nonessential amino acids (Invitrogen, 11140076). The cells were seeded and treated with 0.5 µM or 1.0 µM 5-aza-2-deoxycytidine (5-aza-dC; Tocris, 2624), a DNA methylation inhibitor, every 2 d for 8 d as previously described. Treatment with 0.1% dimethylsulfoxide (DMSO; Sigma D4540) was used as control (Ctrl) in each experimental set. Cell lysates were stored in TRIzol reagent for RNA extraction.

RNA extraction and real-time RT-qPCR

Dorsal prostate tissue was homogenized in TRIzol reagent (Invitrogen, 15596) using Precellys 24 homogenizer (Bertin Technologies, France) and extracted according to the manufacturer's protocol. One microgram total RNA was reversed transcribed using SuperScript III Reverse Transcriptase (Invitrogen, 18080085) in a 20 µl reaction. Primers for real-time RT-qPCR were designed using Primer-BLAST program (Supplemental Table S7). qPCR analysis was performed in a 7500 Fast Real-Time System (ABI, Foster City, CA) in triplicate in a 20 µl reaction. Target gene expression was normalized against the individual endogenous Rpl19 and relative gene expression was calculated using the 2−ΔΔCt method. Similarly, total RNA from treated cells was extracted using TRIzol reagent (Invitrogen, 15596-018) according to the manufacturer's protocol and reversed transcribed for qPCR analyses as previously described.

Functional connectivity analysis

Ingenuity Pathway Analysis (IPA; Qiagen, Redwood City) was performed as previously described for analyzing, integrating, and interpreting the data generated from the methylation array.

Association with the recurrence-free survival of PCa patients

RNAseq data (RNAseqV2, level 3) from 497 PCa subjects as well as their associated clinical data (Biotab) were downloaded from the TCGA database on Feb 1st 2016. Expression levels of the 7 genes were pulled from normalized RNAseq gene results. The data were variance stabilizing transformed before dichotomizing into 2 groups by K-means clustering analysis. Recurrence-free survival was determined based on “days to new tumor event after initial treatment.” Survival analysis with log rank test as well as hazard ratio were calculated using GraphPad Prism software (La Jolla, CA).

Statistical analysis

Data are expressed as mean ± SEM. There were at least 3 individual experimental sets for gene expression analysis. Bisulfite sequencing (BS) analysis of the 7 selected gene candidates was performed on 6 individual samples, with ∼8–12 clones per group. One-way ANOVA and Tukey test were used for comparison among groups. P < 0.05 was considered as statistically significant.