Profiling of histone modifications during transformation

In order to systematically characterize epigenetic changes during transformation, we profiled the four cell lines with available commercial antibodies for histone modifications. Interestingly, we observed a gradual reduction of H3K9me2 and me3 along with transformation (Fig. 2a). On the contrary, H3K9 acetylation showed a gradual elevation (Fig. 2a). No obvious changes were observed for other modifications (Fig. 2a). Next, we studied if our discovery in in vitro cell model is also the same in clinical patient tissues. We performed immunostaining in commercial breast cancer tissue arrays and observed a strong reduction in cancer tissues for both methylations (Fig. 2b). In total, we analyzed 140 cancer and 40 adjacent samples for H3K9me3 and 42 cancer and 66 adjacent samples for H3K9me2. Statistical analyses of fluorescence densities representing H3K9 me2 and me3 show obvious lower methylation levels in cancer tissues (Fig. 2c, ​,dd).

The altered transcription programs of cancer-related genes regulated by H3K9me2

To further confirm our discovery, we performed ChIP-seq studies in the four cell lines of H3K9me2, H3K4me3, and H3K27me3. The data suggested H3K9me2 reduced significantly in HMC-L, -LT, and -LTR; H3K27me3 had no significant change (Fig. 3a, ​,b).b). H3K4me3 enrichment seemed reduced also in transformed cell lines but not significant enough (Additional file 1: Figure S3G). However, we did not see a gradual change by high throughput sequencing among the later three cell lines, which is different from western. The reason might be that two methods used different standards for normalization, histone H3 for western and input for ChIP-seq. Nevertheless, our results indicate that histone H3K9me2 decreased significantly during transformation.

Open in a separate window

Fig. 3

Decreased H3K9me2 alters the transcription program of cancer-related genes at the boundaries of LOCKs. a, b Average enrichment levels of H3K9me2 (a) and H3K27me3 (b) around detected peaks. c The heatmap shows the dynamic H3K9me2 levels on the oncogenes with H3K9me2 reduction. The reads of all H3K9me2 peaks from −10 kb of TSS to TES on each gene were used. The genes verified later by ChIP-PCR and RT-PCR are labelled by arrows. d H3K9me2 levels on several oncogenes were confirmed by ChIP-PCR. e The mRNA expression level of oncogenes with decreased H3K9me2 in transformed cell lines. f The average levels of H3K9me2 LOCKs in the four cell lines. Peaks with length >100 kb were defined as LOCKs. All LOCKs in the four cell lines were extended down- and upstream for 50 kb, and their average H3K9me2 enrichment levels were calculated. g, h The number and average coverage of H3K9me2 LOCK number in four cell lines. i KEGG pathway analysis of genes located at the boundaries of LOCKs with decreased H3K9me2

We then analyzed the dynamic changes of H3K9me2 distribution during tumor transformation. A slight increase of H3K9me2 peaks on distal intergenic regions during transformation was observed (Additional file 1: Figure S3A–C). The length of all H3K9me2 peaks and decreased peaks reduced modestly during the process, while the length of increased peaks remained unchanged (Additional file 1: Figure S3D, E). The biological meanings for these changes still required more studies. Interestingly, we also found that the length of H3K27me3 peaks decreased along with transformation (Additional file 1: Figure S3F), though its total enrichment did not change significantly (Fig. 3b). Various mutations and abnormal expression of EZH2, the major enzyme catalyzing H3K27me3, are associated with cancers [15, 17, 18]. Our discovery provides a clue to the study how H3K27me3 dysregulation is involved in breast cancer transformation.

