Identification of potential therapeutic targets in urothelial bladder carcinoma of Chinese population by targeted next-generation sequencing

ABSTRACT

Patients with urothelial carcinoma (UC) of the bladder have a high risk of death in China. However, a lack of comprehensive molecular profiling in Chinese Han population hinders the development of targeted therapies for bladder cancer. In our present study, we collected fresh bladder tumors from low-grade (T1, N0, M0, G1) non-muscle invasive bladder cancer (NMIBC) patients (n = 16) and high-grade (T2-4, N0, M0, Gx) muscle-invasive bladder cancer (MIBC) patients (n = 16) with their paired normal bladder tissues, and subjected the total genomic DNAs to targeted next-generation sequencing (NGS) for 94 cancer-associated genes. NGS results showed that 30.9% of detected genes (29/94) was mutated in 32 urothelial carcinoma bladder tissues. Furthermore, our results and ICGC database showed that FGFR3, KMT2D, TP53, KDM6A, and ARID1A were the most frequently mutated genes in UC patients. Of note, NMIBC and MIBC displayed distinguishable genomic alterations. FGFR3, KMT2D, AKT1, ARID1A, and STAG2 were the most frequently mutated genes in NMIBC patients, whereas mutations of TP53, CREBBP, FGFR3, KDM6A, KMT2D, and ARID1A were frequently detected in MIBC. Intriguingly, gene ontology and clustering analysis revealed that these frequently mutated genes were highly enriched in signaling pathways responsible for cancer development. Taken together, the mutation frequency of genes associated with UC development in NMIBC and MIBC was screened out in Chinese Han population and elucidation of the related mechanisms provides theoretical basis and technical support for the development of early diagnosis and therapeutic strategies in UC.

KEYWORDS:

• Bladder cancer
• gene mutations
• NGS
• cancer diagnosis and therapeutics

Introduction

Bladder cancer (BCa) is one of the common malignancies worldwide with global prevalence at approximately 1 million each year,¹ which ranked the ninth in the world with sixth in men and tenth in the middle age of women.^2,3 In China, the overall incidence of bladder cancer was about 8.05/10 million in 2015 and continued to increase in recent years.⁴ Epidemiological studies have shown that genetic factors,^5,6 environmental factors such as active smoking and passive smoking,^7,8 occupational exposure in chemical and petrochemical industry,⁹ air pollution,¹⁰ food and drink,¹¹ use of medicines,¹² and sexual distinction¹³ can all affect the development and progression of BCa. In particular, tobacco smoking and gender difference are two important risk factors for BCa incidence.^7,8 The pathogenesis of BCa involves a series of aberrant biological processes such as rapid DNA replication, increased cell proliferation, suppressed cell apoptosis, tumor cell invasion and metastasis, immune escape, and angiogenesis.^14,15 Molecular biological analysis also revealed that alterations in oncogenes or tumor suppressor genes led to the development of BCa.¹⁶

Typically, most BCa (70-80%) are superficial tumors, also called non muscle invasive BCa (NMIBC), with Bacillus of Calmette-Guerin (BCG) or chemotherapy treatment.^17,18 Nevertheless, half of them can relapse and 5% to 20% of NMIBC will eventually become muscle-invasive BCa (MIBC in short) despite of treatments.^19,20 On the contrary, 50% of MIBC patients with systemic chemotherapy or radical cystectomy eventually die due to metastases within 2 years of diagnosis.²¹ Currently, no effective strategies have been applied to cure MIBC, suggesting an emergent development of novel targeted therapies to cure MIBC is needed.

Next generation sequencing (NGS) allows researchers to study biological functions at a massively broad level. NGS has been applied to decipher genomic alterations in BCa by numerous groups.^22–26 Indeed, the active forms of FGFR3 and HRAS due to genomic mutation have been identified in the early stage of BCa, serving as biomarkers for NMIBC.^27–30 Additionally, TP53 mutation and loss of Rb were associated with MIBC incidence.^6,31-33 These facts indicated that a better understanding of the genomic landscape of BCa was essential to improve its therapeutic strategies. Despite of its significance, the molecular profile of BCa at the genomic level in Chinese Han population, however, hasn’t been investigated yet. Here, we performed targeted NGS to test the mutation frequency of 94 well-known cancer-associated genes in 32 bladder cancer tissues. Our aim was to explore the potential biomarkers for BCa diagnosis and prognosis in order to provide novel targets for the treatment of BCa.

