Human BC tissues were selected from 43 patients (aged 50 to 68 years old) who received radical cystectomy between Oct 2015 and Dec 2017 at our institution (Shengjing Hospital of China Medical University, Shenyang, China) without preoperative neoadjuvant chemotherapy or radiotherapy. Three pairs of BC and matched ANNB tissues were used for circRNA microarray analysis, with another 40 pairs for subsequent validation. All samples were obtained within 5 min after resection and kept frozen in liquid nitrogen at −80°C. All tissue samples were pathologically evaluated by two well-experienced pathologists. The inclusion/exclusion criteria were as follows: i) pathologically confirmed bladder urothelial carcinoma; ii) no chemotherapy or radiotherapy before surgery; iii) no history of other malignant tumors; iv) complete clinicopathologic and follow-up data after surgery; and v) no evidence of distant metastasis at the time of surgery. All study participants provided their informed consent for inclusion in this study, and the protocol of the study was ratified and approved by the Research Ethics Committee of Shengjing Hospital of China Medical University (approval no. 2017PS012J).

A NanoDrop ND-1000 (Thermo Fisher Scientific, Inc.) was used to quantify total RNA from each sample. Briefly, circRNAs were enriched by eliminating linear RNAs with Rnase R (Epicentre Biotechnologies), followed by amplification and transcription into fluorescent circRNA using a random priming method (Arraystar Super RNA labeling kit; Arraystar, Inc.). The labeled cRNAs were hybridized onto the Arraystar Human circRNA Array V2 (8×15K; Arraystar, Inc.). After washing the slides, the hybridized arrays were washed by Arraystar Circular RNA microarray v 2.0. and then fixed on the chip scanner for signal value detection. Then, 1 µg of each labeled cRNA was fragmented by adding 5 µl 10X Blocking Agent and 1 µl of 25X fragmentation buffer, then heated the mixture at 60°C for 30 min. Subsequently, 25 µl of 2X hybridization buffer was added to dilute the labeled cRNA. Then, 50 µl of hybridization solution was dispensed into the gasket slide and assembled to the circRNA expression microarray slide. Finally, hybridized arrays were scanned using an Agilent Scanner G2505C (Agilent Technologies, Inc.).

