Oncotarget

Research Papers:

Identification of differentially expressed miRNAs and miRNA-targeted genes in bladder cancer

PDF |  HTML  |  Supplementary Files  |  Order a Reprint

Oncotarget. 2018; 9:27656-27666. https://doi.org/10.18632/oncotarget.24441

Metrics: PDF 368 views  |   HTML 715 views  |   ?  

Jong-Young Lee, Seok Joong Yun, Pildu Jeong, Xuan-Mei Piao, Ye-Hwan Kim, Jihye Kim, Sathiyamoorthy Subramaniyam, Young Joon Byun, Ho Won Kang, Sung Phil Seo, Jayoung Kim, Jung Min Kim, Eun Sang Yoo, Isaac Y. Kim, Sung-Kwon Moon, Yung Hyun Choi & Wun-Jae Kim _

Abstract

Jong-Young Lee1,9,10,*, Seok Joong Yun2,*, Pildu Jeong2, Xuan-Mei Piao2, Ye-Hwan Kim2, Jihye Kim9, Sathiyamoorthy Subramaniyam9, Young Joon Byun2, Ho Won Kang2, Sung Phil Seo2, Jayoung Kim3, Jung Min Kim4, Eun Sang Yoo5, Isaac Y. Kim6, Sung-Kwon Moon7, Yung Hyun Choi8 and Wun-Jae Kim2

1Department of Business Data Convergence, Chungbuk National University, Cheongju, Republic of Korea

2Department of Urology, College of Medicine, Chungbuk National University, Cheongju, Republic of Korea

3Department of Surgery, Department of Biomedical Sciences, Cedars-Sinai Medical Center, University of California Los Angeles, Los Angeles, California, USA

4NAR Center, Inc., Daejeon Oriental Hospital of Daejeon University, Daejeon, Republic of Korea

5Department of Urology, Kyungpook National University Hospital, Kyungpook National University School of Medicine, Daegu, Republic of Korea

6Section of Urologic Oncology and Dean and Betty Gallo Prostate Cancer Center, The Cancer Institute of New Jersey and Robert Wood Johnson Medical School, New Brunswick, New Jersey, USA

7Department of Food Science and Technology, Chung-Ang University, Ansung, Republic of Korea

8Department of Biochemistry, College of Oriental Medicine, Dong-Eui University, Busan, Republic of Korea

9Microarray Division, Theragen Etex Bio Institute, Suwon, Republic of Korea

10Microarray Division, SNP Medicine Co., Ltd, Suwon, Republic of Korea

*These authors made an equal contribution to this article

Correspondence to:

Wun-Jae Kim, email: wjkim@chungbuk.ac.kr

Keywords: bladder cancer tumorigenesis; miRNA; mRNA; miRNA-mRNA interaction; transcription factor

Received: January 31, 2017     Accepted: August 28, 2017     Epub: February 07, 2018     Published: June 12, 2018

ABSTRACT

Background: Differentially expressed genes and their post-transcriptional regulator-microRNAs (miRNAs), are potential keys to pioneering cancer diagnosis and treatment. The aim of this study was to investigate how the miRNA-mRNA interactions might affect the tumorigenesis of bladder cancer (BC) and to identify specific miRNA and mRNA genetic markers in the two BC types: non-muscle invasive bladder cancer (NMIBC) and muscle invasive bladder cancer (MIBC).

Results: We identified 227 genes that interacted with 54 miRNAs in NMIBC, and 14 genes that interacted with 10 miRNAs in MIBC. Based on this data, we found extracellular matrix-related genes are highly enriched in NMIBC while genes related to core nuclear division are highly enriched in MIBC. Furthermore, using a transcriptional regulatory element database, we found indirect regulatory transcription factors (TFs) for enriched genes could regulate tumorigenesis with or without miRNAs.

Materials and methods: Tissue samples from 234 patients histologically diagnosed with BC and 83 individuals without BC were analyzed using microarray and next-generation sequencing technology, and we used different cut-offs to identify differentially expressed mRNAs and miRNAs in NMIBC and MIBC. The selected mRNAs and miRNAs were paired using validated target datasets and according to inverse expression relationships. MiRNA interacted genes were compared with the TF-regulating genes in BC. Meanwhile, pathway enrichment analysis was performed to identify the functions of selected miRNAs and genes.

Conclusions: Identification of differential gene expression in specific tumor types could facilitate development of cancer diagnosis and aid in the early detection of BC.


Identification of differentially expressed miRNAs and miRNA-targeted genes in bladder cancer | Lee | Oncotarget

INTRODUCTION

Bladder cancer (BC) is the second most common urological malignancy in the United States. In 2015, approximately 73,510 new cases of BC were diagnosed and 14,680 individuals died of this disease in the U.S. [1]. In South Korea, BC is the second most common genitourinary tumor and the incidence of BC is about five times greater in men than it is in women [2].

Conventional diagnosis of BC classifies the tumor into two groups based on its clinico-pathologic features: non-muscle invasive (NMIBC) and muscle-invasive (MIBC). [3]. At the time of first diagnosis, a majority of cases are classified as NMIBC. However, approximately 20% of BC is confirmed as MIBC in situ, which is the main cause of cancer-specific deaths in BC patients. In general, NMIBC has a much better prognosis than MIBC and survival rates for NMIBC are higher than those for other malignancies. However, even after transurethral resection (TUR) of the primary tumor, the chance of recurrence is common with NMIBC. Approximately 30-50% of patients with NMIBC experience recurrence in their lifetime. NMIBC also runs the risk of possibly progressing into MIBC but the chances are usually low, ~10-20% [4, 5]. Therefore, frequent recurrence and eventual progression to MIBC are challenges for both patients and urologists. Despite technical improvements in the diagnosis of BC, more efficient diagnostic tests that can facilitate early diagnosis and active surveillance are urgently needed. We currently lack an ideal non-invasive detection tool and there is a need for more accurate and predictive biomarkers.