Our results showed that 30.9% of detected genes were mutated in our study cohort. MIBC and NMIBC displayed a unique mutation signature. Those unique-mutated genes may serve as potential biomarkers or therapeutic targets for BCa treatment.

Materials and methods

Clinical samples collection

Fresh clinical BCa tissues were collected according to the ethical guidelines of the Declaration of Helsinki (1975) and the Institutional Medical Ethics Committee of Xiamen University. Pathological diagnosis of BCa was verified by at least two qualified pathologists. A total of 32 fresh BCa tissues and their corresponding pairs of normal tissues were collected from low stage and low grade (T1, N0, M0, G1) of non-muscle invasive bladder cancer (NMIBC) patients and high stage and high grade (T2-4, N0, M0, G2-G3) of muscle-invasive samples (MIBC) patients. Sixteen pairs of human NMIBC and MIBC specimens were collected from the Department of Urology, The First Affiliated Hospital of Xiamen University from January 2016 to March 2017.

Genomic DNA extraction and target sequencing

Genomic DNA was extracted from RNA Later-preserved bladder cancer samples using DNeasy Blood & Tissue Kit (QIAGEN, Germany, Cat No.69506). The detection of genomic DNA quality was performed on 1% agarose gel, after that the DNA concentration was measured using a Qubit® DNA Assay Kit in Qubit® 2.0 Fluorometer (Invitrogen, USA) according to the manufacturer’s instruction.

