Comprehensive profiling of circular RNA expressions reveals potential diagnostic and prognostic biomarkers in multiple myeloma

Abstract

Background

This study aimed to explore the heterogeneity of circRNA expression pattern via microarray, and further evaluate the potential of 10 specific circRNAs as diagnostic and prognostic biomarkers in multiple myeloma (MM).

Methods

In exploration stage (stage I), circRNA expression profiles were detected by the microarray in bone marrow plasma cells from 4 MM patients and 4 healthy controls (HCs), and bioinformatic analyses were performed. In validation stage (stage II), top 10 upregulated and top 10 downregulated circRNAs identified in stage I were detected in 60 MM patients and 30 HCs for further validation; the diagnostic and prognostic values of these circRNAs in MM patients were analyzed.

Results

In stage I, 122 upregulated and 260 downregulated circRNAs were identified in MM patients compared with HCs. GO, KEGG and pathway enrichment analyses revealed that these circRNAs were implicated in neoplastic pathways such as MAPK and VEGF signaling pathways. In stage II, circ-PTK2, circ-RNF217, circ-RERE, circ-NAGPA and circ-KCNQ5 were validated to be upregulated and circ-AFF2, circ-WWC3, circ-DNAJC5, circ-KLHL2, circ-IQGAP1 and circ-AL137655 were validated to be downregulated in MM compared with controls. Circ-PTK2 and circ-RNF217 were correlated with poor treatment response and survival, while circ-AFF2 predicted good treatment response and survival in MM patients.

Conclusions

This study provides valuable reference for profound understanding about circRNA expression patterns in MM, and validates that circ-PTK2, circ-RNF217 and circ-AFF2 might serve as potential prognostic biomarkers in MM.

Background

Multiple myeloma (MM) is the second most common hematological malignancy derived from long-lived antibody-producing plasma cells in the bone marrow and is characterized by the presence of monoclonal immunoglobins in the serum and/or urine [1]. Over the past half century, the introduction of novel drugs (such as bortezomib) and application of hematopoietic stem cell transplantation have turned the rapid lethal MM into a chronic and manageable disease with extended survival in most of the patients [2, 3]. However, MM lacks symptoms in early stage, and the identification of the disease onset is difficult to be achieved by current examinations [4]. Moreover, the obstacles in treatment such as relapse and multidrug resistance are still common, contributing to poor prognosis in MM patients [5]. Therefore, it is essential to explore novel biomarkers that would help with diagnosis and improve the survival in MM patients.

Circular RNAs (circRNAs) are a class of non-coding RNAs originated by back-splicing of the precursor messenger RNA and forming a covalent loop with no 5′ to 3′ polarity or polyadenylated tail [6]. CircRNAs are stable, abundant and evolutionarily conserved, and mounting studies have proven that they contain target sites for microRNAs (miRNAs), thereby participate in the pathogenesis of various diseases through disturbing miRNAs signal axis [7]. With the development of circRNA microarrays, our knowledge about circRNA expression patterns has been initially uncovered in various diseases. In cancer research, circRNA expression patterns have been studied in some solid tumors including breast cancer, esophageal squamous cell cancer, epithelial ovarian cancer, etc., and a number of circRNAs are disclosed to involve in the pathophysiological progression of these malignancies [8,9,10,11]. As for hematological malignancies, an extensive analysis of circRNA expression profiles reveals a total of 464 dysregulated circRNAs (147 upregulated and 317 downregulated) in acute myeloid leukemia (AML) patients compared with healthy controls, and among these circRNAs, circ_0004277 is validated to be positively associated with prognosis in AML patients [9]. Whereas in MM, the expression profiles of circRNAs are not yet reported. Concidering that circRNAs are differentially expressed and involve in the pathophysiological progression of solid tumors as well as hematological malignancies, we speculated that they might play critical roles in MM as well.

This present study aimed to investigate the heterogeneity of circRNA expression pattern via microarray, and further evaluate the potential of 10 specific circRNAs as diagnostic and prognostic biomarkers in MM.

Methods

Study design

This study consisted of two stages (Fig. 1). In stage I (Exploration Stage), bone marrow samples were collected from 4 MM patients and 4 healthy controls (HCs), and plasma cells were isolated. Then circRNA expression profiles were detected by the microarray, and the bioinformatic analysis of circRNA microarray was performed. CircRNAs with at least 50% samples flagged as “Detected” were selected for further analysis, among which, circRNAs with ≥2.0 fold-changes (FC) and adjusted P values < 0.05 were defined as circRNAs with significant differential expression. In the stage II (Validation Stage), top 10 upregulated and top 10 downregulated circRNAs (based on rank of absolute value for log₂FC) were selected from dysregulated circRNAs identified in the stage I, then were determined by quantitative polymerase chain reaction (qPCR) in the 60 MM patients (including the 4 MM patients in the stage I) and 30 HCs (including the 4 HCs in the stage I) for validation, and the diagnostic and prognostic value of these circRNAs in MM patients were further analyzed.

