Analysis of the Transcriptome: Regulation of Cancer Stemness in Breast Ductal Carcinoma In Situ by Vitamin D Compounds

Breast cancer is the most common cancer and the second leading cause of cancer-related deaths in women worldwide (1). On the basis of the presence or absence of estrogen receptor (ER), progesterone receptor (PR), and HER2, breast cancers are divided into subtypes: luminal A (ER⁺ and/or PR⁺; HER2⁻), luminal B (ER⁺ and/or PR⁺; HER2⁺), basal-like (ER⁻, PR⁻, and HER2⁻), and HER2-enriched (ER⁻, PR⁻, and HER2⁺; ref. 2). Breast cancer development is a multi-step process that involves epigenetic and genetic changes contributing to aberrant cell growth (3). Histologically, breast cancer can be staged into invasive ductal carcinoma (IDC), ductal carcinoma in situ (DCIS), and invasive lobular carcinoma. About 20% of breast cancers newly diagnosed in 2019 among U.S. women will be classified as DCIS, amounting to over 48,000 cases (4). DCIS is an early stage, noninvasive type characterized by proliferation of malignant epithelial cells in the ducts (5). It arises from atypical ductal hyperplasia and may progress to IDC and metastatic cancer (5). It is predicted that up to 50% of the DCIS cases will progress to IDC within 10 years of initial diagnosis (6). Gene expression and miRNA analyses have been performed to elucidate the molecular characteristics of DCIS progression to IDC (7, 8). However, the natural history of progression of DCIS to IDC is yet to be fully determined.

Cancer stem cells (CSC) were first identified in breast cancer using the cell surface markers CD44⁺/CD24⁻ (9). This population is characterized by a stem cell gene expression signatures, drug-resistant phenotype, and self-renewal capacity in vitro and in vivo (10). Human DCIS lesions form spheroids and duct-like structures in ex vivo organoid culture and tumors in immunodeficient mice, suggestive of the presence of CSC-like cells in these early tumors (11). MCF10DCIS.COM was derived from the nontumorigenic MCF10A human breast cell line and exhibits basal-like subtype properties (12). It is similar to human DCIS with bi-potentiality that can give rise to both myoepithelial and luminal cells and spontaneous progression to invasive breast cancer in vivo (13), and it is widely used as a model for IDC development from precursor lesions. MCF10DCIS.COM cultures and tumors contain high ALDH1⁺ and CD44⁺/CD49f⁺/CD24⁻ subpopulations with increased self-renewal and tumor development capabilities, similar to CSCs of fully invasive tumors (14).

Vitamin D signaling is known to be a potential target for breast cancer chemoprevention (15). Our laboratory has shown that vitamin D compounds inhibit triple-negative breast cancer tumorigenesis by reducing expression of CSC-associated genes, including OCT4 and CD44, and by inducing differentiation and upregulating myoepithelial markers (16). These compounds reduce in vitro mammosphere formation and in vivo tumorigenesis, although the molecular mechanisms of these effects are not known. BXL0124 is an analog of calcitriol (1α25(OH)₂D₃) modified with an additional side chain at C21-methyl group, endowing it with more biological activity at lower concentrations without causing hypercalcemia, a limiting side effect of vitamin D (17). BXL0124 also inhibits MCF10DCIS xenograft tumorigenesis more potently than 1α25(OH)₂D₃, apparently through the similar mechanism of suppressing CSCs (18). This study was undertaken to develop global profiles of changes in gene expression and CpG methylation in MCF10DCIS cells induced by vitamin D compounds, to gain an understanding of the pathways involved in their overall effects and the molecular basis of their activities, and to identify potential targets that could be exploited in the chemoprevention of breast cancer progression from DCIS to IDC.