We noticed that decreased H3K9me2 peaks were mainly located on gene bodies (Additional file 1: Figure S3C), therefore, the relationship between gene expression and H3K9me2 during transformation was investigated. Considering the importance of promoters for gene expression, the total reads of H3K9me2 peaks from −10 kb before TTS (usually considered as proximal promoter region) to TES was calculated. We found that H3K9me2 dynamically changed on many oncogenes (Fig. 3c). We then confirmed the change of histone modifications on some of them (Fig. 3d and Additional file 1: Figure S4A). We further analyzed their messenger RNA (mRNA) levels and confirmed with quantitative PCR that many genes exhibited correlated expression (Fig. 3e and Additional file 1: Figure S4B–C). Meanwhile, we found that a few genes did not show the expected elevated expression, such as MYC and BCL2 (Fig. 3e). H3K9me2 is generally believed to be a transcription repression mark and its reduction is required for gene activation. However, H3K9me2 reduction on genes may not directly lead to transcription alteration, and upstream signals sometimes are required.

We also analyzed H3K4me3 and H3K27me3 on genes with decreased H3K9me2. Surprisingly, the average levels of these marks on the above genes are all lower than that in primary cell (Additional file 1: Figure S4D, E). Recently, several groups reported that synergistic decrease of histone H3K4, K9, and K27 methylation is associated with the increase of DNA methylation [36, 37]. It is possible that some genes with H3K9me2 decrease may be also involved in DNA methylation changes.

Boundaries of H3K9me2 LOCKs are enriched with cancer-related genes

H3K9me2 often modifies broad regions on chromatin, which have been previously named as large organized chromatin K9 modifications (LOCKs) [38]. The H3K9me2 on LOCKs has a similar pattern as the total peaks, high in HMC and low in the other three cell lines (Fig. 3f). The total numbers and genome coverages of H3K9me2 LOCKs decreased during transformation (Fig. 3g, ​,h).h). We analyzed the genes overlapped with all decreased H3K9me2 LOCKs but did not find any relationships with cancer. We speculated that change of a broad chromatin region may start from its boundaries, so we analyzed the genes located on the boundaries of the decreased LOCKs. Surprisingly, we found that the genes located in these boundaries are enriched with cancer-related pathways (Fig. 3i and Additional file 2: Tables S8–S10). This suggests that the localization of genes in H3K9me2 LOCKs is related with cellular functional changes and unlocking the oncogenes at the boundaries facilitates the transformation process. The genes at decreased H3K9me2 LOCK boundaries in three transformed cells are very similar, including BTG1, MYC, PIK3R1, and WNT5A for oncogenes and CDKN2A, CDKN2B, FH, INTS6, LINC00032, MITF, PLK2, and WNT5A for tumor suppressors (Additional file 2: Table S11). MYC is the well-known oncogene in breast cancer and many other cancer types [39]. All the other genes may also play important roles in breast tumor transformation.

High expression of H3K9 demethylase KDM3A/JMJD1A in breast cancer cell lines

After determining the reduction of H3K9 methylation as a frequent event in breast cancer, we started to investigate the underlying molecular mechanisms. We firstly examined the mRNA levels of all the known H3K9 methyltransferases and demethylases by RNA-seq and quantitative RT-PCR, but none matched the observed H3K9 methylation pattern (Additional file 1: Figure S5A, B). We then asked if the regulations of these enzymes take place at the protein level. We surveyed the protein levels in transformed cell line with all the available antibodies and found that KDM3A/JMJD1A, a demethylase for H3K9me1 and me2, gradually increased during transformation, inversely matching the decrease of H3K9me2 (Fig. 4a and Additional file 1: Figure S5C). We further found that KDM3A/JMJD1A is much higher in two breast cancer cell lines, MCF and T47D, than that in the primary HMC and other two cancer cell lines, HCT116 and 769-P (Fig. 4b). A commercial breast tissue array containing 48 pairs was stained with KDM3A and the statistical analysis showed that KDM3A significantly increases in breast cancer tissues compared with normal tissues (Fig. 4c, ​,d).d). Another piece of array from the same batch was stained with H3K9me2. Fifteen of 48 pairs (31.3 %) showed both KDM3A increase and H3K9me2 decrease. Taken together, these data show that histone H3K9 demethylase, KDM3A/JMJD1A, increases in breast cancer cell lines.