The feature extraction software (v11.0.1.1; Agilent Technology, Inc.) was use to analyze the collected array images. The limma package from the R software package (https://www.r-project.org/) was used for quantile normalization and subsequent data processing. CircRNA exhibiting fold change (FC) ≥2 and a P<0.05 were considered significant. Volcano-plot (generated by R Software R-3.3.1 gplots) was performed to identify differentially expressed circRNAs with statistical significance between two groups and fold-change filtering was performed to identify differentially expressed circRNAs between two groups. Hierarchical clustering (produced by R software R-3.3.1 gplots, function heatmap2) was performed to display variable circRNA-expression patterns among samples. Data conversion is based on Z-Score. Arraystar miRNA target-prediction software utilizing the TargetScan (http://www.targetscan.org/) (19) and miRanda (http://www.miranda.org) (20) databases was used to investigate circRNA-miRNA relationships, with the top five predicted miRNAs for each differentially expressed circRNA extracted for further analysis. CircRNA-miRNA interactions were annotated in detail.

In order to be more representative of the chip as a whole, 40 sets of BC and ANNB tissues were randomly selected from the microarray for further verification by RT-qPCR analysis. Total RNA was extracted from each frozen sample (the amount of tissue was ~100 mg) in accordance with the manufacturer's protocol (RNAiso PLUS; Takara Bio, Inc.). The 1% agarose electrophoresis determined the integrity of the isolated RNA, and PrimeScript RT reagent kit with gDNA Eraser (Takara Bio, Inc.) was used to synthesize first-strand cDNA. A spectrophotometer was used to measure the optical density ratio at 260 and 280 nm (OD260/280=1.8–2.0) to determine RNA concentration. Primer 5.0 software (Premier Biosoft) was used to design specific primers, which are listed in Table SI. A StepOnePlus real-time PCR system with SYBR-Green I (Takara Bio, Inc.) was used to conduct RT-qPCR. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an internal control. Each sample was tested in triplicate. The thermocycling conditions used were as follows: 95°C for 30 sec, followed by 40 cycles at 95°C for 5 sec and 60°C for 30 sec. Relative mRNA expression levels were determined using the 2^−ΔΔCq method (21).

The mutant and wild-type sequences of hsa-circRNA-100876 were cloned downstream of the firefly luciferase gene pGL3 vector (Promega Corporation). Based on the instructions of the manufacturer, Luciferase Reporter plasmid (Promega Corporation) and miR-136-5p expression plasmid (pmirGLO) were co-transfected instantaneously utilizing Lipofectamine 2000 (Invitrogen; Thermo Fisher Scientific, Inc.). Subsequently, pGL3-hsa_circRNA_100876 vectors and miR-136-5p mimics along with negative control were co-transfected into 293 cells (Bena Culture Collection). Then, 48 h later, Dual-Luciferase reporter gene assay kit (Promega Corporation) was used to detect luciferase activity. The sequences of miR-136-5p mimics were provided by RiboBio Biotechnology Co., Ltd. and were as follows: forward, 5′-ACUCCAUUUGUUUUGAUGAUGGA-3′ and reverse, 5′-UCCAUCAUCAAAACAAAUGGAGU-3′. For each analysis, the Renilla luciferase signal was standardized to the firefly luciferase signal.

The BC cell lines EJ-1 (EJ) and T24 used in the present study were obtained from the Cell Bank of the Chinese Academy of Sciences in Shanghai. In addition, both cells lines were authenticated by STR profiles. Cells were grown at 37°C in a humidified atmosphere of 5% CO₂. All the cell lines were cultured in DMEM medium (Invitrogen; Thermo Fisher Scientific, Inc.) supplemented with 1% penicillin-streptomycin and 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.). Whether hsa_circRNA_100876 is involved in BC cell proliferation ability was measured using Cell Counting Kit-8 (CCK-8; Dojindo Molecular technologies, Inc.) assay. Cells in each group were digested by 0.25% trypsin and seeded into 96-well plates (5,000 cells/well); and three duplicate wells were set for each group of cells. CCK-8 (10 µl) reagent was added to each well and incubated at 37°C for 150 min after transfection for 24, 48, 72 and 96 h. All the experiments were performed in triplicate. The absorbance was detected at 450 nm using a microplate reader instrument (Sunrise™; Tecan Group, Ltd.). The ratio of absorbance difference between the experimental group and blank control group/24-h difference was the cell proliferation rate. The small interfering RNA (siRNA) targeting circRNA_100876 (si-circRNA_100876) was provided by RiboBio Biotechnology Co., Ltd. and was as follows: si-circRNA_100876, 5′-CTCCTACAATGTTGATATG-3′. si-NC (product no. siN0000001-1-5) was also provided by RiboBio Biotechnology Co., Ltd.

