Data Collection

This study utilized the WXS, WGS and RNA-Seq data of NDMM and RRMM patients. The RNA-Seq and variant data was retrieved from National Cancer Institute Genomic Data Commons (NCI-GDC) data Portal utilizing the MMRF CoMMpass study data (https://portal.gdc.cancer.gov/projects.) GDC is a central repository, publicly available to access the vast amount of clinical and genomic data of numerous cancer studies. The MMRF CoMMpass dataset consists of WXS, WGS, RNA-Seq and clinical data of 995 MM patients primarily at the diagnosis and then at the relapse of cancer .

Shortlisting of Candidate Relapse Biomarkers

The genes were selected as candidate biomarkers if meeting any one of the following criteria 1) a gene retrieved as relapse biomarker in literature mining and also upregulated in DGE, 2) a gene carrying SNV and also upregulated in DGE, 3) Frequently mutated genes with higher no of SNVs.

Immune Cells Infiltration

The immune cells infiltration of the RRMM vs NDMM patients was carried out by utilizing the immunedeconv 2.1.0 package in R 4.2.2 [13]. This method estimates the count of immune cells present in the samples by employing deconvolute functions. The gene expression matrix was normalized to obtain fragments per kilobase of transcript per million mapped reads (FPKM) using the count2FPKM function of RNAAgeCalc 1.10.0 package  [14]. The FPKM normalized matrix was further subjected to deconvolution methods to perform immune cells infiltration. Here we utilized two deconvolution methods i.e., quanTIseq and MCPcounter, whereas the ggplot2 version 3.4.1 was used to visualize the results in the form of bar chart and dot plot  [13,15]. The quanTIseq utilizes bulk RNA-seq data to estimate fraction of 10 immune cell types along with the uncharacterized cells. Whereas MCPcounter obtains transcriptomics data to provides the absolute count of immune and stromal cells across all samples with respect to each cell type individually.

Literature Mining

