This article is mentioned in:

Introduction

Long noncoding RNAs (lncRNAs) are a class of RNA (>200 nucleotides in length) that cannot synthesize proteins (1–4). These biomolecules are involved in post-transcriptional regulation (5–8), and are abnormally expressed in multiple types of solid tumor and hematopoietic malignancy; lncRNAs are involved in both carcinogenesis and tumor suppression (9–13).

The expression levels of several lncRNAs have been reported in multiple myeloma (MM); their clinical significance, biological functions and potential molecular mechanisms in the disease have also been investigated (14–16). MM is the second most frequent hematological malignancy, and accounts for ~10% of all such malignancies (17–22). Immunomodulatory drugs (such as lenalidomide and pomamide), proteasome inhibitors (such as bortezomib and carfilzomib) and monoclonal antibodies have significantly increased the survival rate of patients with MM over the past decade (23–25); however, the treatment of relapsed and partially refractory patients remains challenging. The pathogenesis and progression of MM involve complex and heterogeneous genomic alterations (26–30), including modifications that are influenced by lncRNAs.

Certain lncRNAs have been documented to serve an important role in the progression of MM, and can be used as indicators of patient prognosis. For example, metastasis-associated lung adenocarcinoma transcript 1 (MALAT1) is overexpressed in MM tissues and various MM cell lines; upregulation of MALAT1 is significantly associated with poor prognosis, including overall survival (OS) and progression-free survival (PFS) (31–33). Nuclear paraspeckle assembly transcript 1 (NEAT1) has also been reported to serve a pivotal role in promoting MM, and its elevated expression is closely associated with poor prognosis (34,35). The upregulation of urothelial cancer associated 1 (UCA1) (36), protein disulfide isomerase family A member 3 pseudogene 1 (PDIA3P) (37), H19 (38), colon cancer associated transcript 1 (CCAT1) (39) and colorectal neoplasia differentially expressed (CRNDE) (40) are closely associated with poor prognosis in MM; these genes may be used as future indicators in the clinical prognosis of patients with MM. Despite the large numbers of lncRNAs, only a small number have been associated with the prognosis of MM; however, numerous as-yet-undiscovered lncRNAs may also be associated with the progression of MM and patient outcome. Additionally, the predictive ability of a single indicator is limited; a prognostic signature composed of numerous indicators is required to conduct a comprehensive clinical evaluation of tumor prognosis. Prognostic models that combine several indicators have been used in a wide variety of tumors (41–45); however, a prognostic model for MM comprising lncRNAs is yet to be reported.

The present study screened gene chips with expression data from patients with MM and selected prognostic lncRNAs and mRNAs. The associations between the prognostic lncRNAs and mRNAs were mapped, and certain indicators were selected to construct a prognostic model. The MM prognostic model presented in the current study may provide novel insight and directions for the clinical treatment of MM in the future.

Materials and methods

Data acquisition

The microarray gene expression profiling data from the bone marrow of newly diagnosed patients with MM that had not been treated was obtained from the Genome Expression Omnibus (GEO) dataset (46,47) with accession number GSE24080 (48). The data from 559 patients with MM were included for further survival analysis. To separate lncRNAs and mRNAs, probes from the Affymetrix HG-U133_Plus_2.0 array were re-annotated. For genes that matched >1 probe, the expression values of all the measurements were calculated using an average value of the probes. The lncRNAs were extracted according to their Refseq database label (Release 93) (49) and Ensembl annotations (Release version 96) (50).

Survival analysis

Event-free survival (EFS) generally provides more reliable endpoint information for survival analysis (51); thus, it was selected as the survival analysis endpoint in the present study. Univariate Cox analysis was conducted to select prognosis-associated mRNAs and lncRNA using the survival package (version 2.44–1.1) in R (version 3.4.4) (R). P<0.005 was considered to be statistically significant (52). Kaplan-Meier plot was generated to observe the survival status between different survival associated mRNA and lncRNA expression levels.

Gene functional enrichment analyses

To further investigate the potential molecular mechanisms of the top 20 prognosis-associated mRNAs, the biological processes, which were acquired from gene ontology and Kyoto Encyclopedia of Genes and Genomes (Release 87.1) (KEGG) pathways (53–55) were examined based on enrichment analysis using the Clusterprofiler package (version 3.10.1) in R (56). Protein-protein interaction (PPI) networks were developed to explore the associations between each gene using the GeneMANIA plug-in in Cytoscape version 3.6.1 (57,58).