1α25(OH)₂D₃ and a Gemini vitamin D analog [BXL0124; 1α,25-dihydroxy-20R-21(3-hydroxy-3-deuteromethyl-4,4,4-trideuterobutyl)-23-yne-26,27-hexafluoro-cholecalciferol, >95% purity] were provided by BioXCell, Inc. (17). Vitamin D compounds were dissolved in DMSO. MCF10DCIS.com human breast cancer cells (MCF10DCIS, RRID: CVCL_5552) were provided by Dr. Fred Miller at the Barbara Ann Karmanos Cancer Institute (Detroit, MI). Cells were maintained in DMEM/F12 medium supplemented with 5% horse serum, 1% penicillin/streptomycin, and 1% HEPES solution at 37°C, 5% CO₂. The cells were passaged every 3–4 days between passage number p30 and p50. Mycoplasma testing was done every 3 months using Mycoplasma PCR Detection Kit (catalog no. MP0035, Sigma-Aldrich). Cell line authentication was done on December 2019 using short tandem repeat profiling (catalog no. 135-XV) at ATCC.

MCF10DCIS cells were grown to 70%–80% confluence and cells were detached with StemPro Accutase (Life Technologies) cell detachment solution. Cells were then grown in 6-well ultra-low attachment plates at a density of 10,000 cells/mL and maintained in MammoCult Human Medium added with MammoCult proliferation supplement, hydrocortisone solution, and heparin solution (Stemcell Technologies). Spheres were treated with DMSO (0.01%), 1α25(OH)₂D₃ (100 nmol/L), and BXL0124 (10 nmol/L) for 5 days. After 5 days in culture, the culture plates were gently swirled to cluster the spheres in the middle of each well and photographed. Mammosphere forming efficiency (MFE) was calculated by dividing the number of mammospheres (≥100 μm) formed by the number of single cells seeded in individual wells. Three independent experiments were repeated.

MCF10DCIS mammospheres were treated with DMSO (0.01%), 1α25(OH)₂D₃ (100 nmol/L), and BXL0124 (10 nmol/L) for 5 days followed by extraction of RNA and DNA using AllPrep DNA/RNA Mini Kit (Qiagen). The concentrations and quality of the RNA and DNA were determined using NanoDrop spectrophotometer and Agilent 2100 Bioanalyzer separately. We prepared three independent sets of cultures with MCF10DCIS mammospheres treated with DMSO, 1α25(OH)₂D₃, and BXL0124, and extracted each separately to obtain three RNA samples and three DNA samples. We pooled those samples for the RNA-sequencing (RNA-seq) and methyl-sequencing (methyl-seq) analyses, which was carried out at RUCDR Infinite Biologics, a Rutgers University affiliated institution which provides next-generation sequencing (NGS) services. RNA-seq libraries were prepared using Illumina RNA Library Prep Kit v2 according to the manufacturer's user guide with 400 ng of RNA as input. The libraries were then quantified using KAPA Library Quantification Kit according to the manufacturer's user guide and pooled with barcodes. The pooled libraries were sequenced on Illumina NextSeq 550 System, using NextSeq 500/550 Mid Output v2 Kit. The sequencing parameters used were 150 bp, paired-end with around 20 million reads generated per sample. The DNA samples were further processed using an Agilent SureSelect Human Methyl-seq Target Enrichment System (Agilent Technologies) and sequenced on an Illumina NextSeq 500 instrument with 76 bp single-end reads, generating 34–47 million reads per sample.

The RNA-seq reads were mapped to human reference genome with Hisat2 software. Cufflinks was used to assemble the transcript products and calculate the fragment abundance. Cuffdiff was used to quantify transcripts of genes differentially expressed among the three treatment conditions (DMSO, 1α25(OH)₂D₃, and BXL0124). Methyl-seq was aligned using Bismark software (version 0.15.0), CpG sites were counted and clustered into differentially methylated regions (DMR) using DMRfinder (version 0.1; ref. 19). Genomic annotation was performed with ChIPseeker in R (20). R was also used for downstream NGS data analysis and visualization as we have reported previously (21, 22). The RNA-seq and methyl-seq datasets described in this study have been deposited in the NCBI Gene Expression Omnibus with accession number GSE148548.

Analysis of pathways and gene networks of the expression data was performed with Ingenuity Pathway Analysis (IPA) Software from Qiagen (Version 49309495).