Literature mining was conducted to identify biomarkers involved in recurrence of multiple types of cancers using the search terminologies on Google Scholar (https://scholar.google.com)  and PubMed (https://pubmed.ncbi.nlm.nih.gov) as “recurrence in cancer”, “cancer relapse, “recurrent genes involved in lung cancer”, “recurrent genes involved in pancreatic cancer”, “recurrent genes involved in bladder cancer”, “recurrent genes involved in breast cancer”, “recurrent genes involved in colorectal cancer” and “recurrent genes involved in prostate cancer”. Furthermore, genes associated with relapse of MM were also identified through literature review using the terminologies “Relapse” “RRMM”, “cancer relapse” and “multiple myeloma” in Google Scholar and PubMed. Lastly the common genes obtained through literature mining and DGE was identified using grep tool on Linux.

Differential Gene Expression

Differential gene expression of RNA-seq data of CoMMpass was performed by using DESeq2 (1.38.3), an R (4.2.2) package [16]. This method normalizes the counts of each gene and employs shrinkage estimation resulting in more stable and explainable calculations focused on strength based quantitative analysis [16]. The DEGs were selected on the basis of log2 fold change (Log2FC) and p-value. Log2FC and p-value were set at the threshold of 2.0 < Log2FC < -2.0 and < 0.05, respectively to filter out up and downregulated DEGs. Furthermore, the EnhancedVolcano 1.16.0 package was used to retrieve EnhancedVolcano plot of DEGs  [17].

Function Enrichment

Functional enrichment analysis of DEGs and shortlisted candidate biomarkers was performed by GeneCodis 4 (https://genecodis.genyo.es/). GeneCodis 4 is an online enrichment analysis tool that extracts biological information from the experimental data by giving the statistical score for significantly enriched annotations [18].The lists of upregulated and downregulated genes were taken separately as input for annotations specified as Gene Ontology (GO) terms (biological process (BP), cellular component (CC), molecular function (MF)) and KEGG pathway.

After applying the filter of “BM samples” and “non-synonymous SNV” on MMRF- CoMMpass dataset 833 patients were retrieved. Among them 753 were NDMM patients containing 2007 SNVs in 454 genes whereas 80 were RRMM patients along with 226 SNVs affecting 152 genes. The top 5 Missense SNVs maximum affecting cases in cohort of NDMM and RRMM were Q61R in NRAS (7.56%), Q61H in KRAS (6.72%), Q61K in NRAS (5.40%), V640E in BRAF (3.48%) and G12D in KRAS (3.36%). 

Identification of Immune Cells Fraction in MM Samples

The results of immune cell infiltration analysis were assembled in Figure 1a and b. The meticulous analysis of results revealed that, T cell (CD4), natural killer cells (NK), monocytes, macrophages and myeloid dendritic cells (MDC) counts was marginally raised for RRMM in comparison to NDMM. Cytotoxicity score, endothelial cells and cancer associated fibroblasts (CAFs) can only be calculated through MCPcounter showed the similar trend of slightly higher count in RRMM as compared to NDMM. Similarly, regulatory T cell (Tregs) can only be estimated by quanTIseq was somewhat more in RRMM than NDMM. The only discrepancy was observed in the count of B cells as it was more at RRMM than NDMM according to MCPcounter whereas it was higher in NDMM than RRMM according to quanTIseq. Moreover, T cell (CD8+) count declined, and Neutrophils count amplified in RRMM in comparison to NDMM according to both methods (Figure 1a and b).

DGE Analysis

The gene expression profiling of the BM samples of NDMM and RRMM patients revealed a total of 1562 dysregulated genes. Among them 908 genes were significantly upregulated in RRMM (Log2FC > 2, p-value < 0.05) whereas 654 genes were downregulated (Log2FC < -2, p-value < 0.05) as shown in Figure 1c. Furthermore, the top 10 upregulated and downregulated genes were filtered out on the basis of highest and lowest Log2FC values are assembled in Table 1. The top 10 significantly upregulated DEGs include ACTL8, NPTX1, MAGEA4, KLRF2, EGFLAM, UNC13C, GABRA5, NPFFR2, LINC01414, OR7E22P. Whereas the top10 significantly downregulated DEGs were FGFR3, SSTR1, RBM24, LINC03020, HAPLN1, PEG3, IRS4, SLITRK2, RNF128, AGR2.

Literature Mining

Literature mining retrieved 136 genes with significant role in relapse of various cancers as shown in Table 2. Among them 40 genes were found to be involved in RRMM (Table 2). Analysis of results revealed that KRAS, NRAS, TP53, NF1, STK11 and DNMT3A that were found as relapse biomarker in MM also reported for relapse in other cancers. Although these biomarkers were not found to be upregulated in our DGE analysis but mutations in these gene is reported by various research in RRMM.

Shortlisting of Candidate Genes

Genes were shortlisted on the basis of following criteria 1) a gene retrieved as relapse biomarker in literature mining and also upregulated in DGE, 2) a gene carrying SNV and also upregulated in DGE, 3) Top mutated genes retrieved through SNV analysis. According to first criteria 6 DEGs were shortlisted including Colony-Stimulating Factor-1 Receptor (CSF1R), VCAN, Neuropilin 1(NRP1), Collagen Type XXII Alpha 1 Chain (COL22A1), Bactericidal Permeability Increasing Protein (BPI) and Baculoviral inhibitor of apoptosis repeat-containing 5 (BIRC5). These DEGs were upregulated in our analysis but reported as relapse biomarkers for other cancers (Table 3) in literature. Similarly, 5 DEGs were selected in accordance with second criteria included Motor Neuron and Pancreas Homeobox 1 (MNX1), FAT Atypical Cadherin 1 (FAT1) ETS Transcription Factor ERG (ERG), TCL1 Family AKT Coactivator A (TCL1A) and ALF Transcription Elongation Factor 3 (AFF3) (Table 3). These DEGs showed upregulation along with SNV in RRMM patients. These DEGs are reported as diagnostics, prognostics, therapeutic and relapse biomarkers for various other cancers but their role in relapse of MM needs to be explored.   Moreover, two genes (KRAS, NRAS) were selected according to third criteria although, didn’t show any significant DE at relapse but harbor multiple SNVs in both NDMM and RRMM patients (Table 3).

Functional Enrichment and Pathway Analysis of DEGs and Shortlisted Candidate Relapse Biomarker