Prognostic signature construction

As the prognostic value of a single indicator is limited, prognostic signatures were produced that combined multiple indicator candidates. Multivariate Cox analysis was performed on the top 20 most significant prognosis-associated mRNAs and lncRNAs to develop prognostic signatures. The performances of these prognostic signatures were tested using the survivalROC (version 1.0.3) package in R, which provides time-dependent receiver operator characteristic (ROC) curve estimation (59,60). The area under curve (AUC) was calculated at 75 months, as fewer events occurred after this point.

Weighted correlation network analysis (WGCNA)

As lncRNAs cannot be transcribed into proteins, their functional effects are frequently achieved by targeting mRNAs. To investigate the associations between lncRNAs and mRNAs, WGCNA was conducted using the WGCNA package (version 1.63) in R (61,62). The mRNAs were separated into modules, and correlations between the prognostic mRNAs and lncRNAs were calculated. A lncRNA-mRNA axis was identified when an association coefficient >0.4 was obtained. The potential regulatory network was constructed using Cytoscape software.

Results

Prognosis-associated lncRNAs and mRNAs

The present study included 559 patients with MM from the GSE24080 (48) dataset. This dataset was collected and distributed by the Myeloma Institute for Research and Therapy at the University of Arkansas for Medical Sciences. Dichotomized OS and EFS values were examined based on a 2-year milestone cutoff. A univariate Cox regression analysis identified 39 lncRNAs and 1,445 mRNAs that were significantly associated with the EFS of patients with MM (Fig. 1A and B). The top 20 most significant survival-associated lncRNAs and mRNAs are presented in the form of forest plots (Fig. 1C and D; Table I).

Figure 1. Survival-associated lncRNAs and mRNAs in multiple myeloma. (A) Red dots indicate lncRNAs whose expression levels are significantly associated with patient survival, whereas blue dots indicate lncRNAs that do not exhibit an association. P<0.005 was set as the threshold. (B) Red dots indicate mRNAs whose expression levels are significantly associated with patient survival, whereas blue dots indicate mRNAs that do not exhibit an association. P<0.005 was set as the threshold. (C) Top 20 most significantly survival-associated lncRNAs. (D) Top 20 most significant survival-associated mRNAs. LncRNA, long noncoding RNA.

Table I. Top 20 most significant survival-associated mRNAs and lncRNAs.

Table I.

Top 20 most significant survival-associated mRNAs and lncRNAs.

A, mRNAs
Gene symbol	HR	Z-score	P-value
KIF14	1.559189	6.768521	1.30×10−11
FAM72A	1.425627	6.586882	4.49×10−11
CENPL	1.938061	6.233485	4.56×10−10
NEK2	1.468169	6.129078	8.84×10−10
IFI16	1.926689	6.126442	8.99×10−10
DTL	1.443987	6.061892	1.35×10−09
NUF2	1.330607	6.004863	1.91×10−09
SMC4	2.003832	5.957470	2.56×10−09
TPX2	1.458392	5.867954	4.41×10−09
UBE2T	1.469349	5.864287	4.51×10−09
PDE4A	0.548218	−5.729731	1.01×10−08
ABCB10	1.828576	5.704206	1.17×10−08
TIPRL	2.179489	5.665706	1.46×10−08
REEP5	0.500858	−5.643577	1.67×10−08
RBBP8	2.034743	5.624249	1.86×10−08
TOPBP1	2.324429	5.605628	2.08×10−08
MSH2	1.537154	5.602649	2.11×10−08
ANP32E	1.843312	5.600536	2.14×10−08
MCM2	1.529821	5.596079	2.19×10−08
THUMPD2	2.218271	5.567595	2.58×10−08
B, lncRNAs
Gene symbol	HR	Z-score	P-value
NCRNA00201	1.541349	5.562407	2.66×10−08
RP11-164P12.4	0.781350	−4.785449	1.71×10−06
AC116904.1	1.823522	4.696971	2.64×10−06
LOC282997	0.692999	−4.174333	2.99×10−05
HCG26	0.734874	−3.843160	1.21×10−04
CTD-2003C8.1	0.592561	−3.689247	2.25×10−04
RP11-18H21.1	0.832661	−3.620982	2.93×10−04
RP11-875O11.1	0.784430	−3.592916	3.27×10−04
AC022087.1	1.410099	3.472571	5.15×10−04
C9orf130	0.679642	−3.442218	5.77×10−04
RP11-217B7.2	0.765044	−3.388514	7.03×10−04
A1BG-AS	0.783840	−3.372035	7.46×10−04
C21orf34	1.340267	3.260104	1.11×10−03
RP13-15E13.1	0.698713	−3.243736	1.18×10−03
AC073548.1	0.720060	−3.176108	1.49×10−03
AC004383.4	1.289903	3.172203	1.51×10−03
AL356534.1	0.840370	−3.170147	1.52×10−03
RP11-557H15.4	0.799831071	−3.165730	1.55×10−03
CTC-454M9.1	0.731592407	−3.144922	1.66×10−03
RP11-706O15.5	1.263062805	3.144145	1.67×10−03