To validate RNA-seq data, we performed three independent sets of experiments with MCF10DCIS mammospheres treatment with DMSO, 1α25(OH)₂D₃, and BXL0124. After 5 days of treatment, RNA samples are collected and quantitative PCR (qPCR) analysis was performed. RNA was reverse transcribed to cDNA using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems), followed by amplification of the cDNA by using primers of genes of interest and TaqMan 2X Universal PCR Master Mix (Applied Biosystems) and using the ViiA 7 Real-Time PCR System (Applied Biosystems). Primers used were GAPDH (Hs02758991), GDF15 (Hs00171132), LCN2 (Hs01008571), S100A4 (Hs00243202), NGFR (Hs00609976), PPP1R1B (Hs00259967), AGR2 (Hs00356521), KRT6A (Hs04194231), EMP1 (Hs00608055), IGFBP5 (Hs00181213), CAPN6 (Hs00560073), CYP24A1 (Hs00167999), and KRT5 (Hs00361185). GAPDH was used as an internal control. Relative changes of gene expression were calculated using ΔΔC_t method.

The bisulfite-treated DNA was amplified by PCR using Platinum PCR Taq DNA Polymerase (Invitrogen). Primers (5′–3′): forward (GTGGTTTTGTTTTGTTGTTAGAGAG), reverse (biotin-AAAATTCCCTAAAATTAAAAACTTCT), and sequencing (TGTTTTGTTGTTAGAGAGA) are designed with PyroMark Assay Design SW 2.0 software and obtained from Integrated DNA Technologies. Specifically, the reverse primer was biotinylated at the 5′ end. The biotinylated PCR product was captured using streptavidin-coated beads (GE Healthcare). After annealing with the sequencing primer, the single-stranded PCR product was pyrosequenced on a PyroMark Q24 Advanced Instrument (Qiagen). Average methylation index (MI) was calculated by combining the percentage of each CpG peak and dividing with the total number of CpG peaks in the pyrogram as shown previously (23).

Whole-cell lysates (20 μg/lane) were resolved in 10% SDS-PAGE from Bio-Rad. Blots were then probed with the indicated antibodies. Primary antibodies against TP63 Clone 10H7L17 (1:1,000) were from Thermo Fisher Scientific (catalog no. 703809 RRID: AB_2809251); α-SMA (1:1,000) was from Abcam (catalog no. ab5694, RRID: AB_2223021); and β-actin, Clone AC-15 (1:2,000) was from Sigma-Aldrich (catalog no. A1978, RRID:AB_476692). Secondary antibodies were from Cell Signaling Technology. Western blot images are quantified by using GeneGnome XRQ chemiluminescence imaging system and analyzed by GeneTools Analysis Software from Syngene.

The detailed procedure was reported previously (24). MCF10DCIS cells isolated from mammospheres were stained with antibodies against CD44-APC and Clone G44-26 from BD Biosciences (catalog no. 559942, RRID: AB_398683) and CD10-PE from Thermo Fisher Scientific (catalog no. 12-0106-41, RRID: AB_10714985). The stained MCF10DCIS cells were analyzed by flow cytometry using an FC500 Analyzer (Beckman Coulter) to determine the percentage of four different CD44⁻/CD10⁺, CD44⁺/CD10⁺, CD44⁻/CD10⁻, and CD44⁺/CD10⁻subpopulations.

The MFE data are presented as means ± SD. Simple comparisons between two groups were analyzed using Student t test, and comparisons of multiple groups were analyzed using one-way ANOVA. P < 0.05 was considered statistically significant. Statistical analysis was carried out using R statistical software.

To examine the effect of vitamin D compounds on cancer stem–like cells, the mammosphere forming assay was performed. MCF10DCIS cells were grown in mammosphere culture with proliferation supplements as described above for 5 days, as reported previously (25). On the basis of our previous studies, we selected an equivalent effective dose for each compound, 1α25(OH)₂D₃ (100 nmol/L) and BXL0124 (10 nmol/L). Treatment with 1α25(OH)₂D₃ and BXL0124 decreased the MFE from 1.18% in controls to 0.82% (P < 0.05) and 0.87% (P < 0.05), respectively (Fig. 1A and B). While the vitamin D compounds reduced the number of spheres, the proportion of larger spheres 100–200 μm and >200 μm was increased in the treated cultures (Fig. 1C) and the spheres in the vitamin D–treated cultures were rounder in shape (Fig. 1A).

