Clarifying the molecular mechanism associated with carfilzomib resistance in human multiple myeloma using microarray gene expression profile and genetic interaction network

Abstract

Carfilzomib is a Food and Drug Administration-approved selective proteasome inhibitor for patients with multiple myeloma (MM). However, recent studies indicate that MM cells still develop resistance to carfilzomib, and the molecular mechanisms associated with carfilzomib resistance have not been studied in detail. In this study, to better understand its potential resistant effect and its underlying mechanisms in MM, microarray gene expression profile associated with carfilzomib-resistant KMS-11 and its parental cell line was downloaded from Gene Expression Omnibus database. Raw fluorescent signals were normalized and differently expressed genes were identified using Significance Analysis of Microarrays method. Genetic interaction network was expanded using String, a biomolecular interaction network JAVA platform. Meanwhile, molecular function, biological process and signaling pathway enrichment analysis were performed based on Gene Ontology and Kyoto Encyclopedia of Genes and Genomes. Totally, 27 upregulated and 36 downregulated genes were identified and a genetic interaction network associated with the resistant effect was expanded basing on String, which consisted of 100 nodes and 249 edges. In addition, signaling pathway enrichment analysis indicated that cytokine–cytokine receptor interaction, autophagy, ErbB signaling pathway, microRNAs in cancer and fatty acid metabolism pathways were aberrant in carfilzomib-resistant KMS-11 cells. Thus, in this study, we demonstrated that carfilzomib potentially conferred drug resistance to KMS-11 cells by cytokine–cytokine receptor interaction, autophagy, ErbB signaling pathway, microRNAs in cancer and fatty acid metabolism pathways, which may provide some potential molecular therapeutic targets for drug combination therapy against carfilzomib resistance.

Keywords:

• multiple myeloma
• carfilzomib
• drug resistance
• microarray
• interaction network
• compensate pathways

Introduction

Multiple myeloma (MM), also known as plasma cell myeloma, is an incurable cancer formed by malignant plasma cells.¹ As the second most common cancer of the blood next only to non-Hodgkin’s lymphoma, each year, over 20,000 new cases are diagnosed in the USA according to epidemiologic studies from the American Cancer Society.² Over the last 40 years, therapy with melphalan plus prednisone has been recognized as the standard of care for patients with newly diagnosed MM.³ However, older patients and patients with clinically significant coexisting illnesses may not be eligible for high-dose therapy and usually do not tolerate this treatment. For these patients, the proteasome inhibitors (bortezomib and carfilzomib) are active in relapsed or refractory myeloma, which were approved by the Food and Drug Administration for treatment of relapsed/refractory MM in 2003 and 2012, respectively.⁴ In preclinical studies, bortezomib and carfilzomib sensitized melphalan-sensitive and melphalan-resistant myeloma cell lines to melphalan by breaking down enzyme complexes and downregulated cellular responses to genotoxic stress.⁵ However, recent studies revealed that relapse of myeloma developed due to acquisition of resistance to proteasome inhibitors, owing to the mutations of proteasome complex,⁶ upregulation of transporter channels or cytochrome components⁷ and the induction of alternative compensatory pathways.⁸ Although several aspects of the mechanisms associated with acquisition of resistance to proteasome inhibitors have been studied, a systems biological perspective in terms of proteasome inhibitors resistance for MM has not been fully elucidated.

In recent years, with the rapid development of precision medicine, it is possible to analyze high-throughput screening dataset to better understand pathogenesis in terms of disease progression and drug therapeutics.^9–11 To better address this merit, herein, we identified a microarray gene expression profile originating from the carfilzomib-resistant KMS-11 versus parental human myeloma cell line to establish a comprehensive genetic interaction network in order to reveal the molecular mechanisms in carfilzomib resistance in MM, which may provide molecular information or targets for MM clinical interventions in terms of acquisition of resistance to proteasome inhibitors.

Materials and methods

Microarray dataset search strategy

Microarray dataset was downloaded from Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/) with the accession number GSE69078. In this study, Riz et al treated KMS-11 MM cell line with increasing concentrations of carfilzomib over a period of 18 weeks to establish the carfilzomib-resistant MM cell line.⁸ Total RNA was extracted from the KMS-11 cell line with or without carfilzomib treatment, and messenger RNA array was performed based on Affymetrix Human Genome U133 Plus 2.0 platform.

Differently expressed genes identification

