Abstract

Purpose

   To explore the role of senescence-associated secretory phenotype (SASP)
   in ulcerative colitis-related colorectal cancer (UCRCC) precision
   diagnosis and pathogenesis.

Methods

   In this study, we first screened the SASP-related genes (SASPRGs) in
   the ulcerative colitis (UC) and colorectal cancer (CRC) datasets using
   differentially expressed genes (DEGs) and weighted gene co-expression
   network analysis (WGCNA). Subsequently, the core SASPRGs were further
   identified using 113 combinations of 10 machine learning algorithms,
   and evaluating their expression levels, diagnostic efficacy, prognostic
   value, and relationship with immune cells. Expression patterns of core
   SASPRGs were analyzed using scRNA-seq, and a molecular regulatory
   network was constructed. In addition, immunotherapy response and drug
   sensitivity analysis were conducted to screen therapeutic drugs
   targeting core SASPRGs. Finally, their expression levels were validated
   through in vivo experiments.

Results

   We identified 12 core SASPRGs, among which TGFBI and ANXA6 exhibited
   satisfactory diagnostic efficacy and prognostic value for UC and CRC.
   ANXA6 and TGFBI are highly expressed in monocytes and fibroblasts, and
   are associated with the infiltration of multiple immune cells. In
   addition, high TGFBI and ANXA6 groups showed poor responses to
   immunotherapy. The candidate therapeutic drugs screened by drug
   sensitivity analysis showed good binding ability with ANXA6 and TGFBI.
   Finally, the expression levels of TGFBI and ANXA6 were significantly
   increased in UC and CRC mouse models.

Conclusion

   Overall, TGFBI and ANXA6 are crucial in UCRCC, serving as novel
   diagnostic markers. They exhibit robust predictive capabilities on
   patient prognosis and immunotherapy response, offering actionable
   insights to optimize therapeutic decision-making and advance
   personalized treatment paradigms in UCRCC management.

   Keywords: ulcerative colitis, colorectal cancer, ulcerative
   colitis-related colorectal cancer, senescence-associated secretory
   phenotype, machine learning

Introduction

   Ulcerative colitis (UC) is a chronic, recurrent inflammatory bowel
   disease with an unclear etiology, that may be closely related to
   multiple factors such as genetic, environmental, microbial and immune
   factors.[32]^1–4 UC has become a global disease, and in the early 21st
   century, the prevalence of UC is increasing in some newly
   industrialized countries, which leads to a large increase in the burden
   of disease.[33]^5^,[34]^6 UC patients have a high risk of developing
   ulcerative colitis-related colorectal cancer (UCRCC), and a
   meta-analysis showed that the cumulative risk of UCRCC in UC patients
   with 30 years of onset was as high as 13.9%.[35]^7^,[36]^8 Although
   recent research reports indicate that the incidence of colorectal
   cancer (CRC) has been declining for a long time, the risk for patients
   with chronic ulcerative colitis remains higher compared to the general
   population.[37]^9 Furthermore, there are still no effective indicators
   for diagnosing UCRCC at present.

   Cellular senescence is a stable terminal state of cells characterized
   by the loss of proliferative and differentiation capabilities.[38]^10
   Cells can transition to a senescent state due to various stressors,
   such as oncogene-induced senescence, therapy-induced senescence,
   mitochondrial dysfunction-induced senescence, and immune-induced
   senescence.[39]^11–15 The senescence-associated secretory phenotype
   (SASP) is the primary mediator secreted by senescent cells in the
   tissue microenvironment, including chemokines, growth factors,
   cytokines, and stromal metalloproteinases.[40]^16 Increasing evidence
   suggests that senescent cells mainly promote disease progression by
   secreting SASP.[41]^16–18 In UC and CRC, various methods can induce
   cellular senescence, including therapy-induced senescence,
   oncogene-induced senescence, and immune-induced senescence. Among
   these, the senescence induced by changes in the immune microenvironment
   mediated by chronic inflammation has received considerable attention
   from researchers in recent years.[42]^19–22 It is reported that the
   SASP-related genes (SASPRGs) can promote the malignant progression of
   CRC.[43]^23^,[44]^24 However, due to the lack of research on the
   SASPRGs in UC and CRC, it is currently unclear which SASPRGs contribute
   to their progression, so identifying core SASPRGs may provide a
   theoretical basis for the diagnosis, prevention, and formulation of
   personalized treatment strategies for these diseases.

   In this study, we evaluated SASP-related differentially expressed genes
   in UC and CRC patients using bulk RNA transcriptomics. We determined
   the core SASPRGs using 113 combinations of 10 machine learning
   algorithms and identified their diagnostic potential and clinical
   value. Furthermore, we also screened potential therapeutic drugs for
   patients with high expression of TGFBI and ANXA6, and constructed a
   molecular regulatory network, providing new insights for more effective
   diagnostic and therapeutic strategies for UCRCC.

Materials and Methods

