MCF7 cells were purchased from the American Type Culture Collection (ATCC, VA, USA). The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (GIBCO, USA) with 5% fetal bovine serum (GIBCO, USA) and 1% Penstrep (GIBCO, USA) in humidified 5% CO₂ incubator at 37°C. Cell lines were negative for mycoplasma presence (Lonza, Switzerland). Construction of wild-type JMJD6 tagged with V5 and generating stable clones of this construct in MCF7 cells (JOE) has been described earlier [denoted as JOE 1, JOE 2, and so on] (8). Cells transfected with empty vector (Vec) were used as a control. 17β-estradiol (E₂) and 4-Hydroxytamoxifen (4-OHT) (Sigma Aldrich, USA) were used at concentrations indicated in the respective experiments.

MCF7 cells were transiently transfected for 24 hours (h) with V5 tagged JMJD6 using lipofectamine 2000 (Invitrogen, CA, USA) according to the manufacturer’s instructions. For siRNA, the reverse transfection protocol was carried out using JMJD6 siRNA (J6 siRNA; 5’-gcuauggugaacacccuaatt-3’) (8), ESR1 siRNA (Dharmacon: L-003401-00) and scrambled siRNA (ScsiRNA) was used as control (Ambion: 4635). Cells were harvested 48 h after transfection.

MCF7 cells were maintained in DMEM-F12 (GIBCO, USA) without phenol red, 10% Charcoal-dextran treated fetal bovine serum (CDFBS) (GIBCO, USA), 1% Penstrep and 1µM 4-OHT for 6 months (15). Tam resistance was confirmed by increasing 4-OHT up to 5 µM (data not shown). These cells are designated as TAM R. Long term estrogen deprived (LTED) cells were generated from MCF7 cells as described earlier (16).