Fig. 1

Study Flow. MM, multiple myeloma; HCs, healthy controls; PCA, principal component analysis; RT-qPCR, Real Time quantitative polymerase chain reaction

Participants

Between October 2015 and September 2018, 60 de novo MM patients and 30 HCs were consecutively recruited from the Shanghai Jing’an District Zhabei Central Hospital. The inclusion criteria for the MM patients were: (1) newly diagnosed as MM according to International Myeloma Working Group (IMWG) updated criteria for the diagnosis of multiple myeloma (2014); (2) age more than 18 years; (3) life expectancy more than 12 months; (4) able to be regularly followed up. Following MM patients were excluded: (1) relapsed or secondary MM; (2) history of stem cell transplantation (SCT), chemotherapy, radiotherapy or other systematic treatments before enrollment; (3) accompanied with other malignancies; (4) serious infection (e.g. Human Immunodeficiency Virus); (5) pregnant women or lactating women. Besides, all 30 enrolled HCs were healthy bone marrow donors, whose health status was confirmed before donation by appropriate examinations. This study was approved by the Institutional Review Board of Shanghai Jing’an District Zhabei Central Hospital and was conducted according to the Ethical Guidelines for Human Genome/Gene Research issued by the Chinese Government. All participants provided written informed consents before enrollment.

Collection of baseline data

Baseline data were collected after the patients signed the informed consents, including demographic information, such as age and gender, clinical characteristics and laboratory tests, such as immunoglobulin subtype, bone lesion, hemoglobin (Hb), calcium, serum creatinine (Scr), albumin (ALB), Beta-2-microglobulin (β2-MG), Durie-Salmon Stage, the International Staging System (ISS) Stage, lactate dehydrogenase (LDH) and cytogenetics abnormality. Durie-Salmon Stage and ISS Stage were evaluated in accordance with the Durie-Salmon Criteria and ISS Criteria respectively [12, 13], and cytogenetics abnormalities were determined by fluorescence in situ hybridization.

Collection and processing of samples

For enrolled MM patients, bone marrow samples were extracted and collected before any treatment; as for the HCs, bone marrow samples were obtained on the enrollment. Immediately after collection of bone marrow samples, separation of mononuclear cells was performed with gradient density centrifugation, then plasma cells were purified using CD138-coated magnetic beads (Miltenyi Biotec, Germany), and all operations were carried out in strict accordance with the manufacturer’s instructions to ensure greater than 90% plasma cell purity.

RNA extraction and quality control

Total RNAs were extracted from the plasma cells using the Trizol reagent (Invitrogen, USA), according to the manufacturer’s protocol, and RNA integrity was assessed using Agilent 2100 Bioanalyzer (Agilent, USA). Total RNA was quantified using NanoDrop ND-1000 spectrophotometer (Thermo, USA), then linear RNAs were diminished using RNase R (Epicentre, USA).

Microarray detection of circRNAs

After removing linear RNAs, 4 samples from MM patients and 4 samples from HCs were amplified and transcribed into fluorescent cRNA utilizing a random priming method with a Super RNA Labeling Kit (Arraystar, USA), and the labeled cRNAs were purified using RNeasy Mini Kit (Qiagen, Germany). Then, the samples were hybridized with a CapitalBio Technology Human CircRNA Array v1 (Agilent, USA) and Hybridization Kit (Agilent, USA) following the manufacturer’s standard protocols in an Agilent Hybridization Oven (Agilent, USA). After hybridization, the hybridized arrays were washed, fixed and scanned using an Agilent Microarray Scanner (Agilent, USA). Scanned images were imported into Agilent Feature Extraction software (version 11.0.1.1) to obtain raw data. Quantile normalization and low-intensity filtering were carried out with the use of R software package (R version 3.1.2). The circRNAs with at least 50% of samples flagged as “Detected” were selected for further analysis.

Bioinformatics analysis