Datasets and Processing

   Gene expression matrices and related clinical information about UC and
   CRC patients were obtained from the TCGA
   ([45]https://portal.gdc.cancer.gov/) and GEO
   ([46]https://www.ncbi.nlm.nih.gov/geo/) databases. [47]Supplementary
   Table 1 presents the full information about the dataset used in this
   study.

Differential Expression and Functional Enrichment Analysis

   The R package “limma” was used to identify differential expression
   genes (DEGs) in the training datasets. Any gene with adjusted values of
   P < 0.05 and |log2 (FC)| of > 0.8 was regarded as DGEs. Gene Ontology
   (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Set
   Enrichment Analysis (GSEA) analyses were performed using the R package
   “clusterProfiler”. The visualization of this result was plotted by R
   package “ggplot2”.

Weight Gene Co-Expression Network Analysis (WGCNA)

   The weight gene co-expression network analysis (WGCNA) was performed
   using the R package “WGCNA”. Firstly, the hierarchical clustering
   method is used to filter the dataset, and the “pickSoftThreshould”
   function is used to determine the appropriate soft threshold β.
   Finally, the corresponding modules are generated. Among these modules,
   the one most relevant to the clinical phenotype was selected.

Identification of Core SASPRGs Associated with UC and CRC

   Human putative SASP factor genes were extracted from the Gene Ontology
   Consortium and QuickGO database.[48]^25 DEGs and key module genes
   identified by WGCNA were intersected with SASP to obtain SASPRGs
   associated with UC and CRC. To further identify core SASPRGs, we
   calculated the AUC values of 113 combinations of 10 machine learning
   algorithms based on these results, respectively. Ultimately, the genes
   included in the combination with the highest AUC value were considered
   core SASPRGs. The median was used to categorize the core immune-related
   genes into low and high expression groups, and then Kaplan–Meier
   analysis and Cox regression analysis were performed with the TCGA-COAD
   dataset downloaded from TCGA databases.

Immune Cell Infiltrates Analysis

   We evaluated the immune cell infiltration in tissue from UC and CRC
   patients by performing ssGSEA, CIBERSORT and MCPcounter analyses using
   the R package “GSVA”, “CIBERSORT”, and “MCPcounter”. Finally, we
   explored the relationship between the 28 immune cells and core SASPRGs
   using Spearman correlation analysis.

Single-Cell Analysis

   The [49]GSE162335 and [50]GSE231559 datasets were used for single-cell
   analysis. The R package “Seurat” was used for subsequent data
   processing. TSNE analysis was used to identify different clusters.
   Using the “BlueprintEncodeData” dataset from the R package “celldex” as
   a reference, each cluster was initially annotated with “SingleR”. The R
   package “UCell” is used to evaluate the expression of core SASPRGs in
   different cell types, the UCell score of each cell was mapped to the
   TSNE embedding with “ggplot2” for visualization. Bubble plots were used
   to represent the expression levels of core SASPRGs in different cell
   types. In addition, expression patterns of core SASPRGs in different
   patient tissues were validated at the single-cell level.

Prediction of the Immunotherapy Response and Drug Sensitivity Analysis

   CRC patients were divided into two groups according to the average
   expression level of SASPRGs, respectively. The immunotherapy response
   scores of different groups were calculated by TIDE
   ([51]http://tide.dfci.harvard.edu/login/), and the differences in TIDE
   scores between different groups were compared by Wilcoxon test. A
   higher TIDE score indicates a greater likelihood of tumor immune
   escape.

   Drug sensitivity data were obtained from GDSC
   ([52]https://www.cancerrxgene.org/) and CTRP
   ([53]https://portals.broadinstitute.org/ctrp/) databases, and then the
   R package “oncoPredict” was used to calculate the IC50 value of drugs
   in different subgroups of CRC patients. The Wilcoxon test was used to
   compare the difference of IC50 values between different groups.

   To further analyze the correlation between the identified drugs and
   SASP-related proteins, we first obtained the 3D structure of the drugs
   from the PubChem database ([54]https://pubchem.ncbi.nlm.nih.gov/) and
   the 3D structure of SASP-related proteins from the PDB database
   ([55]https://www.ebi.ac.uk/pdbe/) and then used the Autodock tool to
   predict the binding ability of the drugs to the proteins.

Construction of a Molecular Regulator Network

   To predict the regulatory mechanisms of ANXA6 and TGFBI, we constructed
   a molecular regulatory network using the NetworkAnalyst
   ([56]https://www.networkanalyst.ca/NetworkAnalyst/) platform. We
   firstly obtained gene-miRNA interaction information from miRTarBase and
   TF-gene interaction information from JASPAR. Finally, Cytoscape was
   used to visualize the molecular regulatory network.

Construction of the Animal Model

   The UC animal model was established according to the previous
   experimental methods. First, eight-week-old C57BL/6 male mice were
   acclimatized in a specific pathogen-free animal facility for 7 days.
   Normal drinking water was subsequently replaced with a 3% dextran
   sulfate sodium (DSS, meilunbio, MB5535) solution for 7 days to induce
   UC. The mice in the control group were provided with equal volumes of
   water throughout the course of the experiment. Mice were euthanized on
   day 8, and colon tissues were collected for further studies. Body
   weight, fecal occult blood, and fecal viscosity were recorded daily to
   assess disease activity index (DAI). The CRC animal model previously
   constructed by our team was used in this study. The specific methods
   were as follows: Step 1: Eight-week-old C57BL/6 male mice were placed
   in a specific pathogen-free animal facility. The mice in the
   experimental group were intraperitoneally injected with azomethane
   (AOM) at a dose of 10mg/kg and fed with normal drinking water for
   1 week. Step 2: Normal drinking water feeding was replaced with 2.5%
   DSS for one week. Step 3: Feeding with normal drinking water for
   2 weeks. Step 4: Repeat steps 2 and 3 three times. All experiments
   followed the guidelines approved by the Ethical Review Committee for
   Animal Experiments of Nanjing Medical University.

RNA Isolation and qRT-PCR

   Total RNA was extracted from colon tissues by Trizol reagent
   (Invitrogen) and used to synthesize cDNA with the PrimeScript^TM RT
   reagent Kit (TaKaRa, Dalian, China) according to the manufacturer’s
   instructions. qRT-PCR was performed on an ABI 7900HT PCR sequencer
   (Applied Biosystems, Massachusetts, USA) using TB GreenR Premix Ex Taq™
   II (TaKaRa). PCR sequencer (Applied Biosystems, Massachusetts, USA)
   using TB GreenR Premix Ex Taq™ II (TaKaRa). Fold changes in mRNA
   expression were calculated using 2−ΔΔCt and standardized based on
   GAPDH. Primer sequences used for qRT-PCR were as follows: ANXA6,
   forward: 5′-AGAGCTACAAGTCCCTCTACG-3′, reverse:
   5′-CCCACAATCAACCGTTCAAAC-3′. TGFBI, forward:
   5′-CACTCTCAAACCTTTACGAGACC-3′, reverse: 5′-CGTTGCTAGGGGCGAAGATG-3′.
   GAPDH, forward: 5′-TGCACCACCAACTGCTTAGC-3′, reverse:
   5′-GGCATGGACTGTGGTCATGAG-3′.

Histopathological Staining

   Mouse colon tissues were fixed overnight in universal tissue fixative
   solution (Servicebio, G1101). The fixed tissues were embedded in
   paraffin and sliced in a gradient concentration of alcohol after
   dehydration. After staining with hematoxylin and eosin, the slices were
   dehydrated and sealed; observations were made using an upright light
   microscope (NIKON ECLIPSE E100, Nikon, Japan) and imaging system (NIKON
   DS-U3, Nikon, Japan).

Immunostaining

   Antigens of tissue slides were repaired using EDTA antigen repair
   solution (Servicebio, G1203). After completion of repair, the tissues
   were blocked with 10% donkey serum for 30 min. After the addition of
   the primary antibody, the slices were incubated in a wet box overnight
   at 4° C. The next day, tissue slices were washed three times using
   1*PBS (PH7.4), then secondary antibodies were added and incubated for
   50 min at room temperature in the dark. Finally, tissue slices were
   sealed using an anti-fluorescence quench sealant (Servicebio, G1221).
   Images were acquired by an Upright fluorescence microscope (NIKON
   ECLIPSE C1, Nikon, Japan) and scanner (Pannoramic MIDI, 3DHISTECH).

Data Statistics

   All statistical analyses were performed with R version 4.2 software and
   its resource packages. The differences between the two groups were
   compared using t-test (normally distributed data) or Wilcoxon test
   (non-normally distributed data). p < 0.05 was considered statistically
   significant.

Results

Identification of Shared DEGs in UC and CRC

   Firstly, we identified 105 shared genes in UC and CRC through DEGs and
   WGCNA analyses. Next, we identified 35 SASPRGs from the intersection of
   105 genes with human putative SASP factor genes. Using machine learning
   algorithms and K-M survival analyses, we finally determined 2 core
   SASPRGs related to poor prognosis. The flowchart of identifying hub
   genes is shown in [57]Figure 1. In [58]GSE87466 dataset of UC patients,
   a total of 1642 DEGs were identified, including 1061 up-regulated genes
   and 581 down-regulated genes ([59]Figure 2A). In [60]GSE44076 dataset
   of CRC patients, 2950 DEGs were identified, including 1546 up-regulated
   genes and 1404 down-regulated genes ([61]Figure 2B). In order to
   identify shared DEGs, we intersected 1642 DEGs in [62]GSE87211 dataset
   with 2950 DEGs in [63]GSE44076 dataset. Finally, 632 genes were
   identified.

Figure 1.

   [64]Figure 1
   [65]Open in a new tab

   Flowchart of the research.

   Abbreviations: UC, ulcerative colitis; CRC, colorectal cancer; DEGs,
   differentially expressed genes; WGCNA, weighted gene co-expression
   network analysis; SASP, senescence-associated secretory phenotype;
   SASPRGs, SASP-related genes.

Figure 2.

   [66]Figure 2
   [67]Open in a new tab

   Identification of shared DEGs genes in UC and CRC. (A and B) Volcano
   plots of all DEGs in [68]GSE87466 and [69]GSE44076, green indicates
   downregulated DEGs, and red indicates upregulated DEGs. (C and D) The
   selection of a soft threshold in UC ([70]GSE87466) and CRC
   ([71]GSE44076). (E and F) Gene cluster dendrogram for UC and CRC by
   dynamic tree cut algorithm. (G and H) Heatmap of the association
   between modules and clinical traits in UC and CRC. (I) Venn diagram
   demonstrates the intersection of genes obtained by WGCNA and DEGs. (J)
   The enrichment analysis results of the GO and KEGG pathways.

   Abbreviations: BP, Biological Process; CC, Cell Component; MF,
   Molecular Function.

   To explore the biological functions and pathways of 632 shared DEGs. GO
   and KEGG analyses were performed based on these genes. As for GO, these
   genes were significantly enriched in the processes of leukocyte
   migration, cell chemotaxis, extracellular matrix organization,
   leukocyte chemotaxis, and humoral immune response. The KEGG pathway
   enrichment analysis indicated that these genes were significantly
   enriched in pathways, such as viral protein interaction with cytokine
   and cytokine receptor, protein digestion and absorption,
   cytokine-cytokine receptor interaction, and IL-17 signaling pathway
   ([72]Figure 2J).

Identification and Validation of Core SASPRGs Associated with UC and CRC

   WGCNA was conducted to identify biologically significant gene modules,
   enhancing our understanding of genes associated with the pathogenesis
   of UC and CRC in training datasets. We set the optimal soft threshold
   to 2 in [73]GSE87466 dataset and 9 in [74]GSE44076 dataset ([75]Figure
   2C and [76]D). Subsequently, 23 modules in the [77]GSE87466 dataset and
   13 modules in the [78]GSE44076 dataset were identified ([79]Figure 2E
   and [80]F). We selected the modules most relevant to clinical features
   as key modules. The result showed that the turquoise module was
   significantly associated with UC (correlation coefficient=0.72,
   p=3e-18) and the blue module was significantly associated with CRC
   (correlation coefficient=−0.93, p=1e-107) ([81]Figure 2G and [82]H). We
   selected the genes based on MM>0.6 and GS>0.4 as the candidate genes.
   Using Venn diagrams, we identified 105 co-pathogenic genes for UC and
   CRC in the training datasets ([83]Figure 2I).

   Next, we intersected 105 genes with human putative SASP factor genes,
   resulting in 35 SASPRGs ([84]Figure 3A). To further identify core
   SASPRGs, we used 113 combinations of 10 machine learning algorithms and
   calculated their AUC values to evaluate the diagnostic efficacy of each
   combination ([85]Figure 3B and [86]C). We found that the
   Enet[alpha=0.2] and Lasso+GBM models had the highest average AUC value
   across three UC datasets (0.991) and three CRC datasets (0.995).
   Therefore, we selected the genes identified by the Enet[alpha=0.2] and
   Lasso+GBM algorithms as candidate SASPRGs. We obtained 12 core SASPRGs
   through a Venn diagram (MMP7, S100A11, TGFBI, EGFL6, CXCL3, NPY,
   TNFSF11, CTSH, ANXA6, STC2, LRRN2, SRPX2) ([87]Figure 3D).

Figure 3.

   [88]Figure 3
   [89]Open in a new tab

   Identification of core SASPRGs through multiple machine learning
   algorithms. (A) Venn diagram demonstrates the intersection of genes
   obtained by shared genes and SASP. (B and C) A total of 113
   combinations of machine learning algorithms for the core SASPRGs. The
   AUC of each model was calculated across 3 datasets. (D) Venn diagram
   demonstrates the intersection of core SASPRGs obtained by
   Enet[alpha=0.2] and Lasso+GBM algorithms.

   To further evaluate the predictive power of core SASPRGs for CRC
   patient prognosis, Cox regression and K-M survival analyses were
   performed using TCGA-COAD data, which revealed that TGFBI, ANXA6, and
   NPY among the 12 core SASPRGs were potentially associated with poorer
   prognosis (All p values were less than 0.05) ([90]Figure 4A and [91]D).
   Interesting, the results of the Kaplan–Meier Plotter platform
   ([92]http://kmplot.com/analysis/) showed that there was no significant
   difference in recurrence-free survival probability of CRC patients in
   the high NPY group compared with the low NPY group (p=0.26) ([93]Figure
   4E). Consequently, we excluded NPY from further validation.
   Subsequently, we constructed a clinical characteristic map for CRC
   patients. Patients with varying expression levels of ANXA6 and TGFBI
   exhibited distinct patterns of clinical characteristics
   ([94]Supplementary Figure 1). Furthermore, the [95]GSE92415 dataset was
   used to detect the diagnostic efficacy of TGFBI and ANXA6 in assessing
   disease risk for UC patients. The results showed that both TGFBI and
   ANXA6 were significantly higher in the high Mayo score group (>5 score)
   compared to the low Mayo score group (<=5 score), with AUC values of
   0.772 (TGFBI) and 0.773 (ANXA6), respectively ([96]Figure 4B and
   [97]C). These results indicate that TGFBI and ANXA6 have clinical
   significance in stratifying disease risk for UC and CRC patients.

Figure 4.

   [98]Figure 4
   [99]Open in a new tab

   Clinical value of ANXA6 and TGFBI in UC and CRC patients. (A)
   Univariate Cox regression analysis for disease-free survival in
   TCGA-COAD cohorts. (B and C) Clinical value of ANXA6 and TGFBI in UC
   patients. (D) K-M survival analysis of high-expression versus
   low-expression groups of TGFBI, ANXA6, and NPY in TCGA-COAD dataset.
   (E) K-M survival analysis of high-expression versus low-expression
   groups of TGFBI, ANXA6, and NPY in Kaplan–Meier Plotter platform.

   Furthermore, we evaluated expression level and diagnostic efficacy of
   ANXA6 and TGFBI in validation datasets of UC and CRC. The results
   showed that the expression levels of TGFBI and ANXA6 were significantly
   higher in the disease group compared to the normal group. Additionally,
   TGFBI and ANXA6 effectively distinguished between UC patients, CRC
   patients, and healthy controls (AUC>0.75) ([100]Figure 5A and [101]B).
   These findings indicate that TGFBI and ANXA6 have excellent diagnostic
   performance for UC and CRC.

Figure 5.

   [102]Figure 5
   [103]Open in a new tab

   Validation of the core SASPRGs in UC and CRC datasets. (A) Expression
   levels and diagnostic potential of ANXA6 and TGFBI in the validation
   datasets for UC. (B) Expression levels and diagnostic potential of
   ANXA6 and TGFBI in the validation datasets for CRC. *p < 0.05; ** p <
   0.01; *** p < 0.001. **** p < 0.0001.

Biological Functions Associated with ANXA6 and TGFBI in UC and CRC

   To explore the biological functions related to ANXA6 and TGFBI, we
   conducted a GSEA analysis using the genes that are most closely
   associated with them (p<0.05, |r|>0.5). The results show that ANXA6 and
   TGFBI are associated with primary immunodeficiency, NF-kappa B
   signaling pathway, Th1 and Th2 cell differentiation, and B cell
   receptor signaling pathway in UC ([104]Figure 6A). In CRC, ANXA6 is
   associated with focal adhesions, protein digestion and absorption, ECM
   receptor interaction, cell adhesion molecules, and NF-kappa B signaling
   pathway, while TGFBI is related to the hedgehog signaling pathway,
   protein digestion and absorption, ECM receptor interaction, Wnt
   signaling pathway, and carbohydrate digestion and absorption
   ([105]Figure 6B).

Figure 6.

   [106]Figure 6
   [107]Open in a new tab

   Biological functions associated with ANXA6 and TGFBI in UC and CRC. (A)
   GSEA analysis shown the biological processes related to ANXA6 and TGFBI
   in UC. (B) GSEA analysis shown the biological processes related to
   ANXA6 and TGFBI in CRC.

Expression Patterns of ANXA6 and TGFBI in scRNA-Seq

   We analyzed single-cell transcriptome data from the colon tissues of UC
   and CRC patients to further validate the expression patterns of core
   SASPRGs. We identified 35241 and 45270 cells in UC and CRC colonic
   tissues after data quality control, respectively ([108]Figure 7A and
   [109]B). Based on the expression of marker genes, we annotated 6 cell
   types for UC and 10 cell types for CRC patients ([110]Figure 7C and
   [111]D). Subsequently, we evaluated the expression levels of core
   SASPRGs in different cell types, and the results showed that ANXA6 and
   TGFBI were mainly expressed in monocytes and fibroblasts ([112]Figure
   7E–H). Furthermore, we further validated the expression levels of ANXA6
   and TGFBI in both the disease group and the control group at the
   single-cell level. As shown in [113]Figure 7I–J, expression levels of
   ANXA6 and TGFBI were significantly higher in the colon tissues of UC
   and CRC patients (All p values were less than 0.05).

Figure 7.

   [114]Figure 7
   [115]Open in a new tab

   Single-cell transcriptomic characterization of colon tissues from UC
   and CRC patients. (A and B) The tSNE visualization of different cell
   groups in UC ([116]GSE162335) and CRC ([117]GSE231559) patients. (C and
   D) The tSNE visualization after annotation of cell groups in UC and CRC
   patients. (E and F) The tSNE visualization of expression levels of
   ANXA6 and TGFBI in UC and CRC patients. (G and H) The dot plot
   illustrates the expression levels of ANXA6 and TGFBI across different
   cell populations in UC and CRC patients. (I and J) Violin plot showing
   expression levels of ANXA6 and TGFBI in Normal vs UC and Normal vs CRC.

Immune Cell Infiltration and Prediction of the Immunotherapy Response

   Previous studies have shown that the occurrence and development of UC
   and CRC are closely related to immune response pathways. Therefore, we
   performed immune cell infiltration analysis to understand the immune
   infiltration status of UC and CRC. We found a significant difference in
   the level of immune cell infiltration between the disease group and the
   control group ([118]Figure 8A and [119]B and [120]Supplementary Figure
   2A–[121]D). Compared with the control group, the UC group showed higher
   levels of immune cells, including Activated B cell, Natural killer
   cell, Macrophage, Monocyte, T cell CD4 memory activated, Fibroblasts,
   etc. In CRC group, Activated CD8 T cell, Gamma delta T cell, Activated
   B cell, Immature B cell, CD56dim natural killer cell, Plasma cells,
   Macrophages M2, and T cells regulatory showed higher levels, while
   Activated CD4 T cell, Central memory CD4 T cell, Immature dendritic
   cell, Macrophages M1, and Fibroblasts showed lower levels. [122]Figure
   8C and [123]E show the correlation heatmaps among different immune
   cells in UC and CRC, further indicating the immune characteristics
   exhibited by the disease group, as well as highlighting the
   interactions among various types of immune cells. In addition, we
   investigated the relationship between core SASPRGs and immune cell
   components in UC and CRC. The results showed that TGFBI and ANXA6 were
   significantly positively correlated with Activated CD4 T cell, Central
   memory CD4 T cell, Type 2 T helper cell, Regulatory T cell, Natural
   killer cell, Natural killer T cell, and Immature dendritic cell, but
   negatively correlated with Activated CD8 T cell, Activated B cell, and
   CD56dim natural killer cell ([124]Figure 8D and [125]F). In order to
   further explore the relationship between ANXA6 and TGFBI and immune
   cells. We divided them into high-expression and low-expression groups
   based on the median. The results showed that compared with the
   low-expression group, the infiltration level of CD56 bright natural
   killer cells in the high-expression groups was significantly increased
   in UC ([126]Supplementary Figure 2E and [127]F), while in CRC, the
   infiltration levels of T follicular helper cells and Memory B cells
   were significantly increased ([128]Supplementary Figure 2G and [129]H).
   This suggests that TGFBI and ANXA6 may play a key role in the
   pathogenesis of UC and CRC by interfering with immune cell
   infiltration.

Figure 8.

   [130]Figure 8
   [131]Open in a new tab

   Analysis of immune cell infiltration and prediction of the
   immunotherapy response. (A and B) The ssGSEA analysis shown enrichment
   score of immune cells between the UC and CRC disease groups and the
   control group. (C and E) Heatmap of correlations between different
   immune cells in colon tissues of patients with UC and CRC. (D and F)
   Correlation analysis between infiltrating immune cells and core SASPRGs
   in colon tissues of patients with UC and CRC. (G and I) Expression of
   the immune checkpoint genes PDCD1, CD274, and CTLA4 in different
   SASPRGs subgroups. (H and J) TIDE and T-cell exclusion scores in
   different SASPRGs subgroups. *p < 0.05; ** p < 0.01; *** p < 0.001.
   **** p < 0.0001.

   To evaluate the sensitivity of the core SASPRGs to predict
   immunotherapy response in CRC patients, we first evaluated the
   expression levels of PDCD1, CD274, and CTLA4 in different SASPRGs
   subgroups ([132]Figure 8G and [133]I). High expression of immune
   checkpoints often suggests that these patients may be suitable for
   immunotherapy. We then used TIDE to predict immunotherapy efficacy in
   different SASPRGs subgroups. As expected, the high TGFBI group had
   lower expression levels of immune checkpoint and higher TIDE scores and
   T-cell exclusion scores, suggesting that patients did not benefit from
   immunotherapy (All p values were less than 0.05). Notably, the high
   ANXA6 group had higher expression levels of immune checkpoint, while
   the TIDE scores and T-cell exclusion scores were also high ([134]Figure
   8H and [135]J).

Screening of Candidate Drugs for High ANXA6 and TGFBI Patients

   To screen candidate drugs for CRC patients with high expression of
   ANXA6 and TGFBI, we evaluated drug sensitivity in different subgroups
   of patients using drug response data provided by GDSC and CTRP
   databases. Subsequently, we selected the top three drugs or compounds
   with the highest sensitivity as candidates ([136]Figure 9A and [137]B).
   To further analyze the interactions between candidate drugs and
   SASP-associated proteins, we performed molecular docking. The results
   demonstrated that ANXA6 (binding energy: Staurosporine: −7.11 kcal/mol,
   AT7867: −6.31 kcal/mol, BRD_K44224150: −5.93 kcal/mol) and TGFBI
   (binding energy: SCH772984_1564: −7.47 kcal/mol, VX-11e_2096: −6.21
   kcal/mol, ERK_6604_1714: −6.47 kcal/mol) exhibited a certain degree of
   binding affinity for the three drugs ([138]Figure 9C and [139]D). These
   candidate drugs have the potential to become therapeutic drugs for CRC.

Figure 9.

   [140]Figure 9
   [141]Open in a new tab

   Drug sensitivity analysis and Construction of a molecular regulator
   network. (A and B) Outcome of drug sensitivity analysis of three
   potential therapeutic drugs in different SASPRGs subgroups. (C and D)
   3D docking model of three candidate drugs and SASP-associated proteins.
   (E) Molecular regulatory network based on ANXA6 and TGFBI. Green
   hexagons represent miRNAs, and blue arrows signify TFs. *p < 0.05; ** p
   < 0.01; *** p < 0.001. **** p < 0.0001.

Construction of a Molecular Regulator Network

   A thorough understanding of the regulatory mechanisms of core SASPRGs
   is conducive to identifying key biomarkers and enhancing the
   understanding of diseases. Therefore, we constructed a molecular
   regulatory network based on ANXA6 and TGFBI to identify their
   transcriptional and post-transcriptional regulatory molecules.
   Ultimately, we identified 25 miRNAs targeting TGFBI, 23 miRNAs
   targeting ANXA6, 13 transcription factors regulating TGFBI, and 9
   transcription factors regulating ANXA6. Notably, four transcription
   factors, USF2, E2F1, TFAP2A, and FOXC1, jointly regulate both ANXA6 and
   TGFBI ([142]Figure 9E).

Validation of Core SASPRGs Associated with UC and CRC in vivo

   To confirm the previous findings, UC model was constructed using DSS.
   CRC model previously established by our team was used in this study
   ([143]Figure 10A and [144]B). Colon tissues of mice were collected for
   histological staining. Compared with the control group, the UC group
   had disordered colonic mucosa structure, inflammatory cell
   infiltration, and reduced goblet cells, while the CRC group had
   disordered colonic gland arrangement and increased cell atypia
   ([145]Figure 10C). Further evaluation showed that the colon length of
   the UC group mice was significantly shortened, and the intestinal
   mucosa showed erosion, bleeding and edema ([146]Figure 10D). Compared
   with the control group, the body weight of UC group mice decreased
   significantly with the increase of treatment time, and the DAI score
   continued to increase (p<0.05) ([147]Figure 10E and [148]F). In
   addition, we also observed that UC group mice gradually developed
   symptoms such as depression, hair loss, and diarrhea during the
   treatment period. At the same time, we determined the mRNA expression
   levels of ANXA6 and TGFBI in the colon tissues of each group of mice by
   qRT-PCR. The results showed that the mRNA levels of ANXA6 and TGFBI
   were significantly increased in the UC and CRC groups compared with the
   control group (All p values were less than 0.05) ([149]Figure 10H).
   Immunofluorescence staining further revealed the expression of ANXA6
   and TGFBI in the colon tissues of each group of mice ([150]Figure 10G
   and [151]I).

Figure 10.

   [152]Figure 10
   [153]Open in a new tab

   Expression of core SASPRGs in mouse models of UC and CRC. (A and B)
   Flow chart of UC mice and CRC mice model establishment. (C) HE staining
   of colonic tissues from normal, UC, and CRC mice. (D) Colon length of
   normal and UC mice. (E and F) Body weight and DAI score of normal and
   UC mice. (G and I) Immunofluorescence staining and relative
   fluorescence intensity quantification of colonic tissues from normal,
   UC, and CRC mice. (H) The mRNA expression of ANXA6 and TGFBI in colon
   tissue detected by qRT-PCR. *p < 0.05; ** p < 0.01; *** p < 0.001. ****
   p < 0.0001.

   Abbreviations: DSS, dextran sulfate sodium; AOM, azomethane; DAI,
   disease activity index.

Discussion

   UC is a common chronic inflammatory bowel disease.[154]^26 There is
   evidence that active UC is prone to progressing to UCRCC, and the risk
   of UCRCC in UC patients is 2.4 times higher than that in the general
   population.[155]^27 SASP, which has been widely studied in the fields
   of inflammatory diseases and cancer in recent years, plays a key role
   in the development of a variety of diseases.[156]^28–31 According to
   the report, GDF15 promotes the malignant progression of CRC by
   activating the MAPK and PI3K signaling pathways as an essential SASP
   factor.[157]^23 Furthermore, research indicates that SPP1 may be
   associated with the inhibitory immune microenvironment mediated by
   senescent cells in CRC.[158]^24 However, there is still a lack of
   SASP-based studies on the molecular mechanisms of UC and CRC.
   Therefore, exploring the pathogenic genes based on SASP of UC and CRC
   patients not only enhances our understanding of the pathogenesis of
   these diseases, but also helps clinicians to make better decisions and
   improve patient prognosis.

   In this study, we identified 35 overlapping SASPRGs in UC and CRC
   patients by WGCNA as well as differential expression analysis. To
   comprehensively evaluate the processing performance of different
   machine learning algorithms and address the specific characteristics of
   our datasets, we strategically selected 10 representative algorithms.
   These diverse sets of algorithms allow for a robust comparison across
   different models, including Lasso, support vector machine (SVM),
   elastic net (Enet), Ridge, generalized linear model (GLM), random
   forest (RF), generalized boosted regression modeling (GBM), linear
   discriminant analysis (LDA), extreme gradient boosting (XGBoost), Naïve
   Bayes. Subsequently, we identified 12 core SASPRGs using 113
   combinations of 10 machine learning algorithms, among which ANXA6 and
   TGFBI expression were significantly increased in disease groups in
   multiple datasets and were significantly associated with poor prognosis
   in CRC patients. Similar results were obtained in mouse models of UC
   and CRC. Therefore, these two genes were selected for further analysis.
   By analyzing the scRNA-seq data, ANXA6 and TGFBI were found to be
   mainly expressed in fibroblasts and monocytes. In addition, we also
   found that CRC patients with high ANXA6 and TGFBI expression had higher
   TIDE scores and T cell exclusion scores, indicating that these patients
   seem to have a poor response to immunotherapy. Further investigation
   showed that candidate drugs from both pharmacogenomic databases
   exhibited better therapeutic sensitivity in the high ANXA6 group
   (Staurosporine, AT7867, BRD_K44224150) and the high TGFBI group
   (SCH772984_1564, VX-11e_2096, ERK_6604_1714).

   Annexin-a6 (ANXA6), a highly conserved Ca2+ -dependent membrane binding
   protein, is the largest of all annexin families.[159]^32 It is located
   in the cytoplasm and attached to the phospholipid membrane of cells,
   and is involved in a variety of cellular functions, mainly vesicle
   trafficking and membrane repair.[160]^33 Studies have shown that ANXA6
   is differentially expressed in gastric cancer, colorectal cancer,
   melanoma, cervical cancer, and breast cancer, and plays an important
   role in the formation, development, and drug resistance of
   tumors.[161]^34–37 ANXA6 shows dual effects in different tumors.
   Generally, ANXA6 is considered as a tumor suppressor gene mainly
   because ANXA6 negatively regulates the phosphorylation of EGFR and the
   downstream MAPK and PI3K-AKT signaling pathways.[162]^38–40 Some
   studies have found that ANXA6 is detected to be elevated in the feces
   of patients with colorectal cancer, suggesting that ANXA6 is a
   biomarker for early detection of colorectal cancer.[163]^37 This is
   similar to our conclusion. In addition, ANXA6 can promotes the
   secretion of pro-inflammatory cytokines that mediate the malignant
   progression of triple-negative breast cancer.[164]^41 Increasing
   evidence suggests that ANXA6 mediates tumor drug resistance and
   malignant progression by participating in cholesterol metabolism in
   LEs, autophagy regulation, interacting with TGFR, and
   immune-inflammatory responses.[165]^36^,[166]^42–45

   TGFBI (transforming growth factor beta induced) is a secreted protein
   stimulated by TGF-β (transforming growth factor β).[167]^46 As a
   component of extracellular matrix, TGFBI is involved in mediating tumor
   cell adhesion and motility, metabolic reprogramming, angiogenesis, and
   immune microenvironment homeostasis.[168]^47–49 In recent years, many
   studies have found that TGFBI may promote tumor cell proliferation by
   enhancing intercellular signaling, activating the PPARγ or
   PI3K/AKT/mTOR signaling pathway, and is significantly associated with
   poor prognosis of patients, including cholangiocarcinoma, urothelial
   carcinoma, Wilms tumor, and glioma.[169]^50–53 A study on metastatic
   colon cancer found that TGFBI promoted angiogenesis by activating
   spliceosome and lysosome-related pathways.[170]^54 In addition, there
   was a significant correlation between the expression level of TGFBI and
   the degree of immune cell infiltration.[171]^55 TGFBI secreted by
   immune cells promotes tumor progression by creating an inhibitory
   immune microenvironment.[172]^56 Interestingly, high TGFBI expression
   was significantly associated with natalizumab resistance in lung cancer
   patients.[173]^57

   Currently, the most commonly used chemotherapeutic agents for CRC are
   fluorouracil, oxaliplatin, and irinotecan.[174]^58 Our results showed
   that both ANXA6 and TGFBI subgroups showed differences in drug
   sensitivity to fluorouracil, indicating heterogeneity among the
   subgroups. Staurosporine, AT7867 and BRD_K44224150 showed good
   treatment sensitivity in the high ANXA6 group. Staurosporine is an
   ATP-competitive, nonselective protein kinase inhibitor. Staurosporine
   has been shown to inhibit the progression of gastric cancer by
   inhibiting CTNNB1 transcription and Wnt/β-catenin signaling
   activation.[175]^59 AT7867 is a potent AKT inhibitor with demonstrated
   ability to inhibit tumor cell proliferation and induce tumor cell
   apoptosis, especially in colorectal cancer.[176]^60 Interestingly, all
   three therapeutically sensitive drugs screened in the high TGFBI group
   were the ERK inhibitors. Studies have shown that ERK signaling pathway
   can promote tumor cell migration and invasion by regulating the
   expression of TGFBI.[177]^61^,[178]^62 These results provide a solid
   rationale for TGFBI as a novel biomarker for CRC.

   In this study, we utilized multiple UC and CRC datasets from the GEO
   and the TCGA database for a comprehensive and integrated analysis of
   SASPRGs associated with UC and CRC. However, it is important to
   acknowledge that this study contains certain limitations. Firstly, our
   research indicates that ANXA6 and TGFBI, as the functional executors of
   senescent cells, have significant clinical value in patients with UC
   and CRC. However, further in vitro and in vivo experiments are needed
   for verification. Second, the limited sample size may affect the
   generalization ability of the conclusions. Thirdly, the effects of the
   selected drugs have not been verified in the CRC model.

Conclusion

   In conclusion, NXA6 and TGFBI, as core SASPRGs associated with UC and
   CRC, showed good diagnostic efficacy. Notably, ANXA6 and TGFBI appear
   to be effective therapeutic targets in CRC patients. In the future,
   multi-center and large-sample data sets are needed to verify our
   conclusions, and subsequent in vivo and in vitro experiments are needed
   to further explore the possible molecular mechanisms.

Funding Statement

   This study was supported by Scientific Research Fund of Henan
   Provincial Science and Technology Department (No. 182102310302).

Data Sharing Statement

   The datasets supporting the conclusions of this article are available
   in the TCGA ([179]https://portal.gdc.cancer.gov/) and GEO
   ([180]https://www.ncbi.nlm.nih.gov/geo/) databases.
   [[181]https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE87466/GSE10
   7499/GSE59071/GSE92415/GSE44076/GSE17536/GSE87211]. The analysis
   methods and used packages are illustrated in the “Methods” section. All
   other R code and analyses are available from the author upon request.

Ethics Approval

   In accordance with national legislative guidelines on human data, our
   research does not require approval, is derived from Article 32, items
   (1) and (2) of the “Measures for the Ethical Review of Life Science and
   Medical Research Involving Human Subjects” issued by China on February
   18, 2023. The legislation in these two items is that research using
   legally obtained public data, data generated through observation
   without interfering with public behavior, or research using anonymized
   information data can be exempted from ethical approval. This study
   complied with the Declaration of Helsinki and was approved by the
   Ethical Review Committee for Animal Experiments of the Nanjing Medical
   University (2011031).

Author Contributions

   All authors made a significant contribution to the work reported,
   whether that is in the conception, study design, execution, acquisition
   of data, analysis and interpretation, or in all these areas; took part
   in drafting, revising or critically reviewing the article; gave final
   approval of the version to be published; have agreed on the journal to
   which the article has been submitted; and agree to be accountable for
   all aspects of the work.

Disclosure

   The authors report no conflicts of interest in this work.

References