A Novel Defined Necroptosis-Related miRNAs Signature for Predicting the Prognosis of Colon Cancer.

OBJECTIVE: This study aims at exploring the relationship between necroptosis-related miRNAs and colon cancer prognosis.
METHODS: We downloaded the miRNA sequencing data from the TCGA, and eight differentially expressed necroptosis-related miRNAs were screened. Then, we used Cox regression analysis to establish a prediction model of necroptosis-related miRNA. Finally, the prognosis related miRNAs were used to predict the target genes, and functional analysis was used to explore the potential mechanism of these target genes.
RESULTS: The miRNA-seq data of 444 COAD cases were downloaded from TCGA. We identified 8 differentially expressed miRNAs (has-miR-16-5p, has-miR-141-3p, has-miR-148a-3p, has-miR-425-5p, has-miR-7-5p, has-miR-223-3p, has-miR-200a-5p, and has-miR-500a-3p), then Cox analysis was performed for determining eight-miRNA signature prognostic biomarkers with obviously different OS. The area under the curve (AUC) of receiver operating characteristic (ROC) curve for predicting 1-, 3-, and 5-year survival were 0.663, 0.653 and 0.639, respectively. The multivariate analysis also implied that the risk score was an independent prognostic factor considering other confounding factors (HR = 1.847, 95% CI = 1.197-2.848, P = 0.006). According to the Kaplan-Meier analysis, the expression of hsa-miR-500a-3p (P = 0.003), hsa-miR-16-5p (P = 0.004) and hsa-miR-148a-3p (P = 0.035) significantly affected OS outcomes. We predicted the target genes of these three miRNAs and then screened 10 hub genes (CCND1, SMAD3, SMAD2, CDK1, TGFB2, CDC25A, CHEK1, VEGFA, CCNE1, WEE1). In addition, CHEK1 was associated with the survival prognosis.
CONCLUSION: Our study demonstrated that necroptosis is closely associated with colon cancer, and the model of eight necroptosis-related miRNAs are potentially useful prognostic biomarkers and therapeutic targets for colon cancer.

Introduction

Colon cancer is one of the most common cancers in the world and one of the main causes of death around the world. In 2020, there were than 1.1 million new incidents of colon cancer worldwide (6% of all cancer cases), and the death toll was about 0.58 million (5.8% of all cancer-related deaths).1 The poor prognosis is due to the high recurrence rate and early metastasis of colon cancer,2 which is closely related to the limitations of surgery, the lack of specificity of chemotherapy and radiotherapy, and the drug resistance of antitumor drugs.3 However, attributable to the enormous differences in prognosis at the same stage, the most commonly used approach, the tumor-node-metastasis (TNM) classification system for predicting survival, can only provide limited information.4 Therefore, there is an urgent need to develop reliable markers with sound recognition and survival prediction function for colon cancer.

The function of microRNAs (miRNAs) is similar to that of small interfering RNAs (siRNAs). They regulate the expression of tumor suppressor genes or oncogenes at the post-transcriptional level and play a vital role in the occurrence and progression of various cancers, including colon cancer. Given its great potential in cancer screening, diagnosis and prognosis prediction, miRNA is often used as a molecular marker to predict the prognosis of cancer patients.5 Today, thanks to the development of high-throughput microarray technology, we can detect differentially expressed miRNAs between cancer tissues and normal tissues on a large scale.3,6 The network-based screening method can help us identify treatment targets more accurately and efficiently and make up for the differences brought by different individuals. However, when we identify too many target genes, a very complex protein-protein interaction (PPI) network will bring difficulties to our further research, not to mention many false-positive targets. At this time, a method for analyzing and screening the PPI network provided by the systems biology method is fundamental, which principle is that miRNA preferentially regulates proteins that interact more than the average number of interactions with other proteins in the network.7

