VDR regulation of microRNA differs across prostate cell models suggesting extremely flexible control of transcription

Abstract

The Vitamin D Receptor (VDR) is a member of the nuclear receptor superfamily and is of therapeutic interest in cancer and other settings. Regulation of microRNA (miRNA) by the VDR appears to be important to mediate its actions, for example, to control cell growth. To identify if and to what extent VDR-regulated miRNA patterns change in prostate cancer progression, we undertook miRNA microarray analyses in 7 cell models representing non-malignant and malignant prostate cells (RWPE-1, RWPE-2, HPr1, HPr1AR, LNCaP, LNCaP-C4–2, and PC-3). To focus on primary VDR regulatory events, we undertook expression analyses after 30 minutes treatment with 1α,25(OH)₂D₃. Across all models, 111 miRNAs were significantly modulated by 1α,25(OH)₂D₃ treatment. Of these, only 5 miRNAs were modulated in more than one cell model, and of these, only 3 miRNAs were modulated in the same direction. The patterns of miRNA regulation, and the networks they targeted, significantly distinguished the different cell types. Integration of 1α,25(OH)₂D₃-regulated miRNAs with published VDR ChIP-seq data showed significant enrichment of VDR peaks in flanking regions of miRNAs. Furthermore, mRNA and miRNA expression analyses in non-malignant RWPE-1 cells revealed patterns of miRNA and mRNA co-regulation; specifically, 13 significant reciprocal patterns were identified and these patterns were also observed in TCGA prostate cancer data. Lastly, motif search analysis revealed differential motif enrichment within VDR peaks flanking mRNA compared to miRNA genes. Together, this study revealed that miRNAs are rapidly regulated in a highly cell-type specific manner, and are significantly co-integrated with mRNA regulation.

Keywords:

• epigenetic regulation
• microRNA
• prostate cancer
• VDR
• Vitamin D

Introduction

Aberrant transcriptional regulation is evident in many cancers. This arises by genomic, epigenomic, and posttranscriptional mechanisms. Genomic alterations, namely mutations, insertions, deletions, and amplifications, combine with aberrant epigenetic states to manipulate and distort the global cellular transcriptional program. In this manner, the cancer cell transcriptome forms a unique fingerprint of the cell state, which provides for potential therapeutic targeting strategies. Overall, the cancer cell state can be considered to be in a transcriptionally rigid state (reviewed in¹).

This restriction in transcriptional potential is displayed by members of the nuclear receptor (NR) superfamily. Members of this superfamily are expressed in virtually all cell types and are highly integrated functionally.^2-5 These receptors represent some of the most successful examples of targeted therapies^6-9 and, more widely, in human disease they represent the target for approximately 15% of all drugs (reviewed in^10,11). Of the family of 48 NRs, at least 30 are expressed in the prostate² and many are well established in the etiology of prostate cancer (CaP), for example, the androgen receptor (AR) (reviewed in¹²), VDR (reviewed in⁹), and the peroxisome proliferator-activated receptors (PPARs) (reviewed in⁶). Regulation studies on single NR target genes have revealed dynamic mRNA accumulation patterns.^13-16 AR regulation of the TMPRSS2/ETS fusion, a common occurrence in CaP,¹⁷ is critical, precisely because TMPRSS2 is sustained in an AR responsive state in CaP,¹⁸ whereas the regulation of other AR targets, such as the prostate tumor suppressor NKX3.1, becomes attenuated.¹⁹

We have established patterns of transcriptional distortion by the VDR in CaP cell lines and primary cultures, with attenuation of genes, such as GADD45A, whereas the DNA helicase MCM-7 and 1α,25(OH)₂D₃-degrading enzyme CYP24A1 are sustained in a responsive state.^20,21 Similar distortions to PPAR signaling occurs.^2,3,22,23 These de-regulated responses appear particularly amenable to targeted restoration with epigenetic therapies.^2-4,21,24-26

