Multi-Omics analysis elucidates tumor microenvironment and intratumor microbes of angiogenesis subtypes in colon cancer

Abstract

BACKGROUND

Angiogenesis plays an important role in colon cancer (CC) progression.

AIM

To investigate the tumor microenvironment (TME) and intratumor microbes of angiogenesis subtypes (AGSs) and explore potential targets for antiangiogenic therapy in CC.

METHODS

The data were obtained from The Cancer Genome Atlas database and Gene Expression Omnibus database. K-means clustering was used to construct the AGSs. The prognostic model was constructed based on the differential genes between two subtypes. Single-cell analysis was used to analyze the expression level of SLC2A3 on different cells in CC, which was validated by immunofluorescence. Its biological functions were further explored in HUVECs.

RESULTS

CC samples were grouped into two AGSs (AGS-A and AGS-B) groups and patients in the AGS-B group had poor prognosis. Further analysis revealed that the AGS-B group had high infiltration of TME immune cells, but also exhibited high immune escape. The intratumor microbes were also different between the two subtypes. A convenient 6-gene angiogenesis-related signature (ARS), was established to identify AGSs and predict the prognosis in CC patients. SLC2A3 was selected as the representative gene of ARS, which was higher expressed in endothelial cells and promoted the migration of HUVECs.

CONCLUSION

Our study identified two AGSs with distinct prognoses, TME, and intratumor microbial compositions, which could provide potential explanations for the impact on the prognosis of CC. The reliable ARS model was further constructed, which could guide the personalized treatment. The SLC2A3 might be a potential target for antiangiogenic therapy.

Key Words: Angiogenesis, Colon cancer, Immunotherapy, K-means, Single-cell analysis

Core Tip: Angiogenesis plays an important role in colon cancer (CC) progression. This study identified two angiogenesis subtypes (AGSs) with significantly different prognoses, tumor microenvironment, intratumor microbiota, drug sensitivity, and cancer-related pathways single-sample gene set enrichment analysis scores in patients with CC. Based on the two AGSs, a convenient 6-gene angiogenesis-related signature (ARS), was established to predict the prognosis in CC patients. SLC2A3 was selected as the representative gene of ARS, which was higher expressed in endothelial cells and promoted the migration of HUVECs.

INTRODUCTION

Colon cancer (CC) is the third most malignancy worldwide and is the second leading cause of cancer-related deaths[1]. Angiogenesis plays a critical role in the progression and metastasis of CC[2,3]. Antiangiogenic therapy has emerged as an effective means of controlling cancer, and several angiogenesis inhibitors, such as aflibercept and bevacizumab, have demonstrated efficacy in treating CC[4,5].

Angiogenesis and optimizing the efficacy of antiangiogenic therapies have become intense arears of research[6-8]. Systematic studies on angiogenesis-related genes (ARGs) are helpful in understanding the mechanism of cancer, discovering new therapeutic targets for anti-angiogenesis, and uncovering novel combination therapy strategies. Some systematic studies of angiogenesis have been conducted in various cancers, such as gastric cancer and breast cancer, which have constructed the angiogenesis-related prognosis models and revealed the correlation between the tumor microenvironment (TME) and angiogenesis[9-12]. Additionally, the intratumor microbiome is closely associated with the occurrence and development of CC[13,14], and this relationship involves the association with angiogenesis in CC. However, few studies have focused on the systematic study of angiogenesis in CC and the correlation between angiogenesis and TME and intratumor microbial compositions in CC requires further investigation.

In this study, we analyzed the expression and mutational signature of ARGs and identified two angiogenesis subtypes (AGSs). We also differentiated between the two AGSs based on tumor metabolism, TME, and intratumor microbial compositions in CC. Based on the two AGSs, we further developed an angiogenesis-related signature (ARS) model to determine the AGSs, predict the prognosis and assess therapeutic effects, which could guide individualized treatment for CC. Moreover, we revealed the potential mechanisms of SLC2A3 to angiogenesis, which could serve as a potential therapeutic target for anti-angiogenesis in the future.

MATERIALS AND METHODS

Dataset collection

The transcriptome data, somatic mutation data, copy number variation (CNV) files, and corresponding clinicopathological information of CC were retrieved from The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) program, and GSE17536 from the Gene Expression Omnibus (GEO) database was utilized to acquire clinical parameters and normalized gene expression data. The intratumor microbial compositions were retrieved from The Cancer Microbiome Atlas database[15]. single-cell RNA sequencing (scRNA-seq data) (GSE200997) was also retrieved from the GEO database. Samples that lacked significant clinicopathological or survival information were eliminated from further analysis. 36 ARGs in the HALLMARK_ANGIOGENESIS gene set were obtained from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) (Supplementary Table 1)[16]. The overall study flow chart is summarized in Figure 1. Samples without complete survival data were removed from the study.

