Data from laryngeal cancer patients was obtained from The Cancer Genome Atlas (TCGA) and Clinical Proteomic Tumor Analysis Consortium (CPTAC) (Huang et al., 2021). The processed data were accessed via LinkedOmics: http://www.linkedomics.org, including transcriptome, proteome, and corresponding clinical features of all patients from these two datasets. Samples lacking significant clinicopathological or survival information were excluded from further analysis.

Single-variable Cox proportional hazards regression analysis was conducted to identify features that were significantly associated with overall survival (OS) or progression-free survival (PFS) in LSCC cohorts. The results were visualized using a forest plot created with the R package “ggplot2”. The Kaplan-Meier plotter utilizing the R package “survival” was employed to compare the OS or PFS times between specific groups.

The somatic mutation files were obtained from the TCGA and CPTAC databases. We utilized the “Maftool” R package to display a waterfall plot visualizing the top 20 genes with the highest tumor mutation frequency (TMF) in LSCC patients from TCGA and CPTAC (Mayakonda et al., 2018). Different colors were used to annotate specific mutation types. To evaluate the correlation between lymphatic metastasis and mutations, a chi-square test was conducted, and a P-value <0.05 was considered statistically significant.

The “Limma” package of R software was applied to evaluate the differential expression proteins (DEPs) and differential expression genes (DEGs) between lymphatic metastasis and non-lymphatic metastasis LSCC (Ritchie et al., 2015). The DEPs were defined as proteins with P-value <0.05 and log2 (Fold Change) > 0.585. The DEGs were defined as genes with P-value <0.05 and log2 (Fold Change) > 1. Functional enrichment analysis, including gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG), was performed using the DAVID bioinformatics resources (Sherman et al., 2022). The Benjamin–Hochberg adjusted P < 0.05 was regarded as statistically significant.

