Identification of circular RNA signature in bladder cancer.

Circular RNA (circRNA) comprises a class of endogenous species of RNA consisting of a circular loop that is crucial for genetic and epigenetic regulation. The significance of circRNA in bladder cancer (BCa) remains to be investigated. Here we performed genome‑wide circRNA analysis of 5 paired tumour and adjacent normal tissue samples from BCa patients via next generation sequencing (NGS) technology. Next we confirmed NGS data in a separate set of 32 paired BCa samples using quantitative real-time reverse transcription polymerase chain reaction. The results showed that circRNA profile presented a total of 88,732 circRNA in BCa samples. Among them, 14 were upregulated and 42 were downregulated with q-values of <0.001 and fold changes of ≥2 or ≤0.5. The expression level changes of hsa_circ_0091017 and hsa_circ_0002024 in the 32 paired samples were in accord with NGS data. In conclusion, we identified a set of circRNAs that are potentially implicated in the tumorigenesis of BCa and could serve as novel diagnostic markers for BCa.

Introduction

Urothelial bladder cancer (BCa) is the seventh most common cancer among the male population worldwide. Approximately 75% of BCa patients with high recurrence rates present a tumour confined to the mucosa or submucosa, whereas the rest have muscle-invasive BCa with evident aggressive behaviour, which is the main threat to the patients. At present, there is still no effective treatment beyond surgery for bladder cancer. Therefore, understanding the molecular mechanism of BCa development and progression is important to develop novel therapeutic approaches.

Circular RNA (circRNA) was first found in RNA viruses 1. With the advances in next generation sequencing (NGS) technology and bioinformatics, circRNA was discovered as a novel class of abundant endogenous non-coding RNA with regulatory potency in mammalian cells 2 . Unlike normal linear RNA, the 3′ and 5′ ends of circRNA are covalent bonds that are linked to form closed loop structures without polarity nor polyadenylated tails, thereby making them stable 3. CircRNA can function as a microRNA sponge to contribute to competing endogenous RNA (ceRNA) networks. CircRNA can also bind and sequester RNA binding proteins (RBPs) to form RNA protein complexes (RPCs) that can regulate the pool of RBPs and influence gene expression 4,5.

Recently, increasing studies have shown aberrant expression of circRNA in various malignancies, such as laryngeal cancer, basal cell carcinoma, cutaneous squamous cell carcinoma, colorectal cancer, ovarian carcinoma and pancreatic ductal adenocarcinoma 6-10. Moreover, evidence has shown that circRNA plays a critical role in cancer-related pathways. Cir-ITCH can function as a tumour suppressor in oesophageal squamous cell carcinoma and colorectal cancer 11,12. In hepatocellular carcinoma, the expression of circRNA Cdr1as is upregulated while hsa_circ_0001649 is downregulated 13,14. Furthermore, circRNA is proposed as promising biomarker for cancer diagnosis 15.

Huang et al. examined lncRNA/circRNA and mRNA expression profiles in BCa patients using microarray analysis 16. Compared with microarray analysis, NGS is an open technology that can identify undiscovered nucleotide sequences. In this study, we comprehensively investigated aberrant circRNA expression in BCa patients using NGS technology in order to understand the potential role of circRNA in BCa.

Material and Methods

Patients and specimens

The study included five patients with BCa who underwent partial or radical cystectomies at the First Affiliated Hospital of Nanjing Medical University. Patients with a history of cancer and any preoperative chemotherapy or radiotherapy were excluded. After surgery, the matched specimens were immediately preserved in liquid nitrogen until use. All patient samples were confirmed by pathological examination. The pathological diagnosis was in accordance to tumor node-metastasis classification for the stage (2002 International Union Against Cancer) and histopathologic grade (WHO 1973). This study was approved by the Institutional Review Board of the First Affiliated Hospital of Nanjing Medical University, China. Written informed consent was obtained from all the participants in this study.

Total RNA isolation and quality control (QC)

Total RNA was isolated from samples using TRIzol reagent (Life Technologies, CA, USA) following the protocol. The quantity and quality of total RNA samples were measured using NanoDrop ND-1000 (Thermo Scientific, DE, USA). RNA integrity was assessed via electrophoresis using denaturing agarose gel. The spectrophotometric absorbance of RNA samples at 230, 260 and 280 nm and ratios of D260:D280 and D260:D230 were used estimate the purity of total RNA. RNA samples were stored at -80°C before use.