To date, more than 1,100 microRNAs (miRNAs) have been identified in the human genome, far fewer than the total number of mRNAs. However, rather than being strictly one to one, interactions between miRNAs and target genes can be one-to-many or many-to-one, resulting in a very large number of potential regulatory effects. In addition, there is some controversy regarding the specific mode of how miRNA acts [6]. Several studies have shown that miRNAs repress mRNAs by regulating their stability or translation efficiency. Regarding BC, evidence from expression profiling studies using microarray-based approaches has provided abundant information to support the theory of miRNAs. Nevertheless, data about the interactions between miRNAs and mRNAs is inconsistent, and little is known about their expression levels and physiologically relevant functions in cancer cells.

In this study, we sought to identify NMIBC- or MIBC-specific miRNAs, predict their associated genes and miRNA-mRNA interactions based on bioinformatics analyses and validation by next-generation sequencing (NGS), and determine the critical networks and pathways that are dysregulated in BC patients. Our results provide important insight into the mechanisms underlying NMIBC and MIBC in human patients, as well as BC progression.

RESULTS

Identification of differentially expressed genes between BC tissue and normal tissue

Determination of mRNAs was performed as shown in Figure 1. We extracted 522 genes that were differentially expressed in NMIBC versus normal tissue, 105 genes differentially expressed in MIBC versus normal tissue, and 60 genes differentially expressed in NMIBC versus MIBC. All genes were then validated by RNA sequencing from a second cohort. A Venn diagram shows that 402/522 genes in NMIBC and 33/105 genes in MIBC were exclusively dysregulated compared to that of normal tissue, and 11/60 genes were differentially expressed only between NMIBC and MIBC (Figure 2A). To understand the biological functions of these genes, we subjected each of them to network enrichment analysis using GeneMANIA and STRING. The 402/522 differentially expressed genes in the NMIBC were strongly associated with the extracellular matrix (ECM) (false discovery rate (FDR), 1.35E-17) and extracellular structure organization (FDR, 2.98E-13) (Supplementary Table 1). The 33/105 differentially expressed genes in the MIBC were associated with mitosis (FDR, 1.72E-21) and nuclear division (FDR, 1.88E-20) (Supplementary Table 2). The 11/60 genes that were differentially expressed between NMIBC and MIBC were associated with cellular responses to zinc ions (FDR, 4.12E-11) (Supplementary Table 3). Seventy-two genes, whose functions are related to nuclear division and G2/M transition in cell mitosis, showed different expression between BC and normal tissues.

Schematic showing negatively correlated gene-miRNA pairs in NMIBC and MIBC.

Figure 1: Schematic showing negatively correlated gene-miRNA pairs in NMIBC and MIBC. Differentially expressed mRNAs and miRNAs were identified according to -fold changes in expression and p-value criteria. MiRNA-targeted gene pairs were identified using miRNAs and mRNAs differentially expressed in NMIBC and MIBC versus normal tissue according to previously validated informations regarding mRNA-miRNA interactions. NMIBC, non-muscle invasive bladder cancer; MIBC, muscle invasive bladder cancer.

Venn diagram showing mRNAs and miRNAs differentially expressed in BC.

Figure 2: Venn diagram showing mRNAs and miRNAs differentially expressed in BC. (A) mRNA expression. (B) miRNA expression. BC, bladder cancer; NMIBC, non-muscle invasive bladder cancer; MIBC, muscle invasive bladder cancer.

Identification of miRNAs differentially expressed between BC tissue and normal tissue

MiRNAs were determined as shown in Figure 1. We extracted 70 miRNAs differentially expressed in NMIBC versus normal tissue, 25 miRNAs differentially expressed in MIBC versus normal tissue, and 6 miRNAs differentially expressed in NMIBC versus MIBC, all of which were validated by RNA sequencing in a second cohort. The Venn diagram shows that 58 miRNAs were exclusively deregulated in NMIBC versus normal, 17 were specifically deregulated in MIBC versus normal, and 2 miRNAs were differentially expressed only between NMIBC and MIBC (Figure 2B). By comparing our data with 190 miRNAs identified in a previously reported miRNA profiling study of BC [7], we identified 27 novel miRNAs: miR-124, miR-1260a, miR-136, miR-149, miR-191, miR-301b, miR-302a, miR-302c, miR-3065, miR-3195, miR-331, miR-335, miR-34a, miR-363, miR-371a, miR-381, miR-409, miR-4284, miR-4492, miR-4532, miR-4634, miR-483, miR-484, miR-5581, miR-5701, miR-590, and miR-654.

NMIBC-/MIBC-specific gene-miRNA interactions

Next, we attempted to identify gene-miRNA interactions based on the validated mRNAs and miRNAs specific for NMIBC and MIBC. This analysis proceeded by two steps: first, we identified candidate gene target miRNAs based on experimentally validated datasets from the mirTarBase v6 [8] in miRWalk 2.0 [9]. Through this, we obtained 8,542 gene-miRNA interaction pairs for the 522 NMIBC genes derived from the mRNA expression profiles (NMIBC vs normal) and 1,740 gene-miRNA interaction pairs for the 105 MIBC genes (MIBC vs. normal). Next, among the gene-miRNA interaction pairs, only significantly expressed genes and miRNAs that paired with each other were selected to further analyze for differences between normal versus NMIBC and normal versus MIBC. Finally, we obtained 227 genes that paired with 54 miRNAs in NMIBC (Figure 3, Supplementary Table 4), and 14 genes that paired with 10 miRNAs in MIBC (Figure 3, Supplementary Table 5). These miRNAs were compared with190 previously reported miRNAs [10]. This enabled us to identify 14 novel miRNAs: seven (miR-1260a, miR-149, miR-191, miR-335, miR-34a, miR-484, and miR-5701) in NMIBC and seven (miR-124, miR-302c, miR-331, miR-371a, miR-4492, miR-4634, and miR-483) in MIBC (Table 1). To validate biological replication on miRNA-mRNA interacted genes, we were recruited and reanalysis GSE and the cancer genome atlas (TCGA) data with AnaltAnlyzer [11]. Finally, these 14 miRNAs showed negative correlations with 130 genes from differentially expressed analysis that were validated expression in NMIBC using GSE40355. The seven miRNAs showed negative correlations with 10 genes from differentially expressed analysis that were validated its expression in MIBC using GSE40355 [12], GSE66064 [13] and RNA sequences from TCGA (Table 1).

