Exploration of bladder cancer-associated methylated miRNAs by methylated DNA immunoprecipitation sequencing.

BACKGROUND: The current study aimed to explore the association between two epigenomic components, miRNA and DNA methylation, in bladder cancer (BC).
METHODS: Eight paired samples of tumor tissue and matched adjacent normal tissues from BC patients were subjected to methylated DNA immunoprecipitation sequencing and sRNA-Seq for differentially methylated miRNA genes and differential miRNA analysis. The miRNAs regulated by DNA methylation were screened and their functions involved in BC were analyzed using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) as well as a miRNA-mRNA interaction network.
RESULTS: The methylation levels of 212 genes were different between tumors and normal tissues with specific enrichment at transcription initiation and termination sites. Among these genes, 154 were hypermethylated and 58 were hypomethylated. GO and KEGG pathway enrichment analysis indicated that differentially methylated miRNA genes were mainly enriched in tumor-associated GO terms and signaling pathways. Pairwise statistical analysis of MeDIP-Seq and sRNA-Seq data showed that there are 154 and 165 candidate methylation-regulated genes in tumors and normal tissues, respectively. Notably, an interaction network indicated that the miRNAs regulated by methylation regulated a broad range of mRNAs associated with cancer development and progression. In particular, the most differentially expressed miRNAs were validated by qRT-PCR, such that miR-145-5p was downregulated and miR-182-5p was upregulated in patients with bladder cancer.
CONCLUSION: A large number of miRNA genes were modified by methylation in BC. Identification of changes in the expression of these miRNAs provides a great deal of important information for BC diagnosis.

Introduction

Bladder cancer (BC) is a highly prevalent disease of the urinary system in populations worldwide and is ranked as the fifth most common cancer in Western countries.1 In China, the incidence of BC is moderate though it increased year by year between 1998 and 2008 with an average annual growth rate of 4.6%.2 There are multiple risk factors that contribute to the etiology of BC, but smoking and urogenital schistosomiasis are believed to be two of the major risk factors.3–5 Recently studies have indicated that epigenetic changes including DNA methylation, histone modification, and microRNA (miRNA) regulation are associated with an increased risk of BC.6

Increasing evidence has demonstrated that aberrations in methylation affect tumor development and progression in BC.7–9 Methylation of cytosine is a type of genetic annotation which regulates how the genetic code is read and when gene information is transcribed.10 Methylation occurring at CpG sites generally leads to gene suppression or even silencing. Specifically, CpG hypermethylation in promoter regions is associated with a reduction in mRNA levels while CpG methylation in gene body regions often activates transcription.9 Changes in DNA methylation resulting in the silencing of tumor suppressor genes can promote cancer onset and development.11

miRNA is a type of endogenous non-coding RNA that regulates gene expression by inducing direct mRNA degradation or suppressing translation, resulting in abnormal cell behaviors such as cell cycle arrest, apoptosis, proliferation, metastasis, and drug resistance in several cancers inluding BC.12–17 Evidence has been proven that epigenetically regulated miRNAs are another important component of the epigenome and are widely distributed in multiple types of cancers exerting effects on tumor development and progression.18,19 Treatment with methylation inhibitors such as aza-2ʹ-deoxycytidine and 4-phenylbutyric acid can activate the expression of miRNAs and enhance their function in tumorigenesis.20,21

The expression of miRNA can be altered if the methylation occurs in a miRNA-encoding gene.22 Several miRNAs regulated by methylation have been identified in BC and suggested to be associated with the onset and progression of the cancer.14,23,24 However, there are still few comprehensively analyses of the relationship between miRNA expression and DNA methylation in BC. In the present study, we screened interactions between these two epigenomic components by high-throughput deep sequencing and methylated DNA immunoprecipitation sequencing (MeDIP-Seq) with eight paired samples of BC tissue and matched normal adjacent tissues. Further, we suggest a potential regulatory network of methylated miRNA.

Materials and methods

MeDIP-Seq

