Multi-Omics-Based Analysis of the Effect of Longevity Genes on the Immune Relevance of Colorectal Cancer

TL;DR

This study explores how longevity genes affect colorectal cancer and the immune system, offering new insights for treatment.

Contribution

A new prognostic model using five longevity genes for predicting colorectal cancer outcomes and immune status.

Findings

Five longevity genes were identified as an independent prognostic factor for colorectal cancer.

Risk scores correlated with immune cell infiltration and microsatellite instability in patients.

The model revealed distinct molecular subtypes linked to clinicopathological features and tumor microenvironment.

Abstract

Background: Colorectal cancer (CRC) ranks as the third most prevalent cancer globally, with its incidence and recurrence rates steadily rising. To explore the relationship between CRC and longevity-associated genes (LAGs), and to offer new therapeutic avenues for CRC treatment, we developed a prognostic model based on these genes to predict the outcomes for CRC patients. Additionally, we conducted an immune correlation analysis. Methods: We conducted a comprehensive analysis of the effects of 81 LAGs in CRC by integrating multiple omics datasets. This analysis led to the identification of two distinct molecular subtypes and revealed that alterations in LAGs across various layers were linked to clinicopathological features, prognosis, and cell infiltration characteristics within the tumor microenvironment (TME). The training and validation cohorts for the models were derived from the TCGA-COAD, TCGA-READ, and GSE35279 datasets. Subsequently, we developed a risk score model, and the Kaplan–Meier method was employed to estimate overall survival (OS). Ultimately, we established a prognostic model based on five LAGs: BEDN3, EXOC3L2, CDKN2A, IL-13, and CAPN9. Furthermore, we assessed the correlations between the risk score and factors such as immune cell infiltration, microsatellite instability, and the stem cell index. Results: Our comprehensive bioinformatics analysis revealed a strong association between longevity genes and CRC. The risk score derived from the five newly identified LAGs was determined to be an independent prognostic factor for CRC. Patients categorized by this risk score demonstrated significant differences in immune status and microsatellite instability. Conclusions: Our comprehensive multi-omic analysis of LAGs highlighted their potential roles in the tumor immune microenvironment, clinicopathological features, and prognosis, offering new insights for the treatment of CRC.

Figures

Full text

1. Introduction

Colorectal cancer (CRC) is the third most common cancer and the second leading cause of cancer death worldwide [1]. Nearly 30% of new CRC cases and deaths worldwide each year occur in China. According to a previous study, the number of newly diagnosed CRC cases in 2022 will exceed 510,000 [2]. In addition to previous common treatments, such as surgery, radiotherapy and chemotherapy, in recent years, chimeric antigen receptor (CAR) T-cell treatment [3], photodynamic therapy [4], and nanotechnology-based drug/protein delivery systems [5] have gradually become popular and have shown significant efficacy in antitumor treatment. However, the prevention, diagnosis and treatment of cancer have always been major challenges for scientific research and clinical staff, and improving the survival and health status of elderly patients is still challenging.

As the highest primate of humans, the lifespan of C. elegans is no more than tens of days, but humans can live for one hundred years. The longevity genes (anti-aging genes) provide the answer that we want [6]. However, cancer risk usually increases with age, which may be related to long-lived organisms having a longer time to accumulate mutations [7]. A longer lifespan is associated with the evolution of more stringent cell cycle control mechanisms, which increase the sensitivity of cells to growth conditions and affect the tumor suppressor effect [8]. A recent univariate analysis for model construction revealed that age, sex, BMI, family history of CRC, smoking status, etc., were independently associated with CRC risk [9].

Immune senescence is a series of age-related changes that affect the immune system and lead to increased susceptibility to disease over time [10]. As a complex cell network, the tumor microenvironment (TME) is an important factor in regulating the progression of CRC [11]. Immune cells in the TME (cancer-associated fibroblasts, cancer-associated adipocytes, myeloid-derived suppressor cells, etc.) play a role in the development and progression of CRC by producing various cytokines and chemokines and interacting with many chronic inflammatory mediators [12]. The results of immunohistochemistry for mismatch repair proteins (MLH1, PMS2, MSH2, and MSH6) and Polymerase Chain Reaction (PCR) for microsatellite instability have prognostic and therapeutic significance in patients with sporadic CRC [13] and have now become one of the standards guiding clinical medication, but more exploration is needed on individualized treatment options for different types of CRC patients.

