The histone demethylase KDM3A regulates the transcriptional program of the androgen receptor in prostate cancer cells

The lysine demethylase 3A (KDM3A, JMJD1A or JHDM2A) controls transcriptional networks in a variety of biological processes such as spermatogenesis, metabolism, stem cell activity, and tumor progression. We matched transcriptomic and ChIP-Seq profiles to decipher a genome-wide regulatory network of epigenetic control by KDM3A in prostate cancer cells. ChIP-Seq experiments monitoring histone 3 lysine 9 (H3K9) methylation marks show global histone demethylation effects of KDM3A. Combined assessment of histone demethylation events and gene expression changes presented major transcriptional activation suggesting that distinct oncogenic regulators may synergize with the epigenetic patterns by KDM3A. Pathway enrichment analysis of cells with KDM3A knockdown prioritized androgen signaling indicating that KDM3A plays a key role in regulating androgen receptor activity. Matched ChIP-Seq and knockdown experiments of KDM3A in combination with ChIP-Seq of the androgen receptor resulted in a gain of H3K9 methylation marks around androgen receptor binding sites of selected transcriptional targets in androgen signaling including positive regulation of KRT19, NKX3-1, KLK3, NDRG1, MAF, CREB3L4, MYC, INPP4B, PTK2B, MAPK1, MAP2K1, IGF1, E2F1, HSP90AA1, HIF1A, and ACSL3. The cancer systems biology analysis of KDM3A-dependent genes identifies an epigenetic and transcriptional network in androgen response, hypoxia, glycolysis, and lipid metabolism. Genome-wide ChIP-Seq data highlights specific gene targets and the ability of epigenetic master regulators to control oncogenic pathways and cancer progression.


INTRODUCTION
Methylation of histone lysine residues is a significant component of epigenetics and is associated with control of gene expression [1]. Specifically, methylation of lysine 9 of histone H3 (H3K9) has been recognized as hallmark of transcriptionally suppressed genes [2]. KDM3A (lysine demethylase 3A; Gene ID: 55818; also referred to as JMJD1A or JHDM2A) is crucial for gene regulation in a variety of biological activities such as spermatogenesis, metabolism, stem cell activity and tumor progression by demethylating mono-or di-methylated H3K9 [3][4][5]. Although the KDM3A protein regulates a wide array of target genes in tissue-and development-specific settings, chromatin modifiers often lack intrinsic DNA sequence specificity. Therefore, how KDM3A is targeted to specific genes is an area of current research interest and important for understanding epigenetic dysregulation in human disease.

Research Paper
functions as a transcriptional coactivator for the androgen receptor (AR; Gene ID: 367) [3,9]. The ability to cooperate with the AR highlights a potential role of KDM3A as coactivator and driving force for sex-specific tissue development as well as for prostate cancer initiation and progression. In PRAD, androgen-dependent signaling plays a key role in the oncogenesis of prostate epithelial cells and the aggressiveness of the malignancy [10,11]. The AR transcription factor belongs to the nuclear receptor superfamily and contains a C-terminal ligandbinding domain. Upon ligand binding, the AR undergoes a conformational change and dissociates from a cytosolic chaperone protein complex. Its ligand-bound conformation allows the AR to dimerize and to translocate into the nucleus [12]. Once in the nucleus, the activated AR dimer binds to androgen response elements (AREs) present in the promoter or enhancer of AR-regulated target genes and recruits coactivators or co-repressors to regulate gene expression [13]. In addition to the AR, KDM3A has been found to regulate expression and/or activity of several transcription factors such as PPARG, KLF2, ESR1, and HOXA1 [14][15][16][17].
In order to further elucidate the impact of KDM3A on the epigenome, we performed chromatin immunoprecipitation in combination with next generation sequencing (ChIP-Seq) of its binding and demethylation activity. We quantified changes of H3K9me1 or H3K9me2 marks, the two substrates of KDM3A, and mapped ARbinding in the CWR22Rv1 prostate cancer cell line in combination with knockdown of KDM3A. Alteration of H3K9 methylation marks mapped to genomic locations coinciding with AR binding pinpoints target genes and oncogenic pathways cooperatively regulated by KDM3A and AR.