Figure 1 Flowchart for the analytical process in the study. IHC: Immunohistochemical; IF: Immunofluorescence; TCGA-COAD: The Cancer Genome Atlas Colon Adenocarcinoma. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

K-means consensus cluster

Univariate Cox regression analysis was used to extract ARGs via the R package "survival" (P < 0.05). The K-means consensus cluster was performed to explore the potential molecular subtypes between the CC patients using the R package "ConsensuClusterPlus". The cluster count (k) was set from 1 to 9, and the best optimal k value was selected for further investigation.

TME

The "ESTIMATE" package in R was used to evaluate the immune score, the stromal score, and tumor purity. The infiltrating fractions of immune cells were also identified with a single-sample gene set enrichment analysis (ssGSEA) algorithm. The gene set for marking each TME infiltrating immune cell type was obtained from the study, which stored various human immune cell subtypes, including activated CD8 T cell, activated dendritic cell, macrophage, natural killer T cell, regulatory T cell[17,18].

Functional enrichment analysis

Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway analysis were performed using the "clusterProfiler" package in R. Additionally, Gene Set Enrichment Analysis (GSEA) was performed to identify the differences in the set of genes expressed between the two subtypes in the enrichment of the KEGG pathway and HALLMARK pathway from the MSigDB database (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp)[19,20]. The number of permutations were performed 1000 times for each analysis, and the pathways with adjusted P-value < 0.05 and q-value < 0.05 were considered statistically significant.

Tumor mutation burden and mutation

Comprehensive mutation analysis and tumor mutation burden (TMB) scores for each patient were calculated using the "maftool" R package[21].

Microbiota abundance analysis and microbial diversity

Partial least squares discrimination analysis (PLS-DA) was conducted to evaluate the overall differences in microbiota profiles between the two subtypes. The abundance of microbes at the taxonomic levels of phylum, class, order, family, and genus was calculated, and the microbiota profiles of the top 10 most abundant microbes were summarized.

Prediction of immunotherapy and drug response