Bioinformatics analysis was performed using R software package (R version 3.1.2). In brief, principal component analysis (PCA) of circRNA expression profiles was completed by Stats package; differentially expressed circRNAs between MM patients and HCs were analyzed with independent samples t-test using limma package, and circRNAs with a FC ≥2.0 and an adjusted P value (BH multiple test correction) < 0.05 were identified as differentially expressed circRNAs, which were displayed by volcano plots; heatmap plot of differentially expressed circRNAs were completed by pheatmap package. Gene Ontology (GO) enrichment analyses of dysregulated circRNAs were performed based on their located mRNAs and target miRNAs respectively; Kyoko Encyclopedia of Genes and Genomes (KEGG) enrichment analysis and pathway enrichment analysis of dysregulated circRNAs were performed based on their located mRNAs and predicted target miRNAs respectively. In order to investigate the regulation network between circRNAs and their target miRNAs, top 10 upregulated and top 10 downregulated circRNAs in MM sample were selected (based on rank of absolute value for log₂FC) to plot the circRNA-miRNA network using miRanda Database.

Validation of 20 candidate circRNAs by qPCR

A total of 60 MM patients’ samples and 30 HCs’ samples were used for qPCR validation. Top 10 upregulated and top 10 downregulated circRNAs were selected from differentially expressed circRNAs (identified in stage I) by the rank of the absolute value of Log₂FC and were determined by the qPCR, which was performed briefly as follows: after removing linear RNA using RNase R (Epicentre, USA), RNA was reverse transcribed into cDNA with random primer using PrimeScript™ RT reagent Kit (Perfect Real Time) (Takara, Japan) according to the manufacturer’s instructions. Then qPCR was carried out using TB Green™ Fast qPCR Mix (Takara, Japan). The circRNAs relative expression was calculated using 2^-△△Ct method and normalized to GAPDH. All of the quantitative PCR reactions were conducted in triplicate. The primers used in qPCR were listed in Table 1. The expressions of top 10 upregulated and top 10 downregulated circRNAs detected by qPCR between 4 MM patients and 4 HCs from Stage I were shown in Additional file 1: Table S1.

Table 1 Characteristics of MM patients in Stage I and Stage II respectively

Treatment and follow up

All MM patients received appropriate treatments based on disease conditions according to IMWG Recommendations for the Treatment of Multiple Myeloma–Related Bone Disease, and the treatment responses were evaluated referring to the IMWG criteria as well. Response was classified as complete response (CR), very good partial response (VGPR), partial response (PR), and the overall response rate (ORR) was calculated as CR + VGPR+PR. All MM patients were routinely followed up to 2018/12/31 with the median follow-up duration of 24.0 months (range: 5.0–36.0 months). Besides, progression free survival (PFS) was calculated from the date of initiation treatment to the date of disease progression or death; Overall survival (OS) was calculated from the date of initiation treatment to the date of death.

Statistical analysis

Data were displayed as mean and standard deviation (SD), median and interquartile range (IQR) or count (percentage). Comparisons were determined by the independent sample t test, Wilcoxon rank sum test or Chi-square test. Univariate and multivariate logistic regression analyses were performed to screen the circRNAs predicting MM risk. For the independent circRNAs in predicting MM risk, single and combined receiver operating characteristic (ROC) curves were plotted, and the area under the curve (AUC) was calculated to assess the diagnostic value of these circRNAs for MM. Survival profiles were displayed with Kaplan-Meier curves, and the difference in survival was determined by the log-rank test. All statistical analyses were performed using SPSS 24.0 Software (IBM, USA) or R software (Version 3.1.2), and graphs were plotted using GraphPad 7.01 Software (GraphPad, USA). P value < 0.05 was considered as significant.

Results

Baseline characteristics of MM patients in stage I and stage II respectively

In Stage I, 4 MM patients aged 64.5 ± 4.5 years with 2 male and 2 females were included for microarray assay (Table 2). The number of MM patients with Durie-Salmon stage I, II and III were 0 (0.0%), 2 (50.0%) and 2 (50.0%) respectively; and those in ISS stage I, II and III were 1 (25.0%), 2 (50.0%) and 1 (25.0%) respectively. In Stage II, 60 MM patients aged 60.0 ± 9.4 years with 37 males and 23 females were included for qPCR validation. There were 2 (3.3%), 32 (53.4%) and 26 (43.3%) patients in Durie-Salmon stage I, II and III respectively; and 13 (21.7%), 21 (35.0%) and 26 (43.3%) patients in ISS stage I, II and III respectively. See Table 2 for other detailed baseline information of MM patients in Stage I and Stage II.

Table 2 Top 10 upregulated and 10 downregulated circRNAs in MM patients compared to HCs

Identification of differentially expressed circRNAs in MM by microarray

PCA analysis showed clear segregation between 4 MM patients and 4 HCs, which indicated that circRNA expression patterns were able to distinguish MM patients from HCs (Fig. 2a). The valcano analysis was used to determine differentially expressed circRNAs between MM and HCs, which illustrated that 122 circRNAs were upregulated and 260 circRNAs were downregulated in MM compared with HCs (Fig. 2b). The threshold was set to fold change ≥2.0 and adjusted P value < 0.05. Following that, 122 upregulated and 260 downregulated circRNAs were included in heatmap analysis, and were shown to differentiate MM patients from HCs clearly (Fig. 2c).

