Novel estrogen-responsive genes (ERGs) for the evaluation of estrogenic activity.

Estrogenic chemicals are found abundantly in natural and industrial materials, and they exhibit a variety of cellular and physiological activities [ 1 ]. Along with the increasing interest in the estrogenic activity of these chemicals and their beneficial applications, there is growing need to develop new technologies for their detection and characterization. Thus, a variety of assays have been developed, such as ligand-binding assays, reporter-gene assays, yeast two-hybrid assays, transcription assays, protein assays, cell assays, and animal tests, which are based on the mechanisms of estrogen action at the levels of molecules, cells, tissues, and the whole body [ 2 ]. Among the assays, gene expression profiling, such as the DNA microarray assay, is a technology used to detect the alterations of gene expression by monitoring the amount of mRNA or proteins of the estrogen-responsive genes (ERGs), which can detect the estrogenic activity of chemicals and mixtures of chemicals [ 3 – 5 ]. Searching for useful genes was challenging when the technological development was in its infancy. Recently, the development of DNA microarray technology for comprehensive searches of genes that were newly identified by the human genome project accelerated the search for ERGs, where tens to sometimes over a thousand ERGs were identified from more than ten thousand human genes and expressed sequence tags [ 6 – 11 ]. The identified ERGs have been examined for use in various fields such as food, clinical, pharmacological, and environmental applications [ 1 ]. For example, better understanding of the gene-regulation networks based on the ERGs is critical for the development of therapeutics for breast cancer [ 5 ]. Likewise, gene expression profiling based on ERGs has been used for environmental studies [ 12 ]. While the number of ERGs used for gene-expression profiling is important to reliably evaluate the estrogenic activity, the differences in the quality of ERGs due to various factors, such as their amounts in cells, the degrees of their responses, and their different cell functions, can cause unpredictable errors in the assay results. Other factors also need to be considered such as the cost, laboriousness, and the requirements of devices and equipment. The proportion of genes quantitatively showing differences in expression under different stimulations is roughly a few to a few tens of % of all human genes [ 13 ]. When considering statistical stability, the number of the genes that can be used to predict specific results, like clinical outcomes, should be no more than 1,000 and preferably less than 100 [ 13 ]. For example, MammaPrint, a US Food and Drug Administration-cleared molecular diagnostic test used for predicting the risk of breast cancer recurrence, contained 70 genes [ 14 ]. On the other hand, gene expression profiling has been used to predict estrogen activity, where a number of chemicals were examined by DNA microarray-based assays (reviewed by Kiyama & Zhu, 2014 [ 12 ]). A total of 120 ERGs, consisting of six functional groups of genes, were used to evaluate estrogenic activity, with comparable reliability to other estrogen assays [ 15 ]. However, methods for reliably identifying ERGs with statistical stability within satisfactory limits have not been reported, even though more information about genes and their functions has been obtained. Here, we determined the ERGs that have statistical stability based on RNA sequencing (RNA-seq) analysis and examined their usefulness by focusing on their reliability and applicability to predict the estrogen action of chemicals. The recent progress in RNA-seq has enabled its wider use in applications such as food quality control, environmental materials, and medicine [ 16 – 20 ]. For estrogen actions, RNA-seq data obtained from estrogen-stimulated MCF-7 cells were compared with reverse transcription (RT)-PCR and DNA microarray data, suggesting mutually consistent or respectively complimentary gene expression profiles among them [ 21 ]. Furthermore, functional annotations of RNA-seq data with databases, such as the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, have been used to identify targets of antitumor agents against breast cancer such as diallyl trisulfide (a garlic metabolite) [ 22 ], 6-thioguanine [ 23 ], shikonin [ 24 ], tamoxifen [ 25 ], and chickpea isoflavones [ 26 ]. Based on RNA-seq analysis, novel mechanisms of cancer progression and metastasis involving non-coding RNAs [ 27 , 28 ], enhancer RNAs [ 29 ], and circular RNAs [ 30 , 31 ] were demonstrated. Functional regulatory networks examined by comparing chromatin immunoprecipitation-sequencing (ChIP-seq) data and time-course RNA-seq data revealed that transcription factor MYC plays an important role in connecting the networks between promoters and enhancers [ 32 ]. Here, we obtained a set of novel ERGs by RNA-seq that can be used for gene expression profiling of potential estrogenic chemicals and demonstrated their usefulness by examining their statistical reliability and their functional relevance.