Furthermore, the transition to more aggressive CaP is associated with changes to the basal expression of multiple miRNAs, including miR-125b, which targets NCOR2/SMRT,^27,28 and miR-221/miR-222, which target CDKN1B (encodes p27^(kip1)).²⁹ An emergent literature has identified that AR and VDR also regulate miRNAs.^29-33 In parallel to the distorted NR regulation of mRNA in cancer, there is similar altered regulation of miRNA expression, for example, aberrant modulation of miR-21 and others by the AR³⁰ and VDR.³⁴ miRNA and mRNA combine to form feed forward loops. For example, we identified VDR-regulated MCM-7 that encodes the MIR106b cluster³⁵. VDR also regulated CDKN1A,^35-38 which, in turn, is targeted by miR-106b. The net result is fine tuning of p21^(waf1/cip1) expression and the control of the cell cycle arrest.³⁵

Together, these observations support the concept that miRNA regulation by the VDR is distorted in a similar fashion as the mRNA targets. Thus, the transcriptional control of miRNAs appears as precise as that required for mRNA and most likely reflects their function to serve as critical regulatory components in feed forward signaling loops.^39-43

miRNAs hold considerable promise to be exploited as highly accurate and functional prognostic serum markers of CaP stages and drug responses.^27,44-54 They can encapsulate events within the tumor microenvironment, thereby overcoming the limitations of sampling error at biopsy. In parallel, it has become apparent that changes in miRNA regulation in tumors is translated to altered serum expression patterns and, therefore, offers an important diagnostic and prognostic therapeutic window.^{53, 55-57}

Therefore, we reasoned that part of the diversity of cell responses to 1α,25(OH)₂D₃ may be explained by the seemingly large choice of VDR binding sites through the genome and the downstream regulation of miRNA. Specifically, we propose that miRNA expression reflects the epigenetic restriction of VDR transcriptional responsiveness that evolves during CaP progression. In the current study, we applied microarray approaches to examine the VDR regulated expression of miRNA in a wide panel of human prostate non-malignant and malignant cell lines, combined with parallel microarray analyses of mRNA in non-malignant prostate RWPE-1 cells. Combining these data with publically available VDR ChIP-seq data indicates how VDR transcriptional actions evolve in cancer. The findings suggest selective distortion of miRNA regulation in CaP progression.

Results

1α,25(OH)₂D₃ regulates differential expression of miRNA in 7 different prostate cell models

We reasoned that short time exposures of cells to 1α,25(OH)₂D₃ would reveal genes where the VDR was either already bound or could access them as they were in open chromatin context. Therefore, cells were treated with 1α,25(OH)₂D₃ for 30 minutes and miRNA extracted for microarray analyses.

Characteristics and 1α,25(OH)₂D₃ sensitivity of the cell line models are listed in Table 1. We observed distinct patterns of VDR-regulated miRNA expression. In total, 111 miRNA were significantly modulated (P-value < 0.01, Table 1 and Table S1) in different cell line models. Of these, only 5 miRNAs were modulated in more than one model, and of these, only 3 miRNAs were modulated in the same direction (miR-30e, mir-98, and miR-140–5p; Table 2). The microarray expression patterns in the 7 prostate cell lines are represented as a heatmap in Figure 1A. The cell lines with similar clinical features and 1α,25(OH)₂D₃ sensitivity showed related expression changes and clustered together. Generally, non-malignant epithelial cells with the greatest 1α,25(OH)₂D₃ sensitivity (RWPE1 and HPr1AR) clustered separately from androgen-sensitive and partially 1α,25(OH)₂D₃-sensitive malignant models (LNCaP and RWPE2) and ADT-RCaP disease (Fig. 1A). These data suggest that VDR-regulated miRNA expression patterns have the capacity to distinguish between cell lines derived from different stages of CaP. Microarray results were validated using Q-RT-PCR for 5 miRNAs and Q-RT-PCR results showed concordance with microarray data (Fig. 1B).

Table 1. miRNA expression analyses in prostate cell lies. Each cell line was treated with 1α,25(OH)₂D₃ (100 nM) and ethanol as vehicle control for 30 minutes. Log2 changes in expression of miRNAs were calculated in 1α,25(OH)₂D₃ treated samples compared to ethanol treated samples. Number of up- and down-regulated miRNAs in response to 1α,25(OH)₂D₃ treatment in each cell line are shown. Characteristics and reported response to 1α,25(OH)₂D₃ for each cell line is also shown.