Gene Ontology (GO) enrichment analysis (22) was performed to determine functional annotations of miRNA target genes, mainly including three independent ontologies (molecular function, cellular component and biological process). TopGO (https://bioconductor.org/packages/release/bioc/html/topGO.html; version 2.32.0) was used to perform GO analysis of the differentially expressed genes to infer their involvement in molecular functions. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis (23) was used to clarify the interactions and functions among these differentially expressed genes. The R script was used to calculate the significance of differentially expressed genes and KEGG through hypergeometric distribution, and Fisher's exact test was used to calculate the P-value. Data downloads were based on the GO website (http://geneontology.org) and KEGG website (https://www.genome.jp/kegg/). The threshold for the P-value was <0.05 and the count number >2. The top 10 enriched GO items and KEGG pathways of the differentially expressed mRNAs were ranked by enrichment score [-log10 (P-value)].

The selected significantly expressed circRNAs, predicted miRNAs and protein-coding mRNAs were used to construct a circRNA-miRNA-mRNA regulatory network by a software based on TargetScan and miRanda (http://www.miranda.org). The software that was used considered the capacity and number of circRNA-microRNA binding sites, as well as the binding capacity of microRNA and mRNA.

Total proteins (~100 mg of tissue) were extracted from cell lysates with wash buffer (1X PBS, 0.1% SDS, 0.5% nonidet-P-40 and 0.5% sodium deoxycholate) containing 1% (w/v) protease inhibitor (Beyotime Institute of Biotechnology), and a BCA protein assay kit (Beyotime Institute of Biotechnology) was used to quantify the protein concentrations. Protein extract (80 µg) was loaded onto 10% SDS-PAGE gels and then transferred to polyvinylidene fluoride membranes (EMD Millipore). The membranes were then incubated with Tris-buffered saline containing 0.05% Tween-20 and 5% skim milk powder for 2 h at room temperature in the presence of the rabbit anti-CBX4 monoclonal antibody at 4°C overnight (dilution 1:1,000; cat. no. PA5-109482; Thermo Fisher Scientific, Inc.), followed by incubation for 120 min at room temperature with the corresponding secondary antibody HRP-conjugated Affinipure goat anti-rabbit IgG (H+L) (dilution 1:2,000; cat. no. SA00001-2, ProteinTech Group, Inc). GAPDH was used as an internal reference (1:10,000; cat. no. 10494-1-AP; ProteinTech Group, Inc.), and an enhanced chemiluminescence kit was used for visualization (Beyotime Institute of Biotechnology). The relative expression of the target protein was quantified by ImageJ software (version 1.48v; National Institutes of Health).

The 40 pairs of paraffin sections from BC tissues and ANNB tissues were prepared and analyzed. The sections (4-µm thick) were fixed with 10% neutral buffer formalin at room temperature for 12 h. The sections were randomly numbered and double blinded. Immunohistochemistry was performed using an SP kit (ZSGB-BIO; OriGene Technologies, Inc.). After rinsing with phosphate-buffered saline (PBS), 3% H₂O₂ was added at room temperature for 15 min. Then the sections were sealed with goat serum at room temperature for 15 min (ZSGB-BIO; OriGene Technologies, Inc), followed by antigen retrieval in a microwave (total power 800 W, time 10 min). An appropriate dilution of primary antibody (CBX4; 1:200; cat. no. sc-517216; Santa Cruz Biotechnology, Inc.) was then added and the sections were incubated overnight in wet box at 4°C. The sections were then washed and incubated for 30 min at 37°C with a secondary antibody [HRP-conjugated Affinipure goat anti-mouse IgG (H+L); cat. no. SA00001-1; ProteinTech Group, Inc], rinsed with PBS and then incubated with horseradish oxidase-labeled streptomycin avidin working solution (ZSGB-BIO; OriGene Technologies, Inc.) for 30 min at 37°C. Sections were then rinsed again in PBS and antibody binding was developed using a DAB chromogenic kit (Beyotime Institute of Biotechnology). Sections were then counterstained with hematoxylin at room temperature for 90 sec. A light microscope (original magnification, ×100; Nikon Corporation) was used to score the results of immunostaining, by multiplying the score of staining intensity of positive cells and the score of the percentage of positive cells. The final score is the degree of staining multiplied by the intensity.

Quantitative data are presented as the mean ± SEM. One-way ANOVA analysis and paired Student's t-test were used to assess significance among groups. Following one-way ANOVA, Tukey's post hoc test was performed. SPSS software (v19.0; SPSS, Inc.) was used to process all statistical analyses. A receiver operating characteristic (ROC) curve was generated to evaluate its diagnostic value and following calculation of a binomial exact-confidence interval to determine the area under the ROC curve (AUC). The associations between hsa_circRNA_100876 expression and prognosis of BC patients were evaluated using Kaplan-Meier survival analysis and log-rank tests. All tests were performed in triplicate. A P-value <0.05 was considered to indicate a statistically significant difference.

We selected three pairs of primary BC and ANNB tissues with similar clinicopathologic features to perform circRNA microarray for screening dysregulated circRNAs in BC. With a threshold of FC ≥2.0 and P<0.05, 512 differentially expressed circRNAs were identified following scanning and normalization, with 340 circRNAs significantly upregulated and 172 significantly downregulated. As revealed in Table SII, the top 10 regulated circRNAs were presented, along with other detailed information on P-value, FDR, strand, chromosome, circRNA type and gene symbols. The box plot indicated that the circRNA distribution among the six samples was nearly identical after normalization (Fig. 1A), and hierarchical clustering revealed clear variations in the expression profiles between BC and ANNB tissues (Fig. 1B). Differentially expressed circRNAs among samples were assessed by scatter-plot visualization (Fig. 1C) which revealed the variation of circRNA expression between BC tissues and control groups, and volcano map analysis (Fig. 1D) which was applied to visualize differentially expressed circRNAs.

Classification of the dysregulated circRNAs is presented in the bar diagram (Fig. 2A). The results identified 340 upregulated circRNAs comprising 285 exonic, 25 intronic, 19 sense-overlapping, eight antisense, and three intergenic regions, whereas the 172 downregulated circRNAs contained 148 exonic, 13 intronic, nine sense-overlapping, one intergenic, and one antisense in BC-tissue samples. Chromosomal distribution analysis indicated that most were located on chromosome 1, while few located on chromosome 21 and the Y chromosome (Fig. 2B).

Fifteen circRNAs (ten upregulated and five downregulated) were randomly selected from the microarray for further verification by RT-qPCR analysis from 40 sets of BC and ANNB tissues. Eleven of the 15 circRNAs were verified as significantly differentially expressed in cancer tissues, with eight upregulated and three downregulated (Fig. 3A and B). The results revealed upregulation of hsa_circRNA_100086, _038632, _100876, _102241, _004646, _102505, _104600, and _005411 expression in BC tissues relative to their levels in ANNB tissues, whereas hsa_circRNA_001409, _101969, and _103801 were downregulated. These findings were consistent with microarray results.

To investigate how circRNAs regulate the expression of a target gene, nine circRNAs (7 of the nine circRNAs were upregulated and 2 were downregulated) were randomly selected from the validated circRNAs for GO and KEGG pathway analyses. GO analysis annotated genes targeted by the nine differentially expressed circRNAs in domains of molecular functions, biological processes and cellular components. KEGG pathway analysis indicated that several pathways involved in cancer were related to the dysregulated circRNAs. The results were ranked by fold enrichment [(Count/Pop.Hits)/ (List.Total/ Pop.Total)] and revealed enrichment of multiple critical biological functions associated with BC progression (results for hsa_circRNA_100876 and _038632 are presented in Fig. 4A-D). In addition, the supplementary results of GO and pathway analyses are presented in the Figs. S1–S7, identifying multiple items associated with cancer progression. For example, hsa_circRNA_102241 was related to 51 pathways, of which the top 1 pathway was related with ‘proteoglycans in cancer’. As is known, some proteoglycans have pro- and anti-angiogenic activities, whereas other proteoglycans can also directly affect cancer growth by modulating key signaling pathways.

To visualize a circRNA-regulatory network, interactions between the six confirmed differentially expressed circRNAs and their association with miRNAs and target genes were predicted, followed by construction of a unified interaction-network model. The co-expression pattern of circRNA-miRNA-mRNA is presented in Fig. 5A. These results indicated potential roles for the identified circRNAs as endogenous RNAs capable of altering target gene expression. In addition, a ceRNA network of hsa_circRNA_100876-miRNA-mRNA is presented in Fig. 5B. A total of 4 miRNAs (miR-136-5p, miR-217, miR-665 and miR-636) and corresponding target mRNAs were predicted to have an interaction with hsa_circRNA_100876 in the present study. For each of the four identified interacting miRNAs, the intersection of predictions from TargetScan and miRanda were used to find their targeting protein coding genes. A total of 140 miRNA-regulating protein coding genes were identified for hsa_circRNA_100876. It was revealed that hsa_circRNA_100876 may be highly correlated and co-expressed with mRNAs. This provides us with another research strategy to explore the mechanism of hsa_circRNA_100876 by identifying its associated mRNAs and investigating whether it can play a role by regulating the expression of certain associated mRNAs.

To determine possible roles of the dysregulated circRNAs in BC carcinogenesis, hsa_circRNA_100876 was selected for further study. The results revealed that hsa_circRNA_100876 was located on chromosome 11 (71668272–71671937) and differentially overexpressed in BC tissues relative to ANNB tissues (Fig. 6A). ROC curve analysis revealed an AUC value of 0.7003 (Fig. 6B), with a specificity of 52.5% and sensitivity of 90%, suggesting the potential diagnostic efficacy of hsa_circRNA_100876 as a BC biomarker.

To confirm the clinical value of these molecular differences, analyses were performed to evaluate associations between hsa_circRNA_100876 expression and clinicopathologic features. As revealed in Table I, increased levels of hsa_circRNA_100876 were significantly associated with T stage and lymphatic metastasis in BC patients (P<0.05), although the expression level was not related to sex, age, histological grade, or recurrence. Furthermore, univariate analysis indicated that the hsa_circRNA_100876 expression, T stage, and lymphatic metastasis were significantly related to the overall survival (OS) of patients (P<0.05), and multivariate analysis revealed hsa_circRNA_100876 expression and lymphatic metastasis as independent factors affecting BC patient prognosis (Table II). Moreover, Kaplan-Meier survival analysis used to estimate the association between the hsa_circRNA_100876 patient prognosis post-surgery indicated that increased hsa_circRNA_100876 expression was correlated with markedly shorter OS duration (Fig. 6C).

To confirm the role of hsa_circRNA_100876 in cell proliferation, specific siRNA targeting hsa_circRNA_100876 were transfected into BC cells (EJ and T24), resulting in a significant decrease in hsa_circRNA_100876 expression (Fig. 7A and B; P<0.05). The results of CCK-8 assays revealed that, compared to the si-NC group, the proliferation abilities of EJ and T24 cells transfected with si-hsa_circRNA_100876 were significantly reduced (Fig. 7C and D; P<0.05), revealing that knockdown of hsa_circRNA_100876 suppressed proliferation of BC cells.

We constructed a network of circRNA-miRNA-mRNA interactions to visualize potential relationships based on circRNA microarray data. Prediction of miRNA targets revealed that miR-136-5p shared binding sites with hsa_circRNA_100876 (Fig. 8A), and the binding sites were further validated by luciferase reporter assay (Fig. 8C). Furthermore, the same complementary binding sites were revealed between miR-136-5p and CBX4 (Fig. 8A), which was further validated through luciferase reporter assay (Fig. 8E). Moreover, RT-qPCR analysis confirmed significant downregulation of miR-136-5p (Fig. 8B) in BC tissues relative to ANNB tissues, along with its negative correlation with hsa_circRNA_100876 expression (r=−0.7842, P<0.001) (Fig. 8F). Furthermore, it was observed that CBX4 levels were aberrantly upregulated in BC tissues (Fig. 8D), revealing a possible positive correlation with hsa_circRNA_100876 expression and a negative correlation with miR-136-5p expression (r=0.7619 and −0.6404, P<0.001, respectively) (Fig. 8G and H). Western blot analysis indicated increased CBX4 levels in BC tissues as compared with ANNB tissues, which was consistent with the RNA levels (Fig. 9A). Additionally, immunohistochemical staining results confirmed the increased CBX4 levels in BC tissues and relatively low levels in ANNB tissues (Fig. 9B).