Genome-wide identification of potential biomarkers in multiple myeloma using meta-analysis of mRNA and miRNA expression data

Multiple myeloma (MM) is a plasma cell malignancy with diverse clinical phenotypes and molecular heterogeneity not completely understood. Differentially expressed genes (DEGs) and miRNAs (DEMs) in MM may influence disease pathogenesis, clinical presentation / drug sensitivities. But these signatures overlap meagrely plausibly due to complexity of myeloma genome, diversity in primary cells studied, molecular technologies/ analytical tools utilized. This warrants further investigations since DEGs/DEMs can impact clinical outcomes and guide personalized therapy. We have conducted genome-wide meta-analysis of DEGs/DEMs in MM versus Normal Plasma Cells (NPCs) and derived unified putative signatures for MM. 100 DEMs and 1,362 DEGs were found deranged between MM and NPCs. Signatures of 37 DEMs (‘Union 37’) and 154 DEGs (‘Union 154’) were deduced that shared 17 DEMs and 22 DEGs with published prognostic signatures, respectively. Two miRs (miR-16–2-3p, 30d-2-3p) correlated with survival outcomes. PPI analysis identified 5 topmost functionally connected hub genes (UBC, ITGA4, HSP90AB1, VCAM1, VCP). Transcription factor regulatory networks were determined for five seed DEGs with ≥ 4 biomarker applications (CDKN1A, CDKN2A, MMP9, IGF1, MKI67) and three topmost up/ down regulated DEMs (miR-23b, 195, let7b/ miR-20a, 155, 92a). Further studies are warranted to establish and translate prognostic potential of these signatures for MM.