Characteristics	Cell lines	No. of upregulated miRNAs	No. of downregulated miRNAs	Total miRNAs	1α,25(OH)2D3 responsiveness
Non-malignant	RWPE-1	5	5	10	Responsive ^3
Non-malignant	HPr1	11	10	21	Unknown
Non-malignant	HPr1AR	8	5	13	Unknown
Early stage CaP	RWPE-2	1	2	3	Moderately responsive ^66
Androgen sensitive CaP	LNCaP	12	5	17	Partially responsive ^79
ADT recurrent CaP	LNCaP-C4–2	5	4	9	Recalcitrant ^79
Metastatic CaP	PC3	21	22	43	Recalcitrant ^3

Table 2. miRNAs regulated in more than one cell line. Five miRNAs showed differential expression in more than one cell line upon 1α,25(OH)₂D₃ treatment. Table shows the regulated miRNAs, cell lines and their direction of expression change compared to ethanol within each cell line. Two arrows show the direction of expression change in respective cell lines (↑) up-regulated, (↓) downregulated.

miRNA	Cell line	Direction
miR-98	LNCaP & PC3	↑↑
miR-92a	LNCaP & HPr1AR	↑↓
miR-30e	LNCaP & HPr1	↑↑
miR-140–5p	HPr1 & PC3	↑↑
miR-138	HPr1 & HPr1AR	↑↓

Figure 1 (See previous page). Regulation of miRNAs in response to 1α,25(OH)₂D₃ in prostate cell lines. (A) Heatmap of miRNA expression profile across 7 prostate cell lines. The cell lines were treated with either 1α,25(OH)₂D₃ or ethanol for 30 minutes and genome-wide miRNA expression analysis was performed using Exiqon microarrays. Differential expression of miRNAs was calculated in treated samples compared to ethanol vehicle control samples for each cell line. Hierarchical clustering of 7 prostate cell lines was performed using log2 fold changes values of 111 miRNAs across all the cell lines. Hierarchical clustering was performed using Cluster 3.0 and the clusters were visualized by Tree View. Clustering of cell lines showed the separation of cell lines based on their phenotypic characteristics. (B) To validation the microarray results 5 miRNAs across 4 different cell lines were tested using miRNA specific TapMan probes. Validation by Q-PCR showed concordant results for all 5 miRNAs tested. (C) Pathways analysis of miRNAs target genes for significant miRNAs in each cell lines was performed using DIANA-mirPath. Top 20 significantly enriched miRNA targeted pathways for individual cell lines were selected and a unique list of pathways was generated. Pathways were given a score of 1 or 0 if they are present or absent, respectively, in each cell line and hierarchical clustering was performed as described for miRNA expression. Hierarchical clustering of top 20 significantly enriched miRNA targeted pathways for each cell line showed separation of non-malignant and malignant cell lines. Pink color shows the presence of the pathway in the cell lines whereas gray color show absence of the pathway.

Subsequently, we analyzed these 1α,25(OH)₂D₃-regulated miRNAs across the cell models to identify the molecular pathways and networks targeted by these miRNAs using DIANA miRPath v.2.0. Firstly, lists of significantly targeted pathways by miRNAs were generated for each cell line. Individual cell lines showed unique and common pathways targeted by these miRNAs. Subsequently, the top 20 most significant pathways from each cell line were used to find if the common and unique pathway signatures cluster cell lines based on clinical features or 1α,25(OH)₂D₃ sensitivity. The top 20 significant pathways from all the cell lines represent 63 unique pathways. Hierarchical clustering of cell lines based on these 63 pathways showed very unique patterns where non-malignant cell lines (RWPE-1, HPr1 and HPr1AR) and cancer cell lines (LNCaP, LNCaP-C4–2 and PC3) cluster separately (Fig. 1C). This pathway analysis and clustering suggests that 1α,25(OH)₂D₃-regulated miRNAs target pathways that are shared by either malignant or non-malignant cell lines.