Open in a separate window

Fig. 4

Overexpression of KDM3A/JMJD1A in breast cancer cells and tissues. a Analysis of mentioned histone H3K9 methyltransferases and demethylases by Western blotting; KDM3A/JMJD1A levels gradually increase with transformation. b High expression of KDM3A/JMJD1A in breast cancer cell lines MCF7 and T47D. c Tissue array with paired breast cancer and normal tissues was used for KDM3A staining. d The average KDM3A level in breast cancer tissues is significantly higher than that in the paired normal tissues (n=48). e Venn diagrams to show that 15 out of 48 paired samples showed both 1.5-fold KDM3A increase and 1.5-fold H3K9me2 decrease in the cancer tissues

Considering the inconsistency of its mRNA and protein level, KDM3A/JMJD1A is probably regulated at the post-translational level. To further verify it, MG132 (an inhibitor for proteasome) or chloroquine (CQ, an inhibitor for lysosome) was used to treat HMC cell line. Both drugs increased the protein level of KDM3A/JMJD1A (Additional file 1: Figure S5D), suggesting its stability was controlled by both proteasome and lysosome. To verify the function of KDM3A/JMJD1A, we expressed its wild type or catalytic dead mutant (H1180A) and confirmed the expression of wild type decrease H3K9me2 in the cell (Additional file 1: Figure S5E).

H3K9 dimethylation and transcription of cancer-related genes regulated by KDM3A/JMJD1A

To further investigate the role of KDM3A/JMJD1A in regulating transformation, we knocked it down in HMC-LTR using small interfering RNA (siRNA) and found that KDM3A/JMJD1A deficiency rescued the expression of most cancer-related genes in HMC-LTR to the levels in primary cells (Additional file 1: Figure S6A). We performed RNA sequencing and found that KDM3A/JMJD1A deficiency mainly affected the genes involved in cell proliferation, cell cycle, wound healing, and protein transport (Additional file 1: Figure S6B, C). We took the DEGs in HMC-LT (Fig. 5a, Additional file 2: Tables S12 and S13) and HMC-LTR (Fig. 5b, Additional file 2: Tables S14 and S15), respectively, which were rescued by KDM3A/JMJD1A knockdown to the levels in primary cell, and performed GO analysis. The biological processes of these genes are very similar to those changed during transformation (Fig. 5c, ​,d),d), indicating KDM3A/JMJD1A is the key factor for the process. H3K9me2 ChIP-seq analysis in HMC-LT further showed that the increased LOCKs by KDM3A/JMJD1A knockdown largely overlapped with the decreased LOCKs in transformation (Fig. 5e, left). The analysis in HMC-LTR showed similar results (Fig. 5e, right). KDM3A/JMJD1A deficiency also restored a large portion of the decreased H3K9me2 peaks in HMC-LT and -LTR cells (Additional file 1: Figure S6D). The increased genes with KDM3A knockdown were also studied with GO analysis (Additional file 1: Figure S6E).

Open in a separate window

Fig. 5

KDM3A/JMJD1A regulates the transcriptional program of multiple oncogenes during transformation. a, b Expression of DEGs in KDM3A/JMJD1A-deficient cell lines and LT/p6 or LTR/p6 was shown in the heatmap. c, d Biological process analysis of DEGs which were rescued by KDM3A deficiency in transformed cells back to the level in primary cell (marked by arrows in a and b). e Venn diagrams reveal the amount of decreased H3K9me2 LOCKs in LT or LTR restored by KDM3A/JMJD1A knockdown. f KDM3A/JMJD1A knockdown reduced the mRNA level of oncogenes in HMC-LTR. g KDM3A/JMJD1A deficiency increased H3K9me2 level on MYC, BCL2, PAX, PRL, and AGR2 locus. h A cell line stably expressing F-KDM3A/JMJD1A was established from breast cancer cell line T47D. ChIP assay with anti-Flag was performed as indicated. i, j Serum was withdrawn from the medium culturing T47D cells for 24 h, and FBS was added back to induce MYC expression. KDM3A deficiency reduced MYC mRNA and protein induction. *P<0.05, **P<0.01

