Low-Dose Bisphenol A in a Rat Model of Endometrial Cancer: A CLARITY-BPA Study.

BACKGROUND: Bisphenol A (BPA) is known to be biologically active in experimental models even at low levels of exposure. However, its impact on endometrial cancer remains unclear.
OBJECTIVES: This study aimed to investigate whether lifelong exposure to different doses of BPA induced uterine abnormalities and molecular changes in a rat model.
METHODS: Sprague-Dawley rats were exposed to 5 doses of BPA [0, 25, 250, 2,500, or 25,000μg/kg body weight (BW)/d] or 2 doses of 17α-ethynylestradiol (EE2) (0.05 and 0.5μg/kg BW/d) starting from gestational day 6 up to 1 y old according to the CLARITY-BPA consortium protocol. The BW, uterus weight, and histopathology end points of the uteri were analyzed at postnatal (PND) day 21, 90, and 365. Estrous cycling status was evaluated in PND90 and PND365 rats. Transcriptomic analyses of estrus stage uteri were conducted on PND365 rats.
RESULTS: Based on the analysis of the combined effects of all testing outcomes (including immunohistological, morphological, and estrous cycle data) in a semiblinded fashion, using statistical models, 25μg/kg BW/d BPA [BPA(25)], or 250μg/kg BW/d BPA [BPA(250)] exerted effects similar to that of EE2 at 0.5μg/kg BW/d in 1-y-old rats. Transcriptome analyses of estrus stage uteri revealed a set of 710 genes shared only between the BPA(25) and BPA(250) groups, with 115 of them predicted to be regulated by estradiol and 57 associated with female cancers. An interesting finding is that the expression of 476 human orthologous genes in this rat BPA signature robustly predicted the overall survival (p=1.68×10-5, hazard ratio=2.62) of endometrial cancer patients. DISCUSSION: Lifelong exposure of rats to low-dose BPA at 25 and 250μg/kg BW/d altered the estrous cycle and uterine pathology with similarity to EE2. The exposure also disrupted a unique low-dose BPA-gene signature with predictive value for survival outcomes in patients with endometrial cancer. https://doi.org/10.1289/EHP6875.

Introduction

Bisphenol A (BPA) is used extensively in the manufacturing of plastic consumer products. It is a key building block of polycarbonate-based plastics and vinyl ester resin (Lorber et al. 2015). As a result of growing demand among developing countries, the global BPA production is expected to reach 10.6 million metric tons by 2022 (Research and Markets 2016).

It is now known that BPA leaches out from a wide variety of consumer products, making it ubiquitous in the environment (Rubin 2011). Being estrogenic, it may act as an endocrine-disrupting chemical (EDC). Epidemiological studies have linked BPA exposure to a number of chronic human diseases (Rochester 2013). Studies in animal models support a causal relationship between BPA exposure and the pathogenesis of diseases (Ma et al. 2019; Prins et al. 2019; Wazir and Mokbel 2019). As a result, the U.S. Food and Drug Administration banned the use of BPA in food and beverage containers intended for babies and infants in December 2012 (U.S. FDA 2012).

Although urinary levels of BPA in the U.S. population have declined during the last decade (CDC 2019), recent data showed the chemical was still detectable in 95.7% of randomly selected urine samples in the population (Lehmler et al. 2018). Considering its short half-life (6 h), it can be concluded that the general U.S. population is still continuously exposed to BPA today (Arya et al. 2020; Martínez Steele et al. 2020; Teeguarden et al. 2005; Völkel et al. 2005).

The U.S. Environmental Protection Agency (U.S. EPA) has set the “no observed adverse effect level” (NOAEL) and the “lowest observed adverse effect level” (LOAEL) for BPA at body weight (BW)/d and BW/d in mice and rats, respectively (U.S. EPA 2010). Of concern are the findings from multiple “low-dose” studies that reported the development of adverse health effects in rodents exposed to BPA at doses lower than the NOAEL (Prins et al. 2019; Sekizawa 2008; vom Saal and Hughes 2005). When a wide range of BPA doses was tested, the response did not follow a linear dose–response relationship. Instead, BPA showed a nonmonotonic dose–response pattern, having a greater response at the lower-dose range in in vitro (Vandenberg 2014; Villar-Pazos et al. 2017), in vivo (Jenkins et al. 2011; Leung et al. 2017), and human studies (Awada et al. 2019; Derakhshan et al. 2019; Khalil et al. 2014; Li et al. 2020; Pinney et al. 2017). These findings suggest BPA may have multiple modes of action (Wetherill et al. 2007).

Although many studies had shown BPA as an EDC with profound health effects (Alonso-Magdalena et al. 2012), in 2018 the U.S. Food and Drug Administration (FDA) and National Toxicology Program (NTP) reported the data from The Consortium Linking Academic and Regulatory Insights on BPA Toxicity (CLARITY-BPA) core study and deduced that BPA is safe at the current levels occurring in foods (U.S. FDA 2018). The FDA/NTP argued that BPA treatment groups were not dose responsive and did not demonstrate a clear pattern of consistent responses (Camacho et al. 2019; NTP 2018). Subsequently, such claims were criticized by experts in the field, who stated that the FDA failed to recognize the nonmonotonic dose–response relationships of the data, and address other shortcomings of the study (Vandenberg et al. 2019; Vom Saal 2019).

Most BPA-related cancer research has been centered on breast cancer (Binder et al. 2018; Seachrist et al. 2016; Shafei et al. 2018). Studies on other female cancers such as ovarian and cervical cancer are limited (Ma et al. 2015; Oral et al. 2016; Park et al. 2009). For endometrial cancer, only one cohort study has been reported. It found patients with premalignant endometrial hyperplasia and postmenopausal endometrial cancer had significantly lower serum BPA levels than the healthy subjects (Hiroi et al. 2004), which is an unexpected association that needs future follow-up studies to confirm. In contrast, chronic exposure of mice to low-dose BPA led to an aberrant epithelial proliferation in their uteri (Neff et al. 2019). Furthermore, under in vitro settings, BPA induced epithelial-to-mesenchymal transition in a human endometrial carcinoma cell line through augmentation of cyclooxygenase-2 expression (Wang et al. 2015) and promoted cancer cell proliferation via both epidermal growth factor-dependent and -independent pathways (Yaguchi 2019). However, whether lifelong exposure to a specific dose of BPA can increase the risk of endometrial cancer is currently unclear. This question is of particular importance because uterine growth and development continues after birth and throughout life, although Mullerian organogenesis is largely completed at birth in rodents (Cooke et al. 2013; Kobayashi and Behringer 2003). Thus, early-life exposure of the uterus to EDCs has been shown to impair adult uterine functions in mice (McLachlan et al. 1982; Newbold et al. 2004, 2009) and rats (Bitencourt et al. 2019; Bosquiazzo et al. 2013), and induce uterine preneoplastic and neoplastic lesions in mice (Newbold et al. 1990, 2001, 2007).