[i] HR, hazard ratio; lncRNA, long noncoding RNA.

Functional enrichment analyses

Enrichment analyses for the top 20 prognosis-associated mRNAs were conducted to identify risk pathways and biological functions associated with these prognostic genes (Table II). As presented in Fig. 2A, a number of cell cycle-associated biological processes were identified, including signal transduction by a p53-class mediator, cell cycle G2/M transition and mitotic cell cycle G2/M transition. These categories are also closely involved in tumor proliferation. The PPI network revealed that these genes were closely associated with each other (Fig. 2B). It was also observed that these prognostic genes were associated with homologous recombination (Table II). Collectively, the present results suggested that these genes may serve an important role in the pathogenesis of MM.

Figure 2. Biological processes and protein-protein interactions of the top 20 survival-associated mRNAs in multiple myeloma. (A) Biological processes in gene ontology. (B) Interactions between the top 20 most significant survival-associated mRNAs. Purple, blue, red and green connections indicate co-expression, co-localization, physical interaction and pathway, respectively.

Table II. The top 10 most significant biological processes and Kyoto Encyclopedia of Genes and Genomes pathway.

Table II.

The top 10 most significant biological processes and Kyoto Encyclopedia of Genes and Genomes pathway.

Category	ID	Description	P-value	Q-value	Genes	Count
Biological process	GO:0000077	DNA damage checkpoint	5.56×10−07	0.000101	DTL, TIPRL, RBBP8, TOPBP1, MSH2	5
Biological process	GO:0031570	DNA integrity checkpoint	7.35×10−07	0.000101	DTL, TIPRL, RBBP8, TOPBP1, MSH2	5
Biological process	GO:0010389	Regulation of G2/M transition of mitotic cell cycle	1.68×10−06	0.000148	KIF14, NEK2, DTL, TPX2, TOPBP1	5
Biological process	GO:1902749	Regulation of cell cycle G2/M phase transition	2.16×10−06	0.000148	KIF14, NEK2, DTL, TPX2, TOPBP1	5
Biological process	GO:0000075	Cell cycle checkpoint	3.49×10−06	0.000178	DTL, TIPRL, RBBP8, TOPBP1, MSH2	5
Biological process	GO:0000819	Sister chromatid segregation	3.88×10−06	0.000178	KIF14, CENPL, NEK2, NUF2, SMC4	5
Biological process	GO:0000086	G2/M transition of mitotic cell cycle	5.61×10−06	0.00022	KIF14, NEK2, DTL, TPX2, TOPBP1	5
Biological process	GO:0044839	Cell cycle G2/M phase transition	6.67×10−06	0.000229	KIF14, NEK2, DTL, TPX2, TOPBP1	5
Biological process	GO:0031572	G2 DNA damage checkpoint	7.73×10−06	0.000233	DTL, RBBP8, TOPBP1	3
Biological process	GO:0072331	Signal transduction by p53 class mediator	8.48×10−06	0.000233	IFI16, TPX2, RBBP8, TOPBP1, MSH2	5
KEGG pathway	hsa03440	Homologous recombination	8.06×10−04	0.007636	RBBP8, TOPBP1	2