Fig. 2

PCA, Valcano and Heatmap analyses of circRNA expression profiles. a MM patients and HCs were divided into two different groups by circRNAs. b 122 circRNAs were upregulated and 260 circRNAs were downregulated in MM compared with HCs. c The differentially expressed circRNAs were able to differentiate MM patients from HCs. MM, multiple myeloma; circRNAs, circular RNAs; HCs, healthy controls

GO and KEGG enrichment analysis based on located genes and target miRNAs of dysregulated circRNAs in MM

GO enrichment analysis by located genes revealed that the located genes of dysregulated circRNAs in MM were enriched in biological processes (e.g. positive regulation of cytoplasmic mRNA and cellular response to hypoxia), cellular components (e.g. cytosol and membrane), molecular functions (e.g. protein binding and protein kinase activity) (Fig. 3a). And from KEGG enrichment analysis, the located genes of dysregulated circRNAs in MM were enriched in pathways such as VEGF signaling pathway and MAPK signaling pathway, which were well-known neoplastic pathways (Fig. 3b). According to GO enrichment analysis by target miRNAs, the target miRNAs of dysregulated circRNAs in MM were enriched in biological processes (e.g. positive regulation of t cell mediated cytotoxicity and synaptonemal complex assembly), cellular components (e.g. high density lipoprotein particle and MHC class I protein complex), and molecular functions (e.g. adp ribose diphosphatase activities and endodeoxynuclease activity produci) (Fig. 3c). Regarding the pathway enrichment analysis by target miRNAs, the target miRNAs of dysregulated circRNAs in MM were disclosed to be enriched in pathways that underline the ontology of various malignancies such as malignant fibroxanthoma and carcinosarcoma (Fig. 3d).

Fig. 3

GO, KEGG and pathway enrichment analyses for circRNAs. The biological process, cellular component and molecular function that circRNAs were implicated in was revealed by GO enrichment analysis (a, c). In GO category, dot, triangle and square was used to symbolize BP, CC and MF respectively; the size of the symbols presented the gene counts; the gradation of color stood for the value of minus log10 adjusted P value. The cellular signaling pathways that circRNAs were implicated in were assessed by KEGG enrichment analysis (b) and pathway enrichment analysis (d). In KEGG by located genes, the size of the dot presented the gene counts; the gradation of color stood for the value of minus log10 P value. In pathway enrichment analysis by target miRNAs, category, dot, triangle and square was used to symbolize DO, HPO and KEGG respectively; the size of the symbols presented the gene counts; the gradation of color stood for the value of minus log10 P value. BP, biological process; CC, cellular component; MF, molecular function; DO, disease ontology; HPO, human phenotype ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; circRNAs, circular RNAs

Top 10 upregulated and top 10 downregulated circRNAs in MM patient compared with HCs by microarray and regulation network of these circRNAs

Top 10 upregulated and top 10 downregulated circRNAs detected by microarray in MM patients compared to HCs were selected by the rank of the absolute value of Log₂FC, and the detailed information of these circRNAs was listed in Table 3. Besides, the regulation network of these circRNAs with their target miRNAs was shown in Fig. 4.

Table 3 The primers used in qPCR

Fig. 4

Regulation network of top 10 upregulated and top 10 downregulated circRNAs by microarray. The regulation network of circRNAs was as shown. The red squares represented the top 10 upregulated circRNAs; the blue squares represented the top 10 downregulated circRNAs, the gray dots represented the target miRNAs of circRNAs. The circRNA-miRNA network was plotted using miRanda Database. CircRNAs, circular RNAs; miRNA, micro RNA

Expression of candidate circRNAs between MM patients and HCs