The distribution of differential expression genes (DEG) in 1α25(OH)₂D₃ versus control and BXL0124 versus control treatment groups are shown, respectively, in volcano MA plots (Fig. 2A). Of the genes upregulated, 52.8% (371 genes) were common to the two vitamin D compounds, while 49.3% (546 genes) were common in the downregulated group (Fig. 2B). Comparing the gene expression changes of 15,331 genes being sequenced from 1α25(OH)₂D₃ to control group, a list of 12,351 genes was obtained with q value less than 0.01. Of those, 439 genes had a more than 4-fold positive (log₂ > 2) change and 703 genes had greater than 4-fold negative (log₂ < −2) change in normalized RNA expression in 1α25(OH)₂D₃ versus control group. Similarly, when comparing the gene expression changes of BXL0124 treatment group relative to control group, a list of 12,738 genes was obtained at q value of less than 0.01. Of those, 634 genes were upregulated by more than positive 4-fold (log₂ > 2) change and 948 genes were downregulated by more than negative 4-fold (log₂ < −2) change in expression. The fragments per million of differentially expressed genes in response to 1α25(OH)₂D₃ and BXL0124 treatments were normalized by log₂ and shown in a heatmap compared with control (Fig. 2C). The 25 most upregulated and downregulated genes for 1α25(OH)₂D₃ and BXL0124 cultures compared with controls are shown in Table 1 (for all RNA-seq genes, see Supplementary Data S1 and S2).

To validate the findings from RNA-seq data using qPCR analysis, we selected 12 genes of potential interest, on the basis of established relevance to cancer stemness, breast cancer progression, and chemoresistance (Fig. 2D). We first evaluated CYP24A1, a gene that is well established to be involved in vitamin D metabolism and highly induced by vitamin D compounds (26, 27). As expected, the level of CYP24A1 RNA increased with 1α,25(OH)₂D₃ treatment (82,000-fold, P < 0.001) and BXL0124 treatment (177,564-fold, P < 0.01), respectively. Genes that are known to be associated with breast cancer basal–like phenotype (KRT6A and KRT5) and a negative regulator of breast tumorigenesis (EMP1) were also upregulated, in accord with the RNA-seq results. Cytokeratin 6A (KRT6A) gene expression was upregulated 6.2-fold with 1α,25(OH)₂D₃ (P < 0.05) and 9.2-fold with BXL0124 (P < 0.05) treatment. Cytokeratin 5 (KRT5) expression level was increased 3.9-fold and 5.9-fold with 1α,25(OH)₂D₃ (P < 0.01) and BXL0124 (P < 0.001), respectively. EMP1 gene was increased 2.9-fold with 1α,25(OH)₂D₃ (P < 0.05) and 6.1-fold with BXL0124 (P < 0.01), respectively. We validated genes that are involved in epithelial–mesenchymal transition (EMT), invasion, and metastasis (S100A4 and LCN2), chemoresistance (NGFR, PPP1R1B, and AGR2), and basal breast cancer (NGFR). S100A4 gene expression was reduced by 29% with 1α,25(OH)₂D₃ (P < 0.05) and by 35% with BXL0124 (P < 0.05). 1α,25(OH)₂D₃ and BXL0124 decreased LCN2 gene expression by 97% (P < 0.001) and 97% (P < 0.01), respectively. NGFR gene expression was downregulated by 97% with 1α,25(OH)₂D₃ (P < 0.01) and by 99% with BXL0124 (P < 0.001). In our analysis, 1α,25(OH)₂D₃ and BXL0124 decreased PPP1R1B gene expression by 81% (P < 0.01) and 97% (P < 0.001), respectively. AGR2 gene expression was also reduced by 1α,25(OH)₂D₃ (by 91%, P < 0.05) and BXL0124 (by 85%, P < 0.01). GDF15 gene, implicated in the maintenance of breast cancer stem–like cells, was significantly inhibited with 1α,25(OH)₂D₃ by 94% (P < 0.01) and with BXL0124 by 95% (P < 0.01), respectively. Genes that are involved in cell migration and cytokinesis (CAPN6) and mammary gland involution and differentiation (IGFBP-5) were also analyzed. CAPN6 gene expression was decreased by 1α,25(OH)₂D₃ (89%, P < 0.01) and with BXL0124 (97%, P < 0.001) in MCF10DCIS mammospheres. IGFBP5 gene expression was increased 2-folds with 1α,25(OH)₂D₃ (P < 0.05), while reduced by 89% with BXL0124 (P < 0.05). Our qPCR analysis confirmed and validated the expression changes as observed in RNA-seq data (Fig. 2D).