Expression profiles of differentially expressed genes (DEGs) are of paramount importance and have provided critical prognostic insights in MM. Recent transcriptome based studies have reported gene expression prognostic (GEP) signatures associated with tumor classification, survival risk prediction 7,8 , progression of MM 7,[9][10][11][12] , response to drugs 13 , chromosome instbility 14 and others. The DEGs included in GEP signatures are diverse but closely connected to similar pathways. These genes may relate to kinome 15 , autophagy 16 , cell cycle 10,17 , stemness 18 , cytogenetic abnormalities 9,19 , chromosome 1 20 , homozygous deletions, cell death 21 and immune 7 subnetworks. At least 8 to 10 molecular subgroups of MM based on the genomic and transcriptomic patterns have been reported that tend to correlate with different clinical outcomes 8 . Computational and functional analysis of hub genes, nodes, networks and pathways in MM have led to the development of risk scoring systems, relating to the seven genetic subgroups 22 , 70 genes UAMS70 risk signatures 20 , IFM15 risk stratification 17 , 5 gene stemness score 18 , UAMS 17 23 , CINGLEC 214 14 , HOVON-65/GMMG-HD4 EMC 92 9 , HZD 97 21 , M3CN 10 and others.
However, the prognostic scores derived from GEP signatures have low prediction accuracy and limited power to predict risk or response 23 perhaps due to MM heterogeneity and complex interactions between malignant plasma cells and bone marrow environment. A landmark study reported GEP prognostication could be improved when a combination of EMC92 + HZDCD 24 was used. A similar integrative M3CN network study 10 on MMRF-CoMMpass cohort unified eight prognostic gene signatures and demonstrated significantly improved prognostic performance.
Even though a series of MM associated potential DEM/DEG signatures have been identified across several studies over the years, these remain mostly heterogeneous and challenging to interpret in clinics. There are still unresolved questions such as their mutual interdependencies, interactions with microenvironment and their combinatorial synergistic prognostic and therapeutic potentials. There are still lacunae in our knowledge and more studies are needed to understand the signatures that are best valued in clinics for fast and early prognostication of newly diagnosed MM patients. It is thus postulated that a comprehensive analysis of individual GEP identifiers in MM PCs as compared to normal PCs across multiple studies will help unfold common signatures with potential prognostic significance. In this regard, we have performed a meta-analysis of available multiple datasets of DEGs and DEMs in MM patients to derive a unified set of core GEP signatures. We have identified a combination of 'Union 154' DEGs and 'Union 37' DEMs that may aid in achieving improved prognosis and clinical applicability.
Genome-wide miRNA and mRNA expression profiling. Total RNA was isolated from CD138 + plasma cells enriched with MACS beads (Miltenyi Biotech, Germany), collected from 44 newly diagnosed treatment naïve MM patients diagnosed as per IMWG guidelines 45 (Table S2 in Supplementary File 1) and 4 controls (pooled from 10 Hodgkin's disease bone marrow samples). Total RNA was extracted using the miRVana miRNA isolation kit (Thermofisher Scientific, MA, USA).
For the genome wide miRNA expression profiling, RNA was labeled and hybridized to an unrestricted human microRNA v19 Microarray slide (Agilent 046,064, GPL18044) (Agilent Technologies, Santa Clara, CA, USA) according to the manufacturer's protocol. Briefly, 100 ng of total RNA was labeled with Cyanine3 (Cy3) using miRNA Complete Labeling and Hybridization Kit (Agilent Technologies, Santa Clara, CA, USA). The Cy3labeled samples were resuspended in hybridization buffer and hybridized onto Human miRNA 8X60K format microarrays (Agilent Technologies, Santa Clara, CA, USA) at 55 °C for 20 h. After hybridization, microarrays were washed with gene expression wash buffer and the fluorescent signals were scanned using SureScan microarray scanner D (Agilent Technologies, Santa Clara, CA, USA) using one colour scan settings (Scan resolution 3 μm, Dye channel set to Green, Green PMT = 100%). The data generated on miRNA expression in MM using microarrays has been submitted to GEO database with accession no GSE125363.
To correlate whether the miRNA alteration affects gene expression, mRNA expression array analysis was also performed on 44 MM patient samples and 4 controls (pooled from 10 Hodgkin's disease bone marrow samples). Double-stranded cDNA was generated from 200 ng total RNA (isolated with miRVana kit) using the low input quick amp labelling kit (Agilent Technologies, Agilent Technologies, Santa Clara, CA, USA) using T7 primer, dNTPs and affinity script RNase block. Next, cDNA was transcribed to cRNA using T7 RNA polymerase and NTP mix and labeled with Cyanine3 using Cy3-CTP. The labeled cRNA was purified according to manufacturer's protocol using RNAeasy extraction kit (Qiagen, Hilden, Germany). The concentration of Cyanine3 and cRNA was measured using NanoDrop ND1000 spectrophotometer. Samples with specific activity ≥ 6 pmol Cy3/µg cRNA were hybridized onto a SurePrint G3 human GE v3 8 × 60 K microarray slide (Agilent 072,363, GPL20844) (Agilent technologies, Santa Clara, CA, USA), and incubated for 17 h at 65 °C in a hybridization oven. The slides were washed and scanned in SureScan microarray scanner D (Agilent technologies, Santa Clara, CA, USA) with scan settings (Scan resolution 3 μm, Dye channel set to Green, Green PMT = 100%). The data generated on mRNA expression in MM has been submitted to GEO database with accession no. GSE125361.
Preprocessing and mining of DEMs/DEGs from Agilent array. Microarray images (*.tiff) obtained from SureScan scanner were quantified using Agilent Feature Extraction Software (version 11.5.1.1) (Agilent Technologies, Santa Clara, CA, USA). The raw data obtained (tab-delimited text file per hybridisation) was subsequently processed with the Limma R package available in the Bioconductor repository (http:// www. bioco nduct or. org). Limma used linear model statistics to find genes that were differentially expressed between the patients and controls. The raw intensity data were background corrected using normexp method and subsequently normalized using quantile method for one-color. Expression level variations between replicates were analyzed by pairwise comparisons using the lmFit function. The fitted model object was further processed by the eBayes function to produce empirical Bayes test statistics for each gene, including moderated t-statistics, www.nature.com/scientificreports/ p-values and log-odds of differential expression. The t test and Benjamini and Hochberg method were used to calculate the p-values and false discovery rate (FDR) 44 . The adjusted p ≤ 0.05 and |logFC|≥ 1.5 were set as the cut-off criterion for identifying DEGs and DEMs.

Meta-analysis of DEGs/DEMs datasets.
A widely used meta-analysis approach [46][47][48] was applied to integrate the gene/miRNA expression profiles obtained independently from GEO repository as well as datasets generated at our centre following microarray hybridization (Table S1 in Supplementary File 1). The gene and miRNA probes were assigned as per HGNC (HUGO Gene Nomenclature Committee) and miRBase-22.1 identifiers, respectively using the g:Profiler 49 (https:// biit. cs. ut. ee/ gprofi ler/). The differentially expressed genes/miR-NAs obtained through R/Bioconductor limma package 43 from these individual studies were merged by taking the union across them. When multiple probes referred to the same gene/miRNA, the expression values obtained from these probes were minimized to a single value by averaging the expression value (when in the same direction of expression) or were discarded (when had diverse directions of expression). The probes with unknown gene or unknown miRNA identifiers or annotated as antisense RNA, chromosomes, hypothetical loci, non-coding RNAs, non-functional proteins, non-protein coding genes, pseudo-genes and uncharacterized genes were discarded. The DEGs identified were mapped in DisGeNET 50 to determine their known disease associations.
Core analysis using IPA. Ingenuity Pathway Analysis (IPA, Ingenuity Systems, USA; www. qiagen. com/ ingen uity) was used to identify the biological functions, diseases, canonical pathways, and regulatory networks of the functional miRNA-mRNA target interactions. Tab-delimited text files containing gene/miRNA IDs, expression data (fold change), and p-values were uploaded into IPA for their core analysis. The statistical significance of the enrichment was calculated using hypergeometric test and adjusted by FDR method (adj. p-value ≤ 0.05). The top functions (molecular, cellular and biological), diseases, toxicology, and gene signaling networks were calculated using IPA-generated negative logarithm p-values i.e., -log10(p-value) and associated Z-and network scores.

Construction of protein-protein interaction (PPI) network. To examine the interactive associations
among the DEGs at the protein level, MM related genes were mapped on protein-protein interaction (PPI) data using NetworkAnalyst 51 (version 3.0; http:// www. netwo rkana lyst. ca). The network was built based on the original seed proteins through executing the minimum interaction network by trimming the first-order network to keep only those nodes that are necessary to connect the seed nodes. Literature-curated comprehensive PPI data was used to predict interaction network 52 . Network modules containing densely connected group of proteins were predicted using the random walk approach. The significant p-value of a given module was calculated with Wilcoxon rank-sum test 53 . The enriched pathways of DEGs in significant modules (≥ 10 DEGs) were analysed with a threshold of p ≤ 0.05 using DAVID (database for annotation, visualization and integrated discovery) functional annotation tool.
Survival analysis. Sigmaplot 14.0 was used to estimate Kaplan-Meier plots for overall survival (OS) and progression free survival (PFS). The survival analysis was carried out on 35 patients in whom clinical data was available. Comparisons between DEMs were analyzed using means of log-rank test and p ≤ 0.05 as cut-off for statistical significance.
Ethical clearance. The study on 44 MM patients collected from the outpatient department of the All India Institute of Medical Sciences (AIIMS), New Delhi was conducted in compliance with ethical guidelines of the AIIMS and after obtaining approval from the AIIMS ethics committee. Study individuals were enrolled following their voluntary written informed consent.  (Fig. 2). DEGs and DEMs were further investigated for their involvement in most enriched diseases and for their functions in multiple myeloma. On annotation, most of the DEGs were found to be involved in cancer, organismal injury and abnormalities, immunological disease, connective tissue disorder, inflammatory disease (Figure S1a in Supplementary File 2), whereas DEMs were found to be enriched in cancer, organismal injury and abnormalities, reproductive system disease, inflammatory disease, and inflammatory response (Figure S1b in Supplementary File 2). The topmost significant diseases and biofunctions identified for DEMs and DEGs are shown in Table 2. Besides the leading pathways and cellular functions, gene networks were constructed to connect key genes and enriched categories of diseases and functions based on the correlation between DEGs. Core analysis-based network revealed 25 significant networks and each individual network had a maximum of 35 focus genes. Top ranked  Table 3 and Table S11 in Supplementary File 1.

Identification of functional modules in PPI network.
Protein-protein interactions (PPI) network was constructed using aberrantly expressed genes identified in MM to predict biologically significant modules containing a group of proteins that execute similar functions. The minimum interaction network scattered in 1-3 sub-networks including one big network with highest nodes and edges. The network analysis disclosed 1,136 seeds (91.61% of DEGs) associated with 1,937 nodes in the network. The modules containing a group of proteins with identical functions were detected using the random walk approach. A total of 22 significant independent functional modules were observed, whereas 13 modules (module no: 0, 1, 2, 3, 4, 5, 6, 7, 9, 11, 12, 13, and 16) were highly connected with more than 10 nodes and p ≤ 0.05 (Table 4; Table S12 in Supplementary File 1). Out of 1,136 seed nodes, a total of 34.68% (n = 394) nodes were observed with ≥ 10 degrees or connections with other nodes (Table S12 in Supplementary File 1). The betweenness centrality of nodes ranged between of 13.37 to 741,529.3 in the constructed network. All 394 nodes were observed to be targeted by at least one MM-associated DEMs (Table S13 in Supplementary File 1). The top five highly connected hub nodes included UBC, ITGA4, HSP90AB1, VCAM1 and VCP (Table S12 in Supplementary File 1). Module-wise distribution of top three highly connected hub nodes encompassed BRCA1, CDKN1A and PCNA in module 0 ( Figure S2a), UBC, ITGA4 and VCAM1 in module 1 ( Figure S2b)  Identification of biomarker candidates for multiple myeloma. The common molecular biomarker candidates among DEGs and DEMs for diagnosis, disease progression, efficacy, prognosis, response to therapy and safety were identified using the IPA software and HMDD/miRNet database, respectively ( Table 5). The analysis revealed 154 (12.42%) potential biomarkers out of 1,240 observed DEGs that could bear clinical value for MM and were designated as 'Union 154' signature (Fig. 3a). These included common biomarker candidates predominantly with diagnosis (n = 82; 63.25%), efficacy (n = 90; 58.44%), prognosis (n = 56; 36.36%), disease progression (n = 21; 13.64%), response to therapy (n = 23; 14.94%), and safety (n = 9; 5.84%) (Tables S14 and S15 in Supplementary File 1). Among the target gene candidate biomarkers, 42.21% (n = 65) of targets qualified for more than one role (Supplementary Tables S14 and S15). For example, gene CDKN2A was observed to be implicated in six biomarker applications including diagnosis, disease progression, efficacy, prognosis, response to therapy and safety.
In addition, miRNA disease databases such as HMDD and miRNet revealed 37 aberrantly expressed miRNAs as potential biomarkers with clinical utility for MM (Table S16 in Supplementary File 1) and were designated  . 3b). A systematic literature review of 'Union 37' signatures disclosed that 29.73% (n = 11) miRNAs were known circulating biomarkers for diagnostics and prognostics in MM. Some of these miRNAs were identified as epigenetically regulated miRNAs (n = 4), as therapeutic targets (n = 7) and dysregulated miRNAs that resulted in MM disease phenotype (n = 7) and are given in Table 6.

TF-gene/miRNA coregulatory networks.
We further investigated the TF-miRNA-target gene regulatory network for meta-signature gene/miRNAs identified in this study. The gene-TF regulatory network of 5 gene (≥ 4 biomarkers applications) revealed 164 interaction pairs among 5 seed genes (CDKN1A, MMP9, CDKN2A, MKI67, and IGF1) and 139 transcription factors (TFs) (Table S18 in Supplementary File 1). Among them, upregulated gene CDKN1A was found to be regulated by 87 TFs, CDKN2A was regulated by 12 TFs and IGF1 was regulated by 10 TFs (Fig. 5). Similarly, the downregulated gene MMP9 interacts with 20 TFs, and MKI67 interacts with 10 TFs (Fig. 5). TF-gene interactions are shown in Table S17 in Supplementary File 1. TF-miRNA regulatory analysis of top 3 up-and downregulated miRNA biomarkers based on the number of targets showed an association with 339 TFs (Table S19 in Supplementary File 1). From the data of TransmiR, downregulated DEMs such as hsa-miR-20a, hsa-mir-155, and hsa-mir-92a were found to be regulated by 143, 114 and 11 TFs, respectively ( Figure S3a-c). Likewise upregulated DEMs including hsa-mir-23b, hsa-mir-195 and hsa-let-7b were found to be regulated by 140, 63 and 58 TFs, respectively ( Figure S3d-f in Supplementary File 2; Table S17 in Supplementary File 1).

Discussion
In this study, a meta-analysis of mRNA and miRNA expression profiles has been carried out on more than 600 MM patients including 44 Indian myeloma patients (represented in 9 GSE-GEO datasets) in order to compute altered mRNA and miRNA patterns and potential biomarkers of prognostic clinical relevance in multiple     The common loss of miR-15a and miR-16-1 in CLL, as well as the loss of 13q14 in mantle cell lymphoma (50 percent of cases), multiple myeloma (16 to 40 percent) and prostate cancer (60 percent), strongly suggests that these two miRNAs act as tumor suppressor genes. While their full target complement is unknown, they appear to mediate their effects largely by down-regulating the anti-apoptotic protein BCL2. This protein is often found expressed at high levels in CLL and is thought to be important for the survival of the malignant cells    (Fig. 3a and b).
The present study has revealed that 85% (85/100) of DEMs and 91.04% (1,240/1,362) of DEGs were significantly altered, are inversely correlated and involved in regulatory networking in multiple myeloma. The most downregulated miR observed in MM malignant plasma cells as compared to NPCs in our study is miR-155. A reduced expression of this miR in MM PCs vs NPCs suggests a tumor suppressor role as has also been reported previously 55 . A similar study has reported an epigenetic repression of miR-375 in MGUS and MM primary cells as compared to NPCs 56 , which is also concurrent to our findings. Another tumor suppressor miR-144 that can be sponged by lncSOX2OT 57 has been reported to be downregulated in MM plasma cells and cell lines earlier and was found downregulated in plasma cells in our study. Similarly, upregulation of miR-29b in MM PCs in this study is in sync with previous studies, where it has been reported that the overexpression of miR-29b induces apoptosis of multiple myeloma cells by down regulating MCL-1 58 .
Some of the DEMs observed in MM in our study can be extrapolated and categorized on the basis of their previously reported roles relating to pathogenesis, clinical presentation, drug resistance and clinical outcomes. While the deregulated miRs-30d and 181b have been associated with p53 expression 33 , miRs-106/ 181b and miR-181b/ miR-193b are specifically dysregulated in early and late stages of pathogenesis in MGUS and MM respectively 30 . Some of the DEMs have been associated with sensitivities to Bortezomib (e.g., miRs-17-5p, miR-29b-3p, miR-20a-5p) while others with poor survival outcomes (miR-92a, miR-16, let-7e, miR-19b, miR-19a) 25 . Although sample size of inhouse MM subset (n = 44) in our study is small, we observed all the Union37 DEMs in this patient population. Moreover, a significant association of low expression of miR-30d-3p with poor OS and of high expression of miR-16-2-3p with poor OS and PFS (Fig. 4) was also noted. The miR-30d-3p is a   61 its deregulation may be of prognostic significance in MM and needs to be explored further. Coincidentally, IPA analysis has also highlighted importance of WNT pathway in this study. Another integrative study 27 mined two miRNA and two mRNA microarray GEO datasets and identified 39 DEMs and 32 hub genes. Among these DEMs, miR-155 and miR-148 were found to be deregulated in their study 27 as well as in Union 37 profile in the present work. Likewise, another meta-analysis of 7 datasets including MM patients 26 highlighted 13 DEMs, of which hsa-miR-106b, miR-15b, miR-191, miR-19a and miR-20a are also represented in Union 37 profile. A recent meta-analysis by Xu et al 25 reported 7 DEMs of poor prognostic significance among which deregulated miR-92a, miR-16, let-7e and 19b are common to the Union 37 signature.
The IPA core analysis disclosed 12.42% (n = 154) of DEGs as putative biomarkers that could be useful in diagnosis, disease progression, efficacy, prognosis, response to therapy and safety. Further investigation revealed that 42.21% (n = 65) of targets were involved in more than one functional role. It is known that proteins with the highest degree have the highest betweenness in the network. As hub proteins are accountable for holding networks together 62,63 , they are more likely to be master regulators of signaling and transcription and can be used as therapeutic targets or biomarkers 64 . The target genes identified in this study were subjected to PPI network which disclosed a total of 394 nodes with ≥ 10 connections with other nodes and were designated as 'hub' genes. All hub genes were observed to be targeted by MM associated DEMs and could act as possible biomarkers for this disease.
It is noteworthy that IPA based data mining of DEGs and DEMs in this study has revealed five top hub genes lying in the centre of functional networks. These include UBC, ITGA4, HSP90AB1, VCAM1 and VCP. Two genes (UBC and HSP90B1) have been earlier reported to be upregulated and involved in myelomagenesis in malignant plasma cells in other studies 65 as well and may be critically involved in ubiqutin-proteosomal pathway. The HSP90A family members are known to promote anti tumor immunity via their exposure on dying myeloma cells 66 and their interaction with lncRNA MALAT1 is associated with poor prognosis 67 . Gene ITGA4 along with ITGB1 codes for integrin VLA4 that mediates homing of myeloma cells into bone marrow and augment IL6 in the microenvironment. 68 . Similarly, MM cells establish contact with bone marrow stromal cells via adhesion molecules such as VCAM1 and enhance osteoclast stimulating activity that can be reduced by Bortezomib and Lenalidomide 69,70 . The gene VCP is a potential therapeutic target that mediates delivery of ubiquinated misfolded protein aggregates to proteasome 71 and was found to be upregulated in MM plasma cells in this study.

Conclusions
The regulatory crosstalk between DEGs and DEMs in MM is highly complex. This study has identified core putative signatures of DEMs ('Union 37') and DEGs ('Union 154') in MM as compared to normal PCs that may impact clinical outcomes (for instance, miR-16-2 and miR-30d). Further studies on functionally connected hub genes (such as UBC, ITGA4, HSP90AB1, VCAM1, VCP), other potential seed genes (e.g., CDKN1A, CDKN2A, MMP9, IGF1, MKI67), DEMs and their multidimensional networking with regulatory transcription factors are needed for better understanding of their oncogenic/ anti tumor properties and to explore their synergistic prognostic value.

Data availability
Gene expression (GSE125361) and miRNA expression (GSE125363) signatures in multiple myeloma have been submitted to the National Center for Biotechnology Information (NCBI; https:// www. ncbi. nlm. nih. gov/ geo) under BioProject accession number PRJNA515992.