Identification of the anticancer effects of a novel proteasome inhibitor, ixazomib, on colorectal cancer using a combined method of microarray and bioinformatics analysis

Abstract

Purpose

The study aimed to explore the anticancer effects of a novel proteasome inhibitor, ixazomib, on colorectal cancer (CRC) using a combined method of microarray and bioinformatics analysis.

Materials and methods

Cell proliferation was tested by Cell Counting Kit-8 (CCK-8) assay for SW620 cells treated with different concentrations of ixazomib and different treatment times. The microarray analysis was conducted for six samples, including three samples of SW620 cells untreated with ixazomib and three samples of SW620 cells treated with ixazomib. The differentially expressed genes (DEGs) between untreated and treated samples were identified by the Linear Models for Microarray data (LIMMA) package in R language. The Gene Ontology (GO) functional and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed for the DEGs using the Database for Annotation, Visualization and Integrated Discovery (DAVID) and KEGG Orthology-Based Annotation System (KOBAS) online tool. The protein–protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database, and module analysis was performed for the PPI network.

Results

Ixazomib could inhibit the proliferation of SW620 cells in a dose-dependent and time-dependent manner. A total of 743 DEGs, including 203 upregulated DEGs such as HSPA6 and 540 downregulated DEGs such as APCDD1, were identified. Eighty-three GO terms were enriched for DEGs, which were mainly related to protein folding, apoptotic process, transcription factor activity, and proteasome. Thirty-seven KEGG pathways were perturbed, including pathway of apoptosis and cell cycle. Forty-six hub genes, such as TP53, JUN, and ITGA2, were screened out, and three modules with important functions were mined from the PPI network.

Conclusion

The novel proteasome inhibitor ixazomib significantly inhibited the proliferation of human CRC SW620 cells. It exerted anticancer effects through targeting the expression of DEGs, such as HSPA6, APCDD1, TP53, and JUN, and affecting the signaling pathways including apoptosis and cell cycle pathway, which demonstrated the promising potential of ixazomib for CRC therapy.

Video abstract

Point your SmartPhone at the code above. If you have a QR code reader the video abstract will appear. Or use:

http://youtu.be/r5CkNBBaRAQ

Keywords:

• colorectal cancer
• ixazomib
• differentially expressed genes
• functional enrichment analysis
• protein–protein interaction

Introduction

Colorectal cancer (CRC) is a worldwide health burden, ranking the third most common cancer in men and the second in women.¹ Despite progress made in the treatment of CRC, which has increased survival rates of CRC patients, metastasis, recurrence, and chemoresistance remain as obstacles to effective therapy.^2–4 Therefore, discovering innovative therapeutic agents are urgently required to combat these problems.

The ubiquitin–proteasome system (UPS) is the major proteolytic system mediating the degradation of intracellular proteins in eukaryotes.⁵ Misfolded, damaged, or specific proteins are targeted for destruction through the attachment of ubiquitin molecules and degraded by the 26S proteasome.^6–8 These proteasome substrates include regulatory proteins involved in cell proliferation, cell cycle regulation, transcription, signal transduction, apoptosis, and many other biological processes (BPs).^9,10 The stabilization of UPS is necessary for maintaining cellular homeostasis and normal functions.

Malignant cells, with hypermetabolism and uncontrolled proliferation, depend more on high proteasome activity than normal cells. Inhibition of proteasome activity may cause various effects on tumor cells, including anti-inflammatory, antiproliferative, and apoptotic effects, which result in cancer cell death ultimately. It has been reported that malignant cells are highly susceptible to proteasome inhibition.^11–13 Therefore, proteasome inhibitors emerged as attractive anticancer therapeutics.

Ixazomib (MLN2238), a novel proteasome inhibitor, was approved by US Food and Drug Administration in November 2015 for use in combination with lenalidomide and dexamethasone for the treatment of patients with relapsed or refractory multiple myeloma who received at least one previous therapy.¹⁴ It is the first globally approved, orally administered proteasome inhibitor used for the treatment of this patient population.¹⁵ It not only achieved satisfactory clinical outcomes in MM patients but also showed obvious anticancer effects against multiple tumors including chronic lymphocytic leukemia, prostate cancer, breast cancer, and neuroblastoma.^16–19 However, researches focusing on the antitumor effects of ixazomib on CRC are limited, and the anti-CRC mechanisms of ixazomib remain far more to be elucidated.

In this study, we utilized a microarray- and bioinformatics-based method to identify the anticancer effects of ixazomib on CRC cells. Differentially expressed genes (DEGs) were screened out between SW620 cells treated with or without ixazomib. Functional and pathway enrichment analyses were performed for the DEGs to recognize the biological functions and signaling pathways affected by ixazomib. Protein–protein interaction (PPI) network was constructed to identify hub genes. Subnetwork analysis was further performed to mine the molecular complex in the PPI network. To the best of our knowledge, this is the first study using a combined method of microarray and bioinformatics analysis to identify the comprehensive anticancer mechanisms of ixazomib on CRC cells. Our study demonstrated ixazomib as a potential drug against CRC, providing new insights on CRC therapy.

Materials and methods

Cell culture

Human CRC cell line SW620 was obtained from the Cell Bank of the Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). The cells were cultured in RPMI-1640 medium (HyClone Laboratories, Logan, UT, USA) supplemented with 10% fetal bovine serum (Clark Bioscience, Seabrook, MD, USA), 100 U/mL penicillin, and 100 μg/mL streptomycin (Beyotime Institute of Biotechnology, Shanghai, China) at 37°C in a humidified atmosphere containing 5% CO₂.

CCK-8 test for cell proliferation assay

Cell proliferation was measured by CCK-8 assay (Dojindo Molecular Technologies, Inc., Rockville, MD, USA) according to the manufacturer’s instructions. SW620 cells (2×10⁴/mL) suspended in RPMI-1640 medium containing 10% fetal bovine serum were seeded in 96-well plates with 100 μL/well and incubated for 24 h. The cells were then treated with various concentrations of ixazomib (Selleck, Shanghai, China) for 24 and 48 h. Cells incubated without ixazomib for the same periods of time were set as the control groups. After various periods of incubation, 10 μL of CCK-8 solution was added to each well and incubated with cells for 2 h at 37°C. The optical density was then detected with a microplate reader at a wavelength of 450 nm.