The GO term analysis of DEGs revealed that upregulated DEGs were found to be significantly enriched in the following BP, immune response, immune system process, inflammatory response, cell adhesion, positive regulation of T cell activation, antigen processing and presentation of peptides, peptide antigen assembly with MHC class II, immunoglobulin production, antigen processing and presentation of exogenous peptides and neutrophil chemotaxis. However, the downregulated DEGs were notably associated with adaptive immune response, positive regulation of B cell activation, phagocytosis recognition and engulfment, complement activation and development of central nervous system among many others. Similarly, the significant CC terms for upregulated DEGs were plasma membrane, extracellular region and extracellular space whereas, collagen-containing extracellular matrix, external side of plasma membrane, tertiary granule membrane, cell surface, ficolin-1-rich granule lumen and MHC class II protein complex were also notable terminologies. However, the downregulated DEGs were also enriched in plasma membrane, extracellular region and extracellular space, collagen-containing extracellular matrix, external side of plasma membrane along with immunoglobulin complex, synapses, and glutamatergic synapse. Additionally, the significant MF retrieved for the upregulated DEGs through GO analysis were carbohydrate binding, signaling receptor activity, transmembrane signaling receptor activity and cytokine activity. While the antigen binding, sequence-specific DNA binding was the most significant MF related to downregulated DEGs among many others.

The pathway analysis of DEGs revealed their enrichment in various essential pathways (Figure 2). The upregulated DEGs were found significantly enriched in numerous disease pathways including staphylococcus aureus infection, rheumatoid arthritis, malaria, Tuberculosis, asthma, leishmaniasis, graft versus host disease along with phagosome and hematopoietic cell lineage (Figure 2a). The downregulated DEGs, on the other hand, were found to be enriched in diverse signaling pathways including cAMP signaling, calcium signaling, hippo signaling, signaling pathways regulating pluripotency of stem cells along with Neuroactive ligand-receptor interaction and various synapses as depicted in Figure 2b.

Furthermore, the functional enrichment of these shortlisted candidate genes showed significant enrichment of the following BP ontologies,  positive regulation of cell population proliferation, multicellular organism development, positive regulation of serene/threonine kinase activity and positive regulation of endothelial cell proliferation more significant among many others. The shortlisted DEGs showed CC enrichment more significantly in focal adhesion whereas, G protein activity, GDP-binding and cytokine binding were enriched MF. Moreover, KEGG pathways analysis revealed significant enrichment in prostate cancer, colorectal cancer (CRC), acute myeloid leukemia, chemical carcinogenesis receptor activation and PI3K/AKT signaling pathway (Figure 2c).

Figures & Tables

Discussion

The relapse in MM is a constant concern for researchers because of complex multifactorial onset, prognosis and even more intricate recurrence of cancer. Development of effective treatment regimen to minimize relapse and maximize progression free survival of MM patients is conditional with the hunt of new biomarkers for diagnosis, prognosis and relapse along with efficient and precise drug targets and compounds. Literature review of six genes (CSF1R, VCAN, NRP1, COL22A1, BPI and BIRC5) identified by combine analysis of literature mining and DGE revealed their potential in modulating TME and disease progression. CSF1R, a receptor tyrosine kinase, inhibition is reported as potential antitumor strategy. Activated CSF1R are recruited by tumor-associated macrophages and release cytokine that modulated TME to protumoral phenotype. CSF1R activation also initiates many downstream pro survivals signaling cascades including PI3K/AKT, ERK1/2, and JNK [19]. Similarly, BIRC5 is also an immune-related gene that inhibits apoptosis and promotes cellular proliferation. High expression of BIRC5 regulates DNA methylation hence is reported as potential target for developing immunotherapies [20]. Moreover, NRP1 is an independent predictor of relapse and poor survival in non-small cell lung cancer (NSCLC). It is also reported as novel potential therapeutic target in NSCLC because of its critical role in tumorigenesis, cancer invasion, and angiogenesis through VEGF, PI3K, and Akt pathways [21]. Whereas VCAN and COL22A1 are reported as prognostic biomarkers in various cancers. VCAN mRNAs are specifically expressed in cancer-associated fibroblasts and associated with poor relapse free survival of stage II-III patients in CRC. It is a promising biomarker to identify stage II-III patients with high risk of relapse in CRC [22]. While COL22A1 is an integral part of novel prognostic immune related gene signature in CRC [23]. Additionally, COL22A1 is also reported as poor prognosis and relapse biomarker in Head and neck squamous cell carcinoma (HNSCC)   [24]. Furthermore, BPI is associated with human neutrophil as in response to inflammation neutrophils secrets BPI along with many other cytotoxic proteins and  can promote tumor metastasis by the formation of so-called neutrophil extracellular traps (NETs).