Programmed cell death includes apoptosis and non-apoptotic death, just as necroptosis, pyroptosis, and autophagy.8 Necroptosis is a recently identified death mode irrelevant to caspases,9 associated with various human diseases, including inflammatory disorders, ischemic reperfusion injury, degenerative diseases, infectious, autoimmune diseases, and cancer.10,11 On the one hand, its critical regulatory gene can promote the progression and metastasis of cancer; on the other hand, it can eliminate apoptosis tolerant cancer cells, which is why necroptosis has a dual effect on cancer.12 The extravasation and metastasis of tumor cells are caused by endothelial cell necrosis caused by tumor cells, which requires the involvement of amyloid precursor protein and the death receptor 6.13 According to some reports, ever-more regulators have shown anticancer effects by targeting TNF, RIPK1, RIPK3 and MLKL, closely related to necroptosis.14 In recent years, studies have reported that some miRNAs regulated necroptosis in various diseases, and these miRNAs also played an essential role in cancer progression and metastasis, such as miR-143-3p can inhibit RIPK1 mediated inflammation and necroptosis,15 while its mimics might inhibit the progression of bone metastasis of prostate cancer.16

However, no studies have reported using necroptosis-related miRNAs to predict the prognosis of colon cancer patients. In other words, the relationship between them is not clear, which is very necessary to explore molecular markers for predicting the prognosis of colon cancer patients by using necroptosis-related miRNAs. Therefore, we analyzed the data of necroptosis-related miRNAs differentially expressed in normal and colon cancer tissues to explore these miRNAs’ prognostic value and study the role of their target genes.

Materials and Methods

Dataset

The miRNA sequencing data of 444 TCGA-COAD patients were downloaded from the TCGA database on 28 September 2021 (). Meanwhile, the corresponding clinical information composed of 466 files (457 cancer samples and 8 non-cancer samples) was acquired. In addition, the mRNA sequencing data of 456 TCGA-COAD patients for the predictive analysis of hub genes were also obtained from the TCGA database.

Identification of Differentially Expressed Necroptosis-Related miRNAs

We obtained 13 necroptosis-related miRNAs from a previous review,12 presented in . Differentially expressed miRNAs were identified using R language 4.1.1 version “limma” package with a P-value < 0.05. The expression matrix of differentially expressed miRNAs was displayed by heatmap.

Development and Evaluation of the Necroptosis-Related miRNAs Prognostic Model

Patients were chosen whose clinical data included the number of survival days, survival status, age, gender, and TMN classification. Using the Cox regression analysis, it was evaluated whether the miRNAs associated with necroptosis had predictive value in the TCGA cohort. All differentially expressed miRNAs obtained from TCGA were included in the multivariable Cox regression model, and the results were displayed using the “forest plot” package. After that, each miRNA’s product and coefficient were added and calculated into the risk model.

Evaluation of the Risk Score

TCGA-COAD samples were stratified into two subgroups on the basis of their median risk scores, and their overall survival (OS) times were then compared using Kaplan-Meier method. The 1-, 3- and 5-year ROC curves were completed using “survival”, “ROCR”, and “timeROC” packages. After extracting the clinical parameters of the TCGA cohort, such as age, gender and TMN stage (I or II vs III or IV), we analyzed them with the risk score via Cox regression models.

Prediction of Target Genes and Functional Analysis

The target genes of miRNAs associated with OS were predicted by miRDB, miRTarBase and TargetScan. We used R language version 4.1.1 to identify the target genes in all three databases in the meantime, and we also used Cytoscape 3.8.0 to chart the relationship between them. We applied the “clusterProfiler” package to these overlapping genes in order to perform GO enrichment analysis and KEGG analysis.17

Sifting of Hub Genes

The STRING database was used to construct a PPI interaction network based on miRNA targets. Accordingly, the network interrelation file was imported into Cytoscape 3.8.0, and the hub genes were obtained by the cytohubba. Meanwhile, whether the hub genes are related to OS is evaluated by the Kaplan-Meier method.

Statistical Analysis

All statistical analyses were accomplished with R language 4.1.1 version and affiliated packages.

Results

Screening the Differentially Expressed miRNAs

