Abstract Objective To uncover the potential regulatory mechanisms of the relevant genes that contribute to the prognosis and prevention of multiple myeloma (MM). Methods Microarray data ([28]GSE13591) were downloaded, including five plasma cell samples from normal donors and 133 plasma cell samples from MM patients. Differentially expressed genes (DEGs) were identified by Student’s t-test. Functional enrichment analysis was performed for DEGs using the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Transcription factors and tumor-associated genes were also explored by mapping genes in the TRANSFAC, the tumor suppressor gene (TSGene), and tumor-associated gene (TAG) databases. A protein–protein interaction (PPI) network and PPI subnetworks were constructed by Cytoscape software using the Search Tool for the Retrieval of Interacting Genes (STRING) database. Results A total of 63 DEGs (42 downregulated, 21 upregulated) were identified. Functional enrichment analysis showed that HLA-DRB1 and VCAM1 might be involved in the positive regulation of immune system processes, and HLA-DRB1 might be related to the intestinal immune network for IgA production pathway. The genes CEBPD, JUND, and ATF3 were identified as transcription factors. The top ten nodal genes in the PPI network were revealed including HLA-DRB1, VCAM1, and TFRC. In addition, genes in the PPI subnetwork, such as HLA-DRB1 and VCAM1, were enriched in the cell adhesion molecules pathway, whereas CD4 and TFRC were both enriched in the hematopoietic cell pathway. Conclusion Several crucial genes correlated to MM were identified, including CD4, HLA-DRB1, TFRC, and VCAM1, which might exert their roles in MM progression via immune-mediated pathways. There might be certain regulatory correlations between HLA-DRB1, CD4, and TFRC. Keywords: multiple myeloma, functional enrichment, transcription factor analysis, PPI network, pathway enrichment Introduction Multiple myeloma (MM) is a malignant plasma cell (PC) disorder and is characterized by the infiltration of malignant PCs, which are the endpoint of B-cell differentiation, in the bone marrow.[29]^1^–[30]^4 The mortality of MM is relatively high; approximately 63,000 individuals die from MM annually, accounting for 0.9% of all cancer deaths.[31]^5 Therefore, many studies have been conducted to explore the molecular mechanisms of MM etiology.[32]^6^,[33]^7 MM is known to be associated with heterogeneous cytogenetic abnormalities, including hyperdiploid and nonhyperdiploid karyotypes.[34]^8 In patients with a nonhyperdiploid karyotype, many chromosomal translocations are involved, which could lead to the altered expression of related proteins (cyclin D1 and cyclin D3) and genes (WHSC1 and MAF).[35]^4 Additionally, the interaction of genes and proteins are also implicated in the pathogenesis of MM. For instance, a previous study demonstrated that IL-6 contributes to the pathogenesis of MM and that the concentrations of insulin-like growth factor binding protein-1 and the soluble IL-6 receptor are positively associated with the risk of MM.[36]^9 Despite these significant discoveries, due to the difficulty of obtaining the metaphases of malignant PC clone, there is limited knowledge regarding the genetics of MM due to the difficulty of obtaining malignant PC clones in metaphase. Fortunately, the application of advanced molecular techniques such as microarrays and next-generation sequencing facilitate the improvement of our understanding from a genetic level[37]^10 Egan et al[38]^11 and Chapman et al[39]^12 have analyzed the genomic events that initiate MM using the genomic sequencing. Microarray data ([40]GSE13591) were established by Agnelli et al[41]^13 who analyzed these data using a combined FISH and microarray approach and identified near-tetraploidy as a hallmark of the tumor. They also highlighted that loss of heterozygosity is a prominent mechanism in the regulation of mRNA and gene expressions. In their further studies based on these microarray data, the authors provide an elaborate elucidation of the Hedgehog pathway in MM, laying an elaborate foundation for the use of Hedgehog inhibitor detection in clinical trials.[42]^14 Other investigations based on these microarray data ([43]GSE13591) have also been carried out. Lionetti et al[44]^15 defined the microRNA/mRNA regulatory network of MM. HOXB7 expression in MM is identified by both microarray data ([45]GSE13591) analysis and further verification of experiments.[46]^16^,[47]^17 Moreover, another study reconstructed gene regulatory networks by combining gene expression profiles from seven publicly available datasets; however, the emphasis of this study was to identify critical genes that predict overall survival in the prognosis of MM.[48]^18 Nevertheless, none of the aforementioned studies further analyzed the functions, pathways, or the potential correlations of the identified differentially expressed genes (DEGs) between normal and MM PCs. For this study, we downloaded the microarray data ([49]GSE13591) and reanalyzed them using bioinformatics methods including the identification of DEGs in MM, as well as performed the Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), protein–protein interaction (PPI), and PPI subnetwork analyses of these DEGs. These approaches are based on multiple bioinformatic packages, which employ powerful statistical methods,[50]^19^,[51]^20 thus facilitating the identification of critical DEGs and their pathways that are involved in the development of normal PCs into MM. In addition, along with predicting the interactions of the DEGs at the protein level, we aimed to uncover the potential regulatory mechanisms of the relevant genes that contribute to the prognosis and prevention of MM. Methods Microarray data and data preprocessing The gene expression profile [52]GSE13591, deposited by Agnelli et al[53]^13 was downloaded from the Gene Expression Omnibus (GEO) database.[54]^21 The microarray platform was the [55]GPL96 (HG-U133A) Affymetrix Human Genome U133A Array (Affymetrix, Inc., Santa Clara, CA, USA). The expression profile included 138 samples, comprising five PC samples obtained from normal donors (control group) and 133 PC samples derived from MM patients (MM group). The raw expression profile data were preprocessed using the Affy package in Bioconductor[56]^22 and the Affymetrix annotation files from the Brain Array Lab (Microarray Lab, University of Michigan, Ann Arbor, MI, USA). The background correction, quantile data normalization, and probe summarization were performed by the Robust Multi-array Average algorithm to obtain a gene expression matrix.[57]^23 Identification of DEGs The DEGs between the MM and control groups were identified based on Student’s t-test in the linear models for the microarray data (Limma) package.[58]^24 Raw P-values were adjusted into false discovery rates (FDRs) by the Benjamini and Hochberg method.[59]^25 The thresholds for DEG selection were FDR <0.05 and |log[2]FC (fold change)| >1. Functional enrichment analysis of DEGs The GO[60]^26^,[61]^27 and KEGG[62]^28 pathway[63]^29 enrichment analyses were performed to identify the biological functions and pathways of the identified DEGs, based on the hypergeometric distribution algorithm. P<0.05 was chosen as the cut-off value for enriched functions and pathways. Gene function annotation of DEGs By mapping the DEGs using the TRANSFAC database,[64]^30 the transcription factors among the DEGs were identified. Combined with information from the tumor suppressor gene (TSGene)[65]^31 and tumor-associated gene (TAG) databases,[66]^32 all known oncogenes and tumor suppressor genes among the DEGs were extracted. PPI network construction and module analyses The DEGs were mapped using the Search Tool for the Retrieval of Interacting Genes (STRING)[67]^33 database[68]^34 to reveal potential correlations between the DEGs at the protein level. Only the interaction pairs that contained at least one DEG were selected to construct the PPI network, with the cut-off value of a confidence score >0.85. Cytoscape software[69]^35 was used to visualize the PPI network. Furthermore, the PPI subnetwork was identified using BioNet,[70]^36 and FDR <0.01 was chosen as the cut-off value. The KEGG pathway enrichment analysis of the genes in the PPI subnetwork was performed. Results DEGs between the MM and control groups A total of 166 transcripts, corresponding to 63 DEGs, including 42 downregulated DEGs and 21 upregulated DEGs, were identified between the MM and control groups, as shown in the heat map of the cluster analysis of DEGs ([71]Figure 1). Figure 1. [72]Figure 1 [73]Open in a new tab The heat map of the cluster analysis of differentially expressed genes. Note: The X-axis represents samples, while the Y-axis represents genes. Enrichment analysis of DEGs GO enrichment analysis revealed that 262 functional terms for downregulated DEGs and 78 for upregulated DEGs were enriched; the top ten functional terms ranked by statistical significance were listed in [74]Table 1 (for down- and upregulated DEGs, respectively). As presented in [75]Table 1, the top three significantly enriched GO terms for downregulated DEGs were regulation of immune system processes, immune response, and positive regulation of immune system processes (in which HLA-DRB1, VCAM1, A2M, CXCL12, and HMOX1 might be involved; P=9.99E-16), whereas the top three for upregulated DEGs were primary metabolic processes (in which ATF3, CCNB1IP1, DBF4, DUSP5, and EEF1E1 may participate; P=0.0008866), organic substance metabolic processes (P=0.0014942), and response to axon injury (P=0.0016905). Table 1. Top ten enriched GO of downregulated DEGs and upregulated DEGs DEGs GO ID Term Gene counts P-value Genes Downregulated GO:0002682 Regulation of immune system process 21 0 HLA-DRB1, VCAM1, A2M, CXCL12, HMOX1 GO:0006955 Immune response 23 0 HLA-DRB1, VCAM1, A2M, CXCL12, HMOX1 GO:0002684 Positive regulation of immune system process 17 9.99E-16 HLA-DRB1, VCAM1, A2M, CXCL12, HMOX1 GO:0050776 Regulation of immune response 17 8.33E-15 HLA-DRB1, VCAM1, A2M, CXCL12, HMOX1 GO:0002376 Immune system process 23 3.26E-14 HLA-DRB1, VCAM1, A2M, HMOX1, CD79A GO:0002694 Regulation of leukocyte activation 13 7.57E-14 HLA-DRB1, VCAM1, HMOX1, CD27, HLA-DMA GO:0045321 Leukocyte activation 15 9.87E-14 HLA-DRB1, VCAM1, HMOX1, CD79A, CXCL12 GO:0046649 Lymphocyte activation 14 2.09E-13 HLA-DRB1, VCAM1, CD79A, CXCL12, CD27 GO:0050865 Regulation of cell activation 13 2.10E-13 HLA-DRB1, VCAM1, HMOX1, CD27, HLA-DMA GO:0051251 Positive regulation of lymphocyte activation 11 3.15E-13 HLA-DRB1, VCAM1, CD27, HLA-DMA, HLA-DPA1 Upregulated GO:0044238 Primary metabolic process 18 0.0008866 ATF3, CCNB1IP1, DBF4, DUSP5, EEF1E1 GO:0071704 Organic substance metabolic process 18 0.0014942 ATF3, CCNB1IP1, DBF4, DUSP5, EEF1E1 GO:0048678 Response to axon injury 2 0.0016905 TXN, UCK2 GO:0008152 Metabolic process 18 0.0032882 ATF3, CCNB1IP1, DBF4, DUSP5, EEF1E1 GO:0033554 Cellular response to stress 6 0.0039691 ATF3, EEF1E1, GADD45A, MCL1, PPP1R15A, TXN GO:0044237 Cellular metabolic process 17 0.0044393 ATF3, CCNB1IP1, DBF4, DUSP5, EEF1E1 GO:0006259 DNA metabolic process 5 0.0045266 DBF4, GADD45A, HIST1H1C, HIST1H2BK, KPNA2 GO:0030968 Endoplasmic reticulum unfolded protein response 2 0.0054425 ATF3, PPP1R15A GO:0034620 Cellular response to unfolded protein 2 0.0055648 ATF3, PPP1R15A GO:0035967 Cellular response to topologically incorrect protein 2 0.0061954 ATF3, PPP1R15A [76]Open in a new tab Abbreviations: GO, Gene Ontology; DEGs, differentially expressed genes; ID, identifier. The KEGG pathway enrichment analysis revealed that 22 pathways for downregulated DEGs and two pathways for upregulated DEGs were enriched. The top ten significant pathways for downregulated DEGs and the two pathways for upregulated DEGs were shown in [77]Table 2. Among them, the top three KEGG pathways for downregulated DEGs were intestinal immune network for immunoglobulin (Ig) A production (P=5.55E-14), antigen processing and presentation (P=7.19E-14), and asthma (P=7.84E-14). The enrichment result implied that the genes HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, and HLA-DRB4 might be involved in the top three pathways. The two enriched pathways for upregulated genes were the mitogen-activated protein kinase signaling pathway (P=0.0088912) involving DUSP5, GADD45A, and JUND, and the cell cycle (P=0.0178289), involving DBF4 and GADD45A. Table 2. The top ten enriched KEGG pathways of downregulated DEGs and the two enriched KEGG pathways of upregulated DEGs DEGs KEGG pathway Gene counts P-value Genes Downregulated Intestinal immune network for IgA production 9 5.55E-14 CXCL12, HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Antigen processing and presentation 10 7.19E-14 CTSB, HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA- DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4, HSPA1A Asthma 8 7.84E-14 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Staphylococcus aureus infection 9 2.08E-13 C1QA, HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Allograft rejection 8 5.09E-13 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Graft-versus-host disease 8 1.25E-12 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Type 1 diabetes mellitus 8 1.89E-12 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Autoimmune thyroid disease 8 9.59E-12 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Rheumatoid arthritis 9 2.37E-11 CXCL12, HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Viral myocarditis 8 1.15E-10 HLA-DMA, HLA-DPA1, HLA-DPB1, HLA-DQA1, HLA-DQA2, HLA-DRB1, HLA-DRB3, HLA-DRB4 Upregulated MAPK signaling pathway 3 0.0088912 DUSP5, GADD45A, JUND Cell cycle 2 0.0178289 DBF4, GADD45A [78]Open in a new tab Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; DEGs, differentially expressed genes; IgA, immunoglobulin A; MAPK, mitogen-activated protein kinase. Analysis of gene function annotation By integrating the information in the TRANSFAC database, there were three transcription factors identified, including one downregulated gene (CEBPD) and two upregulated genes (JUND and ATF3). Additionally, three TAGs, including one downregulated gene (CEBPD) and two upregulated genes (JUND and EEF1E1), were selected. PPI network and subnetwork analysis The PPI network of genes was constructed, as shown in [79]Figure 2. The top ten genes/proteins with the highest degrees (the numbers of PPI pairs that the gene/protein possessed) in the PPI network were HLA-DRB1, VCAM1, A2M, CXCL12, CEBPD, HMOX1, CD79A, TFRC, CD19, and GADD45A, with the degrees of 113, 88, 82, 70, 67, 46, 44, 42, 37, and 35, respectively. The PPI subnetwork of genes included 15 nodal genes, of which the top two were HLA-DRB1 and VCAM1, both with a degree of 4, while the degree of TFRC was 1 ([80]Figure 3). Figure 2. [81]Figure 2 [82]Open in a new tab Protein–protein interaction network of genes. Note: The red nodes represent upregulated differentially expressed genes and the green nodes represent downregulated differentially expressed genes, while the yellow nodes represent no differentially expressed genes. Figure 3. [83]Figure 3 [84]Open in a new tab Protein–protein interaction in the subnetwork of genes. Notes: The red nodes represent upregulated differentially expressed genes and the green represent downregulated differentially expressed genes, while the white indicate no differentially expressed genes. The fold change of gene expression is shown in color (deeper color indicates higher fold change of gene expression). The square nodes denote genes with lower importance in the subnetwork, and round nodes denote genes with higher importance. For the genes identified in the PPI subnetwork, a total of 20 KEGG pathways were enriched ([85]Table 3). The top three pathways were cell adhesion molecules (P=3.86E-06; in which HLA-DRB1, CD4, HLA-DMA, VCAM1, and HLA-DMB might be involved), the B-cell receptor signaling pathway (P=1.13E-05; in which CD19, CD81, CD79A, and INPP5D might participate), and antigen processing and presentation (P=1.19E-05, in which HLA-DRB1, CD4, HLA-DMA, and HLA-DMB might be involved). Table 3. The 20 enriched KEGG pathways of the protein–protein interaction subnetwork KEGG pathway Gene counts P-value Genes Cell adhesion molecules 5 3.86E-06 HLA-DRB1, CD4, HLA-DMA, VCAM1, HLA-DMB B-cell receptor signaling pathway 4 1.13E-05 CD19, CD81, CD79A, INPP5D Antigen processing and presentation 4 1.19E-05 HLA-DRB1, CD4, HLA-DMA, HLA-DMB Hematopoietic cell lineage 4 2.13E-05 CD19, HLA-DRB1, CD4, TFRC Rheumatoid arthritis 4 2.43E-05 HLA-DRB1, IL18, HLA-DMA, HLA-DMB Asthma 3 2.57E-05 HLA-DRB1, HLA-DMA, HLA-DMB Primary immunodeficiency 3 4.12E-05 CD19, CD79A, CD4 Allograft rejection 3 4.88E-05 HLA-DRB1, HLA-DMA, HLA-DMB Graft-versus-host disease 3 6.66E-05 HLA-DRB1, HLA-DMA, HLA-DMB Type 1 diabetes mellitus 3 7.70E-05 HLA-DRB1, HLA-DMA, HLA-DMB Intestinal immune network for IgA production 3 0.000107226 HLA-DRB1, HLA-DMA, HLA-DMB Malaria 3 0.000128658 CD81, IL18, VCAM1 Autoimmune thyroid disease 3 0.000136377 HLA-DRB1, HLA-DMA, HLA-DMB Staphylococcus aureus infection 3 0.000161334 HLA-DRB1, HLA-DMA, HLA-DMB Phagosome 4 0.000186619 HLA-DRB1, HLA-DMA, TFRC, HLA-DMB Viral myocarditis 3 0.000330845 HLA-DRB1, HLA-DMA, HLA-DMB Leishmaniasis 3 0.000359628 HLA-DRB1, HLA-DMA, HLA-DMB Toxoplasmosis 3 0.002107591 HLA-DRB1, HLA-DMA, HLA-DMB African trypanosomiasis 2 0.002195842 IL18, VCAM1 Systemic lupus erythematosus 3 0.002295939 HLA-DRB1, HLA-DMA, HLA-DMB [86]Open in a new tab Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes; IgA, immunoglobulin A. Discussion In this study, a total of 63 DEGs (42 downregulated, 21 upregulated) between the MM and control groups were identified. Based on functional and pathway enrichment analyses by GO and KEGG, genes such as HLA-DRB1 and VCAM1 might be involved in the positive regulation of immune system processes; meanwhile, HLA-DRB1 might also participate in the intestinal immune network for the IgA production pathway. In the PPI network, genes such as HLA-DRB1 and TFRC were ranked in the top ten nodal genes with high degrees. In the PPI subetwork, critical nodal genes, including HLA-DRB1, VCAM1, CD4, and TFRC, were identified. CD4 encodes a T-lymphocyte membrane glycoprotein, which is a receptor for human immunodeficiency virus and functions to initiate or augment the early phase of T-cell activation.[87]^37 In patients with MM, although impaired T-cell function and viral-specific cytotoxic T-lymphocytes have been reported, the presence of expanded T-cell populations such as clonal CD4^+ T-cells has also been detected.[88]^38 CD4 plays significant roles in immune responses mediated by cells and antibodies.[89]^39^–[90]^41 T-helper 17 cells (T[H]17, a new subset of CD4^+ T-cells) have been suggested as a therapeutic target in MM to improve immune function, as the number of T[H]17 cells is remarkably increased in MM, accompanied by elevated levels of several proinflammatory cytokines.[91]^38 The antigen specificity of T-cells makes great sense for the immune system.[92]^42 In the present study, CD4 was predicted as a crucial nodal gene in the PPI subnetwork and it was implicated in the antigen processing and presentation pathway, indicating that this gene might play a central role in the immune system through the involvement of this pathway in the development of MM cells from normal PCs. HLA-DRB1 (major histocompatibility complex, class II, DR Beta 1), which encodes more than 1,196 alleles, is a polymorphic locus of the HLA system.[93]^43 It is a member of the HLA class II alpha chain paralogs, which comprise an alpha (DPA) and a beta (DPB) chain and are expressed in antigen-presenting cells.[94]^44^,[95]^45 Via presenting peptides from extracellular proteins, HLA-DRB1 plays a pivotal role in the immune system.[96]^46 For instance, it could bind peptides from antigens and present them on the cell surface to be recognized by CD4 T-cells.[97]^47 Additionally, HLA-DPB1 alleles are found to be positively associated with IgA deficiency in German individuals.[98]^48 Although few studies exist regarding the immune function of HLA-DPB1 related to MM, our pathway analysis suggests that HLA-DRB1 might be involved in the intestinal immune network for the IgA production pathway, implying that HLA-DPB1 might exert a regulatory role in MM development via participation in this pathway. Notably, combined with the potential role of CD4 and the predicted correlations in the PPI subnetwork, a regulatory relationship between HLA-DPB1 and CD4 might exist. Transferrin receptor (TFRC) is closely related to cellular iron intake.[99]^49 TFRC is also related to the immune response.[100]^50 The blocking of CD4 results in an increasing peptide dose of CD22 and consequently has an intermediate effect on the response of EK35 (one T-cell clone) response to TFRC,[101]^51 indicating that a relationship might exist between the CD4 and TFRC genes. As revealed in our subnetwork, CD4 was linked to TFRC, and both were predicted to be associated with the hematopoietic cell pathway, providing a clue that TFRC might also be involved in MM progression via mediating CD4. VCAM1 is a member of the Ig superfamily.[102]^52 VLA4, which acts as a cell adhesion receptor of VCAM1, mediates the cell adhesion of myeloma cells to bone marrow stroma in MM.[103]^53 Moreover, the expression of VLA4 in MM cells could mediate VCAM1 binding to mesenchymal stem cells.[104]^4 These results hint that VCAM1 participates in the pathogenesis of MM through its receptor, VLA4. In addition, VCAM1 was implicated in the cell adhesion molecule pathway in our study. Thus, we speculate that VCAM1 might contribute to the development of MM via the cell adhesion molecule pathway. Conclusion In conclusion, several crucial genes involved in MM were identified, including CD4, HLA-DRB1, TFRC, and VCAM1. All of these genes might exert their roles in the progression of normal PCs to MM cells via immune-mediated pathways. Notably, certain regulatory correlations might exist between HLA-DRB1, CD4, and TFRC. However, because the results are based on microarray data with a small sample size, more experimental validations are warranted. Footnotes Disclosure The authors report no conflicts of interest in this work. References