Gene expression of VEGF family members

Using the GTEx database, we initially investigated the expression of VEGF family members in normal human tissues. As shown in Fig. 1A, VEGFA and VEGFB are more abundantly expressed in various normal human tissues than the other members. No family members showed significantly elevated expression levels in normal colon tissue. Using the TCGA-COAD cohort, a pan-cancer analysis revealed that members of the VEGF family exhibited distinct expression patterns in tumor tissues versus normal tissues. VEGFA, VEGFB, VEGFC, VEGFD, and PGF were expressed differently in 13, 17, 16, 19, and 18 tumors and normal samples, respectively (Fig. 1B). The expression of VEGF family members was then analyzed in COAD using paired TCGA-COAD samples (41 cases of the tumor group vs. 41 cases of the normal group). The results indicated that VEGFA, VEGFB, and PGF were highly expressed in tumors, whereas VEGFD expression was low in COAD (Fig. 2A). A total of 98 paired samples from the GSE44076 microarray in the GEO database were utilized for validation analysis, and the results were largely consistent with those of the TCGA-COAD cohort (Fig. 2B, C).

Fig. 1

The mRNA expression of VEGF family in normal and tumor tissues (A) The expression level of VEGF family member in human normal tissues based on GTEx database(B) The expressing levels of VEGF members across diverse cancer types

Fig. 2

Differential mRNA expression of VEGF members in paired samples based on TCGA and GEO database. A Analysis of VEGF members expression in 41 pairs of COAD and adjacent normal tissues. B The heatmap of VEGF members expression in GSE44067 datasets. C Analysis of VEGF members expression in COAD and adjacent normal tissues based on GSE44067 datasets. * P< 0.05; *** P< 0.001

Genetic alterations and gene and protein networks of VEGF family members

Using the cBioportal database, we evaluated the gene mutation status of the VEGF family in 110 colon cancer samples and found that the mutation frequency of VEGF family members varied significantly. VEGFA had the highest mutation frequency (6%), followed by VEGFB and VEGFC (VEGFB 2.9% and VEGFC 2.9%), while VEGFD and PGF were relatively conservative (VEGFD, 1.0% and PGF, 1.0%) (Fig. S1A). Moreover, the VEGF gene mutations did not significantly affect the overall survival of COAD (Fig. S2). Then, we constructed the gene-gene interaction network and protein-protein interaction (PPI) network of the VEGF family using the GeneMANIA and STRING databases, respectively. The gene-gene interaction network uncovered twenty potential VEGF family target genes (Fig. S1B). According to the outcomes of the protein-protein interaction network, there are close interactions between VEGF family genes (PPI enrichment p-value = 4.22e-13) (Fig. S1C). Lastly, correlation analysis based on the TCGA-COAD cohort revealed that all VEGF family members were positively correlated (Fig. S1D).

Survival, prognostic and clinical correlation analysis of the VEGF family in COAD

Using univariate and multivariate Cox regression analyses, the risk factors associated with the prognosis of COAD patients were identified. Age, T stage, N stage, M stage, pathologic stage, CEA level, and VEGFB expression level were identified as risk factors for the prognosis of COAD patients by univariate Cox regression analysis (Fig. 3A, B). The AUC of the ROC curve indicated that VEGFA, VEGFB, VEGFC, and VEGFD had good predictive efficacy in the VEGF family (their respective AUC values were 0.777, 0.736, 0.707, and 0.969). (Fig. 3C). Age, pathologic stage, and the expression level of VEGFB were identified as independent prognostic markers for COAD patients by multivariate Cox analysis (Fig. 3D). High mRNA expression of VEGFB was associated with poor overall survival (OS) (P = 0.008), disease-specific survival (DSS) (P = 0.004), and progression-free survival (PFS) (P = 0.041) in COAD patients as determined by the Kaplan-Meier survival curve (Fig. 3E). In addition, our study demonstrated that high VEGFB expression was associated with adverse clinicopathological characteristics of COAD patients, including advanced tumor stage, increased lymph node and distant metastasis, and poor prognosis (Fig. 3F).

