RUNX2 expression in thyroid and breast cancer requires the cooperation of three non-redundant enhancers under the control of BRD4 and c-JUN.

Aberrant reactivation of embryonic pathways is a common feature of cancer. RUNX2 is a transcription factor crucial during embryogenesis that is aberrantly reactivated in many tumors, including thyroid and breast cancer, where it promotes aggressiveness and metastatic spreading. Currently, the mechanisms driving RUNX2 expression in cancer are still largely unknown. Here we showed that RUNX2 transcription in thyroid and breast cancer requires the cooperation of three distantly located enhancers (ENHs) brought together by chromatin three-dimensional looping. We showed that BRD4 controls RUNX2 by binding to the newly identified ENHs and we demonstrated that the anti-proliferative effects of bromodomain inhibitors (BETi) is associated with RUNX2 transcriptional repression. We demonstrated that each RUNX2 ENH is potentially controlled by a distinct set of TFs and we identified c-JUN as the principal pivot of this regulatory platform. We also observed that accumulation of genetic mutations within these elements correlates with metastatic behavior in human thyroid tumors. Finally, we identified RAINs, a novel family of ENH-associated long non-coding RNAs, transcribed from the identified RUNX2 regulatory unit. Our data provide a new model to explain how RUNX2 expression is reactivated in thyroid and breast cancer and how cancer-driving signaling pathways converge on the regulation of this gene.

INTRODUCTION

Functional genome analysis revealed that gene expression is far more complicated than expected and requires a continuous and widespread regulatory landscape (1). A large part of the expression regulatory function of the genome resides within enhancers (ENHs), small segments of DNA that serve as operational platforms to recruit transcription factors (TFs) (2). Once engaged by their associated TFs, ENHs promote transcription of target genes by interacting with their specific promoter (3,4). However, ENHs are not just collections of TF binding sites for target promoters, but they are also sites of active transcription for many non-coding RNAs and central hubs for the transcriptional machinery and for complexes that control chromatin conformation and function (5,6).

During embryogenesis and in cancer progression, key genes are controlled by stretches of multiple ENHs in close genomic proximity (10–12 kb) called super-ENHs (7). This multiplicity ensures strongest transcriptional activity and increases the precision of gene expression regulation. Beside super-ENHs, evidence exists that multiple discrete ENHs intersperse in the genome may converge on the regulation of the same gene with similar spatio-temporal profiles. It is now established that this multiplexed transcriptional organization is fundamental for the correct execution of developmental pathways. The use of multiple ENHs, may help ensuring the precision of embryonic patterning, contribute to phenotypic robustness and represent functional platforms to support evolution and genetic novelty (8–10).

Many transcription factors that govern tissues and organs morphogenesis are hijacked during cancer progression (11).

RUNX2 is a member of the mammalian RUNT related transcription factor family, necessary during embryogenesis for skeletal development (12–14) and for the morphogenesis of other organs like breast and thyroid (15,16). RUNX2 is increasingly recognized in cancer biology for its oncogenic properties and many studies (17–21) have linked a de-regulation of RUNX2 function with progression and metastasization of different types of human tumors. The RUNX2 gene encodes two major isoforms starting from two alternative promoters (22,23). RUNX2 isoform I, controlled by the proximal P2 promoter is the major RUNX2 isoform in tumor cells (19,20,24). The regulatory mechanisms that control the activity of the P2 promoter and that lead to RUNX2 re-expression in cancer remain widely unknown. We recently showed that the P2 promoter has a limited transcription activity in different cancer models, suggesting that RUNX2 expression in cancer relies on still uncharacterized regulatory elements across the genome (20). In this work, we explored RUNX2 locus chromatin organization, searching for novel ENHs that contribute to RUNX2 expression regulation in cancer cells.

MATERIALS AND METHODS

Cell cultures and treatments

BCPAP (B-Papillary Carcinoma) (25) and TPC1 (Thyroid Papillary Carcinoma-1) (26) cell lines were obtained from Prof. Massimo Santoro (University of Naples, Naples, Italy). MCF7 (Michigan Cancer Foundation-7) (27)cell line were obtained from Dr Massimo Broggini (IRCCS-Istituto di Ricerche Farmacologiche Mario Negri, Milan, Italy). MDA-MB-231 (MD Anderson-MB231) (28) cell line was obtained from Dr Adriana Albini (Scientific and Technology Pole, IRCCS MultiMedica, Milan). All cell lines were authenticated by SNP profiling at Multiplexion GmbH; date of last authentication report is 9 December 2014. All cell lines were grown at 37°C/5% CO2 in DMEM with 10% fetal bovine serum. Unless otherwise specified, cells were treated for 24 h with JQ1 1 μmol/l (purchased from Sigma-Aldrich) or DMSO (Sigma-Aldrich).

Quantitative real time-PCR

Total RNA was purified with Maxwell® RSC simplyRNA Cells (Promega) and retrotranscribed using the iScript cDNA kit (Bio-Rad). Quantitative real time-PCR (qRT-PCR) was conducted using Sso Fast EvaGreen Super Mix (Bio-Rad) in the CFX96 Real Time PCR Detection System (Bio-Rad). See Supplementary Table SI for qRT-PCR primers.

Western blot analysis

Western blot analysis was performed as previously described (29). Antibodies used were rabbit anti-c-JUN (ab31419; Abcam), goat anti-RUNX2 (AF2006, R&D Systems), rabbit anti-BRD4 (A301-985A50, Bethyl), mouse anti-beta-actin (AC-15; Sigma-Aldrich), mouse anti-alpha-tubulin (sc-8035, Santa Cruz), horseradish peroxidase-conjugated anti-rabbit and anti-mouse (GE Healthcare).

Plasmid vectors

All the putative Enhancers and deletion mutants evaluated by Luciferase assays were amplified by PCR and cloned into pGL3-sP2 vector (previously described (20)) downstream the luciferase gene, between BamHI and SalI restriction sites. Point mutations in ENH11.2A1A and ENH13B2 were generated by site-directed mutagenesis, using primers containing the mutated nucleotides to amplify the major fragments. See supplementary Table SIII for primers used.

c-Jun DN expressing vector was a kind gift of Dr Mirko Marabese (IRCCS-Istituto di Ricerche Farmacologiche Mario Negri, Milan, Italy).

Transfection and luciferase assay

Cells were co-transfected with pGL3 reporter vector and pRL-TK vector using Lipofectamine 2000 (Thermo Scientific). Forty-eight hours after transfection, cells were harvested and luciferase activity was measured using the Dual- Luciferase Reporter Assay System (Promega) in a GloMax Discovery Luminometer (Promega), according to the manufacturer's instructions. For each sample, firefly luciferase activity was normalized on Renilla luciferase activity and transactivation of the various reporter constructs was expressed as fold induction on pGL3-basic Empty Vector activity.

Bioinformatic analysis