Following the manufacturer’s instructions, genomic DNA from eight tumors and matched normal adjacent tissues were isolated using a MagMedIP Kit (Diagenode, Denville, NJ, USA). Next, the methylated DNA was separated from the unmethylated fragments by immunoprecipitation using a monoclonal antibody against 5-methylcytidine (Diagenode, Denville, NJ, USA) as previously described.25 Illumina sequencing libraries of the enriched methylated DNA were prepared using NEBNext reagents according to the manufacturer’s recommendations (New England Biolabs, MA, USA). After quantitative evaluation of size distribution and concentration using an Agilent Bioanalyzer 2100 (Agilent Technologies, CA, USA), reads of 100 nucleotides (nt) were sequenced on an Illumina HiSeq 2000 platform (Illumina, CA, USA). All of the samples were obtained with written informed consent and analyzed anonymously. This study was approved by the Medical Ethical Committee of Affiliated Haikou Hospital of Xiangya Medical College, Central South University.

Methylation profiling and differentially methylated gene analysis

The MeDIP-Seq raw data from the Illumina HiSeq 2000 platform were processed according to the standard procedure of the BGI Group (Shenzhen, China). Reads were aligned to a human reference genome (hg19) using the Bowtie alignment algorithm with default parameters. Using a MACS software (MACS 1.4.0), a whole-genome methylation peak scan was performed and sequencing results for tumors and normal tissues were converted into peak values. Distribution trends for peak values across gene bodies and regions 2 kb directly upstream and downstream of the gene bodies were deduced. Genes with a fold change ≥1.5 and P<0.05 were defined as differentially methylated. Differentially methylated genes were further subjected to unsupervised hierarchical cluster analysis to obtain an overview of the differentially methylated genes in BC tumors. Gene ontology (GO) annotation and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis function prediction were performed using the open-source R package. P-values after correction, also termed Q-values were taken as the false discovery rate. GO terms or pathways with a value of Q <0.001 were considered significantly enriched.

sRNA-Seq and differentially expressed miRNA screening

Total RNA was extracted from eight tumors and matched adjacent normal tissues using TRiZol Reagent (Thermo Fisher Scientific, MA, USA) according to the manufacturer’s instructions. RNA quality was evaluated using an Agilent Bioanalyzer 2100 (Agilent). A sRNA library was constructed using TruSeq Small RNA Sample Preparation Kits (Illumina) according to the manufacturer’s protocols. The quality of the library was confirmed with an Agilent Bioanalyzer 2100 and sequencing was performed on an Illumina HiSeq2000 platform (Illumina).

