Identification of epigenetic interactions between miRNA and DNA methylation associated with gene expression as potential prognostic markers in bladder cancer

Identification of interactions between miRNA and methylation associated mRNA levels

Epigenetic interactions associated with survival outcome

To determine whether any of the 227 significant combinations has a significant effect on survival of the patients, an overall survival analysis was performed. Each pair of miRNA and methylation probe from 227 combinations was scatter-plotted, then patients were divided into nine subgroups based on three quantiles of each data set. For example, an interaction between cg13979562 methylation probe and hsa-miR-107 miRNA was significantly associated with gene expression level of NACC2. Out of nine subgroups based on three quantiles of levels of methylation probe intensity at cg13979562 and expression of hsa-miR-107, two groups were chosen to be tested whether survival outcome was significantly different: Group 1 (L/L) having low levels of methylation and low levels of miRNA, and Group 2 (H/H) having high levels of methylation and high levels of miRNA (Fig. 1 (a)). Survival analysis showed that L/L (N = 54) and H/H (N = 58) epigenetic subgroups acting on NACC2 significantly differ in prognosis (p = 3.21e-02), meaning that patients with higher methylation probe intensity at cg13979562 and higher level of hsa-miR-107 expression were associated with a good prognosis (Fig. 1 (b)). Among 227 combinations of gene, miRNA and methylation, 11 combinations acting upon ten genes were found to be significantly associated with survival (p < 0.05): E2F3, CCND1, UTP6, CDADC1, SL35E3, METRNL, TPCN2, NACC2, VGLL4, and PTEN (Table 4 and Additional file 1: Figure S2). Methylation (on two different probes) of METRNL along with the presence of miR-107 was found to be a significant epigenetic interaction impacting survival outcome. To assess whether these interactions are more significant than expected by chance, methylation and miRNA were randomly selected 1,000 times per gene to generate background data set, then, LRT was performed. These 11 epigenetic interactions were still more significant than expected by chance (Additional file 1: Table S2).

Gene	Methylation	miRNA	Bonferroni-corrected LRT p-value	p-value from survival analysis
E2F3	cg21803390	hsa-miR-217	8.31E-08	2.17E-02
CCND1	cg03040489	hsa-miR-944	3.81E-05	1.08E-02
UTP6	cg13453082	hsa-miR-1254	2.44E-04	4.54E-02
CDADC1	cg17226947	hsa-miR-107	7.37E-04	1.03E-02
SLC35E3	cg02006977	hsa-miR-940	1.08E-03	3.58E-02
METRNL	cg01502876	hsa-miR-107	1.49E-03	1.50E-02
TPCN2	cg10490196	hsa-miR-1976	1.69E-03	3.56E-02
NACC2	cg13979562	hsa-miR-107	1.92E-03	3.21E-02
VGLL4	cg25619837	hsa-miR-3662	9.72E-03	1.89E-02
METRNL	cg03155999	hsa-miR-107	1.15E-02	1.19E-02
PTEN	cg166686761	hsa-miR-543	4.97E-02	4.53E-02

Differential gene expression between subgroups defined by epigenetic interactions

To test for differential levels of gene expression between subgroups defined by epigenetic interactions, a T-test was performed for the 11 combinations that were significantly associated with survival (Additional file 1: Figure S3). For example, Fig. 1c shows the boxplot of gene expression levels between the low epigenetic control cohort (L/L) and the high epigenetic control cohort (H/H) for NACC2 at cg13979562 and hsa-miR-107 (p = 4.49E-12). Epigenetic control led to lower expression levels of E2F3, CCND1, CDADC1, SLC35E3, NACC2, and VGLL4, and the epigenetic control of four genes, CCND1, CDADC1, NACC2, and VGLL4, was found to be associated with a worse survival outcome (Additional file 1: Figure S2 and S3). For UTP6, METRNL, TPCN2, and PTEN, high epigenetic control was associated with higher or unchanged expression levels of the target genes (Additional file 1: Figure S3). Of these, high epigenetic control led to a worse survival outcome for UTP6 and PTEN (Additional file 1: Figure S2). In addition, we tested the significance of differential expression levels among four different subgroups, including low methylation & low miRNA (LL), low methylation & high miRNA (LH), high methylation & low miRNA (HL), and high methylation & high miRNA (HH) levels. Levels of gene expression of CDADC1, E2F3, METRNL, NACC2, VGLL4 across four subtypes were significantly different (p < 0.05) and there were consistent patterns that levels of gene expression decreased gradually based on the levels of epigenetic controls, whereas the gene expression levels of PTEN across four subtypes showed the opposite pattern (Additional file 1: Figure S4). Notably, even though NACC2 mRNA levels were significantly lower in HH, HL, and LH subgroups than in LL subgroups, HH subgroup was only associated with good prognosis, indicating that this signal was achieved by only interactions between miRNA and methylation, not by miRNA or methylation alone (Fig. 2). Other results of overall survival analysis for comparing four subgroups can be found here (Additional file 1: Figure S5).