When KDM3A/JMJD1A was knocked down in HMC-LTR, the mRNA of most of these genes decreased, suggesting KDM3A may directly regulate their transcription (Fig. 5f). The ChIP-seq results were verified with quantitative PCR and we found that KDM3A/JMJD1A regulates H3K9me2 on a group of oncogenes, including MYC, PAX3, AGR2, PRL, and BCL2 (Fig. 5g). We then performed Flag ChIP analysis in an F-KDM3A/JMJD1A stable cell lines derived from a breast cancer cell line T47D. The results indicated that KDM3A/JMJD1A binds directly on oncogenes MYC and PAX3 in T47D (Fig. 5h). We speculated that additional signals are probably required to activate MYC expression. We induced MYC in T47D breast cancer cell line with FBS treatment after serum starvation. The elevation of MYC mRNA was impaired in the absence of KDM3A (Fig. 5i), as well as its protein (Fig. 5j). All these results demonstrated that KDM3A/JMJD1A regulates breast tumor transformation through directly binding MYC and PAX3 oncogenes and modulating their transcription.

KDM3A/JMJD1A deficiency impairs the growth and migration of breast cancer cells

To further explore the function of KDM3A/JMJD1A in breast cancer transformation, we made KDM3A stable knockdown cell lines in HMC-L and HMC-LT. We could not get the cell line in HMC-LTR because it was extremely difficult to have four different constructs integrated into one cell line. The results of MTT assay indicated that the growth of HMC-LT was greatly impaired in the absence of KDM3A (Fig. 6a). The cell cycle analysis also indicated that HMC-LT was arrested at G1 phase with KDM3A deficiency (Fig. 6b). The abilities of cell migration and invasion were measured by RTCA real-time monitor or traditional transwell assay with T47D cell line. The results indicated that KDM3A/JMJD1A overexpression enhanced T47D migration and knockdown repressed it (Fig. 6c and Additional file 1: Figure S7A). But cell invasion was not affected (Add Fig. S7B). Colony formation assay indicated that KDM3A/JMJD1A knockdown significantly decreased the colony number in the plate (Fig. 6d). To further test whether MYC and PAX3A really regulate the tumorigenicity ability of the transformed cell, we knock down the two genes respectively in HMC-LTR. The assays performed by RTCA system indicated that MYC and PAX3 are both involved in cell migration and invasion (Additional file 1: Figure S7C). The above data indicated that KDM3A positively regulates the survival and migration of breast cancer cells.

Open in a separate window

Fig. 6

The regulation of tumorigenesis by KDM3A. a KDM3A stable knockdown cell lines were made in HMC-L and -LT. Cell viability was measured by MTT assay. b KDM3A/JMJD1A was knocked down in HMC-LT, and cell cycle was analyzed by flow cytometry with PI staining. c RTCA cell migration assay was performed as described with KDM3A/JMJD1A overexpression and knockdown cells. d KDM3A/JMJD1A was stably knocked down in T47D breast cancer cell line. Cells were seeded in soft agar and colonies were counted 4 weeks later. Experiments were repeated three times, with two duplicates in each experiment. e The protein level of KDM3A and H3K9me2 in the collected breast cancer and paired adjacent tissues is measured with western. f The mRNA levels of MYC and PAX3 in the 15 paired cancer and adjacent tissues are measured by quantitative RT-PCR. g The level of KDM3A protein, H3K9me2, MYC, and PAX3 mRNA was summarized. *P<0.05, **P<0.01

Transformation of human breast tumor cell

