Identification of Hub Gene TIMP1 and Relative ceRNAs Regulatory Network in Colorectal Cancer

Introduction

Colorectal cancer (CRC) has become the predominant cancer worldwide with more than 1.8 million new cases diagnosed annually.^1,2 Furthermore, the five-year survival rate for patients with advanced-stage metastatic cancer is approximately 10%.¹ Like other cancers, CRC is considered to be a heterogeneous disease in which gene aberrations, cellular context, and environmental influences concur with tumor initiation, progression, and metastasis.³ Despite advances in laparoscopic and robotic surgery, more aggressive resection of metastatic disease, radiotherapy, as well as neoadjuvant and palliative chemotherapies, the new treatments had an insignificant effect on long-term survival.⁴ Thus, it is critical to make a thorough inquiry into the underlying biological mechanism of the occurrences and metastases of cancers associated with prognosis so as to discover novel biomarkers for target therapies and prognosis predictions. Although accumulating evidence has demonstrated that multiple genes and cellular pathways participate in the occurrence and development of CRC,^5,6 a paucity of knowledge regarding the potential precise molecules and potential mechanisms underlying CRC progression currently limits the ability to treat this disease.

Bioinformatics analyses, including the analysis of gene interaction networks, microarray expression profiles, and gene annotations are being utilized as powerful tools for identifying potential diagnostic and treatment-relevant biomarkers of cancers.^7,8 For example, by analyzing data from the Gene Expression Omnibus (GEO) database, Cao et al⁹ identified five genes as potential biomarkers and therapeutic targets for gastric cancer. In addition, by analyzing data from GEO and The Cancer Genome Atlas (TCGA), Zhu et al found that high expression of cyclin-dependent kinase 1 (CDK1) is a prognostic factor for hepatocellular carcinoma (HCC), making it a potential therapeutic target and biomarker for the diagnosis of HCC.¹⁰ In particular, the method of integrated bioinformatics analysis can be used to overcome inaccuracies in sequencing arising from small sample sizes. Circular RNAs (circRNAs) are a novel class of endogenous non-coding RNAs that form a covalently closed-continuous loop by back-splicing events via exon or intron circularization.¹¹ Due to the development of high-throughput sequencing, researchers have discovered that thousands of circRNAs are involved in the progression of oncogenesis, invasion, and metastasis by playing the role of “sponges” to microRNAs (miRNAs).¹² For instance, Wang et al¹³ verified that circDLGAP4 regulated lung cancer cell biological processes by sponging miRNA-143 to regulate CDK1 expression and circDLGAP4 may serve as a potential biomarker for the diagnosis and treatment of lung cancer. However, at present, most studies involving circRNAs have been limited to the sequencing of a few samples or exploring the biological function of single circRNAs. To the best of our knowledge, currently, no researchers have used integrated analysis to investigate CRC-related circRNAs.

In this study, differentially expressed mRNAs (DEmRNAs) between human CRC tissues and adjacent non-tumor tissues were identified via analysis of public TCGA datasets. Next, to explore the roles of these DEmRNAs, functional enrichment analyses and pathway enrichment analyses were performed. Then, protein–protein interaction (PPI) networks were successfully constructed. The key genes and significant modules in the networks were identified. Kaplan–Meier analysis was performed to evaluate the prognostic value of these hub genes. Furthermore, three additional circRNA expression profiles were analyzed to identify differentially expressed circRNAs (DEcircRNAs) and differentially expressed miRNAs (DEmiRNAs) between CRC and adjacent non-tumor tissues. Finally, circRNAs–miRNAs–mRNA ceRNAs network was constructed. The research is expected to help to further elucidate the ceRNAs interactions in CRC and generate insight into the potential biomarkers and targets for the diagnosis, prognosis, and therapy of CRC.

Materials and Methods

Data Collection and Preprocessing

CRC gene expression profile data were downloaded from TCGA (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) and standardized, including 41 normal samples and 473 tumor samples and their clinical information. Previous studies have demonstrated that without adjustment, TCGA-COAD READ data set could be generated by merging samples from TCGA-COAD data set and TCGA-READ data set, since principal components analyses and unsupervised hierarchical clustering showed no significant differences.^14,15 CRC miRNA expression profiles from Illumina HiSeqmiRNASeq platforms, including 8 normal samples and 450 tumor samples, were downloaded from TCGA and standardized. In addition, 4 circRNA expression profiles (GSE121895, GSE126094, GSE138589, GSE142837) from Illumina HiSeqRNASeq platforms, including 23 tumor samples and 23 normal samples, were downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo) by searching for the term “CRC” (July 2020), and batch effects were removed using the combat function in the R “sva” package.¹⁶