Fig. 3

Survival, prognostic and clinical correlation analysis of the VEGF family in COAD. A ROC analysis of VEGF family members. B Univariate Cox regression analyses of VEGF family and clinical Parameters in COAD. C Multivariate Cox regression analyses of VEGF family and clinical Parameters in COAD. D Representative images of VEGFB expression in COAD tissues and normal controls. Original magnifications 50× and 100×; E Kaplan-Meier survival curves comparing the high and low expression of VEGFB in COAD patients, including overall survival, disease specific survival and progress free interval. F The mRNA expression of VEGFB in COAD based on pathologic stage, N stage, lymphatic invasion, M stage and DSS event

VEGFB family diagnosis and significance in COAD

We developed a nomogram to predict the 1-, 3-, and 5-year prognosis of COAD patients by integrating the independent risk factors (patient age, pathological stage, and VEGFB expression level) identified by multivariate Cox regression analysis (Fig. 4A). The patient’s prognosis becomes direr as the score calculated by the nomogram increases. The calibration curve demonstrated that the constructed nomogram had excellent predictive performance (Fig. 4B).

Fig. 4

Construction and validation of the nomogram. A The nomogram to predict the 1 -, 3 -, and 5-year prognosis of COAD patients; B The calibration curve of nomogram

Function enrichment analysis of VEGFB in COAD

Using the median VEGFB expression level, samples from the TCGA-COAD cohort were separated into VEGFB high and low expression groups. Between the two groups, a total of 7055 DEGs were identified, including 724 up-regulated genes and 6331 down-regulated genes (Fig. 5A, B). The analyses of GO and KEGG were used to determine the biological functions and pathways associated with DEGs. GO enrichment analysis of up-regulated genes revealed that these genes were mainly related to various protein complexes, extracellular matrix and muscle contraction (Fig. 5C). KEGG enrichment analysis revealed that these genes were primarily involved in the regulation of vascular smooth muscle contraction and the cAMP/Wnt signaling pathway (Fig. 5D). GO enrichment analysis of down-regulated genes revealed that they were primarily involved in mRNA binding, receptor activation, DNA complexes, and chromosome and nucleosome assembly (Fig. 5E), whereas KEGG enrichment analysis revealed that they were primarily involved in RNA transport and viral carcinogenesis (Fig. 5F).

Fig. 5

Function enrichment analysis and DNA methylation of VEGFB in COAD. A Differentially expressed genes (DEGs) for high expression of VEGFB vs. low expression of VEGFB in COAD were shown in the volcano plot, with red dots representing significantly up-regulated genes and blue dots representing significantly down-regulated genes in with high expression of VEGFB. B The heatmap exhibits the expression level. C Enrichment analysis for GO term of up-regulated genes. D Enrichment analysis for GO pathway of down-regulated genes

DNA methylation and and tumor mutational burden analysis of VEGFB in COAD

We investigated the methylation status of VEGFB in the COAD using the GSCA database. The results indicated that the methylation level of VEGFB decreased gradually as VEGFB expression increased (r= -0.33, P 0.001). (Fig. 6A). Utilizing the MethSurv database, we then drew the methylation site expression heat map of VEGFB, yielding a total of seven methylation sites (Fig. 6B). The methylation sites cg05492845 and cg18872604 are significantly associated with colon cancer patients’ prognosis (Fig. 6C). The preceding findings suggest that VEGFB methylation may play a role in the development of COAD. The correlation analysis between VEGFB expression level and tumor mutation burden showed that with the increase of VEGFB expression level, the level of tumor mutation burden decreased, and there was a negative correlation between them (Fig. 6D).

Fig. 6

The DNA methylation analysis of VEGFB in COAD. A Correlation between the VEGFB mRNA expression and DNA methylation levels in COAD. B The heatmap of DNA methylation of VEGFB in COAD obtained from MethSurv database. C The prognostic value of DNA methylation of VEGFB in COAD with different CpG sites. D The relationship between VEGFB expression levels and tumor mutational burden