RNA was isolated from JOE and Vec cells by using RNAeasy Qiagen mini kit (Qiagen, Germany). Libraries were prepared using Illumina True-Seq stranded library kit. Library preparation was done at the National Genomics Core at NIBMG. Paired end sequencing (50 million reads each) was done by NOVASeq. FASTQC check was performed to remove spurious data. FastQ files were initially trimmed by CutAdapt to remove the adapter sequences and were aligned with the human genome (GrCh38.1) by Hisat2.0 algorithm. Aligned sequences were processed by Featurecounts package for obtaining RNA/transcript counts. EdgeR package was used to obtain differentially expressed genes (DEGs) using Galaxy (www.usegalaxy.org). Cluster 2.0 and Tree View programs were used for heatmap generation and visualization. Pathway enrichment analysis was done by Enrichr (https://maayanlab.cloud) and preranked analysis in GSEA (https://gsea-msigdb.github.io/gseapreranked-gpmodule/v6/index.html). Top 10 significant pathways were chosen for further comparisons. Sequencing Data has been uploaded in GEO omnibus (GSE211031). For comparison, microarray data of MCF7 cells treated with siRNA (ESR1 KD) was downloaded from GEO (GSE 27473) and analyzed in GEO2R to get DEGs. Genes with adjusted p value ≤ 0.05 and those that overlapped with DEGs from JOE cells were selected for further analysis in Enrichr. Figures and graphs were plotted using Bioconductor R, ggplot2 and ggvenn packages.

RNA from cell lines was isolated using Trizol (Ambion) as per the manufacturer’s instruction. One microgram of total RNA was reverse transcribed using SuperScript^® III First-Strand Synthesis System (Invitrogen, CA, USA). SYBR green mix (KAPA BIOSYSTEMS, South Africa) was used for quantitative real-time PCR (qRT-PCR), using primers described in Supplementary Table 1 . C_T values of gene specific primers were normalized to C_T value of β-actin. Fold change in gene expression across samples was calculated by the formula 2 ( − Δ Δ C T ) using values from ‘Vec’ cells converted to fold change 1.

Cells were lysed in RIPA buffer containing 1X Protease Inhibitor Cocktail (Sigma Aldrich, USA), quantified by BCA method (Pierce, USA) and equal amounts (50 microgram) were analyzed on 10 – 12% SDS–PAGE gels. Proteins were transferred using PVDF membrane (Sigma Aldrich, USA). Primary antibodies were incubated at 4°C overnight and suitable HRP conjugated secondary antibodies for 1 hour at room temperature. Signals were detected using HI-ECL reagent (Biorad, USA) using Chemidoc (Biorad, CA, USA). Antibodies used were JMJD6 (PSR sc-28348, Santa Cruz Biotechnology; CA, USA; 1:1000), V5 for only exogenously expressed V5-tagged JMJD6 (R960-25, Invitrogen; 1:1000), ER-α (HC-20; F10) (Santa Cruz Biotechnology 1:1000), β Actin (A2228, Sigma Aldrich, USA;1:500), MAPK (ERK1/2) (137F5, Cell Signaling Technologies, 1:1000), Phospho-MAPK (ERK1/2) (Thr202/Tyr204) (197G2, Cell Signaling Technologies 1:1000), RET (3220, Cell Signaling Technologies 1:1000), Phospho-Akt (Ser473) (D9E4060 Cell Signaling Technologies 1:1000), Akt (Cell Signaling Technologies, 1:1000).

Cells were plated in chamber slides (Eppendorf, Germany) for 24 h. 10% Neutral buffered formalin (NBF) was used as a fixative. Cells were washed in 1X PBS followed by permeabilization (1X PBS+ 0.2% Tween-20), washed with Glycine (Sigma Aldrich, USA) and non-specific sites were blocked with 10% normal goat serum for 1 h. Cells were incubated with the primary antibody overnight (JMJD6 1:250, ER 1:250), washed and incubated with the appropriate secondary antibody (1:500) for 1 h. Signals were captured in a confocal microscope (Nikon Eclipse Ti2E). All experiments were performed in triplicates.

Vector and JOE cells were cultured in E₂ deprived media for 72 h before transfection. 10⁵ cells/well were plated in a 24 well plate, p(ERE)₂ luciferase and pRL-TK plasmids were transfected using lipofectamine 2000 (Invitrogen, USA). Cells were treated with either ethanol (Veh), 10 nM E₂, or 1 μM 4-OHT for 24 h and lysed in 1X Passive Lysis buffer (PLB). Reporter gene activity was measured using a Dual Luciferase Assay System (Promega, USA). ‘Veh’ samples were assigned the value of 1 and the remaining treatments were normalized to obtain fold change in luciferase activity.

For TCGA analysis, mRNA expression z-scores relative to all samples (log RNA-seq V2 RSEM) from ER+ Breast cancers from Breast Invasive Carcinoma (Cell, 2015) were used. Correlation plot was generated using www.cbioportal.org. Kaplan Meier plots were generated at www.kmplot.com/analysis using ER+ Breast cancer data. A ratio of JMJD6/ESR1 expression was used to divide the patients into two groups and their association with survival outcome was estimated.

Cell proliferation assay was performed with and without E₂ treatment for both Vec and JOE cells. In brief, cells were starved in E₂ depleted media for 3 days and plated at a density of 10⁴ cells per well (Day 0). 10 nM E₂ was added after 24 h and media was replaced each day. Cells were counted in Hemocytometer at time points mentioned in the figure. Cell viability assay was performed by CCK8 kit (Sigma Aldrich, USA). Briefly, cells were starved for 3 days in E₂ deprived media at a density of 10⁴ cells/well. Different doses of 4-OHT as indicated were added and cells were incubated for 72 h. The media was changed every day and absorbance was measured at 450 nm.

TAM R cells were reverse transfected with JMJD6 siRNA as described above. ScsiRNA was used as a control. After 48 hours, a portion of the cells were subjected to western blot analysis to confirm downmodulation of JMJD6. Remaining TAM R cells were treated with Ethanol (Veh) or 1 µM 4-OHT (Tam). Cell numbers were counted in Hemocytometer at time points mentioned in Figure 7D .

Each experiment was performed independently at least three times. Data are presented as mean ± standard deviation for each experiment except for RNA-Seq analysis. Statistical analysis was performed using GraphPad prism software. Comparisons between two groups were performed using student’s t-test and data was considered statistically significant at p≤ 0.05. Significance in all figures is indicated as follows: * p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001 and ****p ≤ 0.0001.

RNA-seq of JOE and Vec cells identified 4497 significantly differentially expressed genes (adjusted p value ≤ 0.05) following EdgeR analysis. This set of genes is denoted as the “JOE” set. A volcano plot of these genes is shown in Figure 1A . 2208 genes were upregulated and 2289 genes were downregulated by JMJD6. To find the highly regulated genes, a 2-fold change(2FC) cut off was imposed on the significantly changed DEGs and 1131 genes were obtained ( Supplementary Table 2 ). We studied pathway enrichment using both GSEA preranked analysis and Enrichr methods and identified E2F signaling pathway, G2-M checkpoint, mTORC1 signaling, ESR1 mediated early and late gene expression, TGF-β signaling and others as significantly enriched pathways ( Supplementary Figure 1 ). Of these, we have previously described contributions of JMJD6 to regulation of E2F target genes and the TGF-β signaling pathway (7, 8).

We explored estrogen early/late response pathways further and compared the JOE set genes with publicly available dataset for ER regulated genes (MCF7 cells treated with ESR1 siRNA, ESR1 KD set) (17). This dataset was chosen since the JOE dataset was obtained in the absence of E₂ treatment. In the ESR1 KD dataset, 13516 DEGs were obtained (adjusted p value ≤ 0.05). Intersection of these two sets indicated that 3412, that is 76% of JOE set genes, overlapped with ESR1 KD set ( Figure 1B ). Though overlapping, the directionality of gene expression (induced vs repressed) was not conserved, and genes could be divided into 4 subsets based on the treatment and direction of change in gene expression. Figure 1C graphically represents the four quadrants and pathway enrichment analysis from each subset of genes. Details such as number of genes in each quadrant and the Enrichr output for top 10 pathways are represented as a heatmap and bar graph in Supplementary Figure 2 . Both JMJD6 and ESR1 showed concordance in the regulation of the Estrogen early/late gene response pathway. Similarly, they both suppressed the TGFβ signaling, apoptosis and cholesterol biosynthesis pathways. Androgen Receptor (AR) response genes were suppressed by JMJD6, whereas, ER appeared to enhance this pathway. However, certain pathways showed distinct regulation in the JOE and ESR1 KD sets. JMJD6 induced genes (n=1634) showed enrichment for E2F, G2-M transition and mTORC1 signaling pathways ( Figure 1C ). The same set of genes was split into two groups in the ESR1 KD set; 742 genes were less expressed in presence of ER and 892 showed higher expression. For example, the G2-M transition genes were divided into two groups-repressed by ER and induced by ER. Further, while JMJD6 induced Mitotic spindle genes and TNFα signaling pathway, ER suppressed these genes. This suggests that JMJD6 probably induced a larger repertoire of cell proliferation genes as compared to ER ( Supplementary Table 3 ). These data suggest that high JMJD6 may offer ER+ cells an additional advantage via executing an alternate pathway to cell proliferation and this pathway appears to be independent of the E₂ -ER axis.

We validated candidate genes from DEG list that were a) previously known targets of JMJD6 (HOTAIR and MYC); b) genes that are known to be regulated by ER and JMJD6 upon E₂ stimulation (TFF1, NRIP1, SMAD7, FOXC1, GREB1) and c) genes that have JMJD6 binding sites in their regulatory region (AURKA, AURKB) (7, 14, 18, 19). Data for all candidate genes in another JOE clone and western blots to represent successful down regulation of JMJD6 and ER post siRNA treatment is shown in Supplementary Figure 3 . We found an increase in the expression of HOTAIR, MYC and AURKA in JOE cells which could be reversed by JMJD6 siRNA treatment ( Figure 2A ). qRT-PCR data obtained for ESR1 regulated genes was in accordance with microarray ESR1 KD data. ESR1 siRNA decreased the expression of MYC and AURKA. However, change in AURKB levels was not statistically significant. Interestingly FOXC1, GREB1, TFF1, NRIP1 and SMAD7 genes that are induced by and require binding of both JMJD6 and ER to their regulatory sites in the presence of E₂, were repressed in JOE cells. Further, except for FOXC1, the expression of these genes was decreased in MCF7 cells following siRNA mediated knock down of JMJD6 ( Figure 2B ; Supplementary Figure 3 ). This suppression in ER target expression both in presence and absence of JMJD6 was curious, since previously published data suggested that JMJD6 was necessary for ER target gene expression in presence of E_2. It was more likely that ER target genes would display higher expression in the presence of recombinantly expressed JMJD6 or they remain unaltered since JOE cells were not treated with E₂ in our study (14). Unless 1) supraphysiological levels of JMJD6 negatively regulated ER transactivation function in the absence of hormone or 2) JMJD6 suppressed basal expression of ER target genes in absence of E₂ or 3) high JMJD6 expression affected the levels of ER itself. We tested some of the possibilities and checked the expression levels of ESR1 mRNA in JOE cells. We found a decrease in ESR1 expression ( Figure 2C ), and in MCF7 cells transfected with JMJD6 siRNA ESR1 levels appeared to be higher than cells transfected with scrambled control ( Figure 2C ). These data suggested that JMJD6 may negatively affect the ER axis in breast cancer cells by depleting ER expression levels.