Discussion

DNA methylation downregulates gene expression by inhibiting binding of transcription factors to DNA. Downstream, miRNA silencing, mediated by RNA polymerase II, also works to downregulate gene expression by regulating the processing of mRNA transcripts. We expected miRNA activity targeting methylated genes to decrease expression and worsen survival outcome. We found that a common interaction was one in which there was, as expected, lower gene expression with the presence of epigenetic control, and this phenomenon would tend to lead to a worse survival outcome. Canonical CGI methylation is associated with gene expression silencing, but our results seem to support the previous finding that cancer cells seem to activate CGI methylation of hypomethylated genes which were previously lowly expressed in normal tissues; hypermethylation did not increase the expression of the corresponding genes in cancer cells, but transcription factors were overexpressed [27]. In some cases, the target gene maintained same or higher expression even with microRNA and methylation (UTP6, METRNL, TPCN2, and PTEN); for all these except TPCN2, this led to a worse survival outcome (Fig. 3). Out of the 120 genes we found ten genes to be associated with the interaction between miRNA silencing and methylation, the following epigenetic interactions were shown to be especially associated with survival outcome in bladder cancer even though expression patterns largely did not show epigenetically-induced simple downregulation.

Differential miRNA expression patterns in cancer – tumor-suppressing effect

Decreased expression of miR-217 was found to be significantly associated with large tumor size and advanced clinical stage [28], and miR-217 directly suppressed E2F3 and thereby inhibited invasion of hepatocellular carcinoma [29]. MiR-217 was also found to regulate and be regulated by miR-30a-3p, which suppresses p53 [30]. Previously, introduction of synthetic miR-107 suppressed growth of human non-small cell lung cancer cell lines [31] and high levels of miR-107 were associated with a better survival outcome in gastric cancer [32]. MiR-107 was found to target DICER1 and thereby regulate tumor invasion and metastasis (Fig. 3) [33]. Mutations in DICER1 lead to an abnormally short Dicer protein that is unable to aid in the production of miRNA; Dicer acts as an oncogene or tumor suppressor in varying contexts, including varied roles in bladder cancer (Fig. 3) [34]. MIR-940 levels were found to be the highest in invasive and advanced bladder cancer [35] and has previously been found to inhibit the migratory and invasive potential of cells and increase E-cadherin expression by regulating MIEN1. MiR-940 is highly expressed in immortalized normal cells compared to cancer cells and plays a role in mesenchymal-to-epithelial transition (MET) [36]. MiR-543 is known to target SIRT1 in gastric cancer [37], and miR-543-mediated targeting of SIRT1 is known to alleviate insulin resistance [38]. MMP7 (an oncogene) is also targeted by miR-543 in ovarian cancer; downregulation of miR-543 promotes cancer invasion [39]. The function of miR-1976 is poorly characterized although it was identified as aberrantly expressed in lymphoblastic leukemia [40].

Differential miRNA expression patterns in cancer – high levels of expression in cancer

MiR-944 is overexpressed in human cervical cancer cells [41]. MiR-944 is located in the intron of the TP63 gene but has its own promoter; however, miR-944 biogenesis is markedly increased by the binding of a TP63 gene product, ΔNp63 protein. Moreover, miR-944 upregulates p53 expression [42]. In making the case for distinct subtypes of bladder cancer, basal and luminal, Choi et al. found that TP63 knockdown (and inferred from that, lessened miR-944) deceased basal pathway gene expression and also increased PPAR pathway gene expression, associated with luminal-type carcinomas. The evidence suggests a complex pathway for the interaction of miR-944 and p63 in encouraging the development of primary basal MIBC and perhaps discouraging the luminal type of MIBC. Plasma-borne levels of miR-944 and miR-3662 has been suggested as possible biomarkers for lung cancer [43]. In hepatocellular carcinoma, miR-3662 was also found to be upregulated by p53 (Fig. 3) [44]. MiR-1254 was suggested as a serum-based miRNA biomarker for early-stage lung cancer [45] and in contrast to canonical miRNAs, its biogenesis is independent of DROSHA [46]; miR-1254 expression enhancement may re-sensitize tamoxifen-resistant breast cancer cells to tamoxifen [47].

Downregulation of expression of target with high epigenetic control

