Dissection of Estrogen Receptor Alpha Signaling Pathways in Osteoblasts Using RNA-Sequencing

The effects of 17-β-estradiol in osteoblasts are primarily mediated by the nuclear transcription factors, estrogen receptor (ER)α and ERβ. ERs function through three general modes of action: DNA-binding dependent through estrogen response elements (EREs; designated nuclear ERE signaling); nuclear signaling via protein-protein interactions to other transcription factors (nuclear non-ERE signaling); and extra-nuclear signaling (membrane-bound functions of ERs). Identification of the specific transcriptional signatures regulated by each of these modes of action should contribute to an enhanced understanding of estrogen signaling in osteoblasts. To achieve this goal, we utilized specific mutations of ERα that eliminate the ability of the receptor to signal through a specific mode of action. The non-classical ERα knock-in (NERKI) mutation is incapable of signaling through direct DNA binding to EREs and the nuclear only ERα (NOER) mutation eliminates all membrane-localized signaling. Comparison of the gene expression patterns elicited by these mutations with the wild-type ERα (WT) pattern provides mode-specific data concerning transcriptional regulation by ERα. We expressed these constructs in the ER-negative osteoblastic cell line hFOB (−/+ estrogen) and performed global RNA-sequencing. Using a series of pair-wise comparisons, we generated three lists of genes that were regulated either by the nuclear ERE-dependent, nuclear ERE-independent, or extra-nuclear actions of ERα. Pathway and gene ontology analyses revealed that genes regulated through the nuclear ERE and nuclear non-ERE pathways were largely involved in transcriptional regulation, whereas genes regulated through extra-nuclear mechanisms are involved in cytoplasmic signaling transduction pathways. We also intersected our data with genes linked to bone density and fractures from a recent genome-wide association study and found 25 of 72 genes (35%) regulated by estrogen. These data provide a comprehensive list of genes and pathways targeted by these specific modes of ERα action and suggest that “mode-specific” ligands could be developed to modulate specific ERα functionality in bone.


Introduction
Estrogen (E) signaling is mediated by two receptors, estrogen receptor (ER)a and ERb. Although important for many tissue systems, ER signaling is of particular interest with regards to bone biology, as declining E levels during the menopause lead to increased bone turnover and bone loss, as is observed in postmenopausal osteoporosis [1]. Although extensive research has been conducted to understand the role of ERs in bone through use of cellular and animal models, the understanding of gene regulation and cellular pathways regulated by E in osteoblasts is still incomplete.
Mechanistically, ERs modulate gene expression using a number of distinct modes of action. The classical pathway involves direct binding to estrogen response elements (EREs) in the control regions of genes, whereas an alternative nuclear pathway involves protein-protein interactions which are ERE-independent. To facilitate study of these pathways, an ERa mutation [non-classical ERa knock-in (NERKI)] has been made that eliminates ERE but retains non-ERE signaling (both nuclear and extra-nuclear) [2], allowing investigators to study the effects of non-ERE mediated ERa signaling alone. Mice harboring this mutation display an osteoporotic phenotype, demonstrating that nuclear ERE-mediated signaling is important in regulating bone metabolism [3,4,5].
It has become increasingly appreciated that the nuclear actions of ERs cannot completely explain all aspects of ER signaling. Compelling evidence demonstrates that estrogen can act rapidly (within 5-15 minutes of E treatment) through non-genomic mechanisms, to regulate signal transduction through pathways such as ERK and PI3K [6,7,8,9]. This mechanism involves tethering of ERa to the plasma membrane (PM) via palmitoylation of a cysteine residue in its E domain, which facilitates association with caveolin-1 [7]. A mutation of this cysteine has been made (termed ''nuclear-only ER'' or NOER) that completely eliminates membrane localization of ERa [10], and thereby prevents PM ERa signaling. Therefore, in comparison with wild-type ERa (WT), the NOER receptor can be used as a tool to identify those genes regulated through either the nuclear or membraneassociated ERa pathways.
In an effort to identify target genes and pathways regulated by each of these distinct signaling mechanisms, we expressed WT, NERKI, or NOER receptors in an ER-negative osteoblast cell model and measured global gene expression patterns following estrogen treatment using RNAseq. Comparison of estrogendependent gene expression patterns allowed us to compartmentalize these patterns according to the known cellular mechanisms used by ERa, such as dependence on an ERE or whether these genes are regulated through the extra-nuclear (e.g. membraneassociated) function of ERa. These data identify characteristic Eregulated gene signatures from three distinct ERa regulatory mechanisms, which could potentially be used for targeting specific cellular ERa pathways.

