extfieldsofstudy:
Medicine
Biology
Computer Science

text:
Estimating drivers of cell state transitions using gene regulatory network models

Background Specific cellular states are often associated with distinct gene expression patterns. These states are plastic, changing during development, or in the transition from health to disease. One relatively simple extension of this concept is to recognize that we can classify different cell-types by their active gene regulatory networks and that, consequently, transitions between cellular states can be modeled by changes in these underlying regulatory networks. Results Here we describe MONSTER, MOdeling Network State Transitions from Expression and Regulatory data, a regression-based method for inferring transcription factor drivers of cell state conditions at the gene regulatory network level. As a demonstration, we apply MONSTER to four different studies of chronic obstructive pulmonary disease to identify transcription factors that alter the network structure as the cell state progresses toward the disease-state. Conclusions We demonstrate that MONSTER can find strong regulatory signals that persist across studies and tissues of the same disease and that are not detectable using conventional analysis methods based on differential expression. An R package implementing MONSTER is available at github.com/QuackenbushLab/MONSTER. Electronic supplementary material The online version of this article (doi:10.1186/s12918-017-0517-y) contains supplementary material, which is available to authorized users.

Background
Cell state phenotypic transitions, such as those that occur during development, or as healthy tissue transforms into a disease phenotype, are fundamental processes that operate within biological systems. Understanding what drives these transitions, and modeling the processes, is one of the great open challenges in modern biology. One way to conceptualize the state transition problem is to imagine that each phenotype has its own characteristic gene regulatory network, and that there are a set of processes that are either activated or inactivated to transform the network in the initial state into one that characterizes the final state. Identifying those changes could, in principle, help us to understand not only the processes that drive the state change, but also how one might intervene to either promote or inhibit such a transition.
Each distinct cell state consists of a set of characteristic processes, some of which are shared across many cell-states ("housekeeping" functions) and others which are unique to that particular state. These processes are controlled by gene regulatory networks in which transcription factors (and other regulators) moderate the transcription of individual genes whose expression levels, in turn, characterize the state. One can represent these regulatory processes as a directed network graph, in which transcription factors and genes are nodes in the network, and edges represent the regulatory interactions between transcription factors and their target genes. A compact representation of such a network, with interactions between m transcription factors and p target genes, is as a binary p × m "adjacency matrix". In this matrix, a value of 1 represents an active interaction between a transcription factor and a potential target, and 0 represents the lack of a regulatory interaction.
When considering networks, a cell state transition is one that transforms the initial state network to the final state network, adding and deleting edges as appropriate. Using the adjacency matrix formalism, one can think of this as a problem in linear algebra in which we attempt to find an m × m "transition matrix" T, subject to a set of constraints, that approximates the conversion of the initial network's adjacency matrix A into the final network's adjacency matrix B, or B = AT (1) In this model, we describe the differences between cell states with a lower dimensional transition matrix. This matrix allows for the estimation of a relatively smaller number of parameters which focus on larger systemic shifts in regulatory behavior by TFs. Intuitively, one might recognize that the true transition matrix between identical network states is the identity matrix because the diagonal elements of T map network edges to themselves. Deviations from this identity, specifically the observation of meaningful non-zero values off of the diagonal, provide evidence of changes in regulatory network configuration for TFs between states.
While this framework, as depicted in Fig. 1, is intuitive, it is a bit simplistic in that we have cast the initial and final states as discrete. However, the model can be generalized Fig. 1 Overview of the MONSTER approach, as applied to the transition between smokers and those suffering from chronic obstructive pulmonary disease (COPD). MONSTER's approach seeks to find the TF × TF transition matrix that best characterizes the state change in network structure between the initial and final biological conditions. Subjects are first divided into two groups based on whether they have COPD or are smokers that have not yet developed clinical COPD. Network inference is then performed separately on each group, yielding a bipartite adjacency matrix connecting transcription factors to genes. Finally, a transition matrix is computed which characterizes the conversion from the consensus Smokers Network to the COPD Network by recognizing that any phenotype we analyze consists of a collection of individuals, all of whom have a slightly different manifestation of the state, and therefore a slightly different active gene regulatory network. Practically, what that means is that for each state, rather than having a network model with edges that are either "on" or "off, " a phenotype should be represented by a network in which each edge has a weight that represents an estimation of its presence across the population. In other words, the initial and final state adjacency matrices are not comprised of 1's and 0's, but of continuous variables that estimate population-level regulatory network edge-weights. Consequently, the problem of calculating the transition matrix is generalized to solving B = AT + E, where E is an p × m error matrix. In this expanded framework, modeling the cell state transition remains equivalent to estimating the appropriate transition matrix T, and then identifying state transition drivers based on features of that matrix.