For miRNA analysis, low quality, adapter-only, and unclipped reads were discarded. Remaining clean reads 15–30 nt in length were mapped to hg19 using SOAP V2.0. miRNA was identified and annotated using miRBase21 data (http://www.mirbase.org/ftp.shtml). To compare the expression of miRNA between tumors and normal tissues, miRNA reads were normalized to miRNA reads per million reads (RPM) and the average RPM levels of each miRNA were obtained. miRNA with a fold change ≥2 and a P-value <0.01 were considered to be significantly expressed.

miRNA–mRNA interaction network construction

To construct a miRNA–mRNA interaction network differentially expressed miRNAs detected indirectly from the total RNA of the eight paired tissue samples were mapped to the differentially methylated genes obtained by MeDIP-Seq. Overlapping miRNAs were used to generate a miRNA–mRNA interaction network. Candidate miRNA-target mRNA was predicted using miRanda and Target-Scan. An interaction network illustrating the regulatory relationship between the miRNAs and their corresponding mRNAs was drawn using Cytoscape software.

Quantitative real-time polymerase chain reaction (qRT-PCR)

Differentially expressed miRNAs identified by MeDIP-Seq and sRNA-Seq were confirmed by qRT-PCR. Briefly, total RNAs were extracted from tissues and reverse transcribed by an miRNA-specific stem-loop RT primer according to the manufacturer’s protocol. Each cDNA generated was amplified by quantitative PCR using TaqMan® miRNA assays (Thermo Fisher Scientific, MA,USA) using a mature miRNA-specific forward primer and the universal reverse primer. The reactions were incubated at 95 °C for 120 s, followed by 40 cycles of 95 °C for 15 s, and 60 °C for 30 s.

Results

Differentially methylated miRNA genes in BC patients

MeDIP-Seq was carried out to identify differences in the methylation of miRNA genes between BC tumors and normal tissues. As shown in Figure 1A, more than half of the methylated miRNA genes identified were found exclusively in one patient, while only a very small number were found in five of the eight patients studied (Figure 1A). Among these methylated miRNA genes, the number of hypermethylated genes was greater than the number of hypomethylated genes in five out of eight patients, and only three patients had a higher number of hypomethylated miRNA genes (Figure 1B). Analysis of the distribution of DNA methylation across gene bodies and regions 2 kb directly upstream and downstream of the gene bodies reveled two methylation peaks close to the transcription start and termination sites in both tumors and normal tissues (Figure 1C). Most of the changes in methylation between tumors and normal tissues also occurred in these two regions. To gain an overview of the differentially methylated miRNA genes identified, a heatmap was generated and distinct patterns of methylation of miRNA genes were observed in tumors and normal tissues. As shown in Figure 1D, a total of 212 genes were differentially methylated between tumors and normal tissues in BC patients. Of these, 154 genes were hypermethylated and 58 were hypomethylated in the tumor tissues. Higher levels of methylation were observed in tumors tissues than in normal controls. In a volcano plot, hyper- and hypomethylated genes were distributed in quadrants of fold change >1 and P-value <0.05, with larger genes concentrated in the hypermethylated quadrant (Figure 1E).

Figure 1

Determination of DNA methylation by MeDIP-Seq analysis. (A) Distribution of deferentially methylated genes encoding miRNAs in eight patients. (B) Number of hypo- and hypermethylated genes in eight patients. (C) Distribution of reads around gene bodies ±2 kb flanking regions. Upstream 2 kb, region 2 kb upstream of the transcription start site 2K; intragenic, the entire gene from the transcription start site to the end of the transcript; downstream 2 kb, 2 kb downstream of the end of the transcript. (D) Heatmap showing deferentially methylated miRNA-encoding genes. Each gene corresponds to each row. (E) Volcano plot showing the differentially expressed reads in tumors and normal tissues. Red points represents upregulated reads, green points represent downregulated reads, black points are not statistically significant.

Functional prediction for differential methylated genes

Hyper- and hypomethylated genes encoding miRNAs were further subjected to GO and KEGG pathway enrichment analysis for function prediction. KEGG analysis showed that differentially methylated genes were mainly enriched in tumor-associated pathways such as the PI3K-Akt signaling pathway, the MAPK signaling pathway, the Ras signaling pathway, focal adhesion, and miRNAs and proteoglycans associated with cancer (Figure 2A). GO term analysis also determined that most genes were enriched in biological processes, followed by cellular components, with the most enrichment in single-organism cellular processes such as the regulation of nitrogen compound metabolism, cellular and cellular metabolic processes, and intracellular organelles (Figure 2B).

Figure 2

KEGG pathway enrichment and GO analysis of the differentially methylated miRNA genes. (A) Top 20 most enriched signaling pathways. Dot size reflects the number of genes enriched in each signaling pathway. Dot color indicates P-value. (B) GO enrichment analysis of all differentially methylated genes.

Differential expression of miRNAs in BC tumors are associated with methylation of miRNA-encoding genes

To investigate the association between miRNA expression and the methylation of miRNA-encoding genes in tumor tissues, miRNA expression profiles were detected by high throughput sequencing for eight paired tumors and normal tissues. As shown in Figure 3A, cluster analysis revealed distinct profiles of miRNA expression between tumors and matched normal tissues. There were 83 upregulated and 59 downregulated miRNAs in the tumor tissues relative to the normal tissues (Figure 3A). A volcano plot also showed a large number of miRNAs that were up- or downregulated by more than 2-fold in tumor tissues (Figure 3B). Differentially expressed miRNAs were classified into four categories according to the region of methylation in their corresponding genes: both gene body and promoter modified; neither modified; only gene body modified; or only promoter modified. The average expression of miRNAs in each category was calculated. As shown in Figure 3C, there was no significant difference in these four groups of miRNAs between tumors and normal tissues. However, the average expression of unmodified miRNAs was higher, and that of miRNAs with methylation in the promoter region, gene body region, or both was lower (Figure 3C). To determine whether DNA methylation in miRNA-encoding genes was associated with the expression of miRNAs, methylated miRNA genes detected by MeDIP-Seq and miRNAs obtained by sRNA-Seq were merged. A total of 154 and 165 overlapping genes were identified in tumors and normal tissues, respectively (Figure 3D and E).

Figure 3

Differentially expressed miRNAs in BC tumor tissues are associated with methylation of miRNA-encoding genes. (A) Heatmap showing differential miRNA expression profiles in bladder cancer relative to normal tissues. (B) Volcano plot showing the differentially expressed miRNAs in tumor tissues compared with normal tissues. Red points represents up-regulated miRNA, green points represent downregulated miRNAs, black points are not statistically significant. (C) Expression levels of miRNAs corresponding to four genomic regions with different methylation modification levels: both gene body and promoter modified, neither modified, only gene body modified, only promoter modified. (D) Overlap of miRNA-encoding genes identified by MeDIP-Seq and sRNA-Seq in tumor tissues. (E) Overlap of miRNA-encoding genes identified by MeDIP-Seq and sRNA-Seq in normal tissues.

Differentially expressed miRNAs regulated by methylation regulated a large number of mRNAs

Generally, miRNAs regulate gene expression by targeting the 3ʹUTR of mRNA. Changes in the expression of miRNA in tumor tissues may therefore affect a number of mRNAs. To generate a miRNA–mRNA interaction network, candidate mRNA targets of differentially expressed miRNAs that overlapped according to MeDIP-Seq and sRNA-Seq data were predicted. Target prediction revealed that numerous mRNAs were potentially regulated by these miRNAs. A miRNA–mRNA interaction network which contained four subsets and included 12 miRNAs and 278 mRNAs was constructed. The largest regulatory subnetwork was that of miR-145, which included 143 mRNAs, accounting for more than 50% of mRNAs in the entire network (Figure 4).

Figure 4

miRNA–mRNA regulatory network generated using miRNAs of differential methylated genes and mRNAs which containing predicted miRNA targets. Green triangle indicates miRNAs, pink oval represents mRNAs.

To further validate the MeDIP-Seq and sRNA-Seq results, RT-PCR was used to evaluate two up-regulated (miR-1266-5p and miR-182-5p) and three down-regulated (miR-145-5p, miR-126-3p, and miR-152-3p) miRNAs in eight paired tumors and normal tissues. qRT-PCR results showed that miR-1266-5p and miR-182-5p were up-regulated in the BC tissues compared with the normal tissues, whereas miR-145-5p was down-regulated in the tumor tissues (Figure 5), which was consistent with the results of sRNA-Seq. However, miR-126-3p and miR-152-3p were up-regulated in the BC tissues, and qRT-PCR results were inconsistent with the sRNA-Seq results.

Figure 5

Validation of up-regulated [miR-1266-5p (A) and miR-182-5p (B)] and down-regulated [miR-145-5p (C), miR-126-3p (D), and miR-152-3p (E)] miRNAs by qRT-PCR.

Discussion

DNA methylation and miRNAs are two important regulatory factors in human cancers. In the present study, we identified miRNA genes regulated by methylation by combining the profiles of methylated miRNA-encoding genes with the expression profiles of miRNAs in BC patients. Several differentially methylated miRNA-encoding genes were identified and predicted to regulate many miRNAs in BC. Further, we outlined a regulatory network of these differentially methylated miRNAs and highlighted the central role of miR-145 in the methylated miRNA–mRNA interaction network.

Previously, studies have discussed some aspects of the relationship between DNA methylation and miRNA in BC. However, previous studies such as that of Ma et al have mainly focused on the differences in methylation between sensitive and resistant HCC cell lines in response to 5-fluorouracil treatment;19 another study by Shimizu et al identified a set of epigenetically silenced miRNA genes in BC patients, but this was mainly based on a PCR-based method, which may not cover all miRNAs in addition to their counterpart genes information. Therefore, previous works may not have covered all possible methylation changes in miRNA genes.24 Unlike these studies, our current study provides a genome-wide profile of both DNA methylation state and miRNAs expression based on a high-throughput deep sequencing method, which means that the data reflects the global changes in methylation state and miRNA levels in BC tissues.

Previous studies have suggested that methylation levels are higher in gene bodies than in regions flanking the 5ʹ and 3ʹ ends of genes in human cells, bovine placentas, and arabidopsis.26–28 In the current study, partially consistent with previous finding, methylation levels were slightly increased in the gene bodies compared with the flanking regions. However, sharply increased methylation peaks were noted in upstream and downstream regions close to the transcription start and termination sites. These peaks were particular high in the BC tissues relative to normal tissues, suggesting that methylation in these regions may contribute to BC etiology.

Our research shows that methylation does have an effect on miRNA expression levels. Both promoter regions and the gene body regions exhibited the lowest expression of methylated miRNA genes, while unmethylated miRNA had the highest average expression level. Furthermore, integrated analysis using MeDIP-Seq and sRNA-Seq data revealed that DNA methylation in miRNA gene can partially explain the aberrant expression of a proportion of miRNAs in normal and cancer tissues. Previously, several miRNAs have been proven to be regulated by DNA methylation. For example, miR-193, a miRNA which has been shown to be hypermethylated, is downregulated in chemoresistant BC cells.29 Treatment with demethylase resulted in changes in the expression of miR-137, miR-124–2, miR-124–3, and miR-9–3, which are associated with growth and invasion in primary cancers.24 Methylation of the promoters of miR-34a and miR-126 is found in a variety of cancers, including BC.30,31 miR-200b, a miRNA associated with cisplatin treatment sensitivity, is upregulated after treatment with demethylase.32 Meanwhile, in the muscle-invasive bladder tumors, the miR-200b gene is hypermethylated and its expression is specifically silenced.33 Our findings provide more information regarding BC associated methylated miRNAs.

Of these miRNAs, we highlight the tumor suppressor, miR-145, which targets the greatest number of mRNAs in the methylated miRNA–mRNA interaction network produced in the present study. In particular, we confirmed the down-regulation of miR-145 in the bladder tumors with RT-PCR. It has been reported that miR-145 regulates tumor progression and metastasis in a wide range of cancer.34–37 Lower levels of miR-145 induced by DNA methylation have been evidenced to enhance cell behaviors such as proliferation, migration and invasion and could be used as potential prognostic biomarker for cancer therapy.38,39 Removal of the methylation suppresses cancer growth and metastasis.40 In terms of BC, miR-145 has been proven to inhibit BC cell growth, cell invasion and migration, and induce tumor apopotosis.41–44 These findings suggest that miR-145 plays an inhibitory role in tumor development and progress and that DNA methylation occurs with high frequency in cancer. In support of these findings, we observed the downregulation of miR-145 in the present study, and our methylated miRNA–mRNA interaction network shows that miR-145 regulates the biggest subnetwork of cancer-associated mRNAs, including c-Myc, SOX9, SP1, and FSCN1.34,36,45,46 Thus, it will be important for us to identify the role of methylated miR-145 in BC in future work.

Collectively, our data suggest that methylated miRNA genes are involved in the pathogenesis of BC and that drugs that regulate the expression of methylated miRNA genes, especially those related to cancer may be useful in BC therapy. In addition, the differentially methylated miRNAs genes identified in the current study may function as biomarkers for the early detection of BC.