Materials and Methods
Cell culture reagents, adenoviruses and infection of hFOB 1.19 cells The hFOB 1.19 human fetal osteoblastic cell line (hFOB), produced by Harris and colleagues [11], was passaged in phenol red-free aMEM growth medium (Invitrogen, Carlsbad, CA) supplemented with 1X antibiotic/antimycotic (Invitrogen), 10% (v/v) fetal bovine serum (Hyclone, Logan, UT) and 300 ug/mL G418 selection antibiotic (Invitrogen). The Ad-ERa (expressing WT mouse ERa) and Ad-NERKI (expressing mouse ERa containing the double mutation E207A/G208A, which is incapable of binding DNA [2]) adenoviral vectors were produced as previously described [12]. To produce the nuclear-only ERa (NOER), which exhibits solely nuclear signaling [10], the C451A mutation was introduced into WT ERa using the QuikChange II Site-Directed Mutagenesis Kit (Agilent Technologies, Santa Clara, CA) and subsequently used to produce an adenoviral vector (Vector Biolabs, Philadelphia, PA), resulting in ''Ad-NOER''. All ERa molecules were FLAG-epitope tagged at the N-terminus. hFOB cells were plated in 10-cm culture dishes (n = 6) at a final density of 2610 4 cells/cm 2 . Cells were infected at the time of plating at a multiplicity of infection (MOI) of 15 (Ad-ERa), 30 (Ad-NERKI) and 22.5 (Ad-NOER), which is herein shown to result in approximately equivalent protein levels for ERa, NERKI and NOER (Fig. 1). The cells were infected in the presence of 8 mg/ mL hexadimethrine bromide (Sigma-Aldrich, St. Louis, MO) to enhance adenoviral infection. Following another 24 hour incubation to allow for complete infection, the cells were treated with either ethanol vehicle (0.1% v/v) or 10 nM 17-b-estradiol (0.1% v/v; Sigma-Aldrich) in the presence of media containing triplestripped charcoal-treated FBS for an additional 24 hrs.

Western Blotting
Whole cell extracts from adenoviral-infected hFOB cells was prepared and subjected to western blotting analysis as previously described [12]. A total of 10 mg protein, as determined by the Pierce BCA Protein Assay Kit (Thermo Scientific, Rockford, IL) was run from each adenoviral condition. The blots were incubated with primary antibodies directed against mouse-anti-Flag (1:1000) and mouse-anti-b-actin (1:10,000). The mouse anti-IgG secondary antibody linked to horseradish peroxidase was used at 1:5000 dilution. All antibodies were purchased from Sigma-Aldrich. Blots were visualized using enhanced chemiluminscence (GE Lifesciences, Piscataway, NJ) and exposed to X-ray film.

RNA/cDNA isolation
Total RNA was prepared using RNeasy minicolumns (Qiagen, Valencia, CA) and treated with RNase-free DNase (Qiagen) to remove potential contaminating DNA, as previously described [12]. The resulting RNA was used for either RNAseq analysis or in semi-quantitative real-time PCR analysis (QPCR). For the latter analysis, one microgram of total RNA was used in a reverse transcriptase (RT) reaction using the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems by Life Technologies, Foster City, CA). The cDNA was diluted 1:5 with water prior to QPCR analysis.

RNA sequencing (RNAseq)
RNA libraries for RNAseq analysis were prepared from 100 ng of isolated total RNA from each sample using the manufacturer's instructions. Unique indexes were incorporated at the adaptor ligation step for loading multiple samples per flow cell. Three distinct indexed libraries were loaded per flow cell and sequenced on an Illumina HiSeq 2000 using TruSeq SBS sequencing software (version 3) and SCS data collection software (version 1.4.8). Base calling was performed using Illumina RTA (version 1.12.4.2). An average of 130 million reads per sample was achieved resulting in ,96% mapped reads. A more detailed description of the RNAseq methodology used is as previously described [13]. The RNAseq data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus [14] and are accessible through GEO Series accession number GSE55769 (http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc = GSE55769).