MONSTER: MOdeling Network State Transitions from Expression and Regulatory data
The MONSTER algorithm models the regulatory transition between two cellular states in three steps: (1) Inferring state-specific gene regulatory networks, (2) modeling the state transition matrix, and (3) computing the transcription factor involvement.
Inferring state-specific gene regulatory networks: Before estimating the transition matrix, T, we must first estimate a gene regulatory starting point for each state. While there have been many methods developed to infer such networks [1][2][3][4][5][6][7], we have found the bipartite framework used in PANDA [8] to have features that are particularly amenable to interpretation in the context of state transitions. PANDA begins by using genomewide transcription factor binding data to postulate a network "prior", and then uses message-passing to integrate multiple data sources, including state-specific gene co-expression data.
Motivated by PANDA, we developed a highly computationally efficient, classification-based network inference method that uses common patterns between transcription factor targets and gene co-expression to estimate edges and to generate a bipartite gene regulatory network connecting transcription factors to their target genes.
This approach is based on the simple concept that genes affected by a common transcription factor are likely to exhibit correlated patterns of expression. To begin, we combine gene co-expression information with information about transcription factor targeting derived from sources such as ChIP-Seq or sets of known sequence binding motifs found in the vicinity of genes. The process of building prior networks to use as input to MONSTER may be complex, but our tool is agnostic to the source of this input. Users of MONSTER should use domain specific knowledge to generate an appropriate prior network.
We then calculate the direct evidence for a regulatory interaction between a transcription factor and gene, which we define as the squared partial correlation between a given transcription factor's gene expression, g i , and the gene's expression, g j , conditional on all other transcription factors' gene expression: where g i is the gene which encodes the transcription factor TF i , g j is any other gene in the genome, and TF j is the set of gene indices corresponding to known transcription factors with binding site in the promoter region of g j . The correlation is conditioned on the expression of all other potential regulators of g j based on the transcription factor motifs associated with g j . The direct evidence is motivated by the idea that changes in transcription factor expression may lead to similar changes in in target gene expression. The coexpression of co-targeted genes is long established in the literature [9,10], and evidence also points to the coexpression of transcription factor genes with targets of that transcription factor [11,12]. Moreover, studies across multiple tissues have shown widely varying expression of transcription factor genes, indicating that this expression can be used to predict their regulatory involvement [13]. Naturally, transcription factor behavior depends on many factors, including those that occur after translation. However, it makes intuitive sense that the mRNA abundance of a gene for a transcription factor should correlate with target genes to some degree. Next, we fit a logistic regression model which estimates the probability of each gene, indexed j, being a motif target of a transcription factor, indexed i, based on the expression pattern across the n samples across p genes in each phenotypic class: where the response M is a binary p × m matrix indicating the presence of a sequence motif for the i th transcription factor in the vicinity of each of the j th gene. And where g (k) j represents the gene expression measured for sample k at gene j. Thus, the fitted probabilityθ i,j represents our estimated indirect evidence. Combining the scores for the direct evidence,d i,j , and indirect evidence,θ i,j , via weighted sum between each transcription factor-gene pair yields estimated edge-weights for the gene regulatory network. We score each gene according to the strength of indirect evidence for a regulatory response to each of the transcription factors and combine this with the direct evidence of regulation. Combining our measures of direct and indirect evidence presents some challenges. Though both are bounded by [0,1] their interpretations are quite different. The direct evidence can be considered in terms of its conditional gene expression R 2 between nodes, while the indirect evidence is interpreted as an estimated probability. Therefore, we use a non-parametric approach to combine evidence. Specifically, the targets of each transcription factor are ranked and combined as a weighted sum, Our choice of the weight is by default α = 0.5, corresponding to an equal contribution of direct and indirect evidence. This parameter could be adjusted if the context of a study involved reason to prefer one source of evidence over the other (see Supporting Information).
Applying this approach to gene expression data from two distinct phenotypes results in two p × m gene regulatory adjacency matrices, one for each phenotype. These matrices represent estimates of the targeting patterns of the m transcription factors onto the p genes. This network inference algorithm finds validated regulatory interactions in Escherichia coli and Yeast (Saccharomyces cerevisiae) data sets (see Supporting Information).
Modeling the state transition matrix: Many methods have been developed for inferring gene regulatory networks, but more recent work has been proposed for estimating gene regulatory network differentiation [14,15]. Once we have gene regulatory network estimates for each phenotype, we can model the problem of estimating the transition matrix within a regression framework. With this formulation, we solve for the m × m matrix that best describes the transformation between phenotypes (1). More specifically, MONSTER predicts the change in edge-weights for a transcription factor, indexed i, in a network based on all of the edge-weights in the baseline phenotype network.
where b i and a i are column-vectors in B and A that describe the regulatory targeting of transcription factor i in the final and initial networks, respectively.
In the simplest case, this can be solved with normal equations, to generate each of the columns of the transition matrix T such that The regression is performed m times corresponding to each of the transcription factors in the data. In this sense, columns in the transition matrix can be loosely interpreted as the optimal linear combination of columns in the initial state adjacency matrix which predict the column in the final state adjacency matrix. The interpretation of the transition matrix can be best understood by comparing it to the identity matrix. A transcription factor, i, that does not alter its regulatory targets between states will have expected values of 0 for all entries in column i, with the exception of entry i. In the context of discovering changes in network configurations for a transcription factor, we are most interested in evaluating the degree to which each column has non-zero values for all non-i th entries. In essence, we are describing one network as a linear combination of another network. Numerous biological mechanisms, such as the formation of protein complexes, protein inactivation, post-translational modification, epigenetics, etc. allow for the systematic modification of network structures and drive the changes that are detected in the transition matrix (see Supporting Information).
This framework allows for the natural extension of constraints such as L1 and/or L2 regularization (see Supporting Information). For the analysis we present in this manuscript, we use the normal equations and do not impose a penalty on the regression coefficients. The inclusion or exclusion of this feature should depend primarily on assumptions of the underlying network mechanisms. L1 regularization, for example, will tend to infer a sparse transition matrix. This is reasonable in contexts in which the regulatory targeting pattern of TFs is not expected to change for the vast majority of TFs. For our four observational studies of COPD, a highly complex disease, it is not reasonable to assume that transcription factors differ in a sparse manner.
Computing the transcription factor involvement: For a transition between two nearly identical states, we expect that the transition matrix would approximate the identity matrix. However, as initial and final states diverge, there should be increasing differences in their corresponding gene regulatory networks and, consequently, the transition matrix will also increasingly diverge from the identity matrix. In this model, the transcription factors that most significantly alter their regulatory targets will have the greatest "off-diagonal mass" in the transition matrix, meaning that they will have very different targets between states and so are likely to be involved in the state transition process. We define the "differential transcription factor involvement" (dTFI) as the magnitude of the off-diagonal mass associated with each transcription factor, or, where,τ i,j is the value of the element in the i th row and j th column in the transition matrix, corresponding to the i th and j th transcription factors. To estimate the significance of this statistic, we randomly permute sample labels n = 400 times across phenotypes. From these 400 permuted results, we infer the standard error for each estimate under the null. P-values are determined by comparing the observed estimate to this standard error and FDR values is computed from these p-values (see Supporting Information).