Identification of Differentially Expressed RNAs

DEcircRNAs, DEmRNAs, and DEmiRNAs were identified using an R package “DESeq2” (http://www.bioconductor.org/packages/release/bioc/html/DESeq2.html). |log2FC| >2 and FDR <0.05 were set as the cutoff criteria (FC, fold change; FDR, false discovery rate) based on the Benjamini‐Hochberg method for DEmRNAs and DEmiRNAs.¹⁷ DEcircRNAs were screened by |log2FC| >1 and FDR<0.05. R was used to visualize differential genes. For DEcircRNAs, we used Surrogate Variable Analysis to handle multiple GSE profiles as described above. Volcano maps were plotted based on the volcano map of R.

Functional Enrichment Analysis

To identify the biological function of the ceRNAs network, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses are widely used for gene annotation terms and pathway enrichment analysis. GO is a widely used tool for annotating genes with functions, especially molecular function (MF), biological pathways (BP), and cellular components (CC). KEGG Enrichment Analysis is a practical resource for the analytical study of gene functions and associated high-level genome functional information. ClusterProfiler package of R was performed to analyze and visualize functional profiles. A P-value < 0.05 was the threshold for significance for GO and KEGG terms.

Construction of PPI

The PPI network was conducted to analyze the functional interactions between proteins, providing insights into the mechanisms for the development of CRC. The minimum required interaction score is 0.5. The STRING website (https://string-db.org/) was employed to construct the PPI network.

Prognostic Survival Assessment of DEmRNAs

Based on the median expression level of each DEmRNAs, the CRC patients were divided into high and low-expression groups. Kaplan‐Meier analysis and the Log rank test were utilized to paint the survival curves to find the DEmRNAs that were significantly associated with the survival of CRC patients. A P-value < 0.05 was set as the threshold.

Construction of the ceRNAs Network

DEcircRNAs and DEmRNA matched by DEmiRNAs were retrieved using starBase database.¹⁸ Moreover, the prediction results of TargetScan, miRTarBase, and miRDB were integrated by starBase.^19–21 The candidates searched in three databases were associated with the most important DEmiRNAs. Finally, the circRNAs–miRNAs–mRNA ceRNAs network was constructed and visualized using R.

Quantitative Real-Time PCR (qRT-PCR) of DEmiRNAs and TIMP1 in Colon Cancer and Normal Tissue

Total RNA was prepared from colonic tissue using an RNA extraction kit (TIANGEN BIOTECH, Beijing, China), according to the manufacturer’s instructions. The extracted RNA was synthesized to form cDNA using a FastKing one-step kit (TIANGEN BIOTECH, Beijing, China). qRT-PCR was performed using a RealUniversal Color PreMix (SYBR Green) kit (TIANGEN BIOTECH, Beijing, China) to assess the expression of target genes. U6 was used as an internal control for DEmiRNAs. GAPDH was used as internal control for TIMP1. In addition, the relative expression of RNAs was quantified by using the 2^−ΔΔCt method.

Results

Identification of DEcircRNAs, DEmiRNAs, and DEmRNAs in Human CRC

Table 1 shows the clinicopathological data of 473 patients with CRC. According to the cutoff threshold, a total of 412 DEmRNAs (including 82 upregulated and 330 downregulated) were screened out between 473 CRC and 41 normal samples with the standard of logFC> 2 and an adjusted P value (adj.P.Val) <0.05 (Figure 1A). Two hundred and sixty DEcircRNAs (including 253 upregulated and 7 downregulated) were altered significantly between 23 CRC and 23 normal samples by log₂FC > 1 and an adj.P.Val < 0.05 (Figure 1B). To further establish an circRNAs–miRNAs–mRNAs ceRNAs network, we also matched DEmiRNA expression profiles in the 450 CRC and 8 normal samples. As a result, 190 DEmiRNAs reached the inclusion criteria including 82 upregulated and 108 downregulated miRNAs (Figure 1C). The top 10 DEcircRNAs, DEmiRNAs, and DEmRNAs are presented in Table 2.

Table 1 Clinicopathological Characteristics of 473 CRC Patients

Table 2 Top 10 DEcircRNAs, DEmiRNAs and DEGs in Human CRC

Figure 1 The volcano maps of DEGs between CRC samples and normal samples. (A) A total of 412 DEmRNAs including 82 upregulated and 330 downregulated genes. (B) A total of 260 DEcircRNAs including 253 upregulated and 7 downregulated genes. (C) A total of 190 DEmiRNAs including 82 upregulated and 108 downregulated genes. Red represents upregulated genes and green represents downregulated genes.

Functional Enrichment Analysis of GO and KEGG Analysis of DEmRNAs

To further analyze the functional characteristics of DEmRNAs in CRC, GO and KEGG pathway analyses were performed using ClusterProfiler package of R. DEmRNAs were functionally classified into the biological process (BP), cellular component (CC), and molecular function (MF categories). In the BP category, four of the nine most enriched terms were “regulation of protein processing”, “protein activation cascade”, “regulation of acute inflammatory response” and “complement activation”. In the CC category, the four most enriched terms were “extracellular matrix”, “collagen-containing extracellular matrix”, “blood microparticle” and “apical part of cell”. In the MF categories, the three most enriched terms were “antigen binding”, “receptor regulator activity” and “receptor ligand activity” (Figure 2A). In addition, almost 16 KEGG pathways were significantly enriched in our analysis. The three most enriched terms were “cytokine-cytokine receptor interaction”, “kineral absorption” and “steroid hormone biosynthesis” (Figure 2B).

Figure 2 Functional enrichment analysis of DEmRNAs. (A) The top 9 enrichment scores in GO enrichment analysis of the DEmRNAs including biological process enrichment analysis, cellular components enrichment analysis, molecular function enrichment analysis. (B) The top 16 enrichment scores in KEGG enrichment analysis of the DEmRNAs.

PPI Network Construction and Survival Analysis

A total of 412 DEmRNAs (82 upregulated and 330 downregulated) were used to construct the PPI networks, which included 226 nodes and 478 edges. The combined minimum required interaction score >0.5 was considered statistically significant (Figure 3). In addition, the degree distribution of each gene in the PPI network was analyzed, the top five genes [C-X-C chemokine receptor type 8 (CXCL8), TIMP1 (tissue inhibitor of metalloproteinase 1), CXCL1, secreted phosphoprotein 1 (SPP1) and CXCL12] with high connectivity were confirmed as hub genes and next were underwent survival analysis.

Figure 3 The plot of the PPI network of DEmRNAs including 226 nodes and 478 edges by the online database STRING. The combined minimum required interaction score>0.5 was considered statistically significant.

Survival Analysis of Hub Genes

The prognostic values of the five hub genes were assessed in CRC patients using Kaplan‐Meier analysis and Log rank test. The results indicated that CRC patients with high expression of TIMP1 showed worse overall survival (P=0.004). In contrast, the other four hub genes (CXCL8, CXCL1, SPP1, and CXCL12) were not related to the overall survival of CRC patients (P > 0.05) (Figure 4).

Figure 4 Kaplan-Meier survival curves for the top five hub genes including SPP1, CXCL1, TIMP1, CXCL8, and CXCL12. TIMP1 was significantly associated with survival rate of CRC patients.

ceRNAs Regulatory Network Construction

miRNAs-mRNA interactivity was taken into account, in addition to the circRNAs-miRNAs, to construct an integrated ceRNAs network. Based on the starBase database, which masters the function of transcriptome-wide mircoRNA targeting prediction, we matched 61 DEcircRNAs and 3 DEmiRNAs. To clearly show the interaction in ceRNAs, the regulatory network contained some well-described biomarkers, including, hsa-miR-671-5p, hsa-miR-17-3p, hsa-miR-328-3p and TIMP1. This ceRNAs network is particularly informative in locating potential biomarkers for CRC. For instance, hsa-miR-671-5p interacted with TIMP1 and was mediated by hsa-circ-0002191. hsa-miR-17-3p interacted with TIMP1 and was mediated by has-circ-0023397 (Figure 5).

Figure 5 The ceRNAs network of circRNAs-miRNAs-mRNA in CRC. Blue represents DEcircRNAs; Black represents DEmiRNAs; Red represents DEmRNA.

Expression Levels of DEmiRNAs and DEmRNAs in Human Colon Cancer Tissue

To identify the authenticity and feasibility of the ceRNAs regulatory network, some vital DEmiRNAs and DEmRNAs are evaluated in colon cancer tissue and normal tissues. We found that TIMP1 is highly expressed in colon cancer tissue compared to normal tissue (P < 0.001). In contrast, the expression levels of hsa-miR-671-5p, hsa-miR-17-3p, and hsa-miR-328-3p were significantly decreased in colon cancer tissue (Figure 6).

Figure 6 The expression levels of DEmRNA and DEmiRNAs in colon cancer patients compared with those of normal samples. (A) The TIMP1 is highly expressed in colon cancer tissue. (B–D) The miR-671-5p, miR-17-3p and miR-328-3p is low expression in colon cancer tissue, **P <0.01, and ***P <0.001.

Discussion

CRC, the third most commonly diagnosed malignancy and the second leading cause of cancer-related deaths with notably aggressive biological behavior and poor survival rates, has always drawn close attention from researchers.² It is crucial to identify reliable therapeutic targets and biomarkers in order to improve the clinical outcome for CRC patients. “The ceRNAs hypothesis” presents a new pattern of gene expression regulation that cicrRNAs could regulate mRNAs by competing with the corresponding miRNAs.²² Subsequently, benefits from developments in sequencing technology and the applications of bioinformatics confirm the increasingly important biological role in the initiation, progression, and metastasis of tumors.^19–21 CircRNAs differ from other long non-coding RNAs in the structure, which is characterized by covalently linked 5ʹ- and 3ʹ-ends. CircRNAs functionally act as miRNAs sponges, RNA-binding protein sponges, and gene expression regulators. Therefore, circRNAs regulate their target genes expression and proteins network in both transcriptional and post-transcriptional patterns.²³ Increasingly, clinicians consider that circRNAs-miRNAs-mRNAs ceRNAs networks could provide an integrated view of regulatory crosstalk between these CRC-specific RNA transcripts.^24,25

In this study, DEmRNAs were identified between tumor samples and normal control tissues. Then, GO and KEGG analyses were performed to further understand the role of DEmRNAs. The results of GO analyses showed that the DEmRNAs were enriched in regulation of protein processing, protein activation cascade, and acute inflammatory which is confirmed by the knowledge that protein-induced pathology and inflammatory networks underlying CRC are the main cause for tumor development and progression.^26–28 Furthermore, KEGG analyses showed that cytokine–cytokine receptor interaction is a substantial factor in the occurrence of CRC. Cytokines such as TNF-α and IL-6 are classically regarded as central players in CRC by driving activation of the NF-κB and STAT3²⁹. Cytokines including IL-11, IL-17A, and IL-22 have gained attention as facilitators of CRC.²⁹

The top degree hub genes (CXCL8, TIMP1, CXCL1, SPP1 and CXCL12) were presented in the PPI network with DEmRNAs. SPP1, also named Osteopontin, has been proven to be overexpressed in various malignant neoplasms including breast cancer, lung cancer, and gastric cancer.^30–32 Although Seo et al³³ have evaluated the expression of SPP1 in 174 stage II and III CRC specimens and found SPP1 is significantly associated with cell invasion and adherence in CRC, the underlying mechanism was not revealed. Wang et al³⁴ has shown that SPP1 functions as an enhancer of cell growth in hepatocellular carcinoma (HCC) targeted by miR-181c. Further studies have shown that SPP1 promotes the metastasis of CRC by activating epithelial-mesenchymal transition (EMT).³⁵ CXCL8, as a prototypical chemokine, is responsible for the recruitment and activation of neutrophils and granulocytes to the site of inflammation which demonstrated that CXCL8 played a crucial role in facilitating tumor growth and progression in breast cancer, prostate cancer, lung cancer, colorectal carcinoma, and melanoma.³⁶ Phosphorylation of Src-kinases and focal adhesion kinase (FAK) in cancer cells were increased in CXCL8 signaling, which contributed to cell proliferation and chemoresistance.^37,38 The level of CXCL1 are elevated in CRC and increased level of CXCL1 are associated with tumor size, advancing stage, and patient survival.^39,40 It was reported that CXCL1 could promote tumor growth by inducing angiogenesis and the recruitment of neutrophils into the tumor-associated microenvironment.^41,42 CXCL1, the most abundant secreted chemokine by tumor-associated macrophages has been implicated in the promotion of breast cancer growth and metastasis via activating NF-κB/SOX4 signaling.⁴³ The similar phenomenon has been observed in human bladder cancer.⁴⁴ Some researchers have indicated that CXCL1 could increase oncogenes expression in colon cancer, including forkhead box O1 (FOXO1) and transcription factor 4 (TCF4) in CXCL1-treated SW620 cells according to transcriptome analyses.⁴⁵ CXCL1 is also vital for pre-metastatic niche formation and metastasis in CRC.⁴⁶ CXCL12 also known as SDF-1 is widely distributed in human tissues and more than 23 different types of cancers.⁴⁷ Importantly, it has been found that CXCL4 and its ligand CXCL12 are implicated in cell proliferation, angiogenesis, migration, EMT, and tumor metastasis.⁴⁸ TIMP1 belongs to the tissue inhibitor of the metalloproteinases family which includes TIMP1, TIMP2, TIMP3, and TIMP4.⁴⁹ In the present study, TIMP1 has been reported to indicate poor prognosis in CRC (P=0.004), which is consistent with the research of Song et al.⁵⁰ Song et al considers that the expression of TIMP1 was clearly associated with the regional lymph node and distant metastasis. In addition, research by Song et al indicated that TIMP1 was an independent prognostic indicator for the progression and metastasis of colon cancer through FAK-PI3K/AKT and MAPK pathway.⁵⁰ Moreover, TIMP1 could promote receptor tyrosine kinase c-Kit signaling in CRC, while c-Kit is an important oncogene in CRC and plays a role in cell proliferation and migration.⁵¹ For other cancers, TIMP1 inhibited the chemosensitivity of breast cancer cells through the PI3K/AKT/NF-kB pathway.⁵² TIMP1 is in favor of cell survival in melanoma by activating the 3-phosphoinositide dependent kinase-1 signaling pathway.⁵³ TIMP1/CD63/ERK signaling axis induces the formation of neutrophil extracellular traps and facilitates the development of pancreatic cancer.⁵⁴ Clinical studies have demonstrated that the elevated level of TIMP1 was associated with poor prognosis in various tumors, such as breast cancer,⁵⁵ cutaneous melanoma,⁵⁶ and gastric cancer.⁵⁷ The elevated plasma level of TIMP1 predicted a reduced response to second-line hormone therapy and low survival in women with metastatic breast cancer.⁵⁸ Therefore, TIMP1 may be a potential biomarker to predict the prognosis of cancer and play a critical role as a therapeutic target. The TIMP-miRNAs axis has been believed to be a potential therapeutic target against aggressive or drug-resistant variants of human cancers.^59–61 For instance, angiogenesis and tumor growth were increased when TIMP1 banded to CD63 and stimulated miR-210 accumulation by activating PI3K–AKT–HIF1α signaling in the lung adenocarcinoma.⁶⁰ As the hub elements of the ceRNAs network, miRNAs exhibited key roles among different RNA transcripts. In fact, hsa-miR-671-5p have been proven to interact with TIMP1 directly by cross-linking immunoprecipitation.⁶² However, the interaction between hsa-miR-671-5p and TIMP1 still needs to be verified in the occurrence and progress of CRC. The miR-671-5p had a protective role in gastric cancer by targeting upregulator of cell proliferation.⁶³ Meanwhile, miR-671-5p inhibits EMT by directly downregulating FOXM1 in breast cancer.⁶⁴ Interestingly, the levels of miR-671-5p are not only increased in colon cancer tissue but also increased cell proliferation, migration, and invasion by targeting tripartite motif containing 67.⁶⁵ This finding runs against our findings. The same miRNAs can regulate multiple mRNAs molecules and produce different physiological effects. MiR-328-3p was identified in bladder cancer and suppressed cell proliferation, migration, and invasion by targeting integrin subunit alpha 5 as well as by inhibiting EMT and inactivated PI3K/AKT pathway.⁶⁶ Similar tumor suppression effects were observed in colon cancer.⁶⁷

The present study identifies a novel ceRNAs network, which implies that TIMP1 is a potential biomarker underlying the development of CRC, providing new insights into survival predictions and therapeutic targets. However, the limitation of the study still needs investigation. Too many circRNAs were chosen in this study, a more advanced approach to narrow the scope of research is needed. The results of the present CRC-related ceRNAs regulatory network are required to be verified by clinical trials and molecular experiments.

Conclusion

The present study identified a novel circRNAs-miRNAs-mRNA ceRNAs network and provided candidate prognostic biomarkers for predicting the outcome of patients with CRC. Especially, TIMP1 is a potential indicator underlying the development of CRC. This study provided new insights for the survival predictions and therapeutic targets of CRC.