Semi-quantitative Real-time PCR Analysis (QPCR)
The PCR reactions were run in the ABI Prism 7900HT Real-Time System (Applied Biosystems, Carlsbad, CA) using the Quantitect SYBR Green Master Mix (Qiagen), as previously described [15]. The method for data normalization using multiple reference genes and threshold calculations are as previously described [15]. Primer sequences for individual genes were designed using the Primer Express program (Applied Biosystems) and are available on request.

Statistical analyses
For the QPCR analyses, calculations and statistical comparisons were performed using Microsoft Office Excel 2003 (Microsoft Corp., Redmond, WA) and the data are presented as the mean 6 SEM. All values of p,0.05 were considered statistically significant using two-sided t tests. The RNAseq data was analyzed essentially as previously described [13]. Briefly, paired-end reads from the raw RNAseq data were aligned using TopHat (version 2.0.6) [16] against the h19 genome build using the bowtie 1 option [17], and quality control assessments were made using RSeQC software [18]. Gene counts were generated using HTSeq software and gene annotation files were obtained from Illumina. Gene expression data exhibiting a p-value#0.05, a false discovery rate of q#0.05 and median gene counts in at least one group of 10, were used for further investigation. In order to identify which genes are regulated via non-ERE or extra-nuclear actions of ERa, the WT ERa data was compared the NERKI and/or NOER RNAseq datasets. Overlapping genes were identified using the ''Venny'' web resource (http://bioinfogp.cnb.csic.es/tools/venny/). For the genes regulated via each of the three modes of action, pathway analysis was performed using the Ingenuity Pathway Analysis software (Ingenuity Systems, Redwood City, CA). Gene ontology (GO) terms from each dataset were generated using DAVID Bioinformatics Resources Version 6.7 [19].