Immune-related analysis of VEGFB in COAD

Analyzing and comprehending the tumor immune microenvironment will improve immunotherapy’s efficacy. We therefore investigated the connection between VEGFB and the immune microenvironment. There were statistically significant differences between the VEGFB high expression group and the low expression group in terms of immune score, stromal score, and estimate score (Fig. 7A). We hypothesize that VEGFB may exert its biological function by regulating the immune system. Next, we analyzed the difference in immune cell composition between the group with high VEGFB expression and the group with low VEGFB expression. In the group with a high level of VEGFB expression, the proportion of regulatory T cells and M0 macrophages was greater. T cells CD4 memory activated, T cells gamma delta, resting NK cells, activated Dendritic cells, and Eosinophils were even less prevalent in the VEGFB high expression group (Fig. 7B). In addition, correlation analysis between VEGFB and immune cells revealed that VEGFB was positively correlated with regulatory T cells (Tregs) and Macrophages and negatively correlated with eosinophils, dendritic cell activation, T cell gamma delta, NK cells resting, T cell CD4 memory activated, and T cell CD4 memory resting (Fig. 7C, D, E). The abnormal expression and function of immune checkpoint molecules is one of the primary causes of the development and occurrence of numerous tumors. Consequently, we conducted a correlation analysis to identify the immune checkpoints closely associated with VEGFB. Included in the analysis are the following checkpoints: ADORA2A, BTLA, BTNL2, C10orf54, CD160, CD200, CD200R1, CD244, CD27, CD274, CD276, CD28, CD40, CD40LG, CD44, CD48, CD70, CD80, CD86, CTLA4, HAVCR2, HHLA2, ICOS, ICOSLG, IDO1, IDO2, KIR3DL1, LAG3, LAIR1, LGALS9, NRP1, PDCD1, PDCD1LG2, TIGIT, TMIGD2, TNFRSF14, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF14, TNFSF15, TNFSF18, TNFSF4, TNFSF9, and VTCN1. VEGFB was positively correlated with the immune checkpoints TNFRSF4, TNFRSF8, and TNFRSF18 (Fig. 8A–C).

Fig. 7

Tumor microenvironment and immune-related analysis of VEGFB in COAD. A The relationship between VEGFB expression and tumor microenvironment. B The difference of the levels of 22 TIICs between high and low VEGFB expression groups. C The relationship between the abundance of 22 immune cells and VEGFB mRNA expression. The transcription level of VEGFB was significantly positively associated with the levels of Tregs cells and macrophages infiltration in COAD tissues. D Correlation analysis between VEGFB expression level and Tregs infiltration. E Correlation analysis between VEGFB expression level and macrophage cells infiltration

Fig. 8

Correlation between the VEGFB expression and immune checkpoints in COAD. A Radar chart evaluating the relationship of immune checkpoint molecules in COAD. B The co-expression heatmap of VEGFB and the expression of TNFRSF4, TNFRSF8 and TNFRSF18. C Scatter diagrams showed that TRPV3 mRNA expression was significantly positively correlated with the expression of TNFRSF4, TNFRSF8 and TNFRSF18. ***P < 0.001

SHNG17-miR-375-VEGFB regulatory axis construction