To determine whether the gene expression changes were due to changes in gene methylation status and to investigate the global DNA methylation changes induced by vitamin D compounds, we performed single base pair resolution DNA methylation sequencing with the mammospheres treated with DMSO, 1α25(OH)₂D₃, and BXL0124 using Agilent SureSelect Human Methyl-seq library and Illumina NextSeq 500 platform. DNA methylation profiles were characterized using DMRfinder on the basis of a total of 176,900 DMRs for 1α25(OH)₂D₃ and 162,581 DMRs for BXL0124 treatment groups, which were annotated using ChIPseeker. Distribution of DMRs annotated by gene feature are shown in Fig. 3A. With a cut-off criterion of P < 0.05 and methylation difference >10%, we obtained 46 genes with their DMRs located in the promoter regions as shown in the heatmap (Fig. 3B). In addition, we analyzed methyl-seq data of all loci assayed. We found 83 genes with P < 0.05 with a cutoff of methylation difference > 10%, and 145 genes with P < 0.05 without cutoff (Supplementary Fig. S1). Plotting methylation differences for all genes against each other for both 1α25(OH)₂D₃ versus DMSO and BXL0124 versus DMSO showed a linear relationship between the two groups (Supplementary Fig. S2). Notably, MCF10DCIS mammospheres treated with 1α25(OH)₂D₃ and BXL0124 shared similar pattern in heatmap compared with control implying the methylation modification in these DMRs may be a regulatory mechanism for breast cancer chemoprevention.

We further performed correlation between DNA methylation profiles and RNA expression profiles in 1α25(OH)₂D₃ versus control and BXL0124 versus control. We identified a list of 105 DMRs with corresponding RNA expression data in the former group and a list of 165 DMRs with corresponding gene expression data in the latter group (log₂ 2-fold for RNA expression difference and 10% for methylation difference were used as cutoffs). Starburst plots integrating DNA methylation and gene expression are shown in Fig. 3C. Each dot represents a DMR and the corresponding DMR locations are featured by different colors. From the analysis, DMRs can be classified into two major groups: one that has direct association between DNA methylation and RNA expression, and the other group that has inverse association between DNA methylation and RNA expression. The top 10 genes that have inverse relationship between transcript level and promoter DNA methylation are shown in Table 2. Because DNA hypermethylation at CpG island in promoter region is known to silence gene expression, we selected one of the most recognized vitamin D responsive target gene, CYP24A1, from the list of genes that showed higher mRNA expression with lower CpG methylation, and then further validated its methylation status with pyrosequencing. Pyrosequencing analysis of the selected six CpG sites in the promoter region of CYP24A1 gene between 52789045 and 52789434 bps showed that CpG hypermethylation is decreased upon treatment with vitamin D compounds compared with control (average MIs DMSO = 75%, 1α25(OH)₂D₃ = 62%, and BXL0124 = 55%; Fig. 3D). This pyrosequencing finding in CYP24A1 promoter methylation correlates with significant upregulation of RNA expression in mammospheres treated with vitamin D compounds validating the DNA methyl-seq and RNA-seq results.