Cancer cells’ treatment, RNA extraction, and purification

SW620 cells were seeded in a sterile cell culture dish (10×10 cm) for 24 h first to allow for cell attachment and growth. Then, cells were treated with ixazomib at the concentration of 1 μmol/L (1,000 nmol/L), which exerts a moderate inhibition on the proliferation of SW620 cells according to the CCK-8 cell proliferation assay, for 24 h. Cells cultured without ixazomib for the same culture time were set as the controls. Following 24 h incubation at 37°C, the cells were harvested and the total RNA was extracted using the TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). RNA quantity and quality were assessed using the Nanodrop and Agilent 2100 bioanalyzer systems, respectively. The experiment was carried out in triplicate. Then, the qualified RNAs from three samples treated with ixazomib and three control (untreated) samples were further purified using the NucleoSpin^® RNA Clean-up XS kit (Macherey-Nagel, Düren, Germany) and RNase-Free DNase Set (Qiagen NV, Venlo, the Netherlands) according to the manufacturers’ instructions.

Microarray hybridization procedure

Gene expression levels of each sample were measured using the Affymetrix GeneChip Human Genome U133 Plus 2.0 Array (Affymetrix, Inc., Santa Clara, CA, USA). According to the manufacturer’s instructions, total RNA were amplified, labeled using the GeneChip™ 3′ IVT PLUS Reagent Kit (Affymetrix, Inc.) to obtain biotin-labeled cRNA. Microarray hybridization and washing were carried out using the GeneChip™ Hybridization, Wash, and Stain Kit (Affymetrix, Inc.) in Hybridization Oven 645 (Affymetrix, Inc.) and Fluidics Station 450 (Affymetrix, Inc.). Slides were scanned by GeneChip^® Scanner 3000 (Affymetrix, Inc.) and Command Console Software 4.0 (Affymetrix, Inc.) with default settings. Then, a total of six CEL files (three files of ixazomib-treated samples and three control samples) were obtained for further analysis.

Microarray data preprocessing

The raw data in CEL files of the six samples were first converted into probe expression values using the affy package in R language and preprocessed for background adjustment, quantile normalization, final summarization, and log₂ transformation by robust multiarray average (RMA) algorithm in affy package.^20,21 Then, the probe-level data were transformed to the expression values of genes according to the probe set annotation file of Affymetrix GeneChip Human Genome U133 Plus 2.0 Array, which was downloaded from the website of Affymetrix. If more than one probe corresponded to the same gene symbol, the average expression value of these probes was calculated as the expression value of the specific gene. After data preprocessing, we got the expression profile dataset including 21,275 genes for each of the six samples.

DEGs screening and hierarchical clustering analysis

The DEGs between ixazomib-treated samples and control samples were identified using the Linear Models for Microarray data (LIMMA) package in R language.²² The log₂fold change (log₂FC) was calculated. And, the |log₂FC| >1.5 and adjusted P-value <0.05 were used as the threshold for identifying DEGs. Hierarchical clustering analysis of the identified DEGs was further performed using the pheatmap package in R and displayed in a heatmap.²³

Gene Ontology (GO) functional and KEGG pathway enrichment analysis of DEGs

The GO analysis is a commonly used approach providing functional studies of large-scale genomic or transcriptomic data.²⁴ Kyoto Encyclopedia of Genes and Genomes (KEGG) database collects the pathway-related information about how molecules are networked.²⁵ In the present study, the GO functional analysis comprising three ontology categories of BP, cellular component (CC), and molecular function (MF) for the DEGs was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID, Version 6.8; https://david.ncifcrf.gov/), an online tool containing an integrated biological knowledgebase and analytic tools for systematically extracting biological meaning from a large number of genes or proteins.²⁶ The P-value <0.05 was set as the cutoff criterion for GO functional analysis. The KEGG enrichment analysis was carried out by KEGG Orthology-Based Annotation System (KOBAS; http://kobas.cbi.pku.edu.cn/) to identify the most significantly enriched KEGG pathways of the DEGs with corrected P-value <0.05 set as the threshold.²⁷

Construction of PPI network

The PPI network of DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, Version 10.0; http://www.string-db.org/), which contains information regarding >932 million interactions of >9 million proteins.²⁸ Interactions including stable physical associations, transient binding, substrate chaining, and information relay between the DEGs were considered for PPI network construction. All interactions in the PPI network were provided with a confidence score, which represented the existence possibility of the interaction judged by STRING between the two genes. Only interactions with the score of ≥0.4 were selected for PPI network construction, and the PPI network was visualized using the Cytoscape (http://cytoscape.org/) software.²⁹ Genes serve as “nodes” in the PPI network, and edges stand for associated interactions. The “degree” of a node represents the number of interactions between the specific gene and other genes. In the study, the connectivity degree of each network node was calculated by the plugin of CentiScaPe (http://apps.cytoscape.org/apps/centiscape) in Cytoscape.³⁰ Nodes with a degree of more than 10 were deemed as hub nodes, indicating close connections with other nodes and crucial roles in the network.

Module mining and analysis of the PPI network

Genes in the same network module often execute the same BPs together. In this study, the plugin of Molecular Complex Detection (MCODE) in Cytoscape was used to mine and analyze modules from original PPI network.³¹ The significant modules with MCODE score ≥5 and nodes ≥5 were further selected for GO BP analysis to find important BPs where these modules participated in, which was performed using the DAVID online tool with a threshold of P-value <0.05.

Results

Ixazomib inhibits the proliferation of SW620 cells in a dose-dependent and time-dependent manner