As a result of comparing the 13 miRNA expression levels from TCGA data (457 cancer samples and 8 non-cancer samples), 8 miRNAs were found to differ in expression. The expression of all 8 miRNAs was elevated in tumor specimens, and the results were shown by heatmap ().

Building the Prognostic miRNA-Based Signature

We matched samples from the TCGA-COAD study with 427 patients who had complete data on their survival information. A univariate Cox regression revealed 5 miRNAs (P-value < 0.05) were related to patients’ OS (Figure 1A). However, considering that all differentially expressed miRNAs had clinical significance, we included all 8 differentially expressed miRNAs in further analysis (Figure 1B). Thus, the risk score was derived as follows: miRNA signature risk score = (0.307 × has-miR-141-3pexp) + (0.151 × has-miR-148a-3pexp) + (0.382 × has-miR-16-5pexp) - (0.411 × has-miR-200a-5pexp) + (0.034 × has-miR-223-3pexp) + (0.238 × has-miR-425-5pexp) + (0.951 × has-miR-500a-3pexp) - (0.043 × has-miR-7-5pexp). Using the risk score, 427 samples were evenly carved up into different risk subgroups. According to Kaplan-Meier curves, the high-risk group had a lower OS (Log rank test statistics = 8.2, HR = 1.850, 95% CI = 1.220–2.804, P = 0.004, Figure 2A). The sensitivity and specificity of the prognostic model were evaluated by ROC analysis, and we found that the AUC of 1-, 3-, and 5-year survival were 0.663, 0.653 and 0.639, respectively (Figure 2B).

Figure 1

Building the prognostic miRNA-based signature. (A) A univariate Cox analysis revealed that five miRNAs were related to overall survival. (B) A multivariate Cox regression model was constructed for all eight miRNAs with clinical differences.

Figure 2

Evaluation of the necroptosis-related miRNAs risk score. (A) The Kaplan Meier curve showed that the overall survival of the high-risk group was lower than that of the low-risk group. (B) The results of the ROC curve show that the model has a particular prediction ability.

Independence Prognostic Value of the miRNA Signature Model

Cox regression analyses were used to evaluate whether the miRNA signature model could be used as an independent prognostic determinant. Univariate analysis revealed that the risk score was related to OS (HR = 1.856, 95% CI: 1.206–2.856, P = 0.005, Figure 3A). The multivariate analysis also implied that the risk score was an independent prognostic factor considering other confounding factors (HR = 1.847, 95% CI = 1.197–2.848, P = 0.006, Figure 3B), such as age, gender and TMN stage. Therefore, this risk score can be used to evaluate the prognosis of colon cancer. Further, the TMN stage was also considered as an independent indicator of prognosis (HR = 3.511, 95% CI = 2.230–5.528, P < 0.001).

Figure 3

Independence prognostic value of the miRNA signature model. (A) Univariate analysis revealed that risk score was related to overall survival. (B) Multivariate analysis implied that risk score was an independent prognostic factor for colon cancer. The higher the risk score, the worse the prognosis of colon cancer.

miRNA-Target Genes Analysis