We performed IPA to identify upstream transcriptional regulators of the differentially regulated genes by vitamin D compounds (Fig. 4A). Using a cut-off absolute value for the z-score of >2.5, we found 30 genes regulated by 1α25(OH)₂D₃ and 27 genes regulated by BXL0124, among which TP63, VDR, CD24, CST5, and IFNB1 were regulated by both vitamin D compounds. We further identified direct upstream or downstream regulators of vitamin D receptor (VDR) in human mammary gland and breast cancer cell lines in published databases that were also noted in our RNA-seq data (cutoff for P overlap <0.01 and activation z-score >2). Figure 4B shows top candidate genes predicted to be regulated by VDR in the presence of 1α25(OH)₂D₃ or BXL0124. Genes identified in this analysis that were regulated by both vitamin D compounds include TP63, CYP24A1, CD14, TRPV6, STAT4, and FABP6. Analysis of ingenuity canonical pathways regulated by 1α25(OH)₂D₃ versus control and BXL0124 versus control is shown in Supplementary Tables S1 and S2, respectively. Using criteria of P < 0.01 and absolute value z score >2, 45 canonical pathways were regulated by 1α25(OH)₂D₃ and eight canonical pathways were regulated by BXL0124. Among them, IL6, NANOG, GP6, and LXR/RXR pathways were strongly affected.

One target gene of interest that emerged from our analyses is TP63, a member of tumor suppressor p53 family of transcription factors. TP63 is essential for development and maintenance of epithelial stem cells (28). Loss of TP63 expression promotes malignant cells migration, invasion, and distant metastasis in cancers (29). Tp63 is associated with differentiated myoepithelium-specific genes in normal breast tissue and its perturbation is observed in MCF10DCIS cells, which show significant decrease in Tp63⁺ population and impaired differentiation (30). To further examine the effects of vitamin D compounds on TP63, we analyzed protein levels in mammospheres (Fig. 4C; Supplementary Fig. S3) and observed the upregulation of TP63 protein upon treatment with vitamin D compounds. In addition, we analyzed the protein level of α-smooth muscle actin (SMA), a myoepithelial cell differentiation marker in DCIS progression in mammospheres (31). α-SMA level was upregulated in mammosphere samples treated with 1α25(OH)₂D₃ and BXL0124 (Fig. 4C; Supplementary Fig. S3). The disruption of the continuity of myoepithelial basement membrane is a prerequisite to stromal invasion of tumor cells in DCIS progression to IDC (32). Loss of CD10, a myoepithelial surface peptidase, is associated with triple-negative breast cancer and increased risk of DCIS (33, 34). We, therefore, performed flow cytometric analysis of CD10 in MCF10DCIS mammospheres treated with vitamin D compounds. We observed significant upregulation of CD10⁺ in the treated cells (Fig. 4D; Supplementary Fig. S4). Overall, our data suggest that vitamin D compounds prevent DCIS progression to IDC by reducing cancer stem–like cells (or their “stemness”) in the population, and instead inducing myoepithelial differentiation (Fig. 4E).

We initiated this study to gain insights into global changes in gene expression regulated by vitamin D compounds in MCF10DCIS. This was intended to provide clues as to the molecular mechanisms and pathways by which these compounds reduce DCIS progression to IDC. We observed that vitamin D compounds significantly decreased mammosphere formation, an indicator that they reduced the stem cell–like population of the cell's stemness. Among the genes that are commonly and differentially regulated by vitamin D compounds, BXL0124 regulated a higher number of genes compared with 1α25(OH)₂D₃, consistent with the fact that it is a more potent agonist of VDR signaling. From the RNA-seq analysis, we selected a list of 12 candidate genes with which to validate expression changes and for further study. GDF15, growth differentiation factor 15, is a member of the TGFβ superfamily (35). GDF15 increases tumorsphere formation and increases stemness markers Oct4, Sox2, and Nanog, phosphorylation of Smad2, and activation of ERK1/2 by autocrine/paracrine circuit leading to the maintenance of the breast cancer stem-like cell state (36). Both vitamin D compounds decreased GDF15 expression in mammospheres, suggesting inhibitory effects of these compounds on stemness and EMT. LCN2 (lipocalin) is a secreted glycoprotein that transports lipophilic ligands (37). Overexpression of LCN2 promotes cancer cell invasion and motility, upregulates slug, vimentin, and fibronectin, and downregulates E-cadherin indicative of inducing EMT (38). Both 1α25(OH)₂D₃ and BXL0124 reduced LCN2 at mRNA levels.