VDR binding associates with miRNA regulation. In an attempt to establish a role for VDR regulation of these miRNAs, we leveraged publically available VDR ChIP-seq data to identify VDR binding regions in one or more data sets within 5, 10, 50, and 100 kb genomic regions flanking the miRNAs (Table 3). Specifically, we exploited a pre-compiled list of 23,409 non-overlapping VDR binding sites⁵⁸ derived from published VDR ChIP-seq studies. We used this list of VDR binding sites to identify VDR binding regions in the flanking sequences of miRNAs. The VDR binding peaks were compared to the genomic coordinates for primary miRNAs downloaded from miRbase. To test whether 1α,25(OH)₂D₃-regulated miRNAs in CaP cell lines show enrichment of VDR peaks we compared the number of peaks around 1α,25(OH)₂D₃ regulated miRNAs (N = 111) to the number of VDR peaks around all the miRNAs in miRbase (N = 2571). The extent of significant enrichment was tested using hypergeometric distribution.

Table 3. Showing distribution of VDR peaks near miRNAs regulated by 1α,25(OH)₂D₃ in prostate cell lines. Genomic coordinates of all the regulated miRNAs across all the prostate cell lines were downloaded from miRbase and overlapped with combined list of VDR peaks from 5 published VDR ChIP-seq data sets. Number of VDR peak in defined flanking regions of 1α,25(OH)₂D₃ regulated miRNAs were compared to number of VDR peak in similar flanking regions of all the miRNAs reported in miRbase. Statistical significance of enrichment of VDR peaks was calculated using hypergeometric test and P ≤ 0 .05 was considered significant.

Distance	No, of 1α,25(OH)2D3 regulated miRs (N=111 )	All miRs (N=2571)	Hypergeomatric distribution (P-value)
±5 kb	24	399	0.0209
±10 kb	31	747	0.0829
±50 kb	51	1442	0.0070
±100 kb	68	1801	0.0103

The 1α,25(OH)₂D₃-regulated miRNAs showed significant enrichment (hypergeometric distribution P-value < 0.05) of VDR peaks compared to all the miRNAs reported in miRbase database. For example, out of 111 1α,25(OH)₂D₃-regulated miRNAs, 24 miRNAs (21.6%) have at least one VDR peak within ±5 kb, whereas out of 2571 reported miRNAs, 399 miRNAs (15.5%) have at least one VDR peak within ±5 kb (hypergeometric distribution; P-value < 0.021; Table 3). Similarly, distribution of VDR peaks within ±50 and 100 kb of 1α,25(OH)₂D₃-regulated miRNAs also showed significant enrichment compared to all reported miRNAs (Table 3). These results suggest that direct VDR binding events proximal to these miRNAs are functionally relevant in modulation of miRNA expression patterns in response to 1α,25(OH)₂D₃.

Significantly altered mRNA regulation in RWPE-1 cells reflects miRNA regulation

To complement these miRNA analyses, we undertook a parallel mRNA microarray experiment in RWPE-1 cells in response to 1α,25(OH)₂D₃ treatment. This approach, again at a short time point (30 minutes), revealed a restricted transcriptome of 108 targets (P-value < 0.01; Table S2). Of these, 48 genes were upregulated and 60 genes were downregulated. Gene ontology analysis with PANTHER⁵⁹ identified several biological processes to be enriched in this gene list, including cellular and metabolic processes (Fig. 2A).

Figure 2 (See previous page). Functional annotation of 1α,25(OH)₂D₃-regulated mRNA genes in RWPE-1 cells and their co-regulation with miRNAs. (A). Biological processes enriched in the set of genes modulated by 1α,25(OH)₂D₃ in RWPE-1 cells were identified through PANTHER gene list analysis tool. The graph shows significantly enriched biological process and the number of genes in each process which is regulated. (B) The miRNAs and their predicted targets co-regulated in RWPE-1 in response to 1α,25(OH)₂D₃. Gene and miRNAs pairs showing negative correlation and their expression levels are shown. (C) Shows the reciprocal expression of miR-20a and its predicted target gene NOG. miRNA and mRNA expression (RNA-seq) data was downloaded from TCGA data portal. Raw expression values were log2 transformed and t-test was performed between normal tissue and tumors to identify significant difference in expression after FDR correction. Box plot show that miR-20a is upregulated in tumors where and its predicted target gene NOG is down regulated in tumors suggesting reciprocal relation between miRNA and its target mRNA. (D) Shows the enriched DNA binding motifs under the VDR peaks which does not contain DR3-type motif and present within ±5 kb of mRNA and miRNAs genes. Motif search was performed using MEME-ChIP tool of MEME suite.