Venn diagram showing paired gene and miRNA in BC.

Figure 3: Venn diagram showing paired gene and miRNA in BC. (A) Genes associated with gene-miRNA interactions. (B) miRNAs associated with gene-miRNA interactions. BC, bladder cancer; NMIBC, non-muscle invasive bladder cancer; MIBC, muscle invasive bladder cancer.

Table 1: Expression profiling of miRNA-mRNA interacted target genes in NMIBC and MIBC

a. Normal versus NMIBC

miRNAs

miRNA Expression

Target genes

Target Gene Expression

GSE40355 Expression

miRNAs

miRNA Expression

Target genes

Target Gene Expression

GSE40355 Expression

miRNAs

miRNA Expression

Target genes

Target Gene Expression

GSE40355 Expression

miR-1260a

Up

RNF150

Down

NA

miR-335

Up

ATP8B2

Down

NA

miR-335

Up

AKAP12

Down

Down

PRNP

Down

Down

CCL19

Down

Down

CD1D

Down

Down

miR-149

Up

CALD1

Down

Down

CD14

Down

Down

CD37

Down

Down

FLNA

Down

Down

CD8A

Down

Down

CDKN1C

Down

Down

GFPT2

Down

Down

CFI

Down

Down

COL21A1

Down

Down

JUN

Down

Down

CGNL1

Down

Down

DLG2

Down

Down

NFIX

Down

NA

COL3A1

Down

Down

F10

Down

Down

NR2F1

Down

Down

CYGB

Down

NA

GNA14

Down

NA

miR-191

Up

EGR1

Down

Down

CYTL1

Down

Down

GPR162

Down

NA

miR-335

Up

ABCA8

Down

Down

CYYR1

Down

Down

GZMK

Down

NA

ADAMTS8

Down

Down

DACT1

Down

Down

ITPR1

Down

Down

ANTXR2

Down

Down

DNASE1L3

Down

Down

KLRB1

Down

Down

AOC3

Down

Down

DOCK11

Down

NA

MAFB

Down

Down

BIN1

Down

NA

FBLN5

Down

NA

MFGE8

Down

Down

C1S

Down

Down

FILIP1L

Down

Down

MXRA8

Down

Down

CDH11

Down

NA

FZD1

Down

NA

NAALAD2

Down

Down

CLIC4

Down

Down

FZD10

Down

Down

NFIL3

Down

Down

COLEC12

Down

NA

GABARAPL1

Down

Down

RECK

Down

Down

COL6A1

Down

Down

GATA5

Down

NA

RERG

Down

Down

CPE

Down

Down

GSN

Down

Down

TLR1

Down

NA

CRISPLD2

Down

Down

HCST

Down

Down

VNN2

Down

Down

DACH1

Down

Down

ITGB2

Down

Down

miR-34a

Up

AXIN2

Down

Down

FBLN2

Down

NA

KIAA1644

Down

Down

AXL

Down

Down

FGL2

Down

NA

LIFR

Down

Down

BCL2

Down

Down

FHL1

Down

Down

LMOD1

Down

Down

CALD1

Down

Down

FOS

Down

NA

NKG7

Down

Down

CSF1R

Down

Down

GNA14

Down

NA

OTUD1

Down

Down

EPB41L2

Down

Down

HLA-DRB5

Down

Down

P2RX1

Down

Down

GAS1

Down

Down

KLF9

Down

Down

PLSCR4

Down

Down

MAP1B

Down

Down

LMO3

Down

Down

PODN

Down

NA

RTN4

Down

NA

LPPR4

Down

Down

PPAP2B

Down

Down

TGFBR2

Down

Up

LTBP4

Down

Down

PRICKLE2

Down

Down

TGM2

Down

Up

MAOB

Down

Down

PROS1

Down

Down

TPM1

Down

Down

MGP

Down

NA

RARRES2

Down

NA

VASN

Down

NA

NFIB

Down

Down

RASD1

Down

NA

ZEB2

Down

NA

OLFML3

Down

Down

RHOB

Down

Down

miR-484

Up

ACTA2

Down

Down

PDE5A

Down

NA

RORB

Down

NA

ATP1A2

Down

Down

PDK4

Down

Down

SCARA5

Down

Down

CRTAP

Down

Down

RNASE4

Down

Down

SCRG1

Down

NA

DSTN

Down

NA

RNF150

Down

NA

SFRP1

Down

Down

FAM129A

Down

NA

ROR2

Down

Down

SLC46A2

Down

Down

FLNA

Down

Down

SERPINF1

Down

Down

SORBS2

Down

Down

KLF2

Down

Down

SERPING1

Down

Down

TNC

Down

NA

TSC22D3

Down

Down

SMOC2

Down

Down

TNC

Down

NA

TUBB2A

Down

Up

SYNPO2

Down

Down

TSC22D3

Down

Down

miR-5701

Up

RUNX1T1

Down

Down

TGFBR2

Down

Up

VASN

Down

NA

ZFP36L1

Down

Down

ZEB2

Down

NA

b. Normal versus MIBC

miRNAs

miRNA Expression

Target genes

Target Gene Expression

GSE40355 Expression

GSE66064 Expression

TCGA Expression

miRNAs

miRNA Expression

Target genes

Target Gene Expression

GSE40355 Expression

GSE66064 Expression

TCGA Expression

miR-124

Down

EZH2

Up

NA

NA

Up

miR-371a

Down

MCM10

Up

Up

Up

Up

KIF2C

Up

NA

Up

Up

MELK

Up

Up

NA

Up

AURKA