In validation stage, top 10 upregulated and top 10 downregulated circRNAs in MM patients compared to HCs were determined by the qPCR, and compared between MM patients (N = 60) and HCs (N = 30) for validation. Among the top 10 upregulated circRNAs, circ-PTK2 (P < 0.001) (Fig. 5a), circ-RNF217 (P = 0.008) (Fig. 5c), circ-RERE (P = 0.001) (Fig. 5d), circ-NAGPA (P < 0.001) (Fig. 5f) and circ-KCNQ5 (P = 0.002) (Fig. 5g) were validated to be upregulated, while circ-ATIC (P = 0.132) (Fig. 5b), circ-SETD5 (P = 0.329) (Fig. 5e), circ-CSPP1 (P = 0.340) (Fig. H), circ-SFMPT2 (P = 0.918) (Fig. 5i) and circ-UGGT2 (P = 0.221) (Fig. 5j) expression levels were similar in MM patients compared with HCs. As for the validation of the top 10 downregulated circRNAs, circ-AFF2 (P < 0.001) (Fig. 5f), circ-WWC3 (P < 0.001) (Fig. 5g), circ-DNAJC5 (P = 0.008) (Fig. 5p), circ-KLHL2 (P = 0.006) (Fig. 5q), circ-IQGAP1 (P < 0.001) (Fig. 5r) and circ-AL137655 (P = 0.001) (Fig. 5s) expressions were lower in MM patients compared with HCs, while circ-SLAIN1 (P = 0.146) (Fig. 5m), circ-WDR37–1 (P = 0.292) (Fig. 5n), circ-WDR37–2 (P = 0.199) (Fig. 5o) and circ-ASPA1 (P = 0.434) (Fig. 5t) levels were similar between MM patients and HCs.

Fig. 5

Comparing the levels of candidate circRNAs including Circ-PTK2 (a), circ-ATIC (b), circ-RNF217 (c), circ-RERE (d), circ-SETD5 (e), circ-NAGPA (f), circ-KCNQ5 (g), circ-CSPP1 (h), circ-SFMPT2 (i) and circ-UGGT2 (j), circ-AFF2 (k), circ-WWC3 (l), circ-SLAIN1 (m), circ-WDR37-1 (n), circ-WDR37-2 (o), circ-DNAJC5 (p), circ-KLHL2 (q), circ-IQGAP1 (r), circ-AL137655 (s) and circ-ASPA1 (t) in MM patients and HCs. Comparison of CircRNAs expressions between MM patients and HCs were determined by Wilcoxon rank sum test. P < 0.05 was considered significant. MM, multiple myeloma; circRNAs, circular RNAs

Correlation of candidate circRNAs with MM risk

The 5 upregulated (Fig. 6a) and 6 downregulated (Fig. 6b) circRNAs were then included in the ROC analysis, which illuminated that circ-PTK2 (AUC: 0.814, 95% CI: 0.729–0.900), circ-RNF217 (AUC: 0.672, 95% CI: 0.562–0.783), circ-RERE (AUC: 0.725, 95% CI: 0.620–0.830), circ-NAGPA (AUC: 0.804, 95% CI: 0.714–0.895), circ-KCNQ5 (AUC: 0.704, 95% CI: 0.594–0.815) could predict higher MM risk. And circ-AFF2 (AUC: 0.757, 95% CI: 0.641–0.872), circ-WWC3 (AUC: 0.773, 95% CI: 0.673–0.874), circ-DNAJC5 (AUC: 0.672, 95% CI: 0.557–0.787), circ-KLHL2 (AUC: 0.677, 95% CI: 0.564–0.790), circ-IQJAP1 (AUC: 0.758, 95% CI: 0.655–0.860), circ-AL137655 (AUC: 0.708, 95% CI: 0.601–0.816) could predict lower MM risk.

Fig. 6

ROC curves of 5 upregulated (a) and 6 downregulated (b) candidate circRNAs in distinguishing MM patients from HCs. ROC curve, receiver operating characteristics curve; MM, multiple myeloma; circRNAs, circular RNAs; HCs, healthy controls

Correlation of candidate circRNAs with clinical characteristics in MM patients

Among the top 10 upregulated circRNAs, circ-PTK2 was correlated with higher β2-MG level (P = 0.002), advanced ISS stage (P = 0.002) and deletion of 17p (P = 0.020); circ-RERE was associated with higher calcium concentration (P = 0.024) and advanced Durie-Salmon stage (P = 0.037); circ-SETD5 was positively correlated with deletion in 17p (P = 0.024); circ-KCNQ5 was positively correlated with ALB level (P = 0.012); circ-UGGT2 was negatively correlated with Durie-Salmon stage (P = 0.037) and positively correlated with deletion in 17p (P = 0.020) (Table 4). As for the top 10 downregulated circRNAs, circ-AFF2 correlated with lower β2-MG level (P = 0.002) and decreased ISS stage (P = 0.002); circ-WWC3 was associated with lower Durie-Salmon stage (P = 0.037); circ-WDR37–2 was negatively correlated with deletion in 17p (P = 0.016); circ-DNAJC5 was positively correlated with age (P = 0.039) and negatively correlated with LDH level (P = 0.028); circ-KLHL2 was negatively correlated with ALB level (P = 0.012); circ-IQGAP1 was correlated with abundance of IgA (P = 0.038) (Table 4). No correlation of candidate circRNAs with other clinical characteristics was observed, and the detailed information was listed in Table 4.