The effect of ixazomib on the proliferation of CRC SW620 cells was determined by the CCK-8 assay. For ixazomib at the concentration of 100, 500, 1,000, and 2,000 nmol/L with a treatment time of 24 h, the inhibition rate (mean ± standard deviation) was 4.69±8.04, 28.36±4.49, 35.03±3.49, and 39.46%±2.58%, respectively, compared with the control groups. The corresponding inhibition rate with the same ixazomib concentration gradient for a treatment time of 48 h was 13.60±8.74, 41.88±3.34, 59.29±12.70, and 84.07%±2.67%, respectively. As shown in Figure 1, the inhibition rate gradually increased with an increase in concentration and treatment time of ixazomib, indicating that ixazomib inhibits the proliferation of SW620 cells in a dose-dependent and time-dependent manner.

Figure 1 The inhibition of ixazomib on the proliferation of SW620 cells.

Identification of DEGs between ixazomib-treated and untreated samples

To remove systematic bias, the original data were first preprocessed by affy package in R language. Gene expression profile data with an excellent degree of standardization were used for the analysis of DEGs. Based on the cutoff criteria, a total of 743 DEGs were identified between ixazomib-treated and untreated samples, including 203 upregulated and 540 downregulated genes. The genes’ expression variations are shown in volcano plot in Figure 2, with DEGs marked in red. The top 10 significantly upregulated and downregulated DEGs are listed in Table 1.

Figure 2 The alterations of gene expression profile.

Abbreviation: log₂FC, log₂fold change.

Table 1 The identified top 10 upregulated and downregulated DEGs between ixazomib-treated and untreated samples

Upregulated DEGs			Downregulated DEGs
Gene name	Log2FC	P-value	Gene name	Log2FC	P-value
KRT6A	6.2852	1.12E–07	APCDD1	−4.7126	1.69E–07
HSPA6	5.3542	1.42E–07	KRT40	−4.6060	4.98E–08
GABARAPL1	4.3530	1.86E–07	LINC01003	−4.2484	1.65E–07
ADM	4.2515	4.98E–08	FAM3B	−4.2399	1.40E–07
NCF2	4.2303	4.98E–08	LGR5	−3.9729	1.40E–07
MIR22HG	4.2234	3.05E–07	INHBE	−3.9095	0.000192
CD55	4.1233	4.98E–08	ASPSCR1	−3.8382	4.98E–08
DUSP5	4.0788	1.46E–07	TSTD1	−3.7015	4.98E–08
ZFAND2A	3.9474	6.17E–08	FGF3	−3.6730	1.01E–06
CXCL8	3.8177	5.91E–06	TFF3	−3.5758	5.78E–06

Abbreviations: DEGs, differentially expressed genes; log₂FC, log₂fold change.

Hierarchical clustering of DEGs

Hierarchical clustering of the identified DEGs is displayed in Figure 3. As shown in the heatmap, three ixazomib-treated samples owned strong gene expression similarity and were categorized into one cluster, whereas three untreated samples were grouped into the other cluster. The results of hierarchical cluster analysis demonstrated that based on the expression of DEGs, the ixazomib-treated samples could be clearly distinguished from the untreated samples, which indicated good specificity of the DEGs.

Figure 3 Heatmap overview of the hierarchical clustering analysis.

GO functional enrichment analysis

Based on the cutoff criterion, a total of 83 GO terms were enriched by the DEGs, including 67 BP terms, eight CC terms, and eight MF terms. According to the BP terms, DEGs mainly participated in BPs regarding responses to unfolded protein and protein refolding (GO:0006986, GO:0036499, and GO:0042026), apoptosis (GO:0042981, GO:0043065, GO:0006915, and GO:0097193), negative regulation of cell proliferation and growth (GO:0008285 and GO:0030308), and cell cycle arrest (GO:0006977 and GO:0007050). According to the CC terms, DEGs were principally enriched in mitochondrion-related CCs (GO:0005739 and GO:0005759) and proteasome-related CCs (GO:0000502 and GO:0022624). According to the MF terms, DEGs were enriched in MFs including protein binding (GO:0005515 and GO:0051082), transcription factor activity (GO:0000982), hydrolase activity (GO:0016787), and phosphatase activity (GO:0016791). The top 35 BP terms and all the eight CC terms and MF terms are listed in Table 2, and the rest 32 BP terms are listed in Table S1.

Table 2 Top 35 predominant BP terms and all the eight CC terms and MF terms in GO functional enrichment analysis

Category	Term	Count	P-value
GO BP	GO:0006986~response to unfolded protein	9	5.85E–05
	GO:0007568~aging	16	3.79E–04
	GO:0046677~response to antibiotic	7	5.47E–04
	GO:0042981~regulation of apoptotic process	18	7.22E–04
	GO:0090200~positive regulation of release of cytochrome c from mitochondria	6	0.0020
	GO:0042149~cellular response to glucose starvation	6	0.0024
	GO:0007565~female pregnancy	10	0.0026
	GO:0042493~response to drug	21	0.0027
	GO:0090090~negative regulation of canonical Wnt signaling pathway	14	0.0030
	GO:0048661~positive regulation of smooth muscle cell proliferation	8	0.0035
	GO:0090084~negative regulation of inclusion body assembly	4	0.0036
	GO:0006950~response to stress	8	0.0038
	GO:0001666~response to hypoxia	14	0.0047
	GO:0042594~response to starvation	6	0.0055
	GO:2000483~negative regulation of interleukin-8 secretion	3	0.0062
	GO:0036499~PERK-mediated unfolded protein response	4	0.0063
	GO:1990440~positive regulation of transcription from RNA polymerase II promoter in response to endoplasmic reticulum stress	4	0.0063
	GO:0001836~release of cytochrome c from mitochondria	5	0.0063
	GO:0051090~regulation of sequence-specific DNA binding transcription factor activity	5	0.0085
	GO:0002009~morphogenesis of an epithelium	4	0.0099
	GO:0019626~short-chain fatty acid catabolic process	3	0.0102
	GO:0043065~positive regulation of apoptotic process	19	0.0109
	GO:0042026~protein refolding	4	0.0121
	GO:0006915~apoptotic process	30	0.0136
	GO:0050679~positive regulation of epithelial cell proliferation	7	0.0139
	GO:0072655~establishment of protein localization to mitochondrion	3	0.0149
	GO:0006977~DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest	7	0.0162
	GO:0097193~intrinsic apoptotic signaling pathway	5	0.0163
	GO:0044344~cellular response to fibroblast growth factor stimulus	5	0.0163
	GO:0044849~estrous cycle	4	0.0172
	GO:0032436~positive regulation of proteasomal ubiquitin-dependent protein catabolic process	7	0.0187
	GO:0042127~regulation of cell proliferation	13	0.0199
	GO:0032496~response to lipopolysaccharide	12	0.0202
	GO:0050873~brown fat cell differentiation	5	0.0203
	GO:0009408~response to heat	6	0.0205
