Bioinformatics analysis reveals VEGFC’s prognostic significance in head and neck squamous cell carcinoma and its association with immune cell infiltration

Introduction

Head and neck squamous cell carcinoma (HNSCC), an aggressively invasive malignant tumor, ranks as the sixth most common cancer worldwide, with prognoses that are often poor (1,2). This serious situation is primarily due to the intertwined effects of two core challenges: firstly, a scarcity of early and effective diagnostic and prognostic biomarkers, leading to many patients being diagnosed at an advanced stage or with lymphatic or distant metastasis already present (3,4); secondly, the lack of a profound understanding of the molecular mechanisms that influence the survival rate and prognosis of HNSCC patients (5). Therefore, the search for novel biomarkers to enable early diagnosis, risk assessment, and prognosis prediction of HNSCC has emerged as one of the hotspots in current oncology research. To date, researchers have explored the application of various biomarkers in the prognosis assessment of HNSCC, including epidermal growth factor receptor (EGFR) (6), CC chemokine receptor 4 (CCR4) (7), and certain non-coding RNA (ncRNAs) (2,8). These markers have to some extent enhanced our understanding of the biological characteristics of HNSCC. However, in clinical practice, the sensitivity and specificity of many known markers are often insufficient to meet the demands of precision medicine, particularly in terms of accurately predicting individual patient outcomes.

In this context, vascular endothelial growth factor C (VEGFC), as a key factor in promoting angiogenesis and lymphangiogenesis, has become a focus of interest for researchers (9). The protein encoded by VEGFC gene is a member of the platelet-derived growth factor/vascular endothelial growth factor (PDGF/VEGF) family. Studies have shown that the binding of VEGFC to its receptor VEGFR3 can activate the VEGFC/VEGFR3 pathway, which is considered a major inducer of lymphangiogenesis that is involved in lymphatic metastasis (9,10). The prerequisite for tumor lymphatic metastasis is the induction of primary lymphatic vessels and neovascularization by tumor cells invading the tumor stroma. It is evident that VEGFC plays a critical role in promoting tumor immune evasion, enhancing tumor growth (11). At present, multiple studies have found VEGFC to be widely involved in the development of various cancers such as gastric cancer (12), medulloblastoma (13), and pancreatic neuroendocrine tumors (14), and it is closely related to tumor immunity. In HNSCC, some findings have reported VEGFC’s contribution to the growth and motility of HNSCC cells (15), and its involvement in regulating the invasion of HNSCC cells in vitro, as well as its positive correlation with recurrence and lymph node metastasis beyond the primary lesion of HNSCC (16). However, there is a currently lack of research on whether VEGFC indicates the prognosis of HNSCC and deeply explores its association with immune infiltration and immune regulatory molecules.

Therefore, this study aims to analyze comprehensively the expression and prognostic value of VEGFC in HNSCC through multiple cancer databases, including Tumor Immune Estimation Resource 2.0 (TIMER2.0), Gene Expression Profiling Interactive Analysis (GEPIA), University of ALabama at Birmingham CANcer data analysis Portal (UALCAN), The Cancer Genome Atlas (TCGA), Search Tool for the Retrieval of Interacting Genes/Proteins (STRING), Tumor-Immune System Interaction Database (TISIDB), and Gene Expression Omnibus (GEO), and explore its correlation with immune cell infiltration and immune regulation-related genes. This will provide a theoretical basis for clinical diagnosis and treatment of HNSCC. We present this article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-834/rc).

Methods

Database