To explore the coordinated regulation of miRNAs and their target mRNAs we integrated miRNA and mRNA microarray data from RWPE-1 cells. First, we identified predicted mRNA targets for differentially expressed miRNAs (P ≤ 0.05, Table S3) in RWPE-1 cells. For this study, mRNA targets predicted by at least 3 of 4 different prediction algorithms (miRanda,⁶⁰ miRDB,⁶¹ TargetScan⁶², and ExprTargetDB⁶³) were considered for further analysis. Subsequently, predicted mRNA targets that were differentially expressed upon 1α,25(OH)₂D₃ treatment in RWPE-1 cells were used to identify reciprocal expression patterns between miRNAs and their predicted mRNA targets. In this approach, we identified 10 miRNAs showing a statistically significant inverse expression pattern with 13 genes regulated by 1α,25(OH)₂D₃ (Fig. 2B). These findings suggest that in non-malignant RWPE-1 cells, 1α,25(OH)₂D₃-regulated miRNA and mRNA expression are integrated.

We also tested the expression of these miRNAs and mRNAs in TCGA prostate cancer samples using RNA-seq data. Interestingly, out of these 13 genes, 5 genes showed differential expression between tumor and matched normal samples (downregulated genes—WIF1 and NOG; upregulated genes—KIFC1, KLHL6, and SLC34A3). Also miR-20a showed differential expression in TCGA CaP samples. Further, a reciprocal expression pattern was observed between miR-20a and NOG, a predicted target of miR-20a, in TCGA CaP samples (Fig. 2C). These findings in TCGA CaP samples further strengthen the finding of reciprocal expression between mRNA and miRNA in RWPE-1 cells.

VDR binding patterns around miRNA and mRNA