In this study, Sprague-Dawley rats were exposed to a wide dose range of BPA ( BW/d to BW/d) starting from gestational day (GD) 6 until specific terminal end points [postnatal day (PND) 21, 90, and 365], using the same Good Laboratory Practice (GLP) protocol published by the CLARITY-BPA consortium (Heindel et al. 2015). We assessed potential disease risks based on estrous cycling, histological and immunohistochemical (IHC) data for multiple time points, and transcriptomic data at the terminal end points. We compared these data to 1 negative-control group and 2 positive-control groups. The clinical relevance of a unique BPA-specific gene signature was further determined using the survival and RNA-Seq data from The Cancer Genome Atlas (TCGA) endometrial cancer cohort.

Materials and Methods

CLARITY-BPA Consortium Study

This study was one of the CLARITY-BPA Consortium Grantee Studies (NTP n.d.), a research program established between the National Institute of Environmental Health Sciences (NIEHS) and the National Center for Toxicological Research (NCTR) of the FDA to evaluate in vivo effects of chronic exposure to BPA over a broad dose range on different tissues or organs. Detailed descriptions of the test articles used, study material evaluations (diet, chemicals, drinking water, cage and bedding leachates), general study design, animal treatments, and animal allocations to the study conducted at NCTR and other grantee studies can be found in Heindel et al. (Heindel et al. 2015).