RNAseq analysis of hFOB human osteoblastic cells expressing ERa modal variants
The goal of this study was to identify and characterize estrogendependent gene expression patterns elicited by ERa through the nuclear ERE-dependent, nuclear ERE-independent, or through extra-nuclear signaling (e.g. membrane) in the ER-negative, human osteoblastic cell line hFOB [11] using global RNA sequencing (RNAseq). As tools to facilitate this approach, we utilized an ERa mutation that eliminates DNA binding through EREs (NERKI) [2,20] and an ERa mutation which can only signal through the nucleus (NOER) [10]. By comparing the gene expression patterns of these mutant ERa receptor with wild-type ERa (WT), we can identify those patterns elicited by each mode of ERa action. Importantly, our analysis required that estrogenregulated genes in each of the three modes of action were also regulated in WT, thus ensuring that these are physiologically regulated genes by ERa. Therefore, hFOB cells were infected with adenoviruses expressing FLAG-tagged WT, NERKI or NOER. Western blot analysis confirmed that all ERs were expressed at ,equal levels (Figs. 1A/B).
To identify differentially expressed genes with each of the ERa modal variants, hFOB cells were infected with adenoviral vectors expressing WT, NERKI or NOER receptors, treated with estrogen (10 nM E2) for 24 h and RNAseq was performed (see Materials and Methods for details). In a previous study, we had determined that genes exhibiting a median gene count of at least 10 in one group represent an ''expressed'' gene [13]. Therefore genes with a median gene count of less than 10 in both comparison groups were called ''non-expressed'' and excluded from further analysis. In this study, the remaining genes exhibiting a p-value# 0.05 and a false discovery rate of q#0.05 between control and estrogen treatments were used for further analyses.
Using these criteria, an estrogen-regulated gene list for WT was generated. In the WT treatment groups where all modes of estrogen action are preserved, a total of 4353 estrogen-regulated genes were identified (Table S1). Thirty randomly chosen genes from this dataset were analyzed by QPCR to validate the RNAseq data and a high correlation (r = 0.98, P,0.001) was observed (Figure S1), demonstrating a high level of accuracy of the RNAseq technology.
The flow-chart depicted in Figure 2 describes our data analysis strategy in detail. Similarly to WT, the NERKI and NOER datasets were passed through our data filters (median gene count$ 10, p#0.05 and q#0.05) and 2375 estrogen-regulated genes were identified from the NERKI dataset (termed ''ERE-independent'') and 3147 estrogen-regulated genes were identified from the NOER dataset (termed ''Nuclear-only) ( Fig. 2A). Complete lists for all three datasets can be found in Table S1 (tabs 1-3, respectively). To identify the ''ERE-dependent'' genes, the 2375 genes from the NERKI dataset were subtracted from the 4353 genes from the WT dataset ( Fig. 2B; top row). This subtraction results in 1978 genes that are regulated by WT but not regulated by NERKI, therefore by inference these genes are regulated via EREdependent mechanisms. Similarly, to identify those genes whose estrogen-dependent regulation occurs from signals generated outside the nucleus (termed ''Extra-nuclear), the 3147 genes in the Nuclear-only dataset were subtracted from the 4353 genes from the WT dataset (Fig. 2B, bottom row). This subtraction results in 1206 genes that are regulated by WT, but not the NOER, suggesting these genes are regulated in part by mechanisms originating outside of the nucleus. However, our ERE-dependent and ERE-independent gene lists include regulation occurring in both the nucleus and outside of the nucleus. To identify nuclear genes that are ERE-dependent and EREindependent, each of these lists were overlaid with the 1206 genes in the Extra-nuclear list. For the ERE-dependent genes, this resulted in 1167 genes termed ''Nuclear ERE-dependent'' and with 811 genes from the Extra-nuclear list ( Fig. 2C; left side). Similarly, for the ERE-independent genes, this resulted in 1980 genes termed ''Nuclear ERE-independent'' with 395 genes from the Extra-nuclear list ( Fig. 2C; right side). Note that the sum of the Extra-nuclear genes from both comparisons is the complete Extranuclear gene list (811+395 = 1206). In this way, the lists comprising the three modes of ERa action (Nuclear ERE-dependent, Nuclear ERE-independent, Extra-nuclear) were identified. Complete lists for these datasets can be found in Table S2 (tabs 1-3, respectively). The percentage breakdown and number of estrogen regulated genes from each dataset representing these distinct modes of ERa action is shown in Figure 3. The top 20 up-and down-regulated genes for each dataset are shown in Tables 1-3.

Ingenuity pathway analysis (IPA) and gene ontology (GO) analysis of the datasets from the three modes of ERa action
Following global identification of estrogen-regulated genes comprising the three modes of ERa action (described above), we performed IPA on each dataset to identify those cellular pathways mostly like regulated in each of these modes (see Table S3 for complete lists). We focused on those pathways with known regulatory roles in osteoblast biology. Pathways with known roles in bone biology for the Nuclear ERE-dependent (Table 4), Nuclear ERE-independent (Table 5) and Extra-nuclear (membrane) ( Table 6) are indicated. We also performed gene ontology (GO) analysis to categorize the target genes regulated in each mode of ERa action into specific molecular functions (GOTERM_MF_FAT) or cellular components (GOTERM_CC_FAT) (see Table S4). We found that the target genes regulated in the nuclear ERE and non-ERE mode of ERa action were largely enriched for genes involved in transcriptional regulation. Interestingly, the target genes regulated by the extra-nuclear mode of ERa action, presumably occurring through membrane-associated processes, were largely enriched for genes involved in metabolic pathways (i.e. NADH/NADPH oxidoreductase activity and FAD binding) and structural components of the ribosome.

Comparison of the datasets representing the three modes of ERa action with genes with known involvement in bone biology
To place our results in a more physiologically relevant context in terms of bone biology, we intersected our gene lists corresponding to the three modes of ERa action with a recent meta-analysis that identified 72 genomic loci with genome-wide significance to bone mineral density and/or fracture rates [21]. As seen in Figure 4, a total of 25 (,35%) of the GWAS-identified genes are found in our datasets.

Discussion
Estrogen receptor (ER)-a can utilize a number of different mechanisms to modulate gene expression, including nuclear EREdependent, nuclear ERE-independent, and extra-nuclear signaling (e.g. membrane). The goal of this study was to identify the global gene expression patterns regulated by each of these mechanisms in the ER-negative, human osteoblastic cell line hFOB [11] using global RNA sequencing (RNAseq), in an effort to further understand how each mode of ERa action contributes to the overall response of osteoblasts to estrogen. Using specific ERa mutants which eliminate one mode of ERa action, and by comparing these datasets to the WT dataset, we were successful in identifying genes regulated by either the nuclear ERE-dependent, nuclear ERE-independent, and extra-nuclear actions of ERa.
Examination of the relative proportions of the modes of ERa action revealed that a majority of estrogen-regulated genes (72%) are regulated at least partially via nuclear mechanisms. This data is supported by recent report that 82% of estrogen-regulated genes in the uterus are regulated via nuclear ERa pathways [22]. Of these, nearly two-thirds (63%) are from ERE-independent mechanistic pathways and one-third (37%) from ERE-dependent mechanistic pathways. This suggests that a common mechanism for ERa to influence estrogen-dependent transcription in the nucleus is through cooperation with other DNA-bound transcription factor complexes. Our pathway analysis of each mode of ERa action supports this contention, since 135 cellular pathways were significantly regulated by the nuclear ERE-independent dataset, whereas only 48 and 33 pathways were significantly regulated in the nuclear ERE-dependent and extra-nuclear datasets, respectively.
Recent research has demonstrated that ERa is found at the plasma membrane of nearly all ERa-expressing cells [10], and therefore contributes to overall estrogen signaling in most ERpositive cells, including osteoblasts. Indeed, we find that 28% of all estrogen-responsive genes originate from extra-nuclear (e.g. membrane) compartment, since they are regulated by WT ERa but not the nuclear-only ERa (e.g. NOER). In addition, accumulating evidence now suggests that significant cross-talk occurs between the membrane and nuclear ERa signaling [10,23,24,25]. A well-known example of this cross-talk is in the estrogen-dependent regulation of the Ccnd1 gene in breast cancer cells, where membrane-initiated activation of PI3K/AKT and ERK signaling, as well as ERE-dependent mechanisms are necessary for maximal estrogen-stimulated transcription [7,26,27].
An important and interesting aspect of this research is that estrogen-regulated pathways that are either in common or unique to the three modes of ERa action can be identified. Comparison of the pathways from all modes of ERa action revealed a few shared pathways, such as: integrin and integrin linked kinase (ILK) signaling, tight junction signaling, IL8 signaling, MAPK mTOR signaling. These may represent a more generalized function of ERa action, in which multiple modes are involved. However, most of the significant pathways are unique to each mode, suggesting a unique biology is involved dependent on how ERa is functioning. Unique pathways regulated in the nuclear ERE-dependent mode include EGF signaling, sonic hedgehog signaling and many cholesterol and nucleotide biosynthetic pathways. This is in contrast to the more common signaling pathways regulated in the nuclear ERE-independent mode of ERa action such as PPAR, BMP, Wnt, GR, IL6, TGFb, among many others. Unique signaling pathways in the extra-nuclear mode of ERa action include FAK, AMPK, calcium and eNOS signaling. Understanding how each mode regulates a unique aspect of overall ERa signaling in osteoblasts may lead to the generation of ''modespecific'' ligands used to target specific pathways in estrogenresponsive cells. At a physiological level, estrogen signaling is well recognized as an important determinant of bone mineral density (BMD), as declining estrogen levels during the menopause lead to increased bone turnover and bone loss, as is observed in postmenopausal osteoporosis [1]. Therefore it is not surprising that significant overlap exists upon intersection of our estrogen-dependent gene lists with the 72 loci identified in GWAS studies as being related to BMD and/or fracture rates [21]. These include genes with wellknown roles in osteoblast biology, such as Runx2, Spp1 (osteopontin), Esr1 (ERa) and Tnfrsf11b (Opg), as well as other genes with less recognized roles in bone. This demonstrates that all modes of ERa action contribute to gene regulatory patterns important in bone biology. Future studies are needed to examine the role of these genes in bone diseases such as osteoporosis.
In reality, these modes of ERa action do not exist in isolation. It is well known that the rapid effects of membrane estrogen signaling have modulatory effects on nuclear estrogen signaling through phosphorylation of ERa and/or coregulatory transcription factors, as well as induction of numerous signaling cascades  (reviewed in [6,28]). In this way, elimination or modulation of any one mode of ERa action affects the whole system, eventually affecting gene expression and phenotypic expression of the cell, and ultimately the organ and organism [9]. These actions may be very different in various estrogen-responsive tissues, depending on the amount of ERa present and the complement of transcription factors expressed in that particular cell. Since we treated the ER-infected hFOB cells for 24 h, our data most certainly contains not only genes that directly respond to ligand-bound ERa (a primary response), but also those genes whose regulation is dependent on primary response genes (e.g. secondary response genes). Although this is an important aspect to consider when examining the gene lists produced in this study, our primary premise is that these 3 modes of ERa action which are mechanistically distinct, ultimately put into motion a modespecific program that leads to the regulation of distinct subsets of genes. Future studies with shorter treatment periods will be necessary to limit these lists to primary response genes.
We recognize that our study has several limitations. First, we are using an in vitro cell culture system and since the hFOB cells are ER-negative, we expressed these ER mutations using an adenoviral vector. Using an ER-negative cell model was essential to  BMP4, NFATC3, PIK3R1, TCF7, BMP8A, TGFB1, PIK3CG, WNT7B, WNT4, SMAD1, ADAMTS5, BMP8B,  TNFRSF11B, WNT9A, PPP3CC, MAPK12, PIK3R3, IL18, SFRP1, RELA, NFKBIE, FZD1, NFKB1, NFATC1, SMURF1,  JUN, NFKBIA, BMPR1A, NGFR, TNFRSF1B, TNFRSF1A, DVL1, DLX5, IKBKE, NFATC4, TNFRSF11A, CALM1, FZD8,  FOS, TRAF2, FZD4, FOXO1, CSF1, IL1B  eliminate any endogenous ER expression that could confound the results. In addition, during development of this cell culture system we carefully titrated the various adenoviral vectors and used the minimal amount to achieve equal expression of each receptor in . 95% of the cells. Further, while we have not formally compared the levels of ERa expression in our transfected cells with a panel of WT osteoblasts, it is likely that these levels are considerably higher in our transfected cells, despite the fact that we tried to use fairly low viral titers for the transfections. This issue is compounded by the fact that levels of ERa expression in vitro vary considerably with the stage of osteoblast differentiation. Thus, while the strength of our in vitro system is that we were able to study the various estrogen signaling pathways in isolation, we recognize that further validation of these findings in in vivo systems with physiological levels of ERa are needed. A second limitation is that by using hFOB cells, we limited our study to osteoblasts. However, these datasets may be useful in identifying similar ERa action in other estrogen-responsive tissues, such as breast and uterus. Further studies are needed to examine the genes/pathways regulated by these modes of ERa action in other estrogen-responsive tissues. A third limitation is that only ERa, and not ERb was studied. These modes of action also exist for ERb, including the membraneassociated ERb mechanistic pathway that has been extensively studied in cardiomyocytes [29,30]. Future studies will be needed to clarify these modes of ERb action and interaction with ERa, which will be considerably more complex.
In conclusion, we have identified three molecular modes of ERa action in human osteoblastic cells using RNAseq methodology. The identification of estrogen-dependent gene expression patterns and pathways within these modes of action is an important first step into understanding the complete picture of how ERa action contributes to overall ERa signaling in bone, and may provide the data necessary for generation of ''mode-specific'' ligands to preferentially influence a single mode of ERa action. Figure S1 QPCR confirmation of the RNAseq dataset from wild-type ERa. Thirty randomly chosen genes from the wild-type ERa dataset were chosen and QPCR was performed on independent samples. The plot shows a high degree of concordance between the RNAseq dataset and the QPCR analysis, with an r = 0.98 and P,0.001.

(TIF)
Table S1 Complete gene lists of the estrogen-regulated genes from the 3 ERa variants in this study. The estrogen-regulated genes from the wildtype ERa, ERE-dependent, and Nuclear-only ERa variants (from Fig. 2A) are listed in each respective tab. The gene name, accession numbers, fold-change with estrogen (FC), pvalue and false discovery rate (FDR) are indicated. (XLSX)

Table S2
Complete gene lists of the estrogen-regulated genes comprising the 3 modes of ERa action. The ''Nuclear EREdependent'', ''Nuclear ERE-independent'' and ''Extra-nuclear'' estrogen regulated genes are indicated (from Fig. 2C) with the same values as described in Table S1. (XLSX)

Table S3
Complete Ingenuity pathway analysis (IPA) of the genes from the 3 modes of ERa action. The IPA pathways for the genes for each mode of ERa action are indicated. (XLSX)

Table S4
Gene ontology (GO) analysis of the genes from the 3 modes of ERa action. Gene ontology (GO) analysis was performed to categorize the target genes regulated in each mode of ERa action into specific molecular functions (GOTERM_MF_FAT) or cellular components (GOTERM_CC_FAT). (XLSX) Figure 4. Comparison of the RNAseq datasets to gene with known involvement in bone biology from genome-wide association studies (GWAS). The nuclear ERE-dependent, nuclear ERE-independent, and extra-nuclear datasets were intersected with the 72 genes identified by GWAS studies. The overlapped genes for each mode of ERa action are listed and color-coded for clarity. The directionality of the estrogen-dependent regulation is denoted with an arrow following the gene symbol. doi:10.1371/journal.pone.0095987.g004