Next generation RNA sequencing analysis

Total RNA from five matched samples of bladder cancer and adjacent non-tumour tissue were treated with Epicenter Ribo-Zero rRNA Removal Kit (Illumina, CA, USA) and RNAse R (Epicenter, CA, USA) to remove ribosomal and linear RNA. Then, the RNA-seq libraries were constructed using TruSeq® Stranded Total RNA HT/LT Sample Prep Kit (Illumina, CA, USA). Sequencing was determined on Illumina Hiseq 2500 instrument. All the sequencing procedures and analyses were performed in RiboBio (Guangzhou, China).

Polymerase chain reaction and Sanger sequencing

The cDNA templates were PCR-amplified using Veriti® 96-Well Thermal Cycler (Applied Biosystems, CA, USA) with 2×Taq PCR Mastermix (Tiangen Biotech, Beijing, China). A total of 40 PCR cycles were performed, and PCR products were visualised via electrophoresis in 1.5% ethidium bromide-stained agarose gel with D2000 DNA Marker (Tiangen Biotech, Beijing, China). PCR products were purified using TIANgel Midi Purification Kit (Tiangen Biotech, Beijing, China) followed by Sanger sequencing using 3730xl DNA Analyser (Applied Biosystems, CA, USA).

Quantitative real-time PCR (qPCR)

First strand cDNA synthesis was performed using HiScript® II Q RT SuperMix for qPCR with gDNA wiper (Vazyme Biotech, Nanjing, China). Real-time PCR was conducted in triplicate on ROCHE LightCycler® 480 (ROCHE, Basel, Switzerland) with SYBR® Green Master Mix (Vazyme Biotech, Nanjing, China). The primers for qPCR were synthesised in Sango Biotech (Shanghai, China). The convergent primers were used to detect circRNA expression, and the results were normalised against β-actin. The following primers were used: hsa_circ_0091017 F: 5′-ATCATTCGGACAAGACAGG-3′, R: 5′-TCCATATACCAGATCCTCATTG-3′. hsa_circ_0002024 F: 5′-ATGCCCTGCTGACAAG-3′, R 5′-CACGAGTGCCATGATTT-3′, β-actin F: 5′-AGCGAGCATCCCCCAAAGTT-3′ and R 5′-GGGCACGAAGGCTCATCATT-3′. The relative expression was normalised with β-actin using the formula circRNA/β-actin = 2-δCt.

Statistical analysis

