Identification of key genes and miRNA-mRNA regulatory networks associated with bone marrow immune microenvironment regulations in multiple myeloma by integrative bioinformatics analysis

ABSTRACT

The deregulation of microRNAs (miRNAs) and genes in the bone marrow microenvironment have been involved with the pathogenesis of multiple myeloma (MM). However, the exploration of miRNA-mRNA regulatory networks in MM remains lacking. We used GSE125363, GSE125361, GSE47552, GSE2658, GSE136324, GSE16558, and GSE13591 datasets for this bioinformatics study. We identified 156 downregulated and 13 upregulated differentially expressed miRNAs (DEmiRs) in MM. The DEmiRs are associated with the enrichment of pathways mainly involved with cancers, cellular signaling, and immune regulations. We identified 112 hub genes associated with five significant clusters in MM. Moreover, we identified 9 upregulated hub genes (such as IGF1, RPS28, UBA52, CDKN1A, and CDKN2A) and 52 downregulated hub genes (such as TP53, PCNA, BRCA1, CCNB1, and MSH2) in MM that is targeted by DEmiRs. The expression of DEmiRs targeted two hub genes (CDKN2A and TP53) are correlated with the survival prognosis of MM patients. Furthermore, the expression level of CDKN2A is correlated with immune signatures, including CD4+ Regulatory T cells, T cell exhaustion, MHC Class I, immune checkpoint genes, macrophages, neutrophils, and TH2 cells in the TME of MM. Finally, we revealed the consistently deregulated expression level of key gene CDKN2A and its co-regulatory DEmiRs, including hsa-mir-192, hsa-mir-10b, hsa-mir-492, and hsa-mir-24 in the independent cohorts of MM. Identifying key genes and miRNA-mRNA regulatory networks may provide new molecular insights into the tumor immune microenvironment in MM.

KEYWORDS:

• Multiple myeloma
• immune regulations
• immune signatures
• survival prognosis
• immune microenvironment

1. Background

Multiple myeloma (MM) is an incurable hematological malignancy of B lymphocytes characterized by the accumulation of malignant plasma cells in the bone marrow [1]. MM includes 10% of all hematological cancers and is comprised of anemia, hypercalcemia, renal failure, and bone lesions that are ultimately associated with the fatal conditions [2]. The stages of MM included monoclonal gammopathy of undetermined significance, smoldering myeloma, myeloma, and plasma cell leukemia and the end-stage disease is associated with resistance to treatments, increased cellular proliferation, evasion of apoptosis, and an ability to grow independently of the bone marrow microenvironment [2]. The bone marrow microenvironment has a crucial impact on MM. The interaction of plasma cells and bone marrow microenvironments can substantially regulate the disease's cellular proliferation, drug resistance, and prognosis [3]. In MM, the bone marrow niche comprises cellular components including stromal cells, osteoblasts, osteoclasts, endothelial cells, immune cells, and noncellular components including extracellular matrix, cytokines, chemokines, and growth factors [4]. The bone marrow tumor microenvironment influences the expression of MM-related genes associated with the angiogenesis, migration, and metastasis of cancer cells [5]. In the bone marrow microenvironment, targeting the plasma cells and their interaction with other components is a crucial strategy in MM [6].

MiRNAs are novel crossroads between the plasma cells and bone marrow microenvironment [7]. These regulatory miRNAs transfer between the components of the bone marrow microenvironment to establish their oncogenic functions [8]. The aberrant expression levels of miRNAs are found in the bone marrow microenvironment and blood of MM patients. These aberrantly expressed miRNAs are associated with the progression of MM [9]. M Fulciniti et al discovered that the miR-23b targeted Sp1 3'UTR and significantly reduced Sp1-driven nuclear factor-κB activity in MM cells [10]. MiRNAs and their predicted targets are associated with the development and malignancy of MM cells [11]. In response to IL-6, STAT3 is recruited to the miR-21 upstream region in human myeloma cells for contribution to the antiapoptotic function [12]. These studies indicate that the interactions of miRNA-mRNA are associated with the pathogenesis of MM.

Herein, we aimed to identify the differentially expressed genes (DEGs) and differentially expressed miRNAs (DEmiRs) in MM by using two gene expression datasets and a miRNA expression dataset. For these DEmiRs, we identified their associated pathways and target genes. For the hub genes in MM, we identified the significant modular functions in the MM. Furthermore, we investigated the correlations between the expression levels of the significant prognostic hub genes targeted by DEmiRs and antitumor immune signatures in the MM microenvironment. Finally, we identified the miRNA-mRNA interactions associated with immunosuppression in MM.

