- 1Department of Gastroenterology, Affiliated Hospital of North Sichuan Medical College, Nanchong, China
- 2Department of Yunnan Tumor Research Institute, The Third Affiliated Hospital of Kunming Medical University, Yunnan Cancer Hospital, Kunming, China
Background: Colorectal cancer (CRC) is the second leading cause of cancer-related deaths in the world. This study aimed to develop a urea cycle (UC)-related gene signature that provides a theoretical foundation for the prognosis and treatment of patients with CRC.
Methods: Differentially expressed UC-related genes in CRC were confirmed using differential analysis and Venn diagrams. Univariate Cox and least absolute shrinkage and selection operator regression analyses were performed to identify UC-related prognostic genes. A UC-related signature was created and confirmed using distinct datasets. Independent prognostic predictors were authenticated using Cox analysis. The Cell-type Identification by Estimating Relative Subsets of RNA Transcripts algorithm and Spearman method were applied to probe the linkage between UC-related prognostic genes and tumor immune-infiltrating cells. The Human Protein Atlas database was used to determine the protein expression levels of prognostic genes in CRC and normal tissues. Verification of the expression levels of UC-related prognostic genes in clinical tissue samples was performed using real-time quantitative polymerase chain reaction (qPCR).
Results: A total of 49 DEUCRGs in CRC were mined. Eight prognostic genes (TIMP1, FABP4, MMP3, MMP1, CD177, CA2, S100P, and SPP1) were identified to construct a UC-related gene signature. The signature was then affirmed using an external validation set. The risk score was demonstrated to be a credible independent prognostic predictor using Cox regression analysis. Functional enrichment analysis revealed that focal adhesion, ECM-receptor interaction, IL-17 signaling pathway, and nitrogen metabolism were associated with the UC-related gene signature. Immune infiltration and correlation analyses revealed a significant correlation between UC-related prognostic genes and differential immune cells between the two risk subgroups. Finally, the qPCR results of clinical samples further confirmed the results of the public database.
Conclusion: Taken together, this study authenticated UC-related prognostic genes and developed a gene signature for the prognosis of CRC, which will be of great significance in the identification of prognostic molecular biomarkers, clinical prognosis prediction, and development of treatment strategies for patients with CRC.
Introduction
According to global cancer data, colorectal cancer (CRC) ranks third in incidence and second in mortality rate worldwide. GLOBOCAN estimated more than 1.9 million new CRC cases and 935,000 deaths in 2020 (1, 2). The prognosis and survival status of patients with CRC have not improved significantly over the years (3). The prognosis of CRC is predicted using the tumor–node–metastasis (TNM) staging system, histopathologic criteria, and tumor markers, but it cannot accurately predict the clinical prognosis (4). The etiology of CRC is complex, with approximately 10% patients having susceptibility to germline mutations that lead to familial cases. However, the majority of patients with CRC have sporadic cancer caused by a combination of environmental and genetic risk factors (5), although the specific molecular mechanisms remain unknown. Therefore, the molecular mechanisms underlying the development of CRC should be explored further, and novel biomarkers for prognostic assessment should be identified to improve the clinical prognosis of patients. It is challenging to predict the prognosis of CRC owing to its rapid progression and highly heterogeneous nature (6). This necessitates the development of new prognostic models.
The process by which ammonia produced during the metabolism of amino acids in the body is converted into urea via ornithine constitutes the urea cycle (UC). UC eliminates the excess nitrogen and ammonia produced by the breakdown of proteins or the synthesis of nitrogenous compounds in the body. UC enzymes also manipulate nucleotide metabolism in certain types of tumors. UC-related genes (UCRGs) are overexpressed or silenced in different types of cancers, and altered UC gene expression is actively involved in tumorigenesis (7). UC dysregulation (UCD) is observed in various cancer types and is associated with a poor prognosis, but is responsive to immunotherapy (8). Tumor cells reprogram their metabolism to maximize the use of nitrogen and carbon to obtain sufficient energy for tumor proliferation and rapid growth (9, 10). p53 inhibits tumor growth by regulating ammonia metabolism via UC, thereby controlling polyamine synthesis in tumor cells (11). UC enzymes play a critical role in killing cancer cells and suppressing cancer growth. In the human hepatoblastoma cell line, HepG2, deficiency of arginase 1 and ornithine carbamoyltransferase (OTC) results in high ammonia levels and diminished UC function (12). OTC and argininosuccinate synthase (ASS) levels are deficient in acute lymphoblastic leukemia (13). The UC enzyme, carbamoyl-phosphate synthase 1 (CPS1), maintains the pyrimidine pool in non-small cell lung cancer by activating CAD, and silencing CPS1 in KL cells induces cell death and inhibits tumor growth in vivo due to the depletion of pyrimidines (14). ASS1 expression promotes CRC cell proliferation and tumor formation in vitro (15). Whether UCRGs are associated with the prognosis of patients with CRC remains unclear; therefore, systematic studies of prognosis-related UC genes can help to understand their role in CRC development and progression and aid in guiding clinical decisions.
This study aimed to determine the prognostic value of UCRGs in patients with CRC. RNA-sequencing (seq) data for CRC and the corresponding clinical data were extracted from public databases. UC-associated differentially expressed genes (DEGs) closely associated with prognosis were identified to construct a predictive model for CRC prognosis in The Cancer Genome Atlas (TCGA) cohort. We then validated it in the Gene Expression Omnibus (GEO) cohort, determined the expression levels of eight genes in colon cancer tissues using quantitative reverse transcription-polymerase chain reaction (qRT-PCR), and obtained results consistent with our initial prediction.
Methods
Gene and dataset collection
We integrated the transcriptomic data of 673 CRC and 51 normal tissue samples using TCGA database (https://portal.gdc.cancer.gov/repository). Seven CRC datasets were mined from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) (GSE41258, GSE110223, GSE110224, GSE113513, GSE17538, GSE39582, and GSE44076) and used in this study. The number of tissue samples and their usage in each dataset are listed in Table 1. A total of 2,857 UCRGs were derived from the GeneCard database (https://www.genecards.org/) by searching “UC” and are listed in Supplementary Table S1. A flow chart of the present study is shown in Supplementary Figure S1.
Identification of differentially expressed UCRGs (DEUCRGs) in CRC
On the basis of p value < 0.05 and |log2FoldChange(FC)| > 1, we first identified the DEGs (CRC samples vs. normal samples) in the TCGA-CRC, GSE41258, GSE110223, GSE110224, and GSE113513 datasets using the limma package (16). We separately crossed the upregulated and downregulated genes in the five datasets to identify the shared DEGs, which were further intersected with UCRGs to identify DEUCRGs for subsequent prognostic analysis.
Functional annotation analysis
We used the R package clusterProfiler (17) for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. GO was categorized as cellular component (CC), molecular function (MF), and biological process (BP). The significance criterion was adjusted p value ≤ 0.05.
Establishment of UC-related gene signature in CRC
The 590 patients with CRC with survival information in TCGA-CRC dataset served as the training set and 232 patients with survival information from the GSE17538 dataset were used as the external validation set to create and verify the UC-related gene signature for predicting the survival of patients with CRC. Univariate Cox and least absolute shrinkage and selection operator (LASSO) regression analyses were used to select the UC-related prognostic genes in the training set. Depending on the risk score formula (Riskscore = , coef represents the coefficient obtained by LASSO), and cut-off value calculated using the surv_cutpoint function in the cutoff package, patients were classified into two risk subgroups: high- and low-risk groups. Kaplan–Meier curves, receiver operating characteristic (ROC) analysis, and risk curves were used to demonstrate the predictive efficiency of the gene signature.
Relevance analysis of UC-related gene signature and clinical parameters
Risk scores for different clinical factor subgroups were compared using the Wilcoxon or Kruskal–Wallis test. Stratified survival analysis for different clinical subgroups was also performed by developing K–M curves.
Independent prognostic analysis and nomogram creation
The risk score, age, sex, race, stage, pathologic_M, pathologic_N, and pathologic_T were included in Cox analyses (univariate Cox and multivariate Cox) to authenticate independent prognostic predictors. The nomogram comprising the independent prognostic predictors was drawn using R language “rms” to predict survival at 1, 3, and 5 years in patients with CRC. The corresponding calibration and DCA curves were plotted to appraise the precision and reliability of the nomogram model predictions.
Relevance analysis of UC-related prognostic genes and immune infiltration
Discrepancies in 22 types of immune infiltrating cells were assessed and compared between the two risk groups using the Cell-type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm (18) and Wilcoxon test. We calculated the correlation between UC-related prognostic genes and differential immune cells using the Spearman's method.
Analysis and validation of the expression levels of UC-related prognostic genes
We first verified the discrepancy in the expression levels of UC-related prognostic genes in CRC and normal samples in the external datasets, GSE39582 and GSE44076. Corresponding box-line plots were obtained using the ggplot2 package. We used immunohistochemistry images from the Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) to further determine the protein expression levels of UC-related prognostic genes in normal and CRC tissues. To further confirm the results of the public database analysis, we collected six normal tissue samples and six CRC tissue samples from the Affiliated Hospital of North Sichuan Medical College in June and performed RNA isolation and qRT-PCR. This study was approved by the Affiliated Hospital of North Sichuan Medical College (Nanchong, China) (Ethical Application Ref: 2022ER237-1). Total RNA from 12 samples was separated using TRIzol (Ambion, USA), according to the manufacturer's instructions. The inverse transcription of total RNA into cDNA was implemented using the First-strand cDNA synthesis kit (Servicebio, China), according to the manufacturer's instructions. Then, qPCR was carried out using the 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, China), according to the manufacturer's instructions. The primer sequences used for PCR are listed in Supplementary Table S2. Expression was normalized to the internal reference glyceraldehyde 3-phosphate dehydrogenase and computed using the 2−ΔΔCq method (19).
Statistical analysis
All bioinformatics analyses were performed using the R language. Wilcoxon test or Kruskal–Wallis test was used to compare data from different groups. Student's t-test was used to compare the discrepancies in qRT-PCR.
Results
DEUCRGs in CRC
Based on |log2FC| > 1 and p value < 0.05, we first identified the DEGs (CRC samples vs. normal samples) in TCGA-CRC, GSE41258, GSE110223, GSE110224, and GSE113513 datasets. Volcano plots of DEGs for each dataset are shown separately in Figures 1A–E. The specific numbers of DEGs per dataset are presented in Table 2. We intersected the upregulated and downregulated genes in the CRC samples from each of the above five datasets, resulting in 39 crossed upregulated genes and 102 crossed downregulated genes (Figures 1F–G; Supplementary Table S3). We then intersected the 141 crossed differential genes mentioned above with the 2,857 UCRGs derived from the GeneCard database, resulting in 49 crossed genes, namely DEUCRGs in CRC (Figure 1H; Supplementary Table S3). To further probe the function of the DEUCRGs in CRC, functional enrichment analysis was performed. As shown in Supplementary Table S4, 149 GO items (92 BP, 20 CC, and 37 MF items) and 17 KEGG pathways were identified. The top 10 items under each classification are shown in the bubble diagram (Figures 1I–J). These genes were mainly linked to immune-, extracellular matrix-, hormone metabolism-, and nitrogen metabolism-related biological processes. Furthermore, these genes were implicated in “Bile secretion,” “Nitrogen metabolism,” “IL-17 signaling pathway,” “Proximal tubule bicarbonate reclamation,” “PPAR signaling pathway,” “Pyruvate metabolism,” “NF-kappa B signaling pathway,” and “Chemokine signaling pathway.”
Figure 1. Identification of colorectal cancer (CRC)-associated differentially expressed UC-related genes (DEUCRGs) in the cancer genome atlas (TCGA) and gene expression omnibus (GEO) cohorts. Volcano plots of DEGs in (A) TCGA-CRC, (B) GSE110223, (C) GSE41258, (D) GSE110224, and (E) GSE113513 datasets. (F) Venn diagram of upregulated genes in five datasets. (G) Venn diagram of downregulated genes in five datasets. (H) Venn diagram for identifying DEUCRGs. (I) Top 10 Gene Ontology (GO) terms enriched by DEUCRGs under cellular component (CC), molecular function (MF), and biological process (BP) subcategories. (J) Top 10 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched by DEUCRGs.
UC-related gene signature to assess the prognosis of patients with CRC
A total of 590 patients with CRC with survival information in TCGA-CRC dataset served as the training set. To mine the UCRGs relevant to the overall survival (OS) of patients with CRC, we incorporated the 49 DEUCRGs obtained above into a univariate Cox analysis in the training set. Fourteen genes associated with OS of patients with CRC were selected (p < 0.1), with tissue inhibitor of metalloproteinase-1 (TIMP1), fatty acid-binding protein 4 (FABP4), secreted phosphoprotein 1 (SPP1), and butyrylcholinesterase as risk factors for CRC prognosis [hazard ratio (HR) > 1], and matrix metallopeptidase (MMP)-3, MMP1, C-X-C motif chemokine ligand (CXCL)-1, CXCL3, CD177, carbonic anhydrase (CA)-2, S100 calcium-binding protein P (S100P), CA4, nuclear receptor subfamily 3 group C member 2, and CXCL8 as protective factors for CRC prognosis (HR < 1) (Table 3). The above fourteen genes were further integrated into LASSO analysis with a 20-fold cross-validation. As shown in Figure 2A, when lambda min was 0.0140, the corresponding number of genes was eight. Therefore, eight genes (TIMP1, FABP4, MMP3, MMP1, CD177, CA2, S100P, and SPP1) and their corresponding coefficients were determined as UC-related prognostic genes for signature establishment (Table 4). Next, we developed a risk signature based on the following formula: RiskScore = (–0.02030967)×expression (MMP3) + (–0.099075148)×expression (MMP1) + (–0.041645445)×expression(CD177) + (–0.019591304)×expression (CA2) + (–0.074156971)×expression (S100P) + 0.356587004×expression (TIMP1) + 0.06708291×expression (FABP4) + 0.00844555 × expression (SPP1). Based on this formula, we calculated the risk score for each patient with CRC in the training set and classified them into high- and low-risk groups based on the cutoff value. The K–M curve illustrated that patients with higher risk had significantly poorer survival than those with lower risk (Figure 2B). We plotted the ROC curves to check the predictive efficiency of the signature. The area under the curve values for OS in the training set were 0.624 (1 year), 0.658 (3 years), and 0.735 (5 years), indicating decent accuracy (Figure 2C), the cut-off values of the ROC curves were shown in Supplementary Figures S2A–C. Figure 2D shows the distribution of the ranked risk score and survival status for each patient in the training set. Survival status showed that, as the risk score increased, patients had a relatively high risk of death. The expression heatmap showed that FABP4, TIMP1, and SPP1 were highly expressed in patients with high risk scores, whereas S100P, MMP3, MMP1, CD177, and CA2 were highly expressed in patients with low risk scores (Figure 2E). To further prove the applicability and reliability of the risk signature, the above analysis was carried out in the external validation set (GSE17538 dataset). A comparable trend was observed in the external validation set (Figures 2F–I), the cut-off values of the ROC curves were shown in Supplementary Figures S2D–F. The above results confirmed that the UC-related gene signature is a valid survival predictor for patients with CRC.
Figure 2. Establishment of a UC- relevant risk score model. (A) Least absolute shrinkage and selection operator (LASSO) regression analysis to establish an eight-gene signature. (B) Kaplan–Meier curves of the high- and low-risk groups in TCGA-CRC dataset. (C) Receiver operating characteristic (ROC) curves of 1, 3, and 5 years overall survival (OS) in TCGA-CRC dataset. (D) Signature-based distribution of risk scores and survival status in TCGA-CRC dataset. (E) Differences in prognostic gene expression in different risk groups in TCGA-CRC dataset. (F) Kaplan–Meier Curve of the high- and low-risk groups in GSE17538 dataset. (G) ROC curves of 1, 3, and 5 years OS in GSE17538 dataset. (H) Distribution of risk scores and survival status in GSE17538 dataset. (I) Differences in prognostic gene expression in different risk groups in GSE17538 dataset.
Relevance analysis of clinical parameters and UC-related gene signature
To explore the relationship between the UC-related gene signature and clinical factors, we compared the risk scores of different clinical subgroups of patients. As shown in Figures 3A–C, UC-related risk scores were independent of sex and age and correlated with race. Additionally, we found that the UC-related risk scores were notably correlated with stage, pathologic_M, pathologic_N, and pathologic_T stage, and the risk scores tended to increase as the malignancy of the tumor (Figures 3D–G). We further explored the application of the risk score in patients with different clinicopathological characteristics. The results of the stratified survival analysis uncovered significant differences in the survival of the two groups for most clinical subgroups, including age ≤ 65 years, age > 65 years, female sex, white ethnicity, M0, N0, N1/N2, T3/T4, stage I/II, and stage III/IV subgroups (Figure 4).
Figure 3. Comparison of risk scores for different subgroups of clinical parameters. (A) Age > 65 vs. Age ≤ 65. (B) Female vs. Male. (C) Asian vs. Black vs. White. (D) T1 vs. T2 vs. T3 vs. T4. (E) N0 vs. N1 vs. N2. (F) M0 vs. M1. (G) Stage I vs. II vs. III vs. IV.
Figure 4. Stratified survival analysis of subgroups with different clinical parameters. (A) Age > 65. (B) Age ≤ 65. (C) Female. (D) Male. (E) White. (F) T1/T2. (G) T3/T4. (H) N0. (I) N1–N2. (J) M0. (K) M1. (L) Stage I/II. (M) Stage III/IV.
UC-related gene signature is an independent prognostic predictor
Using Cox analyses (univariate and multivariate Cox), we determined that risk score, age, and pathologic_T were independent prognosis predictors in patients with CRC (Figures 5A,B; p value < 0.05). A nomogram containing independent prognostic predictors was generated (Figure 5C), and the calibration curves proved that the performance of the nomogram in predicting the survival at 1, 3, and 5 years in patients with CRC was satisfactory (Figure 5D), whereas the DCA curves revealed that the nomogram had a higher accuracy in predicting the survival of patients at 3 and 5 years than the individual independent predictors (risk score, age, and T) (Figures 5E–F).
Figure 5. Nomogram predicting the OS of patients with CRC. (A) Univariate Cox regression analysis in TCGA cohort. (B) Multivariate Cox regression analysis in TCGA cohort. (C) Nomogram predicting patient survival at 1, 3, and 5 years. (D) Calibration curves of nomogram. (E) DCA curve for predicting 3-year survival of patients. (F) DCA curve for predicting 5-year survival of patients.
Association of UC-related gene signature with immune infiltrating cells
To uncover potential mechanisms for the disparity in the prognosis of the two risk groups, we performed the authentication and functional enrichment analysis of DEGs between the two groups. As shown in Figures 6A–B, a total of 73 DEGs (high-risk vs. low-risk), containing 40 highly expressed and 33 lowly expressed genes, were identified. Accordingly, 24 BP items, 36 CC items, 24 MF items, and seven KEGG pathways were identified (Supplementary Table S5). The top 10 GO items for each classification are shown in Figure 6C. We observed that the above genes were principally linked to the humoral immune response, antibacterial humoral response, and extracellular matrix organization. Furthermore, these genes were implicated in focal adhesion, ECM-receptor interaction, IL-17 signaling pathway, and nitrogen metabolism (Figure 6D). As immune-relevant biological processes and pathways were revealed to be connected to the UC-related gene signature, we investigated the changes in immune infiltration between the two risk groups. Using the CIBERSORT algorithm, 124 samples in the high-risk group and 430 samples in the low-risk group were incorporated to compute the fraction of each immune infiltration cell after excluding samples with p value > 0.05 (Figure 7A). As shown in the violin plot, the fraction of B cell memory, T-regulatory cells, monocytes, and macrophages M0 and M2 was elevated in patients with high risk scores, while the fraction of T cells CD4 memory resting, T cells CD4 memory activated, plasma cells, activated dendritic cells, and neutrophils was higher in patients with low risk scores (Figure 7B). We further calculated the relevance of the UC-related prognostic genes and the 10 differential immune cells mentioned above using the Spearman method (Supplementary Table S6). Based on the thresholds of |cor| > 0.3 and p < 0.05, we discovered that SPP1 was positively correlated with macrophage M2 and neutrophils and negatively correlated with plasma cells (Figure 7C). MMP3 expression was significantly and positively correlated with dendritic cell activation and neutrophils (Figure 7C). MMP1 expression was significantly positively correlated with neutrophil levels (Figure 7C). FABP4 was positively correlated with macrophage M2 (Figure 7C). CA2 expression was significantly negatively correlated with macrophages M0 (Figure 7C).
Figure 6. Functional enrichment of DEGs between high- and low-risk groups. (A) Volcano map of DEGs between high- and low-risk groups. (B) Heat map of DEGs between high- and low-risk groups (C) Top 10 GO terms enriched by DEGs between high- and low-risk groups under BP, CC, and MF subcategories. (D) KEGG pathways enriched by DEGs between high- and low-risk groups.
Figure 7. Differences in immune infiltrating cells between the high- and low-risk patients with CRC. (A) Rrelative percent of immune infiltrating cells in each CRC sample. (B) Violin plot of immune infiltrating cells between high- and low-risk groups. (C) Correlation between UC-related prognostic genes and differential immune infiltrating cells.
Expression levels of UC-related prognostic genes
As illustrated in Figure 8A, CA2, CD177, and FABP4 levels were downregulated, while MMP1, MMP3, S100P, SPP1, and TIMP1 levels were elevated in CRC tissues compared to normal tissues in TCGA-CRC dataset. We further confirmed the same expression trend in the two external validation sets (GSE39582 and GSE44076) (Figures 8C–D). To further determine the changes in the expression of UC-related prognostic genes at the protein level, we obtained the corresponding immunohistochemistry images from HPA database. We did not detect immunohistochemical result of MMP1 CRC. As shown in Figure 9, we found that the protein expression levels of CA2, CD177, and FABP4 were reduced in CRC tissues than in normal tissues. SPP1 was largely unexpressed in normal and CRC tissues at the protein level. Moreover, protein expression levels of MMP3, S100P, and TIMP1 were increased in CRC tissues than in normal tissues. We verified the protein expression in clinical tissue samples via qRT-PCR. In agreement with the results of the public database data analysis, CA2, CD177, and FABP4 were expressed at low levels, while MMP1, MMP3, S100P, SPP1, and TIMP1 were highly expressed in clinical CRC samples compared to normal samples (Figure 10).
Figure 8. Expression levels of UC-related prognostic genes in three datasets. Expression levels of risk model genes in (A) TCGA-CRC, (B) GSE39582, and (C) GSE44076 datasets. **** p < 0.0001.
Figure 9. Immunohistochemistry images of UC-related prognostic genes obtained from the human protein atlas (HPA) database.
Figure 10. Expression levels of eight UC-related prognostic genes in clinical tissues determined using quantitative reverse transcription-polymerase chain reaction (qRT-PCR). (A) CD177 (B) carbonic anhydrase 2 (CA2), (C) fatty acid-binding protein 4 (FABP4), (D) tissue inhibitor of metalloproteinase-1 (TIMP1), (E) secreted phosphoprotein 1 (SPP1), (F) S100 calcium-binding protein P (S100P), (G) matrix metallopeptidase (MMP)-3, and (H) MMP1. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Discussion
CRC is a molecularly heterogeneous malignancy with limited therapeutic strategies for patients with advanced disease. The molecular features of CRC are closely related to its prognosis, and there are no reliable prognostic biomarkers for risk prediction in clinical practice. Therefore, the identification of key biomarkers and therapeutic targets affecting prognosis is essential for improving the clinical outcomes in patients with CRC. As an important factor involved in cellular metabolic reprogramming, UC is involved in the genesis and development of various tumors, including CRC (20), revealing its potential as a prognostic biomarker. For the first time, we constructed and externally validated a novel prognostic model that integrated eight genes associated with UC. The genes constituting the prognostic model were TIMP1, FABP4, MMP3, MMP1, CD177, CA2, S100P, and SPP1 (Table 4).
TIMP1 has been shown to promote the progression of multiple tumors (21–24). Studies in CRC have indicated that upregulation of TIMP1 was associated with poor prognosis and confirmed that TIMP1 can promote tumor cell proliferation and metastasis through the FAK/Akt signaling pathway (22, 25). FABP4 is a mediator of lipid metabolism in adipocytes and can provide fatty acids to tumor cells (26). FABP4 has been found to promote the progression of ovarian cancer, cervical cancer, breast cancer, prostate cancer cell carcinoma, and oral squamous cell carcinoma (26–30). In colorectal cancer, high expression of FABP4 has been found to be closely related to tumor recurrence (31). However, one study showed that FABP4 can inhibit the proliferation and metastasis of CRC cells (32). In this study, high expression of FABP4 was associated with a poor prognosis for CRC patients, so the mechanism of FABP4 in CRC needs further study. The bone bridge protein, SPP1, is a key ECM protein involved in tumor progression and metastasis (33) that is highly expressed in non-small cell lung cancer tissues (34), significantly upregulated in glioma and hepatocellular carcinoma cell lines, and associated with poor prognosis (35, 36). SPP1 expression is upregulated in CRC tissues and is associated with the short OS of patients (37), which is consistent with our experimental and predictive modeling results. MMP3 is a member of the metalloproteinase family along with MMP1, which affects tissue integrity by degrading ECM components (21). MMP3 is a member of the metalloproteinase family along with MMP1, which affects tissue integrity by degrading ECM components (21). Studies have shown that MMP3 can promote the metastasis of CRC (38), melanoma (39), and breast cancer (40). However, studies in CRC also found that the expression of MMP3 in patients without metastasis was significantly higher than that in patients with distant metastasis (41, 42). Our study also found a higher expression of MMP3 in low-risk patients. It is possible that MMP3 triggers tumor metastasis and the expression of MMP3 decreases after the tumor metastasizes. Previous studies have shown that MMP1 has potential as a diagnostic and prognostic marker for CRC (43, 44). Several studies have also confirmed that in CRC, the expression of MMP1 is lower in metastatic CRC than in primary CRC (45, 46). It indicates that MMP1 may have a similar role to MMP3, acting only in the initial stage of tumor metastasis and appearing to be less important once metastasis occurs. Therefore, the specific mechanism of action of MMP3 and MMP1 in colorectal cancer needs to be further studied. CD177 is a glycosylphosphatidylinositol-linked cell surface protein that is heterogeneously expressed by neutrophils, and its expression is associated with good prognosis in breast, prostate, cervical, and lung cancers (47). Studies have shown that CD177+ neutrophils inhibit epithelial tumorigenesis and serve as independent predictors of prognosis in patients with CRC (48). CA2 is a factor that inhibits metastasis and EMT and is associated with good OS in patients with hepatocellular carcinoma (49). CA2 is lowly expressed in ulcerative colitis and CRC tissues (50), and low CA2 expression is associated with a better prognosis. S100P is highly expressed in various solid tumors and associated with poor prognosis in CRC (51, 52), breast cancer (53), pancreatic cancer (54), cholangiocarcinoma (55) lung cancer (56), and ovarian cancer (57). Our qRT-PCR results also indicated that S100P was highly expressed in cancer tissues.
Our results showed that the risk score of the eight DEUCRGs is an independent prognostic marker for CRC patient survival. The survival rate of patients with CRC in the high-risk group was significantly lower than that in the low-risk group. The prognosis of the two risk groups differed in subgroups based on age, sex, clinical stage, T-stage, lymph node metastasis, distant metastasis, and race (Figure 4). Univariate and multivariate Cox regression analyses were used to further investigate the independent prognostic value of the clinicopathological characteristics and risk score, and the results showed that risk score, age, and T stage were significant independent prognostic factors for CRC (Figures 5A, B). An individualized prognostic prediction model was then constructed using a nomogram to quantify the individual risk in the clinical setting by integrating multiple risk factors, including independent prognostic factors, and calibration curves showed a high degree of agreement between the actual and predicted OS rates. Based on the above results, we suggest that our constructed prognostic risk score model is a valid prognostic indicator for patients with CRC.
GO and KEGG enrichment analyses of DEGs in different risk groups were used to determine the roles of UC-related biological processes and classical signaling pathways in different risk groups. As predicted, the results showed significant enrichment of genes in biological processes, such as muscle contraction, antimicrobial humoral response, and muscle system processes (Figure 6C). Changes in ECM components contribute to cancer progression, promote tumor-associated angiogenesis and inflammation, and affect the tumor microenvironment (58). Immune-related signaling pathways and functions, such as the humoral immune response and IL-17 signaling pathway, wernie also sigficantly enriched (Figure 6D). One study showed that the UC is associated with immunity (59). Single-sample GSEA scores showed significant differences in immune cell infiltration scores between the two risk groups. The proportions of B memory cells, T regulatory cells, monocytes, and macrophages M0 and M2 were elevated in the high-risk group, while the proportions of CD4 memory resting and activated T cells, plasma cells, activated dendritic cells, and neutrophils were elevated in the low-risk group (Figure 7B). Different levels of immune cell infiltration exist in different risk groups, and these differences may be important factors affecting the prognosis and treatment response of patients. These results suggest that targeting UC-related genes can potentially improve the immune status of CRC or promote immunotherapy in CRC; however, further studies are needed to confirm this.
Li et al. constructed a prognostic signature that included 25 long non-coding RNAs, but models incorporating too many variables were difficult to implement, limiting their clinical application (60). In contrast, the predictive models that we constructed may be easier to use in clinical practice. However, this study also has some limitations. Firstly, the molecular mechanisms of how urea cycle-associated genes in prognostic models affect the biological behaviour of colorectal cancer cells require further experimental validation. In addition, although our prognostic model was confirmed using an external dataset, further well-designed prospective studies are needed to validate these findings. In future studies, we aim to focus on the roles of these UCRGs in the prognosis of patients with CRC.
Conclusions
In conclusion, we developed and validated a novel prognostic model based on the characteristics of eight UC-related genes in CRC for the first time, which will be of great significance in identifying prognostic molecular biomarkers, clinical prognosis prediction, and treatment strategy decisions for patients with CRC.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The studies involving human participants were reviewed and approved by Affiliated Hospital of North Sichuan Medical College (Nanchong, China). The patients/participants provided their written informed consent to participate in this study.
Author contributions
FW: Research design and manuscript revision. HG: Data analysis and QRT-PCR YW: Drafting of the article. LG, YT, and XW: Collection of clinical samples. All authors contributed to the article and approved the submitted version.
Funding
This study was supported by the Sichuan Provincial Primary Health Development Research Center (China) (grant number: swfz20-z-003).
Acknowledgments
We would like to thank the GEO and TCGA developers.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fsurg.2022.1027655/full#supplementary-material.
References
1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. (2021) 71(3):209–49. doi: 10.3322/caac.21660
2. Keum N, Giovannucci E. Global burden of colorectal cancer: emerging trends, risk factors and prevention strategies. Nat Rev Gastroenterol Hepatol. (2019) 16(12):713–32. doi: 10.1038/s41575-019-0189-8
3. Hubbard J, Grothey A. When less is more: maintenance therapy in colorectal cancer. Lancet (London, England). (2015) 385(9980):1808–10. doi: 10.1016/s0140-6736(14)62350-3
4. Pagès F, Mlecnik B, Marliot F, Bindea G, Ou F, Bifulco C, et al. International validation of the consensus immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet (London, England). (2018) 391(10135):2128–39. doi: 10.1016/s0140-6736(18)30789-x
5. Johnson C, Golla J, Dioletis E, Singh S, Ishii M, Charkoftaki G, et al. Molecular mechanisms of alcohol-induced colorectal carcinogenesis. Cancers (Basel). (2021) 13(17)4404. doi: 10.3390/cancers13174404
6. Roerink S, Sasaki N, Lee-Six H, Young M, Alexandrov L, Behjati S, et al. Intra-tumour diversification in colorectal cancer at the single-cell level. Nature. (2018) 556(7702):457–62. doi: 10.1038/s41586-018-0024-3
7. Hajaj E, Sciacovelli M, Frezza C, Erez A The context-specific roles of urea cycle enzymes in tumorigenesis. Mol Cell. (2021) 81(18):3749–59. doi: 10.1016/j.molcel.2021.08.005
8. Lee J, Adler L, Karathia H, Carmel N, Rabinovich S, Auslander N, et al. Urea cycle dysregulation generates clinically relevant genomic and biochemical signatures. Cell. (2018) 174(6):1559–70.e22. doi: 10.1016/j.cell.2018.07.019
9. Huang D, Li T, Wang L, Zhang L, Yan R, Li K, Xing S, Wu G, Hu L, Jia W, Lin SC, et al. Hepatocellular carcinoma redirects to ketolysis for progression under nutrition deprivation stress. Cell Res. (2016) 26(10):1112–30. doi: 10.1038/cr.2016.109
10. Boroughs L, DeBerardinis R. Metabolic pathways promoting cancer cell survival and growth. Nat Cell Biol. (2015) 17(4):351–9. doi: 10.1038/ncb3124
11. Li L, Mao Y, Zhao L, Li L, Wu J, Zhao M, et al. P53 regulation of ammonia metabolism through urea cycle controls polyamine biosynthesis. Nature. (2019) 567(7747):253–6. doi: 10.1038/s41586-019-0996-7
12. Mavri-Damelin D, Eaton S, Damelin LH, Rees M, Hodgson HJ, Selden C. Ornithine transcarbamylase and arginase I deficiency are responsible for diminished urea cycle function in the human hepatoblastoma cell line Hepg2. Int J Biochem Cell Biol. (2007) 39(3):555–64. doi: 10.1016/j.biocel.2006.10.007
13. De Santo C, Booth S, Vardon A, Cousins A, Tubb V, Perry T, Noyvert B, Beggs A, Ng M, Halsey C, Kearns P The arginine metabolome in acute lymphoblastic leukemia can be targeted by the pegylated-recombinant arginase I Bct-100. Int J Cancer. (2018) 142(7):1490–502. doi: 10.1002/ijc.31170
14. Kim J, Hu Z, Cai L, Li K, Choi E, Faubert B, et al., Cps1 maintains pyrimidine pools and DNA synthesis in Kras/Lkb1-mutant lung cancer cells. Nature. (2017) 546(7656):168–72. doi: 10.1038/nature22359
15. Doubleday P, Fornelli L, Ntai I, Kelleher N. Oncogenic kras creates an aspartate metabolism signature in colorectal cancer cells. FEBS J. (2021) 288(23):6683–99. doi: 10.1111/febs.16111
16. Ritchie M, Phipson B, Wu D, Hu Y, Law C, Shi W, et al. Limma powers differential expression analyses for Rna-sequencing and microarray studies. Nucleic Acids Res. (2015) 43(7):e47. doi: 10.1093/nar/gkv007
17. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. Clusterprofiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Cambridge (Mass). (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141
18. Newman A, Liu C, Green M, Gentles A, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12(5):453–7. doi: 10.1038/nmeth.3337
19. Livak K, Schmittgen T. Analysis of relative gene expression data using real-time quantitative Pcr and the 2[-Delta Delta C(T)] method. Methods (San Diego, Calif). (2001) 25(4):402–8. doi: 10.1006/meth.2001.1262
20. Vander Heiden M. Targeting cancer metabolism: a therapeutic window opens. Nat Rev Drug Discovery. (2011) 10(9):671–84. doi: 10.1038/nrd3504
21. Tian Z, Ou G, Su M, Li R, Pan L, Lin X, et al. Timp1 derived from pancreatic cancer cells stimulates schwann cells and promotes the occurrence of perineural invasion. Cancer Lett. (2022) 546:215863. doi: 10.1016/j.canlet.2022.215863
22. Ma B, Ueda H, Okamoto K, Bando M, Fujimoto S, Okada Y, et al. Timp1 promotes cell proliferation and invasion capability of right-sided colon cancers via the Fak/Akt signaling pathway. Cancer Sci. (2022) 00:1–14. doi: 10.1111/cas.15567
23. Duch P, Díaz-Valdivia N, Ikemori R, Gabasa M, Radisky ES, Arshakyan M, et al. Aberrant timp-1 overexpression in tumor-associated fibroblasts drives tumor progression through Cd63 in lung adenocarcinoma. Matrix Biol. (2022) 111:207–25. doi: 10.1016/j.matbio.2022.06.009
24. Liu H, Xiang Y, Zong QB, Zhang XY, Wang ZW, Fang SQ, et al. Mir-6745-Timp1 axis inhibits cell growth and metastasis in gastric cancer. Aging (Albany NY). (2021) 13(21):24402–16. doi: 10.18632/aging.203688
25. Vočka M, Langer D, Fryba V, Petrtyl J, Hanus T, Kalousova M, et al. Serum levels of Timp-1 and Mmp-7 as potential biomarkers in patients with metastatic colorectal cancer. Int J Biol Markers. (2019) 34(3):292–301. doi: 10.1177/1724600819866202
26. Nieman KM, Kenny HA, Penicka CV, Ladanyi A, Buell-Gutbrod R, Zillhardt MR, et al. Adipocytes promote ovarian cancer metastasis and provide energy for rapid tumor growth. Nat Med. (2011) 17(11):1498–503. doi: 10.1038/nm.2492
27. Jin J, Zhang Z, Zhang S, Chen X, Chen Z, Hu P, et al. Fatty acid binding protein 4 promotes epithelial-mesenchymal transition in cervical squamous cell carcinoma through Akt/Gsk3β/snail signaling pathway. Mol Cell Endocrinol. (2018) 461:155–64. doi: 10.1016/j.mce.2017.09.005
28. Guaita-Esteruelas S, Bosquet A, Saavedra P, Gumà J, Girona J, Lam EW, et al. Exogenous Fabp4 increases breast cancer cell proliferation and activates the expression of fatty acid transport proteins. Mol Carcinog. (2017) 56(1):208–17. doi: 10.1002/mc.22485
29. Uehara H, Takahashi T, Oha M, Ogawa H, Izumi K. Exogenous fatty acid binding protein 4 promotes human prostate cancer cell progression. Int J Cancer. (2014) 135(11):2558–68. doi: 10.1002/ijc.28903
30. Lee D, Wada K, Taniguchi Y, Al-Shareef H, Masuda T, Usami Y, et al. Expression of fatty acid binding protein 4 is involved in the cell growth of oral squamous cell carcinoma. Oncol Rep. (2014) 31(3):1116–20. doi: 10.3892/or.2014.2975
31. Cai JW, Huang XM, Li XL, Qin S, Rong YM, Chen X, et al. An 11-gene signature for the prediction of systemic recurrences in colon adenocarcinoma. Gastroenterol Rep (Oxf). (2021) 9(5):451–60. doi: 10.1093/gastro/goab023
32. Zhao D, Ma Y, Li X, Lu X. Microrna-211 promotes invasion and migration of colorectal cancer cells by targeting Fabp4 via pparγ. J Cell Physiol. (2019) 234:15429–37. doi: 10.1002/jcp.28190
33. Arumugam T, Logsdon C. S100p: a novel therapeutic target for cancer. Amino Acids. (2011) 41(4):893–9. doi: 10.1007/s00726-010-0496-4
34. Zhang Y, Du W, Chen Z, Xiang C. Upregulation of Pd-L1 by Spp1 mediates macrophage polarization and facilitates immune Escape in lung adenocarcinoma. Exp Cell Res. (2017) 359(2):449–57. doi: 10.1016/j.yexcr.2017.08.028
35. Liu L, Zhang R, Deng J, Dai X, Zhu X, Fu Q, et al. Construction of Tme and identification of crosstalk between malignant cells and macrophages by Spp1 in hepatocellular carcinoma. Cancer Immunol Immunother: CII. (2022) 71(1):121–36. doi: 10.1007/s00262-021-02967-8
36. Szulzewsky F, Pelz A, Feng X, Synowitz M, Markovic D, Langmann T, et al. Glioma-associated microglia/macrophages display an expression profile different from M1 and M2 polarization and highly express Gpnmb and Spp1. PloS one. (2015) 10(2):e0116644. doi: 10.1371/journal.pone.0116644
37. Choe E, Yi J, Chai Y, Park K. Upregulation of the adipokine genes Adipor1 and Spp1 is related to poor survival outcomes in colorectal cancer. J Surg Oncol. (2018) 117(8):1833–40. doi: 10.1002/jso.25078
38. Ji Y, Li J, Li P, Wang L, Yang H, Jiang G. C/Ebpβ promotion of Mmp3-dependent tumor cell invasion and association with metastasis in colorectal cancer. Genet Test Mol Biomark. (2018) 22(1):5–10. doi: 10.1089/gtmb.2017.0113
39. Shoshan E, Braeuer RR, Kamiya T, Mobley AK, Huang L, Vasquez ME, et al. Nfat1 directly regulates Il8 and Mmp3 to promote melanoma tumor growth and metastasis. Cancer Res. (2016) 76(11):3145–55. doi: 10.1158/0008-5472.Can-15-2511
40. Chu C, Liu X, Bai X, Zhao T, Wang M, Xu R, et al. Mir-519d suppresses breast cancer tumorigenesis and metastasis via targeting Mmp3. Int J Biol Sci. (2018) 14(2):228–36. doi: 10.7150/ijbs.22849
41. Wang S, Zhang C, Zhang Z, Qian W, Sun Y, Ji B, et al. Transcriptome analysis in primary colorectal cancer tissues from patients with and without liver metastases using next-generation sequencing. Cancer Med. (2017) 6(8):1976–87. doi: 10.1002/cam4.1147
42. Kamal Y, Schmit SL, Hoehn HJ, Amos CI, Frost HR. Transcriptomic differences between primary colorectal adenocarcinomas and distant metastases reveal metastatic colorectal cancer subtypes. Cancer Res. (2019) 79(16):4227–41. doi: 10.1158/0008-5472.Can-18-3945
43. Liang Y, Lv Z, Huang G, Qin J, Li H, Nong F, et al. Prognostic significance of abnormal matrix collagen remodeling in colorectal cancer based on histologic and bioinformatics analysis. Oncol Rep. (2020) 44(4):1671–85. doi: 10.3892/or.2020.7729
44. Sunami E, Tsuno N, Osada T, Saito S, Kitayama J, Tomozawa S, et al. Mmp-1 is a prognostic marker for hematogenous metastasis of colorectal cancer. Oncologist. (2000) 5(2):108–14. doi: 10.1634/theoncologist.5-2-108
45. Bendardaf R, Buhmeida A, Ristamäki R, Syrjänen K, Pyrhönen S. Mmp-1 (collagenase-1) expression in primary colorectal cancer and its metastases. Scand J Gastroenterol. (2007) 42(12):1473–8. doi: 10.1080/00365520701485449
46. Said AH, Raufman JP, Xie G. The role of matrix metalloproteinases in colorectal cancer. Cancers (Basel). (2014) 6(1):366–75. doi: 10.3390/cancers6010366
47. Kluz P, Kolb R, Xie Q, Borcherding N, Liu Q, Luo Y, et al. Cancer cell-intrinsic function of Cd177 in attenuating Β-catenin signaling. Oncogene. (2020) 39(14):2877–89. doi: 10.1038/s41388-020-1203-x
48. Shangkuan W, Lin H, Chang Y, Jian C, Fan H, Chen K, et al. Risk analysis of colorectal cancer incidence by gene expression analysis. PeerJ. (2017) 5:e3003. doi: 10.7717/peerj.3003
49. Zhang C, Wang H, Chen Z, Zhuang L, Xu L, Ning Z, et al. Carbonic anhydrase 2 inhibits epithelial-mesenchymal transition and metastasis in hepatocellular carcinoma. Carcinogenesis. (2018) 39(4):562–70. doi: 10.1093/carcin/bgx148
50. Nakada N, Mikami T, Horie K, Nagashio R, Sakurai Y, Sanoyama I, et al. Expression of Ca2 and Ca9 carbonic anhydrases in ulcerative colitis and ulcerative colitis-associated colorectal cancer. Pathol Int. (2020) 70(8):523–32. doi: 10.1111/pin.12949
51. Lin F, Zhang P, Zuo Z, Wang F, Bi R, Shang W, et al. Thioredoxin-1 promotes colorectal cancer invasion and metastasis through crosstalk with S100p. Cancer Lett. (2017) 401:1–10. doi: 10.1016/j.canlet.2017.04.036
52. Dong L, Wang F, Yin X, Chen L, Li G, Lin F, et al. Overexpression of S100p promotes colorectal cancer metastasis and decreases chemosensitivity to 5-Fu in vitro. Mol Cell Biochem. (2014) 389:257–64. doi: 10.1007/s11010-013-1947-5
53. Kikuchi K, McNamara K, Miki Y, Iwabuchi E, Kanai A, Miyashita M, et al. S100p and ezrin promote trans-endothelial migration of triple negative breast cancer cells. Cellular Oncology (Dordrecht). (2019) 42(1):67–80. doi: 10.1007/s13402-018-0408-2
54. Camara R, Ogbeni D, Gerstmann L, Ostovar M, Hurer E, Scott M, et al. Discovery of novel small molecule inhibitors of S100p with in vitro anti-metastatic effects on pancreatic cancer cells. Eur J Med Chem. (2020) 203:112621. doi: 10.1016/j.ejmech.2020.112621
55. Wu Z, Boonmars T, Nagano I, Boonjaraspinyo S, Srinontong P, Ratasuwan P, et al. Significance of S100p as a biomarker in diagnosis, prognosis and therapy of opisthorchiasis-associated cholangiocarcinoma. Int J Cancer. (2016) 138(2):396–408. doi: 10.1002/ijc.29721
56. Hsu Y, Hung J, Liang Y, Lin Y, Tsai M, Chou S, et al. S100p interacts with integrin Α7 and increases cancer cell migration and invasion in lung cancer. Oncotarget. (2015) 6(30):29585–98. doi: 10.18632/oncotarget.4987
57. Wang X, Tian T, Li X, Zhao M, Lou Y, Qian J, et al. High expression of S100p is associated with unfavorable prognosis and tumor progression in patients with epithelial ovarian cancer. Am J Cancer Res. (2015) 5(8):2409–21.26396916
58. Tzanakakis G, Nikitovic D. Preface of the special issue on the role of extracellular matrix in development and cancer progression. Biomolecules. (2022) 12(3):362. doi: 10.3390/biom12030362
59. Sullivan L. Metabolic frugality marks cancer cells for immune targeting. Cell. (2018) 174(6):1344–6. doi: 10.1016/j.cell.2018.08.023
Keywords: colorectal cancer, urea cycle, Gene signature, prognosis, immune infiltration
Citation: Guo H, Wang Y, Gou L, Wang X, Tang Y and Wang X (2022) A novel prognostic model based on urea cycle-related gene signature for colorectal cancer. Front. Surg. 9:1027655. doi: 10.3389/fsurg.2022.1027655
Received: 25 August 2022; Accepted: 4 October 2022;
Published: 21 October 2022.
Edited by:
Veronica De Simone, Agostino Gemelli University Polyclinic (IRCCS), ItalyReviewed by:
Yuan Yang, First Hospital of Lanzhou University, ChinaGrazisa Scialandrone, Azienda Sanitaria Localedella Provincia di Barletta Andri Trani (ASL BT), Italy
© 2022 Guo, Wang, Gou, Wang, Tang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Xianfei Wang Mjc1MDg1MzQ1OEBxcS5jb20=
†These authors have contributed equally to this work and share first authorship
Specialty Section: This article was submitted to Surgical Oncology, a section of the journal Frontiers in Surgery