GO-CC	GO:0005739~mitochondrion	83	7.17E–08
	GO:0070062~extracellular exosome	143	2.50E–07
	GO:0005759~mitochondrial matrix	31	7.53E–07
	GO:0005615~extracellular space	65	0.0034
	GO:0005777~peroxisome	10	0.0086
	GO:0000502~proteasome complex	7	0.0154
	GO:0022624~proteasome accessory complex	4	0.0183
	GO:0005829~cytosol	130	0.0402
GO-MF	GO:0005515~protein binding	339	4.72E–04
	GO:0000982~transcription factor activity, RNA polymerase II core promoter proximal region sequence-specific binding	5	0.0070
	GO:0016787~hydrolase activity	14	0.0138
	GO:0051087~chaperone binding	8	0.0202
	GO:0008046~axon guidance receptor activity	3	0.0216
	GO:0051082~unfolded protein binding	9	0.0340
	GO:0046703~natural killer cell lectin-like receptor binding	3	0.0355
	GO:0016791~phosphatase activity	5	0.0461

Note: Count refers to the number of DEGs enriched in the GO term.

Abbreviations: BP, biological process; CC, cellular component; DEGs, differentially expressed genes; GO, Gene Ontology; MF, molecular function; PERK, protein kinase R-like endoplasmic reticulum kinase.

KEGG pathway enrichment analysis

There were 37 KEGG pathways enriched in the study, including metabolic pathways (hsa01100), apoptosis (hsa04210), cell cycle (hsa04110), p53 signaling pathway (hsa04115), and MAPK signaling pathway (hsa04010). The top 25 significantly enriched KEGG pathways of DEGs are listed in Table 3, and the KEGG pathway of apoptosis is displayed in Figure 4, with DEGs colored in red. Notably, the pathway of proteasome (hsa03050), which was directly related to the mechanism of ixazomib, was enriched in KEGG analysis, demonstrating the credibility of the results.

Figure 4 The KEGG pathway of apoptosis.

Notes: The red-colored genes are differentially expressed genes enriched in the pathway in the study. Copyright ©2000. Reproduced from KEGG. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.²⁵
Abbreviations: ADP, adenosine diphosphate; CTLs, cytotoxic T lymphocytes; KEGG, Kyoto Encyclopedia of Genes and Genomes; MOMP, mitochondrial outer membrane permeabilization; NF-κB, nuclear factor-κB; NK, natural killer.

Table 3 The top 25 significantly enriched KEGG pathways of DEGs

ID	Term	Input number	Corrected P-value
hsa01100	Metabolic pathways	68	2.07E–12
hsa00280	Valine, leucine, and isoleucine degradation	11	9.53E–07
hsa04210	Apoptosis	13	0.0003
hsa00640	Propanoate metabolism	7	0.0003
hsa04115	p53 signaling pathway	9	0.0006
hsa00561	Glycerolipid metabolism	8	0.0011
hsa04010	MAPK signaling pathway	16	0.0012
hsa04146	Peroxisome	9	0.0012
hsa05169	Epstein–Barr virus infection	14	0.0012
hsa00052	Galactose metabolism	6	0.0012
hsa04110	Cell cycle	10	0.0038
hsa00650	Butanoate metabolism	5	0.0057
hsa04068	FoxO signaling pathway	10	0.0057
hsa04151	PI3K–Akt signaling pathway	17	0.0057
hsa01200	Carbon metabolism	9	0.0061
hsa04141	Protein processing in endoplasmic reticulum	11	0.0061
hsa04310	Wnt signaling pathway	10	0.0071
hsa00051	Fructose and mannose metabolism	5	0.0080
hsa04964	Proximal tubule bicarbonate reclamation	4	0.0174
hsa05132	Salmonella infection	7	0.0185
hsa05200	Pathways in cancer	17	0.0185
hsa03050	Proteasome	5	0.0212
hsa05205	Proteoglycans in cancer	11	0.0214
hsa00010	Glycolysis/gluconeogenesis	6	0.0214
hsa00520	Amino sugar and nucleotide sugar metabolism	5	0.0256

Note: Input number refers to the number of DEGs enriched in the pathway.

Abbreviations: DEGs, differentially expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network construction and hub gene identification

After STRING analysis of the DEGs, the PPI network was constructed with 430 nodes and 1,126 interactions and is shown in Figure 5. After degree calculating, a total of 46 hub genes were identified with the degree of >10 (Table 4). Among these, TP53, TSPO, JUN, FOS, ITGA2, CDKN1A, DECR1, HSPA5, BMP4, and HMOX1 were the top 10 hub genes with the closest connections to other nodes.

Figure 5 Protein–protein interaction network.

Notes: Green nodes stand for downregulated genes, while red nodes stand for upregulated genes. The node size is positively correlated with the degree of the gene. The nodes’ colors from green to red reflect the log₂FC values of the genes from low to high.
Abbreviation: log₂FC, log₂fold change.

Table 4 Hub genes in the PPI network