Development of the prognostic signatures

The top 20 most significant survival-associated lncRNAs (Fig. 3) and mRNAs (Fig. 4) were selected as candidates for analyzing independent MM prognostic factors. Multivariate Cox regression analyses were conducted to develop two prognostic signatures based on lncRNAs and mRNAs (Fig. 5A and B). Both signatures could be used to separate patients into two groups with distinct outcomes. The AUCs of ROC curves were 0.739 for the lncRNA signature and 0.732 for the mRNA signature (Fig. 5C and D). These findings suggested that the two risk scores exhibited a moderate power to predict the survival of patients with MM.

Figure 3. Kaplan-Meier analysis of the top 20 most significantly survival-associated long noncoding RNAs in multiple myeloma. (A) NCRNA00201. (B) RP11-164P12.4. (C) AC116904.1. (D) LOC282997. (E) HCG26. (F) CTD-2003C8.1. (G) RP11-18H21.1. (H) RP11-875O11.1. (I) AC022087.1. (J) C9orf130. (K) RP11-217B7.2. (L) A1BG-AS. (M) C21orf34. (N) RP13-15E13.1. (O) AC073548.1. (P) AC004383.4 (Q) AL356534.1. (R) RP11-557H15.4. (S) CTC-454M9.1. (T) RP11-706O15.5.

Figure 4. Kaplan-Meier analysis of the top 20 most significantly survival-associated mRNAs in multiple myeloma. (A) KIF14. (B) FAM72A. (C) CENPL. (D) NEK2. (E) IFI16. (F) DTL. (G) NUF2. (H) SMC4. (I) TPX2. (J) UBE2T. (K) PDE4A. (L) ABCB10. (M) TIPRL. (N) REEP5. (O) RBBP8. (P) TOPBP1. (Q) MSH2. (R) ANP32E. (S) MCM2. (T) THUMPD2.

Figure 5. Prognostic signatures based on lncRNAs and mRNAs. (A) Kaplan-Meier analysis of the lncRNA-based risk score model predicts MM EFS. (B) Kaplan-Meier analysis of the mRNA-based risk score model predicts MM EFS. (C) ROC curve of the lncRNA-based risk score. (D) ROC curve of the mRNA-based risk score. EFS, event-free survival; lncRNA, long noncoding RNA; MM, multiple myeloma; ROC, receiver operating characteristic.

Construction of the lncRNAs-mRNAs network

WGCNA was used to separate prognosis-associated mRNAs into groups, and to explore the associations between lncRNAs and mRNAs (Fig. 6). The lncRNA-mRNA network provided novel insight into the regulatory mechanisms associated with the progression of MM (Fig. 7). A total of 143 mRNAs were positively or negatively associated with 23 prognosis-related lncRNAs. NCRNA00201, LOC115110 and RP5-968J1.1 appeared to be the most dominant drivers, as they possessed the highest number of connected genes.

Figure 6. Weighted correlation network analysis. (A) Survival-associated genes in multiple myeloma were divided into modules. (B) Relationships between the long noncoding RNAs and mRNAs.

Figure 7. LncRNA-mRNA regulatory network. Red circles indicate risk-associated mRNAs; blue circles indicate protective mRNAs. Red squares indicate risk-associated lncRNAs; blue squares indicate protective lncRNAs. Red edges indicate positive associations; blue edges indicate negative associations. LncRNA, long noncoding RNA.

Discussion

Currently, there is no precise method to assess the prognosis of patients with MM. In the present study, the expression data from a gene chip containing genomic samples from 559 patients with MM were analyzed, and predictive models were constructed based on the lncRNA and mRNA expression profiles. Of note, it was observed that pools of prognostic candidates exhibited greater predictive power than individual indicators. There may also be a targeting relationship between the prognosis-associated lncRNAs and mRNAs. As a previous report has contraindicated the use of a lncRNA prediction model for the prognosis of MM, the present study provides novel insight for the clinical diagnosis and treatment of MM (63).