Statistical analyses were performed using GraphPad Prism 5 (GraphPad Software, CA, USA), R software version 3.2.1 (http://www.r-project.org/) and Microsoft Excel (Microsoft, DC, USA). Significance was determined via Student's t-test for differences in circRNA expression. A two-sided P value < 0.05 was considered statistically significant.

Results

Identification of differentially expressed circRNA in BCa

First we characterised circRNA transcripts via RNA sequencing analysis of ribosomal RNA and linear RNA-depleted total RNA from five matched BCa tissue samples. Each sample was sequenced on Illumina HiSeq and yielded an average of 70 million reads, which were mapped onto the human reference genome (GRCh37/hg19) by TopHat2 17,18. A total of 88,732 distinct circRNA candidates were found in all the samples. Among these, 62,788 contained at least two unique back-spliced reads.

We then defined the statistical criteria for selecting aberrant-expressed circRNA using a q-value of <0.001 with a fold change of >2.0 or <0.5. A total of 14 circRNA were upregulated and 42 circRNA were downregulated significantly (Table 1). The differential expression for circRNA between cancerous and non-cancerous samples was shown in the heat-map and dendrogram generated via unsupervised clustering that was blinded to the clinical annotations (Fig. 1). The Scatter plot was used to assess the variation of circRNA expression between BCa and adjacent normal samples (Fig. 2).

Table 1

Differentially expressed circRNA between bladder cancer samples and adjacent normal samples.

Gene	Gene symbol	Fold Change	p-value	q-value
hsa_circ:chr12:6619445-6619698	NCAPD2	2.443724	6.96E-24	1.28E-19
hsa_circ:chrX:153590053-153593816	FLNA	-2.68072	6.75E-21	9.27E-17
hsa_circ:chr15:99670080-99674092	SYNM	-2.79424	1.36E-14	1.50E-10
hsa_circ:chr6:58776544-58779091	*	1.355364	9.76E-14	8.94E-10
hsa_circ:chr16:15797304-15813543	MYH11	-3.73008	4.89E-13	3.84E-09
hsa_circ:chr11:2151469-2151932	IGF2	4.947237	1.06E-12	7.27E-09
hsa_circ:chr2:220285217-220285660	DES	-3.86951	2.50E-12	1.53E-08
hsa_circ:chr19:861699-861956	CFD	-3.20675	2.72E-12	1.49E-08
hsa_circ:chr11:2150673-2151176	IGF2	4.157646	4.58E-11	2.29E-07
hsa_circ:chr14:20811315-20811566	RPPH1	1.117892	2.28E-10	1.57E-06
hsa_circ:chrX:153588710-153589861	FLNA	-2.58789	2.51E-10	1.06E-06
hsa_circ:chrX:12994905-12995266	TMSB4X	1.911979	3.37E-10	1.32E-06
hsa_circ:chrX:153599241-153599729	FLNA	-2.89218	3.77E-10	1.38E-06
hsa_circ:chr17:48264432-48277270	COL1A1	1.279219	6.19E-10	2.13E-06
hsa_circ:chr11:2018801-2019011	H19	3.397296	8.87E-10	2.87E-06
hsa_circ:chr21:47531364-47532738	COL6A2	-2.35656	1.29E-09	3.93E-06
hsa_circ:chr1:59248113-59248470	JUN	-1.85055	1.29E-08	3.74E-05
hsa_circ:chr2:220290701-220291021	DES	-3.63138	1.36E-08	3.73E-05
hsa_circ:chr19:45971553-45971806	FOSB	-2.38711	1.41E-08	3.69E-05
hsa_circ:chr2:220283167-220283423	DES	-2.1107	1.96E-08	4.89E-05
hsa_circ:chr2:74143893-74146554	ACTG2	-3.31293	1.97E-08	4.72E-05
hsa_circ:chrX:153595100-153596345	FLNA	-2.85275	2.09E-08	4.78E-05
hsa_circ:chr4:119951013-119952102	SYNPO2	-3.34258	3.61E-08	7.94E-05
hsa_circ:chrX:153577183-153578111	FLNA	-2.24401	5.68E-08	0.000116
hsa_circ:chrX:153582998-153583402	FLNA	-3.03816	6.12E-08	0.00012
hsa_circ:chr2:220290916-220291368	DES	-2.96065	1.04E-07	0.000197
hsa_circ:chr6:27833099-27833573	HIST1H2AL	3.389217	1.11E-07	0.000203
hsa_circ:chr11:117073718-117075210	TAGLN	-1.97228	1.25E-07	0.000222
hsa_circ:chr1:59248720-59248891	JUN	-2.03477	2.21E-07	0.000379
hsa_circ:chr2:220286137-220290698	DES	-3.93159	2.25E-07	0.000375
hsa_circ:chr19:39897600-39898320	ZFP36	-1.33434	3.45E-07	0.000557
hsa_circ:chrX:153587675-153588499	FLNA	-2.63766	3.49E-07	0.000548
hsa_circ:chr15:63335003-63353138	TPM1	-2.67037	3.80E-07	0.00058
hsa_circ:chrX:153583193-153585991	FLNA	-3.01628	4.21E-07	0.000626
hsa_circ:chr9:16435553-16437522	BNC2	-4.74179	4.25E-07	0.000615
hsa_circ:chrX:153581764-153582376	FLNA	-2.83011	5.54E-07	0.000781
hsa_circ:chr11:33307959-33309057	HIPK3	-4.57536	8.66E-07	0.00119
hsa_circ:chr2:218666117-218667264	TNS1	-2.92975	9.24E-07	0.001239
hsa_circ:chr11:2151281-2151691	IGF2	4.423465	1.00E-06	0.00131
hsa_circ:chr14:75745740-75747633	FOS	-1.73851	1.04E-06	0.001329
hsa_circ:chr9:35682023-35689885	TPM2	-2.41296	1.55E-06	0.001935
hsa_circ:chrX:12994390-12994949	TMSB4X	2.003411	2.02E-06	0.002411
hsa_circ:chr11:2018630-2018977	H19	3.206163	3.31E-06	0.003867
hsa_circ:chrX:153582010-153582613	FLNA	-3.24423	4.89E-06	0.005596
hsa_circ:chr1:161487765-161565384	FCGR2A	2.791146	5.36E-06	0.006014
hsa_circ:chr21:47546090-47551056	COL6A2	-2.46229	6.15E-06	0.006757
hsa_circ:chr14:75747585-75747857	FOS	-2.39272	8.22E-06	0.008863
hsa_circ:chr19:11657490-11658722	CNN1	-3.61675	1.79E-05	0.01893
hsa_circ:chr1:201459304-201465432	CSRP1	-3.06034	1.98E-05	0.020556
hsa_circ:chr2:217537297-217541342	IGFBP5	-1.66851	2.00E-05	0.020314
hsa_circ:chr21:47402550-47406962	COL6A1	-2.03774	2.55E-05	0.025529
hsa_circ:chrX:153580921-153581292	FLNA	-2.59785	4.22E-05	0.041377
hsa_circ:chr14:75745827-75746235	FOS	-1.39446	4.56E-05	0.043928
hsa_circ:chr2:220283415-220283692	DES	-3.46265	4.88E-05	0.046228
hsa_circ:chr19:14628951-14629188	DNAJB1	2.748679	5.12E-05	0.047735
hsa_circ:chrX:153594676-153595218	FLNA	-2.35286	5.32E-05	0.048747

Figure 1

Heat map and unsupervised hierarchical cluster.

Figure 2

Scatter-plot of circular RNA averaged normalized value variation between the adjacent normal tissues (x-axis), and BCa tissues (y-axis). The circRNAs above the top red dashed line and below the bottom red dashed line indicate more than 2.0 fold change in circRNAs between the two groups of samples.

The bar diagram showed all detected reads of circRNA back-spliced junction on genome elements (Fig. 3A). Most differentially expressed circRNA originated from the exons (Fig. 3B). The box plot showed that the distribution of circRNAs intensities from all the datasets was almost similar in all samples (Fig. 4). Moreover, a volcano plot was drawn to visualise differentially expressed circRNA with statistical significance (Fig. 5).

Figure 3

A: Bar diagram of back-spliced junction reads of circRNA on genome elements; B: The category of most differentially expressed circRNAs. (exon: exon derived circRNAs, intron: intron derived circRNAs, intergenic: intergenic region derived circRNAs, outexon: mixed derived circRNAs incluing exon)

Figure 4

Box plot. X-axis: 10 tissue samples; Y-axis: back-spliced junction counts in circRNAs.

Figure 5

Volcano plots to visualise differential expression between bladder cancer samples and non-cancerous samples, which show the relationship between fold change and statistical significance. The red dots in the plot denote differentially expressed circRNA with a fold change of >2 or <0.5 and a q value of <0.001.

Delineation of gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

For the implication of circRNA in gene regulation 24,25, we performed KEGG pathway enrichment analysis for linear transcripts that matched significantly differentially expressed circRNA. Our results indicated that systemic lupus erythematosus and amphetamine addiction were the top related pathways in upregulated and downregulated linear transcripts, respectively (Fig. 6).

Figure 6

KEGG pathway enrichment analyses of the genes that produce upregulated and downregulated circRNAs. Both upregulated and downregulated circRNAs d circRNAs are significantly changed with fold change >2.0 or < 0.5 and q value < 0.001.

Next we performed GO analyses of the linear counterparts that produced aberrantly expressed circRNA. Our data showed that the linear counterparts of differentially overexpressed circRNA in bladder cancer group favoured platelet alpha granule lumen and platelet activation, compared to non-cancerous group, while the repressed transcripts were more related to actin binding and contractile fibre parts (Fig. 7). These results suggest that these pathways and molecular function-related circRNA may contribute to BCa pathogenesis and development.

Figure 7

GO analyses of the genes that produce differentially expressed circRNA. GO annotations of the linear counterparts of upregulated and downregulated circRNA that cover the domains of biological processes, cellular components and molecular functions. Both upregulated and downregulated circRNA are significantly changed with a fold change of >2.0 or < 0.5 and a q value of <0.001.

Validation of dysregulated circRNA

To validate the screened circRNA candidates, we designed multiple pairs of divergent primers to amplify the back-spliced junctions of circRNA in T24 and 293T cell lines. The amplified products were further confirmed via Sanger sequencing. Two downregulated circRNA (hsa_circ_0091017 and hsa_circ_0002024) were randomly detected in NGS data set via qRT-PCR in 32 pairs of samples. The levels of hsa_circ_0091017 and hsa_circ_0002024 in cancerous tissues were significantly reduced compared with non-cancerous tissues (Fig. 8).

Figure 8

Relative levels of hsa_circ_0002024 and hsa_circ_0091017 in 32 bladder cancer samples. Each point represents the mean result of triplicate samples.

Discussion

As a new member of the ceRNA network, circRNA has drawn extensive attention recently. circRNA are abundant, conserved, stable and cell type-specific in the regulation of cell function 19. Recently, the role of circRNA in tumourigenesis has been demosntrated 22,23. In this study, we detected circRNA in linear and ribosomal RNA-deprived RNA samples from bladder cancer and adjacent normal mucosa using RNA-seq and qPCR technology. The digestion of total RNA using exoribonuclease enzyme RNase R showed overwhelming enrichment for the identified head-to-tail back-spliced junctions compared with linear mRNA, thereby demonstrating the capability of our detection pipeline to recognise legitimate circRNA.

NGS technology has been widely used for mRNA, lncRNA and miRNA. However, the application of NGS in circRNA has only begun. Zheng et al. characterised circHIPK3 in multiple cancers using NGS technology 24. Ahmed et al. showed altered expression patterns of circRNA in epithelial ovarian carcinoma 9. However, systematic studies have not yet been conducted on circRNA profiling in human bladder cancer using NGS technology.

In all the five matched samples, we identified 14 upregulated and 42 downregulated circRNA in cancerous samples compared to non-cancerous samples. The majority of circRNA belonged to exonic, followed by intronic, intergenic and other functional elements (Fig. 3). These dysregulated circRNA are widely distributed in all chromosomes, including sex chromosomes X and Y, suggesting that all chromosomes are involved in the carcinogenesis of BCa.

CircRNAs such as ciRS-7 and Sry have been reported to modulate miRNA activity 5. Xie et al. reported that circRNA_001569 could promote the proliferation and invasion of colorectal cancer by targeting miR-145 25. Circ-ITCH functions as a cancer suppressor in oesophageal squamous cell carcinoma and colorectal cancer 11,12. In bladder cancer, Zhong et al. reported the regulatory role of circTCF25-miR-103a-3p/miR-107-CDK6 pathway 26. However, research on circRNA function is only in the initial phase and requires further exploration in the future.

Several studies have reported transcription regulation of circRNA 19,20. For example, ci-ankrd52, ci-mcm5 and ci-sirt7 enhanced the expression of their parental mRNA 27,28. GO enrichment analyses were performed to classify biological processes, molecular functions and cellular components from differentially expressed mRNA, which could identify the role of differentially regulated circRNA. Our data showed that upregulated mRNA was most associated with kinase activation in biological processes. In contrast, the downregulated transcripts were relevant to muscle contraction in biological processes and actin binding in molecular function. These data indicate that these upregulated cirRNAs would promote kinase activation and stimulate oncogenic signaling pathways to cause BCa.

In addition, KEGG pathway enrichment analysis for significantly differentially expressed mRNA indicated that 2 pathways are associated with upregulated mRNA, whereas 14 pathways are related to downregulated mRNA. The top enriched KEGG pathways in cellular processes for downregulated transcripts included ECM-receptor interaction and focal adhesion. It is well known that reduced ECM-receptor interaction and focal adhesion would impair cell adhesion and promote cell migration and invasion, leading to tumorigenesis and metastasis. Therefore, these pathways may contribute significantly to BCa pathogenesis and development.

In conclusion, this is the first report on aberrant circRNA profile of BCa and our results provide new insights into the role of circRNA in the tumorigenesis of BCa.