This study systematically investigated the characteristics of immune infiltration mediated by longevity-associated genes (LAGs) within the TME of CRC using bioinformatics approaches, including consensus clustering. This work bridges the gap between aging biology and cancer immunology, offering translational implications for precision immunotherapy in CRC and other age-associated malignancies.

2. Methods

2.1. Data Download

In this study, we learned from the official The Cancer Genome Atlas (TCGA) website (https://portal.gdc.cancer.gov/), which integrates RNA sequencing, genomic mutation data and clinical data (including age, sex and survival status) from the TCGA-COAD and TCGA-READ cohorts. The Gene Expression Omnibus database (GEO) (https://www.ncbi.nlm.nih.gov) provides mRNA sequencing data and clinical data (GSE35279) from 74 CRC patients, and this dataset underwent quantile normalization and log2 transformation before analysis. In addition, we extracted the copy number variation data of CRC patients from the Xena UCSC database (https://genome.ucsc.edu/). All public data was accessed on 13 September 2024. This study was performed in accordance with the Declaration of Helsinki (Revised 2013).

2.2. Comprehensive Analysis of Longevity-Related Genes in the TCGA Database

On the basis of a meta-analysis of individuals from 20 cohorts of European, East Asian, or African American descent, gene-level association analysis of tissue-specific gene expression, and a Genome-Wide Association Study (GWAS)of the strict definition of the longevity phenotype, Joris Deelen [14] identified LAGs. The RNA-seq expression data from the TCGA database were first normalized to fragments with a value of millions of bases per thousand, and the RNA-seq expression data from the R language (version: 4.4.1) were first normalized to fragments with a value of millions of bases per thousand [15]. The “limma” package (version 3.56.2) was used to identify differentially expressed genes (DEGs) between normal tissues and tumor tissues; DEGs with screening p values < 0.05 are shown in a box plot. The significance level of the DEGs is expressed as follows: * indicates a p value < 0.05, ** indicates a p value < 0.01, and *** indicates a p value < 0.001. Next, the “survival” package of R and the log-rank test [16] were used. Kaplan-Meier survival analysis of patients with 34 differentially expressed genes at the intersection of COAD and READ was performed, and the purpose was to compare the survival differences among multiple groups. The mutation data associated with COAD and READ patients in the TCGA dataset were downloaded and processed separately, and the “maftools” R [17] package was used to illustrate and visualize the most frequently mutated genes. The “Rcircos” package [16] (version 1.2.2) is helpful for describing the changes in the copy numbers of longevity-related genes on human chromosomes.

2.3. Consensus Clustering Classification

First, RNA expression data and clinical data from the TCGA and GEO databases were integrated using Perl scripts. Next, the “RColorBrewer” package in R was utilized to generate a prognostic network plot of LAGs. Subsequently, the R package “ConsensusClusterPlus” (version 1.64.0) was applied to the data of all 689 CRC patients from the TCGA and GEO databases to investigate the relationships between significantly expressed LAG and CRC subtypes. The optimal K value, representing the clustering variable with the highest correlation within clusters and the lowest correlation between clusters, was selected for classification, and a survival curve was generated. We combined the gene expression and clinical characteristics of all CRC patients in the TCGA and GEO databases, including the tumor stage, sex, age, and sample classification, and visualized these data using a heatmap created with the “pheatmap” software package in R (version 1.0.12). The “limma” package was employed to conduct intercluster variance analysis and principal component analysis (PCA). Based on the PCA results, we identified the common intersection of interacting genes classified as A, C, E, and F to obtain the clustering intersection genes. To facilitate functional enrichment analysis, Weishengxin’s strongest algorithm was used to integrate the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases, allowing for more effective clustering of intersection genes. A significance threshold of p < 0.05 was applied, and results were displayed in histograms and network diagrams. Additionally, the “limma” and “ggpubr” packages in R were used to conduct univariate Cox regression analysis on the cluster intersection genes, enabling the screening and visualization of prognosis-associated genes in a forest plot. Finally, we retyped and screened the samples based on the prognostic genes, selected the most appropriate K value for sample typing, identified DEGs, and generated relevant survival curves.

2.4. Construction of the Prognostic Model

After the CRC prognostic risk genes were incorporated into the model, we used a minimum absolute shrinkage and selection operator (LASSO) [18]. The algorithm was employed for dimensionality reduction to identify prognostic feature genes. The risk score was calculated using the following formula: (expression of each signature gene * model coefficient value). Patients were subsequently divided into high-risk and low-risk groups based on the calculated average risk score. The “ggalluvial” package in R was utilized to construct the prognostic models and create the Sankey diagram. To analyze the differences in risk scores between the initial clustering and gene clustering, the “limma” package was employed, and the results were visualized using box plots generated with the “ggpubr” package. Patients were classified into high-risk or low-risk categories according to the median prognostic risk score. The CRC patients were evenly split into a training group and a test group. Kaplan-Meier survival curves were generated using the “survival” package in R to assess significant differences in overall survival (OS) between these groups. The effectiveness and accuracy of the models were evaluated by measuring the area under the curve (AUC) for 1, 3, and 5 years. The risk curve function was defined in R, and risk curves, survival state maps, and risk heatmaps for both the training and test groups were generated. Additionally, a nomogram model for clinical prediction was constructed by integrating the risk score model with key clinicopathological parameters (sex, age group, tumor stage, risk value, total score, and linear predictive value). “time ROC” [19] package analysis was used to evaluate the performance and accuracy of the predictive features. In addition, the R “include” function was used to draw mutational waterfall plots of the low-risk and high-risk groups.

2.5. Immune, Mismatch Repair and Stem Cell Index Analysis

Analysis of immune cell infiltration in tissues is a key part of disease characterization and prediction of disease prognosis. Using the CIBERSORT algorithm [20], the relative content of each immune cell in the immune microenvironment of high-risk and low-risk samples was evaluated. All immune cells were then analyzed to construct a scatter plot and a correlation heatmap illustrating the relationship between the risk score and immune cell types. The “estimate” package in R was employed to assess the abundance of stromal and immune cells in each sample, while the “ggplot2” package was used to create violin plots for visualization. The results of the differential analysis among the microsatellite stability group, low-frequency microsatellite instability group, and high-frequency microsatellite instability group were presented as percentage histograms and box plots. Furthermore, the R package “ggExtra” was utilized to calculate the correlation between the risk score and the stem cell index score.

2.6. Statistical Analysis

To standardize the processing of transcriptomic, clinical, and gene expression data in this study, PERL software (version 5.30.0) was used. Continuous variables were analyzed using Student’s t-test (or Welch’s t-test when homogeneity of variance was violated), while categorical variables were compared with Pearson’s chi-square test (or Likelihood Ratio chi-square test when assumptions were unmet). All the statistical analyses were performed via R version 4.4.1, and the significance level was set at p < 0.05.

3. Results

3.1. Correlation Analysis of Longevity-Related Genes in the TCGA Database

The longevity-related differential expression genes (DEGs) between colon cancer, rectal cancer, and normal tissues were analyzed, and the intersection of the differential genes between the two groups was obtained by Venn diagram to obtain 36 common differentially expressed longevity-related genes (Figure 1A). Our analysis revealed that 44 longevity-related genes were upregulated in colon cancer, while 19 genes were upregulated in rectal cancer. Survival analysis of differentially expressed longevity-related genes in colorectal cancer found that the expression of 18 genes, including RAD50, ELOVL6, and BEND3, had a greater impact on patient survival (Figure S1). Copy number variations (CNVs) were observed across multiple chromosomes; however, mutations associated with CRC were predominantly located on chromosome one (Figure 1B). Additionally, CNV analysis indicated a significant increase in the proportion of CNVs for CHRNA4 and USP42 in colon cancer, while an increase in CNVs for FURIN and FES was noted in rectal cancer. Conversely, a significant loss of CNVs for LPL was observed in CRC (Figure 1C). The mutation frequencies of these genes were analyzed using a heatmap generated from the TCGA dataset. IGF2R exhibited a high mutation frequency in both colon and rectal cancers, with the majority of mutations being missense mutations (Figure 1D). Furthermore, the network diagram results confirmed the existence of interactions and tight regulatory connections among longevity-related genes (Figure 2A).

3.2. Construction of the Clustering Model

To gain a deeper understanding of the expression characteristics of DEGs in CRC, we employed a consensus clustering method. Initially, we integrated data from patients in the TCGA and GEO databases and conducted consistent cluster analysis. When the clustering variable (K = 6), we observed the highest homogeneity within groups and the lowest correlation between groups, resulting in the division of the 689 CRC patients into six clusters, each exhibiting significant differences in survival (Figure 2B,C). PCA indicated a low correlation among Clusters A, C, E, and F (Figure S2A). A comparison of immune cell infiltration across the six clusters revealed varying levels of activated immune T-cell enrichment among the clusters (Figure S2B). The heatmap further illustrates the clinical characteristics of patients within the different clusters (Figure S2C). This indicates that subtypes classified based on these 36 DEGs can accurately distinguish CRC patients and provide important evidence for personalized treatment of CRC patients. Additionally, we utilized the Venn diagram function to intersect the differentially expressed genes in patients from Clusters A, C, E, and F, resulting in 97 co-expressed clustered differential genes (Figure 2D). GO and KEGG functional enrichment analyses revealed that these genes were primarily associated with neutrophil degranulation, leukocyte chemotaxis, proliferation, differentiation and intercellular adhesion, T-cell activation, leukocyte transendothelial migration, chemokine signaling pathways, and the PI3K-Akt signaling pathway (Figure 2E,F). The clinical significance of the 81 longevity-related genes was further analyzed through univariate Cox regression, identifying 10 genes with notable prognostic value. Specifically, IL-13, ELOVL6, BEND3, PSRC1, ZW10, and CAPN9 were recognized as protective factors, while CDKN2A, FES, EXOC3L2, and MAPK4 were identified as risk factors (Figure 3A). Based on these significant prognostic genes, a second gene cluster analysis was conducted, resulting in the classification of samples into two groups, A and B (Figure 3B). Survival curve analysis demonstrated a significant difference between the two groups (p = 0.006), with patients in group A exhibiting higher survival scores than those in group B (Figure 3C). We also analyzed the expression differences of LAGs between the two genotypes, revealing statistically significant differences in the expression of 38 LAGs between type A and type B patients (Figure 3D). A Sankey diagram illustrates the classification process (Figure 4A). The heatmap presents gene expression profiles along with clinical characteristics (tumor stage, age, etc.), highlighting a significant age difference among patients across the different groups (Figure S2D). From the prognostic genes, five were selected through LASSO regression analysis (Figure S2E,F), leading to the formulation of the risk score calculation = ((IL-13 * −1.82) + (CDKN2A * 0.18) + (BEND3 * −0.50) + (EXOC3L2 * 0.49) + (CAPN9 * −0.12)). This formula was utilized to predict risk scores for samples in the training group, which were subsequently divided into high-risk and low-risk categories based on the median risk score. To evaluate the robustness of the risk score, we conducted a differential analysis of risk scores between the initial and second genotyping groups, confirming reliable results (Figure 4B,C). Based on the calculated mean risk score, we analyzed the differences in LAGs between the high- and low-risk groups, identifying 24 LAGs with significant expression differences (Figure 4D). Furthermore, we examined the mutational landscape of the five prognostic genes in both risk groups, noting that BEND3 exhibited the highest frequency of missense mutations (Figure S2G).

3.3. Construction and Validation of the Clinical Prediction Models

Based on the prediction model derived from the risk score, we performed univariate and multivariate Cox regression analyses on clinical and pathological features, including patients’ sex, age, and tumor stage, to construct a clinical prediction nomogram model (Figure 4E). Encouragingly, the index curve confirmed the model’s accuracy in predicting survival rates, demonstrating superior performance compared to clinical factors such as age, sex, and stage. This indicates that the model is highly accurate in predicting the 1-, 3-, and 5-year survival rates of patients (Figure 4F). To assess the reliability of the risk model, we conducted a differential analysis of survival between high-risk and low-risk groups. Survival curves indicated that patients in the low-risk group exhibited better survival scores compared to those in the high-risk group, both in the training and test groups, as well as in the overall sample (Figure 5A). The area under the curve (AUC) values for the 1-, 3-, and 5-year survival rates were all greater than 0.5, further confirming the robustness of the results (Figure 5B). The risk heatmap illustrated that the expression of high-risk genes, specifically CDKN2A and EXOC3L2, increased alongside rising risk scores (Figure 5C). Additionally, we categorized patients into high-risk and low-risk groups based on risk scores from the training, test, and overall samples (Figure 5D). An increase in risk score corresponded to a heightened risk of mortality, with the probability of survival gradually decreasing as the duration of survival increased (Figure 5E). These results indicate that the prognostic model can effectively predict the prognosis of CRC patients.

3.4. Immune Cell and Mismatch Repair Assays

Our analysis of five prognostic genes revealed distinct immune correlations. CAPN9 expression exhibited a positive association with resting CD4^+^ memory T cells while demonstrating a significant negative correlation with M1 macrophages. Conversely, EXOC3L2 showed positive correlations with regulatory T cells (Tregs) and M0 macrophages but an inverse relationship with eosinophils. IL-3 was positively linked to eosinophil infiltration, whereas CDKN2A displayed a negative association with resting CD4^+^ memory T cells. Elevated risk scores were inversely correlated with infiltration levels of activated dendritic cells, mast cells, resting CD4^+^ memory T cells, neutrophils, and eosinophils (Figure 6B–H), suggesting that low-risk patients may exhibit enhanced immune activation, while high-risk groups display stromal activation patterns (Figure 6I). Notably, the low-risk group demonstrated a higher prevalence of microsatellite instability-high (MSI-H) status, with MSI-H patients consistently showing lower risk scores (Figure 6J,K). Intriguingly, CRC stemness indices were negatively correlated with risk stratification (Figure 6L).

4. Discussion

Human longevity is a phenotype that integrates many aspects of health and the environment into one ultimate quantity: the elapsed time from birth to death [21]. The study of the genetic components of human longevity can provide important insights into the prevention of age-related diseases and the mechanisms of various diseases. Genetic longevity research currently faces a major challenge, namely the lack of a deep understanding of the nature of the genetic component and its size, as well as the lack of optimal strategies for meta-analysis and subject selection in next-generation sequencing (NGS) efforts [22]. Age is an independent prognostic factor for a variety of malignant tumors, and the activation of longevity genes has been confirmed to have strong tumor-promoting effects. However, telomere gene-mediated progressive telomere shortening (replicative senescence) induced cell proliferation arrest is also one of the key anticancer mechanisms in humans [23]. This dual role suggests that longevity genes may be an important bridge between cancer and aging.

The Chinese guidelines for CRC diagnosis and treatment recommend that people after the age of 40 undergo colonoscopy for CRC screening every 5–10 years [24]. However, this will face an enormous cost and resource burden. Therefore, innovative CRC screening indicators, blood-based biomarker detection and stool DNA detection have shown potential for improving the screening participation rate and the sensitivity of precancerous lesions. On the other hand, in response to the heterogeneity of CRC patients, the development of individualized treatment regimens has become a hot topic of current research. Some clinical trials already include combinations of chemotherapy, radiotherapy, hormone therapy and immunotherapy. The combination of cytotoxic compounds and biological therapy (monoclonal or polyclonal antibodies, vaccines, gene therapy, etc.) is also increasingly important [25]. Therefore, the innovative prognostic model constructed in this study may provide new insights for precision immunotherapy.

As an extension of cancer diagnosis, cancer subtype analysis can be considered a consensus clustering problem. This analysis is conducive to providing more accurate treatment for patients [26]. Consistency cluster analysis has been used to combine clinical research data with data from TCGA. The principle is to use the dataset to discover new subtypes and identify biological processes specific to different subtypes [27]. Currently, its application is not limited to genomics research; the consensus clustering analysis model based on CT radiomics is also gradually gaining popularity [28]. Therefore, we applied this method to CRC to identify significant differences among subgroups. Based on multi-omics analysis, this study explored the role of longevity-related genes in CRC. By differential expression analysis, CRC patients were divided into six clusters for the first time, and longevity genes associated with prognosis were screened based on univariate Cox regression analysis. Subsequently, a second consensus cluster analysis of the CRC cohort, based on the expression patterns of these genes, identified two molecular clusters associated with longevity. The risk score was calculated by LASSO regression analysis, and the reliability of the two clusters was verified by the score. In addition, we evaluated patient outcomes using univariate and multifactorial Cox regression analyses combined with clinicopathological features and constructed clinical predictive models. We found that risk scores derived from the expression levels of five specific LAGs (BEDN3, EXOC3L2, CDKN2A, IL-13, and CAPN9) had significant prognostic significance. Survival analysis showed that patients in the high-risk group had a poorer prognosis. These results suggest that longevity-related genetic signatures may provide important guidance for clinical treatment decisions in CRC patients.

Interleukin 13 (IL-13) is an immune regulatory cytokine that is mainly secreted by Th2-like cells and natural killer T cells [29]. Among the CRC samples, IL-13 receptor overexpression was detected in more than 66% of the tumor samples. Moreover, in the immunohistochemistry analysis of stage I-III CRC patients, high IL-13 and IL-13R receptor expression was observed at 50% and 42%, respectively [30]. Saigusa studied 241 CRC patients and reported that serum IL-13 levels were significantly reduced in advanced patients and that low serum IL-13 levels were significantly correlated with poor prognosis [31]. Studies have reported that IL-13 overexpression promotes CRC liver metastasis by activating PI3K/AKT-related signaling pathways and that FAM120A in the IL-13/IL-13R signaling pathway is a key mediator of CRC invasion and liver metastasis [32]. IL-13R silencing inhibits the IL-13-induced proliferation of CRC cell lines through the downregulation of STAT6 activation [33]. These findings reveal the hub role of IL-13 in the occurrence and development of CRC, which provides a theoretical basis for targeting and immunotherapy in CRC patients. The tumor suppressor gene CDKN2A (p16) is a major target of carcinogenesis, with a frequency second only to that of the p53 tumor suppressor gene [34]. A recent study revealed that CDKN2A mRNA expression acts as a biomarker of clinically invasive meningiomas, with potential therapeutic significance [35]. A meta-analysis of CRC revealed a significant association between CDKN2A hypermethylation and lymphovascular invasion, lymph node metastasis, and proximal tumor location [36]. Another observational study revealed that the absence of the tumor suppressor CDKN2A was associated with the risk of CRC recurrence [37]. BEND3 is a four-member BEN domain-containing protein that binds to heterochromatin and functions as a transcription repressor [38]. The lncRNA XLOC_011677-miR-10b-5p-BEND3 potential exosomal ceRNA regulatory network may reveal regulatory pathways and diagnostic biomarkers related to the molecular biology [39] of CRC. EXOC3L2 was previously considered one of the susceptibility genes for Alzheimer’s disease, and its study of cancer may become a new direction for future research [40]. CAPN9 is expressed mainly in the stomach and small intestine and is involved in the regulation of hypertension, heart disease and gastric mucosal defense. It also potentially plays a role in the upregulation of metalloproteinases in the structural changes in the artery wall [41,42]. MiR-585-3p inhibits the proliferation and migration of high-grade serous ovarian cancer cells by targeting CAPN9 [43]. In our study, longevity-related genes were used for the first time to develop a prognostic model of CRC and played an important role in advancing tumor research.

In this study, we integrated data from a total of 689 CRC patients obtained from the TCGA and GEO databases, employing bioinformatics approaches, including multi-omics analyses, to investigate the effects of LAGs on CRC in detail. The construction of a clinically relevant prognostic model offers new avenues for research aimed at potentially improving the prognosis of CRC patients. We also explored the utility of the risk score in predicting responses to immunotherapy and analyzed how varying risk scores influence the differential expression of immune cells within tumors. Despite these advancements, our study has several limitations. First, there is an overreliance on publicly available database data. Second, the model requires validation through large, multicenter, prospective clinical trials to confirm its applicability. Finally, this study lacks experimental validation, and there are no reliable experimental data to support more comprehensive investigations.

5. Conclusions

In this study, we constructed a prognostic model based on five LAGs, which also elucidated the key interactions between LAGs and CRC-specific immune cells. This innovative prognostic model provides new insights for precision immunotherapy, which can be used to fight this disease and improve the prognosis of patients, thus opening a new way for personalized treatment selection.