Chemoresistance is one of the hallmarks of cancer stemness (39). Nerve growth factor receptor (NGFR), a CSC marker in melanoma, along with its downstream target FGF13 (a mediator of chemoresistance), is increased in chemotherapy-sensitive melanoma cells but not in chemo-resistant cells (40). Triple-negative breast cancer cells secrete nerve growth factor that upregulates NGFR expression and inhibits apoptosis induction by chemotherapeutic agents (41). In our study, NGFR expression was decreased in mammospheres treated with vitamin D compounds, suggesting the possible role of vitamin D compounds in reducing chemoresistance by targeting CSCs. KRT5 and KRT6 are members of the keratin gene family coexpressed during differentiation of epithelial tissues (42). In breast cancer, KRT5 and KRT6a are commonly used as IHC markers of the basal-like subtype (43). Immunostaining of tissue microarrays using clinical DCIS samples revealed that cytokeratin 5 and 6 expression was inversely associated with invasion into surrounding normal breast tissue (44). Therefore, increased expression of these markers upon treatment with vitamin D compounds suggests that vitamin D induces myoepithelial differentiation in DCIS, preventing its progression to IDC. CAPN6 (calpain-6), a family member of calcium-dependent cysteine proteases, is downregulated in simple ductal hyperplasia and atypical ductal hyperplasia but upregulated in DCIS breast cancer, suggestive of its function in breast cancer progression (45). The expression level of CAPN6 is significantly decreased with vitamin D treatment, indicating its inhibitory effect in breast cancer progression. To further validate the RNA-seq data from MCF10DCIS mammosphere culture, we performed bioinformatic analysis of our data in relation to publicly available database with MCF10DCIS cells and confirmed a strong positive linear relationship (Supplementary Fig. S5). However, we observed the expression profiles of certain genes (including the ones that we have validated such as NGFR, GDF15, PPP1R1B, KRT5, KRT6A, and LAMA3) in our MCF10DCIS mammospheres are different from those of MCF10DCIS monolayer cells. Treating the mammospheres with vitamin D compounds reversed the expression of these genes. Similarly, several genes highly regulated by vitamin D compounds listed in Table 1 (LACRT, CREB3L1, FCGBP, HMGCS2, BPIFA1, MUC5AC, CLDN10, BPIFB1, and SPINK8) showed different expression profiles in mammospheres compared with monolayer cells, suggesting that these genes might play a role in stem-like–driven mammosphere growth.

Using single base pair resolution methyl-seq, we examined the DNA methylation status of the genes to determine whether their changes in expression was due to changes in promoter CpG island methylation induced by the vitamin D compounds. We found that treatment with 1α25(OH)₂D₃ and BXL0124 have shown similar patterns of changes in CpG methylation profiles. We have also identified that some of the DMRs showed direct relationship between DNA methylation and RNA expression, whereas some showed inverse relationship—while the methylation is increased, gene expression was downregulated and vice versa. Because a large body of evidence shows that CpG island promoter methylation results in gene silencing, we selected a well-established vitamin D metabolizing gene, CYP24A1, that is strongly induced by the compounds to test this hypothesis. We found that when DMRs are hypomethylated by vitamin D compounds, gene expression of CYP24A1 is upregulated. Overall, these findings suggest epigenetic regulation as a potential mechanism of vitamin D compounds for breast cancer chemoprevention.