Up

Up

NA

Up

miR-4492

Down

CDCA8

Up

Up

Up

Up

KIF20A

Up

NA

NA

Up

miR-4634

Down

POLQ

Up

Up

NA

Up

miR-302c

Down

MELK

Up

Up

NA

Up

miR-483

Up

FGF9

Down

Down

NA

NA

miR-331

Up

NELL2

Down

Down

NA

NA

NA: Not available

NMIBC, non-muscle invasive bladder cancer; MIBC, muscle invasive bladder cancer.

To understand the biological functions of the genes that interacted with miRNAs, selected genes were subjected to network enrichment analysis using GeneMANIA and STRING. We found that 221 genes differentially expressed in NMIBC were strongly associated with the ECM (FDR, 4.60E-18) and extracellular structure organization (FDR, 1.72E-11). (Supplementary Table 6), whereas 8 differentially expressed genes in MIBC were involved in spindle formation (FDR, 7.53E-15) and mitosis (FDR, 2.80E-10) (Supplementary Table 7).

The influence of transcription factors (TFs) on BC tumorigenesis

In addition to miRNAs, gene expression is also regulated by TFs. Therefore, we decided to investigate the relationship between TFs and miRNA-interaction genes using the transcriptional regulatory element database (TRED). Twenty-two genes (ACTA2, BCL2, C1S, CAV1, CD14, COL3A1, COL6A1, CTGF, ESR1, F3, FOS, GJA1, ITGB2, JUN, MGP, PRNP, TGFBI, TPM2, TYMS, UCHL1, VIM, and ZFP36L1) are regulated by 62 TFs in 36 TF families [14]. Most target genes of TFs were down-regulated in BC; an exception was the thymidylate synthetase (TYMS) gene, which is regulated by E2F4, MYC, TP53, USF1, and USF2 regulatory elements (Table 2).

Table 2: A lot of TFs interacted with target genes in this study

Transcription factor

Expression

Target genes

Expression

Transcription factor

Expression

Target genes

Expression

Transcription factor

Expression

Target genes

Expression

ELK1

Up

ACTA2

Down

PPARG

Up

ESR1

Down

HIF1A

Up

ITGB2

Down

ATF1

Up

BCL2

Down

RARA

Up

Down

SPI1

Up

Down

BRCA1

Up

Down

SMAD4

Up

Down

ATF1

Up

JUN

Down

CEBPA

Up

Down

SP1

Up

Down

ATF2

Up

Down

CREB1

Up

Down

SP3

Down

Down

BCL6

Up

Down

EGR1

Down

Down

TFAP2C

Up

Down

CREB1

Up

Down

ETS1

Up

Down

USF1

Up

Down

CREM

Up

Down

MYB

Up

Down

AR

Up

F3

Down

EGR1

Down

Down

MYC

Up

Down

CEBPA

Up

Down

ERG

Up

Down

NFKB1

Up

Down

CEBPB

Up

Down

ESR1

Down

Down

NFKB2

Up

Down

CEBPD

Up

Down

ESR2

Up

Down

PPARG

Up

Down

CREB1

Up

Down

ETS2

Up

Down

RARA

Up

Down

E2F4

Up

Down

JUN

Down

Down

RARB

Down

Down

EGR1

Down

Down

MYC

Up

Down

RARG

Up

Down

HIF1A

Up

Down

NFIC

Up

Down

RELA

Up

Down

JUN

Down

Down

PGR

Up

Down

STAT1

Up

Down

JUND

Up

Down

RARA

Up

Down

STAT3

Up

Down

MYB

Up

Down

RARB

Down

Down

WT1

Up

Down

NFIC

Up

Down

RARG

Up

Down

E2F1

Up

C1S

Down

NFKB1

Up

Down

REL

Up

Down

E2F2

Up+

CAV1

Down

REL

Up

Down

SMAD3

Up

Down

E2F4

Up

Down

RELA

Up

Down

SMAD4

Up

Down

ESR1

Down

CD14

Down

SP3

Down

FOS

Down

SP1

Up

Down

SP1

Up

Down

TFAP2A

Up+

Down

STAT1

Up

Down

CEBPA

Up

Down

ARNT

Up

Down

STAT3

Up

Down

CEBPB

Up

Down

ATF1

Up

Down

JUN

Down

MGP

Down

RELB

Up

Down

ATF2

Up

Down

AR

Up

PRNP

Down

SP1

Up

Down

BCL6

Up

Down

E2F4

Up

Down

SP2

Up

Down

CEBPB

Up

Down

JUN

Down

Down

SP3

Down

Down

CREB1

Up

Down

TFAP2A

Up+

Down

JUN

Down

COL3A1

Down

CREM

Up

Down

SMAD1

Up

TGFBI

Down

NFIC

Up

Down

ESR1

Down

Down

MYC

Up

TPM2

Down

SMAD1

Up

Down

ESR2

Up

Down

E2F4

Up

TYMS

Up

SMAD3

Up

Down

ETS1

Up

Down

MYC

Up

Up

AR

Up

COL6A1

Down

ETS2

Up

Down

USF1

Up

Up

SMAD3

Up

Down

ETV4

Up+

Down

USF2

Up

Up

TFAP2A

Up

Down

FLI1

Up

Down

MYB

Up

UCHL1

Down

SP1

Up

CTGF

Down

JUN

Down

Down

MYC

Up

Down

WT1

Up

Down

NFIC

Up

Down

ERG

Up

VIM

Down

AR

Up

ESR1

Down

NFKB1

Up

Down

ETV4

Up

Down

BRCA1

Up

Down

STAT1

Up

Down

LEF1

Up

Down

BRCA2

Up

Down

STAT3

Up

Down

MYB

Up

Down

ERG

Up

Down

JUN

Down

GJA1

Down

NFKB1

Up

Down

HIF1A

Up

Down

SP1

Up

Down

SP1

Up

Down

JUN

Down

Down

CEBPA

Up

ITGB2

Down

MYC