MONSTER finds significantly differentially involved transcription factors in COPD with strong concordance in independent data sets
As a demonstration of the power of MONSTER to identify driving factors in disease, we applied the method to casecontrol gene expression data sets from four independent Chronic Obstructive Pulmonary Disease (COPD) cohorts: Evaluation of COPD Longitudinally to Identify Predictive Surrogate Endpoints (ECLIPSE) [16,17], COPDGene [18][19][20], Lung Genomics Research Consortium (LGRC) [21] and Lung Tissue from Channing Division of Network Medicine (LT-CDNM) [22]. The tissues assayed in ECLIPSE and COPDGene were whole blood and peripheral blood mononuclear cells (PBMCs), respectively, while homogenized lung tissue was sampled for LGRC and LT-CDNM. As a baseline comparison metric, we evaluated the efficacy of applying commonly used network inference methods on these case-control studies. In analyzing phenotypic changes, networks are generally compared directly, with changes in the presence or weight of edges between key genes being of primary interest. It is therefore reasonable to assume that any reliable network results generated from a comparison of disease to controls will be reproducible in independent studies. We investigated whether this is the case for our four COPD data sets using three widely used network inference methods -Algorithm for the Reconstruction of Gene Regulatory Networks (ARACNE) [23], Context Likelihood of Relatedness (CLR) [24], and Weighted Gene Correlation Network Analysis (WGCNA) [25] -computing the difference in edge weights between cases and controls for each of the four studies. We found no meaningful correlation (R 2 < .01) of edge weight difference across any of the studies regardless of network inference method or tissue type (Additional file 1: Supporting Figure S3). Edge weight differences, even when very large in one study, did not reproduce in other studies. This suggests that a simple direct comparison of edges between inferred networks is insufficient for extracting reproducible drivers of network state transitions. This finding may be unsurprising given the difficulty in inferring individual edges in the presence of heterogeneous phenotypic states, technical and biological noise with a limited number of samples.
The lack of replication in edge-weight differences between independent data sets representing similar study designs indicates that we need to rethink how we evaluate network state transitions. MONSTER provides a unique approach for making that comparison. In each of the four COPD data sets, we used MONSTER to calculate the differential transcription factor involvement (dTFI, Eq. 2) for each transcription factor and used permutation analysis to estimate their significance (Fig. 2, Additional file 1: Additional Figures S1-S3). We observed strongly significant (p < 1e − 15) correlation in dTFI values for each pairwise combination of studies. In addition, out of the top 10 most differentially involved transcription factors in the ECLIPSE and COPDGene studies, we found 7 to be in common. Furthermore, three of these seven transcription factors (GABPA, ELK4, ELK1) also appeared as significant in the LGRC results with FDR<0.01 and each of the top five ECLIPSE results were among the top seven in the LT-CDNM results (Additional file 1: Additional Table S1, Additional file 1: Additional Figure S3). This agreement is quite striking considering that the there was almost no correlation in the edge-weight differences across these same studies when we tested the other methods. But it is exactly what we should expect-that the same method applied to independent studies of the same phenotypes should produce largely consistent results.
Many of the top dTFI transcription factors, especially those identified by MONSTER across all four studies, are biologically plausible candidates to be involved in the etiology of COPD (Additional file 1: Additional Table S1, Additional file 1: Additional Figures S1-S3). For example, E2F4 is a transcriptional repressor important in airway development [26] and studies have begun to demonstrate the relevance of developmental pathways in COPD pathogenesis [27].
Some of the greatest effect sizes across all four studies were found for SP1 and SP2. An additional member of the SP transcription factor family, SP3, has been shown to regulate HHIP, a known COPD susceptibility gene [28]. Both SP1 and SP2 form complexes with the E2F family [29,30] and may play a key role in the alteration of E2F4 targeting behavior. Furthermore, E2F4 has been found to form a complex with EGR-1 (a highly significant transcription factor in ECLIPSE and LT-CDNM) in response smoke exposure, which may lead to autophagy, apoptosis and subsequently to development of emphysema [31].
Mitochondrial mechanisms have also been associated with COPD progression [32]. Two of most highly significant transcription factors based on dTFI in ECLIPSE were NRF1 and GABPA (FDR<.001). Indeed, these TFs had highly significant dTFI (FDR<0.1) in all four studies. NRF1 regulates the expression of nuclear encoded mitochondrial proteins [33]. GABPA, also known as human nuclear Fig. 2 MONSTER analysis results in the ECLIPSE study. a Heatmap depicting the transition matrix calculated for smoker controls "transitioning" to COPD by applying MONSTER to ECLIPSE gene expression data. For the purposes of visualization, the magnitude of the diagonal is set to zero. b A network visualization of the 100 largest transitions identified based on the transition matrix in (a). Arrows indicate a change in edges from a transcription factor in the Smoker-Control network to resemble those of a transcription factor in the COPD network. Edge thickness represents the magnitude of the transition and node (TFs) sizes represent the dTFI for that TF. Blue edges represent a gain of targeting features and red represents the loss. c The dTFI score from MONSTER (red) and the background null distribution of dTFI values (blue) as estimated by 400 random sample permutations of the data respiratory factor-2 subunit alpha, may have a similar role in nuclear control of mitochondrial gene expression. Furthermore, GABPA interacts with SP1 [34] providing evidence of a potentially shared regulatory mechanism with E2F4.
Overall, we found a strong correlation across studies in transcription factors identified as significantly differentially involved (Fig. 3a-b). It is reassuring that we find the strongest agreement when comparing studies that assayed similar tissues. However the fact that we see similar dTFI signal across studies involving different tissue types is also notable as it suggests that the transition from smoker control to disease phenotype affects multiple tissues and supports the growing evidence for a role in immune response in COPD pathogenesis.
Gene regulatory networks, and results derived from their comparison, are notoriously difficult to replicate across studies [35]. The four studies we used each has unique aspects, including the choice of microarray platform, study demographics, location, time, and tissue. Nevertheless, MONSTER identified similar sets of transcription factors associated with the transition between cases and controls. This consistency in biologically-relevant transcription factors, associated with the transition from the control phenotype to disease, in four independent studies suggests that MONSTER can Fig. 3 Strong reproducibility in top differential transcription factor involvement found in case-control COPD studies. ECLIPSE and COPDGene profiled gene expression in whole-blood and PBMC while the gene expression data in LGRC and LT-CDNM were assayed in lung tissue. a Results for studies with gene expression data obtained from the same-tissue. Both the blood based (left) and lung tissue studies (right) demonstrate very high Spearman correlation of differential involvement. b Despite using data from different sources we found agreement between studies of different tissues. c Venn diagram depicting the top 20 transcription factors found in each study. The union of all top 20 lists contains 36 transcription factors provide not only robust network models, but also can identify reliable differences between networks.
Despite the overall consistency, some transcription factors had variable dTFI across studies. For example, using the LGRC dataset, we discovered a highly significant (FDR < .0001) differential targeting pattern involving the transcription factors RFX1 and RFX2 (Additional file 1: Additional Table S1). However, these same TFs were not identified as potential drivers of the control to COPD transition in either the ECLIPSE or COPDGene study. This difference is likely due the differences in tissue type as the RFX family transcription factors are known to regulate ciliogenesis [36]. Cilia are critical for clearing mucous from the airways of healthy individuals, but disruption can lead to infection and potentially to chronic airflow obstruction [37][38][39].
The hypothesis behind MONSTER is that each phenotype has a unique gene regulatory network and that a change in phenotypic state is reflected in changes in transcription factor targeting. That hypothesis translates to an expectation that transcription factors driving change in phenotype will have the greatest dTFI scores. One might expect that these "driving transcription factors" would be also be differentially expressed. We compared dTFI to differential expression (ECLIPSE Fig. 4, other studies shown in Additional file 1: Additional Figure S4) and found that many of the transcription factors with high dTFI values were not differentially expressed. This suggests that there are other mechanisms, such as epigenetic modification of the genome or protein modifications, that alter the structure of the regulatory network by changing which genes are targeted by key transcription factors.