Node	Degree	Node	Degree	Node	Degree	Node	Degree
TP53	84	IL8	22	TUBB2A	15	PSMD1	12
TSPO	65	HSPA6	19	SGK1	15	SYK	12
JUN	63	MYB	19	DUSP1	14	TNFSF10	11
FOS	49	CDC25A	18	HSPA2	14	ETS1	11
ITGA2	40	DNAJB1	17	DDIT3	14	PCCA	11
CDKN1A	33	GADD45A	17	RASSF1	14	CD44	11
DECR1	31	PSMD2	16	THNSL1	14	NMU	11
HSPA5	30	HSPH1	16	PSMD12	13	SERPINH1	11
BMP4	27	HADH	16	SFN	13	CA2	11
HMOX1	26	KRT16	16	ITPR1	13	BBS10	11
ATF3	24	FOSL1	15	HDAC4	13
ISG15	24	ANXA1	15	BAG3	13

Abbreviation: PPI, protein–protein interaction.

Module analysis of the PPI network

A total of 17 modules were mined from the PPI network. Among these, three modules (Module 1, Module 2, and Module 3) with both MCODE score ≥5 and nodes ≥5 (Figure 6) were further selected for functional analysis. GO BP functional analysis showed that the functions of Module 1 were mainly associated with transcription (GO:0006366, GO:0045944, GO:1902895, and GO:0061614), the regulation of apoptosis (GO:0043065 and GO:0042981), and the regulation of cell proliferation (GO:0008285 and GO:0042127) (Table S2). Module 2 exhibited a close relationship with BPs about cell cycle arrest (GO:0071850 and GO:0007050) and cellular response to stimuli (GO:0034198 and GO:0006979) (Table S3). The genes in Module 3 were mainly associated with ubiquitin-dependent protein catabolism (GO:0051436 and GO:0043161) and cellular signaling pathway (such as GO:0090090, GO:0038061, and GO:0060071) (Table S4).

Figure 6 The significant modules in the protein–protein interaction network with MCODE score ≥5 and nodes ≥5.

Notes: (A) Module 1 with an MCODE score of 7.565 and the nodes of 24. (B) Module 2 with an MCODE score of 5 and the nodes of 9. (C) Module 3 with an MCODE score of 5 and the nodes of 5. Green nodes stand for downregulated genes, while red nodes stand for upregulated genes. The node size is positively correlated with the degree of the gene. The nodes’ colors from green to red reflect the log₂FC values of the genes from low to high.
Abbreviations: log₂FC, log₂fold change; MCODE, Molecular Complex Detection.

Discussion

CRC is one of the most lethal cancers worldwide.¹ Therapies nowadays can not effectively overwhelm the problems, such as metastasis, recurrence, and chemoresistance.^2–4 Searching for novel strategies is still of great urgency for CRC treatment. Ixazomib, a second-generation oral proteasome inhibitor, has shown excellent therapeutic effects on multiple myeloma patients. Meanwhile, it also displayed obvious anticancer effects against some solid tumors in vitro and in vivo.^17–19 Therefore, it is valuable to explore the antitumor effect of ixazomib on CRC.

In the present study, we first investigated the effect of ixazomib on the proliferation of SW620 cells. It was observed that ixazomib had antiproliferation effect on SW620 cells. CCK-8 assay showed that ixazomib significantly reduced cell proliferation in a dose-dependent and time-dependent manner. To further explore the anticancer activities and the potential mechanisms of ixazomib on CRC cells, we performed a combination analysis of microarray and bioinformatics.

According to the microarray analysis, a total of 743 DEGs were identified between ixazomib-treated and untreated samples. Among these, 203 genes, such as HSPA6, CASP7, and CDKN1A, were upregulated and 540 genes, such as APCDD1, MYB, and TGFBI, were downregulated. Up to 83 GO terms were enriched by the DEGs, including BP terms regarding apoptosis, regulation of proliferation, and cell cycle arrest, MF terms regarding protein binding and transcription factor activity, and CC terms regarding mitochondrion and proteasome. There were 37 KEGG pathways significantly affected by ixazomib, such as pathway of apoptosis and pathway of cell cycle. In addition, 46 hub genes including TP53, JUN, and ITGA2 and three modules with important functions were mined from the PPI network. These results demonstrated that these altered genes, biological functions, and pathways may be important targets of ixazomib acting against CRC.

HSPA6, heat shock 70 kDa protein 6 (also known as HSP70B′), is a member of heat shock protein 70 (HSP70) family. Although the function of HSPA6 remains far more to be elucidated due to limited researches, a recent study has demonstrated that HSPA6 augmented garlic extract-induced inhibition of proliferation, migration, and invasion of bladder cancer cells.³² In that study, HSPA6 was the most highly upregulated gene in bladder cancer EJ cells after garlic extract treatment.³² Overexpression of HSPA6 has been verified to enhance the G2/M phase cell cycle arrest and the inhibition of matrix metallopeptidase 9 enzyme activity induced by garlic extract.³² In the present study, HSPA6 was the second most highly upregulated gene as a result of ixazomib treatment (log₂FC =5.3542, adjusted P-value =1.42E−07), suggesting that HSPA6 may be a key target of ixazomib in preventing tumor proliferation and invasion.

APCDD1 was the most downregulated gene in response to ixazomib in SW620 cells (log₂FC =−4.7126, adjusted P-value =1.69E−07). APCDD1 has been proposed to play an important role in the proliferation of cancer cells.³³ Plasmid-mediated overexpression of APCDD1 has been shown to promote the growth of colon cancer cells, while reducing the expression of APCDD1 by transfecting with antisense S-oligodeoxynucleotides was declared to decrease tumor growth.³³ In the study, APCDD1 was significantly reduced after ixazomib treatment, demonstrating a potential antiproliferation mechanism of ixazomib through downregulating the expression of APCDD1.

An important identified anti-CRC effect of ixazomib is the induction of apoptosis, reflecting in the enrichment of apoptotic GO terms such as positive regulation of apoptotic process (GO:0043065), intrinsic apoptotic signaling pathway (GO:0097193), and the activation of apoptosis KEGG pathway, which is the third significant KEGG pathway enriched by the DEGs.