This study primarily utilized several online databases for direct analysis, including TIMER 2.0 (http://timer.cistrome.org/), UALCAN (http://ualcan.path.uab.edu/), STRING (https://www.string-db.org/), GEPIA (http://gepia2.cancer-pku.cn/), and TISIDB (http://cis.hku.hk/TISIDB/). For prognosis analysis, expression profiles of 503 HNSCC tumor tissues and 44 normal control tissues were retrieved from TCGA (https://cancergenome.nih.gov/), along with corresponding clinical information such as age, gender, race, tumor stage, and nodal status. Clinical information and survival outcomes were extracted for analysis including univariate and multivariate Cox regression analysis, and prognostic nomogram development. The training and validation cohorts for the nomogram were derived from the TCGA dataset, with appropriate stratification to ensure the model’s generalizability and accuracy. The calibration curve was used to assess the accuracy of the nomogram. An ideal calibration curve should be close to the 45-degree diagonal line, indicating that the model’s predicted results are consistent with the observed results. For gene expression validation analysis, the GSE6631 dataset from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) was utilized, which includes whole-genome data from 22 HNSCC tissues and 22 adjacent normal tissues.

Materials and reagents

The reverse transcription quantitative polymerase chain reaction (RT-qPCR) experimental validation was performed using the immortalized human nasopharyngeal epithelial cell line NP69SV40T and the nasopharyngeal carcinoma cell line CNE2. These cell lines were obtained from PuNoSai Life Technology Co., Ltd. (Wuhan, China). The total RNA extraction kit was purchased from Solaibio Technology Co., Ltd. (Shanghai, China). The first-strand cDNA synthesis kit was purchased from BiyunTian Biotechnology Co., Ltd. (Shanghai, China); the Hiff™ qPCR SYBR^® Green Master Mix was purchased from Yisheng Biotechnology Co., Ltd. (Shanghai, China).

Expression analysis of VEGFC gene in various common tumor tissues

The expression levels of VEGFC in various common tumor tissues were analyzed using TIMER 2.0 and GEPIA databases. Within the TIMER 2.0 database, the “Gene_DE” analysis, function under the “Exploration” module, was employed to compare VEGFC expression between various common tumor tissues and adjacent normal tissues. The distribution of gene expression levels were displayed using box plots. Additionally, the GEPIA database was utilized to analyze VEGFC expression across different tumor types relative to control tissues, with results visualized through scatter plots.

Expression analysis of VEGFC gene in HNSCC and normal control tissues

The expression levels of VEGFC in HNSCC and normal control tissues were analyzed using the UALCAN database, followed by validation of expression difference through the GEPIA database. The “TCGA analysis” module within UALCAN facilitated the retrieval of VEGFC expression data in both HNSCC and normal control tissues. In the GEPIA database, the TCGA and GTEX datasets were selected for analysis. Differential expressed genes were identified using the criteria of “|Log₂FC| ≥1” and “P<0.01”. The results of differential expression were visualized using box plots.

Relationship analysis between VEGFC gene and clinical pathological features

The correlations between VEGFC and various clinical pathological features of HNSCC, including tumor stage, gender, age, race, and nodal metastasis status, was examined individually using the “Expression” analysis module of the UALCAN database.

Relationship analysis between VEGFC gene and clinical prognosis of HNSCC

The GEPIA database was utilized to analyze the impact of VEGFC expression on overall survival (OS) and disease-free survival (DFS) in HNSCC patients. Patients were divided into high and low expression groups based on the quartile values of gene expression, and Kaplan-Meier (KM) survival curves were generated. RNA-seq data in transcripts per million (TPM) format and clinical information for HNSCC were obtained from the TCGA database and subjected to univariate and multivariate Cox regression analysis to evaluate VEGFC as a potential independent prognostic factor.

The survival R package was utilized for proportional hazard assumption testing and to create visualizations through univariate and multivariate Cox regression forest plots. The rms package was employed to map the predicted results of the regression model onto scales, allowing for the visualization of the impact of multiple variables and the quantification of a nomogram-related model designed to predict the 1-, 3-, and 5-year OS occurrence rates in HNSCC patients. The ggplot2 package was also utilized to generate risk factor plots, which assessed the influence of risk factors on the survival outcomes of HNSCC patients. Finally, calibration curves were plotted to assess the consistency between the predicted and actual survival rates at different time points (1, 3, and 5 years), thereby evaluating the accuracy and reliability of the predictive model.

Functional enrichment analysis of VEGFC interacting proteins and co-expressed genes

The STRING database was queried to identify proteins interacting with VEGFC, using a highest confidence threshold of 0.9 and limiting the interaction to a maximum of 5 partners, while other parameters were left at their default settings. In the GEPIA database, the “Similar Genes Detection” analysis module was then employed to retrieve genes that demonstrate correlation with VEGFC expression in HNSCC. Subsequent analysis of gene-gene correlation was conducted using GEPIA’s “Correlation Analysis” module with Pearson correlation selected as the method, to assess the correlation between co-expressed genes and VEGFC, with results depicted through scatter plots. Additionally, functional enrichment analysis for the interacting proteins and co-expressed genes of VEGFC was performed using the R package clusterProfiler, which facilitated Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. The top 20 enriched functions were displayed in a bubble plot format.

Correlation analysis between VEGFC expression abundance and immune cell infiltration in HNSCC

The “Gene” module within the TIMER 2.0 database was utilized to examine the relationship between the abundance of VEGFC expression in HNSCC and the infiltration proportions of six distinct immune cells types, including B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells. Additionally, the association between these immune cell types and tumor purity was analyzed. The results of this correlation analysis were visualized using scatter plots.

Correlation analysis between VEGFC expression abundance and tumor-infiltrating lymphocytes (TILs), immune modulators, and inflammatory chemokines in HNSCC

The TISIDB database was utilized to search through the “Lymphocyte”, “Immunomodulator”, and “Chemokine” analysis modules, examining the correlation between VEGFC expression in HNSCC and 28 types of TILs, including immune-related molecules features, three categories of immune modulators [comprising immunoinhibitors, immunostimulators, and major histocompatibility complex (MHC) class-related molecules], as well as the correlation with chemokines and their receptors.

GEO dataset of VEGFC gene expression and RT-qPCR validation

The GSE6631 dataset was obtained from the GEO database, which included files such as “Series Matrix File(s)” and the “GPL8300” platform annotation file. Expression levels and sample groupings were derived from the probe matrix file, with gene names matched to probes based on the annotation file. Following data organization, column plots were generated using GraphPad Prism 8.0 software. Additionally, total RNA was extracted from the in vitro cultured immortalized human nasopharyngeal epithelial cell line NP69SV40T and the nasopharyngeal carcinoma cell line CNE2. The extracted RNA underwent reverse transcription to synthesize the first-strand cDNA, followed by RT-qPCR experimental validation. The PCR amplification protocol consisted of initial denaturation at 95 ℃ for 300 s, followed by 40 cycles of denaturation at 95 ℃ for 5 s, and annealing/extension at 60 ℃ for 30 s. The dissolution curve reaction procedure was carried out according to the recommended protocol of the Roche LightCycler96 RT-qPCR instrument, and the relative quantification was calculated using the 2^-ΔΔCt method. Primer sequences are shown in Table 1.

Statistical analysis

Online databases typically applied default statistical methods for analyses. For instance, the TIMER 2.0 database used the Wilcoxon test to determine the significance of VEGFC between differences between various common tumor tissues and adjacent normal tissues. The GEPIA database employed one-way analysis of variance (ANOVA) to compare VEGFC expression levels between HNSCC and normal control tissues. Pearson correlation was utilized for analyzing the correlation between VEGFC and its co-expressed genes, as well as for immune cell infiltration. In the TISIDB database, Spearman correlation was applied to assess the correlation between VEGFC and TILs, immune modulators, and inflammatory chemokines. The Logrank test was used to evaluate the comparison of survival rates between groups with high and low VEGFC expression. The Welch’s t-test was utilized to compare the expression levels of VEGFC messenger RNA (mRNA) between nasopharyngeal carcinoma cells and normal nasopharyngeal epithelial cells. All statistical analyses were conducted using two-tailed tests, with P<0.05 indicated statistically significant differences.

Ethical statement

The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Results

Expression of VEGFC gene in various common tumor tissues

Analyses conducted in both the TIMER 2.0 and GEPIA database consistently demonstrated high expression levels of VEGFC gene across multiple cancer tissues. Notably, there was a significant upregulation of VEGFC expression (P<0.05) observed in several types of cancer, including diffuse large B-cell lymphoma (DLBC), pancreatic adenocarcinoma (PAAD), stomach adenocarcinoma (STAD), thymoma (THYM), and HNSCC, among others (Figure 1).

Figure 1 Differential expression of VEGFC in tumor and normal tissues. (A) VEGFC gene expression in various common tumor tissues and normal tissues in the TIMER 2.0 database; (B) VEGFC gene expression in various common tumor tissues and normal tissues in the GEPIA database. *, P<0.05; **, P<0.01; ***, P<0.001. “T” and “N” in Figure 1B represent tumor tissue and normal control tissue, respectively. VEGFC, vascular endothelial growth factor C; TPM, transcripts per million; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; DLBC, diffuse large B-cell lymphoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSCC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LAML, acute myeloid leukemia; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; THYM, thymoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma; GEPIA, Gene Expression Profiling Interactive Analysis.

Expression analysis of VEGFC gene in HNSCC and normal control tissues

Data from the UALCAN and GEPIA databases indicated that VEGFC was significantly upregulated compared with normal control tissues (P<0.05) (Figure 2).

Figure 2 Differential expression of VEGFC gene in HNSCC and normal control tissues. (A) Expression differences of VEGFC gene in HNSCC and normal control tissues in the UALCAN database. (B) Expression differences of VEGFC gene in HNSCC and normal control tissues in the GEPIA database. In the figure, “T” represents tumor tissue, “N” represents normal control tissue, and “*” indicates P<0.05. TPM, transcripts per million; HNSCC, head and neck squamous cell carcinoma; VEGFC, vascular endothelial growth factor C; UALCAN, University of ALabama at Birmingham CANcer data analysis Portal.

Relationship between VEGFC gene and clinical pathological features

Analysis from the UALCAN database revealed that there were statistically significant differences (P<0.05) in the expression levels of VEGFC gene among different stages of HNSCC patients, different genders, age groups, ethnicities, numbers of axillary lymph node metastasis, and degrees of differentiation compared to normal control groups (Figure 3).

Figure 3 VEGFC genes expression across clinicopathological features of HNSCC. (A) Relative expression of VEGFC mRNA in HNSCC patients of different stages. (B) Relative expression of VEGFC mRNA in HNSCC patients of different gender. (C) Relative expression of VEGFC mRNA in HNSCC patients of different age. (D) Relative expression of VEGFC mRNA in HNSCC patients of different race. (E) Relative expression of VEGFC mRNA in HNSCC patients with axillary lymph node metastasis. (F) Relative expression of VEGFC mRNA in HNSCC patients of different degrees of differentiation. *, P<0.05; **, P<0.01; ***, P<0.001. VEGFC, vascular endothelial growth factor C; HNSCC, head and neck squamous cell carcinoma.

Correlation of VEGFC gene with clinical prognosis in HNSCC

The effect of VEGFC gene expression on OS and DFS in HNSCC patients, as demonstrated in the GEPIA database, indicated that patients with high VEGFC expression exhibited significantly lower OS and DFS rates compared to those with low expression (P<0.05) (Figure 4A,4B). Univariate Cox regression analysis of clinical characteristics such as age, gender, pathological stage and grade showed that tumor stages T3&T4 [P<0.001, hazard ratio (HR): 1.934, 95% confidence interval (CI): 1.413–2.649], N2&N3 stages (P<0.001, HR: 2.296, 95% CI: 1.686–3.127), Stage III&IV (P=0.003, HR: 1.839, 95% CI: 1.236–2.737), lymph node invasion (P=0.002, HR: 1.708, 95% CI: 1.217–2.397), and high VEGFC expression (P=0.01, HR: 1.403, 95% CI: 1.073–1.835) are prognostic risk factors for HNSCC (Figure 4C). Further multivariate Cox regression analysis confirmed that N2&N3 stages (P=0.006, HR: 1.840, 95% CI: 1.192–2.840) and high VEGFC expression (P=0.04, HR: 1.499, 95% CI: 1.026–2.188) as independent prognostic factors for HNSCC (Figure 4D). The prognostic nomogram, constructed based on clinical features and VEGFC expression levels, showed good predictive ability for 1-, 3-, and 5-year survival rates in HNSCC patients, indicating that the survival probability of the high-risk group is significantly lower than that of the low-risk group at these intervals (Figure 4E). The risk factor graph further demonstrates the relationship between VEGFC expression and survival outcomes in HNSCC patients, indicating a significant increase in disease progression and metastasis risk with increased VEGFC expression (Figure 4F). The calibration curve results show good consistency between the prognostic nomogram, based on clinical features and VEGFC expression levels, and the actual survival rates for 1-, 3-, and 5-year intervals in HNSCC patients, with the calibration curves closely approaching the ideal 45-degree diagonal line (Figure 4G-4I).

Figure 4 Prognostic value of VEGFC gene in HNSCC. (A) Impact of VEGFC gene expression on overall survival of HNSCC patients. (B) Impact of VEGFC gene expression on disease-free survival of HNSCC patients. (C) Univariate Cox regression forest plot of VEGFC gene with clinical characteristics. (D) Multivariate Cox regression forest plot of VEGFC gene with clinical characteristics. (E) Column chart representing the prediction model. (F) Risk factor diagram of the prediction model. (G-I) Calibration plots of the prediction model for 1-, 3-, and 5-year overall survival prediction, respectively. VEGFC, vascular endothelial growth factor C; HR, hazard ratio; CI, confidence interval; HNSCC, head and neck squamous cell carcinoma.

Correlation analysis and functional enrichment of VEGFC interacting proteins and co-expressed genes

Utilizing the STRING database, a protein-protein interaction (PPI) network was constructed, identifying 50 proteins that closely interact with VEGFC (Figure 5A). Additionally, the GEPIA database was employed to identify 100 co-expressed genes in HNSCC. Correlation analysis of these co-expressed genes with VEGFC showed a strong correlation, with genes such as ITGA3, NT5E, PXN, LAMA3, CAV1, TNFRSF12A, LIMA1, and LAMC2 showing particularly high correlation coefficients (P<0.05, and R>0.6) (Figure 5B).

Figure 5 VEGFC interacting proteins and co-expressed genes. (A) Protein-protein interaction network of VEGFC. (B) Scatter plot depicting the correlation between VEGFC gene and select co-expressed genes in HNSCC (representative genes shown). VEGFC, vascular endothelial growth factor C; TPM, transcripts per million.

GO and KEGG enrichment analysis of 50 VEGFC-interacting protein-coding genes and 100 co-expressed genes (Figure 6) revealed that they mainly participate in biological processes such as peptide-tyrosine phosphorylation, epithelial cell migration, positive regulation of the MAPK cascade, and positive regulation of protein kinase B signaling, among others. They also perform molecular functions like protein tyrosine kinase activity, transmembrane receptor protein kinase activity, cell adhesion molecule binding, growth factor binding, and protein heterodimerization activity. Moreover, these genes play roles in several pathways, including the PI3K/Akt signaling pathway, MAPK signaling pathway, focal adhesion, Rap1 signaling pathway, Ras signaling pathway, proteoglycans in cancer, and central carbon metabolism in cancer.

Figure 6 GO and KEGG enrichment analysis of VEGFC interacting and co-expressed genes. (A) Top 20 enriched KEGG pathways. (B) Top 20 biological processes. (C) Top 20 cellular components. (D) Top 20 molecular functions. KEGG, Kyoto Encyclopedia of Genes and Genomes; HIF, hypoxia inducible factor; EGFR, epidermal growth factor receptor; GO, Gene Ontology; VEGFC, vascular endothelial growth factor C.

Correlation of VEGFC gene abundance with immune cell infiltration, TILs, immune modulators, and inflammatory chemokines and their receptors

Analysis from the TIMER 2.0 database regarding the gene expression abundance of VEGFC in HNSCC revealed a significant correlation with the proportions of several immune cell infiltrates. Specifically, VEGFC expression was found to correlate with B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, and dendritic cells (all with P<0.05), while no significant correlation was observed with macrophage infiltration (P>0.05) (Figure 7A). Further Spearman correlation analysis from the TISIDB database between the expression abundance of VEGFC gene in HNSCC and various immune features, including 28 types of TILs, three categories of immune modulators (comprising immunoinhibitors, immunostimulators, and MHC class-related molecules), and chemokines and their receptors, showed several molecules with relatively high correlation to VEGFC expression (R>0.4 for all). These included TILs-related immune feature molecules such as Tcm CD4 and Tgd; immunoinhibitor-related molecules like PDCD1LG2, TGFB1, and TGFBR1; immunostimulator-related molecules such as NT5E, CD276, and PVR; MHC class-related molecules including TAP2, HAL-C, HAL-B, HAL-A, HAL-E, HAL-F, and B2M; and chemokines and receptor-related molecules such as CXCL12 and CXCR1 (Figure 7B-7G).

Figure 7 Correlation of VEGFC gene abundance with immune cell infiltration, TILs, immune modulators, and inflammatory chemokines and their receptor. (A) Correlation of VEGFC gene abundance with immune cell infiltration, depicted through scatter plots for tumor purity, B-cells, CD4⁺ T-cells, CD8⁺ T-cells, macrophages, neutrophils, and dendritic cells infiltration, in that order. (B) Correlation between TILs and VEGFC gene expression; (C-E) Correlation between immunostimulators, including immunostimulators, immunoinhibitors, and MHC class-related molecules with VEGFC gene expression. (F,G) Correlation between chemokines and their receptors with VEGFC gene expression. VEGFC, vascular endothelial growth factor C; TPM, transcripts per million; TILs, tumor-infiltrating lymphocytes; MHC, major histocompatibility complex; ACC, adrenocortical carcinoma; BLCA, bladder urothelial carcinoma; BRCA, breast invasive carcinoma; CESC, cervical squamous cell carcinoma and endocervical adenocarcinoma; CHOL, cholangiocarcinoma; COAD, colon adenocarcinoma; ESCA, esophageal carcinoma; GBM, glioblastoma multiforme; HNSCC, head and neck squamous cell carcinoma; KICH, kidney chromophobe; KIRC, kidney renal clear cell carcinoma; KIRP, kidney renal papillary cell carcinoma; LGG, brain lower grade glioma; LIHC, liver hepatocellular carcinoma; LUAD, lung adenocarcinoma; LUSC, lung squamous cell carcinoma; MESO, mesothelioma; OV, ovarian serous cystadenocarcinoma; PAAD, pancreatic adenocarcinoma; PCPG, pheochromocytoma and paraganglioma; PRAD, prostate adenocarcinoma; READ, rectum adenocarcinoma; SARC, sarcoma; SKCM, skin cutaneous melanoma; STAD, stomach adenocarcinoma; TGCT, testicular germ cell tumors; THCA, thyroid carcinoma; UCEC, uterine corpus endometrial carcinoma; UCS, uterine carcinosarcoma; UVM, uveal melanoma.

Verification of VEGFC gene expression

Analysis of gene chip data from 22 patients, comparing paired HNSCC tumor tissues with normal control tissues using GSE6631 dataset, revealed a significant elevation within the tumor tissues (P<0.01) (Figure 8A). This finding was further substantiated by RT-qPCR validation in the human immortalized nasopharyngeal epithelial cell line NP69SV40T and the nasopharyngeal carcinoma cell line CNE2. The results demonstrated a marked upregulation of VEGFC in the nasopharyngeal carcinoma cells when compared to the normal nasopharyngeal epithelial cells (Figure 8B) (P<0.01). These experimental outcomes corroborate the bioinformatics analysis.

Figure 8 Validation of VEGFC gene expression using GEO dataset and RT-qPCR. (A) Relative expression levels of VEGFC gene in normal control tissues and HNSCC tumor tissues. (B) RT-qPCR validation of relative VEGFC mRNA expression in the human immortalized nasopharyngeal epithelial cell line NP69SV40T and the nasopharyngeal carcinoma cell line CNE2. **, P<0.01. VEGFC, vascular endothelial growth factor C; HNSCC, head and neck squamous cell carcinoma; GEO, Gene Expression Omnibus; RT-qPCR, reverse transcription quantitative polymerase chain reaction.

Discussion

Tumors growth and metastasis are primarily facilitated by lymphangiogenic factor, such as VEGFC and its receptor/co-receptor systems (13). The development of lymphatic networks depends on VEGFC, which acts as a central regulator in this process. It modulates the immune system, aiding tumor cells in evading immune surveillance and thus fostering tumor progression (17). A previous study has identified a cascade effect initiated by molecules within the VEGFC/D-VEGFR3/NRP2 axis during lymphangiogenesis and lymphatic metastasis (18). This cascade mediates the differentiation and maturation of lymphatic endothelial cells (LECs), ultimately facilitating tumor cell chemotaxis, migration, invasion, and metastasis (18). Despite the clear role of VEGFC in lymphatic regulation, the specific molecular mechanisms underlying its involvement in tumorigenesis and metastasis remain under investigation. VEGFC has been implicated in the activation of the Hedgehog signaling pathway and the epithelial-to-mesenchymal transition (EMT) in various human diseases. For instance, VEGFC secreted by breast cancer cells can activate GLI signaling, thereby promoting paracrine-mediated proliferation, migration, and invasion of breast cancer epithelial cells (19). Additionally, the VEGF-C/VEGFR-3 signaling axis has been shown to enhance metastasis by promoting lymphangiogenesis and angiogenesis in a range of tumors such as leukemia, mesothelioma, and Kaposi’s sarcoma (20). VEGFC is considered a potential diagnostic biomarker and a molecular target for therapeutic intervention. This study indicates that VEGFC is highly expressed in HNSCC and is closely associated with HNSCC tumor stage, patient gender, age, race, etc. Moreover, HNSCC patients with elevated VEGFC expression exhibit significantly lower OS and DFS rates compared to those with lower expression levels. Univariate and multivariate Cox regression analyses reinforce VEGFC as an independent risk factor for HNSCC, highlighting its prognostic significance.

Given VEGFC’s crucial role in promoting lymphatic metastasis, this study specifically observed VEGFC expression in patients with different numbers of axillary lymph node metastases. The results show significant differences in VEGFC expression between patients with axillary lymph node metastasis and the normal control group. However, no discernible pattern was observed in VEGFC expression across groups with different extents of axillary lymph node metastases, with the N2 and N3 groups showing an even smaller difference compared to the N0 and N1 groups. This observation underscores the complexity of VEGFC’s role in HNSCC.

The upregulation of VEGFC expression is linked to axillary lymph node metastasis, suggesting a promotive role in metastasis. Yet, the tumor microenvironment, with its myriad of interactive molecules and immune-regulatory genes, may influence the functional impact of VEGFC. Further research is warranted to elucidate these interactions and their implications for VEGFC’s role in HNSCC. Our analysis identified several genes that are closely associated with VEGFC in HNSCC, hinting at a possible cooperative regulation in tumorigenesis and progression. Notably, ITGA3, NT5E, PXN, LAMA3, CAV1, TNFRSF12A, LIMA1, and LAMC2 exhibit high correlation with VEGFC. Integrin subunit alpha 3 (ITGA3), a major extracellular matrix mediator, has been reported to interact with the VEGFR3 receptor, influencing endothelial cell migration and proliferation in Kaposi’s sarcoma (21). ITGA3’s involvement in the development of various cancers, including pancreatic cancer (22), glioblastoma (23), and esophageal squamous cell carcinoma (24), its role in HNSCC radioresistance (25), underscores its multifaceted influence in oncogenesis. Extracellular 5'-nucleotidase (NT5E), an immune modulator, is integral to immune regulation and encodes CD73, a key enzyme in adenosine production from extracellular ATP. Adenosine, a potent factor in immune evasion and angiogenesis, is critical for the growth of HPV-negative HNSCC (26). LIM domain and actin binding 1 (LIMA1) is implicated in cytoskeletal dynamics and cell motility, potentially contributing to tumor proliferation, invasion, and migration (27). The functional enrichment analysis of VEGFC and its related genes in this study revealed their involvement in pathways such as the PI3K/Akt signaling pathway, MAPK signaling pathway, focal adhesion, Rap1 signaling pathway, Ras signaling pathway, proteoglycans in cancer, and central carbon metabolism in cancer. These pathways are well-documented in the context of various cancers. For instance, the PI3K/Akt signaling pathway is implicated in the progression and metastasis of gallbladder cancer (28), while the Rap1 signaling pathway, interconnected with AKT signaling, is known to promote esophageal squamous cell carcinoma metastasis (29). VEGFC’s role in HNSCC may, therefore, be mediated through these pathways.

In recent years, the role of immune cell infiltration in tumor microenvironment has gained considerable attention for its involvement in tumor occurrence, development, and metastasis. It has emerged as a valuable diagnostic and prognostic biomarker across various tumors such as gynecological tumors and osteosarcoma (30,31). A study on HNSCC identified a significant increase in p16 expression in oropharyngeal cancer, which was linked to improved OS and enhanced infiltration of T and B lymphocytes, suggesting a potential correlation between immune cell infiltration and patient prognosis (32). Our study echoes these findings, suggesting that immune cell infiltration may play a role in the pathogenesis and progression of HNSCC. We observed that VEGFC is associated with the infiltration of five major immune cells types in HNSCC: B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, and dendritic cells. To elucidate the relationship between VEGFC and immune cell infiltration, we assessed the impact of VEGFC expression on immune cell levels in HNSCC and investigated its association with various immune feature molecules. Among the 28 TILs-related immune feature molecules, we found that Tcm CD4 and Tgd, immunoinhibitor-related molecule such as PDCD1LG2, TGFB1, and TGFBR1, immunostimulator-related molecules including NT5E, CD276, and PVR, MHC class-related molecules like TAP2, HAL-C, HAL-B, HAL-A, HAL-E, HAL-F, and B2M, as well as chemokine and receptor-related molecules such as CXCL12 and CXCR1, all demonstrated a high correlation with VEGFC expression. This suggests that VEGFC is extensively involved in the regulation of immune molecules in HNSCC, thereby affecting immune infiltration within the tumor microenvironment. These findings position VEGFC as a promising candidate for molecular targeting in immunotherapy for HNSCC.

In summary, this study found that VEGFC is highly expressed in HNSCC and this high expression correlates with clinical staging, patient gender, age, race, etc. Patients with high VEGFC expression exhibit significantly lower OS and DFS rates compared to those with low expression. Univariate and multivariate Cox regression analyses further indicate VEGFC as an independent risk factor for HNSCC. The impact of VEGFC on HNSCC pathogenesis appears to be multifaceted, potentially through the interplay with genes like ITGA3, NT5E, PXN, LAMA3, CAV1, TNFRSF12A, LIMA1, and LAMC2. These genes are implicated in, tumor progression via pathways pivotal to cancer development, including the PI3K/Akt signaling pathway, MAPK signaling pathway, and focal adhesion. Furthermore, VEGFC’s association with the infiltration of five major immune cells types—B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, and dendritic cells—along with its correlation with spectrum of immune regulatory molecules such as TILs, immunostimulators, and chemokines and receptors. This modulation may have profound implications for immune cell infiltration and the overall immune response within the HNSCC tumor context.

Conclusions

Our study consolidates the findings that VEGFC is significantly upregulated in HNSCC and stands out as an independent adverse prognostic factor. The correlation of VEGFC with clinical-pathological variables and its predictive power for OS and DFS underscore its potential utility in patient risk stratification. Furthermore, VEGFC’s association with immune cell infiltration and its role in pivotal pathways suggest it as a promising target for immunotherapeutic approaches in HNSCC. These insights pave the way for the development of personalized treatment strategies, offering a more nuanced understanding of molecular factors affecting patient outcomes.

Acknowledgments

Funding: This research was supported by the research project of Youjiang Medical University for Nationalities in 2018 (grant No. yy2018ky015). The funder had no role in study design, data collection and analysis.

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-834/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-834/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-834/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Ju G, Yao Z, Zhao Y, et al. Data mining on identifying diagnosis and prognosis biomarkers in head and neck squamous carcinoma. Sci Rep 2023;13:10020. [Crossref] [PubMed]
2. Ershadifar S, Young K, Birkeland AC. Utility of miRNA biomarkers for detection of early head and neck squamous cell carcinoma. Transl Cancer Res 2023;12:673-5. [Crossref] [PubMed]
3. Peng Q, Jiang X, Tan S, et al. Clinical significance and integrative analysis of the cuproptosis-associated genes in head and neck squamous cell carcinoma. Aging (Albany NY) 2023;15:1964-76. [Crossref] [PubMed]
4. Bai G, Yue S, You Y, et al. An integrated bioinformatics analysis of the S100 in head and neck squamous cell carcinoma. Transl Cancer Res 2023;12:717-31. [Crossref] [PubMed]
5. Naumov SS, Kulbakin DE, Krakhmal NV, et al. Molecular and biological factors in the prognosis of head and neck squamous cell cancer. Mol Biol Rep 2023;50:7839-49. [Crossref] [PubMed]
6. Umemori K, Ono K, Eguchi T, et al. EpEX, the soluble extracellular domain of EpCAM, resists cetuximab treatment of EGFR-high head and neck squamous cell carcinoma. Oral Oncol 2023;142:106433. [Crossref] [PubMed]
7. Zhang Y, Chen K, Li L, et al. CCR4 is a prognostic biomarker and correlated with immune infiltrates in head and neck squamous cell carcinoma. Ann Transl Med 2021;9:1443. [Crossref] [PubMed]
8. Yang FL, Wei YX, Liao BY, et al. LncRNA HOTAIR regulates the expression of E-cadherin to affect nasopharyngeal carcinoma progression by recruiting histone methylase EZH2 to mediate H3K27 trimethylation. Genomics 2021;113:2276-89. [Crossref] [PubMed]
9. Yang Y, Wang X, Wang P. Signaling mechanisms underlying lymphatic vessel dysfunction in skin aging and possible anti-aging strategies. Biogerontology 2023;24:727-40. [Crossref] [PubMed]
10. Dumond A, Montemagno C, Vial V, et al. Anti-Vascular Endothelial Growth Factor C Antibodies Efficiently Inhibit the Growth of Experimental Clear Cell Renal Cell Carcinomas. Cells 2021;10:1222. [Crossref] [PubMed]
11. Qin T, Xia J, Liu S, et al. Clinical importance of VEGFC and PD-L1 co-expression in lung adenocarcinoma patients. Thorac Cancer 2020;11:1139-48. [Crossref] [PubMed]
12. Li J, Han T. Comprehensive analysis of the oncogenic roles of vascular endothelial growth factors and their receptors in stomach adenocarcinoma. Heliyon 2023;9:e17687. [Crossref] [PubMed]
13. Penco-Campillo M, Comoglio Y, Feliz Morel ÁJ, et al. VEGFC negatively regulates the growth and aggressiveness of medulloblastoma cells. Commun Biol 2020;3:579. [Crossref] [PubMed]
14. Chang TM, Chu PY, Hung WC, et al. c-Myc promotes lymphatic metastasis of pancreatic neuroendocrine tumor through VEGFC upregulation. Cancer Sci 2021;112:243-53. [Crossref] [PubMed]
15. Benke EM, Ji Y, Patel V, et al. VEGF-C contributes to head and neck squamous cell carcinoma growth and motility. Oral Oncol 2010;46:e19-24. [Crossref] [PubMed]
16. Pentheroudakis G, Angouridakis N, Wirtz R, et al. Transcriptional activity of human epidermal growth factor receptor family and angiogenesis effectors in locoregionally recurrent head and neck squamous cell carcinoma and correlation with patient outcome. J Oncol 2009;2009:854127. [Crossref] [PubMed]
17. Wang M, Wisniewski CA, Xiong C, et al. Therapeutic blocking of VEGF binding to neuropilin-2 diminishes PD-L1 expression to activate antitumor immunity in prostate cancer. Sci Transl Med 2023;15:eade5855. [Crossref] [PubMed]
18. Wang J, Huang Y, Zhang J, et al. Pathway-related molecules of VEGFC/D-VEGFR3/NRP2 axis in tumor lymphangiogenesis and lymphatic metastasis. Clin Chim Acta 2016;461:165-71. [Crossref] [PubMed]
19. Kong D, Zhou H, Neelakantan D, et al. VEGF-C mediates tumor growth and metastasis through promoting EMT-epithelial breast cancer cell crosstalk. Oncogene 2021;40:964-79. [Crossref] [PubMed]
20. Vimalraj S, Hariprabu KNG, Rahaman M, et al. Vascular endothelial growth factor-C and its receptor-3 signaling in tumorigenesis. 3 Biotech 2023;13:326.
21. Zhang X, Wang JF, Chandran B, et al. Kaposi's sarcoma-associated herpesvirus activation of vascular endothelial growth factor receptor 3 alters endothelial function and enhances infection. J Biol Chem 2005;280:26216-24. [Crossref] [PubMed]
22. Zheng X, Du Y, Liu M, et al. ITGA3 acts as a purity-independent biomarker of both immunotherapy and chemotherapy resistance in pancreatic cancer: bioinformatics and experimental analysis. Funct Integr Genomics 2023;23:196. [Crossref] [PubMed]
23. Yao J, Wang L. Integrin α3 Mediates Stemness and Invasion of Glioblastoma by Regulating POU3F2. Curr Protein Pept Sci 2023;24:247-56. [Crossref] [PubMed]
24. Du J, Zhao Y, Hu D, et al. Silencing of integrin subunit α3 inhibits the proliferation, invasion, migration and autophagy of esophageal squamous cell carcinoma cells. Oncol Lett 2022;24:271. [Crossref] [PubMed]
25. Steglich A, Vehlow A, Eke I, et al. α integrin targeting for radiosensitization of three-dimensionally grown human head and neck squamous cell carcinoma cells. Cancer Lett 2015;357:542-8. [Crossref] [PubMed]
26. Shi E, Wu Z, Karaoglan BS, et al. 5'-Ectonucleotidase CD73/NT5E supports EGFR-mediated invasion of HPV-negative head and neck carcinoma cells. J Biomed Sci 2023;30:72. [Crossref] [PubMed]
27. Wang X, Zhang C, Song H, et al. Characterization of LIMA1 and its emerging roles and potential therapeutic prospects in cancers. Front Oncol 2023;13:1115943. [Crossref] [PubMed]
28. Wu Z, Yu X, Zhang S, et al. The role of PI3K/AKT signaling pathway in gallbladder carcinoma. Am J Transl Res 2022;14:4426-42. [PubMed]
29. Li Q, Xu A, Chu Y, et al. Rap1A promotes esophageal squamous cell carcinoma metastasis through the AKT signaling pathway. Oncol Rep 2019;42:1815-24. [Crossref] [PubMed]
30. Liu QF, Feng ZY, Jiang LL, et al. Immune Cell Infiltration as Signatures for the Diagnosis and Prognosis of Malignant Gynecological Tumors. Front Cell Dev Biol 2021;9:702451. [Crossref] [PubMed]
31. Yang H, Zhao L, Zhang Y, et al. A comprehensive analysis of immune infiltration in the tumor microenvironment of osteosarcoma. Cancer Med 2021;10:5696-711. [Crossref] [PubMed]
32. Schneider K, Marbaix E, Bouzin C, et al. Immune cell infiltration in head and neck squamous cell carcinoma and patient outcome: a retrospective study. Acta Oncol 2018;57:1165-72. [Crossref] [PubMed]