Previously, two other research groups have analyzed the gene chip data of GSE24080 to obtain MM prognosis-associated lncRNAs using different statistical methods. Zhou et al (64) randomly split the MM cohort into a training dataset (n=280) and a testing dataset (n=279). Univariate regression analysis identified 59 lncRNAs that were associated with the OS of patients. Only four of those lncRNAs (RP4-803J11.2, RP1-43E13.2, RP11-553L6.5 and ZFY-AS1) were reported to exhibit a predictive effect following multivariate regression analysis. These results were inconsistent with a study by Hu et al (63), which identified 176 lncRNAs from the GSE24080 and GSE57317 datasets that appeared to be associated with patient survival. The Hu et al study employed Kaplan-Meier analysis to determine the prognostic influence of lncRNAs, identifying 176 lncRNAs, including RP1-286D6.1, AC008875.2, MTMR9L, AC069360.2 and AL512791.1, as prospective markers for assessing the prognosis of patients with MM. Of note, none of the aforementioned lncRNAs overlapped with the top 20 prognosis-associated lncRNAs identified in the present study, which conducted survival analysis using the survival package in R. Among the top 20 prognosis-associated lncRNAs in the present study, a number were identified as risk factors, including NCRNA00201, AC116904.1, AC022087.1, C21ORF34, AC004383.4 and RP11-706O15.5. The remaining 14 lncRNAs may protect against MM. The use of different statistical tools may partially explain the variations in the lncRNAs identified in each study. In the previously published studies that analyzed GSE24080, Hu et al (63) did not conduct ROC analysis to determine the prognostic values of their lncRNAs, whereas Zhou et al (64) reported four lncRNAs (RP4-803J11.2, RP1-43E13.2, RP11-553L6.5 and ZFY-AS1) together generated an AUC of 0.682 to represent prognostic performance. In the present study, ROC analysis was performed using the survivalROC package in R, resulting in an AUC of 0.739, more favorable than that in Zhou et al (64).

None of the lncRNAs in the presently reported prognostic model have been previously investigated in MM. At present, the majority of the top 20 prognosis-associated lncRNAs have not been reported in any disease. Only three of the lncRNAs have been previously studied, NCRNA00201, HCG26 and C21ORF34.

NCRNA00201, also termed heterogeneous nuclear ribonucleoprotein U processed transcript, has been studied in patients with a submicroscopic deletion at 1q43q44. NCRNA00201 is in a critical deleted region of 1q44 (65). It is expressed in various normal tissues, including brain tissues; the highest expression levels of NCRNA00201 have been observed in the cerebellum (65). As patients with a submicroscopic deletion at 1q43q44 present with intellectual disability, microcephaly, craniofacial anomalies, seizures, limb anomalies and corpus callosum abnormalities, NCRNA00201 may serve a role in these diseases; however, 9 out of 11 patients with a submicroscopic deletion at 1q43q44 did not present with microcephaly or corpus callosum abnormalities, and possessed a small deletion containing NCRNA00201 (65). This study revealed that NCRNA00201 does not affect microcephaly or corpus callosum abnormalities, but that it is a prognostic candidate for intellectual disability and seizures (65).

The potential biological functions of NCRNA00201 have also been documented in pancreatic ductal adenocarcinoma (PDAC); it was significantly overexpressed in PDAC tissues and cell lines compared with noncancerous pancreatic controls. Upregulated levels of NCRNA00201 were associated with poorer prognosis in patients with PDAC. Suppression of NCRNA00201 inhibited cell proliferation, invasion and migration in PDAC cell lines (66). These results indicated that NCRNA00201 served an important role in the tumorigenesis and progression of PDAC (66). NCRNA00201 may serve a similar role in MM, as increased NCRNA00201 levels were closely associated with poor survival in patients with MM. NCRNA00201 may be a risk factor for poor prognosis; however, this finding should be validated in additional cohorts.

The MM prognosis-associated lncRNA HCG26 has also been reported in other diseases. This lncRNA exhibited dysregulated expression in the blood of patients that had suffered an ischemic stroke (67). HCG26 was also reported to be associated with polycystic ovary syndrome (PCOS) (68). Increased HCG26 levels in patients with PCOS were related to antral follicle count. HCG26 knockdown in KGN cells suppressed cell proliferation and cell-cycle progression, and enhanced aromatase gene expression and estradiol production, suggesting that HCG26 may contribute towards the pathogenesis of PCOS (68). HCG26 has also been reported to be associated with nasopharyngeal carcinoma. A genome-wide study of copy number variation associated with nasopharyngeal carcinoma in a Malaysian-Chinese cohort identified candidate loci copy number variations at 11q14.3 and 6p21.3 (including a copy number variant region with HCG26), indicating that HCG26 may serve a role in the development of nasopharyngeal carcinoma (69). In the present study, HCG26 overexpression was associated with improved prognosis in MM, suggesting that it may act to protect against MM. The clinical roles and mechanisms of HCG26 in MM require further investigation.