2. Materials and methods

2.1. Acquisition and preprocessing of microarray data

GSE125363, miRNA expression signatures in multiple myeloma, was downloaded from the NCBI. GSE125363 is based on GPL18044 (Agilent 046064 Unrestricted Human miRNA V19.0 Microarray) that included 44 myeloma samples (plasma cells) and 4 controls [13]. GSE125361, based on GPL20844, is a gene expression signatures database in multiple myeloma (n = 48) [13]. GSE47552 is also gene expression data based on GPL6244, including 41 multiple myeloma samples and 5 controls [14]. We included the plasma cells of multiple myeloma patients and excluded the other samples from the same dataset. Besides, we used the clinical information of the GSE2658 (n = 559) dataset for the survival analysis [15]. Furthermore, we utilized the dataset GSE136324 from 436 newly diagnosed MM patients related to the multiple myeloma (MM) whole tumor microenvironment [16]. Finally, we used the GSE16558 (n = 65, 60 MM patients and 5 control) [17] and GSE13591(n = 138, 133 MM patients and 5 control) [18] for validating the expression level of key genes and miRNAs.

2.2. Identification of differentially expressed miRNAs (DEmiRs) in multiple myeloma and target genes of DEmiRs in bone marrow

The differential expression of miRNAs was carried out using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), an interactive web tool. We utilized the GEO2R tool based on the GEOquery and limma R packages from the Bioconductor project (http://www.bioconductor.org/). The thresholds of adjusted P-value < 0.05 and |logFC| > 2.0 were set to find out the significant level. We used miRNet [19] (miRTarBase v8 and TarBase v8.0 were selected for the prediction of target genes) to identify the target genes of DEmiRs in the bone marrow. The online tool ‘Calculate and draw custom Venn diagrams’ (http://bioinformatics.psb.ugent.be/webtools/Venn/) was utilized for identifying common genes between the target genes of DEmiRs and the hub genes identified in multiple myeloma.

2.3. KEGG pathway enrichment analysis of DEmiRs in multiple myeloma

We identified significant KEGG [20] pathways that were associated with the DEmiRs by using the DIANA-mirPath (version 3) [21]. In this tool, we chose to analyze the validated miRNA-gene interactions using an archive in microT web server v5.0 [22] and DIANA-TarBase v7.0 [23] and identified KEGG pathways. We uploaded 100 downregulated and all upregulated miRNAs into the DIANA-mirPath (version 3). A threshold of adjusted P-value < 0.05 (Fisher's exact test) was used to define the statistical significance. Finally, we identified the common significant pathways enriched between both tools (microT web server v5.0 and DIANA-TarBase v7.0).

2.4. Identification of DEGs in multiple myeloma

We used Network Analyst [24] to identify the DEGs between multiple myeloma and normal samples by a meta-analysis of two multiple myeloma gene expression profiling datasets (GSE125361 and GSE47552). Datasets were normalized by quantile normalization, and the batch effects of multiple datasets were removed using the ComBat method [25]. We used the R package ‘limma’ to identify the DEGs between multiple myeloma and normal samples, and Cochran's combination test performed the meta-analysis [26]. The false discovery rate (FDR) [27] was used to adjust for multiple tests. We identified the DEGs using a threshold of absolute combined effect size (ES) >1.50 and FDR < 0.05.

2.5. Construction of protein–protein interaction (PPI) network of the DEGs and identification of hub genes-associated clusters

We constructed a PPI network of the identified DEGs by search tool for the retrieval of interacting genes (STRING) database (https://string-db.org/cgi/input.pl; STRING-DB v11.0)[24]. The rank of the target proteins based on the degree of interactions in the PPI network was identified using the Cytoscape plugin cytoHubba[25]. The obtained protein–protein interaction data were imported into Cytoscape 3.6.1 software to construct a PPI network[26]. A hub gene was defined as a gene connected to no less than 20 other DEGs in the PPI. We used Cytoscape plugin molecular complex detection (MCODE) to identify the associated modules with hub genes. We identified the significant modules based on the MCODE score and node number. The threshold of the MCODE was Node Score Cutoff: 0.2, Haircut: true, K-Core: 2, and maximum depth from Seed: 100.

2.6. Gene-set enrichment analysis of significant modules

We performed gene-set enrichment analysis of the cluster-specific hub genes by GSEA [27]. We inputted the hub genes associated with a specific cluster into the GSEA. The KEGG pathways significantly associated with the hub genes in an individual cluster were identified (FDR < 0.05).

2.7. Survival analysis of hub genes in multiple myeloma by using the PrognoScan

The prognostic roles of DEmiRs targeted hub genes in multiple myeloma patients were analyzed using the PrognoScan database (http://www.abren.net/ PrognoScan/) [28]. We used the dataset GSE2658 (n = 559) for the survival analysis [15]. Cox P < 0.05 was considered a significant level.

2.8. Evaluation of the immune signature enrichment levels with prognostic hub genes

We used dataset GSE136324 from 436 MM patients for analyzing the correlation of immune infiltrations with prognostic hub genes in the multiple myeloma (MM) tumor microenvironment [16]. We calculated the immune scores, stromal scores, and tumor microenvironment scores using the xCell tool [29]. Also, we identified the enrichment level of various immune signatures in a tumor sample as the single-sample gene-set enrichment analysis (ssGSEA) score [30]. The gene set contains the collection of all marker genes of an individual immune signature. We included 20 immune signatures: B cell, CD8+ T cells, CD4 + regulatory T cells, NK cells, macrophages, neutrophils, pDCs, MHC Class I, HLA genes, Tregs, T cell exhaustion, immune checkpoint genes, T cell activation, cytolytic activity, MDSCs, TH1, TH2, monocyte, CAFs, and Tfh. The marker genes set of individual immune signatures are displayed in Supplementary Table S1. Then, we identified Spearman's correlation between the ssGSEA score and the expression levels of selected hub genes (absolute value of R is not less than 0.20, P < 0.05).

2.9. Statistical and computational analysis

We evaluated Pearson's or Spearman's correlation test to verify the significant levels between the two variables. For analyzing the correlations between the expression levels of hub genes and the enrichment levels (ssGSEA scores) of immune signatures, we used Spearman's correlation test because these data were not normally distributed. For analyzing the correlations between the expression levels of hub genes with the expression levels of other marker genes, we utilized Pearson's correlation test because these data were normally distributed. We employed the Wilcoxon rank-sum test for the expression validation of key genes and miRNAs between the tumor and control groups (P < 0.05). We utilized the NetworkAnalyst [31] tool to calculate the average expression of a gene with multiple probes in the same expression dataset. We used the R package ‘ggplot2’ for plotting the graphs in this study.

3. Results

3.1. Identification of DEmiRs in multiple myeloma

We identified the upregulated and downregulated DEmiRs between the multiple myeloma and control samples. We identified 13 upregulated DEmiRs such as hsa-miR-4301, hsa-miR-548al, hsa-miR-4693-3p, hsa-miR-711, and hsa-miR-4799-3p (Supplementary Table S2). In addition, 156 downregulated DEmiRs were identified in the plasma cells of MM patients (Supplementary Table S2). The top ten (the absolute value of the highest LogFC) upregulated and downregulated DEmiRs are presented in Table 1. Previous studies showed that the DEmiRs obtained from the MM are associated with cancers. For example, hsa-miR-4301 is upregulated in gastric cancer treated with cisplatin [32]. Hsa-miR-299-5p, another upregulated DEmiRs, is associated with the epithelial–mesenchymal transition, proliferation, and tumorigenicity of colorectal cancer cells [33]. Downregulation of miR-23a-5p levels in leukemic cells is associated with their protection of them from chemotherapy-induced apoptosis [34]. In tumors, miR-378a is considered a part of an angiogenic network [35]. Altogether, DEmiRs are associated with the pathogenesis of cancers, including hematological disorders.

Table 1. Top 10 upregulated and top 10 downregulated DEmiRs in multiple myeloma.

ID	Adjusted P value	P value	LogFC	miRNA name
6193	1.52E-02	1.05E-04	2.536	hsa-miR-4301
4589	2.17E-02	1.89E-04	2.452	hsa-miR-548al
35638	3.73E-02	4.62E-04	2.443	hsa-miR-4693-3p
15327	4.67E-03	1.38E-05	2.424	hsa-miR-711
47809	4.51E-02	6.14E-04	2.403	hsa-miR-4799-3p
57759	3.70E-02	4.49E-04	2.331	hsa-miR-5583-5p
6709	1.52E-02	1.05E-04	2.222	hsa-miR-299-5p
50274	2.78E-02	2.82E-04	2.204	hsa-miR-6130
38405	1.83E-02	1.44E-04	2.196	hsa-miR-4481
50878	2.57E-02	2.46E-04	2.188	hsa-miR-758-3p
6629	2.05E-02	1.72E-04	−2.589	hsa-miR-4671-5p
3372	1.42E-02	9.12E-05	−2.607	hsa-miR-4681
15233	8.46E-03	3.76E-05	−2.638	hsa-miR-4755-5p
58185	1.16E-04	3.77E-08	−2.668	hsa-miR-1273a
60712	2.71E-02	2.74E-04	−2.699	hsa-miR-15b-3p
25588	1.16E-04	3.91E-08	−2.706	hsa-miR-3177-3p
25883	4.82E-02	6.78E-04	−2.762	hsa-miR-2355-3p
59150	9.95E-04	1.21E-06	−3.046	hsa-miR-4744
49869	3.20E-02	3.60E-04	−3.117	hsa-miR-378a-3p
49708	1.54E-02	1.09E-04	−3.296	hsa-miR-23a-5p

3.2. Identification of pathways significantly associated with DEmiRs

We used microT web server v5.0 [22] and DIANA-TarBase v7.0 [23], which are built-in DIANA-miRPath [21] to identify the KEGG [20] pathways. We identified common significant 48 KEGG [20] pathways associated with the downregulated DEmiRs in MM (Table 2). These pathways are mainly involved in cancers, immune-regulation, stromal regulations, metabolisms, cellular signaling, and cellular developments (Table 2). Previous studies revealed that the enriched pathways are associated with hematological cancers. For example, it was reported that proteoglycans are required for multiple myeloma cell adhesion, in vivo growth, and vascularization [36]. The hippo pathway is associated with resistance in MM and other hematologic malignancies [37]. In MM, the TGF-beta signaling pathway can control proliferation, differentiation, migration, and cell survival [38]. In addition, we also identified 27 KEGG pathways associated with upregulated DEmiRs in MM (Supplementary Table S3).

Table 2. Downregulated DemiRs are associated with the enrichment of KEGG pathways. MicroT web server v5.0 [22], and DIANA-TarBase v7.0 [23] are used for identifying common pathways.

KEGG pathways	P-value (MicroT-CDS v5.0)	P-value (TarBase v7.0)
Proteoglycans in cancer	3.22E-09	3.73E-14
Pathways in cancer	5.33E-09	3.69E-04
Hippo signaling pathway	7.22E-09	6.80E-07
Glioma	2.48E-07	8.62E-06
Phosphatidylinositol signaling system	2.48E-07	2.29E-02
TGF-beta signaling pathway	7.16E-07	1.92E-04
ErbB signaling pathway	1.21E-06	3.69E-04
Signaling pathways regulating pluripotency of stem cells	1.31E-06	4.57E-03
Fatty acid metabolism	3.02E-06	1.22E-11
Thyroid hormone signaling pathway	4.32E-05	1.92E-04
Wnt signaling pathway	5.35E-05	1.39E-02
Renal cell carcinoma	6.02E-05	8.33E-07
Fatty acid biosynthesis	6.74E-05	1.69E-03
Chronic myeloid leukemia	9.04E-05	8.62E-06
Focal adhesion	9.04E-05	3.62E-04
Endocytosis	1.99E-04	3.62E-04
Melanoma	2.48E-04	1.76E-02
Prostate cancer	2.50E-04	3.13E-05
Adherens junction	4.38E-04	1.30E-06
Colorectal cancer	6.91E-04	9.67E-04
Pancreatic cancer	9.28E-04	1.47E-03
Endometrial cancer	1.04E-03	5.51E-04
Neurotrophin signaling pathway	1.06E-03	5.67E-04
Acute myeloid leukemia	1.22E-03	6.69E-03
PI3K-Akt signaling pathway	1.22E-03	2.59E-02
Non-small cell lung cancer	1.38E-03	1.09E-04
FoxO signaling pathway	2.93E-03	1.26E-04
Oocyte meiosis	2.93E-03	3.32E-06
Ubiquitin mediated proteolysis	3.30E-03	3.32E-06
mTOR signaling pathway	3.79E-03	1.41E-04
Fatty acid degradation	5.89E-03	2.50E-05
p53 signaling pathway	6.28E-03	6.93E-04
Viral carcinogenesis	6.83E-03	1.06E-07
Prion diseases	7.62E-03	1.30E-06
Central carbon metabolism in cancer	8.51E-03	5.73E-04
Hepatitis B	1.02E-02	1.92E-04
Prolactin signaling pathway	1.14E-02	8.72E-03
Transcriptional misregulation in cancer	1.14E-02	4.52E-04
Regulation of actin cytoskeleton	1.17E-02	3.44E-02
MAPK signaling pathway	1.66E-02	3.35E-02
Bacterial invasion of epithelial cells	1.82E-02	2.50E-05
AMPK signaling pathway	2.05E-02	5.51E-04
Lysine degradation	2.12E-02	1.70E-06
Thyroid cancer	2.77E-02	1.44E-03
TNF signaling pathway	2.87E-02	7.92E-04
HIF-1 signaling pathway	3.06E-02	2.59E-02
Estrogen signaling pathway	4.05E-02	2.48E-07
N-Glycan biosynthesis	4.17E-02	1.25E-02

3.3. Identification of differentially expressed hub genes and their modular analysis

We identified 208 upregulated and 550 downregulated genes (Supplementary Table S4) by a meta-analysis of two GEO datasets (GSE125361and GSE47552) (absolute value of combined ES > 1.5, P < 0.05). The top ten upregulated genes (the highest combined ES) are ARPC5L, CDKN1A, IGF1, UBALD2, RPS28, RABAC1, SURF1, PELI1, TBCC, and YBEY in MM. It was stated that the expression of ARPC5 was significantly higher in relapsed MM cells than in baseline MM cells [39]. IGF1, one of the top upregulated genes, is a major growth factor for MM cell lines and is associated with prognostic relevance in MM [40]. In addition, NCAPD2, BORA, MSH2, TF, CLUAP1, CTPS2, ZNF382, ZMYM1, DDIAS, and SF3B3 are the top ten (the lowest combined ES) downregulated genes in MM. In MM, MSH2 has interacted with RECQ1 helicase, which is ultimately associated with DNA repair pathways [41]. A study has reported that a transcription inhibitor ZNF382 is downregulated in MM due to hypermethylation in the promoter region[42].

We inputted all DEGs into the STRING v11 to construct protein–protein interaction networks for identifying hub genes in MM. We identified 112 hub genes (degree of interaction ≥ 20) in MM that included 10 upregulated DEGs, and 102 downregulated DEGs (Figure 2A and Supplementary Table S5). The upregulated hub genes in MM included IGF1, IRF1, RPS28, UBA52, RPL32, RPS19, CDKN1A, RNF7, PSMB8, and CDKN2A. Some of the downregulated hub genes in MM are TP53, BRCA1, CCNB1, CDC20, PCNA, CDC6, BPTF, POLR2B, HDAC2, SMC1A, and MSH2. Furthermore, we identified the significant clusters of identified hub genes by using the Cytoscape plugin MCODE. The hub gene-associated five significant clusters in MM (Table 3, Figure 1(A, C, E, G, and H)). The functional enrichment analysis revealed that cluster 1 (Figure 1A) is associated with the enrichment of 20 KEGG pathways, including cell cycle, pathways in cancer, p53 signaling pathway, mismatch repair, DNA replication, bladder cancer, glioma, melanoma, chronic myeloid leukemia, small cell lung cancer, prostate cancer, homologous recombination, nucleotide excision repair, non-small cell lung cancer, pancreatic cancer, and TGF-beta signaling pathway (Figure 1B). Hub gene-associated cluster 2 (Figure 1C) is associated with the enrichment of 3 KEGG pathways (ubiquitin-mediated proteolysis, oocyte meiosis, and proteasome (Figure 1D)), and cluster 3 (Figure 1E) is related to the enrichment of 3 other KEGG pathways (cell cycle, nucleotide excision repair, and pancreatic cancer (Figure 1F)). It indicated that the hub genes-associated clusters are associated with the enrichment of cancer-related pathways in MM.

Figure 1. Construction of cluster-specific PPI of hub genes and cluster-specific pathway enrichment analysis. A. Hub genes of cluster 1 involved with PPI. B. KEGG pathways that are significantly associated with hub genes of cluster 1. C. Hub genes of cluster 2 involved with PPI. D. Three KEGG pathways are enriched in Cluster 2. E. PPI of hub genes in cluster 3. F. Three KEGG pathways are enriched associated with the hub genes of cluster 3. G-H. PPI of hub genes in cluster 4 and cluster 5.

Table 3. MCODE identified hub gene-associated five significant clusters in MM.

Cluster	Score	Nodes	Edges	Hub genes
1	13.2	31	198	UBA52, TOP1, CENPE, ANAPC7, RBL1, CDKN1A, TFDP1, RAD50, SMC6, CDKN2A, ORC3, TOP2A, BORA, SKP2, RRM1, MCM8, PRIM1, ANAPC4, RBBP4, RFC5, NCAPD2, WRN, MSH2, CDC16, CCNB1, RB1, CCNE1, CDC23, PCNA, MDC1, RAD51C
2	10.32	26	129	BTRC, SPDL1, KIF20A, CENPF, SLBP, FBXW8, NUP37, SMC1A, WDHD1, HUWE1, CDC20, NUP107, FBXO5, UBE2O, NCAPG, XPO1, RNF7, ESCO2, FBXL20, PSMD1, FBXO7, NUP160, PSMB8, DSN1, MKI67, CKAP5
3	5.091	12	28	BPTF, BRCA2, GTF2H1, CDC6, POLR2B, RFC1, ATF2, ORC4, ORC2, BRCA1, TP53, TOP2B
4	4.667	7	14	HDAC2, CTCF, TBP, SMARCC2, KDM1A, SUPT16H, SMARCC1
5	2.727	12	15	CXCL12, NOC3L, CAT, CEP290, WDR36, NAT10, PRKAR2B, BBS10, IGF1, NEDD1, STAT5B, GART

3.4. DEmiRs target the differentially expressed hub genes in MM

We predicted the target genes of DEmiRs in bone marrow by using the miRNet [19] (miRTarBase v8 and TarBase v8.0 were selected for predicting the genes). We predicted 10,995 genes for downregulated DEmiRs (Figure 2B) and 3781 genes for upregulated DEmiRs (Figure 2C) using miRNet, respectively. We found that the upregulated IGF1, IRF1, RPS28, UBA52, RPS19, CDKN1A, RNF7, PSMB8, and CDKN2A hub genes interact with the downregulated DEmiRs and the downregulated 52 hub genes (such as TP53, PCNA, BRCA1, CCNB1, POLR2B, TOP2A, MSH2, CDC6, BPTF, and CDC20) are interacted with the upregulated DEmiRs in bone marrow (Supplementary Table S6). Haimeng Yan et al. revealed that UBA52 and RPS19 are differentially expressed in MM and associated with the enrichment of signaling pathways [43]. It was reported that TP53-inducible miRNAs are related to the development of MM [44].

Figure 2. Differentially expressed hub genes were targeted by DEmiRs. A. Identification of upregulated and downregulated hub genes. B. The upregulated 9 hub genes interact with downregulated DEmiRs. C. The downregulated 52 hub genes interact with upregulated DEmiRs.

3.5. DEmiRs targeting hubs genes are associated with survival in MM

We investigated the survival of DEmiRs targeted all hub genes (9 upregulated and 52 downregulated). We used the Prognoscan [28] for survival analysis and found that two hub genes (CDKN2A and TP53) had a negative expression correlation with survival prognosis (Disease-specific survival) in MM (GSE2658, n = 559). We found that the elevated level of upregulated hub gene CDKN2A is associated with poor prognosis, and the lower level of downregulated hub gene TP53 is correlated with a poor survival prognosis in MM patients (Figure 3).

Figure 3. DEmiRs targeting hubs genes are associated with survival prognosis in multiple myeloma patients. Kaplan-Meier curves show that the expression of two hub genes (CDKN2A and TP53) correlates with survival prognosis (Disease-specific survival) in MM patients (log-rank test, P < 0.05).

3.6. CDKN2A is associated with the regulation of tumor microenvironment in MM

The level of tumor-infiltrating lymphocytes (TILs) is an independent predictor of survival in cancer [45]. We investigated the correlations between the expression levels of two hub genes (CDKN2A and TP53) and the levels of immune signatures in the TME of MM. Interestingly, we found that the expression level of CDKN2A is negatively correlated with immune scores (R = −0.27, P < 2.2e-16) and microenvironment scores (R = −0.31, P < 2.2e-16) (Figure 4A), but TP53 is not correlated with immune scores and tumor microenvironment scores. Furthermore, the expression level of CDKN2A is positively correlated with CD4+ Regulatory T cells, T cell exhaustion, MHC Class I, and immune checkpoint genes (Figure 4B) and negatively correlated with the macrophages, Neutrophils, and TH2 cells in the TME of MM (Figure 4B). It was indicated that the CD4 T cells are responsible for controlling autoimmune phenomena in MM [46]. T cell exhaustion in multiple myeloma is associated with the relapse of the disease [47]. Immune checkpoint pathways are associated with impeding antitumor immune responses and contributing to malignancy's persistence and/or relapse [48]. Macrophages and neutrophils are now considered critical regulatory components of the hematopoietic niche [49]. Altogether, these results indicate that the expression level of CDKN2A is associated with the immunosuppressive TME in MM.

Figure 4. The expression level of CDKN2A is associated with the regulation of the tumor microenvironment in MM. A. Immune scores and microenvironment scores are negatively correlated with the expression levels of CDKN2A. B. Immune signatures are significantly correlated with the expression levels of CDKN2A.

3.7. CDKN2A is positively correlated with immunosuppressive markers in MM

Since CD4+ Regulatory T cells, T cell exhaustion, and immune checkpoint genes are positively correlated with the expression levels of CDKN2A (Figure B), we calculated the correlation of marker genes with the expression of CDKN2A. Interestingly, we found that the immunosuppressive markers including PDCD1 (PD-1), CTLA4, FOXP3, HHLA2, CD70, TMIGD2, TNFSF9, and IL5 are positively correlated with the expression levels of CDKN2A in the TME of MM (Figure 5). The Akt pathway is activated in MM through PD-1-bound PD-L1, which ultimately plays a substantial role in the progression of MM disease [50]. Walter et al. discovered that the expression levels of FOXP3 and CTLA4 presented a sixfold and 30-fold higher, respectively, in MM patients than in controls [51]. Another immunosuppressive marker, HHLA2, significantly reduces cytokine production by T cells that ultimately inhibit the functions of CD8 T-cells [52]. It suggests that the expression level of CDKN2A is associated with the immunosuppressive TME in MM.

Figure 5. Positive correlation of CDKN2A with immunosuppressive markers in MM microenvironment. Immunosuppressive markers PDCD1 (PD-1), CTLA4, FOXP3, HHLA2, CD70, TMIGD2, TNFSF9, and IL5 are positively correlated with the expression levels of CDKN2A in the microenvironment of MM.

3.8. The expression validation of CDKN2A and its co-regulatory DEmiRs in the independent MM cohorts

Since CDKN2A is associated with the survival prognosis and immune infiltrations in MM, we validated the expression level of CDKN2A and its co-regulatory miRNAs (Figure 6). We utilized the GSE16558 [17] and GSE13591[18] for the expression level validation of CDKN2A and its co-regulatory miRNAs, including hsa-mir-492, hsa-mir-24, hsa-mir-192, and hsa-mir-10b (Figure 7). Interestingly, we found that these four co-regulatory miRNAs (hsa-mir-492, hsa-mir-24, hsa-mir-192, and hsa-mir-10b) are consistently downregulated in the GPL8695 of GSE16558 (Wilcoxon rank-sum test, P < 0.05) (Figure 6A). Moreover, we revealed that the expression level of CDKN2A is consistently upregulated in the GSE16558 and GSE13591 (Wilcoxon rank-sum test, P < 0.05) (Figure 6B). It indicates that the expression level CDKN2A and its co-regulatory miRNAs are consistently deregulated in MM patients.

Figure 6. The expression validation of CDKN2A and its co-regulatory DEmiRs in the independent MM cohorts. A. The expression level of CDKN2A co-regulatory miRNAs: hsa-mir-492, hsa-mir-24, hsa-mir-192, and hsa-mir-10b are consistently downregulated in the GSE16558. B. The expression level of CDKN2A is consistently upregulated in the GSE16558 [17] and GSE13591[18].

Figure 7. MiRNA-CDKN2A interactions in MM. The downregulation of miRNAs interacts with the upregulation of CDKN2A. This co-regulatory network is associated with the immunosuppressive TME of MM.

4. Discussion

We identified DEGs in MM by a meta-analysis of two gene expression datasets. In addition, DEmiRs were identified in MM relative to normal cells. We found that several enriched pathways are associated with DEmiRs (Table 2). The pathways mainly involved cancers, cellular signaling, and immune regulations. PPI-based analysis revealed that the hub genes are involved in significant clusters associated with KEGG pathway enrichments (Figure 1). Based on the differentially expressed genes and miRNA, we build mRNA-miRNA regulatory networks in MM. Two DEmiRs targeting hub genes (CDKN2A and TP53) are linked with the survival prognosis of MM patients (Figure 3). The CDKN2A is upregulated (LogFC = 2.0529 and adjusted P value = 7.82E-05) in MM patients, and the survival probability of the upregulated-CDKN2A group is lower in MM patients. In contrast, the lower expression group of downregulated-TP53 group (LogFC = −1.7835and adjusted P value = 0.0005) is correlated with a shorter survival time in MM patients. It indicates that the higher expression group of CDKN2A and the lower expression group of TP53 are associated with poor survival prognosis in MM. Since the level of immune infiltration is an independent predictor of survival in cancer [45], we investigated the correlation of CDKN2A and TP53 with the immune infiltration levels in the bone marrow microenvironment of MM. First, the expression level of CDKN2A is negatively correlated with immune scores and microenvironment scores (Figure 4A), indicating that the higher expression level of CDKN2A is associated with lowering the immunity in MM patients. Second, the expression level of CDKN2A is positively correlated with various cells, including CD4+ Regulatory T cells, T cell exhaustion, MHC Class I, and immune checkpoint genes, and negatively correlated with macrophages, neutrophils, and TH2 cells in the TME of MM (Figure 4B). Third, the immunosuppressive markers including PDCD1, CTLA4, FOXP3, and IL5 are positively correlated with the expression levels of CDKN2A (Figure 5). It suggests that the higher expression level of CDKN2A is associated with the suppression of TME by depressing the immune-modulatory cellular components and activating the immunosuppressive markers and immune-inhibitory cellular components. Altogether, CDKN2A mediated the MM pathogenesis by regulating the bone immune microenvironment.

Based on the hub genes and their regulatory DEmiRs, we constructed miRNA-mRNA regulatory networks in MM that may be associated with the regulation of the bone microenvironment. These networks included the hsa-mir-192-5p-CDKN2A, hsa-mir-10b-5p-CDKN2A, hsa-mir-492-CDKN2A, and hsa-mir-24-1-5p-CDKN2A (Figure 7). In newly diagnosed MM patients, the overexpression of CDKN2A is associated with a shorter survival time [53]. In the human tumor, a lower expression level of miR-192 could predict poor clinical outcomes in several cancer types [54]. The deregulated expression level of hsa-mir-10b is a fundamental miRNA in various metastatic cancer [55]. The aberrant expression of miR-24 contributes to tumorigenesis, tumor progression, and poor clinical outcomes in multiple cancers [56]. Moreover, our validation results revealed that the expression of the level of CDKN2A, hsa-mir-492, hsa-mir-24, hsa-mir-192, and hsa-mir-10b is consistently deregulated in the MM cohorts (Figure 6). This result suggests that the hsa-mir-192-5p-CDKN2A, hsa-mir-10b-5p-CDKN2A, hsa-mir-492-CDKN2A, and hsa-mir-24-1-5p-CDKN2A networks may be contribute to the immunosuppression in MM (Figure 7).

Some previous studies tried to identify the microRNA-mRNA regulatory networks in multiple myeloma. For example, Hongyu Gao et al. identified the miR-135b-GADD45A and miR-148a-USPL1 pairs in MM [57]. Yang Yang et al. identified DEmRs and predicted their targets [58]. Also, Marta Lionetti et al. identified microRNA expression patterns and predicted microRNA/mRNA regulatory networks in distinct molecular groups of multiple myeloma[59]. However, these authors did not discover how these co-regulatory networks are associated with the regulation of TME in MM. For the first time, we identified the hsa-mir-192-5p-CDKN2A, hsa-mir-10b-5p-CDKN2A, hsa-mir-492-CDKN2A, and hsa-mir-24-1-5p-CDKN2A pairs in MM, which ultimately associated with immunosuppressive TME (Figure 7). The major drawback of this study is that the MM-associated miRNA–mRNA regulatory network identified by bioinformatics analysis has not been validated by experimental research. Thus, although our findings could provide potential biomarkers for multiple myeloma diagnosis and therapeutic targets, further experimental and clinical validation is necessary to transform these results into practical application in MM treatments.

5. Conclusions

Identifying key genes and miRNA-mRNA regulatory networks in MM may provide new molecular insights into the miRNA-mRNA regulatory networks mediated mechanism of the bone marrow tumor immune microenvironment in the MM.

Authors’ contributions

All authors have contributed to the research concept, design, and interpretation of the data. All authors reviewed and approved the final version of the document.

Data availability

The GEO datasets GSE125363, GSE125361, GSE47552, GSE2658, GSE136324, GSE16558, and GSE13591 used in this study are available in The National Biotechnology Information Center Gene (NCBI-Gene) database (https://www.ncbi.nlm.nih.gov/gene).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by National Natural Science Foundation of China: [Grant Number 81560027].