We subsequently examined the role of VDR binding associated with mRNA and miRNA targets in RWPE-1 cells. In the first instance, we compared all the VDR peaks within ±5 kb of 1α,25(OH)₂D₃-regulated miRNAs to identify peaks containing the canonical DR3-type VDR binding motif. For 1α,25(OH)₂D_3-regulated miRNAs, out of 81 VDR peaks within ±5 kb, 14 (17.3%) peaks contained the DR3-type motif. However, the same analyses of 1α,25(OH)₂D_3-regulated mRNA genes revealed 183 VDR peaks within ±5 kb and that 66 (36.1%) peaks contained the DR3-type motif. This difference of DR3-type containing motifs between 1α,25(OH)₂D₃-regulated miRNAs and mRNAs was statistically significant (Fisher's exact test P-value = 0.0022).

To identify other transcription factors which may be interacting with VDR to regulate these miRNAs and mRNA genes we searched for binding motif enrichment in VDR peaks that do not contain DR3-types motif using Multiple Em for Motif Elicitation (MEME).⁶⁴ The binding motif of ZNF263 was significantly enriched in VDR peaks flanking both mRNA and miRNA genes. Whereas the binding motif of IRF1/STAT1,2 was only significantly enriched under the VDR peaks flanking mRNA genes (Fig. 2D).⁶⁵ This difference suggests that mRNA and miRNA genes are under control by unique VDR-containing complexes.

Discussion

The current study was undertaken to define the extent and flexibility of the VDR-regulated miRNAs in a wide range of different prostate cell models from non-malignant to malignant (androgen sensitive and androgen deprivation therapy recurrent) and metastatic models. Remarkably, the expression of VDR-regulated miRNAs was strikingly flexible and suggests that the VDR actions in this disease context are highly modulated by other cellular forces. These data also suggest that patterns of VDR-regulated miRNA and mRNA are integrated as identified in RWPE-1 cells. Furthermore, analyses of TCGA CaP data provide evidence that this integration remains in the diseased state.

We examined the expression of miRNA in response to 1α,25(OH)₂D₃ at an early time point (30 minutes) with the goal to identify initial and direct responses toward 1α,25(OH)₂D₃ and direct VDR targets in prostate cell lines. We identified differential miRNA expression with minimal overlap between cell lines supporting the concept that VDR regulation of miRNA is highly variable. As with other NRs, the genomic actions of VDR oscillate between on and off states and result in cyclical gene regulation that is apparent at early time points and often continues to modulate for at least several hours, subsequently. These findings suggest that the differentially expressed miRNA and mRNA signatures identified in this study after 30 minutes treatment are strongly suggestive of genes that will continue to modulate at later time points and have true biological consequences.^{35, 66}

This flexibility of VDR signaling is also reflected by the 5 published VDR ChIP-seq data sets, which show very little conservation of the choice of VDR binding site across the genome in different cell types. These five studies identified roughly a similar number of VDR binding sites (approximately 1500 – 2500) despite examining different cell lines at different time points. Namely, lymphoblastoid cells after 36 h treatment, ⁶⁷ monocytes (THP-1) after 40 minutes treatments,⁶⁸ and colorectal cancer cells (LS180) after 180 min with 1α,25(OH)₂D₃.⁶⁹ Remarkably, only 20% of these genomic sites are identical in all 3 investigated cell lines. This indicates that the majority of VDR binding is both time and cell background specific. In general, this suggests that for the VDR, and probably also for other nuclear receptors, there is a high number of binding sites available within the genome and that these potential binding sites are distilled within a cell. The basis for this unrealized potential remains ambiguous but probably reflects a diversity of receptor interactions. These findings are consistent with the concept that VDR interacts with other “pioneer factors” or transcription factors that bind regulatory genomic regions and facilitate the opening of these loci via interactions with chromatin modifying enzymes. This then allows “following factors” to bind a subset of these accessible regions and to execute their regulatory actions. Analyses of the enrichment of the DR-3 motifs in miRNA compared to mRNA support the concept of differential targeting of protein coding versus non-coding genes. Similarly the enrichment of the IRF1 motif in the mRNA further supports the concept of unique regulation of genes depending on cell function.

As to what, in more aggressive cancer cells, controls the ability of the VDR to regulate miRNA, it is worth noting that in the absence of ligand, the VDR is resident in the nucleus. This in turn may support a mechanism by which inappropriate interactions with co-repressors and H3K9 methylation favor a more stable epigenetic event, namely, DNA methylation. We have revealed recently that the increased recruitment of co-repressor association with the VDR appears to convert a transient epigenetic silencing process, which is part of the normal nuclear receptor-transcriptional cycle, into a stable and heritable epigenetic event.⁶⁶ The consequences of this could therefore be the targeted methylation of genes where NCOR1 and other co-repressors are recruited by the VDR, or other nuclear resident transcription factors. Therefore, these receptors may provide a route for the silencing of critical miRNA transcriptional control by selective NCOR1 recruitment and thereby allow CaP cells to escape mitotic restraint.

Tumor expression patterns of miRNAs have been shown to distinguish either cancer subtypes or predict tumor biological and clinical behavior, and recently it has been demonstrated that altered levels of miRNA as a consequence of tumor expression are detectable in serum.²⁷ Thus, an emergent and highly attractive area in cancer biology is exploiting these non-invasive biomarkers to characterize clinical cancer phenotypes and disease progression more accurately. The current study addressed how miRNAs are regulated by the VDR pathway and how these actions are distorted in CaP phenotypes. Thus, by considering miRNA expression behavior in the presence of different vitamin D levels potentially gives an accurate and practical route to stratify patient responses toward VDR-centered therapies.

Materials and Methods

Cell Lines and reagents: Prostate epithelial cell lines representing a spectrum of CaP phenotypes were used, as described previously.²⁷ Cell lines derived from normal prostate epithelial cells (RWPE-1, HPr1, HPr1AR), early stage CaP cells (RWPE-2), androgen sensitive CaP (LNCaP), ADT recurrent CaP (LNCaP-C4–2), and metastatic CaP (PC3) were used. RWPE-1, RWPE-2, HPr1, and HPr1AR cells were maintained in KSF media supplemented with EGF and BPE (Invitrogen). Other cells were maintained in RPMI-1640 supplemented with 10% v/v FBS (Invitrogen). 1α25(OH)₂D₃ [gift of Dr. Milan Uskokovic (BioXell S.p.A., Italy)] was used at 100 nM concentration for all experiments and compared to ethanol treated cells.

Microarray analyses: Total mRNA and miRNA isolation was performed using Trizol (Invitrogen). MiRNA expression analysis in all the cell lines was performed using Exiqon miRNA 5th generation expression arrays using single color hybridization. For mRNA expression Illumina microarray (Illumina HT12v4) analysis was performed in RWPE-1 cells. All the microarrays were performed in biological triplicates for each cell line treated with 1α,25(OH)₂D₃ (100 nM) or ethanol (vehicle control) for 30 minutes. All analyses were performed using the R platform for statistical computing,⁷⁰ and a range of library packages implemented in Bioconductor, including the LIMMA package.⁷¹ Normalized expression of miRNA and mRNA microarrays was compared across cell lines using Linear Models for Microarray Data in LIMMA. Hierarchical clustering was performed using Cluster 3.0⁷² and the clusters were visualized by Tree View.⁷³ Gene Ontology analysis of 1α,25(OH)₂D₃-regulated genes in RWPE-1 was performed using PANTHER (Protein ANalysis THrough Evolutionary Relationships) classification system.⁵⁹ Network analysis of differentially expressed miRNAs was undertaken to reveal common targeting of pathways. To understand the functional significance of these miRNAs with CaP disease stages, we undertook pathway enrichment analyses using Diana-mirPath.⁷⁴

Quantitative real time PCR (Q-RT-PCR): For validation of microarray results, TaqMan probes for selected miRNAs were used and Q-RT-PCR analysis was performed on ABI 7300.

Integration with VDR ChIP-seq data: Tuoresmaki et al.⁵⁸ compiled data sets from 5 published VDR ChIP-seq studies including immortalized lymphoblastoids,⁶⁷ hepatic stellate cells,⁷⁵ and cancer cell lines representing monocytic leukemia,^58,68 and colon cancer.⁶⁹ These studies differed in the time of cell exposure to 1α,25(OH)₂D₃ and identified binding sites, on the order of hundreds to several thousands of different binding sites per cell type.

From these 5 VDR ChIP-seq data sets a total of 23,409 non-overlapping genomic VDR binding sites were identified.⁵⁸ We used this list of VDR peaks to map VDR binding near miRNA and mRNA genes. Genomic coordinates of all miRNA primary transcripts were downloaded from miRbase (www.mirbase.org). Presence of VDR binding peaks was identified using Galaxy^76-78 within ±5, 10, 50 and 100 kb of miRNA primary transcripts for all 1α,25(OH)₂D₃-regulated miRNAs. Statistical significance of overlap was calculated by hypergeometric tests using R. The MEME Suite web server was used to identify enriched motif under VDR peaks. Specifically, MEME-ChIP⁶⁴ tool with JASPAR Vertebrates and UniPROBE Mouse motif database and default settings were used.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Funding

MJC acknowledges support in part from National Institute of Health grants [R01 CA095367–06 and 2R01-CA-095045–06] and the Prostate program of the Department of Defense Congressionally Directed Medical Research Programs [GRANT11496202], MJC also acknowledges support, in part, of the NCI Cancer Center Support Grant to the Roswell Park Cancer Institute [CA016056]. MDL acknowledges support of Molecular Pharmacology and Experimental Therapeutics NRSA T32 program [T32CA009072]. MJC and LES acknowledge support of the Pilot Studies program of the University at Buffalo Clinical and Translational Research Center and the Buffalo Translational Consortium.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.