Gene expression profiling is a method to evaluate the estrogenic activity of chemicals, and a variety of assays have been developed based on transcriptomic and proteomic technologies [ 1 ]. Here, we screened ERGs for evaluating estrogenic activity by transcriptomic assays such as RNA-seq and RT-PCR. To obtain a set of ERGs, we adopted a two-step process involving a correlation analysis and a coefficient of variation (CV) analysis ( Fig 1A ). MCF-7 cells were used as a cell system for the evaluation, which was performed first by two assays, a cell proliferation assay (SRB assay; Fig 1B ) and protein assay (Western blotting to detect the phosphorylation of two marker proteins, Erk1/2 and Akt; Fig 1C ), which have frequently been used to evaluate estrogenic activity [ 42 – 44 ]. The cells showed the response to estrogen (17β-estradiol; E 2 ), and the response was inhibited by treatment with an ER antagonist ICI in both assays in the same way as previously reported [ 15 ], suggesting that the cells responded to E 2 through the ER. Then, we performed the RNA-seq analysis to evaluate the cell response at the transcription level ( Fig 1D ). A total of six sets of gene expression profiles for MCF-7 cells treated with E 2 were obtained by RNA-seq, where a set of the 120 ERGs, which were previously used in DNA microarray assays [ 15 ], were used to compare the profiles (first five panels in Fig 1D ). The six combinations of profiles showed R -values of 0.83 to 0.95, indicating strong and statistically significant correlations among the profiles (see the data for all combinations in S2 Table ). In contrast, there was no correlation between the profiles for the E 2 treatment and E 2 +ICI treatment (the last panel in Fig 1D ), suggesting the response to be mediated by the ER. Therefore, the ERGs that we obtained using gene expression profiling methods and by RNA-seq techniques can be used as markers for the evaluation of estrogenic chemicals. (A) A strategy adopted here. (B) Cell proliferation assay. MCF-7 cells treated with E 2 (10 nM), E 2 (10 nM) + ICI (1 μM), or vehicle (0.1% DMSO) were subjected to cell proliferation assays with SRB. The relative proliferation index for MCF-7 cells treated with E 2 (10 nM), E 2 (10 nM) + ICI (1 μM), or ICI (1 μM) alone are shown: * p < 0.05, vs. control (0.1% DMSO). (C) Western blotting. MCF-7 cells were treated with E 2 (10 nM), E 2 (10 nM) + ICI (1 μM), or vehicle (0.1% DMSO) for indicated times and cell extracts were subjected to Western blotting for the evaluation of phosphorylated proteins (P-Erk1/2 and P-Akt) (upper images) and total proteins (T-Erk1/2 and P-Akt) (lower images). (D) Correlation analysis of RNA-seq data. The gene expression profiles for MCF-7 cells treated with E 2 obtained by RNA-seq were compared using a set of the 120 ERGs, and the results are visualized in scatter-plot graphs. The vertical and horizontal axes indicate log 2 -values of FPKMs. R - and p -values were calculated for each graph based on linear regression between two profiles. Next, we examined the quality of the obtained ERGs by CV analysis ( Fig 2 ). Among the six sets of ERGs obtained by RNA-seq ( Fig 1D ), the sets of 203 ERGs that were used to evaluate estrogenic activity by DNA microarray assays [ 15 ] were again examined by CV analysis (indicated by the bar graph in Fig 2 ), where the 203 genes included expression standards and 174 original ERGs, and the sets of 150, 120, 90, 60, or 30 genes were selected from the top of the list of genes ranked by CV values (mean of SD/Av). The top-ranked ERGs exhibited low CV values, 0.17 +/- 0.04 for the top 30 ERGs, indicating that gene expression profiling by RNA-seq can provide mean CV values that are lower than those by DNA microarray assays (mean CV value of 0.22; [ 15 ]). Furthermore, the correlation coefficients ( R -values) for the top 30 ERGs in the five data sets obtained by RNA-seq were 0.97 to 0.99, which are also more reliable than those with nine data sets obtained by DNA microarray assays (0.93 to 0.97; [ 15 ]), confirming that the statistical stability of data is greater in RNA-seq than in DNA microarray assays. The RNA-seq for E 2 was repeated 6 times and gene expression profiles were obtained (E 2 -1 to E 2 -6). The correlation coefficients ( R -values) between one profile (E 2 -1) and the others (E 2 -2 to E 2 -6), giving five sets of data, were calculated for respective numbers of gene sets. The R -values of the total 203 genes, or the indicated numbers (174, 150, 120, 90, 60, or 30) of genes taken from the top of the list in the order of the genes ranked by the values of SD/Av indicating the statistical stability, are shown for the five sets of data in a line graph. The mean of SD/Av for respective numbers of genes are shown in a bar graph; SD/Av for 203 genes is 2.31 ± 6.24 and that for 174 genes is 2.34 ± 6.6. We selected new ERGs for gene expression profiling based on the approach shown in the previous section. A total of 300 ERGs were selected from the ranks created by CV analysis, from all human genes available for RNA-seq (a total of approx. 26,000 genes), which consist of 150 up-regulated and 150 down-regulated genes (see S3 Table for the complete list). Among the 300 ERGs, the top 30 ERGs ( Table 1 ) were examined by real-time RT-PCR ( Fig 3A ), and the data were compared with those obtained by RNA-seq ( Fig 3B ). Both datasets were identical, except for minor differences due to slight variation in the locations of the transcripts that were detected by these assays. Notably, the statistical significance levels are lower for some genes, probably due to the lack of expression levels in the statistical consideration. (A) Real-time RT-PCR analysis for 30 ERGs. The 30 ERGs (listed in Table 1 ) that were selected by RNA-seq analysis in Fig 2 were analyzed again by real-time RT-PCR. The expression level of each gene was normalized by that of β-actin and indicated in a bar graph, where the bars show the mean ± SD (n = 3) of log 2 -transformed ratios of E 2 + data and E 2 - data (E 2 +/E 2 -). (B) RNA-seq analysis for 30 ERGs. The values of FPKMs (n = 6) obtained by RNA-seq were log 2 -transformed and indicated in a bar graph, where the bars show the mean ± SD (n = 6) of the ratios, as shown in panel A. The top 30 genes selected according to SD/Av and expression ratios have stable gene expression. * p < 0.05: between E 2 and the control (0.1% DMSO). The top 30 ERGs among the 300 genes stably responding to E 2 are listed. The newly identified ERGs were examined by signaling pathway analysis using the GO and KEGG databases ( Fig 4 ). The characteristics of ERGs were analyzed first by the statistical evaluation of the relationship of the significance ( p -values; Fig 4A ) or stability (CV values; Fig 4B ) of expressional alterations with the degrees of the alterations (log 2 -transformed fold changes) and then by the functional searches in the GO ( Fig 4C and 4D ) and KEGG ( Fig 4E and 4F ) databases. We confirmed that 30 ERGs belong to the groups of genes showing statistically significant alterations of their expression due to the treatment with estrogen (red or green circles in Fig 4A and 4B ). From the database searches, significant levels of association were found between the genes of significant alterations of their expression and specific cell functions and pathways such as the cellular macromolecule catabolic process (a 262-gene group showing up-regulation), protein localization to organelles (a 220-gene group showing up-regulation), mRNA metabolic process (a 201-gene group showing up-regulation), negative regulation of gene expression (a 290-gene group showing down-regulation), negative regulation of biosynthetic process (a 283-gene group showing down-regulation), negative regulation of cellular biosynthetic process (a 281-gene group showing down-regulation), and negative regulation of macromolecule biosynthetic process (a 276-gene group showing down-regulation) in the GO database ( Fig 4C and 4D ; S4 Table ), and metabolic pathways (a 279-gene group showing up-regulation), protein processing in endoplasmic reticulum (a 58-gene group showing up-regulation), cell cycle (a 55-gene group showing up-regulation), pathways in cancer (a 90-gene group showing down-regulation), human papillomavirus infection (a 66-gene group showing down-regulation), endocytosis (a 46-gene group showing down-regulation), and breast cancer (a 32-gene group showing down-regulation) in the KEGG database ( Fig 4E and 4F ; S4 Table ). Furthermore, significant numbers of the 30 and 300 ERGs that were characterized in Fig 2 were included in the groups of genes found by the database searches ( S4 Table ), suggesting the association of these ERGs with specific cell functions. (A) Volcano plots for ERGs in MCF-7 cells. The genes showing significant ( p < 0.05) up-regulation (red circles) or down-regulation (green circles) are indicated, while those that are not significant are shown by gray circles. The novel 30 ERGs are shown by black dots. The graph shows log 2 -transformed fold-changes (X-axis) and -log 10 of p -values (Y-axis). (B) Volcano plots for the 30 genes that show stable expression in response to E 2 . The graph shows log 2 -transformed fold-changes (X-axis) and mean of SD/Av values (Y-axis). (C, D) GO analysis for the biological processes (BPs). The top 10 significantly affected ( Q < 0.05) BPs for the ERGs that exhibit up-regulation (panel C) or down-regulation (panel D) are shown. The number of genes in each group is visualized by the size of the circles, and the Q -values of the group by color. (E, F) KEGG analysis for metabolic pathways. The top 10 significantly affected ( Q < 0.05) metabolic pathways for the ERGs that exhibit up-regulation (panel C) or down-regulation (panel D) are shown. Since ERGs are likely controlled by ERs through the interaction at their promoter and enhancer regions, we examined the ChIP-seq data that are available in the Gene Expression Omnibus database to determine whether the regions containing these ERGs were identified as ER subtype-specific binding sites ( Table 2 ). The results indicate that 29 ERGs (except LINC02593 ) were classified into four groups; ERα-specific ( SUSD3 , RAPGEFL1 , PKIB , INSYN1 , SPOCD1 , BARX2 , SYNE1 , SRGAP3 , CSTA , and FRY ), ERβ-specific ( ACOX2 and LOXL2 ), ERα/β-specific ( EGR3 , GATA4 , IL20 , PDLIM3 , CTSD , TMPRSS3 , FDFT1 , DOK7 , B4GALT1 , RAB26 , RAP1GAP , MATN2 , CYP1A1 , and CCDC68 ), and other types ( ZNF521 , IGSF1 , and CRISP3 ) ( Fig 5 ). Gene No. is based on the order of the data from real-time RT-PCR ( Fig 3 ). Gene No. 20 ( LINC02593 ) is not listed here because it is a long non-coding RNA and inappropriate for the analysis. a The original data analyzed in this table are available in the NCBI’s Gene Expression Omnibus with accession numbers GSE117569 (for ERα; [ 135 ]) and GSE149979 (for ERβ; [ 151 ]). b The Q -value was determined using the ChIP-Atlas database ( https://chip-atlas.org/ ; [ 152 ]), where each gene was analyzed for the presence of peaks within ± 10 kb from the transcription start site (TSS). c ERα binding is positive (+) when the Q -value of E 2 + is bigger than that of E 2 -. d The integrative genomics viewer (IGV) was used to visualize the results of ChIP-seq, where each gene was analyzed for the presence of peaks within ± 10 kb from the TSS. Then, the S/N ratio was calculated, where the average of E 2 +/ERβ ChIP-seq values, rep 1 and 2, was used as S, and the E 2 +/input value was used as N. e ERβ binding is positive (+) when the S/N ratio is bigger than 1.0.

In this paper, we identified the ERGs that show statistical stability based on the CV values. We evaluated their usefulness as markers for estrogenic activity by examining the stability of the data when compared with correlation analysis, RT-PCR, and functional association with estrogen action through database searches. The ERGs identified here have been known for (1) the enzymes, such as metabolic enzymes, proteases, and protein kinases, (2) the genes with specific cell functions, such as cell-signaling mediators, tumor-suppressors, and roles in breast cancer, (3) the association with transcriptional regulation, and (4) estrogen-responsiveness. The top 30 ERGs were further classified as ERα-binding, ERβ-binding, ERα/β-binding, or ER-independent types based on ChIP-seq data. Therefore, the ERGs identified here represent various cell functions and cell signaling pathways, including estrogen signaling, and thus, may be useful for evaluations of estrogenic activity.

Human breast cancer MCF-7 cells were obtained from the Japanese Collection of Research Bioresources Cell Bank. Antibodies used for Western blotting were those against total Erk1/2 (T-Erk; #9102, Cell Signaling Technology, Ipswich, MA, USA), phospho-Erk1/2 (P-Erk; #9101), total Akt (T-Akt; #4691), and phospho-Akt (P-Akt; #4060). MCF-7 cells were cultured in a phenol red-free RPMI 1640 medium (Gibco, Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco, Thermo Fisher Scientific) and maintained at 37°C with 5% CO 2 in an incubator (Thermo Fisher Scientific). The sulforhodamine B (SRB) assay was performed as described by Dong et al. [ 33 , 34 ]. Before stimulation of MCF-7 cells with chemicals, the cells were cultured at a density of 1.5×10 4 cells per well in RPMI 1640 containing 10% (v/v) dextran-coated charcoal-treated FBS (DCC-FBS) (Gibco, Thermo Fisher Scientific) for 3 days in 24-well plates. After treatment with 10 nM 17β-estradiol (E 2 ; Sigma-Aldrich, St. Louis, MO, USA), 1 μM ICI 182,780 (ICI; Sigma-Aldrich), or 0.1% dimethyl sulfoxide (vehicle) (DMSO; FUJIFILM Wako Pure Chemical, Osaka, Japan) for 3 more days, MCF-7 cells were fixed with trichloroacetic acid (Sigma-Aldrich) at 4°C for 30 min and washed with ultrapure water. After washing, samples were stained with acetic acid containing 0.4% SRB (Sigma-Aldrich) at room temperature for 20 min and then washed with acetic acid. The bound protein was dissolved with 10 mM unbuffered Tris-base (pH = 10.5) at room temperature for 10 min and transferred into 96-well plates to measure the absorbance at 490 nm using the Modular-Designed Multimode Reader SH-9000 (Corona Electric, Ibaraki, Japan). Three independent assays were performed for each chemical, and the data were analyzed by t -test. Western blotting was performed as described by Dong et al. [ 33 , 34 ]. MCF-7 cells at a density of 1.0×10 5 cell per well were cultured for 2 days in DCC-FBS in 6-well plates and then one more day in serum-free medium. The cells were pretreated with ICI for 1 h and then treated with 10 nM E 2 or vehicle (0.1% v/v DMSO). The total protein was extracted from the cells and examined by SDS-PAGE using e-PAGEL (ATTO, Tokyo, Japan) and electro-transferred onto nitrocellulose membranes (Millipore, Billerica, MA, USA) using a semi-dry transfer cell (Bio-Rad Laboratories, Benicia, CA, USA). The membranes were blocked with EzBlock BSA (ATTO) and incubated with the antibodies against signal proteins (T-Erk, P-Erk, T-Akt or P-Akt) overnight at 4°C after appropriate dilution (1:500 or 1:1000). The antibody-antigen complexes were detected with a horseradish peroxidase-coupled goat antibody against rabbit IgG (Cell Signaling Technology) after dilution (1:1000) and visualized using WSE-6100 LuminoGraph I (ATTO). MCF-7 cells at a density of 1.0×10 6 cells per well were cultured in DCC-FBS for 3 days at 37°C in 5% CO 2 . The cells were treated with 10 nM E 2 , 1 μM ICI, or 0.1% DMSO (vehicle) for 2 more days. Total RNA was extracted using a RNeasy Mini-kit (QIAGEN, Venlo, the Netherlands). The DNA libraries for sequencing were constructed using the Truseq Stranded mRNA Library Prep (Illumina, San Diego, CA) and MGIEasy Universal Library Conversion Kit (MGI Tech, Shenzhen, PRC). Then, the libraries were subjected to sequencing using the DNBSEQ-G400 (MGISEQ-2000RS, MGI Tech) instrument, and the obtained reads were aligned using splitBarcode ( https://github.com/MGI-tech-bioinformatics/splitBarcode ). The read data were downsized to 20 M read pairs using seqkit v0.13.0 [ 35 ]. After the adopter sequences were cleaned with cutadapt ver.2.10 [ 36 ], the read data were mapped with HISAT v2.2.1 [ 37 ]. The value of FPKMs (reads per kilobase of transcript per million reads mapped) for each gene was estimated with Cuffdiff v2.2.1 [ 38 ], and the gene expression was analyzed as described previously [ 15 ]. After RNA-seq was repeated 6 times for E 2 , the standard deviation (SD) and the average (Av) for the 203 genes, which includes ERGs and control genes [ 15 ], were determined and used for calculating the log 2 of FPKMs. These genes were evaluated based on the absolute value of SD/Av and used in the correlation analysis with the correlation coefficient ( R value) between one of the assays (E 2 -1) and each of the other five assays (E 2 -2~ E 2 -6). The FPKM values obtained with RNA-seq were log 2 -transformed and a total of 300 genes were selected based on the log 2 value (gene expression ratio). Among the 300 genes, 30 genes were selected by ranking the SD/Av values and further analyzed by real-time RT-PCR [ 8 ]. Real-time RT-PCR was conducted with the CFX Connect real-time PCR detection system using the iTaq Universal SYBR Green One-Step Kit (Bio-Rad). The conditions for real-time RT-PCR were as follows: reverse transcription for cDNA at 42°C for 10 min and denaturation at 95°C for 1 min, followed by 46 cycles of denaturation at 94°C for 10 sec, annealing at 57°C for 30 sec, and extension at 72°C for 20 sec. After PCR, a melting curve was constructed by increasing the temperature from 65 to 95°C to confirm the validity of the reaction. Real-time RT-PCR was repeated 3 times for each gene and the Av and SD were calculated using CFX Maestro Software (Bio-Rad). The nucleotide sequences of PCR primers ( S1 Table ) were obtained from the references for EGR3 [ 8 ], LOXL2 [ 39 ], and SYNE1 [ 40 ], or otherwise by web-based calculation. The FPKM values obtained by RNA-seq were log 2 -transformed and p -values were calculated. The genes showing statistically significant ( p < 0.05) up-regulation or down-regulation, or those without significant changes, were visualized by Volcano plotting [ 26 ], where log 2 -transformed expressional changes (X-axis) and -log 10 -transformed p -values or SD/Av values (Y-axis) are indicated. The pathway analysis was performed with WebGestalt software (WEB-based Gene SeT AnaLysis Toolkit) ( www.webgestalt.org/ ; [ 41 ]) using the gene sets showing statistical significance ( p < 0.05) from the GO ( www.geneontology.org/ ) and KEGG ( www.genome.jp/kegg/ ) databases. False discovery rates (FDRs: Q -values) were estimated using the Benjamini-Hochberg method, and the top ten categories of FDR < 0.05 were visualized.

(XLSX) Click here for additional data file. (DOCX) Click here for additional data file. A total of 300 ERGs containing 150 up- and 150 down-regulated genes among ERGs stably responding to E 2 . (XLSX) Click here for additional data file. (A) up-regulated genes in GO analysis; (B) down-regulated genes in GO analysis; (C) up-regulated genes in KEGG analysis; (D) down-regulated genes in KEGG analysis. (XLSX) Click here for additional data file. (DOCX) Click here for additional data file. (PDF) Click here for additional data file.