For external validation of the ARS score, we downloaded an independent dataset, IMvigor210, from a freely available data package. To externally validate the predicted response to immunotherapy, we calculated ARS scores for each patient in IMvigor-210. (http://research-pub.gene.com/IMvigor210CoreBiologies). In the IMvigor-210 cohort, anti-PD-L1 immunotherapy responses were also compared between groups with high or low ARS scores.

The immunopheno score (IPS) of CC patients was obtained from The Cancer Immunome Database (https://www.tcia.at/home). Moreover, the immune-checkpoint inhibitors response was assessed through the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/Login/)[22,23].

In order to estimate the chemotherapeutic response based on the half-maximal inhibitory concentration (IC₅₀), common clinical drugs were selected. The R package "pRRophetic" was used to estimate the chemotherapeutic response.

Differential gene identification and Gene Set Variation Analysis

Differential genes (DEGs) between AGSs were identified by the "limma" package. Furthermore, genes with |log2 value (FC)| > 1 and P < 0.05 after adjusting for false discovery rate (FDR) were considered angiogenesis-related DEGs. GO and KEGG enrichment analyses were further performed using these DEGs.

Gene Set Variation Analysis (GSVA) was conducted on the gene profile through the "GSVA" R package to compare the distinctions of the biological processes between low and high-score groups. GSVA, a non-parametric and unsupervised method, can evaluate the pathway variations or biological processes through an expression matrix sample. The "c5.go.v7.4.symbols.gmt" gene set from the molecular signatures database (https://www.gsea-msigdb.org/gsea/msigdb) was used as the reference gene set. FDR < 0.05 indicated a statistically significant enrichment pathway.

Development and validation of the angiogenesis-related risk and nomogram

We analyzed the prognostic significance of these DEGs between two subtypes by univariate Cox regression using the R package 'survival'. Subsequently, we selected the most robust prognostic gene signatures in the Least absolute shrinkage and selection operator (LASSO) regression model using the R package "glmet"[23]. Then, an angiogenesis-related risk (ARR) score was calculated by the corresponding coefficients of selected signatures. The ARR formula was established as follows: Score = Σi Coefficient (mRNA) × Expression (mRNA).

According to the median value of ARR, samples were divided into low and high groups. We then evaluated the association between the two subsets on the risk score.

Nomograms were developed to predict 1-, 3-, and 5-year overall survival (OS) for CC patients with their risk scores and other clinicopathological characteristics. Next, we performed receiver operator characteristic (ROC) and calibration curve analysis to verify the clinical reliability of the established nomogram in TCGA-COAD patients and external validation cohort (GSE17536).

Single-cell RNA-seq analysis

The “Seurat” package explored the transcriptional heterogeneity of cells within the normal and CC tissue. To cluster cells, we used the FindCluster function. To find markers, we used the Findmarkers function. According to the marker genes, we annotated different cell types.

Cell culture and transient transfection

HUVECs were gifted from Beijing Anzhen Hospital, which were then incubated in ECM (ScienCell, United States, 1001) supplemented with 5% FBS and 1% penicillin-streptomycin in a humidified incubator and maintained at 37°C with 5% carbon dioxide. Plasmids encoding SLC2A3 were constructed by Changsha Youbio Co, Ltd (NM_006931). Transfection was performed with Lipofectamine 2000 (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's guidelines. Western Blot assessed the transfection efficiency of SLC2A3 overexpression after 48 h post-transfection.

Immunohistochemical and immunofluorescence staining

Immunohistochemical (IHC) and immunofluorescence (IF) staining were performed on normal and paired tumor tissue. Immunohistochemistry staining SLC2A3 protein was detected by immunohistochemistry. Colocalization of SLC2A3 and CD31 was observed using IF. The tissues, which were incubated with rabbit polyclonalanti-GUT3 (proteintech 20403-1-AP, IHC: 1:200, IF: 1:100) and rabbit polyclonalanti-CD31 (proteintech, IF: 1:100). SLC2A3 expression was evaluated by using a system considering the staining intensity (0 means negative 1 means weak; 2 means moderate; and 3 means strong) and the percentage of positively stained cells (< 5% = 0, 5% to < 25% = 1, 25% to 50% = 2, > 50% to < 75% = 3, > 75% = 4). A staining intensity scale (0 means negative, 1 means weak, 2 means moderate, and 3 means strong) and percentage of positively stained cells (< 5% = 0, 5% to < 25% = 1, 25% to 50% = 2, > 50% to < 75% = 3, > 75% = 4) were used to evaluate SLC2A3 expression. The final scores were calculated by multiplying the extent scores by the intensity scores. The study was approved by the Medical Ethics Committee of Beijing Shijitan Hospital, Capital Medical University (Approval number: 2020-11).

Cell proliferation, migration, and tubule formation

Cell proliferation was determined using the EdU assays. Each group of HUVECs was seeded into 24-well plates at 2 × 10⁴ cells per well. Then, the cells were incubated with fresh media containing 50 nM EdU reagent (RiboBioInc, cat. C00003, China) for 2 h at 37 °C, then fixed with 4% paraformaldehyde solution, followed by Apollo and Hoechst staining. Cell proliferation was assessed by the relative EdU-positive cell index. Relative EdU-positive cells index = the number of proliferating cells stained in Apollo/ total number of cells stained in Hoechst.

The migration abilities of HUVECs were assessed using the transwell assay. HUVECs (2 × 10⁴ cells per well) were seeded in fresh medium (without serum) in the upper chamber. The lower chamber was filled with a complete medium. After 5 h, the upper chamber was washed, fixed, dyed, and photographed.

To assess the ability of HUVECs to form tube-like structures, HUVECs (1.5 × 10⁴ cells per well) were resuspended in media and seeded on growth factor-reduced Matrigel (Corning 356231) for 4 h. The number of branches and total segment lengths were quantified using Image J software (NIH, United States).

Statistical analysis

All statistical analyses were performed with R software (version 4.1.2). Statistical significance was set at a probability value of a two-sided P < 0.05. Wilcoxon test, One-way ANOVA, Student's t-test, and chi-square test were used to compare different groups. The survival differences between the groups were estimated with Kaplan-Meier, which was compared by performing a log-rank test. Univariate, multivariate Cox regression analyses and LASSO analyses were performed to identify independent prognostic factors and establish the risk-score model. The ROC curves and calibration curves were compared to determine the predictive accuracy of the prognostic models.

RESULTS

General information on ARGs in CC

We obtained the gene set of HALLMARK ANGIOGENESIS with the 36 ARGs from the MSigDB database[16]. The expression levels of those 36 ARGs in The Cancer Genome Atlas Colon Adenocarcinoma tumor samples (n = 399) and normal samples (n = 39) were shown in a heatmap (Figure 2A). Generally, the GSVA angiogenesis score was significantly higher in CC samples as compared to the normal samples (Figure 2B). A total of 8 DEGs were identified between those two groups (Figure 2C), among which 6 were up-regulated (TIMP1, STC1, MSX1, VEGFA, SPP1), and the other 2 were down-regulated (SLCO2A1 and FGFR1) (Figure 2D). Univariate COX regression analysis determined 7 ARGs (JAG2, NRP1, PGLYRP1, POSTN, SERPINA5, STC1, and TIMP1) associated with prognosis (Figure 2E). For mutation analysis, VCAN was identified as the most commonly mutated ARGs (9%), followed by COL5A2 (6%), POSTN (6%), COL3A1(5%), and JAG (5%) (Figure 2F). CNV alterations were identified in all 36 ARGs (Figure 2G), and their site on chromosomes was displayed in a circle plot (Figure 2H).

Figure 2 Expression landscape of angiogenesis-related genes in colon cancer. A: Heatmap (green: low expression level; brown: high expression level) of the angiogenesis-related genes (ARGs) between the nontumor (blue) and tumor samples (red); B: Differences in angiogenesis pathway single-sample gene set enrichment analysis score between normal and colon cancer (CC) tissues; C: Volcano plot of the differential ARGs; D: The expression of 8 differentially expressed angiogenesis-related genes between normal and CC tissues; E: Univariate Cox regression analysis of 36 angiogenesis-related genes and seven genes with P < 0.05; F: Mutation frequencies and types of ARGs; G: The location of copy number variation (CNV) alteration of 36 ARGs on 23 chromosomes; H: Chromosomal localization of ARGs with CNV. ARGs: Angiogenesis-related genes. ^cP < 0.001.

Identification of AGSs

We clustered all CC samples into two subtypes via the K-means consensus cluster according to the expression levels of the ARGs (Figure 3A). The t-SNE analysis was used to confirm the intergroup difference between the subtypes (Figure 3B). A significant difference in OS was also identified between the subtypes (Figure 3C). Generally, AGS A (AGS-A, n = 249) showed a higher expression of ARGs compared to AGS B (AGS-B, n = 115) (Figure 3D). For GSEA analysis, AGS-A and AGS-B exhibited differences in many go terms, such as actin filament organization, activation of protein kinase activity, and actomyosin structure organization (Figure 3E and F).

Figure 3 Tumor consensus clustering based on angiogenesis-related genes. A: Two angiogenesis clusters were identified according to the K-means consensus clustering matrix (k = 2); B. The t-SNE analysis for the transcriptome profiles of two angiogenesis subtypes; C: Kaplan-Meier curves between two angiogenesis subtypes; D: The heatmap showed the differences in expression levels of angiogenesis-related genes between the two angiogenesis subgroups; E and F: Gene Set Enrichment Analysis for Go terms in two angiogenesis subtypes.

Tumor characteristics between two AGSs

The ssGSEA scores of many KEGG and hallmark pathways, including cytokine receptor interaction, allograft rejection, and cell adhesion molecules cams, were different between the AGS-A and AGS-B groups (Figure 4A and B). The estimated, immune, and stromal scores significantly differed between these two groups (Figure 4C). The percentages of 23 immune cells were also showed significant differences between the two groups (Figure 4D). The expression of PD-1-related checkpoint genes (PD1, PDL1, PDL2), CTLA4-related checkpoint genes (CTLA4, CD80, CD86), and agonists of T cell-activation-related checkpoint genes (OX40, GITR, TNFRSF9, ICOS, CD40, CD27, CD70) in the AGS-B group were higher than in the AGS-A group (Figure 4E-G). The IPS scores of the samples with negative expression of PD1 in the AGS-B group were lower than in the AGS-A group (Figure 4H-K). The TIDE scores in the AGS-B group were higher than in the AGS-A group (Figure 4L). Compared with the AGS-A group, the AGS-B group had a higher percentage of high degree microsatellite instability (MSI), a lower percentage of microsatellite stability (MSS), and a lower degree MSI (Figure 4M). The tumor purity in the AGS-A group was higher than in the AGS-B group (Figure 4N). There were no significant differences in TMB between the AGS-A and AGS-B groups (Figure 4O).

Figure 4 Identification of characterization in angiogenesis subtypes. A and B: The differences in Kyoto Encyclopedia of Genes and Genomes and HALLMARKER pathway single-sample gene set enrichment analysis scores between the two angiogenesis subgroups; C: Correlations between the two angiogenesis subgroups and tumor microenvironment score; D: The abundance of immune cells in two angiogenesis subtypes; E: Differential expression of PD-1, PD-L1, and PD-L2 between two angiogenesis subtypes (P < 0.001); F: Differential CTLA4, CD80, and CD86 expression between two angiogenesis subtypes (P < 0.001); G: Differential expression of OX40, GITR, TNFRSF9, ICOS, CD40, CD27 and CD70 between two angiogenesis subtypes (P < 0.001); H-K: Immunopheno score in two angiogenesis subtypes; L: Tumor Immune Dysfunction and Exclusion scores in two angiogenesis subtypes; M: The proportion of angiogenesis subtypes in the three modification patterns [microsatellite stability (MSI), blue; MSI-low red; MSI-high orange]; N: Differences in tumor purity between two angiogenesis subtypes in colon cancer (CC) patients (P < 0.001); O: Differences in tumor mutational burden between two angiogenesis subtypes in CC patients (P < 0.001). ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

There were 86 samples of CC in the TCMA database. A total of 11, 22, 38, 67, and 221 taxa were obtained for each sample at the phylum, class, order, family, and genus levels, respectively. There was no difference in alpha diversity (Shannon index, Simpson index, richness index, and Chao index) among different subtypes at the phylum, class, order, family, and genus levels, respectively (Figure 5A-E). We used PLS-DA maps to compare microbiome profile landscapes in the two subtypes. From the phylum to genus level, microbial composition features are increasingly differentiated (Figure 5F-J). Figure 5K-O showed the microbiota profiles of the top 10 most abundant microbes in the two subtypes at the phylum, class, order, family, and genus levels, respectively.

Figure 5 Identification of intratumoral microbial compositions in angiogenesis subtypes. A-E: The alpha diversity (Shannon index, Simpson index, richness index, and Chao index) among different subtypes at the phylum, class, order, family, and genus levels; F-J: The Partial least squares discrimination analysis maps among different subtypes at the phylum, class, order, family, and genus levels; K-O: The microbiota profiles of the top 10 most abundant microbe in two subtypes at the phylum, class, order, family, and genus levels.

Construction of an ARS model

A total of 279 DEGs were identified between AGS-A and AGS-B. Among these, 276 were higher in AGS-A, while three were higher in AGS-B (Figure 6A). These DEGs were enriched in some KEGG pathways and Go terms, including phagosome and staphylococcus aureus infection (Figure 6B and C). Among these DEGs, TIMP1, CERCAM, SERPINE1, FSTL3, ZNF385A, SCG2, CALB2, TPM2, SLC2A3, UCHL1, GPX3, FABP4, CRABP2, BST2 and TNNT1 were confirmed to be prognostic (Figure 6D). Then LASSO regression and multivariate Cox analysis were used to establish a predictive risk model based on these DEGs (Figure 6E). The risk model was accessed as follows: ARS score = (0.0532 × expression of TIMP1) + (0.0232 × expression of FSTL3) + (0.0260 × expression of CALB2) + (0.0715 × expression of SLC2A3) + (0.0479 × expression of FABP4) + (0.1456 × expression of TNNT1).

Figure 6 Identification of differentially expressed genes between angiogenesis subtypes and construction of the prognostic model. A: Differentially expressed genes (DEGs) between two angiogenesis subtypes; B and C: Kyoto Encyclopedia of Genes and Genomes and Gene Ontology enrichment analyses of DEGs between two angiogenesis subtypes; D: Univariate Cox regression analysis of each DEG between two angiogenesis subtypes and 13 genes with P < 0.05; E: Cross-validation for tuning the parameter selection in the Least absolute shrinkage and selection operator regression; F: Sankey diagram of subtype distributions in angiogenesis subtypes and angiogenesis-related score (ARS) groups; G: Differences in ARS among two angiogenesis subtypes in colon cancer patients (P < 0.001); H: Heatmap (green: low expression level; brown: high expression level) of the angiogenesis-related genes between the low and high ARS groups; I: Correlations between the expression of 36 Angiogenesis related-genes and selected 6 genes in the ARS model. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

All CC patients were divided into two subgroups based on the median ARS score. The relationship between these risk groups and the aforementioned AGSs was shown in a Sankey diagram (Figure 6F). In general, AGS-B had a higher ARS score than AGS-A (Figure 6G). The expression of 36 ARGs in the two ARS score subgroups was shown in a heatmap (Figure 6H). Figure 6I showed that most ARGs were correlated with the selected genes in the ARS score.

SLC2A3 is a representative gene of ARS for angiogenesis

Among the six genes (TIMP1, FSTL3, CALB2, SLC2A3, FABP4, and TNNT1) in the ARS score, SLC2A3 showed the highest correlation with the ssGSEA score of the angiogenesis pathway (Figure 7A). Therefore, SLC2A3 was selected as the representative gene of ARS for the subsequent analysis. SLC2A3 also exhibited the highest correlation with endothelial cell markers among all ARS genes (Figure 7B). Kaplan-Meier curve indicated that the SLC2A3 was prognostic in CC patients (Figure 7C).

Figure 7 Single-cell analysis for SLC2A3 in colon cancer. A: Correlations between selected 6 genes in the angiogenesis-related score (ARS) model and angiogenesis pathway single-sample gene set enrichment analysis score; B: Correlations between selected 6 genes in the ARS model and marker of the endothelial cell; C: Kaplan-Meier survival curves between low and high expression of SLC2A3 groups; D: Different distinct clusters in normal and colon cancer (CC) tissue by t-SNE analysis; E: The expression levels of SLC2A3 in different cell types; F: The expression levels of SLC2A3 in normal and CC tissue; G: The expression levels of SLC2A3 in different cell types between normal and CC tissue; H: Differences in Gene Set Variation Analysis analysis of metabolism-related pathways in endothelial cells between normal and CC tissues. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Single-cell analysis was conducted to identify the specific cell type with high SLC2A3 expression. There were eight cell types, including epithelial cells, endothelial cells, monocyte cells, T cells, B cells, plasma cells, fibroblast cells, and myeloid cells in CC and normal colon tissue (Figure 7D). Expression of SLC2A3 in different cells, arranged from the highest expression to the lowest expression, was shown in a density plot (Figure 7E and F). Notably, SLC2A3 was nearly absent in the endothelial cells of normal colon tissue but overexpressed in the endothelial cells of CC (Figure 7G). GSVA analysis also indicated the changes in glucose metabolic pathways in CC endothelial cells, supporting the role of SLC2A3, a glucose transporter, in linking metabolism with angiogenesis in CC (Figure 7H).

The overexpression of SLC2A3 in CC tissues was further confirmed by IHC staining (Figure 8A). Colocalization of IF signals of SLC2A3 and CD31 suggested that most endothelial cells (CD31+) in CC tissues were SLC2A3+ (Figure 8B). The biological functions of SLC2A3 in endothelial cells were also revealed by a series of gain-of-function assays. Overexpression of SLC2A3 did not alter cell proliferation (Figure 8C) and tubule formation (Figure 8D). However, overexpression of SLC2A3 promoted endothelial cell migration (Figure 8E).

Figure 8 Immunohistochemistry, and immunofluorescent for SLC2A3 and cell migration, proliferation, and tube formation for SLC2A3 overexpression cells. A: The immunohistochemistry of SLC2A3 in normal and colon cancer (CC) tissue; B: Immunofluorescence experiment demonstrates colocalization of SLC2A3 and CD31 in CC tissue; C: Cell proliferation ability of SLC2A3 overexpression cells was measured using the EdU experiment; D: Migratory ability of SLC2A3 overexpression cells was measured using the transwell migration assay; E: Tube formation assay was used to determine the effect of SLC2A3 overexpression on angiogenesis, and the number of nodes and total length were measured using Image J software. ^bP < 0.01; ^cP < 0.001.

The ARS and tumor characteristics

The ssGSEA scores of many Hallmark pathways, including DNA repair, mTORC1 signaling, and unfolded protein response, were significantly different between low and high ARS score subgroups (Figure 9A). The heatmap showed the predicted abundance of 28 immune cell subtypes in those subgroups (Figure 9B). The expression levels of the six genes in the ARS score were positively associated with the percentages of different immune cells (Figure 9C). The expression of PD-1-related checkpoint genes (PD1, PDL1, PDL2), CTLA4-related checkpoint genes (CTLA4, CD80, CD86), and agonists of T cell activation-related checkpoint genes (OX40, GITR, TNFRSF9, ICOS, CD40, CD27, CD70) in the high ARS score subgroup were higher than in the low ARS score subgroup (Figure 9D-F). The IPS score of samples with a negative expression of PD1 in the high ARS score subgroup was lower than in the low ARS score subgroup, whether the expression of CTLA4 was positive or negative (Figure 9G-J). The TIDE score in the high ARS score subgroup was higher than in the low ARS score subgroup (Figure 9K).

Figure 9 Identification of characterization, the abundance of immune cells, immune checkpoint, and immunotherapy effect between low and high-risk angiogenesis-related signature groups. A: The heatmap showed the differences in HALLMARKER pathway single-sample gene set enrichment analysis scores between the two angiogenesis-related signature (ARS) groups; B: Differences in the abundance of immune cells between two ARS groups; C: Correlations between the abundance of immune cells and selected 6 genes in the ARS model; D: Differences in PD-1, PD-L1, and PD-L2 expression between two ARS groups (P < 0.001); E: Differences in CTLA4, CD80, and CD86 expression between two ARS groups (P < 0.001); F: Differences in OX40, GITR, TNFRSF9, ICOS, CD40, CD27, CD70 expression between two ARS groups (P < 0.001); G-J: Immunopheno score in two groups; K: Tumor immune dysfunction and Exclusion in two ARS groups; L: Survival analyses for two ARS groups in the anti-PD-L1 immunotherapy cohort using Kaplan-Meier curves (IMvigor210 cohort, P = 0.060); M-Q: Relationships between ARS groups and chemotherapeutic sensitivity; R: The proportion of high and low ARS groups in the three modification patterns. Microsatellite stability, blue; microsatellite instability-low red; microsatellite instability-high orange. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Kaplan-Meier analysis showed that the low ARS score subgroup had a marginally significantly better OS than other patients in the IMvigor210 cohort (Figure 9L). For chemotherapeutic response prediction, the high ARS score subgroup generally had a lower IC₅₀ compared to the low ARS score subgroup in most drugs (Figure 9M-Q). Additionally, the high ARS score subgroup had a higher percentage of MSI-H and a lower percentage of MSS and MSI-L compared to other patients (Figure 9R).

The ARS and clinical outcome

To evaluate the prognostic value of the ARS score, we subjected two independent cohorts (the TCGA-COAD cohort and the GSE17536 cohort) to prognostic analysis. The clinical characteristics of patients in the TCGA-COAD and GSE17536 cohorts are shown in Supplementary Table 2. For the TCGA-COAD cohort, Kaplan-Meier analysis showed that patients with a low ARS score had a better OS than those with a high ARS score (P = 0.007, Figure 10A). The area under the curves (AUCs) of 1-, 3-, and 5-year OS of CC patients were 0.667, 0.631, and 0.683 in the TCGA-COAD cohort with good calibration fitness (Figure 10B and C). For the GSE17536 cohort, Kaplan-Meier analysis also verified that patients with a low ARS score had a better OS (P = 0.004, Figure 10D-F).

Figure 10  Construction, Internal, and validation of the angiogenesis-related signature model and the Nomogram. A: Survival analyses for low and high angiogenesis-related signature (ARS) groups in The Cancer Genome Atlas Colon Adenocarcinoma (TCGA-COAD) cohort using Kaplan-Meier curves (P = 0.007); B: Time-dependent receiver operator characteristic (ROC) curves for ARS in the TCGA-COAD cohort; C: Calibration plots for ARS in the TCGA-COAD cohort; D: Survival analyses for low and high ARS groups in the GSE17536 cohort (P = 0.004); E: Time-dependent ROC curves for the ARS in the GSE17536 cohort; F: Calibration plots for the ARS in the GSE17536 cohort; G and H: Univariate and Multivariate Cox regression analysis of clinicopathological features and the ARS; I: Nomogram based on ARS, age, and tumor stage for predicting overall survival; J and K: Time-dependent ROC curves and calibration plot for nomogram in TCGA cohort; L and M: Time-dependent ROC and calibration plot curves for nomogram in GSE17536 cohort. AUC: Area under the curve. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

Univariate and Multivariate Cox regression analyses suggested that the ARS score was an independent risk factor for OS (Figure 10G and H). Other independent prognostic factors, including age, risk score, and stage, were used to construct a nomogram for OS prediction (Figure 10I). The ROC curves showed that the AUCs of the nomogram model at 1, 3, and 5 years were 0.755, 0.775, and 0.747 in the TCGA-COAD cohort, respectively (Figure 10J). The calibration curve showed that the 1, 3, and 5 year nomogram model was consistent with the ideal model, indicating good accuracy (Figure 10K). Similar results were obtained from the GSE17536 cohort (Figure 10L and M). Together, the ARS score is a prognostic factor for CC patients. ^aP < 0.05; ^bP < 0.01; ^cP < 0.001.

DISCUSSION

Angiogenesis and TME have been reported to be associated with the growth rate and biological characteristics of the tumor[24]. We identified two AGSs (AGS-A and AGS-B) with significantly different prognoses, TME, intratumor microbiota, drug sensitivity, and cancer-related pathways ssGSEA scores in patients with CC. Notably, patients in AGS-B, which had a worse prognosis, had a higher immune score and stronger immune cell infiltration but also exhibited high immune escape. This indicates that angiogenesis may be related to the involvement of the TME[25,26]. Our study is consistent with the previous studies that have also reported a correlation between angiogenesis and immune TME was also reported, including gastric cancer, kidney cancer, and breast cancer[9-12]. Immunotherapy in combination with anti-angiogenesis therapy may be considered a promising synergistic strategy in CC treatment. Current clinical studies have shown that immunotherapy combined with antiangiogenic therapy is effective for some cancers, including non-small cell lung cancer[7,27-29]. This article provides additional evidence for the combination of anti-angiogenesis therapy and immunotherapy.

The gut microbiome environment affects tumor outcomes, which might be related to angiogenesis in CC[30,31]. Several studies have shown that the expression of the AF-1 operon in E. coli upregulates the expression of HIF-1α in intestinal epithelial cells, leading to increased expression of VEGF and IL-8. This promotes tumor growth through the up-regulation of angiogenesis[32-34]. However, there have been few systematic studies on the relationship between microbes and angiogenesis. Our study further explored the relationship between angiogenesis and intratumor microorganisms in CC. We found that some pathogenic microorganisms had a larger proportion in AGS-B than AGS-A, including Alistipes, Fusobacterium, and Prevotella. Alistipes, as a potential pathogen, may be closely related to the occurrence and development of CC[35]. A.finegoldii, a species in the genus Alistipes, promotes CC via the IL-6/STAT 3 pathway[35]. Alistipes produce sulfonates, which can act as antagonists of Willebrand factor receptors and inhibit TNF-α and are closely related to the function of the vascular endothelium[36]. They can also play a beneficial role in cancer immunotherapy by regulating the immune microenvironment[37]. Previous studies have suggested that the most significant alteration in CC is Fusobacterium[38-41]. Fusobacterium nucleatum (F. nucleatum) has been reported to induce a proinflammatory immune response and promote carcinogenesis, which was positively associated with IBD and CC[42-44]. Kostic et al[44] reported that F. nucleatum increases tumor multiplicity and selectively recruits tumor-infiltrating myeloid cells, which can promote tumor progression. Borowsky reported that the amount of tissue F. nucleatum is associated with the lower density of stromal memory helper T cells[45]. F. nucleatum also affects the polarization of macrophages[46,47]. Prevotella was reported to be highly enriched in proximal CC[48]. Prevotella can promote the progression of CC and affect the therapeutic effect of FOLFOX[49]. Although intratumor microorganisms have been reported to have a pathogenic role in CC, further investigation is warranted to determine their involvement in angiogenesis.

In our study, we established a risk score (ARS score) based on the AGSs This model consisted of six genes (TIMP1, FSTL3, CALB2, SLC2A3, FABP4, and TNNT1), and was successfully validated as a prognostic factor. Additionally, this model was confirmed as a predictor of therapeutic efficacy in CC. Among these genes, TIMP1[50], FSTL3[51], SLC2A3[52], and TNNT1[53] have been reported to be closely associated with the progression and prognosis of CC. The role of FABP4 in regulating the progression of CC is still controversial[54,55]. Tian et al[54] reported that FABP4 promoted invasion and metastasis of CC by regulating fatty acid transport, while Zhao et al[55] reported that FABP4 inhibited the migration, invasion, and EMT of colon cells. Although these genes have been reported to be related to the progression of CC, the mechanism of their association with angiogenesis in CC needs to be further investigated.

Of these six genes, SLC2A3 showed the highest correlation with the ssGSEA score of angiogenesis. SLC2A3 is involved in glucose transport and is closely associated with tumor metabolic reprogramming[56-58]. Kuang et al[59] reported that targeting SLC2A3 in resistant tumors could enhance the potential of antiangiogenic treatments. However, the relationship between SLC2A3 and angiogenesis in CC is still unknown. We found that SLC2A3 was highly expressed in endothelial cells in CC tissue and overexpression of SLC2A3 promoted the migration of HUVECs. The initial step in angiogenesis is endothelial cell migration[60]. Cellular glycolysis promotes endothelial cell activation, migration, and contraction[61]. Our study also indicated that the glycolytic pathway was significantly enriched in endothelial cells in CC through single-cell analysis, suggesting that SLC2A3 might affect endothelial cell migration by regulating endothelial cell glycolysis, thus affecting the angiogenesis in CC.

While our study provides a novel perspective on angiogenesis in CC, there are some potential limitations that should be acknowledged. Firstly, the data from public databases were obtained retrospectively, and inherent selection bias might affect may eir robustness. Second, due to the lack of microbial information on the species level in the TCMA database, we were unable to further explain the correlation of different microbial species with angiogenesis. Moreover, large-scale prospective cohort studies and systemically functional experimental assays are further needed for an in-depth understanding of the TME and microbiome in relation to the two AGSs.

CONCLUSION

Our study identified two AGSs for CC patients with distinct clinical outcomes, TME, and intratumor microbial compositions. Based on the DEGs between the two subtypes, we further established a reliable ARS score to determine the subtypes and predict the prognosis and therapeutic response in CC patients, which could assist in personalized treatment. In addition, our study further explored the mechanism of SLC2A3 in promoting angiogenesis, which could support the development of potential anti-angiogenic therapeutic targets in the future.