IPA of transcriptomes revealed potential signaling pathways and downstream targets regulated by vitamin D compounds in DCIS. Among these, we selected the TP63 pathway for further analysis. Tp63, a member of tumor suppressor p53 family of transcription factors, is essential for development and maintenance of epithelial stem cells (29). Six different isoforms of Tp63 have been identified on the basis of the presence of transactivating (TAp63) or dominant-negative (ΔNp63) domains (29). Tp63-knockout mice showed impaired mammary, epithelial, and craniofacial developments, indicating its role in regulating proliferation and differentiation of these structures (46). IHC studies in normal and diseased breast tissues show that Tp63 is immunoreactive in myoepithelial cells of normal and benign lesions but not in invasive breast carcinoma (47). Recent analysis based on METABRIC cohort studies showed that the best overall survival is associated with myoepithelial mammary cell phenotype with Tp63 expression, similar to normal breast–like class (48). In addition, Tp63 is found to be highly expressed in metaplastic carcinomas of the breast with spindle cell or squamous differentiation (49). Loss of TP63 expression promotes malignant cells migration, invasion, and distant metastasis in cancers (13, 50). ΔNp63α is found to mediate CXCL12–CXCR4 pathway that activates diverse oncogenic downstream signaling pathways including PI3K/AKT, MAPK, JAK/STAT, and NF-κB. By activating CXCR4 as its transcriptional target, ΔNp63α promotes prostem cell activity and chemotaxis of breast cancer cells to metastasis sites (51). In irradiated human epithelial cells, a high expression level of key stem cell factor OCT4, a downstream target of miR34a, activates p63 by cooperating with p63 isoform (TAp63α) to promote oncogenic transformation (52). Our data showed that in DCIS mammospheres treated with 1α25(OH)₂D₃ and BXL0124, protein levels of TP63 increased, consistent with our hypothesis that the mechanism by which the vitamin D compounds inhibit progression of DCIS to IDC involves a reduction in the CSC-like population or the stemness of the cells. To explore the possibility that these effects are manifest in a more highly differentiated state of these treated populations, we examined markers of myoepithelial differentiation. We found that the vitamin D compounds increased levels of α-SMA and CD10, myoepithelial and basal progenitor cell surface markers, respectively.

In summary, we have identified epigenetic and gene expression changes regulated by vitamin D compounds in MCF10DCIS mammospheres using RNA-seq and methyl-seq techniques. RNA-seq data provided a global view of genes differentially regulated by 1α25(OH)₂D₃ and BXL0124. A list of potential genes affected included those involved in breast cancer stemness and progression. Methyl-seq data revealed epigenetic modification of gene expression in MCF10DCIS mammospheres treated with vitamin D compounds. This observation is further supported by validation of CYP24A1 gene using pyrosequencing technique. IPA identified top upstream regulators and downstream targets of the VDR including, notably, TP63. These findings identify potential key pathways that could play a significant role in DCIS progression to IDC and cancer stemness and offer targets that might be exploited for inhibition of this progression.

D. Sargsyan reports personal fees from Johnson & Johnson and personal fees from Rutgers University outside the submitted work. No potential conflicts of interest were disclosed by the other authors.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.

N.L. Shan: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing-original draft, writing-review and editing. A. Minden: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing-original draft, writing-review and editing. P. Furmanski: Conceptualization, data curation, formal analysis, supervision, investigation, writing-original draft, writing-review and editing. M.J. Bak: Formal analysis, validation, investigation, methodology, writing-original draft. L. Cai: Conceptualization, resources, data curation, software, formal analysis, methodology, writing-original draft. R. Wernyj: Conceptualization, software, formal analysis, visualization, methodology, writing-original draft. D. Sargsyan: Conceptualization, data curation, software, formal analysis, investigation, methodology, writing-original draft, writing-review and editing. D. Cheng: Software, formal analysis, investigation, methodology. R. Wu: Software, formal analysis, methodology. H.-C.D. Kuo: Software, formal analysis, methodology, writing-original draft, writing-review and editing. S.N. Li: Data curation, software, formal analysis, validation, investigation, methodology. M. Fang: Validation, investigation, methodology, writing-original draft, writing-Review and editing. H. Maehr: Conceptualization, resources. A.-N. Kong: Conceptualization, resources, data curation, supervision, funding acquisition, visualization, methodology, writing-original draft, writing-review and editing. N. Suh: Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing-original draft, project administration, writing-review and editing.

This research was supported by the National Center for Complementary and Integrative Health of the NIH grant R01 AT007036 (to N. Suh) and R01 AT009152 (to A.-N. Kong), the National Institute of Environmental Health Sciences grant ES005022, Charles and Johanna Busch Memorial Fund at Rutgers University, and the New Jersey Health Foundation (to N. Suh).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.