Discussion
One of the fundamental problems in biology is modeling the transition between biological states such as that which occurs during development or as a healthy tissue transforms into a disease state. As our ability to generate large-scale, integrative multi-omic data sets has grown, there has been an increased interest in using those data to infer gene regulatory networks to model fundamental biological processes. There have been many network inference methods published, each of which uses a different approach to estimating the "strength" of interactions between genes (or between transcription factors and their Fig. 4 Differentially involved transcription factors are not necessarily differentially expressed. A plot of the differential expression versus the differential involvement for transcription factors based on our analysis of the ECLIPSE data. MONSTER commonly finds transcription factors which are differentially involved but are expressed at similar levels across cases and controls. Importantly, these transcription factors would not have been identified using conventional differential expression methods. This demonstrates the unique potential MONSTER has for discovery beyond standard gene expression analysis targets). But all suffer from the same fundamental limitation: every method relies on estimating weights that represent the likelihood of an interaction between two genes to identify "real" (high confidence) edges. In comparing phenotypes, most methods then subtract discretized edges in one phenotype from those in the other to search for differences.
MONSTER represents a new way of looking at phenotypic transitions, but one that captures many aspects of what we should expect. First, we have to recognize that there is no single network that represents a phenotype, but that each phenotype is represented by a family of networks that all vary slightly from each other, yet which have essential features that are consistent with the phenotype. What this means is that each regulatory edge in a network representation has to be represented by continuous, rather than discrete, variables. This captures the fact that regulatory interactions are stronger in certain individuals and weaker in others, or present in some and absent in others, but that, on average, they represent a distribution.
Second, when we consider a change in phenotype, that will be reflected in altered patterns of gene expression, and ultimately in the networks that represent the phenotype. In a transition, some individuals will experience a greater change while others will experience a smaller change. But overall, regulatory patterns in the network will shift as the phenotype changes.
Third, the change in the gene regulatory network structure between phenotypes will be driven by changes in the connectivity of the regulators-the transcription factors that alter when, how, and how strongly genes are expressed. A natural hypothesis in this model is that the transition between phenotype is likely associated with the transcription factors that experience the greatest change in their regulatory patterns between states, and that the activation or inactivation of their target genes, and the functions carried out by those genes, likely reflect the phenotypic differences between states.
MONSTER captures these features, creating initial and final state network representations and estimating the change in transcription factor regulatory patterns by estimating a transition matrix. For each transcription factor, the "off diagonal mass" calculated as the differential transcription factor involvement (dTFI), identifies those transcription factors that are ultimately likely to drive the phenotypic state transition.
This approach has several limitations. For example, MONSTER does not attempt to infer specific regulatory mechanisms. Rather, by focusing on transcription factors, the goal of the method is to identify which transcription factors change their behavior between study groups. A deeper and more targeted investigation of the specific regulatory mechanisms that may be underlying these changes in transcription factor involvement would be needed to fully interpret the changes in targeting patterns that MONSTER is able to identify.
In applying MONSTER to four independent COPD gene expression data sets surveying both COPD and smoker controls, a highly consistent picture of the transcription factors associated with disease development emerges. This consistency is, to some, surprising as gene expression data is notoriously noisy, with each study finding sets of differentially expressed genes that often are not concordant. By focusing on transcriptional regulators, MONSTER seems to be able to separate a cleaner signal from the noise and one that makes some biological sense. Indeed, when one looks at the transcription factors found by MONSTER as associated with the transition, all are biologically plausible candidates which provide new and important opportunities for future molecular studies of COPD pathogenesis. It is also noteworthy that many of these transcription factors could not have been found through a simple differential expression analysis as their transcriptional levels do not change significantly between disease and control populations. Rather, it is the regulatory patterns of these transcription factors, possibly driven by epigenetic or other changes, that shifts with the phenotype.

Conclusion
The systems biology research community has long framed the discussion of the transitions between phenotypes in terms of gene expression state space, in which a change in phenotype corresponds to a transition between one characteristic expression profile to another. Here we extend that framework to the gene regulatory network state space, recognizing that patterns of gene expression are driven by alterations in patterns of gene regulation-and therefore, changes in the corresponding gene regulatory network. This shift allows us to ask how phenotypic alterations, and the corresponding regulatory changes, are effected through transcription factor "rewiring, " and to identify those transcription factors that are altering their regulatory targets. Thus MONSTER represents a novel approach to identifying and understanding the factors that drive changes in biological states and one with broad potential for application in a range of systems.


extfieldsofstudy: Chemistry Medicine

s2fieldsofstudy: Biology Medicine Chemistry

text:
Loss of the vitamin D receptor in human breast and prostate cancers strongly induces cell apoptosis through downregulation of Wnt/β-catenin signaling

Vitamin D co-regulates cell proliferation, differentiation and apoptosis in numerous tissues, including cancers. The known anti-proliferative and pro-apoptotic actions of the active metabolite of vitamin D, 1,25-dihydroxy-vitamin D [1,25(OH)2D] are mediated through binding to the vitamin D receptor (VDR). Here, we report on the unexpected finding that stable knockdown of VDR expression in the human breast and prostate cancer cell lines, MDA-MB-231 and PC3, strongly induces cell apoptosis and inhibits cell proliferation in vitro. Implantation of these VDR knockdown cells into the mammary fat pad (MDA-MB-231), subcutaneously (PC3) or intra-tibially (both cell lines) in immune-incompetent nude mice resulted in reduced tumor growth associated with increased apoptosis and reduced cell proliferation compared with controls. These growth-retarding effects of VDR knockdown occur in the presence and absence of vitamin D and are independent of whether cells were grown in bone or soft tissues. Transcriptome analysis of VDR knockdown and non-target control cell lines demonstrated that loss of the VDR was associated with significant attenuation in the Wnt/β-catenin signaling pathway. In particular, cytoplasmic and nuclear β-catenin protein levels were reduced with a corresponding downregulation of downstream genes such as Axin2, Cyclin D1, interleukin-6 (IL-6), and IL-8. Stabilization of β-catenin using the GSK-3β inhibitor BIO partly reversed the growth-retarding effects of VDR knockdown. Our results indicate that the unliganded VDR possesses hitherto unknown functions to promote breast and prostate cancer growth, which appear to be operational not only within but also outside the bone environment. These novel functions contrast with the known anti-proliferative nuclear actions of the liganded VDR and may represent targets for new diagnostic and therapeutic approaches in breast and prostate cancer.