We analyzed the expression of eight miRNAs on patient survival using the Kaplan-Meier curve. Significant effects on OS outcomes were observed from the presentation of hsa-miR-500a-3p (Log rank test statistics = 9.1, HR = 1.869, 95% CI = 1.207–2.895, P = 0.003), hsa-miR-16-5p (Log rank test statistics = 8.5, HR = 1.869, 95% CI = 1.233–2.833, P = 0.004), hsa-miR-148a-3p (Log rank test statistics = 4.4, HR = 1.555, 95% CI = 1.015–2.382, P = 0.035) (Figure 4A–C). According to the data from three databases, each miRNA regulated different target genes. In total, 539 genes were obtained as the target genes of the three miRNAs (Figure 4D). To explore the association between miRNAs and target genes, we constructed this action network using Cytoscape 3.8.0, as shown in . The network consists of 542 nodes and 591 edges, of which hsa-miR-16-5p, hsa-miR-148a-3p and hsa-miR-500a-3p correspond to 421, 137 and 33 target genes, respectively. The average of node degree and clustering coefficient in this network were 2.2 and 0.0, respectively. Based on GO enrichment and KEGG analysis, we determined how these target genes function biologically. A total of 338 terms about the target genes associated with colon cancer were obtained by GO analysis (). As a result, the dot plot of the top 5 terms of biological process (BP), cellular component (CC), and molecular function (MF) terms were displayed (Figure 5A). In the three categories, changes in the BP of target genes were enriched in histone modification, gland development, and dephosphorylation; changes in MF were mainly enriched in protein serine/threonine kinase activity, protein serine kinase activity, and SMAD binding, while CC were mainly enriched in late endosome, transferase complex, and phosphorus−containing groups recycling endosome. KEGG pathway analysis of target genes associated with colon cancer revealed 44 results (), of which counts ≥ 20 were mainly enriched in the PI3K-Akt signaling pathway, Human papillomavirus infection, MicroRNAs in cancer, MAPK signaling pathway, and mTOR signaling pathway (Figure 5B).

Figure 4

Screening miRNAs related to the prognosis of colon cancer and predicting their target genes. (A) The higher the expression of has-miR-500a-3p, the worse the prognosis of colon cancer. (B) The high expression of has-miR-16-5p also leads to a poor prognosis of colon cancer. (C) The has-miR-148a-3p is also an adverse factor for the prognosis of colorectal cancer. (D) Venn diagram showed that 539 genes were obtained as the target genes of the three miRNAs.

Figure 5

The bubble diagram showed the biological functions of these 539 target genes. (A) The results of gene ontology analysis of target genes. (B) The results of KEGG pathway analysis of target genes.

Screening Hub Genes

The PPI network complex containing 498 nodes and 1931 edges consists of 498 target genes. This network’s average degree and clustering coefficient were 7.8 and 0.2, respectively. Ten hub genes (CCND1, SMAD3, SMAD2, CDK1, TGFB2, CDC25A, CHEK1, VEGFA, CCNE1, WEE1) were screened by Cytoscape 3.8.0 and its plug-in (). The node attributes information of these hub genes, such as degree, betweenness centrality, closeness centrality, and clustering coefficient, was shown in Table 1. A Kaplan-Meier analysis was performed to determine how hub genes expression affected patient survival, revealing that only the CHEK1 (Log rank test statistics = 6.7, HR = 0.600, 95% CI = 0.406–0.886, P = 0.010) exerted significant effects on OS ().

Table 1

The Node Attributes Information of Ten Hub Genes Based on PPI Network Analysis

Hub Genes	Degree	Betweenness Centrality	Closeness Centrality	Clustering Coefficient
CCND1	59	0.069544	0.411519	0.149036
SMAD3	41	0.036435	0.386667	0.178049
SMAD2	33	0.016226	0.369288	0.242424
CDK1	42	0.060424	0.404098	0.14518
TGFB2	15	0.002259	0.325198	0.495238
CDC25A	26	0.01503	0.371235	0.292308
CHEK1	27	0.016257	0.380108	0.262108
VEGFA	42	0.043709	0.39127	0.159117
CCNE1	22	0.008628	0.36846	0.324675
WEE1	17	0.001974	0.354167	0.485294

Notes: Discussion. Degree is the number of connections between nodes. Betweenness centrality refers to the number of times one node acts as the shortest bridge between the other two nodes. Closeness centrality represents the average length of the shortest path from one node to all other nodes. Clustering coefficient is the cluster tendency of nodes in the network.

Discussion

Research has shown that necroptosis plays a significant part in cancer progression and metastasis, and cancer treatment based on targeting necroptosis-inducing factors is a strategy for cancer treatment.18 In the current study, this is the first time to explore the use of miRNAs related to necroptosis to predict the prognosis of colon cancer. This study first analyzed the expression of 13 miRNAs associated with necroptosis in colon cancer and normal tissues, then constructed a prognostic risk score composed of 8 miRNAs to evaluate the predictive value of these necroptosis-related miRNAs. To clarify the mechanism of these necroptosis-related miRNAs that clearly impacted patient survival, we predicted their target genes through three databases. The functional analyses indicated that the differentially expressed miRNAs are highly connected with the PI3K-Akt and MAPK signaling pathways.