Genomic annotation and transcriptional regulation of KDM3A specific demethylase activity
Knockdown of KDM3A in CWR22Rv1 cells showed minor effects on global levels of H3K9me1 or H3K9me2 by Western blot analysis ( Figure 1A), suggesting that KDM3A demethylates a small pool of methylated histone marks and regulates a specific set of gene targets. In order to establish the genome-wide impact of the epigenetic regulator KDM3A, we conducted a matched ChIP-Seq experiment using antibodies specific for histone marks H3K9me1 and H3K9me2 in combination with small hairpin RNA  Table 1). Histone lysine demethylation (KDM) events mediated by KDM3A were defined by gain of methylation ChIP-Seq signals following knockdown of KDM3A ( Figure 1B-1E). Alterations of H3K9me1 and H3K9me2 histone marks upon KDM3A knockdown were evaluated in comparison to reference genomic DNA input or control non-coding shRNA samples. Overall, the peak counts of both H3K9me1 and H3K9me2 ChIP-Seq experiments showed a gain of signal (32244 to 34162 and 23353 to 46599, respectively). Since H3K9 methylation is a mark associated with the highly condensed heterochromatin state, we characterized the specific genomic regions associated with both H3K9 histone methylation marks. KDM events were functionally annotated by mapping bound regions to the human genome and by classifying them according to the nearest gene locus and relative position within coding regions. Promoter or transcription start sites (TSS) and transcription termination sites (TTS) genomic annotations are defined as being within ± 1000 bp of the gene-coding body. Intergenic regions were defined as the remaining regions outside the gene body. In the H3K9me1 ChIP-Seq experiment the intergenic regions were the most frequently found region with 21112 peaks (46.5%) followed by 21495 (47.7%) as intronic regions, 822 (1.8%) as exonic regions, 606 (1.3%) as promoter-TSS regions, 549 (1.2%) as TTS regions, 424 (0.9%) as 3'UTR regions, and 46 (0.1%) as 5' untranscribed (UTR) regions ( Figure  1C). Similarly, the H3K9me2 ChIP-Seq had the intergenic region as its most frequent region with 18195 (55.9%) followed by 13167 (55.9%) as intronic regions, 373 (1.1%) as exonic regions, 355 (1.1%) as promoter-TSS regions, 246 (0.8%) as TTS regions, 204 (0.6%) as 3'UTR regions, and 25 (0.1%) as 5' UTR regions ( Figure 1E). Taken together, ChIP-Seq profiles monitoring histone 3 lysine 9 methylation marks following KDM3A knockdown revealed selective histone demethylation effects of this epigenetic modifier.
Following genomic annotation, we were curious if there was a specific gene expression program underlying demethylation of H3K9me1 and H3K9me2. Half of the detected genes, 8841 (55.1%), contain both H3K9me1/2 marks ( Figure 1F-1G). While the ChIP-Seq data shows that H3K9me1 has more annotated genes compared to H3K9me2 (5121 and 2089 respectively), both histone marks showed an equal fraction of genes being transcriptionally responsive to KDM3A knockdown according to the transcriptomic dataset. Overall, from the transcriptomic experiments, 1408 (58.4%; using a significance cutoff with adjusted p-values below 0.05) genes are reported as differentially downregulated upon shRNA knockdown of KDM3A while 1002 genes are reported as up-regulated by KDM3A knockdown (Figure 1H-1I). The combination of ChIP-Seq histone demethylation events and transcriptomic assessment showed major transcriptional activation by KDM3A, suggesting that KDM3A may synergize with distinct transcriptional regulators for epigenetic control of gene expression.