In JOE cells, immunoblot analysis indicated that expression of ER protein was lower than in Vec cells ( Figure 3A ). This observation was consistent in several clones stably expressing recombinant JMJD6 tagged with V5 ( Supplementary Figure 4A ). In addition, decrease in ER protein was evident after transient transfection of JMJD6 in MCF7 cells for a period of 24 h ( Supplementary Figure 4B ). Using immunofluorescence, we found that nuclear ER staining in JOE cells was lower as compared to Vec cells, and no ER specific fluorescence was detected in the cytoplasm ( Figure 3B ). Therefore, JMJD6 decreased ER but did not appear to affect its cellular localization. Knock down of JMJD6 by siRNA resulted in a slight but consistent increase in ER levels in MCF7 cells ( Figure 3C ).

Next, we tested ERE luciferase activity in JOE and Vec cells to determine if the residual ER in JOE cells could be activated by hormone stimulation and if high JMJD6 suppressed basal ER activity. As expected from the reduction in ER levels, the basal ERE-Luc activity was lower in JOE cells than Vec cells, however, it increased by at least 10-fold in the presence of E_2, and 4-OHT treatment negated this induced activity ( Figure 3D ). A similar 10-fold induction was evident in Vec cells following E₂ treatment though the basal ERE-luc activity in these cells was far higher due to presence of substantial ER protein levels. We conclude that thought elevated levels of JMJD6 impinge upon total ER RNA/protein levels, decrease basal activity of ERE-luc, they appear not to interfere with the capacity of ER to transactivate gene expression in the presence of E_2.