Our study generated a signature risk score composed of eight necroptosis-related miRNAs and detected that it could be forecasted OS in colon cancer. In tumor patients, tumor-promoting miRNA was upregulated, and tumor inhibitory miRNA was downregulated.19 The progression and poor prognosis of metastatic colorectal cancer are linked to increased circulating levels of miR-141-3p. According to reports, miR-141-3p targets PTEN, which in turn causes phosphoinositide to activate the PI3K/Akt pathway.20 What’s more, in Caco-2 cells treated with LPS, miR-141-3p inhibited the effect of RIPK1, which activated downstream molecules and promoted necrotic body formation, thereby reversing intestinal epithelial cell injury.15 Clinically, the serum levels of exosome miR-425-5p were joined to progression and metastasis of CRC, which can regulate PTEN by activating PI3K/Akt signaling pathway, inducing M2 polarization of macrophages, enhancing epithelial-mesenchymal transition and secret VEGF to promote tumor metastasis.21,22

MiR-7-5p mediates mitochondrial damage by targeting mitochondrial proteins SLC25A37 and TIMM50, thereby inducing apoptosis and necroptosis of rhabdomyosarcoma cells.23 Recently, lncRNA TINCR could serve as a miR-7-5p sponge was proposed to be associated with colorectal cancer progression, which may promote the development of colorectal cancer mediated by miR-7-5p.24 MiR-223-3p negatively manages the expression of RIP3, and it significantly affects the inhibition of RIP3-mediated necroptosis and inflammation to reduce spinal cord neuron damage.25 LINC00961, as a tumor suppressor gene, up-regulates the expression of its target gene SOX11 by acting as a miR-223-3p sponge to inhibit the migration and invasion of colorectal cancer.26 Compared with normal colonic mucosa, miR-200a-5p was significantly up-regulated in colorectal cancer, but colon adenoma is not substantially different from normal mucosa.27,28 It was reported that miR-200a-5p and its target gene RNF11 were differentially expressed in cardiac tissue, and overexpressed miR-200a-5p induces RIP3-dependent necroptosis.29

In our study, the expression of hsa-miR-500a-3p, hsa-miR-16-5p, and hsa-miR-148a-3p significantly affected the OS outcomes. Studies have shown that miR-500a-3p can bind to the 3ʹUTR of MLKL,30 promote the characteristics of tumor stem cells and drive the activation of STAT3 signaling by targeting a variety of negative regulators of the JAK/STAT3 pathway.31 Although MiR-16-5p is considered a tumor suppressor,32,33 its mechanism of action with colorectal cancer is controversial.34 The content of miR-16-5p in serum and plasma of patients with colorectal cancer was significantly more abundant than that of healthy individuals,35 and a similar phenomenon was also observed in the feces of patients with precancerous lesions.36 In colorectal cancer, the inhibition of tumor promoter miR-148a-3p can restore the surface level of MHC-I and significantly enhance the effect of CD8+ T cell-mediated immune attack by promoting the expression of calnexin.37 However, a study reported that the low expression of miR-148a is associated with the significant reduction of disease-free survival and overall survival,38 indicating that the mechanism between miR-148a-3p and colorectal cancer is not precise.39

By analyzing accessible online databases, we argued the potential target genes and functions of the above three miRNAs. Target genes mainly took part in histone modification, protein serine/threonine kinase activity, and SMAD binding, targeting the PI3K-AKT and MAPK pathway, which played an essential role in the progression and metastasis of cancer. A report revealed that the inhibition of TIMP1 could reduce proliferation and metastasis by regulating FAK-PI3K/AKT and MAPK pathway but increase apoptosis.40 It was reported that emodin induces high phosphorylation of RIP1 and MLKL through ROS mediated PI3K/AKT, thereby inducing necroptosis of renal cells carcinoma.41 In glioblastoma, RIP1 and RIP3 as upstream of ROS-JNK-p53 loops, mediate the mammalian target of PI3K/AKT/mTOR/P70S6K, thereby inducing necroptosis.42