Identification of an epigenetic and transcriptional network in androgen receptor signaling regulated by KDM3A
Following characterization of histone H3K9me1/2 marks we determined enrichment of transcriptional motifs associated with these histone marks controlled by KDM3A in prostate cancer cells. The goal of this analysis is to identify potential transcription factors that cooperate with KDM3A to regulate gene expression. Using the Jaspar motif database, we conducted an unbiased search for significant enrichment of transcription factor families (analysis of motif enrichment search with p-values below 0.05). Top hits included the androgen receptor, sterol regulatory element binding factor (SREBF), hypoxia inducible factor (HIF), activator protein 1 (AP1) complex of JUN/FOS, Krüppel-like factors (KLF), v-myc avian myelocytomatosis viral oncogene homolog (MYC), and forkhead box (FOX) families of transcription factors with significant enrichments and p-values below 1.0E-04 (Table 1). In addition, we enhanced simple ChIP-Seq-based searches with position-specific matrices to determine which transcription factor motifs were enriched compared to shuffled background sequences. The enrichment analysis showed the androgen response element (ARE) with 2915 incidences as one of the most frequent motifs detected (Table 1). Next, we analyzed significantly altered expression levels upon KDM3A shRNA knockdown (Figure 2A). At the transcriptional level, the androgen response gene set was the most enriched with a p-value and false discovery rate q-value each below 1.0E-20 ( Figure 2B). Next, we sought inferred transcriptional regulators by comparing transcriptional targets to datasets that outline targets of transcription factors through the use of Ingenuity Pathway Analysis. The transcription factors AR, HIF, MYC, and AP1 complex were significantly enriched with p-values below 1.0E-07. Lastly, we merged ChIP-Seq profiles of H3K9Me1/2 and KDM3A transcriptional data focusing on 1408 annotated genes (overlap of H3K9ME1/2 ChIP-Seq with transcriptomic data that were down-regulated upon shRNA KDM3A knockdown) ( Figure 1I, Supplementary Table 2). The data contained the highest enrichment ratio (26.7%) in a significantly enriched set of 27 genes in androgen signaling with p-values below 1.0E-17 and q-values below 1.0E-15 (Supplementary Table 3). In detail, putative KDM3A-regulated genes included pathways involved in androgen response, androgen receptor signaling, androgen biosynthesis, prostate cancer, pathways in cancer, cholesterol homeostasis, bile acid metabolism, aldosterone-regulated reabsorption, and progesterone regulation, hinting at the possibility of hormone nuclear steroid receptor involvement (enrichment with p-values below 2.62E-02 and q-values below 9.35E-02 correcting for multiple hypotheses testing) ( Table 2). Interestingly for the concept of cooperative control, the pathways of Adjusted p-values accounting for multiple hypotheses testing using the false discovery rate controlling procedure of Benjamini and Hochberg of transcriptional motifs below 0.05 in H3K9me1/2-KDM3A ChIP-Seq data were considered significant.  The analysis is based on a query set of 1408 genes positively regulated by KDM3A with coincidence of demethylation and transcriptomic regulation (compare Figure 1). P-values below 0.05 from the hypergeometric distribution for K, the number of elements in the pathway of interest, k, the number of elements in the intersection of the input gene set, N=45956, all known human genes, and n=1408, the number of elements in the input gene set and false discovery rate q-values below 0.10 for multiple hypotheses testing according to the controlling procedure of Benjamini and Hochberg were considered significant.
by all of these different enrichment approaches, suggesting a key role for KDM3A in regulating AR activity.

Matched KDM and AR ChIP-Seq experiments reveal coincidence of demethylase binding, demethylation and AR binding events
Knockdown of KDM3A in CWR22Rv1 cells resulted in loss of KDM3A ChIP-Seq binding accompanied by specific, matched gain of histone lysine 9 demethylation ( Figure 4A-4C). Knockdown of KDM3A had little effect on the protein level of AR [3,18]. We examined the alteration of AR binding by ChIP-Seq with an AR antibody following KDM3A knockdown and quantified the overlap of AR ChIP-Seq events with KDM3A binding and changes in epigenetic H3K9me1/2 marks ( Figure  4A-4D). The activity-based ChIP-Seq array matched with knockdown of KDM3A resulted in 37525 peaks associated with KDM3A binding, 45246 and 32665 H3K9 mono-and di-demethylation (H3K9me1/2-KDM) events, respectively, and 34614 peaks for KDM3A-matched AR binding. Such an experimental design allows one to distinguish between histone demethylase binding ( Figure 4E), epigenetic activity ( Figure 4F), and coactivator binding events ( Figure  4G). Gain of H3K9me1/2 was coupled to specific changes in AR binding in the KDM3A knockdown experiments ( Figure 1A-1E). Overall 37.0% of the AR ChIP-Seq peaks with altered H3K9me1/2 signal were suppressed upon knockdown of KDM3A, while the remaining fraction was not affected. The genome-wide ChIP-Seq analysis is consistent with the biochemical data, demonstrating that KDM3A in effect recruited AR to target genes [3,18]. KDM3A ChIP-Seq and H3K9me1/2-KDM ChIP-Seq in combination with matched knockdown of KDM3A produced an epigenetic network that overlaid with the AR ChIP-Seq data ( Figure 4H-4J). In the case of matched and merged datasets of AR ChIP-Seq in combination with KDM3A knockdown, we identified in total 77911 H3K9me peaks ( Figure 4F) and directly overlaid them with 34614 AR peaks containing 121700 ARE motifs ( Figure 4I). Importantly, using such matched ChIP-Seq analyses, a set of 1912 genes was identified that showed an overlap of demethylation and AR binding events (KDM3A/AR ChIP, 2381 peaks, 1912 genes) ( Figure 4I, Supplementary Table  1). Epigenetic events identified by ChIP-Seq were overlaid with transcriptomic data, defining a set of 421 genes that   Table 4). Similar to the initial gene set based exclusively on H3K9me1/2 marks (57.2%) ( Figure  1), 60.2% of genes scored as activated upon KDM3A/AR coactivation while 39.8% were silenced. Merged ChIP-Seq data of KDM3A/AR coactivation with transcriptomic data of KDM3A knockdown defined a set of 260 genes (whereas activation is interpreted as down-regulation upon shRNA knockdown) ( Figure 4K-4L, Supplementary Table 4).