We attempted to determine the regulatory network upstream of VEGFB’s miRNA. Through the “miRNA-mRNA” module of the StarBase database, 57 potential miRNA genes were identified. Given that miRNAs typically inhibit mRNA expression, the target miRNA should have a negative correlation with VEGFB, be lowly expressed in tumor tissues, and be associated with a favorable prognosis. The correlation analysis revealed that, among the 57 candidate miRNAs, 17 were significantly negatively correlated with VEGFB (Table S1). Expression analysis revealed that six microRNAs were significantly overexpressed in normal colon tissues (Fig. S3). We then conducted a survival analysis to narrow down the candidate miRNAs. Only miR-375 was significantly associated with a favorable prognosis for COAD patients, according to the results. (Fig. S4). We concluded that miR-375 is the most likely upstream miRNA of VEGFB in COAD based on correlation, expression, and survival analyses (Fig. 9A, B). We analyzed its association with immune infiltrating cells and immune checkpoints in greater detail. There is a significant positive correlation between hsa-miR-375 and VEGFB (TNFRSF4, TNFRSF8, and TNFRSF18) and a negative correlation between it and macrophages and immune checkpoints (Fig. 9C-E). The aforementioned results confirmed miR-375 as the upstream miRNA of VEGFB. We utilized the same method to examine the upstream lncRNAs of miR-375 in COAD. Through the “miRNA-lncRNA” module, eight lncRNAs were extracted initially. lncRNAs should be negatively correlated with miRNAs and positively correlated with mRNAs, given that they compete with miRNAs for mRNA binding. Only SNHG17 satisfies the aforementioned two requirements (Fig. 10A). According to the correlation analysis, SNHG17 had a negative correlation with miR-375 and a positive correlation with VEGFB (Fig. 10B). Expression and survival analysis confirmed the reliability of SNHG17 as an upstream lcnRNA. The results demonstrated that SNHG17 was highly expressed in COAD tissues and was associated with a dismal prognosis for patients (Fig. 10C, D). Finally, we constructed a predictive model for the SNHG17-miR-375-VEGFB axis in relation to COAD prognosis and progression (Fig. 10E).

Fig. 9

Identification of miR-375 as a potential upstream miRNA of VEGFB in COAD. SHNG17-miR-375-VEGFB Regulatory Axis Construction. A miR-375 was downregulated in COAD. ***P < 0.001. B COAD patients had a better survival with high expression of miR-375. C The relationship between the abundance of 24 immune cells and miR-375 expression. D The level of Macrophages was significantly down-regulated in the miR-375 high-expression group compared to the low-expression group. ***P < 0.001. E Scatter diagrams showed that miR-375 expression was significantly negatively correlated with the expression of TNFRSF4, TNFRSF8 and TNFRSF18

Fig. 10

Identification of SNHG17 as potential upstream lncRNAs of has-miR-375 in COAD (A) Correlation of VEGFB and miR-375 with candidate lncRNAs. B SNHG17 was negatively correlated with miR-375 while positively correlated with VEGFB. C SNHG17 were upregulated in COAD. ***P < 0.001. D High expression of SNHG17 predicted a worse survival of COAD patients. E The model of SNHG17-miR-375-VEGFB axis in carcinogenesis of COAD

VEGFB enhance the proliferation and migration of colon cells

The expression level of VEGFB in colon cancer tissues was significantly higher than that in adjacent normal tissues, both RT-qPCR and immunohistochemistry confirmed (Figs. 11A, B). We further transfected SW480 cells with siRNA and successfully down-regulated the expression of VEGFB, and SiRNA NC was also transfected into SW480 cells (Fig. 11C). The results of CCK-8 assay and clonal formation assay showed that down-regulation of VEGFB expression reduced cell viability and proliferation (Figs. 11D, E). Scratch assay was used to examine the effect of VEGFB on the migration ability of sw480 cells (Fig. 11F), and the results showed that down-regulation of VEGFB expression reduced the migration of SW480 cells.

Fig. 11

VEGFB is highly expressed in colon cancer and its activity is necessary for tumor cell proliferation and migration. A The expression of VEGFB in tumor tissues was significantly higher than that in adjacent normal tissues verified by RT-q-PCR; B Immunohistochemistry confirmed that the expression of VEGFB in tumor tissues was significantly higher than that in adjacent normal tissues; C Treatment of SW480 cells with control siRNA or VEGFB-targeting siRNA confirmed that VEGFB was successfully inhibited by RT-q-PCR analysis. D The CCK-8 showed that inhibition of VEGFB expression could reduce cell viability. E The colony arrangement assay confirmed that the expression of VEGFB was inhibited and the proliferation ability of cells was reduced. F The scratch assay confirmed that the expression of VEGFB was decreased and the cell migration ability was decreased