For the identification of putative enhancers, we analyzed a 335 kb genomic region spanning the RUNX2 gene using ENCODE project annotations, integrated in Genome Browser (http:// genome.ucsc.edu/). We used ENCODE regulation super-track displaying enrichment in H3K27Ac, H3K4me1 and H3K4me3 histone marks in a panel of seven cell lines of different origin (GM12878, H1-hESC, HSMM, HUVEC, K562, NHEK, NHLF) and DNase I sensitivity in a panel of 125 cell lines. To analyze sequence conservation across species of the RUNX2 locus we used the Placental Mammal Basewise Conservation by PhyloP track. The source of the H3K27Ac tracks for MCF7, A549 and HCT-116 cells was ENCODE, while for the H3K27Ac binding in MDA-MB231 cells we retrieved the GSM1693018 track from the GEO database (30). For BRD4 ChIP-Seq profiles in MCF7, SUM159 and HCC1395 cells we retrieved the following tracks from GEO database (31,32): GSE55921, GSM1842711 and GSM1842697. For c-JUN binding in MDA-MB231 cells we used GSM1700784 track from the GEO database (33). Putative enhancer regions were identified on the basis of conservation in mammals, the presence of DNAse I hypersensitivity clusters, and the presence of specific histone modifications associated with enhancer regions (H3K27Ac and H3K4Me1). Transcription factors–binding sites analysis on the core of ENH11 and ENH13 was performed by quering the TRANSFAC Professional database (http://www.biobase-international.com/) using the Match algorithm and applying the matrices’ vertebrate_non_redundant_minSUM’.

For pathway analysis, TFs list obtained by TRANSFAC analysis of ENH11 and ENH13 binding sites were included in a enrichment analysis performed through the use of IPA (Ingenuity® Systems, www.ingenuity.com). The four most significant TOP canonical pathway were considered.

Protein-protein interaction analysis was performed using STRING database (http://string-db.org/) (34). For both ENH 11 and ENH 13, all TFs obtained by TRANSFAC were included in the analysis. Only the experimentally determined interaction and the protein homologies were considered.

Chromatin immunoprecipitation (ChIP)

After cross-linking with 1% formaldehyde, cells were lysed and nuclei extracts were sonicated using a Bioruptor® Pico sonicator (Diagenode). Chromatin was precipitated with Magna ChIP™ Protein G Magnetic Beads (16-662, Millipore) and the appropriate antibody (Supplementary Table SIV). The immunoprecipitated DNA fragments were quantified by qPCR. For each experiment, a chromatin amount corresponding to 1% of chromatin used for immunoprecipitation was kept as input control. Each qPCR value was normalized over the appropriate input control and reported in graphs as input %. (qPCR value/input value × 100).

siRNA transfections

c-JUN and BRD4 Silencer Select RNAi (20 nmol/l) and control oligos (Thermo Scientific) were transfected using RNAiMax Lipofectamine (Thermo Scientific). Cells were harvested 48 h after transfection for qRT-PCR.

3C-Chromosome conformation capture

3C assays were performed as described in Naumova et al. (35) with some modifications. As restriction enzyme we used HindIII (New England Biolabs, Ipswich, MA, USA). The anchor fragment used to analyze the loops between Runx2-P2 promoter and the Enhancers regions spans from chr6:45 388 380 to 45 401 628 (13.249 bp). ∼107 TPC1 cells were harvested, fixed with 1% formaldehyde (Sigma Aldrich, St. Louis, MO, USA) for 5 min at room temperature. Formaldehyde crosslinking was quenched by 0.15 M glycine for 5 min at room temperature. Nuclei were obtained by 5 min incubation at 65°C in Lysis Buffer (10 mM Tris–HCl pH 8, 10 mM NaCl, 0.2% NP-40) with Protease Inhibitor Cocktail EDTA-Free (Biotool, Munich, Germany). Nuclei were treated with 0.3% SDS for 1 h at 37°C shaking (900 rpm), the reaction was stopped by 2% Triton, for 1 h at 37°C shaking. Nuclei were collected and, for each sample, digestion was performed in two separate tubes, with 900 U HindIII each at 37°C overnight with constant shaking. The enzyme reaction was halted by an incubation at 65°C for 25 min in the presence of 1.6% SDS. Samples were diluted 10-fold in Ligation Buffer (New England Biolabs, Ipswich, MA, USA), One percent triton was added and the samples were incubated 1 h at 37°C, shaking. Proximity-mediated ligation was performed with 10 000 U of T4 DNA Ligase (New England Biolabs, Ipswich, MA) dividing the samples in five tubes, 4 h at 16°C, 30 min at room temperature. The samples were reverse crosslinked by overnight incubation with 300 μg Proteinase K (Promega, Madison, WI, USA) at 65°C. Ligated chromatin was extracted and purified by phenol–chloroform followed by ethanol precipitation. Digestion efficiency was determined by qPCR, using a series of primers flanking the restriction sites (C) and a control internal primer set (R). Comparing an aliquot of sample retrieved before restriction enzyme digestion (UNDIG), and an aliquot retrieved after the overnight digestion (DIG) we calculated the %digestion = 100 – 100/[(CtR-CtC)DIG – (CtR-CtC)UNDIG]. Only samples digested for >60% were further analyzed. Samples concentration was determined by qPCR with a set of internal primers that amplify far from any restriction site. Samples purity was checked by a titration curve using 3C anchor primer and a primer of a nearby region. 3C primers and probes were designed by Primer3 software and are available in Supplementary Table SII.

Interaction frequencies were normalized by assessing fold change of 3C PCR amplification product (SsoAdvanced, Biorad) of samples compared to randomly ligated HindIII-digested Bacterial Artificial Chromosome (BAC) library that span the RUNX2 locus. The BAC library is constituted by the following clones: rp11-1019c24, rp11-1141H14, ctd-3210G7, rp11-1019C24, ctd-2346D15 (Thermo Fisher Scientific, Waltham, MA, USA). The clones were mixed in equimolar concentrations, digested and ligated. To determine statistical significance of the observed interactions, the signal of fragments corresponding to RUNX2 ENHs was compared to the signal of proximal regions.

CRISPR/Cas9

For the deletion of the three RUNX2 functional ENHs we first generated a BCPAP cell line overexpressing the Cas9 enzyme under an inducible promoter by transfecting BCPAP cells with the pCW9-Cas9 plasmid. pCW-Cas9 was a gift from Eric Lander & David Sabatini (Addgene plasmid # 50661) (36). After puromycin selection we isolated a set of clones and we chose the one with maximum Cas9 expression in >95% of cells (BCPAP clone #1). For each ENH we designed flanking sgRNAs using the CRISPR design tool (http://crispr.mit.edu/) and each sgRNA has been cloned in the pX330-U6-Chimeric_BB-CBh-hSpCas9 plasmid (from which we removed the Cas9 coding sequence) in the BbsI clonig site. pX330-U6-Chimeric_BB-CBh-hSpCas9 was a gift from Feng Zhang (Addgene plasmid # 42230) (37). See Supplementary Table SI for sgRNAs sequence. pX330 plasmids containing sgRNAs flanking each ENH have been co-trasfected in the BCPAP #1 clone using Lipofectamine2000 (Thermo Scientific). Cas9 expression was induced by 5 ng/ml Doxycicline (Sigma-Aldrich) for 5 days. The cells have been seeded at single cell density in 96-well plates and grown for 3 weeks. Each clone has been screened by standard PCR with primers flanking the expected deletion (see Supplementary Table SI for the sequence of primers). The clones showing a band corresponding to the deleted allele have been further analyzed by droplet digital PCR (ddPCR) to establish the copy number of the wild-type allele. Inside the deleted region of each ENH, we designed a ddPCR assay comprising primers and a TaqMan probe conjugated with the FAM fluorophore. Outside the deleted regions we designed a reference assay comprising primers and a HEX-conjugated probe. Each clone was analyzed by QX200 ddPCR System (Biorad) using the appropriate assay and copy number was established in comparison with the reference assay.

RACE

To determine the sequence of the full RAIN transcripts we perfomed 5΄ RACE and 3΄ RACE using the SMARTer RACE 5′/3′ (Clontech) kit following the manufacturer's instructions. Briefly, 1 μg of DNAse-treated RNA extracted from TPC1 cells has been retro-transcribed to generate 5′ RACE-ready and 3′ RACE-ready cDNA. 5′ ends and 3′ ends of RAIN transcripts have been amplified using the RAIN_RACE_R and RAIN_RACE_F primers respectively (see Supplementary Table SI for sequence) and a touchdown PCR program. Amplified fragments have been extracted from agarose gel, cloned into the pRACE plasmid and sequenced.

Patients selection and RUNX2 ENHs mutational analysis

This study was approved by the local Ethical Committee (protocol no.: 2014/0014425 of 06/05/2014).

Fifty-four well-differentiated PTCs were retrieved from the archive of the Pathology Unit of the Azienda USL-IRCCS of Reggio Emilia, comprising over 2700 thyroid cancers collected in the past 30 years. Tumors were classified according to the criteria recommended by the World Health Organization Classification of Tumors to exclude high-grade carcinomas, and staged according to the American Joint Committee on Cancer Staging Manual, seventh edition. Among the 54 selected PTCs, 26 developed distant metastases outside the neck (DM-PTC) while the remaining 28 did not develop distant metastases, based on a follow-up of 7 years, and were selected to match DM-PTCs by age. All patients underwent total thyroidectomy, and in 80% of the cases, they also underwent ipsilateral central neck dissection. Follow-up for all cases ranged from 13 to 372 months.

DNA from formalin-fixed paraffin-embedded (FFPE) tumor tissue was extracted using the FFPE Plus LEV DNA Purification Kit (Promega, Madison, WI) and the Maxwell 16 instrument (Promega).

DNA quality and quantity were assessed using Nanodrop and Qubit (Thermo Fisher Scientific, Waltham, MA). A TruSeq custom panel (Illumina, San Diego, CA, USA) was designed for amplicon multiplexed targeted resequencing of P2, ENH 3, ENH 11, ENH 13. In the panel design a dual strand approach was implemented to eliminate frequent erroneous mutation call in FFPE samples owing to the fixation-derived DNA deamination. The coordinates of the regions included in the custom panel design are summarized in Supplementary Table SV. After DNA quality control test libraries were generated, quantified and normalized as indicated in TruSeq custom panel low input library prep reference guide. Libraries were sequenced with the MiSeq Desktop Sequencer instrument (Illumina) and data were analyzed with the AmpliconDS protocol, VariantStudio (Illumina) and Integrated Genomics Viewer 2.3 software (https://www.broadinstitute.org/igv/). Somatic variants were identified by selecting mutations with no reported minor allele frequency in germline mutation repositories according to the VariantStudio software.

ChIP on human thyroid tumor samples

Fresh frozen thyroid tumor and normal thyroid samples for this analysis were retrieved from the Research Biobank of Arcispedale S. Maria Nuova-IRCCS of Reggio Emilia. Chromatin from patient tissues was prepared following a modified protocol. Frozen tissue was minced with a scalpel blade and resuspended in LB1 (50 mM HEPES, 1 mM EDTA, 140 mM NaCl, 0.25% Triton X-100, 0.5% NP-40, 10% glycerol). A Dounce Homogenizer was used to obtain a single cell suspension. Cells have been centrifuged and resuspended in LB2 to isolate nuclei (10 mM Tris pH 8, 1 mM EDTA, 200 mM NaCl, 0.5 mM EGTA). Nuclei have been crosslinked in 1% Formaldehyde for 10′ at RT. Cross-linking has been stopped by incubating with glycine 0.125 M. Crosslinked nuclei have been lysed in nuclei lysis buffer and chromatin immunoprecipitation has been performed as described for cell lines.

Each value corresponding to H3K27Ac enrichment in RUNX2 P2 and ENHs in normal and tumor thyroid tissues has been first normalized over input and then over an intergenic negative control region to allow a better comparison between different samples showing different immunoprecipitation efficiencies. In parallel RNA was collected for qRT-PCR expression analysis of RUNX2, c-JUN and TEAD1.

Statistical analysis

Statistical analysis was performed using GraphPad Prism Software (GraphPad). Statistical significance was determined using the Student's t-test. Each experiment was replicated two to five times.

RESULTS

Mapping active ENHs in the RUNX2 genetic locus in cancer cells

We searched the annotation data of the ENCODE project to identify elements with features of potential ENHs in the genomic regions surrounding the RUNX2 locus. Based on H3K27Ac, H3K4Me1, H3K4Me3, DNAse sensitivity and conservation tracks from seven cell lines of different origin, we identified 14 putative ENHs that we named ENH1–14 (Figure 1A, Supplementary Figure S1A). We also analyzed H3K27Ac tracks from four epithelial cancer cell lines and osteoblast (Supplementary Figure S1B). As we previously reported, A549, MCF7 and HCT116 express low RUNX2 levels and, accordingly, the H3K27Ac enrichment profile is low along the entire RUNX2 locus. By contrast, the MDA-MB231 cell line, in which RUNX2 expression is higher, displayed a H3K27Ac enrichment profile similar to the one resulting from the merge of the seven cell lines. ENH1–14 were enriched in H3K27Ac also in osteoblasts, confirming that these regions may play a relevant role in RUNX2 regulation in many contexts. We previously showed that ENH3, differently from ENH1 and ENH2, is transcriptionally active and cooperates with the P2 promoter to drive RUNX2 expression in cancer cells (20).

Figure 1.

Identification of active ENHs in the RUNX2 locus. (A) Schematic representation of RUNX2 genomic locus showing the position of the two promoters (P1 and P2) and the putative ENHs. Tracks corresponding to the histone modifications used for the identification of the ENHs are displayed. The diagram was obtained by modification of the genome browser view (http://genome.ucsc.edu). (B) qRT-PCR and Western blot analysis of RUNX2 expression in TPC1, BCPAP, MDA-MB231 and MCF7 cell lines. (C–E) ChIP analysis of RUNX2 P2 promoter and ENH3–14 regions with anti-H3K4Me3 (C), anti-H3K4Me1 (D), anti-RNA-PolII and anti-H3K27Ac (E) antibodies in TPC1, BCPAP, MDA-MB231 and MCF7 cells. An unrelated DNA region upstream of the P2 promoter was used as negative control. The bars represent the average enrichment of the indicated genomic regions in the immunoprecipitated DNA expressed as percentage of the input. All data are expressed as mean values ±SEM. N = 3. (F) Luciferase analysis of ENH3–14 regions activity in TPC1, BCPAP, MDA-MB231 and MCF7 cells. Cells were transfected with the indicated pGL3 constructs. The values represent the average fold change of luciferase activity in cells transfected with the pGL3-P2 or pGL3-P2/ENHs vectors, normalized to Renilla luciferase activity for transfection efficiency control and to empty vector activity. All data are expressed as mean values ± SEM. N = 5. For luciferase experiments statistical significance was calculated for each construct containing an ENH fragment compared to the construct containing the P2 promoter alone. ENH3: TPC1 P = 0.02, BCPAP P = 0.04, MDA-MB231 P = 0.01, MCF7 P = 0.006; ENH11.2: TPC1 P = 0.04, BCPAP P = 0.0008, MDA-MB231 P = 0.03, MCF7 P = 0.004; ENH13: TPC1 P = 0.0003, BCPAP P = 0.003, MDA-MB231 P = 0.008, MCF7 P = 0.03.

Since our goal was to identify regulatory regions driving RUNX2 expression in thyroid and breast cancer, we first confirmed that these regions had features of active ENHs in these cells. We performed ChIP experiments using chromatin markers in two thyroid cancer cell lines (TPC1 and BCPAP) and two breast cancer cell lines (MDA-MB231 and MCF7) expressing different levels of RUNX2 (Figure 1B). In MDA-MB231 and MCF7 cells, as we previously reported for TPC1 and BCPAP, the RUNX2 isoform I transcribed from the P2 promoter was the major RUNX2 isoform, while the expression of Isoform II (transcribed from the P1 promoter), was not detectable (Supplementary Figure S2A). All tested ENHs, with the exception of ENH6 and ENH14, showed a significant enrichment of H3K4Me1 as compared with the P2 promoter (Figure 1D), while H3K4Me3 was significantly enriched on the P2 promoter and barely present on the ENHs (Figure 1C). ChIP analysis of H3K27Ac profiles showed that ENH3, ENH10, ENH11 and ENH13 were significantly enriched in this modification as compared with the other ENHs (Figure 1E). In MDA-MB231 ENHs 4, 8 and 9 were also enriched in H3K27Ac, while in MCF7 ENH10 was the most enriched in H3K27Ac. Annotation data for osteoblasts showed an additional intragenic H3K27Ac peak (between exons 3 and 4). We did not observe H3K27Ac enrichment in this region in thyroid and breast cancer cell lines indicating that this may represent an osteoblast-specific RUNX2 regulatory element (Supplementary Figure S1B and C). Next, we evaluated RNAPolII recruitment on RUNX2 locus by ChIP analysis (Figure 1E). A significant enrichment of RNAPolII binding on ENH3, ENH10, ENH11 and ENH13 was observed in all cell lines tested, reflecting the H3K27Ac binding profile in these regions. RNA-PolII binding was also visible on ENH12 in all cell lines. ENH3, ENH10, ENH11, ENH13 showed higher levels of RNAPolII. In particular, RNA-PolII binding levels in ENH11 were similar to that detected on the promoter. Noticeably, the extent of RNAPolII recruitment on P2 promoter and on active ENHs correlated with the level of RUNX2 expression. Next, we assessed the ability of ENH3–14 to transactivate the RUNX2-P2 promoter in a luciferase assay (Supplementary Figure S2B).

As shown in Figure 1F and Supplementary Figure S2C-F, ENH3, ENH11 and ENH13 induced a significant activation of the P2 promoter in all cell lines even if the activity of ENH11 and ENH13 was significantly stronger than that of ENH3. The effect of ENH13 was stronger in TPC1 and MDA-MB231 than in BCPAP and MCF7 in accordance with the higher RUNX2 expression in these cell lines. By contrast ENH11 activity was significantly higher in TPC1 as compared to the rest of cell lines analyzed suggesting a possible cell-type specific activity of these elements.

Taken together, these data indicate that ENH3, ENH11 and ENH13 are transcriptionally active and suggest the possibility that these elements are involved in the regulation of RUNX2 in thyroid and breast cancer.

ENH3, ENH11 and ENH13 interact with the P2 promoter and are required for RUNX2 expression

Three-dimensional chromatin organization is fundamental for gene expression regulation and is necessary for the functional interaction of promoters and ENHs. To assess the physical association between RUNX2-P2 promoter and the predicted ENHs we performed Chromosome Conformation Capture (3C) experiments.

We queried the interaction profile of RUNX2-P2 promoter with sequences flanking ∼95 kb 5′ and 330 kb 3′ in TPC1 cells, spanning from P1 promoter to ENH14. As shown in Figure 2A, the RUNX2-P2 anchor displayed high-interaction frequency with fragments encompassing ENH3, ENH11 and ENH13, but not with the other ENHs tested. An additional region upstream of ENH6 showed interaction with P2-anchor, but this did not reach significance. This region does not display features of active ENH. Recently, it has been shown that RUNX2-P2 promoter is in the same topologically associating domains (TADs) with RUNX2-P1 promoter and with the promoter of the syntenic SUPT3H gene in osteoblast precursors cells (38). In our model, we observed a significant interaction between RUNX2-P2 anchor and the SUPT3H promoter region, but not with RUNX2-P1 promoter, suggesting a different structural organization of the RUNX2 genomic locus in the tested cancer cell model. To confirm the relevance of ENH3, ENH11 and ENH13 in RUNX2 regulation we used a CRISPR-Cas9 genome editing approach to delete ENH3, ENH11 and ENH13 in BCPAP cells.

Figure 2.

Cooperation between ENH3, ENH11 and ENH13 is required for RUNX2 expression. (A) Chromosome Conformation Capture (3C) of RUNX2 genomic locus. The positions of the RUNX2 P1, P2 and SUPT3H promoters and of the active ENHs are displayed. The interaction frequency for each fragment is expressed as a fraction of the fragment nearest to the anchor. N = 3. ENH3 P = 0.02, ENH10 P = 0.0007, ENH11 P = 0.029, ENH13 P = 0.014. Peak corresponding to region upstream of ENH6 P = 0.06737. (B) Schematic representation of the CRISPR/Cas9 approach employed to generate deletions of the three active ENHs. For each ENH, positions of sgRNAs and extent of deletion are indicated. (C and D) qRT-PCR analysis of RUNX2 (C) and SUPT3H (D) expression in 20 control clones, 9 ENH3 deleted clones, 7 ENH11 deleted clones and 9 ENH13 deleted clones. All data are expressed as mean values ± SEM. * P < 0.05.

We designed two single guide RNA (sgRNA) flanking each ENH (Figure 2B). The sgRNA guides were transfected into BCPAP cells overexpressing Cas9 under a doxycycline inducible promoter. Cas9 expression resulted in a site-specific deletion of the fragments defined by the sgRNAs (922 bp for ENH3, 1905 bp for ENH11 and 1501 bp for ENH13). We derived single clones for each sgRNA couple and for control (no sgRNAs). Copy number for each ENH in each clone was determined by digital PCR (Supplementary Figure S3A and B). Figure 2C shows the levels of RUNX2 expression in ENHs-deleted and CTRL clones. Noticeably, deletion of each one of the three ENHs resulted in a significant down regulation of RUNX2 expression. By contrast, no relevant effect was observed in the expression of the syntenic gene SUPT3H following ENHs deletion (Figure 2D). We observed no differences in RUNX2-Isoform II (P1-driven) expression which is unaltered and barely detectable in all clones tested (data not shown). These observations indicate that ENH3, ENH11 and ENH13 are specific and non-redundant regulatory elements for RUNX2 expression.

ENH3, ENH11 and ENH13 are active in vivo in human thyroid cancer patients

We investigated the activity of ENH3, ENH11 and ENH13 in human normal and thyroid cancer samples.

In these samples, RUNX2 expression was strongly induced in thyroid tumor as compared with matched normal tissue, in accordance with our previous work (19) (Figure 3A). To establish whether ENH3, ENH11 and ENH13 are activated during thyroid cells oncogenic transformation, we used ChIP analysis to evaluate the levels of H3K27Ac in these regions. As shown in Figure 3B, the levels of H3K27Ac in ENH3, ENH11 and ENH13 were higher in cancer samples as compared to normal thyroid tissues. This observation is in accordance with the increased transcriptional activity of the P2 promoter (Figure 3B) and with the higher RUNX2 expression (Figure 3A) observed in the same samples and confirms that changes in the epigenetic status of ENH3, ENH11 and ENH13 correlates with the activation of RUNX2 expression in vivo.

Figure 3.

ENH3, ENH11 and ENH13 are active in human cancer patients. (A) qRT-PCR analysis of RUNX2 expression in normal thyroid and papillary thyroid carcinoma (PTC) tissue of human patients. (B) ChIP analysis of RUNX2 P2 promoter and ENH regions with anti-H3K27Ac antibodies in PTC patients normal or tumoral tissue (N = 14). C) Average number of mutations per patient in RUNX2 promoter and active ENHs in a cohort of 26 papillary thyroid carcinomas that developed distant metastasis (DM-PTCs) and 28 control PTCs that did not develop metastasis (CTRL-PTCs). All the identified mutations were single nucleotide substitutions. Data are expressed as mean values ± SEM. * P < 0.05. (D) Distribution of mutations in the analyzed regions (P2, ENH3, ENH11 and ENH13) in all PTC samples and in the two subgroups (DM-PTC and CTRL-PTC).

To obtain an indirect measure of RUNX2 ENHs activity in vivo, we investigated whether mutational load in RUNX2 functional regulatory elements correlates with cancer aggressiveness. It has been recently demonstrated that, in cancer cells, active ENHs are characterized by a higher mutation frequency as compared to non-active elements. This is likely determined by the fact that active ENHs are generally not accessible to DNA repair machinery since stably occupied by transcription factors. Consequently, active ENHs accumulate mutations faster than other non-transcriptionally active regions (39). We analyzed by NGS the mutational profile of P2, ENH3, ENH11 and ENH13 in a unique cohort of 26 highly aggressive papillary thyroid carcinomas that developed distant metastasis (DM-PTCs). As control, we analyzed the mutational profile of the same regions in a set of 28 control PTCs (CTRL-PTCs) that did not develop metastasis with a minimum follow-up of 8 years. Noticeably, the average number of mutations per patient in the entire panel was significantly higher in DM-PTCs than in CTRL-PTCs (Figure 3C). This was a specific feature of RUNX2 regulatory elements and not just the consequence of a higher mutational burden of DM-PTCs as compared with Ctrl-PTCs. Indeed, NGS analysis in a 174-amplicon panel corresponding to 26 coding genes, on the same cohort of tumors, did not show significant difference in the overall mutation frequency between DM-PTCs and Ctrl-PTCs (data not shown). Furthermore, mutations distribution within RUNX2 ENHs showed that P2 promoter and ENH13 were the regions with the highest incidence of mutations, while ENH11 and ENH13 were the regions with the highest difference in the rate of mutations between DM-PTCs and CTRL-PTCs (Figure 3D). Likely, these differences reflect the actual activation status of these elements in thyroid cancer. All detected mutations were private nucleotide substitutions. This type of genetic mutations may alter ENHs activity by introducing new binding sites for TFs. However, the mutations heterogeneity detected within RUNX2 ENHs suggests that this is not a driving mechanism for the aberrant activation of RUNX2 ENHs in thyroid cancer.

BRD4 controls RUNX2 expression in cancer by binding to ENH3, ENH11 and ENH13

Cooperation between distal ENHs and target promoters requires the binding of multiple proteins that modify the tridimensional structure of chromatin bringing together DNA regions far apart in the linear sequence of the genome. Among these proteins, BRD4 member of the Bromodomain and Extra-Terminal (BET) motif family, plays a primary role in the regulation of distal ENHs, in particular in cancer where its activity is required to sustain the expression of key oncogenes (e.g c-MYC) (40–42). Analysis of publically available BRD4 tracks in MCF7 and in two models of triple negative breast cancer cell lines (HCC1395, SUM159) confirms that BRD4 is enriched within RUNX2 putative regulatory regions and that its binding is higher in triple negative breast cancer cell as compared to MCF7 (Supplementary Figure S4A). We analyzed BRD4 expression in thyroid and breast cancer cell lines to establish potential correlation with RUNX2 levels (Figure 4A). Western blot analysis showed higher BRD4 expression in MCF7 cell line compared to TPC1, BCPAP and MDA-MB231 that showed comparable levels. BRD4 protein expression seems to be inversely correlated with RUNX2 levels since MCF7 are the cells with the lowest RUNX2 expression among the panel tested. ChIP analysis of BRD4 occupancy on RUNX2 locus showed enrichment on RUNX2-P2 promoter and ENH3, ENH11 and ENH13 in all tested cell lines (Figure 4B). Noticeably, BRD4 accumulation was mild on the rest of predicted but non active ENHs (Supplementary Figure S4B).

Figure 4.

RUNX2 is a target of BRD4 and of BRD4 inhibitor JQ1. (A) Western blot analysis of BRD4 expression in TPC1, BCPAP, MDA-MB231 and MCF7 cell lines. (B) ChIP analysis of RUNX2 P2 promoter and active ENH regions with anti-BRD4 antibodies in TPC1, BCPAP, MDA-MB231 and MCF7 cells. An unrelated DNA region upstream of the P2 promoter was used as a negative control. The bars represent the average enrichment of the indicated genomic regions in the immunoprecipitated DNA expressed as percentage of the input. (C) Proliferation curve of TPC1, BCPAP, MDA-MB231 and MCF7 cells treated with the indicated concentrations of JQ1 or mock (DMSO). D) qRT-PCR analysis of RUNX2 in TPC1, BCPAP, MDA-MB231 and MCF7 cells treated with the indicated concentrations of JQ1 or mock (DMSO). (E and F) qRT-PCR and Western blot analysis of BRD4 expression in the indicated cells lines treated with siRNA against BRD4 or control siRNA. G) qRT-PCR analysis of RUNX2 expression in the indicated cell lines treated with siRNA against BRD4 or control siRNA. All data are expressed as mean values ± SEM. * P < 0.05. N = 2.

To investigate the functional relevance of BRD4 binding in the regulation of RUNX2 expression we inhibited BRD4 using the pan-BET protein inhibitor (BETi) JQ1.

We tested the sensitivity of TPC1, BCPAP, MDA-MB-231 and MCF7 cell lines to JQ1 (Figure 4C, Supplementary Figure S4C). TPC1, BCPAP and MDA-MB231 cells showed a medium sensitivity to JQ1 (IC50 250–500 nM) while MCF7, which express the lowest RUNX2 levels, were also the less sensitive to JQ1 (IC50 > 10μM) among the tested cell lines, as already reported (43,44). Next, we analyzed RUNX2 expression 6 and 24 h after JQ1 treatment (Figure 4D). As control, we analyzed c-MYC expression, the best characterized BETi target in cancer and a mediator of the anti-proliferative effect of the drug (Supplementary Figure S4D). Noticeably, RUNX2 expression was significantly down-regulated in all cell lines upon JQ1 treatment at IC50 compatible dosage. Down-regulation was already detectable after 6 h, suggesting that RUNX2 is directly affected by JQ1 inhibition. The strength of repression seemed to be dose dependent at least in thyroid cancer cells. By contrast, c-MYC expression was not significantly affected, suggesting that the anti-proliferative effect of JQ1 in these cancer cell lines mainly relies on c-MYC independent mechanisms. To confirm that BRD4 participate to RUNX2 regulation, we silenced BRD4 in thyroid and breast cancer cells by siRNA (Figure 4E–G). BRD4 silencing resulted in a significant repression of RUNX2 expression in TPC1, BCPAP and MDA-MB231. The effect is less evident in MCF7 cells, in accordance with their lower RUNX2 expression and reduced sensitivity to JQ1 mediated BRD4 inhibition.

To get a deeper insight into the role of BRD4 in regulating RUNX2 expression, we investigated by ChIP analysis, the effect of BRD4 inhibition on the chromatin features of RUNX2 locus (Figure 5A–D). JQ1 treatment induced a marked reduction of BRD4 binding to all active RUNX2 regulatory elements in TPC1, BCPAP and MDA-MB23. In MCF7, JQ1 treatment induced BRD4 binding inhibition only on ENH13, in line with the low sensitivity to this drug. Noticeably, BRD4 inhibition also affects the interaction of the transcriptional co-activator MED1 to RUNX2 ENHs (Figure 5B). Indeed, JQ1 treatment induced a significant decrease of MED1 binding to ENH11 and ENH13 in TPC1 and MDA-MB-231, which are the cell lines with the highest RUNX2 expression. A significant reduction in MED1 binding to ENH11 was also observed in BCPAP, while in MCF7 JQ1 treatment led to a partial inhibition of MED1 recruitment on ENH3 and ENH13.

Figure 5.

BRD4 inhibition modifies chromatin status of RUNX2 locus. (A–D) ChIP analysis of RUNX2 P2 promoter and active ENH regions with anti-BRD4 (D), anti- MED1 (E), anti-RNA-PolII (F) or anti-H3K27Ac (G) antibodies in TPC1, BCPAP, MDA-MB231 and MCF7 cells treated with 1 uM JQ1 or mock (DMSO). An unrelated DNA region upstream of the P2 promoter was used as a negative control. The bars represent the average enrichment of the indicated genomic regions in the immunoprecipitated DNA expressed as percentage of the input. All data are expressed as mean values ± SEM. * P < 0.05. N = 3.

To define how BRD4 affects RUNX2 transcription, we performed ChIP analysis of RNAPolII upon JQ1 treatment (Figure 5C). In thyroid cancer cell lines (TPC1 and BCPAP) BRD4 inhibition led to a dramatic drop of RNAPolII binding on P2 promoter, that likely reflects a decreased initiation complex formation. On the contrary, in breast cancer models (MDA-MB-231 and MCF7), JQ1 treatment determined a strong accumulation of RNAPolII on P2 promoter, suggesting a block of transcription elongation. As previously described (45), H3K27Ac levels were only marginally affected by BRD4 inhibition (Figure 5D).

Taken together, these observations demonstrate that, as for other key pro-oncogenic transcription factors, RUNX2 expression in cancer is controlled by BRD4 and provides a proof-of-principle for targeting RUNX2 expression in cancer using BETi.

Distinctive signaling pathways potentially converge on ENH11 and ENH13 regulation

Binding of specific TFs tunes the activation of ENHs. To map TFs binding to ENH11 and ENH13, we generated sequential deletion mutants that were tested in luciferase assays for their ability to transactivate the P2 promoter in thyroid and breast cancer cell lines. We identified fragments of 232 bp for ENH11 (ENH11.2A1) and 258 bp for ENH13 (ENH13B2) that retained a transactivation activity comparable to the entire ENH sequence in all tested cell lines (Figure 6A–F). Further deletions completely abolished the activity of these elements (data not shown). These fragments likely represent the minimal functional core of these ENHs. Analysis of the ENCODE annotated chromatin features showed that both ENH11.2A1 and ENH13B2 corresponded to regions devoid of H3K27Ac and with a high density of TFs, as previously reported for ENH functional cores. (Supplementary Figure S5A). Next, we used the TRANSFAC database to predict which transcription factors were likely to bind ENH11 and ENH13 cores. This search retrieved 58 putative binding sites for 37 TFs within ENH11.2A1 and 54 putative binding sites for 37 TFs within ENH13B2 (Supplementary Table SVI). To establish whether these TFs were expressed in our cell lines we searched the expression data from the Cancer Cell Line Encyclopedia (CCLE) project. This analysis showed that members of these families of TFs are expressed in at least three of the analyzed cell lines (MDA-MB231, MCF7 and BCPAP while TPC1 are not present in this database) (Supplementary Table SVII).

Figure 6.

ENH11 and ENH13 respond to distinct signaling pathways. (A and B) Schematic representation of ENH11 (A) and ENH13 (B) genomic loci, showing H3K27Ac binding profile and deletion mutants used in luciferase assay to identify the core of each ENH. The fragments retaining transcriptional activity are displayed in black. (C and F) Luciferase analysis of subsequent rounds of deletion mutants of ENH11 (C-E) and ENH13 (D-F) in TPC1, BCPAP, MDA-MB231 and MCF7 cells. Cells were transfected with the indicated pGL3 constructs. The bars represent the average fold change of luciferase activity in cells transfected with pGL3-P2 or pGL3-P2/ENHs vectors normalized to Renilla luciferase activity for transfection efficiency control, and to empty vector activity. Data are expressed as mean values ± SEM. * P < 0.05. N = 4. (G and H) Workflow diagram of the analysis for the identification of the pathways converging on ENH11(G) and ENH13 (H) cores. (I and J) Protein–protein interaction analysis of the relationships among the TFs predicted to bind ENH11 (I) and ENH13 (J). Pink lines represent experimentally validated interactions, blue lines indicate protein homology.

We performed enrichment analysis on the list of predicted TFs, using the Ingenuity Pathway Analysis (IPA) software.

Figure 6G and H report the top scoring canonical pathways for ENH11 and ENH13, obtained by this analysis. With the exception of TGFβ mediated signaling, which is predicted to impact on both, ENH11 and ENH13 are under well distinct regulatory networks. ENH11 was primarily associated with cancer related signaling pathways, while ENH13 was related with development-related signaling pathways. This is interesting since cancer and embryonic development are the two main biological processes in which RUNX2 is engaged and during which its expression must be precisely controlled. Furthermore, the fact that signaling pathways converging on ENH13 are related to developmental control suggest a possible role of this element during embryogenesis.

c-JUN is the master regulator of ENH3, ENH11 and ENH13

To define the potential relationships among the TFs predicted to bind each ENH, we performed protein–protein interaction analysis. This analysis revealed a network of interactions among the TFs predicted to bind to ENH11 (Figure 6I). By contrast, very few connections were detected among the TFs predicted to bind to ENH13 (Figure 6J). In both cases, the heterodimeric complex AP1 resulted as a central node of the interaction map. Two and three AP1 binding sites were predicted within ENH11 and ENH13 respectively (Figure 7A). Furthermore, we recently showed that c-JUN (subunit of AP1) binds to ENH3 and controls RUNX2 expression in cancer cells (20). Thus, we hypothesized a model in which AP1 is central for the activity of the identified ENHs and for their cooperation in the regulation of the P2 promoter.

Figure 7.

c-JUN is the master regulator of RUNX2 ENHs. (A) Sequence and localization of the predicted TFs binding sites on ENH11 and ENH13 core fragments. Consensus diagram for AP-1, TEF-1 and RUNX2 are reported. (B) Luciferase analysis of the effect of targeted point mutations of the indicated binding sites in the ENH11 and ENH13 core fragments in TPC1, BCPAP, MDA-MB231 and MCF7 cells. Cells were transfected with the indicated pGL3 constructs. The bars represent the average luciferase activity of each mutant construct expressed as a fraction of the WT fragment. Values have been normalized to Renilla luciferase activity as transfection efficiency control. All data are expressed as mean values ± SEM. * P < 0.05. N = 3. (C) qRT-PCR analysis of RUNX2 expression in TPC1, BCPAP, MDA-MB231 and MCF7 cells transfected with an empty vector (NT) or dominant negative c-JUN mutant (c-JUN DN). In each cell line, c-JUN DN expression has been verified by immunoblot with anti-c-JUN antibodies. (D) Luciferase analysis of the P2 promoter, ENH11 and ENH13 cores in TPC1 cells cotransfected with increasing amount of c-JUN DN vector. For each sample c-JUN DN expression has been verified by immunoblot with anti-c-JUN antibodies. Values are representative of two independent experiments. (E–H) ChIP analysis of RUNX2 P2 promoter and active ENH regions with anti-c-JUN antibodies in TPC1 (E), BCPAP (F), MDA-MB231 (G) and MCF7 (H) cells. An unrelated DNA region upstream of the P2 promoter was used as a negative control. The bars represent the average enrichment of the indicated genomic regions in the immunoprecipitated DNA expressed as percentage of the input. All data are expressed as mean values ± SEM. * P < 0.05. N = 3. (I and J) qRT-PCR analysis of c-JUN and TEAD1 expression in normal thyroid and papillary thyroid carcinoma (PTC) tissue of human patients. N = 12. (K) Western Blot analysis for c-JUN in a subset of normal thyroid and PTC samples analyzed by qRT-PCR. Quantification indicates the levels of c-JUN normalized on the Actin levels in normal vs tumor tissues.

To test this hypothesis, we first investigated the relevance of the predicted AP1 binding sites in ENH11 and ENH13. To this end we introduced point mutations in the AP1 sites within ENH11.2 A1 and ENH13B2 and tested their functional effects by luciferase assay (Figure 7B). Noticeably, mutation of the adjacent AP1 sites within ENH11 significantly impaired the activity of ENH11.2A1 in both thyroid and breast cancer cell lines. Similarly, mutation of the central AP1 binding site (AP1-II) within ENH13 was sufficient to induce profound repression of ENH13B2 fragment activity. To further confirm these data, we overexpressed in thyroid and breast cancer cells a mutant form of c-JUN that, lacking the transactivation domain, acts as dominant negative (DN) and inhibitor of the c-JUN transcription activity (46). Over-expression of c-JUN DN strongly impairs RUNX2 expression (Figure 7C, Supplementary Figure S5B) in TPC1, BCPAP, MDA-MB231 and MCF7. Furthermore, luciferase analysis demonstrated that c-JUN DN induces repression of RUNX2 ENHs, confirming the centrality of AP1 in the regulation of these elements (Figure 7D). ChIP confirmed that c-JUN is recruited to ENH3, ENH11 and ENH13 in TPC1, BCPAP and MDA-MB231 while its binding on the RUNX2 P2 promoter is limited (Figures 7E–H). In MCF7, c-JUN is enriched exclusively on ENH3. These results are confirmed by the analysis of c-JUN ChIP-Seq track available for MDA-MB231 (33) that shows the significant enrichment of c-JUN on ENH3, ENH11 and ENH13 (Supplementary Figure S5C). TFs rarely act alone, often working in a combinatorial manner to control regulatory elements. Thus, we searched for potential AP1 partners in the regulation of ENH11 and ENH13. Recently, it has been reported that ENHs that are under BRD4 regulation are characterized by the enrichment of consensus sequences for the AP1 family and for the Hippo/YAP-regulated Transcriptional enhancer factor TEF1 (TEAD1) (41). Noticeably, TEF1 is among the TFs predicted to bind both ENH11 and ENH13. We also noticed that ENH13 contains a putative RUNX2 binding site. The ability of RUNX2 to control its own transcription has already been described (47). We used site direct mutagenesis to abolish TEF1 and RUNX2 binding sites within ENH11.2A1 and ENH13B2 and we used luciferase assay to measure the effect of these alterations on the activity of RUNX2 ENHs alone or in combination with AP1 mutations.

Figure 7B shows the results of these experiments. Mutation of the TEF1 binding site within ENH11.2A1 strongly repressed ENH11 and further abolished the transcriptional activity of the AP1 mutant, demonstrating a functional cooperation between TEF1 and AP1 in the regulation of this element. By contrast, mutation of TEF1 binding site within ENH13B2 had no effect on its activity. Mutation of the RUNX2 binding site inhibited the activity of ENH13B2 and further limited the activity of the AP1 mutated ENH13B2, suggesting a cooperation between AP1 and RUNX2 in the regulation of this element. In conclusion, these data confirm that AP1 is a master regulator of RUNX2 and suggest that it orchestrates the complex network of regulatory elements required to sustain RUNX2 expression in thyroid and breast cancer. The activation of transcriptional regulatory elements in cancer is primarily driven by the aberrant up-regulation of TFs. Thus, we analyzed the expression of c-JUN and TEAD1 (member of the TEF family) in normal and thyroid tumor samples (Figure 7I–K). We observed a slight but not significant trend of c-JUN upregulation in tumor samples as compared to normal thyroid both at mRNA and protein level (Figure 7I–K). Conversely, we observed a robust upregulation of TEAD1 mRNA expression in tumor samples, comparable to the one observed for RUNX2 (Figure 7J).

Based on this observation, we may hypothesize that upregulation of c-JUN co-activators in thyroid cancer increases c-JUN recruitment to RUNX2 regulatory elements driving its expression. Otherwise it is possible that post-transcriptional modifications increase c-JUN transcriptional activity in thyroid cancer leading to RUNX2 up-regulation. Further analyses are needed to improve the characterization of the molecular mechanisms driving activation of RUNX2 ENHs in cancer.

RAIN is a novel family of RUNX2 Associated Intergenic Long Non-coding RNAs

Active ENHs are often transcribed into short RNAs (eRNA) or represent transcription starting site for the expression of long non coding RNAs (lncRNA) that, in turn, cooperate to the regulation of gene expression (48). We searched publicly available RNA-Seq annotation data to map uncharacterized transcripts within the RUNX2 locus. We identified a large set of intergenic lncRNAs downstream of RUNX2 3′-UTR and we assessed their expression by qRT-PCR (Supplementary Figure S5D). Only the TCONS00011820 was expressed in thyroid and breast cancer cell lines. The annotated sequence of this transcript was located downstream ENH11 and comprised only two short exons. We used a 5′- and 3′-end Rapid Amplification of cDNA ends (RACE) approach to map the full-length sequence of this transcript in TPC1 cells. We characterized a family of lncRNAs that we named RAIN (RUNX2 Associated Intergenic Long Non-coding RNAs) (Figure 8A). This family comprises four major isoforms characterized by an iper-variable common central region and two alternative 5′- and 3′- ends. The common central region is characterized by a wide sequence variability, while the two alternative 3′-ends only differ in length, with the long RAINs isoforms, exceeding of 3010 bp the short one. Noticeably, the two alternative 5′-ends identify different transcription starting sites located within ENH10 and ENH11 respectively. Analysis of RAINs expression in thyroid and breast cancer cell lines demonstrated that the levels of these lncRNAs strongly correlate with RUNX2 expression (Figure 8B). Together with the fact that RAINs are transcribed from RUNX2 functional ENHs, this observation indicates that RUNX2 and RAINs are co-regulated in cancer cells. To further explore this relationship, we analyzed RAINs expression in CTRL and ENHs- deleted CRISPR/Cas9 clones previously described (Figure 8C). As expected, deletion of ENH11, containing one of the RAINs transcription starting sites, strongly reduced RAINs expression. Surprisingly, deletion of ENH3 and ENH13 also strongly impairs the expression of these lncRNAs as already observed for RUNX2. We also demonstrated that c-JUN silencing or inhibition by overexpression of the c-JUN-DN significantly repress RAINs expression (Figure 8D and E). Finally, we observed that JQ1 treatment profoundly inhibits RAINs expression with a stronger effect than the one observed on RUNX2 expression (Figure 8F). This is coherent with the observation that BRD4 controls the transcription of non-coding RNAs from regulatory ENHs (49,50). These data describe a new family of RUNX2 associated intergenic lncRNAs and envisage a functional relationship between RUNX2 and its cognate lncRNAs in cancer.

Figure 8.

Identification of a novel family of RUNX2 Associated Intergenic Long Non coding RNAs (RAIN). (A) Schematic representation of RAIN isoforms and of the genomic locus from which they are transcribed. Arrows indicate the positions of specific primers used in the RACE experiment. Positions of ENH10 and ENH11 are indicated above the diagram. The variable region of RAINs is indicated in gray. (B) qRT-PCR analysis of the expression of RUNX2 and RAIN transcripts in TPC1, BCPAP, MDA-MB231 and MCF7 cell lines. (C) qRT-PCR analysis of the expression of RAIN in 20 control clones, 9 ENH3 deleted clones, 7 ENH11 deleted clones and 9 ENH13 deleted clones. (D) qRT-PCR analysis of c-JUN, RUNX2 and RAIN expression in TPC1 cells transfected with siRNA oligos against c-JUN or control siRNA. (E) qRT-PCR analysis of the expression of RAIN in TPC1 cells transfected with empty vector (NT) or increasing amounts of a plasmid expressing c-JUN DN. (F) qRT-PCR analysis of RAIN expression in TPC1, BCPAP, MDA-MB231 and MCF7 cells treated with 1 μM JQ1 or mock (DMSO). All expression analysis on RAIN transcripts have been conducted with a common primer pair recognizing all isoforms. All data are expressed as mean values ± SEM. * P < 0.05. N = 3. (G) Schematic representation of RUNX2 ENHs cooperation in controlling RUNX2 P2 promoter.

DISCUSSION

Tumor development and progression are the consequence of the deregulation of many pathways, which are orchestrated by the aberrant function of transcription factors. RUNX2 is re-activated in cancer cells and its expression has a negative prognostic value in several cancer types including breast, thyroid, prostate, colon, and lung carcinoma (17,19,21,51–57). In spite of its fundamental role, the mechanisms that regulate the expression of RUNX2 are still largely unknown.

In this work, we provide evidence that RUNX2 expression is regulated by a network of non-redundant ENHs that cooperate with RUNX2 P2 promoter through the selective binding of cancer associated TFs and through chromatin topological conformation (Figure 7G). Other ENHs were previously described to cooperate with the P1 promoter in osteoblasts (58).

Together with RUNX1 and RUNX3, RUNX2 belongs to the mammalian RUNT-related transcription factor family. All three RUNX proteins serve as master regulators of development and are aberrant reactivated in cancer. Even if their expression is restricted to specific tissues, all three RUNX genes conserve a dual P1/P2 promoter gene organization, suggesting that the localization of other regulatory elements may be conserved among these three genes. Each RUNX gene is followed by a CLIC gene in an inverted orientation (Supplementary Figure S5E–G). In these regions, that could resemble the localization of ENH11 and ENH13, several H3K27Ac peaks are visible, suggesting that also RUNX1 and RUNX3 may be regulated by the crosstalk between their promoters and distal ENHs.

The discovery of super-ENHs demonstrated that ENHs may cooperate to form higher-order and more efficient regulatory units (11,59). Increasing the density of TF binding sites and chromatin binding proteins within cluster of closely spaced ENHs improves cooperation, leading to maximal transcriptional activity at lower TFs concentration.

We showed that ENH3 ENH11 and ENH13, even if located at great distance on the genome, are brought together by chromatin folding. This structural organization creates a platform for the high density binding of TFs and chromatin remodeling proteins that likely support the high levels of RUNX2 expression required during cancer progression. Our data seem to indicate that tridimensional chromatin looping can organize long range cooperative ENHs within super-ENHs-like structures to potentiate and finely regulate the expression of key genes. In line with our hypothesis, we demonstrated that RUNX2 ENHs (and in particular ENH11 and ENH13) are characterized by the significant accumulation of BRD4 and MED1, which are well known super-ENHs associated proteins (41,59). Furthermore, differently from super-ENHs, where deletion of single elements has been shown to only partially affect the expression of target genes, we observed that CRISPR/Cas9 deletion of just one of RUNX2 ENHs results in a dramatic decrease of RUNX2 mRNA levels (60,61).

ENHs act as the logic gates of many regulatory circuits. We predict that many relevant non-overlapping pathways converge on ENH11 and ENH13 regulation. Cancer related pathways like MAPK, already known to activate RUNX2 expression (62,63), are predicted to affect ENH11. Constitutive activation of these pathways in cancer is usually determined by genetic mutations in components of the signaling cascade. Thus, ENH11 could serve as the final receptor of cell intrinsic de-regulated signaling. By contrast, ENH13 appears to be primarily affected by developmental pathways, confirming that cancer overturns silent embryonic mechanisms to fulfill survival and progression needs.

We demonstrated that c-JUN, member of the AP1 family of TFs, binds to each of the RUNX2 active ENHs and its binding is required for their transcription activity (Figure 7G). Similar to RUNX2, activated c-JUN plays an important role in carcinogenesis and cancer progression (64).

RUNX2 and c-JUN are downstream targets of many oncogenic signals and have been shown to cooperate in the regulation of target gene expression in cancer cells. However, we have been the first to demonstrate the requirement of c-JUN for RUNX2 expression in cancer, unveiling a previously unknown functional hierarchy between these TFs (20). A recent paper demonstrated that, during osteoblast differentiation, BRD4 drives the activation of a specific subset of lineage specific ENHs (41). Analysis of TF predicted binding sites identified an enrichment of consensus sequences for the AP1 and TEF-1 families within the BRD4-depedent lineage specific ENHs (41). Our data, in line with these observations, consolidate the existence of a functional relationship between BRD4 and AP1 in the regulation of cancer genes. Two possible models to explain this cooperation can be imagined. BRD4, binding to highly acetylated elements, functions as docking site for the sequence specific recruitment of c-JUN on ENHs that can be triggered. Viceversa, c-JUN binding is a prerequisite for the accumulation of BRD4 on ENHs favoring their interaction with target promoters. Further experiments are needed to disclose the molecular details of this cooperation and to define whether acetylation of c-JUN can impinge on this mechanism.

One of the most exciting and unexpected findings in the genomic era is the extensive transcription activity of the non-coding genome (1). Here, we report the identification of RAINs, a previously unknown family of RUNX2 associated lncRNA. Like RUNX2, which expression is driven by two context specific alternative promoters, RAINs are transcribed starting from two alternative ENHs. ENHs-associated lncRNAs often act in cis to regulate the expression of neighbor genes. We showed that RAINs expression is strongly associated with RUNX2 expression, suggesting a possible involvement of RAINs in the complex network of events that controls RUNX2 expression in cancer.

In cancer, the effect of BETi is mediated by oncogenes repression (42,65). We demonstrated that RUNX2 expression is dependent on BRD4 activity and can be targeted by BETi in cancer cells. The best characterized target of BETi is c-MYC. Noticeably, all thyroid and breast cancer cell lines tested in this work showed a significant sensitivity to JQ1. However, differently from RUNX2 which was profoundly inhibited by JQ1, c-MYC was not significantly affected. These observations indicate that the anti-proliferative effect of BETi in these tumors is primarily mediated by c-MYC independent mechanisms, among which RUNX2 repression could be of crucial importance. We already demonstrated that targeting RUNX2 expression by epigenetic drugs like HDACi is an effective strategy to restrain tumor growth (20). One of the major effects of HDACi is to promote histone acetylation leading to re-expression of silenced oncogenes. However, this mechanism cannot explain how HDACi promote the simultaneous repression of key oncogenes (66–68). BRD4 binds to active, highly acetylated regulatory elements. It is possible that, increasing acetylation levels across the genome, HDACi promote BRD4 delocalization from active ENHs to other sites reducing the expression of highly expressed genes. Our observation that JQ1 represses RUNX2 by blocking BRD4 recruitment to RUNX2 promoter and ENHs, seems to strengthen and confirm this hypothesis.

Furthermore, based on this model, we may speculate that HDACi and BETi cooperate to the repression of RUNX2 and other cancer driver genes, supporting the rational for a combinatorial use of this drugs in cancer (69).

In conclusion, our data identify a comprehensive model to explain how RUNX2 expression is aberrantly reactivated in cancer and how cancer driving signaling pathways converge on its activation.

Click here for additional data file.