This article is mentioned in:

Introduction

Colorectal cancer (CRC) is the third most common cancer type worldwide, with high morbidity and mortality rates (1). Annually, the global incidence of CRC is estimated to be ~1 million, with ~500,000 mortalities (2). Obesity, smoking, diet and a lack of exercise are risk factors associated with CRC (3). Despite advanced detection approaches, including colonoscopy and fecal immunochemical testing in early stage and precancerous lesions (4), the incidence of CRC remains high. In a previous study, in the United States in 2014, a cohort of 136,830 individuals was estimated to be diagnosed with CRC and 50,310 patient (36.8%) succumbed (5). In China, rapidly increasing incidence and mortality rates of CRC have been detected in past decades (6). Therefore, extensive studies have been conducted to investigate more effective biological therapies for CRC management. The accumulation of mutations in a large number of oncogenes and tumor suppressor genes, which could active or inhibit the pathways critical for the initiation and progression of CRC, were detected (7). Several biomarkers have been established for the detection of metastatic CRC, including KRAS and RAS mutations (8,9). Additionally, the crucial pathways were also observed. Smith et al showed that tumor protein p53 promoted the progression of CRC through the alteration of genetic pathways (10). The nuclear factor-κB signaling pathway was reported to contribute to the carcinogenesis of CRC (11). MicroRNAs (miRNAs/miRs) are small RNAs that play central roles in cancer development via the regulation of its target genes. The altered expression of miR-21, miR-31, miR-143 and miR-145 was implicated in CRC progression (12). A recent study recruiting a genome-wide screening method identified 16 vital genes in CRC, such as SCARA5, which was affected by methylation (13). However, the comprehensive regulatory mechanisms of CRC, particularly the interplayed associations between miRNAs and genes, remain obscure. The present study utilized the expression profile data in the study by Khamas et al (13) to identify the differentially-expressed genes (DEGs) between CRC tissues and paired normal control tissues. In addition, the interactions amongst the DEGs were further investigated through protein-protein interaction (PPI) network analysis. Furthermore, the miRNAs that targeted the DEGs were also predicted. As a whole, all these bioinformatical analyses were aimed to identify potential biomarkers for the prognosis and prevention of CRC, and to uncover the underlying regulatory mechanism of CRC progression.

Materials and methods

Gene expression profile data

The gene expression profile data GSE32323, which was deposited by Khamas et al (13), was used. The public Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/), was utilized in the study. The platform used was GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array; Agilent Technologies, Palo Alto, CA, USA). In the expression profile, there were 34 samples derived from the CRC patients, consisting of 17 from cancerous tissues (CRC samples) and 17 from paired normal tissues (control samples).

Identification of DEGs