Transformed breast cell lines were generated as previous described. 293FRT cells were co-transfected with packaging plasmid ZV77 (psPAX and pMD2G for lentivirus) and pBabe retroviral plasmids containing desired complementary DNA (cDNA). Supernatants containing virus were harvested 48 h later and HMC was infected together with 8 μg/ml polybrene. Typically, more than 80 % of cells were infected as measured by parallel infections with a GFP-expressing construct. Drug selection was performed with 200 μg/ml G418 for neomycin, 50 μg/ml hygromycin, or 0.5 μg/ml puromycin. pBabe with large T antigen, hTERT, or RAS (V12) were purchased from Addgene.

Immunofluorescent staining of cancer tissues and cultured cells

Tumors tissue array (Alenabio, www.alenabio.com) were fixed in 10 % formalin, embedded in paraffin, and followed by standard dewaxing procedures. Cells were cultured on the cover slips and fixed with freezing methanol after washing twice in PBS. The cover slips or tumors tissue array were then washed three times by PBS and blocked in PBS with 1 % BSA for 10 min or 1 h. The cover slips or tumor tissue arrays were hybridized with first and second antibodies for 1 h, respectively. Then, the slips were mounted with prolong anti-fade kit (Invitrogen) and observed with fluorescent microscopy. The arrays used for staining are as follows: H3K9me2 - BR243K, BR243L, BR243M, BR724, BR725, BCN963a; H3K9me3 - BR243B, BR243D, BR243K, BR243L, BR243A, BC081120, BCN963a; and KDM3A - BR724, BR725.

Reverse transcription and quantitative PCR

Cells were scraped down and collected by centrifugation. Total RNA was extracted with RNA extraction kit (Yuanpinghao) according to manufacturer’s manual. Approximately 1 μg of total RNA was used for reverse transcription with a first-strand cDNA synthesis kit (Toyobo). The amount of mRNA was assayed by quantitative PCR. β-Actin was used to normalize the amount of each sample. Assays were repeated at least three times. Data shown were average values ± SD of one representative experiment. All primer sequences are presented in Additional file 2: Table S16.

ChIP assay

ChIP assay was performed as previously described [40]. Briefly, approximately 1×10⁷ cells were fixed with 1 % formaldehyde and quenched by glycine. The cells were washed three times with PBS and then harvested in ChIP lysis buffer (50 mM Tris-HCl, pH 7.6, 1 mM CaCl₂, 0.2 % Triton X-100). DNA was digested to 150–300 bp by MNase (Sigma) before extensive centrifugation. Four volumes of ChIP dilution buffer (20 mM Tris-HCl, pH 8.0, 150 mM NaCl, 2 mM EDTA, 1 % Triton X-100, 0.1 % SDS) was added to the supernatant. The resulted lysate was then incubated with protein G beads and antibodies at 4 °C over night. The beads were washed five times and DNA was eluted by Chip elution buffer (0.1 M NaHCO₃, 1 % SDS, 20 μg/ml proteinase K). The elution was incubated at 65 °C over night and DNA was extracted with DNA purification kit (TIANGEN). The purified DNA was assayed by quantitative PCR with Biorad MyIQ. Assays were repeated at least three times. Data shown were average values ± SD of representative experiments. The sequences of primers are in Additional file 2: Table S16.

Pipeline of RNA-seq analysis