To get the target gene regions, we designed probes on the website of Agilent (http://www.agilent.com/) for the 94 genes according to the design description. Initially, DNA fragmentation was performed by hydrodynamic shearing machine (Covaris, Massachusetts, USA) to generate 180–280 bp fragments. Extracted DNA was then subjected to ligation-mediated PCR (LM-PCR), followed by the product purification and hybridization with the probe for enrichment. All LM-PCR products (non-captured or captured) were used to estimate the magnitude of enrichment by real-time PCR. The captured library was subjected to cluster generation in an Illumina flow cell. Then, sequencing was performed to produce 150 bp paired-end reads by using an Illumina Hiseq platform. Finally, to ensure each sample meets the desired average fold coverage, the captured libraries were sequenced independently.

Sequence data quality control

The shortreads (Raw data) were transformed from the original fluorescence image files obtained from HiSeq platform by base calling. They were recorded in FASTQ format, including sequence data and their corresponding quality scores. Reads containing adapter contamination and low-quality/unrecognizable nucleotides were excluded through quality control. Then, downstream bioinformatics analyses were performed on clean data. At the same time, percentage of reads with average quality >20 and with average quality >30, sequencing error rate, the total reads number, and GC content distribution were calculated.

Reads mapping and somatic genetic alteration detection

Original mapping results were obtained and stored in BAM format by mapping valid sequencing reads to the human reference genome (UCSC hg19) using Burrows-Wheeler Aligner (BWA) software. After that, BAM files were sorted with SAM tools and Picard (http://broadinstitute.github.io/picard/). The resulting BAM files were obtained and used for calculating the sequence coverage and depth after duplicate marking, local realignment, and base quality recalibration.

MuTect and Strelka were used to call somatic single nucleotide variations (SNVs), small insertions, and deletions (InDels) from paired tumor-normal samples. In addition, polymorphisms of somatic SNVs and InDels referenced in the 1000 Genomes Project or Exome Aggregation Consortium (ExAC) with a minor allele frequency over 1% were removed to default filters. Finally, ANNOVAR was used for functional annotation of genetic variants.³⁴

Control-FREEC was used to identify somatic copy number variations (CNVs) and GISTIC algorithm was used to infer recurrently amplified or deleted genomic regions. Based on the frequency and amplitude of amplification or deletion affecting each gene, G-scores were calculated for genomic and gene-coding regions. Amplification or deletion region with a G-score > 0.1 and a corresponding p-value less than 0.05, based on a permutation-derived null distribution, was defined as a significant CNV region.

Results

Patient characteristics

Among these 32 identified BCa patients, 5 patients were female and 27 were male, with age of 33 y to 82 y (Table 1). Pathologically, 16 patients were classified as NMIBC and others were MIBC. The maximum diameter of mass from 27 patients (84.4%) was less than 5 cm while only 5 patients have their mass larger than 5 cm in diameter. Fifty percent of patients (n = 16) were in TNM stage I and 8 patients (25%) were in stage II, while only 3 and 5 patients had stage III and IV disease (9.375% and 15.625%), respectively. Sixteen (50%) patients were G1 degree and 16 (50%) had G2-G3 degree.

Table 1. Clinical information of BCa patients in our study.

Clinicopathological features		n	%
Gender	M	27	84.375
	F	5	15.625
Age	<60	14	43.75
	≥60	18	56.25
Maximum diameter of mass	<5 cm	27	84.4
	≥5 cm	5	15.6
TNM stage	I	16	50
	II	8	25
	III	3	9.375
	IV	5	15.625
Lymph node status	≥1	2	6.25
	0	30	93.75
Degree of differentiation	G1	16	50
	G2-G3	14	50
classification	NMIBC	16	50
	MIBC	16	50

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

Genomic characterization of NMIBC and MIBC cases

To characterize the landscape of genomic mutations in Chinese Han population, we collected 32 BCa tissues (the detailed clinical information is shown in Table 1) and their paired normal bladder controls, then targeted next generation sequencing was applied to these samples for 94 well defined UC-associated genes (Table S1). The average sequencing depth for the Normal and tumor samples from 32 UBC patients analyzed by NGS was 474.5×, with 99.98% of the reference human genome covered as shown in Supplementary Figure 1. In the combined discovery and validation cohorts, we identified an average of 5595.4 point mutations and 806.3 inDels per sample. According to the NGS results, we noticed that 53.1% of UC patients (17/32) had different degree of gene mutation (Figure 1a). In addition, 30.9% of detected genes (29/94, Table 2 & Table S2) were mutated in 32 urothelial carcinoma bladder tissues. Of note, FGFR3 (31.3%), KMT2D (18.8%), TP53 (12.5%), ARID1A (12.5%), AKT1 (12.5%), KDM6A (9.4%), STAG2 (9.4%), and LRP1B (9.4%) were the most frequently mutated genes in UC patients (Figure 1a & Table 2). We also analysis the sequencing data of 103 China bladder cancer samples download from ICGC database and found that TP53 (27.2%), KDM6A (24.3%), FGFR3 (16.5%), HRAS (15.5%), ERBB2 (13.6%), CREBBP (12.6%), ERCC2 (12.6%), ARID1A (11.7%) and XIRP2 (10.7%) were the top-mutated genes (Figure 1b & Table S3). Consistent with our results, we all found coding genes of histone modification-related genes (KMT2D, KDM6A, EP300, and CREBBP), FGFR3, TP53, ARID1A, AKT1, STAG2, LRP1B, HRAS2, ERBB2, ERCC2, and FAT1 were frequently mutated genes in UC patients (Figure 1). However, different from their results, FGFR3 is the most frequently mutated gene in our results other than TP53 in their results (Figure 1). This may be due to sample capacity, different regions, diets and living habits, however, it needs more samples to further verify. Previous studies showed that the FGFR3 mutations occured predominantly in Urothelial cell carcinoma, whereas TP53 mutations were found in over 50% of human cancers.^35,36 Data from The Cancer Genome Atlas (TCGA) Research Network showed that TP53 (49%), MLL2 (49%), ARID1A (49%), KDM6A (49%), PIK3CA (49%), EP300 (49%) and CDKN1A (49%) were the most frequently mutated genes.³⁷ These results are consistent with our results and complement each other, elucidate the molecular mutation basis for BCa development.

Table 2. Summary of gene mutation frequency in 32 BCa tissues.

Gene	Gene ID	Description	Mutation frequency (Total/NMIBC/MIBC, %)
FGFR3	NM_000142.4	fibroblast growth factor receptor 3	31.3/50/12.5
KMT2D	NM_003482.3	lysine methyltransferase 2D	18.8/31.3/6.3
TP53	NM_000546.5	tumor protein p53	12.5/0/25
ARID1A	NM_006015.5	AT-rich interaction domain 1A	12.5/18.8/6.3
AKT1	NM_005163.2	AKT serine/threonine kinase 1	12.5/25/0
KDM6A	NM_021140.3	lysine demethylase 6A	9.4/6.3/12.5
STAG2	NM_001042749.2	stromal antigen 2	9.4/18.8/0
LRP1B	NM_018557.2	LDL receptor related protein 1B	9.4/12.5/6.3
HRAS	NM_005343.3	HRas proto-oncogene, GTPase	6.3/12.5/0
CREBBP	NM_004380.2	CREB binding protein	6.3/0/12.5
ERCC2	NM_000400.3	ERCC excision repair 2, TFIIH core complex helicase subunit	6.3/0/12.5
EP300	NM_001429.3	E1A binding protein p300	6.3/0/12.5
ATM	NM_000051.3	ATM serine/threonine kinase	6.3/12.5/0
FAT1	NM_005245.3	FAT atypical cadherin 1	6.3/6.3/6.3
ERBB2	NM_004448.3	erb-b2 receptor tyrosine kinase 2	3.1/6.3/0
XIRP2	NM_001024618.2	xin actin-binding repeat containing 2	3.1/6.3/0
ERBB3	NM_001982.4	erb-b2 receptor tyrosine kinase 3	3.1/12.5/6.25
DCHS2	NM_001358235.1	dachsous cadherin-related 2	3.1/12.5/12.5
CACNA1D	NM_000720.3	calcium voltage-gated channel subunit alpha1 D	3.1/0/6.3
RANBP2	NM_006267.4	RAN binding protein 2	3.1/6.3/0
FBXW7	NM_033632.3	F-box and WD repeat domain containing 7	3.1/0/6.3
BRCA2	NM_000059.3	BRCA2, DNA repair associated	3.1/6.3/0
EGFR	NM_005228.5	epidermal growth factor receptor	3.1/6.3/0
RARG	NM_000966.6	retinoic acid receptor gamma	3.1/6.3/0
OBSCN	NM_052843.3	obscurin, cytoskeletal calmodulin and titin-interacting RhoGEF	3.1/0/6.3
RBM10	NM_005676.4	RNA binding motif protein 10	3.1/6.3/0
ZFP36L1	NM_004926.3	ZFP36 ring finger protein like 1	3.1/6.3/0
CCND1	NM_053056.2	cyclin D1	3.1/0/6.3
RAF1	NM_002880.3	Raf-1 proto-oncogene, serine/threonine kinase	3.1/6.3/0

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

Figure 1. The genomic landscape of bladder cancer. (a) Targeted NGS analysis of the mutated genes in 32 bladder cancer samples. (b) ICGC data showing the mutated genes in 103 bladder cancer samples of Chinese Han population.

Furthermore, subgroup analysis showed that FGFR3 (50.0%, 8/16), KMT2D (31.3%, 5/16), AKT1 (25.0%, 4/16), ARID1A (12.5%, 2/16) and STAG2 (18.8%, 3/16) were the most frequently mutated genes in NMIBC patients (Figure 2a & Table 2). FGFR3, one member of fibroblast growth factor receptor family, mediates cell proliferation in the response to FGF stimulation. The activation of FGFR3 due to gene amplification or gene mutation triggers numerous signaling pathways which are involved in the carcinogenesis of UC.^38,39 Consistent with a previous report documenting that FGFR3 mutation (49%) was frequently detected in NMIBC,²⁶ our sequencing result also confirmed this finding. Another prevalently mutated gene was KMT2D (50%). As a histone methyltransferase, KMT2D mediates epigenetic modification by methylating lysine-4 position of histone 3,⁴⁰ the mutation of which may lead to aberrant epigenetic profile occurred in NMIBC. Therefore, our sequencing result not only strengthens the prevalence of FGFR3 and KMT2D in NMIBC but also directs further scientific investigations.

Figure 2. Twenty-nine mutated genes in NMIBC and MIBC tissues assayed by Targeted NGS. (a) Gene mutation frequency in NMIBC. (b) Gene mutation frequency in MIBC. Abbreviations: NMIBC, non-muscle invasive bladder cancer; MIBC, muscle-invasive bladder cancer.

In addition, NSG results also showed that TP53 (25.0%, 4/16), CREBBP (12.5%, 2/16) and FGFR3 (12.5%, 2/16) were the most frequently mutated genes in MIBC patients (Figure 2b & Table 2). Of note, TP53, CREBBP, CACNA1D, and CCND1 were only frequently mutated in MIBC but not in NMIBC patients (Figure 2 & Table 2). We have followed up the operation patients and found there are 25% (8/32) of patients had relapse or metastasis (Data not shown). Among the genes with high mutation frequency, 75% (3/4) of patients with TP53 mutation, 66.7% (2/3) of patients with LRP1B or ARID1A and 33.3% (2/6) of patients with KMT2D mutation had relapse and metastasis, suggesting their mutations were correlated with recurrence, metastasis, or associated with other tumors such as gastric cancer and ovarian cancer.

Among these mutated genes, CACNA1D (Calcium Voltage-Gated Channel Subunit Alpha1 D), which is located in cell membrane and mediates calcium homeostasis, is highly associated with neuron diseases.^41,42 However, the detailed function of CACNA1D in the development of MIBC is still elusive. We speculate that the CACNA1D mutation disrupts calcium homeostasis and alters calcium-related signaling. In comparison with other reports, CACNA1D seems to have a higher mutation rate in Chinese Han population, suggesting the geographic disparity exists in the carcinogenesis of BCa. Interestingly, another mutated gene is CREBBP (CREB binding protein), which complexes with phosphorylated CREB and enhances the transcription of cAMP-responsive genes. CREBBP also functions to regulate gene expression via acetylating histone and allowing chromatin to a relaxed state.^43,44 Accordingly, an other group also demonstrated that CREBBP was commonly mutated in high grade BCa (13%).²⁴ Taking together, scientists have compelling rationale to develop these mutated genes as potential biomarkers or therapeutic targets for MIBC treatment.

Signaling pathway alterations

To investigate which signaling pathways these mutated genes are correlated with, we performed gene ontology and clustering analysis by using KEGG and Reactome program. Data revealed that these frequently mutated genes were tightly associated with signaling pathways which were implicated into cancer development, such as PI3 K/AKT, MAPK, Jak-STAT, FOXO and ERBB signaling pathways (Figure 3 & Figure 4). Previous studies have reported that dysfunction of these signaling pathways drives cancer progression and targeting them with their specific inhibitors slows down tumor growth in various cancer types.^45,46 Collectively, our NGS result suggests that activation of oncogenic signaling pathways (PI3 K/AKT, MAPK, Jak-STAT, FOXO and ERBB) in BCa development may partially result from genomic alterations.

Figure 3. KEGG gene enrichment analysis showed that these mutated genes are highly correlated with 30 physiological and pathological pathways.

Figure 4. Reactome program revealed that these mutated genes were enriched in signaling pathways which were implicated into cancer development.

Discussion

BCa is one of the leading causes of cancer deaths worldwide and continues to threaten human health. Although current treatments show promising therapeutic benefits, BCa can recur and advance to a high stage with limited curable options. Therefore, understanding of the molecular profiling of BCa not only provides rationale to develop targeted therapies but also helps identify potential diagnostic or prognostic biomarkers. Although genomic alterations in BCa have been comprehensively investigated by analyzing TCGA database or performing next negation sequencing by other groups,^24,37 however, geographic disparity in BCa led us to screen gene mutations in Chinese Han population. Our results showed that 29/94 genes were mutated in 32 bladder cancer tissues, indicating a higher mutation burden in Chinese BCa patients. According to our analysis, a differential profile of mutated genes existed between NMIBC and MIBC, which may serve as the rationale to develop promising biomarkers to distinguish these two types of BCa. However, further investigation is required to strengthen our results given the relatively small sample size of our study.

Previous investigations into genomic alteration in Chinese Han population have identified some mutated genes. For instance, Cai et al. found that three genes (STAG2, ESPL1, TACC3), were involved in sister chromatid cohesion and segregation (SCCS) process and-, displayed a mutation in bladder cancer with transitional cell carcinoma.⁴⁷ Wu S. et al. performed whole-genome sequencing in 65 urothelial bladder carcinomas (UBCs) and targeted sequencing in an additional 196 UBCs. Besides FGFR3, TP53, and PIK3CA, they also found ZFP36L1, HRAS, and KDM6A were frequently mutated in UBCs. Notably, they find recurrent ADGRG6 enhancer mutations and FRS2 duplications which are associated with higher protein expression in the tumor and poor prognosis.⁴⁸ In our study, we also identified a mutated STAG2 and ZFP36L1 in BCa patients. This consistent data suggest that dysfunctional SCCS may be implicated into BCa carcinogenesis. Cai et al. also demonstrated that some chromatin remodeling genes (UTX, MLL-MLL3, EP300, NCOR1, and ARID1A) were frequently mutated in transitional cell carcinoma of BCa.⁴⁹ Our result revealed that three genes (KMT2D, EP300 and ARID1A) involved in chromatin remodeling were genetically changed in UC patients, implying that aberration in chromatin remodeling is a hallmark of BCa patients. Novel genetic alterations in BCa have continuously been identified in Chinese Han population,⁵⁰ however, none of these genes was genetically changed on the basis of our sequencing result, which may be explained by the fact that evident difference on genetic background exists in BCa patients.

Cystoscopy is the primary diagnostic tool for BCa patients.⁵¹ Nevertheless, it is an invasive detection method because it needs urologists to take biopsy from patients. Therefore, a noninvasive diagnostic tool is emergently required to replace cystoscopy. Circulating tumor cells (CTCs) are tumor cells that move to the bloodstream, detection of biomarkers in CTCs represents a promising noninvasive tool in numerous cancers.^52,53 Here, a higher mutation frequency of FGFR3 (50.0%), KMT2D (31.3%), and AKT1 (25.0%) was observed in NMIBC but not in MIBC, indicating these three mutation events occur in the early stage of urothelial carcinogenesis. Therefore, it may be practical to detect the mutated form of FGFR3, KMT2D, and AKT1 in CTCs for the diagnosis of early stage of BCa.

In addition, we found that several mutants (TP53, CREBBP, CACNA1D, CCND1, FAT1, FBXW7, and OBSCN) were only detected in MIBC but not in NMIBC. As to MIBC patients, few therapeutic options were effectively utilized to cure them. Consequently, investigating the biological functions of these mutated genes may provide us identify novel targets for MIBC. Actually, studies have also shown that mutation of TP53 has been identified to be associated with MIBC.³¹ Here, in addition to TP53, we identified more mutated genes specific to MIBC. Our NGS results and ICGC data analysis also showed that FGFR3 and TP53 are the two most frequently mutated genes (Figure 1). Interestingly, patients with FGFR3 mutation generally do not share with TP53 mutation (Figure 1). So we can divide bladder cancer patients into FGFR3-mutation or TP53-mutation types. These mutated genes may give us a variety of directions to elucidate the underlying mechanisms responsible for MIBC development.

For the first time, we showed that CACNA1D has a high mutation rate in MIBC, suggesting the calcium pathway may play important roles in BCa progression of in Chinese Han population. Actually, calcium ATPase inhibitors have been pre-clinically used for cancers such as prostate cancer and breast cancer.^54,55 However, these inhibitors have been limitedly tested in BCa. These results provide insights into potential effects of ATPase inhibitors in MIBC treatment.

Disclosure of potential conflicts of interest

No potential conflicts of interest were disclosed.

Competing interests

The authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

This work was granted by Fujian Provincial Department of Science and Technology (^#2016D009), Xiamen science and technology plan (^#3502Z20184014), Natural Science Foundation of Fujian Province (^#2017D0010 and ^#2017J01355), Young and Middle-aged Backbone Foundation of the Fujian Provincial Health and Family Planning Commission (^#2019-ZQNB-26) and The Science Fund founded by the First Affiliated Hospital of Xiamen University for Young Scholars (^#XYY2017004).