Up

ZFP36L1

Down

PAX5

Up

Down

CEBPB

Up

Down

TF, transcription factor; BC, bladder cancer

DISCUSSION

In this study, we used two different methods to identify global miRNA-mRNA interactions and their enriched functional networks that are specific to NMIBC and MIBC. Microarrays were used to obtain differential expression profiles, while RNA-Seq was used for experimental validation of both miRNA and mRNA levels. Additionally, those gene expression pattern were validated using another cohort which were previously reported GSE and TCGA [1113]. These results showed that expression pattern of gene-miRNA interacted genes were conserved in NMIBC and MIBC on tumorigenesis. We uncovered detailed patterns of miRNA and gene expression, allowing us to propose molecular markers for specific subtypes of BC. The 402 genes differentially expressed in NMIBC were related to ECM function, platelet activation, and multicellular organismal metabolism. Moreover, we found that the miRNAs interacted with the 221 genes in NMIBC were mainly related to the ECM, organization of extracellular structure, and muscle contraction. In contrast, the 33 genes differentially expressed in MIBC were mainly related to mitosis, antigen processing and presentation via MHC class II, and G2/M transition. In addition, we found that 8 miRNA-interaction genes in MIBC were related to spindle formation, mitosis, and nuclear division.

The results of gene-network analysis of gene-miRNA pairs revealed that ADAMTS8, COL21A1, COL3A1, COL6A1, COL6A2, COL18A1, FBLN2, FBLN5, F3, LTBP4, MFGE8, MGP, PARRES2, SERPINF1, SFRP1, SNCA, TGFBI, TGFB3, TNC, and VEGFA genes, which were down-regulated in NMIBC, are involved in the formation or function of the ECM, including the release of new growth factors and ECM molecules. The ECM serves as a structural scaffold that provides the support necessary to maintain tissue integrity and sustainability [15]. Therefore, ECM regulation is important for the generation of new tissue structures, as well as maintaining the architecture and homeostasis of adult tissues [15]. In addition to our discovery of a potential panel of genes for BC, the genes we found that modulate smooth muscular contraction (ACTA2, ATP1A2, CALD1, CAV1, EDNRA, GJA1, LMOD1, MYH11, MYLK, TPM1, TPM2, and VIM) could represent potential therapeutic targets for the treatment of diseases related to bladder contraction [15].

Cell cycle-related genes (AURKA, CDCA8, KIF20A, and KIF2C) were overexpressed in MIBC compared to NMIBC. Moreover, genes that play a role in cell cycle division are over-expressed in BC when compared with normal tissues. Therefore, since cell division, chromatin replication, and chromosome segregation are key targets for drug development, miRNA-targeted genes involved in nuclear division and G2/M transition may be candidate targets for drugs designed to treat BC [16].

Previously published studies propose 440 up- or down-regulated miRNAs as BC biomarkers, which comprise 190 union miRNAs. Among them, miR-145 was down-regulated in multiple different types of cancer tissue, helping to maintain the differentiation status of smooth muscle cells. This observation was reported consistently by 11 other groups [7], suggesting that miR-145 is likely to be a critical factor in BC tumorigenesis. Finally, our comprehensive analyses identified 14 novel miRNAs specific to BC (miR-124, miR-1260a, miR-149, miR-191, miR-302c, miR-331, miR-335, miR-34a, miR-371a, miR-4492, miR-4634, miR-483, miR-484, and miR-5701).

Members of the miR-183 cluster are located within a 5 kb region on human chromosome 7q32.2, and are transcribed in the same direction (from telomere to centromere) [17]. Previous studies have shown that the miR-183 cluster was abnormally expressed in a variety of tumors and could be directly involved in human cancers [18]. Notably, miR-183-3p, miR-182-5p, and miR-96-5p have oncogenic functions in BC [19, 20], and we found here that these miRNAs were over-expressed in both NMIBC and MIBC. In particular, the miR-183 cluster regulates genes that overlap in both NMIBC and MIBC, i.e., DKK3 (a target of miR-183-5p), FGF9 (a target of miR-182-5p), and monooxygenase, DBH-Like 1 (MOXD1, a target of miR-96-5p), suggesting that these miRNAs play important roles in bladder tumorigenesis. Our present study is also one of the first to reveal that MOXD1 is regulated by miR-96-5p in BC tissue. MOXD1 encodes a member of the copper monooxygenase/dopamine/β-hydroxylase family, which is localized to the endoplasmic reticulum. Members of this protein family are mainly involved in biosynthesis of neurotransmitters and hormones.

Furthermore, several lines of evidence indicate that miR-200b and miR-200c are over-expressed to a greater extent in NMIBC than in MIBC, consistent with previous studies [21]. The miR-200 family is thought to play an essential role in tumor suppression by inhibiting the epithelial-mesenchymal transition (EMT), an initiating step of metastasis [22]; because cells lose adhesion during the EMT, the rate of cell mortality is elevated during this process. In particular, miR-200 promotes the final step of metastasis in which migrating cancer cells undergo EMT during colonization of distant tissues [23]. Interestingly, our data shows that miR-200b and miR-200c were more highly expressed in NMIBC than in MIBC. Moreover, our study is the first to report the BCL2 (a target of miR-200b-3p and miR-200c-3p), CDH11 (a target of miR-200c-3p), FBLN5 (a target of miR-200c-3p), and TIMP2 (a target of miR-200c-3p) interaction pairs in NMIBC (Supplementary Table 4a).

Through analysis of the gene-miRNA pairs after linking the expression profiles of TFs and miRNAs to their target genes, we found that the expression of almost all TFs was moderately or severely repressed while miRNAs were over-expressed in BC tissues versus normal tissues. The exceptions to this general trend were E2F2, TFAP2A, and ETV4; TFs that were over-expressed in BC tissue. These results suggest that these TFs and miRNAs may regulate target genes with oncogenic properties in BC tumorigenesis. While almost all of the target genes in the gene-miRNA pairs showed lower expression, the TYMS gene interestingly showed over-expression, indicating possible oncogenic properties. This suggests that miR-99a-5p, which interacts with TYMS, and TFs E2F4, MYC, USF1, and USF2 regulate the expression of the TYMS gene. Further work is needed to gain a better understanding of how TFs and miRNAs act as both oncogenes and tumor suppressors. Advancements could potentially lead to the development of therapeutic targets from the expression profiling data of TFs, miRNAs, and the target genes in BC.