Since JMJD6 appeared to suppress ER, we explored the clinical ramifications of this observation in publicly available datasets. In ER+ patient samples from TCGA breast cancer data, JMJD6 and ESR1 showed a trend towards negative correlation (Spearman, r² = -0.23) ( Figure 4A ). We studied if this negative correlation had any bearing on patient outcome. We used the ratio of JMJD6 and ESR1 mRNA and JMJD6 and ER protein to test their association with survival. A high ratio would indicate high JMJD6, low ER samples and a low ratio would identify samples having low JMJD6, and hence higher level of ER. Interestingly, in more than 2500 ER+ samples, a high JMJD6 to ESR1 (RNA level) and a high JMJD6 to ER ratio (protein level) significantly associated with poorer prognosis and patient survival ( Figures 4B, C ).

Reduction in ER amounts is a well-known indicator of tamoxifen insensitivity (6, 20, 21). Other pathways involved in lower responsiveness include higher expression of RET (REarranged during Transfection), a receptor tyrosine kinase, and upregulation of the Mitogen Activated protein kinase (MAPK) pathway (22–25). In JOE cells, RET was induced at both RNA and protein level ( Figures 5A, B ). Knockdown of either JMJD6 or ESR1 in MCF7 cells using gene specific siRNAs suppressed the expression of RET ( Figure 5A ). This is not surprising since it is a well-known transcriptional target of ER (26). In addition, Tam resistant tumors display higher MAPK (ERK1/2) levels (27). Immunoblot analysis showed higher protein expression for both total ERK1 and p-ERK1 in JOE cells as compared to Vec cells ( Figure 5C ). Based on these data, we propose that JMJD6 may be involved in predisposing cells to Tam insensitivity by decreasing ER, increasing RET expression and total ERK1 levels in ER+ breast cancer cells.