Table 4 Correlation of candidate circRNAs relative expression with clinical characteristics

Correlation of candidate circRNAs with treatment response in MM patients

The correlation of candidate circRNAs with treatment response in MM patients was assessed and we observed that, in the top 10 upregulated circRNAs, circ-PTK2 (P = 0.015) was associated with reduced CR; circ-RNF217 (P = 0.020) and circ-SETD5 (P = 0.029) were correlated with lower ORR (Table 5). As for the top 10 downregulated circRNAs, circ-AFF2 (P = 0.002) was positively correlated with CR. No correlation of other candidate circRNAs with treatment response was observed.

Table 5 Correlation of circRNAs relative expression with treatment response

Correlation of candidate circRNAs with survival profiles in MM patients

Circ-PTK2 (P = 0.035) (Fig. 7a), circ-RNF217 (P = 0.011) (Fig. 7c) and circ-DNAJC5 (P = 0.027) (Fig. 7p) were correlated with lower PFS, but circ-AFF2 (P = 0.003) (Fig. 7k) predicted longer PFS. Whereas the other candidate circRNAs, including circ-ATIC (P = 0.261) (Fig. 7b), circ-RERE (P = 0.277) (Fig. 7d), circ-SETD5 (P = 0.293) (Fig. 7e), circ-NAGPA (P = 0.541) (Fig. 7f), circ-KCNQ5 (P = 0.147) (Fig. 7g), circ-CSPP1 (P = 0.870) (Fig. 7h), circ-SFMBT2 (P = 0.251) (Fig. 7i), circ-UGGT2 (P = 0.351) (Fig. 7j), circ-WWC3 (P = 0.226) (Fig. 7l), circ-SLAIN1 (P = 0.919) (Fig. 7m), circ-WDR37–1 (P = 0.334) (Fig. 7n), circ-WDR37–2 (P = 0.468) (Fig. 7o), circ-KLHL2 (P = 0.823) (Fig. 7q), circ-IQGAP1 (P = 0.995) (Fig. 7r), circ-AL137655 (P = 0.082) (Fig. 7s) and circ-ASAP1 (P = 0.316) (Fig. 7t) were not correlated with PFS in MM patients.

Fig. 7

Correlation of candidate circRNAs including Circ-PTK2 (a), circ-ATIC (b), circ-RNF217 (c), circ-RERE (d), circ-SETD5 (e), circ-NAGPA (f), circ-KCNQ5 (g), circ-CSPP1 (h), circ-SFMPT2 (i), circ-UGGT2 (j), circ-AFF2 (k), circ-WWC3 (l), circ-SLAIN1 (m), circ-WDR37-1 (n), circ-WDR37-2 (o), circ-DNAJC5 (p), circ-KLHL2 (q), circ-IQGAP1 (r), circ-AL137655 (s) and circ-ASPA1 (t) with PFS in MM patients. PFS was displayed with Kaplan-Meier curves, and the difference in survival was determined by the log-rank test. P < 0.05 was considered significant. PFS, progression free survival; MM, multiple myeloma; circRNAs, circular RNAs

Regarding OS, circ-PTK2 (P = 0.004) (Fig. 8a) and circ-RNF217 (P = 0.022) (Fig. 8c) were associated with lower OS, but circ-AFF2 (P = 0.015) (Fig. 8k) was associated with longer PFS. Other candidate circRNAs including circ-ATIC (P = 0.823) (Fig. 8b), circ-RERE (P = 0.350) (Fig. 8d), circ-SETD5 (P = 0.460) (Fig. 8e), circ-NAGPA (P = 0.841) (Fig. 8f) circ-KCNQ5 (P = 0.219) (Fig. 8g), circ-CSPP1 (P = 0.301) (Fig. 8h), circ-SFMBT2 (P = 0.430) (Fig. 8i), circ-UGGT2 (P = 0.848) (Fig. 8j), circ-WWC3 (P = 0.760) (Fig. 8l), circ-SLAIN1 (P = 0.274) (Fig. 8m), circ-WDR37–1 (P = 0.485) (Fig. 8n), circ-WDR37–2 (P = 0.328) (Fig. 8o), circ-DNAJC5 (P = 0.228) (Fig. 8p), circ-KLHL2 (P = 0.889) (Fig. 8q), circ-IQGAP1 (P = 0.772) (Fig. 8r), circ-AL137655 (P = 0.085) (Fig. 8s) and circ-ASAP1 (P = 0.871) (Fig. 8t) were not correlated with OS in MM patients.

Fig. 8