INTRODUCTION
Breast and prostate cancers are among the most prevalent malignancies in industrialized countries. Although mortality has steadily declined over the past 20 years, a significant proportion of patients eventually develop metastatic disease, most frequently to the skeleton. 1 Skeletal related events due to bone metastasis are frequent and a major cause of mortality and morbidity. [2][3] We have previously shown in rodent models that reduced bone turnover inhibits, while increased bone turnover accelerates, breast and prostate cancer growth in bone. [4][5][6][7][8][9][10] These experimental findings provide a logical explanation for the clinical observation that accelerated bone turnover is associated with higher rates of skeletal related events and poorer prognosis in patients with breast or prostate cancer. 11 They also offer a rationale for the use 1 of anti-resorptive agents in these patients. 12 We have further reported that in rodent models, vitamin D deficiency promotes the growth of breast and prostate cancer cells in bone. 5,8,10,13 These effects appear to be mediated mainly through an increase in bone remodeling, that is, indirectly through changes in the bone microenvironment. However, inhibition of bone remodeling with potent anti-resorptive treatments (for example, osteoprotegerin, zoledronic acid) fails to completely reverse the pro-proliferative effects of vitamin D deficiency, 8,10 suggesting that vitamin D deficiency may also promote cancer cell growth by an additional and possibly direct mechanism.
Apart from its function in calcium and phosphate homeostasis, vitamin D is known to exert strong antiproliferative, pro-differentiation and pro-apoptotic actions in numerous cell types and tissues, including cancers. 14 The biologically active metabolite of vitamin D, 1,25-dihydroxyvitamin D [1,25(OH) 2 D], acts through binding to the vitamin D receptor (VDR), a member of the nuclear steroid hormone receptor superfamily. In the absence of vitamin D, the VDR remains in the cytoplasm. Ligand binding causes the VDR to form a heterodimer with the retinoid X receptor, which facilitates movement of the VDR-ligand complex out of the cytoplasm and into the nucleus. Within the nucleus this complex binds to vitamin D-responsive elements in the regulatory region of target genes. 15 In keeping with the classical function of vitamin D in regulating calcium and phosphate homeostasis, the VDR is expressed at high levels in tissues such as the intestine, bone and kidney. 15 However, research into the pleiotropic actions of vitamin D has revealed that the VDR is also expressed in numerous other tissues, including malignant tumors. 14 Limited clinical studies in patients with breast and prostate cancer demonstrated that VDR expression in these tumors is negatively associated with tumor size and lymph node involvement. [16][17][18] Furthermore, mice with global VDR knock-out show increased sensitivity to carcinogen challenges. 14,[19][20] These and the findings of our previous studies 8, 10 point to a broader role of the VDR in the regulation of cell growth, which may go beyond its classical function as a ligand-specific nuclear receptor. Using a VDR knockdown approach, the current in vitro and in vivo studies aimed to further define the role of the VDR in the regulation of breast and prostate cancer growth.