The DEPs and DEG encoded proteins were annotated based on the cellular component of the Gene Ontology. The frequency of proteins found in different locations, such as the cytosol, membrane, nucleus, cytoskeleton, mitochondria, Golgi apparatus, and lysosome, is illustrated in the bubble diagram. For epitope prediction, the PDB file of the protein intended for analysis was downloaded from the AlphaFold Protein Structure Database. Afterwards, the files were used as input for Ellipro (http://tools.iedb.org/ellipro/) to estimate the frequency of discontinuous epitopes for each structure (Ponomarenko et al., 2008).

The plasma MPs that were overexpressed in lymphatic metastasis LSCC were analyzed using the algorithm developed by Anderson et al. for prioritizing immunotherapy targets (Anderson et al., 2022). However, some slight modifications were made to adapt the algorithm to our specific study. The following were determined: fold-change of expression in tumor compared to normal; significance of the difference in expression between the lymphatic metastasis group and the non-lymphatic metastasis group; and off-tumor expression according to Human Proteome Map (Kim et al., 2014). To ensure that each of the three features had a range of 0–1, they were individually scaled based on their maximum value. For each protein, the scaled values were converted into a three-dimensional vector in R3. The final score was then calculated based on the magnitude of this vector, which was capped at a maximum value of approximately 1.732.

For gene set enrichment analysis (GSEA), the patients were divided into high- and low-expression groups according to the median protein abundance. GSEA was performed to identify the primarily enriched pathways using the “clusterProfiler” package of R software (Wu et al., 2021). Pathways with the nominal P < 0.01 were considered statistically significant. Protein-protein interaction (PPI) networks were constructed using the STRING database (v.11.5) (Szklarczyk et al., 2019). The minimum required interaction score was set to a threshold of 0.400.

LSCC cell lines (TU212 and TU686) were purchased from the Cell Center of Life Science of Chinese Academy of Science (Shanghai, China). TU212 and TU686 cells were maintained in RPMI-1640 medium (Gibco, Carlsbad, CA), supplemented with 10% fetal bovine serum (FBS, Gibco, United States) and 1% penicillin/streptomycin. All cell lines were incubated and cultured at 37°C and 5% CO₂. EpCAM siRNA (siEpCAM), MGST1 siRNA (siMGST1) and their respective negative controls (siNC) were purchased from Sangon (Shanghai, China). For transfection, cells were seeded onto plates and transfected with specific siRNA using Lipofectamine 3,000 reagent (Invitrogen, Carlsbad, CA, United States) following the manufacturer’s instructions.

LSCC cells were seeded in 150 µL of medium in 96-well plates and incubated at 37°C and 5% CO₂ for 24 h. Next, cells were incubated in 3-[4, 5-dime-thylthiazol-2-yl]-2, 5-diphenyl tetrazolium bromide (MTT) solution (5 mg/mL) for 4 h in a humidified incubator. After dimethyl sulfoxide (DMSO) was added to each well, absorbance was measured with a microplate reader at 570 nm. The experiment was performed at least three times (Bio-Rad, Hercules, CA).

Tu212 and Tu686 cells were collected and washed with PBS, lysed in RIPA buffer containing PMSF (phenylmethylsulfonyl fluoride) and protease and phosphatase inhibitors on ice for 30 min. Proteins were separated by SDS-PAGE and transferred to a polyvinylidene fluoride (PVDF) membrane. Next, the membrane was blocked with a 5% nonfat milk, then incubated with primary antibodies at 4°C overnight. The primary antibodies were used: GAPDH (AC001), EpCAM (A19301), and MGST1 (A0880) from ABclonal. After washing with PBS + Tween 20, cells were incubated with the appropriate secondary antibodies for 2 h at room temperature. Protein bands were detected with ECL (enhanced chemiluminescence) reagents and the ChemiDoc MP imaging system (Bio-Rad). ImageJ software was utilized to perform protein quantification.

The migration and invasion of TU212 and TU686 cells were evaluated by transwell assay. For the migration assay, cells were cultured in serum-free RPMI 1640 media in the upper chamber (8-μm pore size, Corning Inc.; New York, United States) and 600 µL medium containing 20% FBS was placed in the lower chamber. After incubation at 37°C and 5% CO₂ for 24 h, cells were fixed with a 4% PFA solution and stained with 0.5% crystal violet for 20 min. For the invasion assay, the protocol was similar to the migration assay but for cells were grown in Matrigel-coated chambers. The number of migrated and invasive cells was photographed and counted across 5 random fields under an inverted microscope.

Continuous variables were compared by Student’s t-test, while the categorical variables were compared using the chi-square test. The Kaplan-Meier method was used to generate OS and PFS curves, and the log-rank test was performed. Statistical analyses were performed using R version 4.3.2. A P-value <0.05 (two sides) was considered statistically significant.

Lymph node metastasis is widely recognized as a crucial factor contributing to a poor prognosis among tumor patients. The baseline characteristics of LSCC patients from the TCGA database and CPTAC are shown in Supplementary Table S1. The univariate independent prognostic analysis was conducted on CPTAC-LSCC and TCGA-LSCC cohorts, which showcased a significant correlation between lymph node metastasis and prognosis in LSCC (Figures 1A,D). The patients were then divided into two groups based on the presence or absence of lymph node metastasis. The group without lymph node metastasis demonstrated a longer PFS in comparison to the group with lymph node metastasis, with P-values of 0.01898 for CPTAC cohort and 0.03245 for TCGA cohort. (Figures 1B,E). Consistently, the OS of the group without lymph node metastasis significantly surpassed that of the group with lymph node metastasis, with P-values of 0.033 for CPTAC cohort and 0.030 for TCGA cohort (Figures 1C,F). In line with the current consensus, both TCGA-LSCC and CPTAC-LSCC cohorts displayed the poor prognosis of patients with lymph node metastasis. This suggests that these cohorts are appropriate for investigating the mechanisms of lymph node metastasis and identifying novel therapeutic targets for LSCC. In light of this, gene expression data from 111 LSCC cases was obtained from the TCGA database, while protein expression data from 47 LSCC cases was collected from the CPTAC database. As shown in Figure 1G, a proteogenomic analysis workflow was conducted to examine gene mutations and protein expressions associated with lymph node metastasis. The cell MPs associated with lymph node metastasis were identified and assessed for their potential as therapeutic targets. Additionally, in vitro functional experiments were carried out to investigated the underlying mechanisms through which these proteins may regulate lymph node metastasis.

It has been reported that certain gene mutations are correlated with the risk of lymph node metastasis in some types of tumors. To investigate the mutation events related to lymph node metastasis in laryngeal cancer, we acquired whole-exome sequencing data of the TCGA-LSCC and CPTAC-LSCC cohorts. Then the correlation between the top 20 genes with the highest mutation frequencies and lymph node metastasis were analyzed. The results revealed that certain genes had significantly higher mutation frequencies in the non-lymphatic metastasis group compared to the lymphatic metastasis group. These genes include LRP1B (P = 0.000360), NSD1 (P = 0.033), USH2A (P = 0.007), FAT1 (P = 0.026), PAPPA2 (P = 0.046), and RYR2 (P = 0.048) (Supplementary Figures 1A,B). Among these, the LRP1B and USH2A, which are common mutant genes across various cancer types (Brown et al., 2021; López et al., 2023; Möhrmann et al., 2022; Yang et al., 2023), were found to be highly expressed in the lymphatic metastasis group (Supplementary Figure 1C) (Han et al., 2021). The low expression or mutation of LRP1B has been reported to be potentially responsible for a lower risk of lymphatic metastasis, which aligns with our findings in this study (Han et al., 2021).

To gain a comprehensive understanding of the regulatory factors associated with lymph node metastasis in LSCC, we performed an integrative analysis on the proteome and corresponding transcriptome of the CPTAC-LSCC cohort. By conducting differential expression analysis between the non-lymphatic metastasis group and the lymphatic metastasis group, we discovered 228 differentially expressed proteins (DEPs) including 54 upregulated proteins and 174 downregulated proteins (with a fold change >1.5 and p < 0.05) (Figure 2A). Additionally, 535 differentially expressed genes (DEGs), consisting of 277 upregulated genes and 258 downregulated genes (with a mRNA fold change >2 and p < 0.05) were identified (Figure 2C).

Functional enrichment analysis was conducted using the DAVID service to generate representative biological processes involved in lymph node metastasis. As depicted in Figures 2B,D, the DEPs and DEGs were primarily enriched among terms such as “muscle contraction,” “intermediate filament organization,” “sarcomere organization,” “keratinization,” and “epidermis development.” These processes are closely associated with cell growth and motility. The DEPs and DEGs were then classified based on their subcellular localizations, using the Cell component (CC) catalog of the Gene Ontology (GO) database. Majority of DEPs were localized in the cytosol, cell membrane, and nucleus (Figure 2E). Given the accessibility of MPs in the targeted therapy or immunotherapy, our interest has been centered around the highly expressed MPs found in the lymph node metastasis group. These proteins have the potential of serving as therapeutic targets for lymph node metastatic LSCC.

We have previously discovered that some of the DEPs related to lymph node metastasis are located on the cell membrane (Figure 2E). Differential analysis was conducted on all quantified cell MPs within the proteome data, resulting in the identification of 15 cell MPs that were overexpressed in the tumors with lymph node metastasis and 39 downregulated MPs (Figure 3A). The segment of the MP that is exposed on the cell surface is a highly promising drug target, as approximately 60% of drug targets are MPs (Overington et al., 2006). Figure 3B displays the expression level of 15 metastasis-specific MPs, the majority of which exhibited overexpression in tumor tissues in comparison to normal tissues (Figure 3C). To verify their potency as therapeutic targets, we first performed antibody epitope prediction using three-dimensional protein structure models generated by the Alphafold database, utilizing the ElliPro server. The results demonstrated that these proteins possess antigenic epitopes of at least 70 amino acids or more, making them susceptible to targeting by monoclonal antibodies (Figure 3D). Subsequently, a prioritization of these 15 MPs as immunotherapy targets was conducted, taking into account their distribution across various organs, abundance in laryngeal carcinoma, and correlation with lymph node metastasis (Figures 3E,F). Epithelial cell adhesion molecule (EpCAM) and Microsomal glutathione S-transferase 1 (MGST1) are the top-scoring proteins, exhibiting significantly higher scores than others (Figure 3F).

Since EpCAM and MGST1 were found to be overexpressed in lymphatic metastatic LSCC, we investigated the prognostic relevance of these proteins. Our findings showed that patients with overexpression of EpCAM or MGST1 had a worse prognosis (Figures 4A,B). EpCAM and MGST1 showed significant overexpression at both the protein and mRNA levels in lymphatic metastatic LSCC (Figure 4C). It is also noteworthy that expression of EpCAM was higher in cancer tissue compared to adjacent tissue. On the other hand, MGST1 was only upregulated in lymphatic metastatic LSCC, with no significant difference between cancer and adjacent tissue (Figure 4D). Furthermore, the risk curve illustrates a positive correlation between EpCAM and MGST1 with disease progression, which aligns with the findings of the survival curve (Figure 4E).

To explore the role of EpCAM and MGST1 in LSCC cells, TU212 and TU686 cells were separately transfected with siRNA-targeting EpCAM, siRNA-targeting MGST1, or non-specific RNAi control (Figure 5A). Western blot showed more than 50% knockdown efficiency in TU212 and TU686 cells (Figure 5B) for both EpCAM and MGST1. MTT assays showed that silencing EpCAM or MGST1 significantly reduced the proliferation of TU686 and TU212 cells (Figures 5C,D). The transwell assay demonstrated that the reduction of EpCAM leads to a significant decrease in the migration and invasion of LSCC cells (Figure 5E), and similar effects were observed with MGST1 knock-down (Figure 5F). These results collectively provided evidence that the deficiency of EpCAM or MGST1 suppresses the proliferation, migration, and invasion of LSCC cells. These proteins might therefore represent promising therapeutic opportunities in attempts to improve the poor prognoses of lymphatic metastatic LSCC.

EpCAM is regarded as a carcinoma cell-surface marker involved in cell adhesion, proliferation, migration, stemness, and epithelial-to-mesenchymal transition. MGST1 is a membrane-bound glutathione transferase with the ability to detoxify reactive intermediates. Our findings demonstrated that overexpression of EpCAM and MGST1 are closely correlated with lymphatic metastasis and play crucial role in driving cell growth, migration, and invasion. However, there was no significant correlation between the expression of these two proteins (Figure 6A), implying distinct mechanisms by which EpCAM and MGST1 are likely involved in lymphatic metastasis. Here we aimed to examine the regulatory mechanisms behind the disordered expression of these proteins in LSCC.

Patients with high EpCAM featured more CDKN2A mutations and tumors of higher histologic grade (p = 0.002) (Fisher’s exact test) (Figures 6B,C). Tumors with high EpCAM showed specific elevation in pathways involving spliceosome and mismatch repair, and downregulation in pathways related to endocytosis and regulation of the cytoskeleton (Figures 6C,D). It can be inferred from Figure 4 that the overexpression of EpCAM, along with its facilitation of lymph node metastasis, could be associated with a protein-level regulatory mechanism, such as PPI. Given this, PPI analyses were carried out by STRING server (Figure 6E), the EpCAM-interacting proteins were mainly enriched in cancer-related biological processes, including the Wnt signaling pathway, adhesion junctions, and the Hippo signaling pathway (Figure 6F). Patients with low MGST1 were characterized by harboring more RYR2, LRP1B, or DNAH5 mutations (Figure 6B). Tumors with high MGST1 showed specific upregulation of pathways relevant to ribosome assembly and oxidative phosphorylation and downregulation of pathways involving extracellular matrix, proteoglycans of binding, protein digestion and absorption (Supplementary Figures S2A,B). The top three proteins showing the highest correlation with MGST1 were SUPT16H, NOP2, and APEX1. Both APEX1 and NOP2 are recognized regulators of tumor metastasis, suggesting a potential mechanism wherein the lymphatic metastasis associated with high MGST1 could be attributed to the overexpression of either APEX1 or NOP2 (Supplementary Figures S2C–E). Taken together, these in silico analyses may provide guidance for further research on the regulatory mechanisms of EpCAM and MGST1 in cancers.

It has been reported that EpCAM and MGST1 participate in cell signaling, including phosphaditylinositol-3 kinase and wnt/β-catenin pathway (Dai and Lu, 2022; Li, et al., 2023c; Munz et al., 2009; Ni et al., 2013). This is significant due to the support of these pathways in tumor cell proliferation, survival, and anti-apoptotic responses. To investigate the mechanism by which EpCAM or MGST1 regulates the progression of LSCC, we examined the effects of EpCAM and MGST1 knockdown on the wnt/β-catenin and PI3K signaling pathways in TU686 and TU212 cells. As anticipated, silencing EpCAM resulted in downregulation of PIK3CA, NFKB, and phosphorylated NFKB, as well as the inhibition of β-catenin and downstream phosphorylated GSK3β (Figures 7A,B). Of note, knockdown of MGST1 resulted in a similar signaling interference (Figures 7C,D). These findings support the mechanism whereby the two biomarkers, enriched in lymphatic metastatic LSCC, are consistently associated with the Wnt/β-catenin or PI3K signaling pathways, thus leading to tumor-promoting effects.

While the integrated proteomic approach identified EpCAM and MGST1 as promising therapeutic targets in metastatic LSCC, this study has limitations that merit consideration: Functional characterization of EpCAM and MGST1 was confined to in vitro systems, and in vivo validation using appropriate animal models would be valuable to substantiate translational potential. Although integration of TCGA and CPTAC datasets provided a cohort, larger independent cohorts incorporating additional proteomics datasets would strengthen generalizability. Furthermore, while the data suggest modulation of the Wnt/β-catenin-PI3K pathway, elucidation of the precise molecular mechanisms governing EpCAM/MGST1 interactions with these pathways represents an important direction for future investigation.