In summary, our validated data obtained from NGS-based experiments, combined with comprehensive and unbiased computational analyses of published data, identified 14 novel miRNAs in BC. Furthermore, genes differentially expressed in NMIBC were associated with the ECM, muscle contraction, and nuclear division. Gene enrichment analysis of the gene-miRNA pairs also suggest that the organization of the ECM is important in NMIBC, whereas DNA replication and G2/M transition are important in MIBC. This report provides the first evidence for novel gene-miRNA interaction pairs involving 227 genes and 54 miRNAs in NMIBC and 14 genes and 10 miRNAs in MIBC.

The promising findings from this study may lead to potentially novel diagnostic and therapeutic interventions in BC and other cancers.

MATERIALS AND METHODS

Patients

The study cohort included 234 patients histologically diagnosed with bladder urothelial carcinoma and 83 individuals without BC. Tumor levels were assessed according to standard criteria [3]. The BC samples were further subcategorized into two groups, NMIBC and MIBC. NMIBC patients, who were usually treated with TUR, were periodically assessed by cystoscopy and urinary cytology every 3 months for the first 2 years, every 6 months for the next 3 years, and annually thereafter. Patients with MIBC underwent radical cystectomy and complete pelvic lymph node dissection. These patients were also subjected to urinary diversion, including creation of a conduit, continent cutaneous reservoir, and orthotopic ileal neobladder. Moreover, MIBC patients with pT3/pT4/lymph node–positive disease were subjected to four to six cycles of cisplatin-based adjuvant chemotherapy. These MIBC patients were subjected to physical examination, urine cytology, serum chemistry, chest x-ray, and abdominal and pelvic computerized tomography every 3 months for the first 2 years, every 6 months for the following 2 years, and then annually thereafter. Clinical data was assessed for all patients retrospectively. The mean follow-up of these patients was 71 months (median 61 months; range, 15–115 months). Normal muscosae was collected far from tumors and real normal mucosae from benign diseases. All specimens were rapidly frozen in liquid nitrogen and stored at -80°C until use. Collection and analysis of samples was approved by the Institutional Review Board of Chungbuk National University (IRB approval number 2006-01-001 and GR2010-12-010), and informed consent was obtained from each subject.

Datasets

High-throughput molecular datasets were generated for miRNAs using two different technologies. Microarray and RNA sequencing mRNA data were generated by RNA sequencing. The analysis also included microarray datasets published previously by our group (Gene Expression Omnibus accession number GSE13507) [24]. Sample sizes are provided in Figure 1.

Library preparation and data generation

RNA extraction

Total RNA was isolated from the indicated tissues using the TRIzol reagent (Life Technologies, Carlsbad, CA) and purification was performed using phenol based extraction methods according to the manufacturer’s instructions. The RNA concentration was determined using NanoDrop ND-1000 spectrometer and the RNA integrity number was evaluated on a 2100 Bioanalyzer (Agilent, Santa Clara, CA) using the RNA 6000 Nano Kit (Agilent, Santa Clara, CA).

miRNA microarrays

Total RNA (100 ng) from each sample was dephosphorylated, 3' end–labeled with Cy3-pCp, purified on Micro Bio-Spin columns, dried, and hybridized using the miRNA Microarray System labeling kit ( Illumina, San Diego, CA) and the Agilent Human miRNA Microarray Release 16.0 platform, which contains 1,205 human and 144 viral miRNAs [25]. The protocol used to generate microarray gene expression datasets is provided in reference [26].

Stranded mRNA library construction