The pattern of epigenetic control leading to lower expression held for E2F3, CCND1, CDADC1, SLC35E3, NACC2, and VGLL4, and for some of these (CCND1, CDADC1, NACC2, VGLL4), the epigenetic control was found to be associated with a worse survival outcome. CDADC1 is a domain of cytidine and dCMP deaminase, also called NYD-SP15. CDADC1 was found to dynamically shuttle between nucleus and cytoplasm and overexpression of CDADC1 was found to reduce cell growth and block G1 to S phase transition in the cell cycle [48]. It was also dysregulated in liver cancer [49] and involved in regulating testicular development [50]. High epigenetic control of CDADC1 was associated with lower expression and with a worse survival outcome; we posit that lower expression of CDADC1 would encourage cell cycle deregulation (including less breakdown of cytidine, a component of RNA) and thereby worsen survival outcome. NACC2, also known as RBB, is a transcription repressor that is an important regulator of the p53 pathway: NACC2 inhibits the expression of MDM2, which stabilizes p53 expression [51]. VGLL4 is a Hippo pathway member and acts as a YAP agonist [52]; it is said to function as a tumor suppressor in gastric cancer [53], lung cancer [54] and was also included in a smoking cessation quit-success genotype score calculation [55]. As esophageal squamous cell carcinoma progresses, VGLL4 expression is downregulated [56].

CCND1 gene amplification was previously found to be correlated to histopathological tumor characteristics, cancer-specific survival and response to chemotherapy in bladder cancer [57, 58]. CCND1 protein regulates the cell cycle during the G(1)/S transition, is a substrate for SMAD1, and also phosphorylates and inhibits members of the RB protein family, including RB1. When RB1 is phosphorylated, E2F, a transcription factor, disassociates from the RB/E2F complex and allows the E2F target genes to be transcribed. Our model also implicated epigenetic control of a component of E2F, E2F3, as associated with a better survival outcome. MiRNA targeting (miR-577) of E2F3 was found to inhibit gastric cancer cell progression [59], and miRNA targeting of E2F3 and CCND1 (miR-449b) was found to inhibit the proliferation of SW1116 colon cancer stem cells [60].

In a study of metastatic urothelial carcinoma, significant copy number amplifications were found in E2F3 (30% vs. 7% amplification) and CCND1 [61]. For E2F3 and CCND1, we found that the interaction of methylation of the target gene (E2F3/CCND1) and presence of a targeting miRNA was associated with lower expression of the target gene; high epigenetic control was associated with a worse survival outcome for CCND1, but a better survival outcome for E2F3.

For SLC35E3, the lowered expression resulting from the interaction between miRNA targeting and methylation of target, as for E2F3, led to a better survival outcome. SLC35E3, also called Bladder Cancer-Overexpressed gene 1 (BLOV1), acts in transmembrane transport and belongs to the drug/metabolite transporter protein superfamily (Luscombe, MJ, A novel gene which is overexpressed in advanced bladder cancer May 1999, EMBL/GenBank/DDBJ databases). SLC35E3 was also bioinformatically identified as a potential secreted or transmembrane protein [62].

Maintenance or increase of expression of target with high epigenetic control

For UTP6, METRNL, TPCN2, and PTEN, high epigenetic control was associated with higher expression of the target gene. For UTP6 and PTEN, high epigenetic control led to a worse survival outcome. UTP6 (HCA66) is required for both centriole duplication and ribosome synthesis [63]. Haploinsufficiency-derived UTP6 under-expression in neurofibromatosis type 1 (NF1) cells resulted in those cells being less susceptible to apoptosis [64]. UTP6 was also implicated as element of the interactome of the human histone deacetylase family [65]. PTEN has been shown to be hypermethylated in ovarian cancer cell lines and also highly regulated at the translational level [66], and aberrantly expressed in many forms of cancer [67]. In MIBC squamous epithelia, the expression of PTEN was reduced or lost, and mTOR expression was negatively correlated with PTEN expression only in urothelial squamous cell carcinoma, not schistosomal bladder squamous cell carcinoma [68].

METRNL expression was reduced in the high methylation/high miR-107 (H/H) cohort as compared to the low methylation/low miR-107 (L/L) group for one methylation probe (cg01502876) and unchanged for the other methylation probe (cg03155999), suggesting that methylation of METRNL in cancer will reduce its expression in some cases (Fig. 3). We found that the cohort of patients with high levels of miR-107 and high methylation of its target METRNL at both probes had a better survival outcome. METRNL promotes glucose tolerance and energy expenditure, and blocking METRNL in vivo causes reduces thermogenic gene response and may play a role in tissue inflammation [69]. TPCN2 encodes an NAADP (nicotinic acid adenine dinucleotide phosphate)-induced two-pore Ca(2+) ion channel which is ubiquitously expressed but has elevated expression in liver and kidney and operates as a sensor of both luminal pH and Ca(2+) [70]. TPCN2 is highly sensitive to pH, and is localized intracellularly to endolysosomal organelles [71]. The fact that TPCN2 is induced by nicotinic acid raises questions about potentially novel mechanisms that may underlie well-established link between bladder cancer and smoking.