In this study, the animal dosing experiment was carried out at the Association for Assessment and Accreditation of Laboratory Animal Care (AALAC)-accredited NCTR facility according to FDA GLP regulations. All procedures were approved by the NCTR Institutional Animal Care and Use Committee. In brief, BPA [CAS # 80–05–7, TCI America; catalog # B0494, Lot # 111,909/AOHOK (air-milled)] was administered to Sprague-Dawley rats (from the NCTR breeding colony, Sprague-Dawley/CD23/NctrBR) by oral gavage from gestational day (GD) 6 through the start of labor and then by oral gavage to pups from PND 1 ( 0) until the experimental end point at PND21, PND90, or PND365. BPA doses were 2.5; 25; 250; 2,500; and BW/d. One pup per litter was used as the experimental unit. Ten pups from were used for any outcome analyses. A vehicle (0.3% carboxymethylcellulose, Sigma-Aldrich; catalog # C5013, Lot # 041M0105V) control group (CTL) and 2 positive-control doses (0.05 and BW/d) of (EE2) were included (CAS # 57-63-6, Sigma-Aldrich; catalog # E4876, Lot # 071M1492V, purity).

Sample Collection

All uterus samples were collected by the NCTR laboratory according to their published protocol (Camacho et al. 2019). Briefly, each uterus was weighed and divided longitudinally into 2 halves. The right half of the uterus was fixed in 10% normal buffered formalin for 24 h and then stored in 70% ethanol at room temperature. The left half was cut transversally into 3 segments: the one-third from the proximal end closest to the left ovary was flash-frozen in liquid nitrogen, whereas the remaining 2 segments from the middle and distal portions were individually embedded in OCT embedding medium (Fisher Scientific) and frozen on dry ice. Additionally, when rats reached PND90 or PND365, the estrous cycling status of each rat was determined via examination of vaginal cytology (Cora et al. 2015), and uteri were collected only at estrus. If a rat failed to reach estrus at the end of the fifth day of daily checking, it was euthanized by carbon dioxide asphyxiation and its uterus was harvested. The final cycling stage was recorded at the sacrifice.

Group Data Decoding

When the samples were received from the CLARITY-BPA program, the dosing information for each treatment group was not disclosed to us (the grantees) before completion of data collection and uploading of the data to a secured NTP Chemical Effects in Biological Systems database. All experiments were conducted in a single-blind fashion. After all raw data, including histology data, IHC data, and RNA-sequencing counts, were uploaded to the NTP database and were locked without further modifications, the dose data were decoded and further analyses were carried out accordingly. All methods and results presented in this study were based on the decoded values that represent concentrations of all compounds used.

Histological and IHC Analyses

Five () sections were cut from the formalin-fixed paraffin-embedded (FFPE) samples, and the fifth section was stained for hematoxylin and eosin (H&E) staining and scored for hyperplasia, squamous metaplasia, and other endometrial abnormalities such as endometrial breakdown and inflammatory response. H&E staining was done for sections from three time points (PND 21, 90, and 365). All samples were reviewed by a certified pathologist (A.K.) in a blinded fashion based on representative H&E-stained sections.

Similarly, IHC analyses were performed on the fourth FFPE section with an antibody for proliferating cell nuclear antigen (PCNA, sc-56, Santa Cruz Biotechnology) and antimouse secondary antibody (PI-2000-1, Vector Labs) at a dilution ratio of 1:1,000, and positive signals were developed via 3,3′-diaminobenzidine oxidation reaction according to our published protocol (Leav et al. 2001). Apoptotic cells were stained in the third section using the terminal deoxynucleotidyl transferase dUTP nick-end labeling (TUNEL) assay (Millipore, ApopTag Peroxidase In Situ Apoptosis Detection Kit S7100). The IHC and the TUNEL signals were visualized and captured using LEICA DM300 microscope with DMC2900 camera and subjectively scored (based on the percentage of positivity) by a certified pathologist (A.K.) in a blinded manner.

Whole-Genome RNA-Seq Analysis

To elucidate the long-term effects of BPA on the uterus, flash-frozen uterus samples () from 1-y-old rats in estrus were selected for genome-wide RNA-Seq profiling. Uterus samples () were lysed in QIAzol buffer (Qiagen) and homogenized using Precellys 24 tissue homogenizer (Bertin Corp.); the total RNA was extracted using a Qiagen RNeasy Lipid Tissue Kit (Qiagen). The RNA quality and quantity were assessed using an Agilent Bioanalyzer and NanoDrop ND-1000 spectrophotometer (Thermo Scientific), respectively. Samples with an RNA integrity number (RIN) greater than 9 were used for sequencing. RNA libraries were prepared with the TruSeq RNA Library Preparation Kit according to the manufacturer’s protocol (Illumina) and were sequenced with the HiSeq 1000 sequencing system in a sequencing core at the University of Cincinnati.

RNA-Seq data processing and analysis were performed by following the published protocol (Huber et al. 2015). More specifically, sequence reads were aligned to the reference Rattus norvegicus genome (rn5) using the TopHat2 aligner (Trapnell et al. 2009) followed by assessing the sequence’s quality with FastQC and RNA-SeQC based on nucleotide composition bias, PCR bias and GC bias, sequencing saturation, mapped reads distribution, coverage uniformity, strand specificity, and transcript level RNA integrity. (Andrews 2010; DeLuca et al. 2012). Aligned reads were sorted according to genome coordinates, and samples were assigned to specific ID/group with Picard from Broad Institute (http://broadinstitute.github.io/picard), as a preprocessing step for performing quality control (QC) with RNA-SeQC. Reads aligning to each known transcript were counted using R packages GenomicFeatures and GenomicAlignments (version 3.4.4, R Development Core Team) (Anders et al. 2013). Only genes with more than 5 reads in at least n samples, where n is the number of samples in the smallest groups, were analyzed. -Values were adjusted for multiple comparisons testing using the false discovery rates (FDR) (Storey and Tibshirani 2003). Gene expression profiles () in the heat map were clustered using the Bayesian infinite mixture model (Freudenberg et al. 2010; Medvedovic and Sivaganesan 2002), and the groups were clustered using hierarchical clustering with the average method. Multidimensional scaling analyses were conducted to visualize the level of similarity between individual samples. The enrichment analysis was performed with knowledge-based Ingenuity Pathways Analysis (IPA) (Qiagen; www.qiagen.com/ingenuity). Fisher’s exact test was used to measure whether there was a statistically significant overlap between the 710 unique genes common to both groups and the genes regulated by a known transcriptional regulator listed in the Ingenuity Knowledge Base (from IPA). ToppFun (ToppGene Suite; https://toppgene.cchmc.org/enrichment.jsp), was used to identify any genes specific to a certain disease (Chen et al. 2009). RNA-seq data were deposited in the National Center for Biotechnology Information Gene Expression Omnibus database with accession number GSE139308.

TCGA Survival Analyses

All analyses were carried out using packages in R (R Development Core Team; http://www.r-project.org). Harmonized RNA-Seq and clinical data for TCGA uterine corpus endometrial carcinoma samples were downloaded from TCGAbiolinks on 26 October 2018. The data for 548 patient samples were then normalized with ComBat (Johnson et al. 2007) to remove batch effects. Rat to human gene conversion was done with the orthologous conversion function in biomaRt (Durinck et al. 2009). Human transcripts with low variability across the TCGA cohort were excluded from further analysis. The exclusion process was arbitrary to keep most variable transcripts but not the transcripts with low variance. As a result, 50,000 most variable transcripts were selected for further analysis. Hierarchical clustering was performed to dichotomize the cohort into two groups based on the difference of expression level of the human orthologous genes without an a priori cutoff. Overall and recurrence-free survival between the groups were compared with a Kaplan–Meier plot with a log-rank test and with Cox proportional hazard models (with adjustments for race and age at initial pathologic diagnosis). Recurrence-free survival was determined based on “days to new tumor event after initial treatment.” To establish the statistical significance of the selected group of genes used for clustering human samples, the p-value for difference in survival was compared to the distribution of p-values obtained by clustering samples based on randomly selected sets of genes of the same size. The empirical p-value was calculated as the proportion of 10,000 random samples that resulted in more significant p-values than the p-value observed for the original set of genes.

To further define the characteristics between the two groups, clinical features including age, race (Caucasian vs. African American), menopausal status (pre- vs. postmenopausal, as defined by TCGA consortium), obesity [yes/no obese, obese body mass index ], estrogen receptor level (high/low, ESR1 transcript level), cancer grade (high/low, 3 or presented as “high” in the original TCGA data), cancer status (yes/no cancer-free), and cancer histology (yes/no serous type) were compared using a Student’s t-test for continuous variables or Pearson’s chi-square test with contingency tables for categorical variables using GraphPad Prism 6.0.

Statistical Analyses

Numerical biomarkers were summarized using the mean and standard deviation, and binary or categorical biomarkers were summarized using frequency (%). The association of a biomarker to the treatment group was assessed using an F test from the ANOVA model. Post hoc means were compared between groups under the ANOVA model and adjusted for multiple comparisons using Bonferroni’s method. For a binary or categorical biomarker, its association with the group effect was assessed using Fisher’s exact test. The same Fisher’s exact test was also used to compare frequencies between treatment groups after the p-value was adjusted for multiple comparisons using Bonferroni’s method.

To investigate which BPA group was more similar to EE2 or CTL based on observations (i.e., predictors) of estrous status, histology (yes/no proliferation), and IHC staining for nuclear PCNA and apoptosis (percentage of positive staining), we used two predictive modeling methods, the generalized logistical regression model (GLRM) and Fisher’s linear discriminant analysis (LDA). To minimize potential analytical bias in the process, the BPA doses were blinded to the statistician.

For the GLRM method, the treatment group was defined as a categorical dependent variable, and estrous status, histology, PCNA, and apoptosis were treated as predictors of independent variables. From the GLRM, a propensity score was predicted for each subject, or more specifically, the probability of being categorized as EE2 (or CTL) was predicted. The summarized propensity scores from a treatment group were then used to determine if the group was more likely to be similar to EE2 or CTL. Intuitively, a treatment group would be considered closer to the EE2 (or CTL) group if its propensity score of being EE2 (or CTL) was higher than that of the CTL (or EE2).

For the LDA method, the linear combination of the predictors that could maximize the distance (or discrimination) among the groups was determined. Using the LDA method, each of the treatment groups was developed and used to assign a score from the first canonical component of LDA. The difference in scores relative to EE2 (or CTL) could be used to identify if the group were more “close” or “similar” to EE2 (or CTL). In particular, a treatment group was considered more similar to EE2 (or CTL) if its difference was smaller than that of the CTL (or EE2). All statistical computation was performed using SAS (version 9.3; SAS Institute Inc.). were considered statistically significant.

Results

All data presented in figures or tables are now fully decoded.

Continuous Exposure to BPA on Body and Uterus Weight

To examine whether BPA at any dose could be toxic, we first compared BW in the BPA-treated rats with the CTL. As shown in Figure S1, there was no difference in BW in any of the rats treated with BPA or EE2 at any time points measured (). Similar results were observed for uterus weight, except for a significantly higher observed uterus weight () in the high-dose EE2(0.5) group in comparison with all other groups at PND21 (Figure S1). Similar results were obtained when the uterus weight data were normalized with the BW (Figure S1).

Histological Analysis of Uteri of 1-Y-Old Rats Treated with EE2 or BPA

We next examined the histology of the uterus samples collected at PND21, PND90, and PND365. None of the samples at PND21 showed abnormal histology (data not shown). However, half (5 out of 10) of the uterus samples from the EE2(0.5)-treated group showed squamous metaplasia at PND90 (, Figure S2). Interestingly, after a 1-y exposure, 2 out of 9 samples (22.2%) in the EE2(0.5) group; 1 out of 9 (14.3%) in BPA(25); and 1 out of 10 (10%) in BPA(25,000) showed squamous metaplasia, even though these differences were not statistically significant (, Figure S3), and none was observed in the CTL group.

Proliferation and Apoptosis in Uteri of 1-Y-Old Rats Treated with EE2 or BPA

Estrogens are known to induce uterine epithelial cell proliferation (O’Brien et al. 2006; Winuthayanon et al. 2010). We therefore immunostained the uterus samples with an anti-PCNA antibody as an indicator of cell proliferation, quantified the staining, and compared the results between groups. In PND21 samples, we observed slightly but nonsignificantly higher numbers () in PCNA-positive cells in both the EE2(0.5) (high-dose) and BPA(250) groups in comparison with the CTL group (Figure S4). For all other time points, there were no significant differences in PCNA positivity ( and 0.35, Figure S4).

Similarly, because apoptosis is normally observed in cycling endometrium (Dharma et al. 2001), we examined whether BPA could affect this process. In PND21 rats, high-dose EE2 or either doses of BPA(2.5) or BPA(2,500) exhibited more apoptotic cells than the CTL or any other EE2 and BPA-treated rats (Figure S5). However, in PND90 and 365 rats, we did not observe any significant difference in apoptosis between the CTL and any BPA-exposed groups () (Figure S5).

Estrous Cycling Analysis in PND90 Rats Treated with EE2 or BPA

Because we aimed to collect the uterus at the estrus stage for gene expression analysis, we checked the cycling stage of each rat before it was euthanized at approximately PND90. As shown in Figure 1A, more than 50% of the rats in almost all groups were in the estrus stage except the EE2(0.5) group, but we observed a mixed transition between diestrus and estrus phase in this group ( vs. other groups). A minor abnormal cycling transition was also observed in the BPA(250) group, but it was not statistically significant [against other groups except EE2(0.5), ]. When compared with the CTL, EE2(0.05), EE2(0.5), BPA(250), BPA(2,500), and BPA(25,000) showed different degrees of cycle disruption, but the differences were not significant ().

Figure 1.

Estrous cycling of PND90 and PND365 rats. Female rats were treated with following treatment from gestation day 6 till (A) postnatal day (PND) 90 or (B) PND365. CTL: vehicle only; EE2(0.05): body weight/day ; EE2(0.5): body weight/day EE2; BPA(2.5): BW/d bisphenol A; BPA(25): BW/d bisphenol A; BPA(250): BW/d bisphenol A; BPA(2,500): BW/d bisphenol A; BPA(25,000): BW/d bisphenol A. When rats reached PND90 or PND365, the estrous cycling status of each rat was determined by examining the cellular composition of the vaginal smear, and uteri were collected only at estrus. If a rat failed to reach estrus at the end of the fifth day of daily checking, it was euthanized, and its uterus harvested. The final cycling stage was recorded at the sacrifice. The percentage of each stage was calculated in each group. Fisher’s exact test was used to compare frequencies between treatment groups. * EE2(0.5) vs. any group; @; # EE2(0.5) vs. EE2 (0.05). Please refer to Excel Figure 1 for the raw data. BW, body weight.

Estrous Cycling Analysis in 1-Y-Old Rats Treated with EE2 or BPA

Similar analyses were performed for the PND365 rats. As shown in Figure 1B, the distribution of estrous status in the EE2(0.5) group differed from that of the CTL group and the EE2(0.05) group (). Other comparisons were not statistically significant, even though we noticed that the ratio was not dramatically different in the EE2(0.05), BPA(2.5), BPA(2,500), and BPA(25,000) groups; the ratio was not maintained in the EE2(0.5), BPA(25), and BPA(250) groups; in addition, exposure to BPA(25) and BPA(250), like EE2(0.5), seemed to prolong the estrus stage and affect the normal cycling of 1-y-old rats.

Logistic Regression Analysis of Combined Immunohistological, Histological, and Estrous Cycle Data in 1-Y-Old Rats Treated with EE2 or BPA

To identify which dose of BPA acted in a manner similar to that of EE2, we performed the GLRM method combining all measured parameters collected (including IHC, histology, and estrous cycle data) at PND365 as predictors. As shown in Table 1, both BPA(250) and BPA (25) showed higher propensity scores for being similar to EE2(0.5). We also applied the LDA method to the same data set (Table 2). Similarly, BPA(250) and BPA(25) had smaller distances from EE2(0.5), suggesting that they were more similar to EE2(0.5) than to the CTL. Other groups, including BPA(2.5), BPA(2,500), BPA(25,000), and EE2(0.05) were more likely to be similar to the CTL.

Table 1

Predictive probability of being classified as high-dose EE2/CTL using logistic regression analyses.

	Propensity score
Group	CTL	EE2(0.5)
EE2 (0.05)	0.16	0.06
BPA (2.5)	0.13	0.10
BPA (25)	0.11	0.12
BPA (250)	0.10	0.19
BPA (2500)	0.15	0.10
BPA (25000)	0.13	0.11

Note: Rats were treated as follows from gestational day 6 until postnatal day 365. CTL: vehicle only; EE2(0.05): BW/d (EE2); EE2(0.5): BW/d EE2; BPA(2.5): BW/d bisphenol A (BPA); BPA(25): BW/d BPA; BPA(250): BW/d BPA; BPA(2500): BW/d BPA; BPA(25000): BW/d BPA. In the generalized logistic regression model, the treatment group (BPA, EE2, or CTL) was defined as a categorical dependent variable, and estrous status, histology, PCNA, and apoptosis were treated as predictors of independent variables. A propensity score [the probability of being categorized as EE2(0.5) or CTL] was predicted for each subject. Summarized propensity scores from a treatment group were used to determine if the group was more or less likely to be similar to EE2(0.5) or CTL. BPA, bisphenol A; BW, body weight.

Table 2

Predictive probability of being classified as high-dose EE2/CTL using linear discriminant analyses.

Group	Average canonical score	Diff. against CTL	Diff. against EE2(0.5)
CTL	0.62	—	—
EE2 (0.05)	−1.01	0.02	1.61
EE2 (0.5)	0.60	—	—
BPA (2.5)	0.13	0.49	1.14
BPA (25)	−0.10	0.72	0.90
BPA (250)	−0.56	1.18	0.44
BPA (2500)	0.23	0.38	1.24
BPA (25000)	0.12	0.50	1.12

Note: Rats were treated as follows from gestational day 6 until postnatal day 365. Estrous status, histology, PCNA, and apoptosis were treated as predictors of independent variables. For the LDA method, the linear combination of the predictors that could maximize the distance (or discrimination) among the groups was determined. Using LDA, each of the treatment groups was developed and used to assign an “average canonical score” for the first canonical component of LDA. The difference between scores for a specific treatment group and EE2(0.5) (or CTL) determined if the group was more “close” or “similar” to EE2(0.5) (or CTL). A treatment group was considered more similar to EE2(0.5) (or CTL) if its difference against EE2(0.5) (or CTL) was smaller than its difference against CTL [or EE2(0.5)]. —, no data; BPA(2.5), BW/d bisphenol A (BPA); BPA(25), BW/d BPA; BPA(250), BW/d BPA; BPA(2500), BW/d BPA; BPA(25000), BW/d BPA; BW, body weight; CTL, vehicle only; EE2(0.05), BW/d (EE2); Diff against CTL/EE2(0.5), difference in scores between a particular group against CTL or EE2(0.5), respectively; EE2(0.5), BW/d EE2.

Transcriptomic Analysis of Differential Gene Expression in Uteri of 1-Y-Old Rats Treated with EE2 or BPA

We performed RNA-Seq analysis on selected samples from 1-y-old rats that were at estrus stage ( per group). This analysis allowed us to investigate the effect of BPA under the same hormonal background. Using pairwise comparison, we found that more than 400 genes were differentially expressed relative to the CTL in each of the BPA-treated groups () (Figure S6).

We then performed content-specific Bayesian clustering analysis of the top 1,000 significant genes to determine the similarity among the groups based on the transcriptomic signature. Genes affected by BPA(25) and BPA(250) formed tighter clusters and shared similar nodes in the dendrogram (Figure 2A). On the other hand, EE2 groups (both low- and high-dose) shared more similar gene expression patterns with other groups, including the CTL and other doses of BPA. Using multidimensional scaling analysis to visualize the level of similarity between individual samples based on the expression level of all genes, we observed that genes affected in the BPA(25) and BPA(250) groups were closely related to each other but were distinctively different from the CTL and both EE2 groups (Figure 2B).

Figure 2.

Transcriptomic RNA-Seq analysis of uterus samples at the estrus stage collected after exposure to CTL, EE2, or BPA for 1 y. Female rats were treated as follows from gestation day 6 until postnatal day (PND) 365. CTL: vehicle only; EE2005: BW/d ; EE205: BW/d EE2; BPA2_5: BW/d bisphenol A; BPA25: BW/d BPA; BPA250: BW/d BPA; BPA2500: BW/d BPA; BPA25000: BW/d BPA. Transcriptomic analysis was conducted on uterus samples at estrus () using RNA-Seq. (A) The top 1,000 genes with significant differences (), compared to the CTL, were analyzed with a Bayesian infinite mixture model, and the samples were analyzed with average-linkage hierarchical clustering and presented as a heatmap; rows indicate individual genes and columns indicate individual samples. Treatment groups in columns are color-coded. (B) All genes were subjected to multidimensional scaling analyses to visualize the similarity between individual samples. Samples are color-coded by groups; some groups are circled to show the proximity of each sample in a particular group. BW, body weight.

Knowledge-Based Gene Ontology Analysis of the Transcriptomic Data Derived from BPA(25)- and BPA(250)-Treated Rats

A Venn diagram (Figure 3A) showed that 710 genes were common to both the BPA(25) and BPA(250) groups but did not overlap with others. It is interesting to note that all of the genes were concordantly up- or down-regulated by both BPA doses, and more than one-third of those genes (262 genes) showed a greater fold difference in the BPA(25) group than in the BPA(250) group (Table S1). Using IPA, we found that more than 16% of the 710 genes (115 genes) were predicted to be estrogen-associated genes because “beta-estradiol” was the key upstream regulator ( Table 3 and Figure 3B). Using ToppFun, we identified a subset of genes (57 out of 710 genes) that was significantly () associated with “malignant tumor of the cervix” and “cervical cancer” (Table 4).

Figure 3.

(A) Venn diagram of differentially expressed genes. Female rats were treated as follows from gestation day 6 until postnatal day (PND) 365. CTL: vehicle only; BPA2.5: BW/d bisphenol A; BPA25: BW/d BPA; BPA250: BW/d BPA; BPA2500: BW/d BPA; BPA25000: BW/d BPA. Transcriptomic analysis was conducted on uterus samples at estrus stage using RNA-Seq. Significantly differentially expressed genes () are shown in the Venn diagram. (B) Upstream regulator prediction with Ingenuity Pathway Analysis for a subset of differentially expressed genes that overlapped only between BPA(25) and BPA(250) groups. Dotted arrows represent indirect interactions; solid arrows represent direct interactions. Note: BW, body weight; FDR, false discovery rate.

Table 3

Upstream regulators predicted for 710 unique genes common to BPA(25) and BPA(250) groups.

Name	# Molecules in data set	p-Value of overlap
Beta-estradiol	115	1.84×10^−11
TP53	97	1.02×10^−8
Dexamethasone	105	7.29×10^−8
TGFB1	89	3.89×10^−7

Note: RNA-Seq data indicated that 710 genes were differentially expressed in both the BPA(25) and BPA(250) treatment groups but did not overlap with others. Ingenuity pathway analysis (IPA) was used to predict upstream regulators of those genes. Only IPA-annotated genes (677) were used for the analysis. The overlap -value was calculated using Fisher’s exact test and determined whether there was a statistically significant overlap between the 677 unique genes common to both treatment groups and the genes regulated by a known transcriptional regulator in the Ingenuity Knowledge Base (from IPA). # Molecules in data set, number of genes that overlapped with the submitted data. TP53, tumor protein p53; TGFB1, transforming growth-factor beta 1.

Table 4

Disease prediction by ToppFun.

	ID	Name	Source	p-Value	FDR B&H	Genes from input	Genes in annotation
1	C0700095	Central neuroblastoma	DisGeNET BeFree	1.083×10^−6	5.767×10^−3	93	1655
2	C0007847	Malignant tumor of cervix	DisGeNET Curated	3.185×10^−6	8.481×10^−3	57	883
3	C4048328	Cervical cancer	DisGeNET Curated	5.459×10^−6	9.692×10^−3	60	964

Note: ToppFun performs enrichment analysis based on a hypergeometric model and was used to determine whether the 710 genes differentially expressed in the both the BW/d bisphenol A and BW/d bisphenol A groups were enriched in sets of genes associated with human diseases. Only ToppFun-annotated genes (650) were used for the analysis. The values for enriched annotations were corrected for multiple testing using the Benjamini–Hochberg method. The enriched annotations with corrected p were identified as overrepresentative annotations for the submitted gene set. BW, body weight; DisGeNET, discovery platform containing one of the largest publicly available collections of genes and variants associated with human diseases; FDR B&H, false-discovery rate calculated based on the Benjamini–Hochberg method; Genes from input, number of submitted genes matched with disease genes; Genes in annotation, number of genes associated with a particular disease; Source, database for the disease gene set.

Survival Analysis of the Human TCGA Endometrial Cancer Patients Using BPA-Responsive Genes Identified in BPA25- and BPA250-Treated Rats

To define the clinical relevance of these 710 genes, we performed survival analyses with TCGA data. We first mapped the 710 rat genes into 674 human genes using biomaRt and then removed the genes with the least variation in the TCGA cohort to give 478 genes (Table S2). We then performed modeling-based hierarchical clustering to dichotomize the cohort into two groups based on the expression level of 478 genes (Figure 4A). Survival analyses revealed that patients in group 2 had significantly poorer overall survival [, ], as well as poorer recurrence-free survival (, ) (Figure 4B). When the data were adjusted for age and race, these associations remained significant; specifically, the HR of overall survival decreased to 2.62 (), and the HR of recurrence-free survival increased to 1.91 ().

Figure 4.

Survival analyses of human The Cancer Genome Atlas (TCGA) Uterine Corpus Endometrial Carcinoma (UCEC) data with low dose BPA-targeted genes. (A) Hierarchical clustering of 548 human uterine cancer samples from TCGA UCEC cohort using the 710 genes differentially expressed during estrus in 1-y-old rats treated with 25 or BW /d BPA. The 710 rat genes were mapped into 674 human genes in the UCEC cohort of TCGA using biomaRt. Furthermore, 196 genes were removed due to low variation among samples; 478 genes were used for modeling-based hierarchical clustering and dichotomizing the cohort into two groups based on the difference of expression level of the human orthologous genes without a priori cutoff. The results were presented as a heatmap; rows indicate individual genes, and columns indicate individual samples. Sample features in columns are color-coded; details in (C). (B) Survival analysis of the TCGA UCEC cohort using unique genes identified in BPA(25) and BPA(250)-treated rat uteri. Modeling-based clustering (GIMM) was applied to the human genes to dichotomize the cohort into Groups 1 or 2; and Kaplan Meier curves of (left panel) overall survival () and (right panel) recurrence-free survival () were compared between the two groups using a log-rank test. (C) Detailed column headings from (A). Each row represents one parameter examined. Based on the expression of the 478 human genes, the cohort was divided into 2 groups [group 1 (red) vs. group 2 (black)] using average-linkage hierarchical clustering. American, American Indian or Alaska native; Black, Black or African American; BW, body weight; Endometrioid EA, endometrioid endometrial adenocarcinoma; ESR, estrogen receptor alpha; Free, tumor-free; GA, Genome Analyzer sequencing platform; G1, Grade 1; G2, Grade 2; G3, Grade 3; HiSeq, HiSeq sequencing platform; Mixed EA, mixed serous and endometrioid; Native, Native Hawaiian or other Pacific Islander; Neg, Negative; Pos, Positive; Serous EA, serous endometrial adenocarcinoma; Tumor, with tumor; Un, unknown.

To test that similar results could not be obtained using a random set of genes, we performed 10,000 random samplings of 478 human genes, repeated the survival analyses and calculated the empirical -value. The empirical -value was 0.026 for overall survival and 0.0414 for overall survival when adjusted for age and race. This finding indicated that only 2.6% and 4.14% of the random sampling could yield a significant difference in overall survival without and with covariate adjustment, respectively. However, the empirical -value went up to 0.4749 for recurrence-free survival, indicating that a significant difference could be obtained by random samplings in almost half of the cases (47.5%).

Because the TCGA RNA-Seq data were obtained with two different sequencing platforms (HiSeq vs. Genome Analyzer system), we examined whether the group assignment was confounded by this batch difference. As shown in Figure 4C, both groups showed a similar number of samples sequenced by either platform, and the batch effect was successfully corrected with the ComBat pipeline. To further characterize the patients in group 2 (with poor overall survival), we investigated several clinical parameters. As shown in Figure S7, patients in group 2 were significantly older (average age: 66.4 years of age in group 2 vs. 62.5 in group 1, ), and there was a significantly higher proportion of African Americans in group 2 than in group 1 (27%, 57 out of 298, African Americans in group 2 vs. 19%, 49 out of 179 in group 1, ). When we considered pre- and postmenopausal status, we found that the proportion of premenopausal women in group 2 was significantly smaller than that in group 1 (4%, 8 out of 186, premenopausal in group 2 vs. 9%, 27 out of 294, in group 1, ). Based on BMI (obese, ), patients in group 2 were significantly less obese (50%, 96 out of 191, obese in group 2 vs. 65%, 207 out of 319, in group 1, ). Consistent with the survival analysis, patients in group 2 had more aggressive disease, with fewer patients surviving; a higher number of patients with tumor even after initial treatment (68%, 133 out of 196, tumor-free in group 2 vs. 89%, 290 out of 327, in group 1, ); higher neoplasm grade (92%, 191 out of 208, high-grade in group 2 vs. 40%, 134 out of 335, in group 1, ); and more aggressive histological cancer type (58%, 121 out of 208, serous in group 2 vs. 4%, 15 out of 335, in group 1, ) in their cancer (Figure 4C and Figure S7). Serous cancer was predominantly found in group 2, and gene expression data for ESR1 reflected this (i.e., 89%, 185 out of 208, of the patient samples expressed a low level of ESR1 in group 2 vs. 36%,122 out of 335, in group 1; ) (Figure 4C and Figure S7).

As an auxiliary analysis, we investigated how many genes were highly differentially expressed between the two groups. We first conducted IPA analysis on the 486 human orthologous genes and found that 78 genes shared “” as the top upstream regulator. ToppGene analysis on transcription factor binding-site prediction on these 78 genes did not identify any classical ER-binding elements or known tethering binding sites on their putative promoters (Table 5). Of these 78 genes, 37 of them showed a significant difference in terms of expression level () between the two groups, but only 4 of the 37 genes (TRHDE, CDKN1A, COLEC12, and SLC8A1) showed fold differences greater than 2 (Figure 5A). An interesting aspect is that the same four genes in rats were selectively higher or lower in rats treated with the two low doses of BPA (25 and ), and two of them, Colec12 and Slc8a1, followed a nonmonotonic dose–response pattern (Figure 5B). The expression of Colec12 in both the BPA(25) and BPA(250) groups differed from the CTL, EE2(0.05), EE2(0.5), BPA(2.5) and BPA(2,500) groups (). The expression of Slc8a1 in the BPA(250) group differed from the CTL, BPA(2.5), BPA(2,500) and BPA(25,000) groups (); and the BPA(25) group differed from the CTL and BPA(2.5) and 0.007, respectively.

Table 5

Transcription factor binding site by ToppFun.

	ID	Name	p-Value	FDR B&H	Genes from input	Genes in annotation
1	V$AP2ALPHA 01	V$AP2ALPHA 01	4.451×10^−6	2.698×10^−3	24	195
2	V$HNF3ALPHA Q6	V$HNF3ALPHA Q6	6.889×10^−5	2.087×10^−2	20	173
3	V$AP2GAMMA 01	V$AP2GAMMA 01	1.961×10^−4	3.461×10^−2	21	201

Note: ToppFun performs enrichment analysis based on a hypergeometric model and was used to determine whether the 710 genes differentially expressed in both the BW/d bisphenol A and BW/d bisphenol A groups were enriched in sets of genes associated with a transcription factor binding site. Only ToppFun-annotated genes (650) were used for the analysis. The -values for enriched annotations were corrected for multiple testing using the Benjamini–Hochberg method. The enriched annotations with corrected were identified as overrepresentative annotations for the submitted gene set. FDR B&H, false-discovery rate calculated based on the Benjamini–Hochberg method; Genes from input, number of submitted genes matched with the disease genes; Genes in annotation, number of genes associated with a particular transcription factor binding site; V$AP2ALPHA 01, AP2 alpha; V$HNF3ALPHA Q6, FOXA1; V$AP2GAMMA 01, AP2 gamma.

Figure 5.

Expression of genes in two sets of RNA-Seq data. Comparison of the expression of (A) between two groups in the Cancer Genome Atlas (TCGA) cohort; (B) in different BPA dose groups in the rat study. Female rats were treated as follows from gestation day 6 until postnatal day (PND) 365. CTL: vehicle only; EE2005: BW/d ; EE205: BW/d EE2; BPA2_5: BW/d bisphenol A; BPA25: BW/d BPA; BPA250: BW/d BPA; BPA2500: BW/d BPA; BPA25000: BW/d BPA. Transcriptomic analysis was performed on the uterus samples at estrus stage using RNA-Seq; 710 genes were differentially expressed only in the uteri of 1-y-old rats exposed to BPA25 and BPA250 compared with CTLs. The 710 rat genes were mapped into 674 human genes in the Uterine Corpus Endometrial Carcinoma cohort of TCGA using biomaRt. A total of 196 genes were removed due to low variation among samples; 478 human orthologous genes were submitted to Ingenuity Pathway Analysis (IPA) for predicting upstream regulators; 78 genes shared “” as the top upstream regulator. Within 78 genes, 37 showed significant differences in expression () between the two groups, as defined by hierarchical clustering, but only four (TRHDE, CDKN1A, COLEC12, and SLC8A1) showed fold differences . The expression level of those four genes was compared between the RNA-Seq data from the TCGA cohort and the RNA-Seq data from the rat–BPA study. Data were plotted as the error. Student’s -test was applied to compare two groups (***; ****) for the human TCGA study. For the rat study, ANOVA was applied to compare the gene expression level of a particular BPA dose group with the CTL (*; **). Trhde: BPA250 and CTL were different with ; Cdkn1a: BPA2.5 was different from BPA25, BPA250, BPA25000, and EE2 0.05 with ; Colec12: Both BPA25 and BPA250 were different from CTL, EE2 0.05, EE2 0.5, BPA2.5, and BPA2500 with ; Slc8a1: BPA250 was different from CTL, BPA2.5, BPA2500, and BPA25000 with ; and BPA25 was different from CTL and BPA2.5 with and 0.007, respectively. BW, body weight; CDKN1A/Cdkn1a, cyclin-dependent kinase inhibitor 1A; COLEC12/Colec12, collectin subfamily member 12; SLC8A1/Slc8a1, solute carrier family 8 member A1; TRHDE/Trhde, thyrotropin-releasing hormone degrading enzyme.

Discussion

In this study, we found that lifelong exposure (GD6 to PND365) to low doses of BPA (25 or BW/d), but not higher doses, altered the estrous cycle, associated transcriptomes, and uterine pathology. This conclusion was based on an analysis of the combined effects of all testing outcomes in a semiblinded fashion using two statistical models (logistic regression and linear discriminant analysis). These two low doses [BPA(25) and BPA(250)] shared a similar gene expression signature, which was distinct from other dose groups in unbiased transcriptome analysis. Further bioinformatics analyses uncovered that a subset of those BPA-responsive genes was linked to “,” exhibited a nonmonotonic dose–response pattern, and could be oncogenic in cervical cancer. Differences in the expression of the orthologous human genes were associated with poor overall survival in a TCGA endometrial cancer cohort. Overall, our data suggest that the unique gene signature induced by lifelong exposure to specific doses of BPA might be involved in endometrial cancer growth and progression.

To our knowledge, the endocrine-disrupting effect of BPA on nonreproductive end points of the uterus in a female rat model had not been evaluated thoroughly until this study. The model we report here is unique because it covers five log doses of BPA, ranging from 2.5 to BW/d, with two doses lower than U.S. EPA’s safe human dose (e.g., BW/d) (U.S. EPA 2010). Also, ours is a lifelong exposure model in which the rats were chronically exposed to BPA daily starting from GD6 to PND365. Most of the rat studies reported the effect of BPA on the uterus either within a fixed window of exposure or within a limited dose range (Suvorov and Waxman 2015). Only two of the models conducted similar lifelong exposure in rats but with a different dose range ( BW/d to BW/d or BW/d to BW/d) and with much shorter time points for final assessment (PND70–117) (Ema et al. 2001; Tyl et al. 2002). Similar to our results, neither study found a significant difference in cyclicity, uterus weight, or histopathology. This finding was also in agreement with the data reported by the NTP that BPA exposure groups showed no difference in the estrous cycling of approximately 16-wk-old rats (NTP 2018). However, when evaluating all measured parameters in statistical models, we found that both the BW/d and BW/d doses of BPA were more “estrogenic” than other doses because they produced effects similar to high-dose EE2 ( BW/d).

When we analyzed the transcriptome data for the samples collected at the estrus stage, we identified a specific gene signature, containing 710 genes, that tied both the BPA(25) and BPA(250) groups together and formed unique gene clusters, which were distinct from the EE2 groups. Even though both groups behaved like a high-dose EE2 group in a combined estrous cycling and histopathology analysis, the molecular pathways underlying the effects of BPA (25 and 250) and EE2 in the estrus-stage uterus were quite different. In other words, at low doses, BPA acted similarly to an estrogen to disturb estrous cycling and tissue histopathology but could exert extraestrogenic effects in regulating specific gene expression in the uterus during estrus stage. Notably, more than one-third of the 710 genes (262 genes) showed the greatest fold change in the lower BPA dose group, indicating that those genes were sensitive to minimal BPA exposure (Table S1). Further analyses revealed that a subset of the signature (78 genes) was shared estradiol as the top upstream regulator in the IPA analysis, and another subset (57 genes) was associated with cervical cancer in a disease-prediction analysis, suggesting that some BPA-responsive genes were also estrogen-responsive and/or might be involved in cervical carcinogenesis. Of the 78 genes, 4 genes (CDKN1A, TRHDE, COLEC12, and SLC8A1) were highly differentially expressed in only two low-dose–exposed groups in the rat data and also were presented in the gene signature that predicts the survival of endometrial cancer patients in the human data. An interesting finding is that exposure to EE2 did not alter the expression level of all four genes although they shared as the common upstream regulator. This suggests that the estrogenic signaling exerted by EE2 could be different from the signaling elicited by BPA, especially when the concentration of endogenous estrogens is high, but the detailed mechanism remains to be uncovered.

CDKN1A was shown to be an estrogen-regulated gene in a study with a breast cancer cell line (MCF-7) (Mandal and Davie 2010), and its elevated expression level correlated with better prognosis in endometrial adenocarcinoma patients (Tang et al. 2019). TRHDE and COLEC12 were reported to be estrogen-regulated in in vitro (Zhao et al. 2009) and in vivo studies (Bauer 1988; Schomburg and Bauer 1997), but its role in endometrial cancer remains largely unknown. On the other hand, SLC8A1 expression level did not seem to be affected by estrogen in a study using human cell lines (Yang et al. 2011). In our current rat–BPA study, two of these genes, Colec12 and Slc8a1, were induced by only two doses of BPA and exhibited a typical nonmonotonic dose–response pattern as supported by the statistical analyses.

BPA is a weak estrogen and known to interact with both ER and (ESR1 and ESR2) and with G protein-coupled estrogen-receptor 1 for downstream signaling (Bergeron et al. 1999; Kurosawa et al. 2002; Li et al. 2013; Thomas and Dong 2006). Transcription factor binding-site prediction of the 78 genes, in which they shared estradiol as the common regulator, did not identify any classical ER-binding elements or known tethering binding sites. However, there was significant enrichment for AP-2–binding sites (including both alpha and gamma subunits) in the 78 genes. Both subunits were able to transactivate a cloned ESR1 promoter in a breast cancer cell line (Tan et al. 2011), but served as a negative regulator of ESR1 in endometrial cancer cell lines (Lin et al. 2016). How BPA regulates the expression of those genes via AP-2 signaling in endometrial cancer requires further investigation.

In this study, we further investigated whether those specific BPA-responsive genes were associated with endometrial cancer progression. A total of 476 human orthologous genes robustly predicted the overall survival of endometrial cancer patients, and this was independent of age and race. Nevertheless, the two groups of patients dichotomized by the expression level of the genes showed a significant difference in terms of age, race, menopausal status, obesity status, cancer histology/grade, and ESR1 expression level. An interesting aspect is that the same characteristics were recently reported for type 2 endometrial cancer patients (i.e., being older, less obese, nonwhite, menopausal, and with an advanced-stage disease) (Feinberg et al. 2019). More importantly, the patients identified in this study to have a shorter time to recurrence-free survival had lower expression of (fewer transcripts) in their cancer tissue, which is a key characteristic of type 2 endometrial cancer (i.e., estrogen-independent)(Jongen et al. 2009; Mendivil et al. 2009; Shen et al. 2017). Serous endometrial cancer is estrogen-independent and was reported to express a low level of ESR1 (Kounelis et al. 2000). Therefore, this could explain why this particular BPA-responsive gene set was associated with in our gene ontology analysis. It is possible that these genes could be used to accurately differentiate between type 1 and 2 endometrial cancer.

The core study findings published by the CLARITY-BPA consortium included more detailed vaginal cytology and pathology data (Camacho et al. 2019). Based on their vaginal cytology data, the BPA effect was not found. A similar conclusion was made in our study, but the time points investigated and the evaluation protocol used were different. In our transcriptome study, we identified a gene signature that was significantly associated with the overall survival of endometrial cancer patients. Based on the pathological data published in the core study (Camacho, et al. 2019), the incidence rate of endometrial hyperplasia (noncystic) appeared to be higher in the interim sacrifice experiment (1 y) for all “continuous-dose BPA” groups compared with the CTL group, despite not being statistically significant. However, this was not observed in the terminal sacrifice (2 y) experiment (Camacho et al. 2019). This may imply that the effect of BPA on the uterus could be age-dependent—having less influence in aged animals—possibly due to the diminished effect of the hormone in aging animals. Because we did not observe any gross endometrial tumors in 1-y-old rats in the current study of in a study of 2-y-old rats by the NCTR, the gene set we discovered is unlikely to be tumor-initiating, but it may promote cancer progression. This possibility could be tested in the future with a 2-hit animal model using a known cancer inducer followed by BPA treatment.

In summary, we identified two low, estrogenic doses (25 and of BW/d) of BPA that differentially affect gene expression in the uterus in a rat model of lifelong BPA exposure. This gene signature could potentially be used to explain the molecular subtype of endometrial cancer in humans and to predict survival. Because type 1 and 2 designations have never been part of the formal staging or risk stratification in clinics, understanding their molecular signatures would provide an opportunity for developing early diagnostic and prognostic biomarkers for endometrial cancer.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.