Correlation of candidate circRNAs including Circ-PTK2 (a), circ-ATIC (b), circ-RNF217 (c), circ-RERE (d), circ-SETD5 (e), circ-NAGPA (f), circ-KCNQ5 (g), circ-CSPP1 (h), circ-SFMPT2 (i), circ-UGGT2 (j), circ-AFF2 (k), circ-WWC3 (l), circ-SLAIN1 (m), circ-WDR37-1 (n), circ-WDR37-2 (o), circ-DNAJC5 (p), circ-KLHL2 (q), circ-IQGAP1 (r), circ-AL137655 (s) and circ-ASPA1 (t) with OS in MM patients. OS was displayed with Kaplan-Meier curves, and the difference in survival was determined by the log-rank test. P < 0.05 was considered significant. OS, overall survival; MM, multiple myeloma; circRNAs, circular RNAs

Discussion

From this comprehensive analysis of circRNA expression profiles in MM, (1) we found that circRNA expression patterns were able to distinguish MM patients from HCs, and there were 122 upregulated and 260 downregulated circRNAs in MM compared with HCs, which were implicated in neoplastic signaling pathways such as MAPK signaling pathways and VEGF signaling pathway. (2) In validation stage, 5 out of 10 upregulated and 6 out of 10 downregulated circRNAs by microarray were confirmed by qPCR, and these circRNAs could distinguish MM patients from healthy controls. (3) Circ-PTK2 and circ-RNF217 were correlated with poor treatment response and survival, while circ-AFF2 predicted favorable treatment response and survival in MM patients.

With the rapid progression and wide application of high-throughput sequencing and microarray, the expression pattern of circRNAs in various human diseases are increasingly reported, and the dysregulated expressions of circRNAs are shown to contribute to the pathogenesis of various cancers [7, 9,10,11, 14,15,16,17,18,19]. For instance, in hepatitis B-related HCC, 189 upregulated circRNAs and 37 downregulated circRNAs were found by circRNA microarray, and circRNA_100,338 is further validated to be associated with metastatic progression by acting as an endogenous sponge for miR-141-3p [7]. Another study identifies 2556 upregulated and 1832 downregulated circRNAs in epithelial ovarian cancer tissues compared with normal ovarian tissues by microarray and bioinformatic analysis [11]. As for hematological malignancies, one previous study reveals the expression patterns of circRNAs in AML and exhibits 147 upregulated and 317 downregulated circRNAs in AML patients compared with healthy controls [10]. These previous studies uncover the expression patterns of circRNAs in several cancers including hematological malignancy, however, there is currently no study on heterogenicity of circRNA expression profiles in MM yet. In our study, we performed circRNA microarray and identified 122 upregulated and 260 downregulated circRNAs in bone marrow plasma cells of MM patients compared with HCs. In addition, these dysregulated circRNAs were shown to be involved in neoplastic signaling pathways including MAPK signaling pathway and VEGF signaling pathway. To our knowledge, this was the first study that investigates the expression patterns of circRNAs in MM, which might serve as valuable reference for further investigation of circRNAs functions in MM.