mRNA sequencing libraries were prepared using the TruSeq Stranded mRNA Sample Preparation Kit (RS-122-2101) (Illumina, San Diego, CA). Oligo dT attached magnetic beads were used to purify poly-A–containing mRNA from 1 μg of total RNA. Next, the purified mRNA was disrupted into short fragments, and first-strand cDNAs were synthesized using SuperScript II reverse transcriptase (Invitrogen, Carlsbad, CA) and random hexamers. cDNA with adapters ligated to both ends were enriched by PCR. cDNA library size and quality were evaluated electrophoretically using the Agilent DNA 1000 Kit (part # 5067-1504) on a 2100 BioAnalyzer. Subsequently, the libraries were sequenced on an Illumina HiSeq 2500. Image analysis was performed using the HiSeq control software version 2.2.58. Raw data were processed and base calling was performed using the standard Illumina pipeline (CASAVA version 1.8.2 and RTA version 1.18.64).

Small RNA library construction

Small RNA sequencing libraries were constructed using the TruSeq Small RNA Sample Preparation protocol (RS-200-0012) (Illumina, San Diego, CA). Illumina adapters were directly and specifically ligated to microRNA molecules with a 3'-hydroxyl group and a 5'-phosphate. The quality and size distribution of the adapter-ligated RNAs and amplified libraries were confirmed using the High Sensitivity DNA Analysis Kit (Cat. #5067-4626) (Agilent, Santa Clara, CA). Libraries were quantitated using the Library Quantification Kit for NGS (KK4824) (Kapa Biosystems, Wilmington, MA). Subsequently, libraries were sequenced on an Illumina HiSeq 2500. Real-time image analysis and base calling were performed on the instrument using the HiSeq Sequencing Control Software version 2.2.58. CASAVA software version 1.8.2 and RTA version 1.18.64 were used for de-multiplexing and generation of FASTQ sequence files.

Data analysis

Microarrays

The Robust Multiarray Average, in the R package [27], was used to perform global correction, quantile normalization, and median Polish summarization. P-values (t test) were calculated from bead mRNA signal intensities [27].

mRNA sequencing

Total sequencing reads were subjected to preprocessing as follows: adapter trimming was performed using cutadapt with default parameters, and quality trimming (Q30) was performed using FastQC with default parameters. Processed reads were mapped to the human reference genome (Ensembl 72 [GRCh37: hg19]) using tophat and cufflink with default parameters [28]. Fragments Per Kilobase of exon per million fragments Mapped (FPKM) values were normalized and quantitated using R package Tag Count Comparison (TCC) [29] to determine statistical significance (e.g., P and Q values) and differential expression (e.g., -fold changes).

miRNA sequencing

Total sequencing reads were subjected to preprocessing: adapter trimming, quality trimming, size selection (17–24 nt), and clustering using mirDeep2 with default parameters [30]. The preprocessed cluster representatives were mapped to the human reference genome (Ensembl 72 [GRCh37; hg19]), and miRNA regions were mined with respect to miRBase reference co-ordinates. Total read counts were subjected to Trimmed mean m-values normalization and quantitation using R package edgeR, and statistical significance (e.g., P and Q values) and differential expression (e.g. -fold changes) were determined.

Integrated analysis of differentially expressed miRNAs and mRNAs

For data obtained using each technology, we used different cut-offs to identify differentially expressed mRNAs and miRNAs. Initially, thresholds for p-values (P ≤ 0.05) and -fold changes (log2FC ≥ 1.0) were applied to the microarray dataset. Thresholds for the RNA-sequencing datasets were as follows: mRNAs (FPKM ≥ 0.3, P ≤ 0.05, FDR ≤ 0.05, and log2FC ≥ 1.0) and miRNAs (read count ≥ 10, P ≤ 0.05, FDR ≤ 0.05, and log2FC ≥ 1.0). Finally, mRNAs common to the microarray dataset and mRNA-Seq datasets were selected as validated mRNAs. MiRNAs common to the microarray dataset and miRNA-Seq datasets were selected as validated miRNAs. The selected mRNAs and miRNAs were paired using validated target datasets and according to inverse expression relationships (i.e., miRNAs whose level increased as that of their target decreased and vice versa).

Analysis of the relation between transcription factors (TFs) and miRNA-interacting genes

TRED was used to identify the TFs related to the miRNA targeting genes in BC. There are curated 36 cancer-related TF families in TRED [14]. And we compared the miRNA-mRNA interacted pairs with the known TF families.

Pathway enrichment analysis

Pathway enrichment analysis was performed by submitting selected miRNAs to DIANA mirPath V3.0 [31]. Significantly correlated and enriched pathways were identified by calculating significant p-values using Fisher’s meta-analysis method for specific miRNAs. Finally, the list of genes was examined using GeneMANIA [32] to identify interactive functional networks from existing databases.

Abbreviations

BC: Bladder cancer; ECM: Extracellular matrix; ETM: Epithelial-mesenchymal transition; FDR: False discovery rate; MIBC: Muscle invasive bladder cancer; miRNA: microRNA; MOXD1: Monooxygenase DBH-like 1; NMIBC: Non-muscle invasive bladder cancer; TCC: Tag count comparison; TCGA: The cancer genome atlas; TF: Transcription factor; TRED: Transcriptional regulatory element database; TUR: Transurethral resection; TYMS: Thymidylate synthetaset.

Author contributions

WJ Kim led the project and the manuscript preparation.

Isaac Y. Kim, SK Moon, YH Choi designed and managed the project and together with WJ Kim.

P Jeong, XM Piao, YH Kim, YJ Byun, HW Kang, SP Seo obtained the list of BC samples to align and performed the alignment.

S Subramaniyam and JM Kim performed the microarray and sequencing data analysis.

JY Lee and SJ Yun performed integration analysis and wrote a first draft of the paper.

J Kim and P Jeong contributed to revisions.

ACKNOWLEDGMENTS

The authors wish to thank Ms. Eun-Ju Shim from the National Biobank of Korea at Chungbuk National University Hospital for the sample preparations and her excellent technical assistance.

CONFLICTS OF INTEREST

The authors declare no conflicts of interest.

FUNDING

This research was supported by the International Science and Business Belt Program through the Ministry of Science, ICT and Future Planning (2017K000490) and Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2017R1A6A3A11028303). The authors also acknowledge support from Department of Defense (W81XWH-15-1-0415) and the U.S.-Egypt Science and Technology Development Fund by the National Academies of Sciences, Engineering, and Medicine (to J.K.).

REFERENCES

1. Siegel R, Naishadham D, Jemal A. Cancer statistics, 2012. CA Cancer J Clin. 2012; 62:10-29.

2. Jung KW, Park S, Won YJ, Kong HJ, Lee JY, Seo HG, Lee JS. Prediction of cancer incidence and mortality in Korea, 2012. Cancer Res Treat. 2012; 44:25-31.

3. Babjuk M, Oosterlinck W, Sylvester R, Kaasinen E, Bohle A, Palou-Redorta J, European Association of Urology (EAU). EAU guidelines on non-muscle-invasive urothelial carcinoma of the bladder. Eur Urol. 2008; 54:303-314.

4. Messing E. Markers of detection. Urol Oncol. 2007; 25:344-347.

5. Quan C, Cha EJ, Lee HL, Han KH, Lee KM, Kim WJ. Enhanced expression of peroxiredoxin I and VI correlates with development, recurrence and progression of human bladder cancer. J Urol. 2006; 175:1512-1516.

6. Wilczynska A, Bushell M. The complexity of miRNA-mediated repression. Cell Death Differ. 2015; 22:22-33.

7. Lee JY, Ryu DS, Kim WJ, Kim SJ. Aberrantly expressed microRNAs in the context of bladder tumorigenesis. Investig Clin Urol. 2016; 57:S52-59.

8. Chou CH, Chang NW, Shrestha S, Hsu SD, Lin YL, Lee WH, Yang CD, Hong HC, Wei TY, Tu SJ, Tsai TR, Ho SY, Jian TY, et al. miRTarBase 2016: updates to the experimentally validated miRNA-target interactions database. Nucleic Acids Res. 2016; 44:D239-247.

9. Dweep H, Gretz N. miRWalk2.0: a comprehensive atlas of microRNA-target interactions. Nat Methods. 2015; 12:697.

10. Lee JY, Ryu DS, Kim WJ, Kim SJ. Aberrantly expressed microRNAs in the context of bladder tumorigenesis. Investigative and Clinical Urology. 2016; 57:S52-S59.

11. Emig D, Salomonis N, Baumbach J, Lengauer T, Conklin BR, Albrecht M. AltAnalyze and DomainGraph: analyzing and visualizing exon expression data. Nucleic Acids Res. 2010; 38:W755-W762.

12. Hecker N, Stephan C, Mollenkopf HJ, Jung K, Preissner R, Meyer HA. A new algorithm for integrated analysis of miRNA-mRNA interactions based on individual classification reveals insights into bladder cancer. PLoS One. 2013; 8:e64543.

13. Dhawan D, Paoloni M, Shukradas S, Choudhury DR, Craig BA, Ramos-Vara JA, Hahn N, Bonney PL, Khanna C, Knapp DW. Comparative gene expression analyses identify luminal and basal subtypes of canine invasive urothelial carcinoma that mimic patterns in human invasive bladder cancer. PLoS One. 2015; 10:e0136688.

14. Jiang C, Xuan Z, Zhao F, Zhang MQ. TRED: a transcriptional regulatory element database, new entries and other development. Nucleic Acids Res. 2007; 35:D137-140.

15. Andersson KE, Wein AJ. Pharmacology of the lower urinary tract: basis for current and future treatments of urinary incontinence. Pharmacol Rev. 2004; 56:581-631.

16. Choschzick M, Lebeau A, Marx AH, Tharun L, Terracciano L, Heilenkotter U, Jaenicke F, Bokemeyer C, Simon R, Sauter G, Schwarz J. Overexpression of cell division cycle 7 homolog is associated with gene amplification frequency in breast cancer. Hum Pathol. 2010; 41:358-365.

17. Xu S, Witmer PD, Lumayag S, Kovacs B, Valle D. MicroRNA (miRNA) transcriptome of mouse retina and identification of a sensory organ-specific miRNA cluster. J Biol Chem. 2007; 282:25053-25066.

18. Li P, Sheng C, Huang L, Zhang H, Huang L, Cheng Z, Zhu Q. MiR-183/-96/-182 cluster is up-regulated in most breast cancers and increases cell proliferation and migration. Breast Cancer Res. 2014; 16:473.

19. Han Y, Chen J, Zhao X, Liang C, Wang Y, Sun L, Jiang Z, Zhang Z, Yang R, Li Z, Tang A, Li X, Ye J, et al. MicroRNA expression signatures of bladder cancer revealed by deep sequencing. PLoS One. 2011; 6:e18286.

20. Liu Y, Han Y, Zhang H, Nie L, Jiang Z, Fa P, Gui Y, Cai Z. Synthetic miRNA-mowers targeting miR-183-96-182 cluster or miR-210 inhibit growth and migration and induce apoptosis in bladder cancer cells. PLoS One. 2012; 7:e52280.

21. Wszolek MF, Rieger-Christ KM, Kenney PA, Gould JJ, Silva Neto B, Lavoie AK, Logvinenko T, Libertino JA, Summerhayes IC. A MicroRNA expression profile defining the invasive bladder tumor phenotype. Urol Oncol. 2011; 29:794-801 e791.

22. Korpal M, Kang Y. The emerging role of miR-200 family of microRNAs in epithelial-mesenchymal transition and cancer metastasis. RNA Biol. 2008; 5:115-119.

23. Dykxhoorn DM, Wu Y, Xie H, Yu F, Lal A, Petrocca F, Martinvalet D, Song E, Lim B, Lieberman J. miR-200 enhances mouse breast cancer cell colonization to form distant metastases. PLoS One. 2009; 4:e7181.

24. Lee JS, Leem SH, Lee SY, Kim SC, Park ES, Kim SB, Kim SK, Kim YJ, Kim WJ, Chu IS. Expression signature of E2F1 and its associated genes predict superficial to invasive progression of bladder tumors. J Clin Oncol. 2010; 28:2660-2667.

25. Wang H, Ach RA, Curry B. Direct and sensitive miRNA profiling from low-input total RNA. RNA. 2007; 13:151-159.

26. Takeshita N, Hoshino I, Mori M, Akutsu Y, Hanari N, Yoneyama Y, Ikeda N, Isozaki Y, Maruyama T, Akanuma N, Komatsu A, Jitsukawa M, Matsubara H. Serum microRNA expression profile: miR-1246 as a novel diagnostic and prognostic biomarker for oesophageal squamous cell carcinoma. Br J Cancer. 2013; 108:644-652.

27. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP. Exploration, normalization, and summaries of high density oligonucleotide array probe level data. Biostatistics. 2003; 4:249-264.

28. Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012; 7:562-578.

29. Sun J, Nishiyama T, Shimizu K, Kadota K. TCC: an R package for comparing tag count data with robust normalization strategies. BMC Bioinformatics. 2013; 14:1-14.

30. Friedlander MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012; 40:37-52.

31. Vlachos IS, Zagganas K, Paraskevopoulou MD, Georgakilas G, Karagkouni D, Vergoulis T, Dalamagas T, Hatzigeorgiou AG. DIANA-miRPath v3.0: deciphering microRNA function with experimental support. Nucleic Acids Res. 2015; 43:W460-466.

32. Zuberi K, Franz M, Rodriguez H, Montojo J, Lopes CT, Bader GD, Morris Q. GeneMANIA prediction server 2013 update. Nucleic Acids Res. 2013; 41:W115-122.


Creative Commons License All site content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 License.
PII: 24441