We studied the expression of JMJD6 in two models of Tamoxifen resistance that are derived from parental MCF7 cells; TAM R cells and LTED cells (28). LTED cells have different sensitivity to estrogen (LTED-Q- quiescent; LTED-H- hypersensitive; LTED-I- Independent) (29). Our premise was that if JMJD6 associated with Tam insensitivity, TAM R and LTED-I cells would show higher JMJD6 expression. Accordingly, TAM R cells displayed elevated levels of JMJD6 and pAKT but not ERK1/2, when compared to parental MCF7 cells ( Figure 6A ). LTED-Q and LTED-H cells did not possess high JMJD6 protein levels. Interestingly, the property of Tamoxifen’s agonist behavior in LTED-I cells, coincided with an increase in JMJD6 protein expression ( Figure 6B ).

Further, comparison of DEGs from JOE cells with published tamoxifen resistance gene signatures was undertaken. Analysis of gene expression data from 136 tamoxifen resistant patients had reported a signature of 36 genes (30). Of these, 18 genes appeared in the JOE dataset. Of note is the compete conservation of the direction of gene regulation (high expression vs low expression) in JOE cells and patient samples ( Figure 6C ). A second 30 gene signature has been derived from Tam R cells originating from MCF7, T47D and ZR75.1 (31, 32). Of these 30 genes, 20 genes were present in the JOE dataset with complete conservation of the direction of gene expression. These data clearly indicate that high JMJD6 levels may impose Tam insensitivity in ER+ breast cancer cells.

The ERE luciferase activity increased in JOE cells following E₂ treatment indicating that the ability of lower levels of ER expression did not affect its ability to transactivate gene expression. Since ER induces cell proliferation in MCF7 cells, we assessed if E₂ could do so in JOE cells. Further, we determined if E₂ was capable of increasing ER levels to mediate this change. Our data reveal that E₂ not only increased JOE cell proliferation but also recovered the expression of ER in these cells ( Figures 7A, B ). Since JOE cells were responsive to E_2, we studied the effect of 4-OHT on cell viability in Vec and JOE cells. At each dose of 4-OHT, JOE cells consistently showed a trend towards better survival than Vec cells, except at the highest dose of 10 μM Tam ( Figure 7C ). As shown, while only 30-45% of Vec cells survived at 1 and 2.5 μM doses, 60-75% of JOE cells remained alive. Therefore, JOE cells display lower sensitivity to Tam induced cytotoxicity

To check if JMJD6 directly contributed to Tam resistance, we treated TAM R cells with JMJD6 siRNA. 48 hours are plating, reduction in JMJD6 levels was confirmed by immunoblotting ( Figure 7D ). 1 μM Tam was added at this time and cell proliferation was monitored for the next 24, 48 and 72 h ( Figure 7D ). JMJD6 depletion led to a significant decrease in TAM R proliferation as compared to ScsiRNA treatment. More importantly, treatment of cells with both JMJD6 siRNA and Tam decreased cell proliferation further. Together these data raise the possibility that higher JMJD6 expression predisposes ER+ cancers to endocrine therapy resistance.