High levels of circulating and intratumoral neutrophils have been shown to correlate with poor survival in pancreatic cancer and pancreatic ductal adenocarcinoma [25]. These results were also consistent with our immune cell infiltration analysis which revealed the presences of higher neutrophils count in RRMM. The immune cells infiltration analysis in our study corroborated the previously reported role of immune dysfunction in MM invasion and progression [19]. Development of effective therapeutic strategies against neutrophils for treating cancer and various other diseases has already been in consideration. According to recent studies neutrophils have multiple phenotypes that perform diverse functions, particularly modulation of inflammation and immune response [26]. Similarly, another study mentioned the impairment of immune response in MM due to reduced phagocytic activity of neutrophils in comparison to healthy control [27]. Moreover, the role of neutrophils in facilitating tumor progression through immune deregulation and increased vulnerability to infection in MUGS and MM has already been established [28]. CD8+ T cells count was found to be decreased in NDMM, are imperative for immune defense against intracellular pathogens (viruses, bacteria) and for tumor surveillance. The trafficking or transporting of CD8+ T cells into the TME is crucial to exert its anti-tumor function. The elevated levels of CD8+ T cells in the TME are linked with positive anti-tumor effects and good prognosis in breast, colorectal, glioblastoma, and cervical cancers [29]. Hence the lower count of CD8+ T in TME is not only linked with cancer progression/relapse or poor prognosis but also with the increased risk of secondary viral or bacterial infection. The lower count of CD8+ T cells retrieved through immune infiltration investigation also supports our findings of pathway analysis which showed enrichment of various other cancer and bacterial infections pathways in RRMM (Figure 2).

The DGE and SNV analysis of dataset retrieved five genes (MNX1, FAT1, ERG, TCL1A, AFF3) harboring mutation and upregulation in RRMM. These mutations include MNX1 (P392L), FAT1 (N3716K), ERG (E353Q), TCL1A (T38I), and AFF3 (P1129L).  KRAS and NRAS, although not found in DGE but are the most mutated genes in both NDMM and RRMM patients according to SNV data. Both KRAS and NRAS are proto-onco genes, belonging to the Ras family of proteins, small GTPase, involved in regulation of biological processes particularly cell growth, proliferation and apoptosis  [30]. It has been reported that newly acquired mutations and pre-treatment sub-clonal mutations of KRAS and NRAS in MM possibly induce chemo resistance and relapse [31]. Mutation in KRAS and NRAS have been reported to be more common in relapsed patients (>70% patients) [32,33]. The codons G12, Q61 and G13 are mutation hotspots for both KRAS and NRAS but various substitutions in each codon elicit different signaling pathway, hence express distinct pathophysiology [34]. In our study all the three hotspots were found substituted with various codon (NRAS: Q61H, Q61R, Q61E, G12V, G12D, G12R, G13D KRAS: Q61K, Q61H, G12D, G13R, G13D) along with few new codons (NRAS: K117N, A59E KRAS: Y64D, E153Q) in RRMM. The impairment of PI3K pathway due to KRAS (G12R) mutation is reported in literature [35] and is also consistent with our results of pathway analysis (Figure 4).  Additionally, MNX1 encodes a transcription factor (HB9) that contains a homeodomain. The overexpression of MNX1 has been reported in many cancers (prostate, colon, liver, breast and bladder cancer, glioma and pancreatic progenitor tumors and acute myeloid leukemia) and has been suggested as potential diagnostic and prognostic biomarker [36, 37]. Upregulation of MNX1 stimulates the Wnt/β‐catenin signaling and via expression of downstream genes c‐Myc and CCND1, hence plays a vital role in CRC progression [38].  FAT1 is among the most frequently mutated genes in many types of cancer.  The role of FAT1 in cancer progression is highly dependent upon cancer type. In some cancers epithelial-mesenchymal transition (EMT) and the formation of cancer initiation/stem-like cells is promoted by loss of FAT1 function promotes whereas overexpression of FAT1 leads to EMT in others. The paired analysis of diagnosis and relapse sample in B-cell acute lymphoblastic leukemia overexpression of FAT1 was found correlated with shorter relapse-free and overall survival [39]. Several studies have reported a correlation of FAT1 mutation or expression with prognosis in various cancers, such as breast cancer, NSCLC, gastric cancer and T-cell lymphoma [40].  Positive correlation of FAT1 overexpression with proliferation and WNT/β‐catenin signaling pathway in T-cell acute leukemia (T-ALL) is recently reported in a study [41]. ERG encodes a transcription factor involved in development and differentiation affecting vasculogenesis, haematopoiesis, angiogenesis and embryogenesis, and is associated with regulation of cellular processes  [42]. ERG over expression is associated with poor prognosis and oncogenesis promotion, in prostate cancer, Ewing’s sarcoma, acute myeloid leukemia, acute T-lymphoblastic leukemia [43]. ERG high expression stimulates gene fusion event (ERG-TMPRSS2) that leads to early relapse in AML [43]. ERG upregulation is associated with upregulation of PI3K/AKT pathway [44]. TCL1A is a proto-oncogene expressed in embryonic stem cells, activated T and B lymphocytes and coactivator of kinases and interacting partners crucial in signaling pathways (PI3K and NF-κB) and cellular activities  [45]. Dysregulated TCL1A has a well-documented role in hematopoietic malignancies i.e., development of T-cell leukemia, correlation of overexpression with aggressiveness, deregulation of the cell cycle and genomic instability in chronic lymphocytic leukemia [46]. It is also suggested as prognostic biomarker for stage II/III CRC [47], therefore, proposed as potential biomarker for colorectal and hematological malignancies [45, 47]. AFF3 is primarily expressed in B cells and encodes a protein involved in transcription regulation [48]. AFF3 upregulation has been found in many cancers (gastric, breast, Adrenocortical and AML), also involved in modulating TME thus suggested as target for immune therapy [48-50]. In gastric cancer, dysregulated AFF3 is a potential marker for diagnosis and prognosis as well as correlated with immune checkpoints response whereas in breast cancer  overexpression is associated with drug (tamoxifen) resistance therefore suggested as  predictive marker for ER+ breast cancer [48, 50]. Moreover, AFF3 meditates oncogenic effects of β-catenin as constitutive activation of Wnt/β-catenin signaling pathway in mice leads to formation of malignant adrenocortical tumors [51].