The third previously reported lncRNA that exhibited potential prognostic value in MM was C21orf34. C21orf34 has been studied for its role in blood pressure by the Hypertension Genetic Epidemiology Network; African Americans and European Americans exhibited associations between blood pressure and intronic single nucleotide polymorphisms on chromosome 21q21.1 (70). The C21orf34 gene was linked to African American patients, improving understanding of the pathophysiology of hypertension (70). C21orf34 has also been studied in malignancy; C21orf34, which is the host gene of microRNA-125b, was reported to be downregulated in human metastatic melanoma (71). The role of C21orf34 in MM is yet to be determined. This study is the first to identify C21orf34 as an MM prognostic indicator. Its increased expression levels may predict the improved survival of patients with MM, suggesting that it may act as a protective factor against MM. As the prognostic value of C21orf34 was only determined via gene chip data mining, little is known regarding the functional role and mechanism of C21orf34 in MM. Therefore, further investigation is required.

To investigate the potential functional implications of prognostic markers for the onset and progression of MM, various bioinformatics computational methods were combined. The most reliable prognostic biomarkers identified in the present study were actively involved in cell cycle-associated processes. Sustained and proliferative signaling has been increasingly acknowledged as a fundamental trait of cancer cells, so the present findings are not unexpected (72). Previous studies reported that cell cycle interference may exert an antitumor function in MM (73,74). As a result of the complex mechanisms of tumorigenesis and tumor progression, a single gene is unlikely to underpin poor prognosis in MM. Therefore, a lncRNA-mRNA network was proposed to comprehensively explore the molecular characteristics of MM. The WGCNA results indicated that there may also be targeting relationships between the prognosis-associated lncRNAs and mRNAs. These lncRNAs may exert their prognostic effects by targeting closely associated mRNAs. Investigations of lncRNA-based regulatory networks are limited, particularly from the perspective of prognosis. Ronchetti et al (75) previously proposed a network constructed by lncRNAs and miRNAs. Further research should be conducted into the functional relationships between mRNAs and lncRNAs.

Certain shortcomings in the present study should be stated. A total of 559 cases were included in this study; however, the reported findings should be confirmed in additional independent cohorts. Furthermore, the prognostic values of the lncRNAs in this study were investigated using a gene chip; this single detection method should be verified by other methods, such as reverse transcription-quantitative PCR. Additionally, the majority of the lncRNAs identified in our prognostic model have not been previously reported. Their specific clinical significance, biological functions and potential mechanisms of action should be studied in further experiments. Finally, the molecular associations between identified lncRNAs and mRNAs in the expression network should be further investigated. Additional experiments are required to determine whether the prognosis-associated lncRNAs serve a role in MM via their corresponding mRNA targets.

In conclusion, the present study constructed a model that is capable of predicting prognosis in MM and generated a network with corresponding prognosis-associated mRNAs. Of note, the clinical significance and function of the majority of the lncRNAs identified in the present study remain unknown. These results offer novel perspective for the clinical diagnosis and treatment of MM and suggest novel directions for investigating the mechanisms underlying the development of MM.

Acknowledgements

Not applicable.

Funding

The present study was supported in part by the National Nature Science Foundation of China (grant no. 81560024), Program of Scientific and Technology Project, Guilin Science Research and Technology Development (grant no. 2016012706-2), National Natural Science Foundation of China (grant no. 81460038) and Guangxi Natural Science Foundation of China (grant no. 2017GXNSFAA198178).

Availability of data and materials

The datasets used and/or analyzed during the present study are available from the corresponding author on reasonable request.

Authors' contributions

YRL, ZPG and ZZY conceived and designed the study, as well as desiged the figures and tables. FXZ, XTW and ZZY contributed to the statistical analysis, as well as writing and correcting the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References