Induction of apoptosis is an important mechanism of anti-cancer agents. There are two mainly independent apoptotic signaling cascades: the extrinsic and intrinsic pathways. The former is usually initiated by the binding of extracellular death receptor ligands, such as tumor necrosis factor (TNF), TNF-related apoptosis-inducing ligand (TRAIL), and Fas ligand (FasL) to their transmenbrane receptors, which results in the activation of caspase-8 and caspase-10. The latter, in contrast, is often triggered by internal sensors, such as DNA damage, metabolic stress, and the presence of unfolded proteins, which leads to the release of mitochondrial cytochrome c into the cytosol and the subsequent activation of caspase-9. Caspase-8, caspase-9, and caspase-10 are called initiator caspases. Once activated, they can spark the downstream executor caspases, such as caspase-3, caspase-6, and caspase-7 (CASP7), which cleave a series of substrates, leading to characteristic apoptotic signs including nuclear condensation and DNA fragmentation.³⁴

In the present study, the important executor caspase, CASP7, was significantly upregulated in ixazomib-treated samples (log₂FC =2.1502, adjusted P-value =3.05E−07). Meanwhile, apoptotic BPs such as intrinsic apoptotic signaling pathway (GO:0097193) and positive regulation of release of cytochrome c from mitochondria (GO:0090200) together with apoptosis-related CCs such as mitochondrion (GO:0005739) and mitochondrial matrix (GO:0005759) were enriched by DEGs, demonstrating that ixazomib may execute its anti-CRC effects by upregulating apoptotic proteins and inducing cytochrome c releasing from mitochondria into cytosol, which resulted in the apoptosis of SW620 cells eventually.

Another important anticancer mechanism of ixazomib is the perturbation of cell cycle progression, reflecting in GO cell cycle arrest terms (such as GO:0007050~cell cycle arrest and GO:0006977~DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest) and the KEGG cell cycle pathway (hsa04110).

CDKN1A, also known as p21, was enriched in all the cell cycle-related GO terms and the KEGG cell cycle pathway in the study. The most known function of p21 is serving as a negative cell cycle regulator. It can block the cell cycle progression through binding to and inhibiting cell cycle promoters such as cyclin-dependent kinases (CDKs) and proliferating cell nuclear antigen (PCNA).³⁵ Overexpression of p21 has been shown to induce G1 phase arrest and inhibit tumor growth in vitro and in vivo.³⁶ It is also essential to sustain G2 arrest after DNA damage.³⁷

In the study, p21 (CDKN1A) was significantly upregulated after ixazomib treatment (log₂FC =2.7978, adjusted P-value =1.01E−05), together with cell cycle-related GO terms and cell cycle KEGG pathway, suggesting a possible anticancer mechanism of cell cycle arrest induction by ixazomib on CRC cells.