Cell culture
The human breast cancer cell line, MDA-MB-231, 21 was obtained from ATCC. Cells were cultured in Dulbecco's Modified Eagle's medium, supplemented with 10% fetal bovine serum (FBS, JRH Biosciences, KS, USA) and 1% penicillin-streptomycin. The human prostate cancer cell line, PC3, was obtained from ATCC. Cells were cultured in RPMI, supplemented with 10% fetal bovine serum (FBS, JRH Biosciences, Lenexa, Kansas, USA) and 1% penicillinstreptomycin solution. Unless otherwise stated, tissue culture media and supplements were from Life Technologies (Carlsbad, CA, USA).
Knockdown of VDR expression in cancer cells VDR expression was silenced in MDA-MB-231 and PC3 cells via a lentiviral-based expression system driving the production of short hairpin RNA species (shRNAs, Sigma, St. Louis, MO, USA). The clones selected were VDR (referred as "VDR-KD", transfected with shRNA, TRCN000019506, Sigma) and non-target control (referred as "NT", transfected with non-target RNA, SHC002V, Sigma). VDR-KD and NT cells from both MDA-MB-231 and PC3 cells were selected using 2 μg·mL − 1 puromycin (Sigma, USA) for two weeks. The cells resulting from this transfection will be referred to as MDA-NT, MDA-VDR-KD, PC3-NT, and PC3-VDR-KD.
Cell growth assays NT and VDR-KD cells were treated with 10 -8 mol · L −1 1,25dihydroxy vitamin D 3 [1,25(OH) 2 D 3 ] or vehicle (0.1% ethanol) in medium supplemented with 2% charcoal stripped heat-inactivated fetal bovine serum (FBS). In our hands, untreated FBS contains 1,25(OH) 2 D at a concentration of 360 pmol·L − 1 when measured by RIA (Diasorin, Stillwater, MN, USA). In contrast, 1,25(OH) 2 D is undetectable in media supplemented with 2% charcoal stripped heat-inactivated FBS. We therefore consider culture media supplemented with 2% charcoal stripped heat-inactivated serum to be free of 1,25(OH) 2 D ("ligand-free media"). Cell growth assays were performed as described previously. 6,22 Briefly, cells were seeded in media containing 2% charcoal stripped FBS and treated with either 10 − 8 mol · L −1 , 1,25 (OH) 2 D 3 or vehicle. For the rescue experiments, cells were treated in the same manner with 1 μmol · L −1 6-bromoindirubin-3'-oxime (BIO, Sigma), a compound that binds to, and inactivates the enzyme that normally phosphorylates and degrades β-catenin (that is, GSK-3β). Addition of BIO to the cell culture inhibits this process and therefore stabilizes β-catenin signaling, 23  Western blot analysis For protein expression, cells were harvested at 24 h post plating. Whole cell lysates were prepared using RIPA buffer and cell fractions were prepared using NE-PER Nuclear and Cytoplasmic Extraction kit (Thermo Scientific, Rockford, IL, USA), supplemented with protease inhibitor cocktail (Roche Pharmaceutical). Protein concentrations were determined using a Pierce 660 nm protein assay (Thermo Scientific), proteins were separated on Bolt 4%-12% mini gels and then transferred onto nitrocellulose membrane by means of the iBlot2 rapid transfer (Life Technologies). Membranes were blocked in TBS-T 5% skim milk for 1 h, incubated overnight with primary antibody and then with secondary antibody for 1 h. Proteins were visualized with Amersham ECL Western Blotting Detection Reagents (GE Healthcare, Amersham, Buckinghamshire, UK) and detected with the Bio-Rad Image Lab system (Bio-Rad, Hercules, California, USA). Primary antibodies used were: anti-VDR (Clone 9A7γ.E10.E4, Thermo Scientific, Rockford, IL, USA), anti-β-catenin, anti-β-actin, (Cell Signaling Technology, Danvers, MA, USA), GSK-3β (BD Biosciences, Franklin Lakes, NJ, USA), and anti-caspase 3 (Merck Millipore, Frenchs Forest, NSW, Australia).
Animal models Five-week-old female and male BALB/c nu/nu mice were purchased from the Animal Resources Centre, Canning Vale, WA, Australia. Previously studies in our laboratory confirmed these mice to be fully replete for 25(OH)D levels. 8,10 Murine xenograft models were employed to male mice) cancer cell growth in bone and, orthotopically, in either the mammary fat pads (female mice) or the subcutaneous tissue of the inferior-lateral flank (male mice), as described previously. 6,22 In accordance with Animal Welfare Guidelines and approved protocols, all mouse manipulations were performed inside a laminar-flow hood under aseptic conditions whilst maintaining general anesthesia following intra-peritoneal injection of ketamine/ xylazine (75/10 mg·kg − 1 ). The analgesic carprofen (Rimadyl) (5 mg·kg − 1 subcutaneously) was administered before inoculation to minimize post-surgical pain.
For the orthotopic and subcutaneous models, cells were prepared as above and suspended in a mixture of cold Matrigel and PBS (1:1) at a concentration of 1 × 10 7 cells per mL. 100 μL of the suspension was either implanted orthotopically into the mammary fad pad or subcutaneously at the flank (n = 8-10 per group). Tumor growth was measured by calipers every second day from day 6 onwards and tumor volume was determined using the modified ellipsoidal formula [volume = (length × width 2 )/2]. For the intratibial xenograft model, 10 μL of cell suspension (for MDA, 5 × 10 4 cells per injection, for PC3 cells 10 5 cells per injection) were slowly injected through the tibial plateau of the left proximal tibia using a Hamilton syringe. 9 Contralateral tibiae were injected similarly with vehicle alone (PBS). For the monitoring of lytic bone lesions in vivo, mice were anesthetized as above and assessed by digital radiography: for MDA cells, images were taken on days 10, 17 and 21 post cell injection; for PC3 cells, images were taken on days 17, 24 and 31 post cell injection. Lytic bone areas within the tibiae were measured using interactive image analysis software (ImageJ, NIH, Bethesda, Maryland, USA) and the size of osteolytic lesions was calculated. The growth of lytic lesions was compared among experimental groups. The intratibial experiments were terminated before mice showed adverse effects of the tumor growth.
Histological staining and immunohistochemistry Tibiae were fixed for 36 h in 4% paraformaldehyde buffered with 0.1 mol · L −1 phosphate buffer (pH 7.4) and decalcified in 10% EDTA at 4°C for 2 weeks. The tissues were then processed and embedded in paraffin. Four-micron sections were cut from each specimen and stained with haematoxylin and eosin for routine histological examination. Osteoclasts were detected via histochemical analysis of tartrate-resistant acid phosphatase (TRAcP), using naphthol AS-MX phosphate as substrate and fast red violet Luria-Bertani salt (both from Sigma Chemical, St. Louis, MA, USA) as a stain for the reaction product; incubation was performed at room temperature for 30 min.
Micro-computed tomography imaging After tissue harvest, representative micro-computed tomography (CT) images of tibiae were obtained using a SkyScan 1172 scanner (SkyScan). Scanning was done at 100 kV, 100 μA using a 1-mm aluminum filter. In total, 1 800 projections were collected at a resolution of 6.93 μm per pixel. Reconstruction of sections was done using a modified Feldkamp cone-beam algorithm with beam hardening correction set to 50%. VGStudio MAX software (Volume Graphics GmbH) was used to obtain three-dimensional visualization of tibiae from reconstructed sections.
Bone histomorphometry Histomorphometric analysis of the proximal tibial metaphysis was conducted to evaluate bone volume, tumor burden and static measures of bone resorption. Measurements were performed in longitudinal 4 μm thick sections stained with haematoxylin and eosin (×12.5 magnification) using the OsteoMeasure System (Osteometrics, Atlanta, GA, USA). Multiple sagittal sections were cut through each tibial sample and stained. The sagittal section with largest tumor area was selected and the tumor area was measured and used as an index of tumor burden. Cortical bone area was measured in the same sections. Osteoclasts were identified as TRAcP-positive multinucleated cells at the tumor-bone interface. The number of osteoclasts per millimeter of tumor-bone interface was calculated at a magnification of × 200.
Statistical analysis All data were presented as the mean ± SE. Statistical analyses were conducted using one-way ANOVA followed   (Figure 1g). Surprisingly, however, when cultured under ligand-free conditions the growth of MDA-VDR-KD cells was also significantly reduced, with growth curves similar to those of 1,25(OH) 2 D 3 treated MDA-NT cells. As expected, treatment of VDR knockdown cells with 1,25(OH) 2 D 3 had no further effect on growth (Figure 1d). Compared to NT controls, cell proliferation was significantly reduced in MDA-VDR-KD cells, again of similar magnitude to that seen in 1,25(OH) 2 D 3 -treated MDA-NT cells (Figure 1e and f). In contrast, VDR knockdown was associated with a pronounced rise in cell apoptosis: Thus, the proportion of TUNEL-positive cells was 12-fold higher in MDA-VDR-KD cells compared to MDA-NT controls (Figure 1g). In keeping with these observations, Caspase 3 mRNA expression and cleaved Caspase 3 protein levels were significantly higher in MDA-VDR-KD than NT cells (Figure 1h   (d-j) Intra-tibial implantation: At endpoint (day 21 post implantation) lytic lesion size on X-ray, micro-CT (d,e) and histological tumor area (f) were significantly smaller in mice implanted with MDA-VDR-KD than with MDA-NT cells (n = 15). Compared to NT controls, tumors derived from VDR knockdown cells were characterized by an increased proportion of apoptotic cancer cells (g) and lower mitotic activity (h). At the bone/tumor interface, tumors derived from VDR knockdown cells exhibited increased total and cortical bone area (i,j), and reduced osteoclast number (k). Asterisks denote significant difference from MDA-NT (*P o0.05; **P o0.01). Data are mean ± s.e.m.
Bone Research (2017) 17023 VDR controls cancer cell apoptosis and proliferation Y Zheng et al into the mammary fat pad or by inoculation into the left tibia. 22 Following orthotopic implantation of MDA-MB-321 cells into the mammary soft tissue, VDR ablation was associated with significantly reduced tumor growth compared to VDR-expressing NT controls (Figure 2a). At endpoint (day 33 post implantation (p.i.)), tumor mass was 40% lower and apoptosis was 36% higher in tumors derived from MDA-VDR-KD cells compared to those grown from MDA-NT cells (n = 10 per group, Po0.05; Figure 2b and c).
When MDA-NT and MDA-VDR-KD cells were implanted into the tibiae, tumors derived from VDR knockdown cells produced significantly smaller lytic lesions (Figure 2d and e) and tumor areas (Figure 2f) at endpoint (day 21 p.i.) than tumors generated by NT cells (n = 15, Po0.05). In keeping with these phenotypes, tumors derived from VDR knockdown cells contained a higher proportion of apoptotic cells (+80%) and demonstrated lower mitotic activity (−39%) (Figure 2g and h). Implantation of MDA-VDR-KD cells also resulted in significantly reduced bone destruction (as demonstrated by significantly greater total and cortical bone area compared to MDA-NT; Figure 2i and j) and reduced osteoclast number at the bone/tumor interface (Figure 2k). These results indicate that knockdown of VDR in breast cancer cells changes tumor growth characteristics in vivo regardless of the tissue environment.

VDR knockdown in PC3 cells-in vitro and in vivo studies
In order to investigate whether the effects of VDR ablation on cancer cell behavior are reproducible across different types of cancer, we next studied the human prostate cancer cell line, PC3, in analogous experiments. Compared to non-target controls (PC3-NT, cells transfected with non-target RNA), VDR mRNA and protein expression was reduced by 80% in PC3 cells at 24 h post plating (PC3-VDR-KD) (Figure 3a and b). Functionally effective ablation of the VDR was again confirmed by a significant reduction in CYP24 expression following stimulation with 10 − 8 mol · L −1 1,25(OH) 2 D 3 (Figure 3c). The in vitro growth curves of PC3 cells exhibited similar results to those seen with MDA-MB-231 cells: Knockdown of VDR expression significantly reduced PC3 cell growth, while at the same time stimulating cancer cell apoptosis. These changes were independent of the effects of 1,25(OH) 2 D 3 and of similar magnitude to those  (Figure 5a). Pathway enrichment analysis demonstrated that genes within the Wnt/β-catenin pathway were all downregulated in VDR knockdown cells (Figure 5b). Expression of these genes was further validated by real-time RT-PCR, establishing significant downregulation of relevant genes most of which are downstream of β-catenin (Figure 5c). Western blot analysis indicated that both cytoplasmic and nuclear β-catenin protein levels were reduced in MDA-VDR-KD compared to MDA-NT cells. Treatment with 1,25(OH) 2 D 3 induced nuclear β-catenin protein levels in MDA-NT but not in MDA-VDR-KD cells (Figure 5d). This reduced β-catenin protein levels could was reduced by 40% and the proportion of apoptotic cells (c) was increased by 125% in PC3-VDR-KD compared to PC3-NT tumors (n = 10). Asterisks denote significant difference from controls (*P o0.05). Data are mean ± s.e.m. (d-j) Intra-tibial implantation: Lytic lesion size (d,e) and tumor area (f) at endpoint (day 31 p. i.) were significantly smaller in mice implanted with PC3-VDR-KD cells compared to those implanted with PC3-NT cells (n = 12). Compared to NT controls, tumors derived from PC3-VDR-KD cells were characterized by increased apoptosis (g) and lower mitotic activity (h), and at the bone/tumor interface exhibited increased total and cortical bone area (i,j) and reduced osteoclast number (k). Asterisks denote significance difference from controls (*P o 0.05; **Po 0.01). Data are mean ± s.e.m.
Bone Research (2017) 17023 VDR controls cancer cell apoptosis and proliferation Y Zheng et al 8 be due to, at least in part, the reductions of Wnt3a and FZDR, the upstream initiative genes (Figure 5c) but were not due to GSK-3β as its expression is not altered in MDA-VDR-KD cells (Figure 5d).
Ex vivo analyses of tumor tissue confirmed that β-catenin protein levels were reduced in tumors derived from VDR knockdown cells relative to NT controls (Figure 5e). There is solid evidence to indicate that aberrant activation of Wnt/β-catenin signaling inhibits cell apoptosis while promoting cell proliferation. [26][27][28] Thus, the slower growth of tumors generated by VDR knockdown cells may be due to a reduction in Wnt/β-catenin activation. To test this hypothesis, we treated the cells with the GSK-3β inhibitor, 6-bromoindirubin-3'-oxime (BIO, Sigma), to stabilize and therefore enhance β-catenin levels. 23 At concentrations of 1 and 10 μmol · L −1 , BIO elevated cytoplasmic β-catenin and increased nuclear β-catenin protein levels in both MDA-NT and MDA-VDR-KD cells (Figure 5f). Consequently, Axin2 mRNA levels were increased 6-and 3-fold in BIO-treated NT and VDR-KD cells, respectively (Figure 5g). In contrast, BIO had no effect on mRNA expression of Cyclin D1 or Ki67 (Figure 5g). Growth curve analyses showed that treatment of VDR knockdown cells with BIO restored cell growth to approximately 70% of NT controls (Figure 5h). Consistent with the results for mRNA expression (Figure 5g), immunohistochemistry for Ki67 confirmed that treatment of NT or VDR-KD cells with BIO had no effect on proliferation (Figure 5j). In contrast, treatment of MDA-VDR-KD cells with BIO resulted in a significant reduction in the proportion of apoptotic cells, similar to the levels observed in NT cells (Figure 5i).

DISCUSSION
Our in vitro and in vivo studies demonstrate that the VDR has a function in promoting the growth of breast and prostate cancer cells, independent of its cognate ligand, 1,25(OH) 2 D. This previously unrecognized function of the VDR is operational both within and outside the bone microenvironment. Mechanistically, we established that knockdown of the VDR in cancer cells interferes with the Wnt/β-catenin signaling pathway. Thus, VDR knockdown is associated with a robust reduction in cytoplasmic βcatenin protein levels, resulting in enhanced cell apoptosis and reduced proliferation rates. Of note, stabilization of β-catenin and hence re-activation of β-catenin signaling in VDR knockdown cells partly restored cell growth, mostly through inhibition of cell apoptosis.
We previously reported that vitamin D deficiency promotes the growth of human breast and prostate cancer cells in bone. 5,8,10,13 This acceleration in cancer cell growth was, at least in part, attributable to the effects of vitamin D deficiency on the bone microenvironment, as hypovitaminosis D-via an increase in parathyroid hormone secretion-induces a significant increase in osteoclast-mediated bone resorption. 5,8,10,13 In addition, treatment with 1,25(OH) 2 D 3 inhibits breast and prostate cancer cell growth in vitro and in vivo. 8,10 This and the fact that inhibition of bone remodeling failed to completely reverse the pro-proliferative effects of vitamin D deficiency provided evidence that vitamin D signaling may be directly involved in the control of tumor growth in bone. Disrupting vitamin D signaling in cancer cells would therefore be expected to promote cell proliferation and tumor growth, similar to the effects of vitamin D deficiency.
However, contrary to this assumption and our expect ations, we found that silencing of the VDR in breast and prostate cancer cell lines consistently resulted in reduced cell growth and increased apoptosis in vitro. Importantly, these effects were seen not only in the presence but also in the absence of 1,25(OH) 2 D 3 . Furthermore, the effects of VDR knockdown on tumor growth, cancer cell proliferation and apoptosis observed in vitro were fully reproducible in vivo, evident not only in bone but also in soft tissues, and independent of whether human breast or prostate cancer cell lines were used. The fact that the VDR knockdown cells are insensitive to treatment with 1,25(OH) 2 D 3 is particularly relevant to the in vivo experiments, as these were performed in vitamin D-replete mice, thus providing evidence that the reduction in tumor growth seen in vivo with VDR-KD cells must be independent of vitamin D. Taken together, our findings indirectly suggest that the VDR holds hitherto unknown functions that promote cancer cell growth independent of its ligand.
Most of the known biological effects of 1,25(OH) 2 D are mediated through binding to the VDR. Ligand binding causes the VDR to form a heterodimer with the retinoid X receptor and to translocate from the cytoplasm to the nucleus, where the complex binds to vitamin D-responsive elements in the regulatory region of its target genes. 15 In the absence of vitamin D, the VDR is mostly retained in the cytoplasm, 29-30 however, recent evidence suggests that shuttling of the unliganded VDR to the nucleus produces basal transcriptional activity via binding to co-repressor or co-activators. [30][31][32][33][34][35] Global chromatin immunoprecipitation (ChIP) analyses have also indicated that a significant number of DNA-binding sites are occupied by the VDR in the absence of its cognate ligand. 36 There is little data available regarding a potential unliganded function of the VDR. Previous studies in skin suggest that the unliganded VDR plays a role in the context of normal keratinocyte stem cell function (through interaction with the canonical Wnt pathway). 20,[37][38] The fact that ablation of VDR expression in cancer cells results in a marked inhibition of cell and tumor growth raises the possibility of the unliganded VDR being involved in the regulation of cancer cell proliferation and apoptosis. However we acknowledge that this conclusion cannot be seen as definitive due to incomplete VDR ablation. Our data indicate that the proapoptotic and anti-proliferative effects of VDR knockdown are mediated through inhibition of the Wnt/β-catenin signaling pathway. Indeed, stabilizing β-catenin in VDR-KD cells partially restored the growth of these cancer cells. In this context, it was interesting to note that in vitro, BIO-induced stabilization of β-catenin reversed the effect of VDR-KD on cell apoptosis but not on cell proliferation (for example, Ki67, cyclin D1 expression). Also, compared to NT controls, the expression of interleukin-6 (IL-6) and IL-8 mRNA was significantly reduced in breast cancer cells with disrupted VDR signaling. We have previously shown that IL-6 mediates direct paracrine-autocrine signaling between cells of the osteoblast lineage and cancer cells, thereby enhancing the growth of breast cancers in bone. 22 The reduced growth of VDR knockdown cells and tumors derived from such cells may, at least in part, be due to their diminished expression and secretion of interleukins like IL-6.
Bone Research (2017) 17023 VDR controls cancer cell apoptosis and proliferation Y Zheng et al Wnt/β-catenin signaling plays a crucial role in breast and prostate cancer cell growth. [26][27][28] Aberrant activation of the Wnt/β-catenin signaling pathway in cancer cells is associated with the modification of pivotal processes such as Wnt expression, β-catenin stabilization or degradation, and nuclear translocation of β-catenin. In the nucleus, β-catenin complexes with T-cell factor and facilitates transcription of target genes activating cell proliferation and inhibiting apoptosis. Currently, it remains to be seen how the VDR interacts with β-catenin in cancer cells. For example, we did not detect changes in GSK-3β protein levels, suggesting that the low β-catenin levels observed in VDR knockdown cells may not be due to its degradation. Instead, in cells with disrupted VDR signaling we found significant changes in the expression of genes upstream of β-catenin (for example, Wnt3a, FZD4), although it appears that these changes would not fully explain the profound reduction in β-catenin concentrations and signaling.
Clinical research into the role of vitamin D and VDR status in breast and prostate cancer has remained controversial. While it appears that the VDR is present in differentiated breast cancers, [39][40] the associations between VDR expression, its downstream signaling and patient prognosis or survival still require further clarification. 16,18,[41][42] VDR expression levels may be associated with the transition from benign to invasive breast tumor. 40,[43][44] Some studies showed cancer progression in breast tumors was characterized by reduced CYP27B1 and increased CYP24A1 expression, indicating that during the process of malignant transformation the reduction in VDR expression and higher CYP24A1 levels renders the tumor less sensitive to vitamin D. 40,43 Our results may help to explain some of the discrepancies between vitamin D status and clinical outcome [45][46] as tumor behaviour and patient outcome may likely depend on both vitamin D and VDR status. Thus, in the case of vitamin D deficiency, cancer cells may lose the ligand-mediated nuclear signaling effects of the VDR and as a consequence, the unliganded VDR may become the dominant functional form. Hence, in a ligand-depleted context (that is, vitamin D deficiency), the unliganded VDR may be responsible for increased breast cancer cell growth. In the presence of adequate vitamin D, increased nuclear VDR localization and reduced unliganded VDR actions could together reduce breast cancer proliferation. It is possible that local production of 1,25(OH) 2 D, dependent on adequate serum 25(OH)D levels and tissue CYP27B1 activity, may be more important than levels of serum 1,25(OH) 2 D which do not correlate closely with vitamin D status. 47 These possibilities need to be further investigated by measuring nuclear and unliganded VDR expression in human tumors at different stages of progression and correlating these results with the patient's vitamin D status in a prospective clinical trial. 46 In conclusion, the current study demonstrates that the unliganded VDR possesses intrinsic functions that control cancer cell growth in a ligand-independent manner through interference with the Wnt/β-catenin signaling pathway. Targeting the unliganded VDR may open up new avenues for the development of novel diagnostic and therapeutic approaches in breast and prostate cancer.