Following the data preprocessing, including background correction and the transformation from probe level to gene symbol using the Affy package (14) in R language (http://www.bioconductor.org/packages/release/bioc/html/affy.html), the data was subjected to normalization with the preprocessCore package (version 1.28.0; http://www.bioconductor.org/packages/3.0/bioc/html/preprocessCore.html) (15). Subsequently, the DEGs between CRC and normal samples were selected basing on a t-test of Linear Models for Microarray Analysis package in R (version 3.22.7; http://www.bioconductor.org/packages/release/bioc/html/limma.html) (16). The fold-change (FC) of the gene expression was also calculated. The threshold criteria for the DEG selection were P<0.05 and |log2FC| ≥1.

Functional enrichment analysis of the DEGs

To investigate the functions and processes that may be altered by the identified DEGs, the Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed, using the online tool of the Database for Annotation Visualization and Integrated Discovery (version 6.7; http://david.abcc.Ncifcrf.gov/) (17), a potent program integrating the gene or protein functional annotations with graphical summary. The cut-off value for the screening of significant functions and pathways was P<0.05.

Establishment of the PPI network

The Search Tool for the Retrieval of Interacting Genes (STRING) database (version 9.1; http://string-db.org/) (18) was recruited to predict the potential interactions amongst the identified DEGs from the protein level. Only the interactions containing at least one DEG were filtered out to build the PPI network, with the criterion of a combined score of >0.4, as visualized by Cytoscape (version 3.2.1; http://cytoscape.org/) software (19).

Prediction of targets of microRNAs

Using the web-based gene set analysis toolkit (WebGestalt; Vanderbilt University, TN, USA; http://bioinfo.vanderbilt.edu/webgestalt/) (20), the regulatory miRNAs of the DEGs were selected.

Results

DEGs between CRC and normal samples

According to the aforementioned selection criteria, a set of 1,347 DEGs, including 659 upregulated genes and 688 downregulated genes, were identified.

Altered functions and pathways by the DEGs

As indicated in the results of the enrichment analysis (Table I), the upregulated DEGs were significantly enriched in biological processes (BPs) that included the mitotic cell cycle (GO:0000278), nuclear division (GO:0000280) and the cell cycle (GO:0007049), and pathways such as the cell cycle (Hsa04110) and DNA replication (Hsa03030). For the downregulated DEGs, the over-represented functional GO terms were cellular response to zinc ion (GO:0071294), cellular response to chemical stimulus (GO:0070887) and cellular response to chemical stimulus (GO:0070887), while the prominent pathways were metabolic pathways (Hsa01100) and pancreatic secretion (Hsa04972) (Table II).

Table I. GO and pathway enrichment analysis of the upregulated DEGs (top 5 in each category, as ranked by the P-value).

Table I.

GO and pathway enrichment analysis of the upregulated DEGs (top 5 in each category, as ranked by the P-value).

Category	ID	Term	Count	P-value
BP	GO:0000278	Mitotic cell cycle	102	2.63×10−24
BP	GO:0000280	Nuclear division	55	2.26×10−22
BP	GO:0007049	Cell cycle	128	3.47×10−21
BP	GO:0007067	Mitosis	55	1.25×10−18
BP	GO:0022402	Cell cycle process	117	1.15×10−18
CC	GO:0031981	Nuclear lumen	140	1.68×10−17
CC	GO:0044428	Nuclear region	158	2.32×10−16
CC	GO:0043233	Organelle lumen	164	1.44×10−15
CC	GO:0031974	Membrane-enclosed lumen	166	1.55×10−15
CC	GO:0070013	Intracellular organelle lumen	161	2.22×10−15
MF	GO:0005515	Protein binding	319	2.45×10−8
MF	GO:0005488	Binding	450	1.92×10−6
MF	GO:0003678	DNA helicase activity	9	2.10×10−5
MF	GO:0004386	Helicase activity	15	1.42×10−4
MF	GO:0008009	Chemokine activity	8	1.70×10−4
KEGG pathway	Hsa04110	Cell cycle	24	1.21×10−11
KEGG pathway	Hsa03030	DNA replication	11	3.64×10−8
KEGG pathway	Hsa03013	RNA transport	21	1.19×10−7
KEGG pathway	Hsa03008	Ribosome biogenesis in eukaryotes	15	1.60×10−7
KEGG pathway	Hsa04115	p53 signaling pathway	10	1.72×10−4

[i] GO, gene ontology; DEGs, differentially-expressed genes; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological process; CC, cell component; MF, molecular function; Count, numbers of DEGs enriched in each term.

Table II. GO and pathway enrichment analysis of the downregulated DEGs (top 5 in each category, as ranked by the P-value).

Table II.

GO and pathway enrichment analysis of the downregulated DEGs (top 5 in each category, as ranked by the P-value).

Category	ID	Term	Count	P-value
BP	GO:0071294	Cellular response to zinc ion	7	2.45×10−8
BP	GO:0070887	Cellular response to chemical stimulus	112	2.90×10−7
BP	GO:0010035	Response to inorganic substance	32	3.18×10−7
BP	GO:0006629	Lipid metabolic process	77	4.91×10−7
BP	GO:0050896	Response to stimulus	303	1.34×10−6
CC	GO:0005615	Extracellular space	69	6.65×10−11
CC	GO:0005576	Extracellular region	131	3.24×10−10
CC	GO:0044421	Extracellular region part	81	1.11×10−9
CC	GO:0071944	Cell periphery	224	1.50×10−9
CC	GO:0016020	Membrane	346	6.02×10−9
MF	GO:0019955	Cytokine binding	12	1.47×10−6
MF	GO:0097367	Carbohydrate derivative binding	20	9.75×10−6
MF	GO:0008201	Heparin binding	16	1.03×10−5
MF	GO:0005539	Glycosaminoglycan binding	18	2.74×10−5
MF	GO:0016616	Oxidoreductase activity, acting on the CH-OH group of donors, NAD or NADP as acceptor	14	3.79×10−5
KEGG pathway	Hsa01100	Metabolic pathways	69	1.21×10−4
KEGG pathway	Hsa04972	Pancreatic secretion	12	6.96×10−4
KEGG pathway	Hsa04960	Aldosterone-regulated sodium reabsorption	7	1.29×10−3
KEGG pathway	Hsa00910	Nitrogen metabolism	5	1.90×10−3
KEGG pathway	Hsa00232	Caffeine metabolism	3	2.02×10−3

[i] GO, gene ontology; DEGS, differentially-expressed genes; KEGG, kyoto encyclopedia of genes and genomes; BP, biological process; CC, cell component; MF, molecular function; Count, numbers of DEGs enriched in each term.

PPI network of the DEGs

By mapping the DEGs into the STRING database, the potential interactions of the DEGs from the protein level were predicted. As a result, a PPI network comprising 1,478 edges and 462 nodes were established. A protein in the network serves as a ‘node’, and the ‘degree’ of a node represents the number of the interactions between two nodes. Based on this definition, the top ten nodes with high degrees in the PPI network were cyclin-dependent kinase 1 (CDK1; degree=59), cyclin B1 (CCNB1; degree=48), NDC80 kinetochore complex component (degree=45), non-SMC condensin I complex, subunit G (degree=45), MAD2 mitotic arrest deficient-like 1 (MAD2L1; degree=44), centromere protein F (degree=41), BUB1 mitotic checkpoint serine/threonine kinase B (BUB1B; degree=39), centromere protein A (degree=37), PDZ-binding kinase (degree=36) and TPX2, microtubule nucleation factor (degree=36) (Fig. 1).

Figure 1. Protein-protein interaction network of the DEGs. Red nodes represent protein products of upregulated DEGs, green nodes represent protein products of downregulated DEGs and the lines between two nodes denote the interactions between them. DEGs, differentially-expressed genes.

Integrated miRNA-target regulatory network

Using the WebGestalt software, the integrated miRNA-target network was built, consisting of 459 nodes (305 miRNAs and 154 DEGs) and 646 edges (Fig. 2). In this network, the notable genes that were targeted by multiple miRNAs included SRY (sex determining region Y)-box 4 (SOX4; targeted by 27 miRs, including hsa-mir-129, hsa-mir-133a/b and hsa-mir-204), CCND1 (cyclin D1; targeted by 21 miRs, including hsa-let-7b, hsa-mir-155, hsa-mir-16 and hsa-mir-195) and v-myc avian myelocytomatosis viral oncogene homolog (MYC; targeted by 10 miRs, including hsa-mir-34a, hsa-let-7c, hsa-mir-145 and hsa-mir-24.

Figure 2. Integrated miRNA-target regulatory network. Red circle nodes represent protein products of upregulated DEGs, green circle nodes represent protein products of downregulated DEGs, pink triangular nodes represent miRNAs and the lines between two nodes denote the interactions between them. DEGs, differentially-expressed genes; miRNA/miR, microRNA.

Discussion

CRC is one of the most lethal cancers in the world (3). Biomarker therapeutic methods may be the most effective approaches for the management of CRC. In the present study, a total of 1,347 DEGs (659 upregulated and 688 downregulated) were identified between CRC and normal tissues. Among them, CDK1, CCNB1, MAD2L1 and BUB1B, which are mainly enriched in cell cycle-related BPs and pathways, were also the predominant nodes in the PPI network. The integrated miRNA-target network identified crucial genes, including SOX4 (targeted by hsa-mir-129, hsa-mir-133a/b and hsa-mir-204), MYC (targeted by hsa-mir-34a, hsa-let-7c, hsa-mir-145 and hsa-mir-24) and CCND1 (targeted by hsa-let-7b, hsa-mir-155, hsa-mir-16 and hsa-mir-195), which were all enriched in cell cycle-related pathways. CDK1, CCNB1 and CCND1 were also associated with the p53 signaling pathways.

Cell cycle-related genes that promote the proliferation of endothelial cells contribute to the progression of tumor growth and metastasis of CRC (21). CDK1 encodes for a serine/threonine kinase that controls the eukaryotic cell cycle by regulating mitotic onset, as well as the centrosome cycle (22). CDK1 promotes cell proliferation via phosphorylation and inhibition of forkhead box O1 transcription factor (23). The alteration of CDK1 has been found in numerous cancer types, including breast cancer (24), esophageal adenocarcinoma (25) and oral squamous cell carcinoma (26). Deregulated CDK1 has been found in CRC (27), and it has been demonstrated that cantharidin, the traditional Chinese medicine that could induce cell cycle arrest and apoptosis in various cancers, exerted the anticancer function via the inhibition of CDK1 activity (28).

CCNB1 is a regulatory protein involved in mitosis (29). The increased expression of CCNB1 has also been observed in multiple cancer types, including non-small cell lung cancer (29) and gastrointestinal stromal tumors (30). Moreover, CCNB1 serves as a biomarker for the prognosis of estrogen receptor-positive breast cancer (31). CCNB1 plays important roles in the cell proliferation at the G2 phase. It was previously verified that the suppression of CCNB1 by miR-93 resulted in the inhibition of tumor growth in CRC (32).

MAD2L1 and BUB1B are two major mitotic spindle checkpoints. Previous studies considered that the mutation or deficiency in checkpoint proteins may contribute to enhancing the tumor development in breast cancer (33), and the mutation of BUB1, the paralog of BUB1B, was first reported in CRC (34). However, in contrast with these findings, Yuan et al validated the overexpression of MAD2L1 and BUB1B by reverse transcription-quantitative polymerase chain reaction in breast cancer and proposed that it may alternatively be the overexpression of checkpoint genes that account for genomic instability (35).

The high expression level of SOX4, the transcription factor responsible for the regulation of embryonic development and cell control, was significantly associated with the recurrence of CRC (36). Notably, it was reported that the oncogene SOX4 was regulated by miR-129-2 in endometrial cancer, and that the overexpression of SOX4 was partly caused by the suppression of miR-129-2 (37).

MYC is a central gene that plays important regulatory roles in cell cycle progression. The deficiency of c-MYC inhibited the proliferation of tumor cells in numerous cancer types during the cell cycle through G1 into S phase (38), while the upregulation of MYC transcription by the SNP rs6983267 was demonstrated to promote the development of CRC (39). Moreover, a spectrum of studies has reported the suppression of MYC by miRNAs, including let-7a (40), miR-23a/b (41) and miR-145 (42), in various cancer types. Furthermore, the overexpression of stromal genes, such as collagen type I α2 chain (COL1A2), was also detected in CRC (43).

In the present study, the aforementioned 7 genes were upregulated in CRC samples, and the genes were all enriched in cell cycle-related BP terms and pathways, implying that these genes mediated cell cycle pathways that may play a crucial role in the tumorigenesis and progression of CRC. Combining the previous confirmations with the present predicted miRNA-target interactions, it can be speculated that SOX4 may be the target of miR-129, while MYC may be targeted by hsa-mir-145 and hsa-let-7c.

The p53 protein acts as a tumor suppressor, as it could prevent DNA damage by promoting cell cycle arrest in the G1 phase or by apoptosis. The alteration of genes in the p53 signaling pathway is tightly correlated with cancer development (44) CCND1 is a cyclin protein that functions as a regulator of CDKs, such as CDK4 or CDK6, during the cell cycle G1/S transition. Amplification of CCND1 has been observed in CRC (45) and the association between increased CCND1 and the activation of the p53 pathway has been established (46). Besides, the involvement of CDK1 and CCNB1 in the p53 signaling pathway have also been implied (47,48). The present findings indicated that CDK1, CCNB1 and CCND1 were all enriched in the p53 signaling pathway, providing a hint that the three genes may have vital roles in the progression of CRC by the regulation of the p53 signaling pathway. An extensive number of miRNAs downregulated the expression of CCND1, including miR-193b (49), miR-200b (50), miR-138b (51) and let-7b (52). Based on the correlations in the integrated miRNA-target network, CCND1 was regulated by 21 miRNAs, including hsa-let-7c, suggesting that CCND1 may be the target of hsa-let-7c.

In conclusion, the cell cycle-related pathways mediated by the CDK1, CCNB1, MAD2L1, BUB1B, SOX4, COL1A2 and MYC genes, and the p53 signaling pathway regulated by the CDK1, CCNB1 and CCND1 genes may play important roles in the progression of CRC. All these genes may be used as biomarkers for the prognosis of CRC. Furthermore, SOX4 may be targeted by miR-129 and MYC by hsa-mir-145 and hsa-let-7c, while CCND1 may be the target of hsa-let-7c. However, further experimental validation is warranted to confirm these putative regulatory correlations.

References