Tumor protein p53 (TP53), a hub gene with the highest degree in the PPI network (degree =84), was downregulated in response to ixazomib treatment (log₂FC =−2.0571, adjusted P-value =1.96E−06). According to the International Agency for Research on Cancer (IARC) TP53 database (version R18, http://p53.iarc.fr/), SW620 is a cell line bearing mutant p53.³⁸ Compared to wild-type p53, mutations in TP53 not only result in the loss of tumor suppressive functions but also lead to oncogenic properties, which is called gain oncogenic function (GOF) of the gene.^39,40 The GOF of mutant p53 results in tumor progression and resistance to therapy.⁴⁰ Therefore, targeting mutant p53 is an attractive strategy of cancer therapy. In the present study, the expression of TP53 was decreased, indicating that the oncogenic mutation p53 could be inhibited by ixazomib.

JUN (also known as c-Jun), is a hub gene with the third highest degree (63) in the PPI network. Despite the oncogenic role of JUN reported by some researches,^41,42 recent studies have revealed the tumor suppressive function of the gene.^43–45 It has been shown that the compound TRX-E-002-1 could induce c-Jun-dependent apoptosis in ovarian cancer stem cells through activating and enhancing the expression of c-Jun.⁴³ Prolonged activation of c-Jun promotes apoptosis induced by α-tocopherol ether-linked acetic acid in human prostate cancer cells.⁴⁴ Moreover, proteasome inhibitors, including bortezomib, carfilzomib, and ixazomib, were demonstrated to trigger multiple myeloma cell death via upregulation of c-Jun.⁴⁵

Consistent with previous researches, c-Jun was significantly upregulated in ixazomib-treated samples in the study (log₂FC =2.5028, adjusted P-value =2.05E−06). c-Jun was also involved in GO BPs such as response to drug (GO:0042493), release of cytochrome c from mitochondria (GO:0001836), and negative regulation of cell proliferation (GO:0008285). Combining these results, ixazomib was indicated to inhibit the proliferation and induce the apoptosis of SW620 cells through upregulating the expression of c-Jun.

ITGA2, a hub gene in Module 1, has a connectivity degree of 40, which is the fifth highest degree in the PPI network. Integrin α₂, the protein encoded by ITGA2, is an important subunit of integrin α₂β₁, which is a receptor for collagens.⁴⁶ Integrin α₂β₁ has been demonstrated to suppress metastasis of breast cancer.⁴⁷ Loss of expression of integrin α₂ subunit resulted in more advanced disease and reduced the survival of CRC patients.⁴⁸ These researches suggested a tumor-suppressing function of ITGA2. At present, the expression of ITGA2 was increased in response to the treatment of ixazomib (log₂FC =1.5901, adjusted P-value =1.48E−04), indicating a potential ITGA2-mediated anticancer effect enhanced by ixazomib.

Activating transcription factor 3 (ATF3), a member of the activation transcription factor/cAMP responsive element-binding family of transcription factors, is a hub gene in Module 2 with a degree of 24. ATF3 has significant effects on the cell cycle progression and apoptosis during carcinogenesis.^49–52 It is reported to induce S phase arrest and apoptosis in diarylheptanoid-induced neuroblastoma cell death.⁵³ It also acts as a proapoptotic protein in response to retigeric acid B-induced apoptosis in prostate cancer cells through activating downstream genes involved in cell cycle arrest and cell death, such as DDIT4, GADD45A, and CDC25A.⁵⁴ Moreover, ATF3 plays a driving role in various drugs induced apoptosis in CRC.^55–57 According to the microarray analysis, ATF3 was significantly upregulated as a result of ixazomib treatment (log₂FC =2.0248, adjusted P-value =2.32E−05), indicating that ixazomib exerts an antitumor function by elevating the expression of ATF3, inducing cell cycle arrest, and apoptosis in CRC.

The functions of Module 3 were mainly related with ubiquitin-dependent protein catabolism and proteasome. The genes of PSMD1, PSMD2, PSMC4, and PSMD12 in Module 3 encode separate subunits of proteasome (gene database at The National Center for Biotechnology Information (NCBI), https://www.ncbi.nlm.nih.gov/gene). They all were found to be upregulated after ixazomib treatment (log₂FC =1.5276, 2.2345, 2.0028, and 1.7703, respectively, and adjusted P-value =1.42E−05, 7.15E−07, 1.30E−06, and 7.24E−06, respectively). Since ixazomib potently inhibits the function of proteasome, upregulation of these genes were considered to be a compensate feedback of CRC cells trying to remedy the decreased performance of proteasome.

Conclusion

Our study provides a comprehensive microarray and bioinformatics analysis of the anticancer effects of ixazomib on CRC. Ixazomib was demonstrated to significantly inhibit the proliferation of SW620 cells in a dose-dependent and time-dependent manner. DEGs, such as HSPA6, APCDD1, and p21, hub genes, including TP53, JUN, and ITGA2, and pathways of apoptosis and cell cycle may be the therapeutic targets of ixazomib against CRC. Our results indicated the promising potential of ixazomib for CRC therapy, and further experimental studies are needed to confirm our findings.

Supplementary materials

Table S1 The rest 32 BP terms in GO functional enrichment analysis

GO BP term	Count	P-value
GO:0007165~signal transduction	52	0.0223
GO:0055114~oxidation-reduction process	30	0.0226
GO:0009116~nucleoside metabolic process	4	0.0234
GO:0001822~kidney development	8	0.0234
GO:0071277~cellular response to calcium ion	6	0.0259
GO:0071493~cellular response to UV-B	3	0.0267
GO:0030178~negative regulation of Wnt signaling pathway	6	0.0279
GO:0071356~cellular response to tumor necrosis factor	9	0.0291
GO:0005975~carbohydrate metabolic process	12	0.0297
GO:0043407~negative regulation of MAP kinase activity	5	0.0300
GO:0071243~cellular response to arsenic-containing substance	3	0.0335
GO:0010165~response to X-ray	4	0.0345
GO:0000079~regulation of cyclin-dependent protein serine/threonine kinase activity	5	0.0388
GO:0006869~lipid transport	7	0.0393
GO:0008285~negative regulation of cell proliferation	21	0.0394
GO:0071456~cellular response to hypoxia	8	0.0394
GO:0003334~keratinocyte development	3	0.0410
GO:0070365~hepatocyte differentiation	3	0.0410
GO:0006105~succinate metabolic process	3	0.0410
GO:0031953~negative regulation of protein autophosphorylation	3	0.0410
GO:0006983~ER overload response	3	0.0410
GO:0045071~negative regulation of viral genome replication	5	0.0420
GO:0030819~positive regulation of cAMP biosynthetic process	5	0.0420
GO:0006974~cellular response to DNA damage stimulus	13	0.0434
GO:0007050~cell cycle arrest	10	0.0441
GO:0006954~inflammatory response	20	0.0470
GO:0030308~negative regulation of cell growth	9	0.0471
GO:0043433~negative regulation of sequence-specific DNA binding transcription factor activity	6	0.0475
GO:0030433~ER-associated ubiquitin-dependent protein catabolic process	6	0.0475
GO:0016266~O-glycan processing	6	0.0475
GO:0006935~chemotaxis	9	0.0490
GO:0006081~cellular aldehyde metabolic process	3	0.0491

Note: Count refers to the number of DEGs enriched in the GO term.

Abbreviations: BP, biological process; DEGs, differentially expressed genes; ER, endoplasmic reticulum; GO, Gene Ontology; MAP, mitogen-activated protein.

Table S2 The GO BP functional analysis of Module 1

GO BP term	Count	P-value
GO:0007565~female pregnancy	5	5.01E–06
GO:0042493~response to drug	6	3.85E–05
GO:0006366~transcription from RNA polymerase II promoter	7	3.88E–05
GO:0045944~positive regulation of transcription from RNA polymerase II promoter	8	1.78E–04
GO:1902895~positive regulation of pri-miRNA transcription from RNA polymerase II promoter	3	3.07E–04
GO:0045893~positive regulation of transcription, DNA-templated	6	4.55E–04
GO:0051090~regulation of sequence-specific DNA binding transcription factor activity	3	4.83E–04
GO:0043065~positive regulation of apoptotic process	5	5.67E–04
GO:0032570~response to progesterone	3	0.0012
GO:0007568~aging	4	0.0013
GO:0006954~inflammatory response	5	0.0014
GO:0008285~negative regulation of cell proliferation	5	0.0016
GO:0051591~response to cAMP	3	0.0016
GO:0042127~regulation of cell proliferation	4	0.0017
GO:0034097~response to cytokine	3	0.0021
GO:0061614~pri-miRNA transcription from RNA polymerase II promoter	2	0.0026
GO:0009612~response to mechanical stimulus	3	0.0027
GO:0048661~positive regulation of smooth muscle cell proliferation	3	0.0028
GO:0006950~response to stress	3	0.0029
GO:0060395~SMAD protein signal transduction	3	0.0030
GO:2000144~positive regulation of DNA-templated transcription, initiation	2	0.0091
GO:0009629~response to gravity	2	0.0104
GO:0008283~cell proliferation	4	0.0116
GO:0000733~DNA strand renaturation	2	0.0117
GO:0070365~hepatocyte differentiation	2	0.0130
GO:0032496~response to lipopolysaccharide	3	0.0193
GO:0035994~response to muscle stretch	2	0.0208
GO:0006940~regulation of smooth muscle contraction	2	0.0208
GO:0044849~estrous cycle	2	0.0221
GO:0030224~monocyte differentiation	2	0.0221
GO:0038095~Fc-epsilon receptor signaling pathway	3	0.0225
GO:0031668~cellular response to extracellular stimulus	2	0.0233
GO:0051412~response to corticosterone	2	0.0233
GO:0001701~in utero embryonic development	3	0.0246
GO:0045892~negative regulation of transcription, DNA-templated	4	0.0264
GO:0045987~positive regulation of smooth muscle contraction	2	0.0272
GO:0007186~G-protein coupled receptor signaling pathway	5	0.0276
GO:0010165~response to X-ray	2	0.0284
GO:0032967~positive regulation of collagen biosynthetic process	2	0.0297
GO:0045648~positive regulation of erythrocyte differentiation	2	0.0310
GO:0042981~regulation of apoptotic process	3	0.0313
GO:0021983~pituitary gland development	2	0.0361
GO:0006915~apoptotic process	4	0.0366
GO:2000379~positive regulation of reactive oxygen species metabolic process	2	0.0386
GO:0046677~response to antibiotic	2	0.0411
GO:0070555~response to interleukin-1	2	0.0424
GO:0010634~positive regulation of epithelial cell migration	2	0.0424
GO:0030501~positive regulation of bone mineralization	2	0.0449
GO:0045786~negative regulation of cell cycle	2	0.0474

Note: Count refers to the number of DEGs enriched in the GO BP term.

Abbreviations: BP, biological process; DEGs, differentially expressed genes; GO, Gene Ontology.

Table S3 The GO BP functional analysis of Module 2

GO BP term	Count	P-value
GO:0036499~PERK-mediated unfolded protein response	3	9.81E–06
GO:1990440~positive regulation of transcription from RNA polymerase II promoter in response to endoplasmic reticulum stress	3	9.81E–06
GO:0071850~mitotic cell cycle arrest	3	1.16E–05
GO:0007050~cell cycle arrest	3	0.0014
GO:0036500~ATF6-mediated unfolded protein response	2	0.0037
GO:0006983~ER overload response	2	0.0042
GO:0043065~positive regulation of apoptotic process	3	0.0063
GO:0034198~cellular response to amino acid starvation	2	0.0104
GO:2000379~positive regulation of reactive oxygen species metabolic process	2	0.0124
GO:0071479~cellular response to ionizing radiation	2	0.0129
GO:0043066~negative regulation of apoptotic process	3	0.0141
GO:0034605~cellular response to heat	2	0.0153
GO:0000079~regulation of cyclin-dependent protein serine/threonine kinase activity	2	0.0161
GO:0042542~response to hydrogen peroxide	2	0.0211
GO:0070373~negative regulation of ERK1 and ERK2 cascade	2	0.0239
GO:0043433~negative regulation of sequence-specific DNA binding transcription factor activity	2	0.0247
GO:0006977~DNA damage response, signal transduction by p53 class mediator resulting in cell cycle arrest	2	0.0256
GO:0006979~response to oxidative stress	2	0.0450

Note: Count refers to the number of DEGs enriched in the GO BP term.

Abbreviations: ATF6, activating transcription factor 6; BP, biological process; DEGs, differentially expressed genes; ER, endoplasmic reticulum; ERK1, extracellular signal-regulated kinases 1; ERK2, extracellular signal-regulated kinases 2; GO, Gene Ontology; PERK, protein kinase R-like endoplasmic reticulum kinase.

Table S4 The GO BP functional analysis of Module 3

GO BP term	Count	P-value
GO:0090090~negative regulation of canonical Wnt signaling pathway	5	8.56E–09
GO:0006521~regulation of cellular amino acid metabolic process	4	1.05E–07
GO:0002479~antigen processing and presentation of exogenous peptide antigen via MHC class I, TAP-dependent	4	2.01E–07
GO:0038061~NIK/NF-kappaB signaling	4	2.31E–07
GO:0051436~negative regulation of ubiquitin-protein ligase activity involved in mitotic cell cycle	4	2.89E–07
GO:0051437~positive regulation of ubiquitin-protein ligase activity involved in regulation of mitotic cell cycle transition	4	3.55E–07
GO:0031145~anaphase-promoting complex-dependent catabolic process	4	4.00E–07
GO:0060071~Wnt signaling pathway, planar cell polarity pathway	4	6.34E–07
GO:0043488~regulation of mRNA stability	4	8.93E–07
GO:0002223~stimulatory C-type lectin receptor signaling pathway	4	9.46E–07
GO:0033209~tumor necrosis factor-mediated signaling pathway	4	1.35E–06
GO:0090263~positive regulation of canonical Wnt signaling pathway	4	1.42E–06
GO:0050852~T-cell receptor signaling pathway	4	2.67E–06
GO:0038095~Fc-epsilon receptor signaling pathway	4	4.65E–06
GO:0000209~protein polyubiquitination	4	5.14E–06
GO:0043161~proteasome-mediated ubiquitin-dependent protein catabolic process	4	6.90E–06
GO:0000165~MAPK cascade	4	1.48E–05

Note: Count refers to the number of DEGs enriched in the GO BP term.

Abbreviations: BP, biological process; DEGs, differentially expressed genes; GO, Gene Ontology; MAPK, mitogen-activated protein kinase; MHC, major histocompatibility complex; NF-kappaB, nuclear factor-kappaB; NIK, nuclear factor-kappaB-inducing kinase; TAP, transporter associated with antigen processing.

Disclosure

The authors report no conflicts of interest in this work.