Benefiting from the stable nature and RNA degradation resistance, circRNAs are considered as prominent and novel biomarkers for many diseases, especially cancer, and there are specific circRNAs whose clinical values have been highlighted in several cancers. For instance, circ-LDLRAD3 is upregulated in pancreatic cancer, and is disclosed to be a potential biomarker in disease diagnosis [20]. In addition, circ_0014130 is positively correlated with TNM stage as well as lymphatic metastasis, and is of good diagnostic potential for NSCLC [21]. As for hematological malignancies, circ_0004277 is downregulated and offers a diagnostic biomarker in AML [10]. These aforementioned studies emphasized the potential of several specific circRNAs as biomarkers in diagnosis of solid tumors and hematological malignancies to a certain extent, however, the diagnostic value of circRNAs in MM is still misty. In addition, circRNA expression profile is a novel concept developed in recent years, and the current comprehensive screening of circRNA expression such as microarray is still limited by the accuracy. In addition, the sample size for microarray was far smaller than that of q-PCR, therefore, it was highly possible that the two analyses yielded deviation in results. CircRNA expression profiles by microarray aimed to give us a macroscopic view about the expression patterns of circRNAs, however, a larger sample size and more accurate tool were needed for a more refined understanding. Therefore, we carried out the Stage II analysis using q-PCR in a larger sample size. We selected the top 10 upregulated and top 10 downregulated circRNAs from the previous bioinformatic analyses, and validated regarding their diagnostic potential in MM with a larger sample size by qPCR. Our analyses revealed that, 5 out of 10 upregulated and 6 out of 10 downregulated circRNAs by microarray were confirmed by qPCR, and these circRNAs could distinguish MM patients from healthy controls. The possible explanations could be: (1) These circRNAs might influence the transcription of their parental genes by acting as restoration pools. For instance, circ-PTK2 serves as restoration pool for its paternal gene PTK2, which was identified as oncogene in MM, and increases the expression of PTK2 gene, thereby increases MM risk [22]. (2) These circRNAs might influence the pathogenesis of MM by sponging their target miRNAs. For example, circ-AFF2 might sponge miR-638 and inhibit the oncogenic function of miR-638 in MM (as shown in circRNA regulation network (Fig. 4)). Additionally, circ-PTK2 might act as sponge for anti-oncogenic miR-1298-5p and promote the neoplastic progression in MM (retrieved from tissue specific circRNA database: http://gb.whu.edu.cn/TSCD/). Although further studies were needed to analyze and demonstrate the detailed mechanisms of these circRNAs in MM, our study still illuminated that circ-PTK2, circ-RERE, circ-AFF2 and circ-WWC3 could serve as novel diagnostic biomarkers in MM.

Although rarely shown, it is still evident from the existing studies that some specific circRNAs are closely correlated with treatment response and may have potential prognostic value in cancer patients. For example, circ_0000285 expression is lower in cisplatin-resistant bladder cancer patients compared to cisplatin-sensitive patients and is independently correlated with poor treatment outcomes in bladder cancer patients [23]. Regarding survival, circ-RAD23B, an oncogene in NSCLC, predicts shorter OS in NSCLC patients [24]. Additionally, circRNA expression profiles display that circ_0001017 and circ_0061276 are correlated with longer OS in gastric cancer patients [25]. As for hematological malignancies, only one study exhibits that circ_100053 contributes to leukemogenesis in chronic myeloid leukemia (CML) and predicts increased resistance to imatinib as well as poor survival in CML patients [26]. Although the correlation of several specific circRNAs with patients’ prognosis in solid tumors as well as hematological malignancy has been reported, the correlation of circRNAs with prognosis in MM is still unknown [27]. In order to get a more profound understanding of the correlation of circRNAs with prognosis in MM, we evaluated the correlation of the top 10 upregulated and top 10 downregulated circRNAs with treatment response as well as survival, and disclosed that circ-PTK2 and circ-RNF217 were correlated with poor treatment response and survival, while circ-AFF2 predicted good treatment response and survival in MM patients. The possible reasons were: (1) According to our analyses, these circRNAs were closely correlated with clinicopathological features in MM patients, therefore, they would affect the prognosis of MM patients via influencing the clinicopathological features such as Durie-Salmon stage and deletion at 17p. (2) These circRNAs might change the cell sensitivity to chemotherapy and develop drug resistance via targeting miRNAs, thereby influence prognosis in cancer patients. For instance, circ-AFF2 might sponge miR-638, which was previously shown to induce drug resistance in human breast cancer, thereby reduced drug resistance and improved prognosis in MM patients [28]. (3) As explained above, these circRNAs might impact the normal function of miRNAs by serving as miRNA sponges in MM (see Fig. 4 for regulation network of candidate circRNAs and the detailed miRNAs), thereby influenced prognosis in MM patients. In addition, the tubular form of potential miRNA targets of the top 10 upregulated and top 10 downregulated circRNAs was shown in Additional file 2: Table S2, and the potential target miRNAs of all the 122 upregulated and 260 downregulated circRNAs were listed in Additional file 3.

This study first revealed the differential expressions of circRNAs and determined circRNAs with diagnostic and prognostic potential in MM, whereas there were still some shortcomings. Firstly, although several circRNAs with potential as diagnostic and prognostic biomarkers for MM were identified, the molecular mechanisms of these circRNAs in MM pathology were not investigated. Secondly, due to the budget, the sample size was relatively small for stage I, and MM patients from different clinical stage were not included for analysis, which could be improved in further studies. Thirdly, we explored the prognostic value of circRNAs but not in a logical approach because it was not the main goal in this study. However, it would be of great clinical significance to further detect the correlation of these circRNAs with prognosis in MM patients in a more logical way in the future. Moreover, the use of circRNAs as biomarkers for cancers is still in the early stage of research, and thorough practical proofs and standards were needed for clinical application. Studies that further validate the feasibility of circRNAs as diagnostic and prognostic biomarkers in cancer are needed to lead the bench side findings to real-life application.

Conclusion

In conclusion, this study provides valuable reference for profound understanding about expression patterns of circRNAs in MM, and validates that circ-PTK2, circ-RNF217 and circ-AFF2 might serve as potential prognostic biomarkers in MM.