Using Cytoscape 3.8.0 and its plug-in to screen out 10 hub genes (CCND1, SMAD3, SMAD2, CDK1, TGFB2, CDC25A, CHEK1, VEGFA, CCNE1, WEE1) that regulate colon cancer miRNA. Cyclins are a group of proteins that regulate the cell cycle by binding to cyclin-dependent kinases (CDKs).43 Studies have shown that CCND1 and CCNE1 were differentially expressed between normal and cancer colon tissues, and the latter was significantly increased in the early stage of TNM and was significantly related to OS.44 In epithelial ovarian cancer, miR-490-3p regulated and affected the expression of SMARCD1 and CCND1 by targeting CDK1.45 CDK1 could form CCNB/CDK1 complex with CCNB to control the G2/M phase transition and was crucial for initiating mitosis.46 CDC25A was also a cell cycle-related gene, which was overexpressed in many types of cancer and promotes tumor progression by regulating the cancer cell cycle and may be a key target for effective cancer therapy.47 SMAD protein was TGF-β Signal transduction molecules downstream of family receptors. In cancer cells from different sources, acidosis activated TGF-β2, promoting partial epithelial-mesenchymal transition, which supported SMAD2 acetylation.48 There is a lot of evidence that TGF-β2 was abundant in the glioma microenvironment, which induced autophagy played a vital role in glioma invasion. This autophagy also initiated the feedback of TGF-β2 in glioma by affecting the expression of SMAD2/3/7.49 Angiogenesis was a necessary condition for tumor progression. VEGFA secreted by tumor cells and the surrounding matrix was overexpressed in most human tumors, which stimulated endothelial cell proliferation and survival, closely related to tumor invasion, metastasis, and prognosis.50 Like other DNA damage response (DDR) related kinases, the primary biological function of WEE1 as a tumor suppressor was to prevent DNA altered cell replication. However, more and more studies have shown that WEE1 was overexpressed in various solid tumors and was associated with a poor prognosis.51

It is worth noting that CHEK1 is an essential gene and associated with the OS of patients. Kaplan-Meier showed that high CHEK1 expression is linked to a higher survival rate. CHEK1 was a serine/threonine protein kinase, a significant regulator of DNA damage response.52 Because it played a regulatory role in DNA damage, CHEK1 was previously considered as a tumor suppressor gene. However, a study found that the differential expression of CHEK1 and the decrease of CHEK1 mRNA are poor prognostic factors for gastric cancer and colorectal cancer.53 Another study showed that treatment with specific miR inhibitors in colon cancer reduced the high level of endogenous DNA damage sensitive miRs (DDSMs), increasing CHEK1 transcription level, which revealed the inhibitory effect of CHEK1 on colon cancer.54

We found a novel miRNA signature risk score and analyzed its predictive effect on the prognosis of colon cancer. We have no conditions to collect our own clinical data to verify our model, which is a limitation of our study. In addition, we have no experimental conditions to verify the action mechanism of these miRNAs and targets, which is also the limitation of our research.

Conclusions

In conclusion, our study shows that necroptosis is closely related to colon cancer because there are differences in miRNA expression between normal and colon cancer tissues. In addition, our model of eight necroptosis-related miRNAs can be used as independent prognostic molecular markers of colon cancer. In our study, we have discovered a new miRNA marker that could predict the prognosis of colon cancer and help identify reliable therapeutic targets for colon cancer. The ten hub genes we screened may be effective targets for targeted cancer therapies, especially CHEK1. In addition, they have the potential to become biomarkers of CRC. However, more experimental studies and clinical trials are needed to verify these targets.