mRNA-seq library was performed by using Illumina TruSeq library construction kit. A 5 μg of total RNA was used as initiation and then prepared according to the manufacturer’s instruction. mRNA-seq libraries were sequenced using HiSeq2000 for 100-bp paired-end sequencing. Quality control of mRNA-seq data was performed using Fatsqc, and then low quality bases were trimmed. After quality control, data were mapped to hg19 genome reference by Tophat2 and allow maximum 2 mismatch. Cufflinks were used to find out differential expression genes. Gene ontology analysis was performed using DAVID (http://david.abcc.ncifcrf.gov) [44, 45].

Pipeline of ChIP-seq analysis

ChIP was performed using desired antibodies. Library was prepared using Illumina TruSeq kit according to the manufacturer’s procedure. Briefly, DNA was prepared for end repair and “A” tailing, adaptor ligation, and library amplification. ChIP-seq libraries were sequenced on HiSeq 2000 for 100-bp paired-end sequencing.

Quality control of ChIP-seq data was performed using Fastqc, and then low quality bases and adaptor contamination were deleted. After quality control and data filtering, data were mapped to hg19 using BWA aln algorithm. Since H3K9me2 and H3K9me3 appear large scale in chromatin, SICER software was used for peaks calling with window size 1000 and gap size 10,000. H3K9me2 and H3K9me3 enrichment region gene annotation was performed using RefSeq gene reference [46]. Gene ontology analysis and KEGG pathway analysis were performed using DAVID.

TCGA breast cancer differential expression gene analysis

The gene expression data pf 100 paired breast cancer and normal tissues were downloaded from TCGA data portal for analysis of differential expressed genes (https://tcga-data.nci.nih.gov/tcga/, the data of total 101 pairs were downloaded but that of one cancer tissue were not readable, so only 100 pairs were used). In order to figure out significant differential expression genes between cancer sample and normal tissue, genes’ expression level less than 5 FPKM in all 201 samples were deleted, and then ANOVA analysis was performed for the rest genes with P value cutoff 0.001. After ANOVA analysis, genes’ average expression level between cancer and normal tissues less than twofold change was deleted. Gene ontology analysis of differential expression genes was performed using DAVID.

Cell cycle analysis with flow cytometry

Cells were harvested after digestion with 0.05 % Trypsin-EDTA. The cells were then washed twice with PBS and fixed in ice-cold 70 % ethanol overnight. Fixed cells were washed twice with PBS and stained in PBS containing propidium iodide (PI, 50 μg/mL) and RNase (100 μg/mL) for 30 min at 37 °C. Cell cycle analysis was performed on an Epics XL-MCL flow cytometer (Beckman Coulter) with System II (version 3.0) software (Beckman Coulter). Additional analysis of cell cycle distribution was determined using Flowjo software.

Cell viability assay

Cell viability was performed by the MTT assay as previously described [47]. Briefly, cells were split at 1×10³ per well in 96-well plates. Next every 24 h, the cells were added with MTT (0.25 μg) in each well for 4 h at 37 °C; the medium with the formazan sediment was dissolved in 50 % DMF and 30 % SDS (pH 4.7). The absorption was read at 570 nm.

Colony formation assay

The bottom layer of 0.6 % agar noble in medium was first placed onto 6-well plate. Cells were seeded in 0.35 % top agar containing medium. Fresh top agar was added 1.5 weeks later, and colonies were counted 8 weeks later. For HMCs, the 50,000 cells were seeded while 5000 or 10,000 cells were seeded for T47D.

Cell migration and invasion assay

The RTCA assay was done as the manufacturer’s protocol. Cells were cultured at 6000 per well in CIM-Plate wells coated with (invasion) or without (migration) matrigel. The cell index signals were read by xCELLigence RTCA DP Analyzer (ACEA Bioscience Inc.). Invasion and migration are monitored continuously over a 48-h period.

The transwell assay was done as follows. Briefly, cells were split at 1×10⁵ per well in 24-well transwell plates coated with (invasion) or without (migration) matrigel. The cells were fixed in 4 % PFA and stained by crystal violet after 48 h. The positive cells were counted under microscope.

Each experiment was repeated three times and results were presented as mean ± SD.

Cancer tissue collection

All the cancer tissues are collected after obtaining the consents of the patients. All the experiments are carried out in accordance with the approved guidelines and protocols by Medical Ethics Committee of Zhongnan Hospital, Wuhan University.

We thank Dr. Man Mohan (Shanghai Jiao Tong University, China), Dr. Hong-Bing Shu (Wuhan University, China), and Dr. Wei Xie (Tsinghua University, China) for discussion and revision of the manuscript.

This work was supported by grants from the National Basic Research Program of China (973 Program, 2011CB504206 and 2012CB518700) and the National Natural Science Foundation of China to Min Wu (31470771 and 91019013) and Lianyun Li (31221061, 31200653, and 31370866).