The functional enrichment analysis of all shortlisted genes revealed their significant enrichment in G protein activity, GDP binding, positive regulation of cell population proliferation, serene/threonine kinase activity, endothelial cell proliferation. However, pathway enrichment for PI3K-Akt signaling and prostate cancer, acute myeloid leukemia and CRC pathways along with many cancer pathways was also significant (Figure 2). It can be inferred that KRAS and NRAS mutants may affect the activation and hydrolysis due to dysregulation of G protein activity, GDP and GTP binding, and GTPase activity. The aberrant GTPase activity due to KRAS mutations has been reported to affect the GTP-hydrolysis [52]. Furthermore, enrichment of PI3K-Akt signaling pathway was also consistent with our other results as upregulation of CSF1R, NRP1, TCL1A and ERG activate PI3K-Akt signaling [5, 44].  This pathway is crucial to many cellular processes and plays a significant role in cancer proliferation and multidrug resistance. Decrease in cellular apoptosis is mediated by continuous phosphorylation of various transcription factors by AKT, thus the promotion of the proliferation, angiogenesis, and survival of cell [53]. However, the PI3K-Akt signaling pathway is crucial for pathophysiology of MM and is associated with therapy resistance  [54, 55]. Moreover, Wnt/β‐catenin signaling plays a dual and disease stage-specific role in the pathogenesis of MM. Wnt/β‐catenin pathway activation during MM disease progression is mediated through epigenetic silencing by antagonists which facilitates pathway activation and proliferation of MM cells [56]. The pathway although not found enriched in pathway analysis but the 3 DEGs (MNX1, AFF3, and FAT1) among the 5 selected candidate genes are found to have direct relationship in upregulation of wnt/ β‐catenin pathway [38, 51]. This study proposed that upregulation of the following genes (CSF1R, VCAN, NRP1, COL22A1, BPI,  BIRC5, MNX1, FAT1, ERG, TCL1A, AFF3) might be the cause of drug resistance and relapse in RRMM. These genes exert their effect by modulating apoptosis, TME and PI3K-Akt signaling pathway mainly along with various other pathways. Their role as potential biomarker for diagnosis, prognosis and drug resistance in various cancers has also been reported. Thorough investigation on role of these genes in drug resistance might be helpful in understanding the mechanism of relapse in MM.  Furthermore, the neutrophils, P13K-Akt signaling pathway along with aforementioned genes are potential biomarkers of relapse in MM and targets of new effective treatment regimens. 

Conflict of Interest

The authors declare that there is no conflict of interest.

Author Contributions

All authors contribute equally to designing study, acquisition of data, analysis and interpretation of data, and drafting of the manuscript.

References