KDM3A and AR coactivation results in oncogenic pathway activation of AR signaling
H3K9 demethylation is known for stimulating gene expression [3,5,19,20]. We used gene set enrichment analysis to identify signaling networks or functional clusters of genes controlled by both KDM3A and AR. At the gene level we studied pathway enrichment for the following three sets: A) 1912 genes defined by overlapping H3K9 demethylation/AR ChIP ( Figure 4I, Supplementary Table 1-2); B) 421 genes defined by overlapping H3K9 demethylation/AR ChIP with differential expression upon KDM3A/AR knockdown ( Figure 4K, Supplementary Table 2); and C) 260 genes that are down-regulated upon KDM3A or AR knockdown and having overlapping H3K9 demethylation/AR ChIP ( Figure 4K, Supplementary Table  2). The gene sets were designed in a hierarchical fashion such that the parental set A of 1912 genes includes subset B of 421 genes, and B includes subset C of 260 genes (Supplementary Table 2). Such a hierarchical structure of gene sets allows one to monitor and identify requirements and conservation of a functional outcome. The gene set enrichment analysis revealed oncogenic activation of androgen signaling and metabolic pathways with p-values below 1.40E-02 and q-values below 4.44E-02 in hypoxic response, glycolysis, and lipogenic metabolism ( Table 3, Supplementary Table 5). Enrichments of androgen response, metabolic pathways, hypoxia, aldosterone-regulated sodium reabsorption, glycolytic and glycerophospholipid metabolism pathways were conserved and enriched from the overlapping H3K9 demethylation/ AR ChIP (1912 genes) to the genes displaying downregulation upon KDM3A/AR knockdown and overlapping H3K9 demethylation/AR ChIP signal (260 genes). For example, androgen response and metabolic pathways showed incrementally higher enrichment with better defined input gene sets (significant enrichment of genes in androgen response with p-value of 7.51E-07 for 1912 gene set and decreased p-value of 5.23E-09 for 260 gene set). Similarly, pathways of glycolytic and glycerophospholipid metabolism showed consistently higher enrichment with the lowest p-values in the set of 260 genes, indicating that KDM3A demethylation targets these pathways and causes an up-regulation of gene expression. Of all genes detected by both H3K9 demethylation and AR ChIP-Seq experiments as well as KDM3A and AR differential expression following knockdown (core set of 45 genes; 28 genes up-regulated; 17 genes down-regulated), the AR response was the only significantly enriched pathway with a p-value of 5.29E-08. 7 genes, NDRG1, PTK2B, ACSL3, KRT19, INPP4B, NKX3-1, and MAF, with direct implication in AR signaling were present in all gene sets of the hierarchical enrichment analyses, thereby subject to KDM3A control by H3K9 demethylation, AR binding, and differential expression following knockdown (Supplementary Table 1 -3, 5). In the androgen response pathway KDM3A and AR had strong ChIP-Seq activity ( Table 4). The majority of epigenetic H3K9 demethylation events resulted in up-regulation of target gene activity and occurred at the TSSs of target genes. In addition, multiple introns, exons, and TTS of target genes were also implicated in KDM3A control. A network upstream regulator analysis implicated the AR to coordinate with epigenetic and transcriptional responses observed upon knockdown of KDM3A (Supplementary

DISCUSSION
Epigenetic regulators like KDM3A specifically demethylate histone marks H3K9me1 and H3K9me2, thereby playing a central role in the histone code. In cancer, demethylation and decondensation of chromatin can lead to dysregulated gene expression and transcriptional activation of gene targets [21]. Previous cell biological studies have suggested that KDM3A may stimulate transcription mediated by nuclear receptors and/or that KDM3A may be involved in activation of forkhead proteins during cell differentiation [17,22]. While histone modifiers are not limited to specific DNA cognate sites, transcriptional specificity can be accomplished by cooperation with transcription factors recognizing distinct DNA motifs [13].
The chosen bioinformatics approach of quantifying changes of matched epigenetic remodeling in combination with transcription factor ChIP-Seq and transcriptomic analysis following knockdown of an epigenetic master regulator offers insight from two angles: it can identify in an unbiased way genome-wide cooperation of epigenetic remodeling with other members of the transcriptional machinery, and it can elucidate details of the interaction. In prostate cancer, steroid ligand dependency or independency can influence the prognostic outcome. The CWR22Rv1 cell line offers the benefit of being able to probe a static, aggressive end point of the disease, while it has limitations due to lacking the dynamic ligand-dependent aspect of AR signaling. Since the CWR22Rv1 cell line expressed permanently activated AR without the ligand binding domain, it provides a stable model for studying the dynamic response to epigenetic regulation by KDM3A knockdown in combination with ChIP-Seq analysis. KDM3A knockdown abolished tumor formation in an orthotopic prostate tumor model using CWR22Rv1 cells [18]. Interestingly, KDM3A knockdown in other prostate cancer cell lines, including the androgen dependent LNCaP cells, blocked cell proliferation [18].
Androgen signaling is subject to multilevel control. In addition to agonist and antagonist ligand chaperones, intracellular localization and interactions with other transcription factors or histone modifiers can influence the transcriptional outcome. Multiple epigenetic regulators have been described to interact with the AR. KDM4C co-localizes with the androgen receptor in prostate carcinomas, and knockdown of KDM4C inhibits transcriptional activation and tumor cell proliferation [23]. KDM3A is involved in spermatogenesis by regulating expression of target genes such as PRM1 and TMP1, which are required for packaging and condensation of The analysis is focused on coactivation by KDM3A/AR ChIP and transcriptomic change upon knockdown. Input gene sets were defined by coactivation by KDM/AR ChIP-Seq (1912 genes; marked with *), coactivation by KDM/AR ChIP-Seq and transcriptomic change upon KDM3A/AR knockdown (421 genes; marked with **), oncogenic gene up-regulation by KDM3A/AR and KDM/AR ChIP coactivation (260 genes; marked with **) (compare Figure 4). Comprehensive lists of gene sets and enrichment ratios are deposited in the Supplementary Information. To determine significance of pathway enrichment thresholds of 0.05 and 0.10 were used for p-values and q-values, respectively according to multiple hypotheses testing using the controlling procedure of Benjamini and Hochberg. sperm chromatin [3]. Furthermore, its involvement in obesity resistance through regulation of metabolic genes such as PPARA and UCP1 highlight a transcriptional network focused on lipid modifiers. In this epigenetic and transcriptomic study we aimed to outline specific pathways of KDM3A demethylase action and enriched transcriptional networks under its control. Previous studies have outlined KDM3A expression levels in prostate cancer phenotypes [24], but KDM3A-regulated target pathways by ChIP-Seq analysis were unknown [24][25][26]. Motif-guided searches for cooperating transcription factors can link transcriptional programs with genome-wide histone modifications. Motifs of SREBP, HIF, AP1, KLF, MYC, and FOX families are enriched in the H3K9me1/2-KDM ChIP-Seq data and were described to play a role in prostate cancer progression [3,9,[27][28][29][30][31]. However, of the enriched transcription factors characterized so far, only AR and MYC have strong biochemical links to KDM3A [3,18]. Solely based on KDM3A regulation of ChIP-Seq and transcriptomic data, a gene set of 1408 genes revealed androgen-related signaling as top hit. Within the androgen response, transcriptionally validated genes, NKX3-1 and KLK3, have been shown to have dynamically regulated histone modification states [32].
In a second step, matched ChIP-Seq studies of KDM3A and the AR allowed us to focus on a distinct transcriptional network. Aside from the androgen response, cellular metabolism is highly enriched in the executioner program of KDM3A and the AR. In particular, lipogenic and hypoxic metabolism stand out for their regulation in the combined KDM3A/AR ChIP-Seq and transcriptomic data. Previously established AR targets involved in cellular metabolism include ACSL3, which is involved in converting free long chain fatty acids into fatty acyl-CoA esters [33,34]. Genome-wide ChIP-Seq profiles show that KDM3A regulates the expression of FKBP5, CHUK, HSP90AA1, and VHL. These genes were previously implicated as being regulated by the AR and their proteomic function is involved in AR folding, transactivation, and translocation in the nucleus [35][36][37][38]. HSP90AA1 is described to be under the control of KDM3A [39][40][41]. Prominent shared targets of KDM3A and AR include gene targets associated with hypoxic metabolism. Transcriptional regulation in response to hypoxia is regulated by the actions of HIF1A and controls glycolytic Transcriptomic validation of gene expression change of target gene following knockdown of KDM3A or AR is listed as LOG2 (fold change).
metabolism [42]. Previously, KDM3A gene expression was identified as one of the genes under control of HIF1A [43]. In addition, there is increased transcriptional activity of AR within castration-resistant prostate cancer by hypoxia [44]. A similar overlapping network of KDM3A demethylation, nuclear hormone signaling, and hypoxia is described in estrogen independent breast cancer models [16]. KDM3A's ability to regulate metabolic gene expression by controlling AR binding site availability in hypoxic cells may be the molecular action KDM3A utilizes to stimulate tumor progression [45]. The additional enriched pathways with shared KDM3A and AR regulation appear to the hypoxic cell response. Cytokine and metabolic signaling is able to induce expression of HIF1A and promote its activation in an oncogenic fashion [46][47][48][49][50][51].
The functional impact of coordinated action between a lysine demethylase and transcription factors may lend to its target specificity, or at the very least, create accessibility for DNA binding [52][53][54][55][56]. KDM3A controlled H3K9me1/2 ChIP-Seq data shows a strong enrichment of AR binding sites within the CWR22Rv1 castration-resistant prostate cancer cell line with continued expression of genes involved in the androgen response. However, about a third of the AR ChIP-Seq peaks were suppressed in response to increased H3K9me1/2 signal in the KDM3A experiments. It remains to be determined which other coactivators and corepressors take charge of AR binding sites in genomic locations where KDM3A has no effect on AR binding or transcriptional response. Notably, selected KDM3A-or AR-dependent genes show no clear association of modulated H3K9me1/2 marks within 50,000 bp around the gene body, suggesting that KDM3A can regulate gene expression either by accessing distant enhancers or by physical interaction with the transcriptional machinery independent of H3K9 demethylation. For example, KDM3A was identified to erase lysine 372 monomethylation of TP53, a protein methylation site crucial for the stability and pro-apoptotic function of chromatin-bound tumor suppressor [57]. Significantly, next to the AR, SREBF, and HIF, the conducted transcription factor motif analysis confirmed a strong presence of TP53 target sequences controlled by KDM3A. Extensive future ChIP-Seq studies of KDM3A as well as other candidate transcription factors associated with KDM3A will be necessary to further characterize the full spectrum of epigenetic and transcriptional control of the master regulator KDM3A.
In conclusion, the ChIP-Seq study refined the genomic sites of KDM3A-mediated H3K9me1/2 histone demethylation within the CWR22Rv1 prostate cancer cell line. The cancer systems biology analysis expanded underlying transcription factor motifs associated with oncogenic KDM3A demethylation, suggesting an underlying transcriptional network that directs transcriptional activation. Future experimental verification of epigenetic hotspots is needed to determine when detected response elements are functional in gene regulation. The matched transcriptomics and epigenomics approach identified an overlap between androgen receptor ChIP-Seq and KDM3A-regulated H3K9me1/2 ChIP-Seq. The comprehensive genome-wide mapping of matched ChIP-Seq profiles highlighted mechanistic details of how an epigenetic master regulator can exhibit control over selected transcriptional programs, such as metabolic pathways and hypoxia response in cancer.

Tissue culture of prostate cancer cell lines
CWR22Rv1 is a human prostate carcinoma epithelial cell line derived from a xenograft that was serially propagated in mice after castration-induced regression and relapse of the parental, androgen-dependent CWR22 xenograft [58,59] (CRL-2505, American Type Culture Collection, Manassas, VA). The CWR22Rv1 prostate cancer cell line was kindly provided by Dr. James Jacobberger (Case Western Reserve University, Cleveland, OH), and are maintained in RPMI 1640 medium supplemented with 10% fetal bovine serum and antibiotics. Cell cultures are regularly tested to ensure absence of mycoplasma. The CWR22Rv1 prostate cancer cell line expresses AR full-length with a duplicated DNA binding domain in exon 3 and AR splice variants, for example AR-v7, lacking a ligand binding domain. Thus, the CWR22Rv1 cell line displays constitutively active AR even in the absence of androgen [59]. In the ChIP-Seq study an AR antibody (PG21, 06-680, Sigma EMD Millipore, Darmstadt, Germany) was used recognizing both forms of the AR. All experimental protocols were approved by the Institutional Review Board at the University of California Merced. The study was carried out as part of IRB UCM13-0025 of the University of California Merced and as part of dbGap ID 5094 on somatic mutations in cancer and conducted in accordance with the Helsinki Declaration of 1975.

qRT-PCR analysis
Total RNA from prostate cancer cells was extracted using a mammalian RNA mini preparation kit (RTN10-1KT, GenElute, Sigma EMD Millipore, Darmstadt, Germany) and then digested with deoxyribonuclease I (AMPD1-1KT, Sigma EMD Millipore, Darmstadt, Germany). Complementary DNA (cDNA) was synthesized using random hexamers. Triple replicate samples were subjected to SYBR green (SYBR green master mix, Qiagen SABiosciences) qRT-PCR analysis in an Eco system (Illumina, San Diego). Gene expression profiles were analyzed using the ΔΔCT method. qRT-PCR threshold cycle (CT) values were normalized using the housekeeping gene cyclophilin A (PPIA; peptidylprolyl isomerase A; Gene ID: 5478). The following primers served for ChIP-qRT-PCR validation of ChIP-Seq signal of H3K9me1/2-KDM around the AREs of AR target genes: KLK3 (kallikrein related peptidase 3; Gene ID: 354; also known as PSA, prostate specific antigen):

Next generation sequencing and ChIP-Seq analysis
The purified DNA library was sequenced using an Illumina HighSeq2000 at the Sanford-Burnham Medical Research Institute at Lake Nona, National Genome Library Core Facility. Sequenced regions were aligned to the reference human genome 19 using the Bowtie alignment program that utilizes an extended Burrows-Wheeler indexing for an ultrafast memory efficient alignment [60]. Peak calling utilized a model-based analysis of ChIP-Seq (MACS) algorithm [61,62]. The overlap analysis, plot of genomic location, sequence extraction, motif identification, and peak filtering were performed using ChIPseek: a web-based analysis for ChIP data [63].
ChIPseek also employs scripts from BEDtools [64] using a genome binning algorithm used by the UCSC genome browser to sort genomic regions into groups along the length of chromosome [65]. Data visualization was carried out using the integrative genomics viewer [66].

Motif analysis based on position site-specific matrix models
Computational response element searching algorithms are able to estimate a sequence's likelihood in belonging to the response element of the query transcription factor using position site-specific matrix (PSSM) models where each position in the query transcription factor model gives each of the four letters in the DNA alphabet a score based on the probability of that nucleotide being found at that position [67]. Summation into a logs-odd score is converted into a p-value assuming a zero-order background model, and all response elements less than the threshold are reported [68]. Motif discovery, motif enrichment, and motif scanning used the multiple expectation maximization for motif elicitation (MEME) and discriminative regular expression motif elicitation (DREME) suite software toolkits from a set of user supplied unaligned sequences for ChIP-Seq regions [69]. De novo motif analysis programs MEME and DREME identify frequently detected DNA sequences patterns and similarity matches of recurring ChIP-Seq sequences with DNA motifs www.impactjournals.com/oncotarget of deposited studies in genomic sequence databases [68,70]. After a motif of interest is discovered the genomic sequences of the ChIP sequenced data is scanned using the MEME suite software find individual motif occurrences (FIMO) [68] for individual motif occurrences using PSSMs to compute a log-likelihood ratio score for each submitted sequence. Sequencespecific matrix models are further used to analyze the next generation sequencing data for motif enrichment and potential coactivators [13,71].

Microarray analysis
In order to quantify the transcriptomic effect of KDM3A or AR knockdown, a microarray profiling analysis was conducted on CWR22Rv1 knockdown cells. CWR22v1 cells were transduced with lentiviral pLKO.1 control shRNA vector, AR, or KDM3A shRNA for 48 h. Total RNA was isolated from cells, and 500 ng was used for synthesis of biotin-labeled cRNA using an RNA amplification kit (Ambion, Thermo Fisher Scientific, Waltham, MA). Biotinylated cRNA was labeled by incubation with streptavidin-Cy3 to generate a probe for hybridization with the GeneChip Human Transcriptome Array 2.0 (Affymetrix Inc, Santa Clara, CA) or the genome-wide transcriptome Human HT-12 V4.0 (Illumina Inc, San Diego, CA). Four samples from two experimental groups (n=2 per group) were hybridized to the chip to obtain raw gene expression data, which was processed to obtain raw data in the form of expression intensities. Raw data was then exported for further processing and analysis using R statistical software version 3.1 in combination with the BioConductor oligo, affy and genefilter packages [72]. The raw signal intensities were background corrected, LOG2 transformed, and quantile normalized to generate robust multi-array average (RMA) normalized intensities [73]. Quality control analyses did not reveal any outlier samples. Differential expression between experimental groups was assessed by generating relevant contrasts corresponding to the two-group comparison and was evaluated using the Linear Models for Microarray Analysis (LIMMA) package [72,74]. Raw p-values were corrected for multiple hypotheses testing using the false discovery rate controlling procedure of Benjamini and Hochberg, and adjusted p-values below 0.05 were considered significant [75]. Genes with significantly altered expression levels with adjusted p-values below 0.05 following KDM3A knockdown were selected and analyzed through the use of Ingenuity Pathway Analysis (IPA, Qiagen, Redwood City, CA). Pathway enrichment of differentially expressed genes was determined by gene set enrichment analysis (GSEA) by pairing each gene with its expression value and ranking genes in descending order (Broad Institute, Cambridge, MA) [76,77]. As we are testing multiple gene sets at the same time, the p-values need to be controlled for false positives. The false discovery rate estimation for the pathway enrichment is as summarized in q-values with a threshold of 0.10 controlling the error rate and correcting for multiple hypotheses testing according to Benjamini and Hochberg [75]. Acquired data of transcriptome profiling microarray analysis of CWR22Rv1 cells with AR knockdown using GeneChip Human Transcriptome Array 2.0, platform GPL16686 (Affymetrix Inc, Santa Clara, CA), is deposited in the NCBI GEO database under accession number GSE86547. Acquired data of CWR22Rv1 cells with KDM3A knockdown using hybridization with genome-wide transcriptome Human HT-12 V4.0, platform GPL10558 (Illumina Inc, San Diego, CA), is deposited under accession number GSE70498.

Availability of supporting data
The Supplementary Information contains tables on genome-wide mapping, annotation, and overlap of H3K9me1/me2 demethylation ChIP-Seq and AR ChIP-Seq (Supplementary Table 1), gene sets based on ChIP-Seq and transcriptomic data (Supplementary Table 2), pathway enrichment analysis based on H3K9me1/me2 demethylation ChIP-Seq gene set (Supplementary Table  3), transcriptomic response upon shRNA(KDM3A) and shRNA(AR) knockdown (Supplementary Table 4), hierarchical gene set enrichment analysis of identified KDM3A target genes (Supplementary Table 5), and an upstream regulator analysis based on ingenuity pathway analysis (Supplementary Table 6).

Author contributions
F.V.F. and J.Q. designed the study and directed next generation sequencing. S.W. and F.V.F. conducted the ChIP-Seq data analysis. L.F. and N.S. conducted transcriptomic validation of ChIP-Seq targets. S.W., J.Q., and F.V.F. performed data interpretation, wrote the text, and reviewed the final manuscript. ACKNOWLEDGMENTS F.V.F. is grateful for the support of grants CA154887 and CA176114 from the National Institutes of Health, National Cancer Institute. J.Q. is supported by grants CA154888 and CA207118 from the National Institutes of Health, National Cancer Institute. This work is supported by UC CRCC CRN-17-427258, GCR, and HSRI program grants.

CONFLICTS OF INTEREST
The authors declare that there is no competing interest as part of the submission process of this manuscript.