Comparison of the gene expression profiles of carfilzomib-resistant derivatives versus parental human KMS-11 MM cell line was normalized using log₂ transformation after normalization. Significance Analysis of Microarrays (SAM, http://statweb.stanford.edu/~tibs/SAM/), a statistical technique for finding significant genes in a set of microarray experiments, was applied according to a previous publication.¹²

Genetic interaction network construction

To better understand how these significant genes identified by SAM interacted with each other, genetic interaction network was expanded using String JAVA consortium (http://string-db.org/). String, a website-based biomolecular interaction network database, has an application programming interface which enables the user to get the data without using the graphical user interface of the web page.

To better understand the potential drug-resistant mechanisms in MM, Gene Ontology consortium (GO; http://www.geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/) functional enrichment were also applied through Database for Annotation, Visualization and Integrated Discovery¹³ (https://david.ncifcrf.gov/) plug-in in String database.

Statistical analysis

For differently expressed genes identification, gene expression was considered to be significant if the threshold of false discovery rate was ≤5% and fold change was ≥2. For GO and Kyoto Encyclopedia of Genes and Genomes enrichment analysis, biological process, molecular function and signaling pathways were identified as different if the P-value was ≤5%.

Results

Sixty-three genes were found to be significantly expressed in carfilzomib-resistant KMS-11 cells

To better understand which regulators contribute to carfilzomib resistance in KMS-11 cells, differently expressed genes were screened out using SAM plug-in in Excel frame. As shown in Figure 1, after performing SAM, 63 genes were found to be differently expressed in carfilzomib-resistant KMS-11 cell line compared to its parental one, with a false discovery rate ≤5% and a fold change ≥2. Figure 2 reveals the heatmap representation of these 63 genes, which indicates that 27 genes were upregulated and 36 genes decreased dramatically. The detailed information of these genes could be found in Table 1.

Figure 1 SAM plot result output of the gene expression profiling of the microarray dataset from GSE69078.

Note: In this plot, red and green dots represent the gene sets that were up- and downregulated, respectively.
Abbreviation: SAM, Significance Analysis of Microarray.

Figure 2 Heatmap visualization of the differently expressed genes identified by SAM in carfilzomib-resistant KMS-11 (GSM1692587, GSM1692588 and GSM1692589) versus parental human myeloma cell line (GSM1692593, GSM1692594 and GSM1692595).

Note: In this picture, red represents upregulated genes, while green represents downregulated genes.
Abbreviation: SAM, Significance Analysis of Microarray.

Table 1 Significant genes identified by SAM in carfilzomib-resistant KMS-11 versus parental human myeloma cell line

Gene ID	Gene name	Fold change	Gene regulation
202201_at	BLVRB	3.652113	Up
219332_at	MICALL2	2.988681	Up
208792_s_at	CLU	2.831521	Up
208791_at	CLU	2.87382	Up
205943_at	TDO2	2.881487	Up
235343_at	VASH2	2.793228	Up
207469_s_at	PIR	2.619397	Up
205081_at	CRIP1	2.318266	Up
244407_at	CYP39A1	2.900208	Up
206140_at	LHX2	2.505247	Up
205348_s_at	DYNC1I1	2.19834	Up
211458_s_at	GABARAPL1	2.346701	Up
223464_at	OSBPL5	2.171483	Up
206435_at	B4GALNT1	2.188087	Up
226884_at	LRRN1	2.112008	Up
227307_at	TSPAN18	2.121718	Up
219740_at	VASH2	2.261368	Up
223633_s_at	BC005081	2.245276	Up
208869_s_at	GABARAPL1	2.316713	Up
203729_at	EMP3	2.063944	Up
217728_at	S100A6	2.149586	Up
232549_at	RBM11	2.141973	Up
219489_s_at	NXN	2.168232	Up
222742_s_at	IFT22	2.030835	Up
214453_s_at	IFI44	2.126185	Up
223434_at	GBP3	2.010582	Up
220432_s_at	CYP39A1	2.042243	Up
202983_at	HLTF	0.24513	Down
204273_at	EDNRB	0.261948	Down
228167_at	KLHL6	0.301736	Down
213478_at	KAZN	0.367101	Down
231202_at	ALDH1L2	0.318751	Down
209723_at	SERPINB9	0.293765	Down
204271_s_at	EDNRB	0.359747	Down
206701_x_at	EDNRB	0.320717	Down
205549_at	PCP4	0.398058	Down
210644_s_at	LAIR1	0.395248	Down
47069_at	PRR5	0.447339	Down
229830_at	Unknown	0.35393	Down
205402_x_at	PRSS2	0.44509	Down
215071_s_at	HIST1H2AC	0.455068	Down
219259_at	SEMA4A	0.412331	Down
213725_x_at	XYLT1	0.444801	Down
205016_at	TGFA	0.447057	Down
219168_s_at	PRR5	0.444994	Down
206691_s_at	PDIA2	0.466426	Down
205822_s_at	HMGCS1	0.411619	Down
219255_x_at	IL17RB	0.456308	Down
205506_at	VIL1	0.472652	Down
212816_s_at	CBS	0.459518	Down
218280_x_at	HIST2H2AA3	0.499175	Down
236451_at	LOC100996579	0.431235	Down
225502_at	DOCK8	0.456397	Down
220565_at	CCR10	0.470778	Down
228821_at	ST6GAL2	0.394775	Down
214455_at	HIST1H2BC	0.487997	Down
205463_s_at	PDGFA	0.469407	Down
205898_at	CX3CR1	0.433747	Down
209598_at	PNMA2	0.454474	Down
216470_x_at	PRSS2	0.46771	Down
224156_x_at	IL17RB	0.498466	Down
208962_s_at	FADS1	0.484027	Down
225846_at	ESRP1	0.482879	Down

Abbreviation: SAM, Significant Analysis of Microarray.

Carfilzomib-resistant genetic interaction network

To address the merit of systems biology and deepen our understanding toward how these genes regulated carfilzomib resistance in MM in a system perspective, all these significant genes were submitted to String bioinformatics platform future analysis. As shown in Figure 3, the interaction network involved in carfilzomib resistance consists of 100 nodes (genes) and 249 edges (molecular interaction), with the average node degree (the number of edges connected to the node) being 4.98. Besides, network analysis also indicated that the clustering coefficient and protein–protein interaction enrichment P-value were 0.788 and 5.41e⁻¹², respectively, which means the network has a reliable robustness.

Figure 3 Genetic interaction network associated with carfilzomib resistance in multiple myeloma based on String platform. In this picture, each circle represents a gene (node) and each connection represents a direct or indirect connection (edge).

Note: Line color indicates the type of interaction evidence and line thickness indicates the strength of data support.

GO analysis

To assess the protein–protein interaction network involved in carfilzomib resistance in the context of GO, all the nodes were submitted to Database for Annotation, Visualization and Integrated Discovery bioinformatics platform for further functional annotation. As shown in Table 2, molecular function analysis indicated that most of these genes regulated protein or enzyme binding and activities. Besides, we also evaluated the biological processes involved in this carfilzomib-resistant network (Table 3). Table 3 summarizes all the potential biological processes for carfilzomib resistance. Among them, immune response, mitopahgy/macroautophagy and cellular stress ranked as top candidates.

Table 2 Molecular function analysis of the genetic interaction network associated with carfilzomib resistance in KMS-11 cell line in terms of GO

GO ID	Molecular function	Observed gene count	FDR
GO.0003988	Acetyl-CoA C-acyltransferase activity	5	6.20E-08
GO.0005515	Protein binding	42	0.00197
GO.0046983	Protein dimerization activity	14	0.00466
GO.0042802	Identical protein binding	16	0.0063
GO.0005102	Receptor binding	17	0.0073
GO.0048407	Platelet-derived growth factor binding	3	0.0107
GO.0005161	Platelet-derived growth factor receptor binding	3	0.0124
GO.0003774	Motor activity	6	0.0164
GO.0042803	Protein homodimerization activity	11	0.0202
GO.0003985	Acetyl-CoA C-acetyltransferase activity	2	0.0243
GO.0005017	Platelet-derived growth factor-activated receptor activity	2	0.0243
GO.0003824	Catalytic activity	41	0.0257
GO.0005125	Cytokine activity	6	0.0257
GO.0016740	Transferase activity	22	0.0264
GO.0019899	Enzyme binding	17	0.0298
GO.0038085	Vascular endothelial growth factor binding	2	0.0322
GO.0004714	Transmembrane receptor protein tyrosine kinase activity	4	0.049

Abbreviations: FDR, false discovery rate; GO, Gene Ontology.

Table 3 Biological process analysis of the genetic interaction network associated with carfilzomib resistance in KMS-11 cell line in terms of GO

GO ID	Biological process	Observed gene count	FDR
GO.0009605	Response to external stimulus	38	3.22E-12
GO.0002376	Immune system process	34	1.27E-08
GO.0006955	Immune response	26	3.71E-07
GO.0009991	Response to extracellular stimulus	16	8.57E-07
GO.0009628	Response to abiotic stimulus	23	1.42E-06
GO.0060548	Negative regulation of cell death	21	1.42E-06
GO.0031667	Response to nutrient levels	15	1.79E-06
GO.0006950	Response to stress	40	2.12E-06
GO.0051716	Cellular response to stimulus	54	2.41E-06
GO.0007173	Epidermal growth factor receptor signaling pathway	11	5.14E-06
GO.0010941	Regulation of cell death	25	5.37E-06
GO.0043066	Negative regulation of apoptotic process	19	6.83E-06
GO.0000422	Mitophagy	6	1.42E-05
GO.0001934	Positive regulation of protein phosphorylation	18	1.42E-05
GO.0008284	Positive regulation of cell proliferation	18	1.42E-05
GO.0033554	Cellular response to stress	25	2.09E-05
GO.0042981	Regulation of apoptotic process	23	2.09E-05
GO.0044710	Single-organism metabolic process	42	2.09E-05
GO.0050896	Response to stimulus	56	2.38E-05
GO.0016236	Macroautophagy	7	3.06E-05
GO.0043410	Positive regulation of MAPK cascade	13	3.62E-05
GO.0016049	Cell growth	8	3.78E-05
GO.0044712	Single-organism catabolic process	18	3.78E-05
GO.0031668	Cellular response to extracellular stimulus	10	4.04E-05
GO.0030334	Regulation of cell migration	15	4.43E-05
GO.0044804	Nucleophagy	5	4.43E-05

Abbreviations: FDR, false discovery rate; GO, Gene Ontology; MAPK, mitogen-activated protein kinase.

Pathway enrichment analysis

To assess the relationship between the significantly expressed genes and carfilzomib resistance, we also evaluated the signaling pathways involved in this pathogenesis (Table 4). Notably, cytokine–cytokine receptor interaction, autophagy, ErbB signaling pathway, microRNAs in cancer and fatty acid metabolism pathways seem to confer carfilzomib resistance in human KMS-11 MM cell line.

Table 4 Signaling pathway analysis of the genetic interaction network associated with carfilzomib resistance in KMS-11 cell line in terms of GO

Pathway ID	Signaling pathway	Observed gene count	FDR
4060	Cytokine–cytokine receptor interaction	13	1.28E-07
4140	Regulation of autophagy	6	4.82E-06
280	Valine, leucine and isoleucine degradation	6	6.81E-06
1212	Fatty acid metabolism	6	8.74E-06
5215	Prostate cancer	7	1.12E-05
5214	Glioma	6	2.51E-05
71	Fatty acid degradation	5	8.59E-05
900	Terpenoid backbone biosynthesis	4	0.000108
72	Synthesis and degradation of ketone bodies	3	0.000297
270	Cysteine and methionine metabolism	4	0.000733
5200	Pathways in cancer	9	0.000733
1100	Metabolic pathways	17	0.00106
4962	Vasopressin-regulated water reabsorption	4	0.00139
5206	MicroRNAs in cancer	6	0.0015
5212	Pancreatic cancer	4	0.00433
650	Butanoate metabolism	3	0.00487
5218	Melanoma	4	0.00615
4012	ErbB signaling pathway	4	0.0111
4540	Gap junction	4	0.0111

Abbreviations: FDR, false discovery rate; GO, Gene Ontology.

Discussion

Combined with bioinformatics, high-throughput screening has become a convenient assay for drug-resistance or off-target identification.^14,15 As early as 2003, a glass-based microarray suitable for detecting multiple tetracycline (tet) resistance genes was developed and applied.¹⁶ Then, Hongisto et al developed a high-throughput three-dimensional (3D) screening method that revealed drug sensitivities between the culture models of JIMT1 breast cancer cells. Compared with the traditional method for studying cancer in vitro, the anchorage-independent three-dimensional models allowed cells to grow in two dimensions and resulted in screening out 102 compounds with multiple concentrations and biological replicates for their effects on breast cancer cell proliferation.¹⁷ Using a similar method, in the present study, we also established a genetic interaction network using the publicly available microarray dataset and the functional protein interaction platform – String. Our results revealed that cytokine–cytokine receptor interaction, autophagy, ErbB signaling pathway, microRNAs (miRNAs) in cancer and fatty acid metabolism pathways were highly associated with carfilzomib resistance in MM.

A previous study indicated that autophagy contributed to carfilzomib resistance in MM by KLF4-SQSTM1/p62, which proved our bioinformatics prediction between carfilzomib resistance and autophagy.⁸ In this study, Riz et al identified high levels of KLF4 expression often occurring in MM patients carrying the t(4;14) translocation, and acquisition of carfilzomib resistance in both t(4;14)-positive MM cell line models was associated with reduced cell proliferation, decreased plasma cell maturation and activation of prosurvival autophagy by regulation of KLF4 expression.⁸ Meanwhile, basing on the proteostasis network analysis by Acosta-Alvear et al,¹⁸ inhibition of proteasome resulted in the compensatory mechanisms through inhibition of translation and induction of autophagy, which also confirmed our prediction regarding the role of autophagy in the acquisition of resistance to carfilzomib in MM.¹⁸

miRNAs, a group of noncoding RNA molecules composed of 19–25 nucleotides, can posttranscriptionally regulate target gene expression, which results in cell development, differentiation, apoptosis and proliferation.^19,20 Besides, miRNAs are also involved in the development of drug resistance by miRNA dysregulation.²¹ By far, several labs have already focused on exploring the role of miRNAs in drug resistance using microarrays. They discovered that the epigenetic modulations of miRNAs contributed to cancer drug resistance.²² As to carfilzomib resistance, miRNA also plays a major role in regulating the fundamental cellular processes that control MM resistance to proteasome inhibitors.²³ Malek et al identified that the expression of miR29 family and Let-7A1 increased in response to bortezomib, carfilzomib and ixazomib. However, Let-7A2, Let-7D, Let-7E and Let-7F2 were downregulated in bortezomib-, carfilzomib-and ixazomib-resistant cells, compared to drug-sensitive parental cells. According to our bioinformatics analysis, MTOR, EGFR, ERBB2, PDGFA, PDGFRA and PDGFRB were involved in the subnetwork of miRNAs in cancer pathways. Since mammalian target of rapamycin (mTOR) inhibition can also induce autophagy,^24,25 previous results also support the protective role of autophagy during proteasome inhibition, indicating that mTOR inhibition may desensitize carfilzomib both through inhibition of translation and induction of autophagy by regulation by miRNAs.¹⁸

As to the ErbB signaling pathway, the relation between drug resistance and ErbB pathway has already been predicted by Azad et al.²⁶ Using the Bayesian modeling framework, potential cross-talks between epidermal growth factor receptor (EGFR)/ErbB signaling and six other signaling pathways (Notch, Wnt, G protein coupled receptor [GPCR], hedgehog, insulin receptor/insulin-like growth factor 1 receptor [IGF1R] and transforming growth factor-beta [TGF-b] receptor signaling) contributed to drug resistance in breast cancer cell lines. However, limited information regarding carfilzomib resistance in MM is available.

Besides the signaling pathways mentioned above, we also discovered many pathways like valine, leucine and isoleucine degradation,²⁷ fatty acid metabolism, fatty acid degradation,²⁸ cysteine and methionine metabolism,²⁹ and terpenoid backbone biosynthesis, which are also involved in carfilzomib resistance in MM. However, detailed information regarding the association between these pathways and carfil-zomib resistance is not available. Notably, all these pathways seem to participate in cancer energy/nutrition metabolism. Whether there are any cross-talks between cancer metabolism and MM resistance is still unknown.

Conclusion

In conclusion, using the integrated microarray gene expression profile and genetic interaction network, we explored the molecular mechanisms underlying carfilzomib resistance in MM cell line and highlighted some potential signaling pathways such as cytokine–cytokine receptor interaction, autophagy, ErbB signaling pathway, miRNAs in cancer and fatty acid metabolism pathways which may be involved in this process.

Authors contribution

All authors contributed toward data analysis, drafting and revising the paper and agree to be accountable for all aspects of the work.

Acknowledgments

We would like to thank Gene Expression Omnibus (GEO), Significance Analysis of Microarrays (SAM), Database for Annotation, Visualization and Integrated Discovery and String databases for making their data readily available to the scientific community.

This work was supported by the National Natural Science Foundation of China (No 81141052). The funders had no role in the study design, data collection and analysis, decision to publish or preparation of the manuscript.

Disclosure

The authors report no conflicts of interest in this work.