Metabolic expression profiling analysis reveals pyruvate-mediated EPHB2 upregulation promotes lymphatic metastasis in head and neck squamous cell carcinomas

Abstract

Lymphatic metastasis is a well-known factor for initiating distant metastasis of head and neck squamous cell carcinoma (HNSCC), which caused major death in most patients with cancer. Meanwhile, metabolic reprogramming to support metastasis is regarded as a prominent hallmark of cancers. However, how metabolic disorders drive in HNSCC remains unclear. We firstly established a new classification of HNSCC patients based on metabolism gene expression profiles from the TCGA and GEO database, and identified an enriched carbohydrate metabolism subgroup which was significantly associated with lymphatic metastasis and worse clinical outcome. Moreover, we found that highly activated pyruvate metabolism endowed tumors with EPHB2 upregulation and promoted tumor lymphangiogenesis independently of VEGF-C/VEGFR3 signaling pathway. Mechanically, high nuclear acetyl-CoA production from pyruvate metabolism promoted histone acetylation, which in turn transcriptionally upregulated EPHB2 expression and secretion in tumor cells. EPHB2 bound with EFNB1 in lymphatic endothelial cells promoted YAP/TAZ cytoplasmic retention, which alleviated YAP/TAZ-mediated prospero homeobox protein 1 (PROX1) transcriptional repression, and then triggered tumor lymphangiogenesis. Importantly, combined treatment with EFNB1-Fc and VEGFR3 inhibitor synergistic abrogated lymphangiogenesis in vitro and in vivo, suggesting that targeting EPHB2 might be a potential strategy to patients with no or slight response to VEGFR3 inhibitor. These findings uncover the mechanism by which pyruvate metabolism is linked to lymphatic metastasis of tumor and provides a promising therapeutic strategy for the prevention of HNSCC metastasis.

Introduction

Approximately 20–30% of patients with head and neck squamous cell carcinoma (HNSCC) suffer from local and regional lymphatic metastasis, which can lead to distant metastasis and cause poor clinical outcomes [1,2,3]. Lymphangiogenesis is considered to be the predominant mechanism for the expansion of tumor-associated lymphatic network, which mediated by lymphatic endothelial cells (LECs) with expressing characteristic molecular markers, including lymphatic vessel endothelial receptor 1 (LYVE-1), prospero homeobox protein 1 (PROX1) and podoplanin (PDPN) [4]. Given that approximately 95% of tumor-associated vessels infiltrated by tumor cells are lymphatic vessels, which play a major role in tumor metastasis, targeting lymphangiogenesis has been considered as a potent and direct strategy for controlling metastatic diseases.

VEGF-C, one of the vascular endothelial growth factors (VEGF) family, has long been recognized as a major molecular driver of lymphangiogenesis by activation of VEGFR3 expressed on LECs [5, 6]. The VEGF-C/VEGFR3 interaction promotes proliferation, migration and tube formation of LECs to induce tumor lymphangiogenesis via maintenance of the MAPK or PI3K/Akt signaling activity [5, 6]. Based on its essential role in lymphangiogenesis, numerous drugs targeting VEGF-C/VEGFR3 pathway have been established and applied to clinical practice or under the clinical testing and presented partial efficacy against lymphatic metastasis in various cancers [7], especially in renal cell carcinoma [8] and colorectal cancer [9]. However, numerous clinical trials demonstrated that the magnitude of the clinical benefit of targeting VEGF-C/VEGFR3 pathway varies among different patient populations [10,11,12]. Moreover, some patients with lymphatic metastatic cancers exhibit low expression of VEGF-C, suggesting the occurrence of VEGF-C-independent lymphangiogenesis. Therefore, further research is needed to identify the comprehensive mechanism of lymphangiogenesis in HNSCC.

Metabolic reprogramming is one of the well-established hallmarks of cancer and is responsible for the metastatic spread of tumors [13]. In order to adapt the changes in energy requirements at each step during the metastatic cascade, tumor cells selectively and dynamically modify their metabolic programs [14, 15]. Recently, more and more studies showed that tumor cells have heterogeneous metabolic preferences and dependencies. Within the same tumor, tumor cells with higher carbohydrate metabolism, which provided fuel of energy production, presented a more aggressive and malignant phenotype via increasing oxidative phosphorylation [16, 17]. Compared with non-metastatic cancer, the serum pyruvate concentration was higher in patients with metastatic breast cancer [18] and highly invasive ovarian cancer [19]. Pyruvate uptake inducing the production of α-ketoglutarate activated collagen hydroxylation by increasing the activity of the enzyme collagen prolyl-4-hydroxylase (P4HA), which remodeling the extracellular matrix in the lung metastatic niche for breast cancer [20]. Multiple studies have provided evidence that metabolic changes in cancer cells are not just consequential bystander effects because of metabolites, but also could directly modulate activity of signaling pathways and global gene expression programmes [21,22,23]. Hence, deciphering how metabolic disorders drive lymphatic metastasis may provide potential therapeutic vulnerabilities in HNSCC.

Here, we firstly classified HNSCC patients into four subtypes based on unsupervised clustering of metabolism gene expression profiles using TCGA public data and evaluated in GEO cohort. We found that an enriched carbohydrate metabolism subtype, Cluster D, was significantly associated with lymphatic metastasis and worse clinical outcome. Meanwhile, we also investigated the underlying mechanism of carbohydrate metabolism triggering lymphatic metastasis in laryngeal cancer. Importantly, blocking EPHB2 with recombinant EFNB1-Fc synergised with VEGFR3 inhibitor treatment to potently suppress lymphatic metastasis in laryngeal cancer. Our findings bring new insight into the links between pyruvate metabolism and lymphatic metastasis, and explore strategies to optimize the treatment for lymphatic metastasis in HNSCC.

Methods and materials

Patients and datasets

This study collected 725 HNSCC patients from two databases: TCGA (n = 532) and GEO (n = 193). For the TCGA cohort, the RNA-seq data and corresponding clinical information were downloaded from the TCGA database (http://cancergenome.nih.gov/). The molecular and immune features were also retrieved [24]. For the GEO cohort, the RNA-seq accompanied with clinical information were obtained from the GEO database (GSE142083, GSE127165, GSE112026 and GSE74927).

All human samples used in this study were collected from the Guangdong Provincial People’s Hospital. Twenty-eight cases of fresh tissues, including 8 of nasopharyngeal carcinoma tissues, 6 of oropharyngeal cancer tissues, 6 of hypopharyngeal cancer tissues and 8 of laryngeal cancer tissues with or without lymphatic metastasis were used for western blot analysis. Prior patients’ consents and approval from the Institutional Research Ethics Committee (2020180 H) were obtained for the use of these clinical specimens in research.

Identification and validation of metabolic subtypes

The 1064 metabolism-related genes from MSigDB were selected for consensus clustering (R package “consensusClusterplus”) to identify robust clusters in the training cohort from TCGA and the genes were detailed in Supplementary Table. S1. We utilized two methods of elbow method and gap statics to identify the optimized K categories, meanwhile considering the associations between metabolic clusters and survival outcomes. The cumulative distribution function (CDF) and consensus heatmap were adopted to assess the optimal K.

To validate the metabolic subtypes in GEO cohort, we first used the R package “limma” to remove batch effects on the 4 GEO-derived HNSCC datasets. Subsequently, we trained a classifier based on the subtypes in TCGA cohort via “Random Forest” algorithm with R package “randomForest”. Then, the classifier was used to predict metabolic subtype for patients in the GEO cohort and each sample was assigned to a metabolic subtype. The in-group-proportion (IGP) statistic with R package “clusterRepro” was performed to assess the reproducibility and similarity of the acquired metabolic subtypes between the training and validation cohorts.

Gene set variation analysis (GSVA)

GSVA was utilized to evaluate the most significantly enriched molecular pathways of the stemness subtypes by using the “GSVA” package in R [25]. Differential analysis of the enrichment scores of molecular pathways between the different groups was performed by the limma package in R [26]. the signatures with an absolute log2 fold change (FC) > 0.2 (adjusted P < 0.05) were defined as differentially expressed signatures.

Connective mapping analysis

The CMap database (https://clue.io/) was used to explore potential compounds targeting the molecular pathways and key genes associated with the metabolic subtype of HNSCC patients [27]. The DEGs between the different metabolic subtypes were employed to query the CMap database, and the most significantly highly expressed genes of each metabolic subtype were considered as key regulatory gene.

Cell culture

The State Key Laboratory of Oncology in South China (Sun Yat-sen University Cancer Center, Guangzhou, China) provided the cells lines. The normal nasopharyngeal epithelial cell line NP69 and laryngeal epithelial cell line HuLa-PC were cultivated in Keratinocyte Growth Medium 2 (C-39011, Sigma-Aldrich). The nasopharyngeal carcinoma cell lines HK1 and SUNE1, oral squamous carcinoma cell line SCC9 were grown in RPMI 1640 (Invitrogen, Carlsbad, CA, USA) medium supplemented with 10% fetal bovine serum (FBS, Life Technologies, Carlsbad, CA). The oral squamous carcinoma cell line UM-SCC-1, hypopharyngeal squamous carcinoma cell line FaDu, and laryngeal squamous cell carcinoma cell line Tu212, Tu686 and AMC-HN-8 were cultivated in DMEM (Invitrogen, Carlsbad, CA, USA) medium supplemented with 10% FBS. Human lymphatic endothelial cells (LECs) were cultured in endothelial cell medium (Sciencell), containing 5% FBS, 1% endothelial cell growth supplement and 1% penicillin/streptomycin solution. All cells maintained in an incubator at 37 ℃ with 5% CO₂ and 100% humidity with the medium changed every other day.

Plasmids, infection, and transfection

The human EPHB2, YAP (Full, 5SA), EFNB1 (Full, Y317A/Y331A), TAZ and shRNAs targeting EPHB2 or VEGFR3 were cloned into the pLVX-TRE3G-IRES vector. Transfection of siRNAs or plasmids was performed using the Lipofectamine 3000 reagent (Thermo Fisher Scientific, Waltham) according to the manufacturer’s instruction. All primers and oligonucleotides are listed in Supplementary Table 2. Stable cell lines were generated via lentiviral infection using 293T cells and selected for 10 days with 0.5 µg/mL puromycin 48 h after infection.

RNA isolation and qRT-PCR

RNA from the indicated cells were extracted using the Trizol reagent (15596026, Invitrogen). cDNA synthesis was carried out using M-MLV Reverse Transcriptase and then the cDNA was used as the template for amplification on a CFX384 Real-Time System (Bio-Rad Laboratories, Singapore). The data were analyzed using the ΔΔCt method. The results from each experiment were normalized to the expression of ACTB. The sequences of the primers used for all qRT-PCR assays are listed in Supplementary Table. S1.

Western blot analysis

Western blot analyses were performed according to a standard protocol using primary antibodies, including anti-EPHB2 (Cell Signaling Technology, 14389), anti-H3 (Cell Signaling Technology, 9715), anti-H3ac (pan-acetyl) (Invitrogen, PA5-114693), anti-H4 (Cell Signaling Technology, 2592), anti-H4ac (pan-acetyl) (Invitrogen, PA5-32029), anti-EFNB1 (Santa Cruz Biotechnology, sc-515264), anti-p-YAP (Cell Signaling Technology, 4911), anti-YAP (Cell Signaling Technology, 4912), anti-p-TAZ (Cell Signaling Technology, 59971), anti-TAZ (abcam, ab119373), anti-PROX1 (Proteintech, 67438-1-Ig), anti-TEAD (Cell Signaling Technology, 13295), anti-VEGFR3 (Invitrogen, MA5-15651), anti-flag (Invitrogen, MA1-91878), anti-HA (Cell Signaling Technology, 3724) and anti-a-Tubulin (Cell Signaling Technology, 2144) antibodies. a-Tubulin were used as loading controls.

Acetyl-CoA measurement

To measure the nucleus acetyl-CoA concentrations, 50 ml of exponentially growing cells was collected in SB buffer (1 M sorbitol, 20 mM Tris (pH 7.4) and 10 mg/ml zymolyase-20T) and lysed with EBX (20 mM Tris (pH 7.4), 100 mM NaCl, 0.5% Triton X-100, 15 mM 2-mercaptoethanol and protease inhibitors). One aliquot of the lysate was collected to measure the total acetyl-CoA and the remaining lysate was layered over NIB (20 Mm Tris (pH 7.4), 100 mM NaCl, 1.2 M sucrose, 15 mM 2-mercaptoethanol and protease inhibitors). After centrifugation, the pellet was lysed with 1% Triton X-100, centrifuged and collected as the nuclear fraction. The Acetyl-CoA quantity was measured using an Acetyl-CoA kit (Suzhou Comin Biotechnology) according to the manufacturer’s instructions. The assay was performed in 96-well clear-bottom plates, and the fluorescence was quantified using the SpectraMax M2 (Molecular Devices). The acetyl-CoA levels were normalized to the corresponding protein concentrations. For the total and nuclear fractions, the quantity of acetyl-CoA levels was normalized to the corresponding protein concentrations. Histone H3 were analysed by western blot to ensure that there was no cross-contamination between compartments during subcellular fractionation.

Streptavidin affinity purification of dCas9-captured DNA and proteins

Streptavidin affinity purification of dCas9-Captured DNA and proteins was performed following the published protocol [28]. Briefly, Tu212 (5 × 10⁷) cells transfected with FB-dCas9 plasmid and EPHB2 promoter sgRNAs were harvested and treated with 1% formaldehyde to cross-link the proteins to DNA. The genomic locus-associated proteins were isolated using streptavidin purification and further analyzed by Western blot.

Enzyme-linked immunosorbent assay (ELISA)

Quantitative detection of EPHB2 expression in the cell culture medium was examined by a Human Ephrin Type-B Receptor 2 (EPHB2) ELISA Kit (CUSABIO, Wuhan, China), according to the manufacturer’s instructions. Appropriate sample dilutions were made if needed. Each sample, standard, and control were assayed in triplicate. The optical density (OD value) of each well was determined using a microplate reader set to 450 nm.

Tube formation of LECs

The LEC tube formation assay was performed by first pipetting 200 µl 50% Matrigel (BD Biosciences, Bedford, MA) into a 24-well plate, which was then polymerized for 30 min at 37 °C. LECs (2 × 10⁴) in 200 µl of CM were added to each well and incubated at 37 °C, 5% CO₂ for 12 h. Images were taken using a bright-field microscope (Eclipse 80i, Nikon, Tokyo, Japan) at ×100 magnification. The capillary tubes were quantified by measuring the total length of completed tubule structures.

8×GTIIC luciferase reporter activity assay

The 8×GTIIC luciferase reporter, which harbors eight TEAD-binding sites (5’-ACATTCCA-3’), was utilized to indicate the transactivation YAP/TAZ. After transfection with indicated siRNAs per well for 12 h, LECs cells were further transfected with 8×GTIIC-luciferase reporter and pRL-CMV using Lipofectamine 3000. Cells were lysed at 36 h and subjected to a plate-reading luminometers for measurement and calculation the Firefly/Renilla ratios.

Co-immunoprecipitation (Co-IP) assay

Cell lysates were generated from the indicated cells utilizing lysis buffer (150 mM NaCl, 10 mM HEPES, 0.1% NP-40, pH = 7.4). These lysates were subsequently subjected to overnight incubation at 4 °C with indicated antibodies or conjugated beads. Following incubation, beads containing proteins bound to the antibodies were thoroughly washed six times with IP wash buffer (150 mM NaCl, 10 mM HEPES, 0.1% NP-40, pH = 7.4), followed by elution using 1 M glycine at pH 3.0. The resulting eluates were mixed with sample buffer, denatured, separated using SDS-PAGE, and utilized for subsequent western blotting analysis. The antibodies and conjugated beads used in this study were as follows: anti-EFBN1 (Santa Cruz Biotechnology, sc-515264), anti-flag M1 agarose affinity gel (sigma, MA1-91878) and anti-HA immunoprecipitation kit (sigma, IP0010).

Patient-derived organoids (PDOs) culture

Freshly excised laryngeal tumor tissue of patients with lymphatic metastasis were harvested and dissociated enzymatically for 30–40 min in digestion media containing 2 mg/ml collagenase I (Sigma-Aldrich). Supernatant was washed and filtered through a 100 μm cell strainer. Cells were then embedded in basement membrane extract (BME, #E0282, Sigma Aldrich) and solidified in 48-well plates. According to the manufacturer’s instructions, patient-derived tumor organoid culture medium (Accurate international biotech), containing with 10 µM Y‑27,632, 100 ng/mL FGF‑10, 5% FBS and 2% B-27 supplement, was administered immediately and refreshed every 2–3 days. Organoids were collected for western blot and ELISA analyses after 3 to 4 weeks of culture. Subsequently, organoids were subcultured approximately every 14–21 days. All tissues were collected from patients with laryngeal cancer who have undergone surgical treatment in Guangdong Provincial People’s Hospital. The use of human specimens was carried out in accordance with the Helsinki declaration and approved by the ethics committee of Guangdong Provincial People’s Hospital.

Xenograft tumor models and treatments

BALB/c-nude mice (Male, 4–5 weeks old) and NOD-SCID (Male, 4–5 weeks old) obtained from Gempharmatech Co., Ltd. were used for examination of the LN metastasis of Tu212 tumors and the growth of PDOs, respectively. Mice were housed in barrier facilities on a 12 h light/dark cycle. The inguinal LN metastasis model was performed as previously reported [29]. Recombinant human EFNB1-Fc was purchased from R&D systems (7654-EB) and SAR131675 was purchased from MedChemExpress (MCE, HY-15458). EFNB1-Fc (2.5 mg/kg, twice per week) [30] or SAR131675 (100 mg/kg, once per day) [31] were administered 10 days after tumor inoculation. For the inguinal LN metastasis model, after a month of inoculation, mice were euthanized. Inguinal LNs were first measured to calculate volumes. LNs were used for quantitative reverse transcription polymerase chain reaction (qRT-PCR) analysis of hCK and mACTB to determine the proportion of LN-spread Tu212 cells. To construct the organoid xenograft model, the dome solidified by BME-cell solution was digested by collagenase to get a single-cell suspension. Then, mice were subcutaneously injected with 5 × 10⁵ indicated single-cell suspension respectively (n = 6 per group), and tumor volumes were measured weekly. Mice were euthanized and tumors were collected after 7 weeks of inoculation, and IHC staining of lymphatic endothelial cell markers LYVE-1 to indicated the intra-tumoral lymphangiogenesis. All experimental procedures were approved by the Institutional Animal Care and Use Committee of Sun Yat-sen University and performed following the Declaration of Helsinki.

Statistical analysis

All statistical analyses were conducted using R software, GraphPad Prism 6.0 (GraphPad Inc., San Diego, CA, USA), or SPSS 16.0 (IBM, Chicago, IL, USA). Kaplan-Meier analysis with log-rank test was performed to assess survival difference between subtypes. Chi-square test was carried out to determine the difference of clinical and molecular features between subtypes. The variables with P < 0.05 were select to further perform Chi-square test between Cluster D and non-Cluster D group. For comparisons of three groups, one-way analysis and ANOVA test of variance were used. P < 0.05 was considered statistically significant.

Results

Consensus clustering identifies four subtypes in HNSCC based on the metabolic features

The workflow was shown in Fig. 1A. A total of 725 patients were enrolled in our study: 532 from TCGA were assigned to the discovery cohort, 193 from the GEO database were assigned to the validation cohort. To explore a novel molecular classification of HNSCC patients, we performed unsupervised consensus clustering in the TCGA cohort based on the expression patterns of 1064 metabolism-relevant genes from MSigDB (Supplementary Table. S2). To select the optimal and stable cluster number, cumulative distribution function (CDF) curve was calculated to determine the optimal k value, and k = 4 was eventually chosen as an optimal number of clusters after comprehensive consideration (Supplementary Fig. 1A-G), four subclasses were designated Cluster A, Cluster B, Cluster C, and Cluster D. Following this, tSNE (t-Distributed Stochastic Neighbor Embedding) analysis also demonstrated that HNSCC patients under the 4 classifications had significantly different expression profiles (Fig. 1B). To further validate the assignments of subtypes, we first corrected the batch effects across datasets from GEO database and examined the expression patterns via principal component analysis (Supplementary Fig. 1H-I). Subsequently, leveraging TCGA cohort-based grouping data, we established a random forest model to classify patients from GEO into the same four metabolic subgroups (Fig. 1C). Of note, a significant prognostic difference in recurrence-free survival (RFS) and overall survival (OS) was observed among subtypes, especially with shorter RFS for Cluster D than others and better RFS for Cluster A than others (RFS, P = 0.013; OS, P = 0.0036, Fig. 1D-E).

Fig. 1

Identification distinct metabolic subtypes in HNSCC through metabolism gene profiling. A. a Flow chart shows the computational method to classify tumor samples into metabolic subtypes. TCGA cohort was used as a discovery set, GEO cohort was collected as validation sets. B-C. Principal component analysis (PCA) supported the stratification into four HNSCC subclasses in TCGA cohort (B) and GEO cohort (C). D. Recurrence-free survival (RFS) of four subclasses (Cluster A, Cluster B, Cluster C and Cluster D) in TCGA cohort. E. Overall survival (OS) of four subclasses (Cluster A, Cluster B, Cluster C and Cluster D) in TCGA cohort. F. Heatmaps of differential enrichments of metabolism-related signatures in the TCGA cohort. Lipid, carbohydrate, amino acid, nucleotide and vitamin metabolism signatures were presented

To further explored the different metabolic characteristics in distinct subclasses, GSVA was conducted to quantify the enrichment levels of multiple metabolism processes which were divided into five main categories (Lipid metabolism, Carbohydrate metabolism, Amino acid metabolism, Nucleotide metabolism, and Antioxidant metabolism) and the enrichment patterns were demonstrated as a heatmap in Fig. 1F. We found that higher levels of metabolism processes in the Cluster D subtype were mainly carbohydrate, amino acid, nucleotide and antioxidant, which indicated that tumors of this subtype might have enhanced energy metabolism as well as proliferation and viability under oxidative stress. Cluster B displayed higher levels of metabolism processes that were related to lipid, nucleotide and antioxidant, whereas amino acid and carbohydrate metabolism were relatively enriched in Cluster C. Cluster A subtype, by contrast, showed few enrichments of metabolism signatures, which implied low metabolic activities.

Cluster D with upregulated carbohydrate metabolism signatures was associated with lymphatic metastasis in HNSCC

To assess the clinical relevance of the metabolic subtypes identified above, we next determined the component ratio of clinicopathologic features in four subtypes. In the TCGA cohort, these four subtypes had significantly different component in histopathology grade (P < 0.001), HPV infection (P < 0.001), smoking habit (P < 0.001), gender (P < 0.001), and age (P = 0.08) (Supplementary Fig. 2A and Supplementary Table. S3). Similarly, the different component in advanced clinical stage (P < 0.001), smoking habit (P < 0.001) and age (P < 0.001) was also observed in GEO cohort (Supplementary Fig. 2A and Supplementary Table. S4). In particular, we found that the Cluster D subtype was significantly associated with advanced N-stage in both cohorts, which is generally associated with a poorer prognosis in head and neck cancer patients (TCGA Cohort, P = 0.0001; GEO Cohort, P < 0.0001, Fig. 2A-B and Supplementary Table. S5). In addition, the GSVA score of lymphangiogenesis-related gene sets showed that Cluster D possessed molecular features which highly related to lymphatic metastasis, while Cluster A was relatively deficient in these features (all P < 0.05, Fig. 2C-D).

Fig. 2

Connective mapping found that EPHB2 was the key gene of Cluster D and associated with lymphatic metastasis. A-B. Correlation of the HNSCC subclasses with N-stage in the TCGA cohort (A) and GEO cohort (B). C-D. Boxplot of the signature score for lymphangiogenesis-related signatures distinguished by different subclasses in the TCGA cohort (C) and GEO cohort (D). E. A Flow chart shows the connective mapping method to find the key gene of different metabolic subtypes. F. The relative expression level of EPHB2 in different metabolic subtypes. G. Correlation analyses between Cluster D and EPHB2 expression in HNSCC patient specimens. H. Correlation analyses between lymphatic metastasis and EPHB2 expression in HNSCC patient specimens. I. Boxplot of the signature score for lymphangiogenesis-related signatures distinguished by different EPHB2 expression levels. J. Heatmaps of differential enrichments of metabolism-related signatures in high- or low-EPHB2 groups

EPHB2 was the key driven gene of the cluster D subtype in HNSCC

Focusing on the subcluster with the most potential for lymphatic metastasis, we further conducted the connective mapping to identify the key gene for the subcluster in HNSCC (Fig. 2E). The key differentially expressed genes (DEGs) were exhibited in Supplementary Fig. 2B-C, and EPHB2 was the top upregulated-gene of DEGs, which was specifically higher expressed in Cluster D (Mean ± SD, 50.32 ± 18.76) than others (33.74 ± 18.92 of Cluster A, 35.65 ± 17.75 of Cluster B, and 28.09 ± 18.45 of Cluster C, P < 0.0001, Fig. 2F). Notably, a significantly high expression of EPHB2 was also observed on the HNSCC patients in Cluster D compared to those patients in non-Cluster D (75.0% vs. 37.9%, P < 0.0001, Fig. 2G). Moreover, correlation analysis indicated that the high EPHB2 expression was significantly associated with lymphatic metastasis of HNSCC (LMs vs. non-LMs, 57.6% vs. 42.5%, P = 0.0005, Fig. 2H). The GSVA score of lymphangiogenesis-related gene sets also showed that high EPHB2 expression closely related to lymphatic metastasis (all P < 0.05, Fig. 2I). In addition, the specific metabolism signatures of high-expressed EPHB2 group were similar to the metabolic patterns of previously reported in Cluster D subtype, involving gene signatures of carbohydrate, amino acid and nucleotide metabolism (Fig. 2J). These findings indicated that EPHB2 was defined as the key gene of the Cluster D subtype and may played a predominant role in lymphatic metastasis of HNSCC patients.

Enhanced carbohydrate metabolism upregulated EPHB2 in HNSCC

Indeed, compared with normal epithelial cell from nasopharynx and larynx, we found that EPHB2 was widely upregulated in HNSCC cell lines (Fig. 3A). Consistently, a marked upregulation of EPHB2 in the HNSCC tissues with lymphatic metastasis, corroborating our earlier preclinical observations (Fig. 3B). These results suggested that EPHB2 may act as a potential key factor involved in mediating lymphatic metastasis in HNSCC. Previous studies have been reported that carbohydrate metabolism feeds into the regulation of gene expression via metabolic enzymes and metabolites, of which pyruvate metabolism was regarded as a central part. Therefore, we detected whether pyruvate contributed to EPHB2 upregulation. Indeed, exogenously adding pyruvate upregulated the levels of EPHB2 expression in Tu212 and Tu686 cells (Fig. 3C). In addition, TM-1, a blocker of pyruvate dehydrogenase kinase (PDHK), treatment abolished pyruvate-induced increase of EPHB2 expression (Fig. 3C), suggesting that pyruvate metabolism played an essential role in EPHB2 upregulation, which may result from the accumulation of acetyl-CoA, the major metabolite of pyruvate. Acetyl-CoA was a substrate for acetylation reactions and directly affects the function of all acetyltransferases [32,33,34]. Acetyl-CoA produced by pyruvate was localized inside mitochondria and needed to be exported for nuclear histone acetylation as citrate or acetate, which was converted back to acetyl-CoA in the cytosol by ATP citrate lyase (ACLY) or acetyl-CoA synthase (ACSS2) (Fig. 3D). Given the upregulated pyruvate metabolism signatures in EPHB2-high cells and Cluster D subtype, we found that pyruvate treatment increased the acetyl-CoA levels in nucleus, and eventually enhanced histone H3 and H4 acetylation, which resulted in the augment of EPHB2 transcription levels. Whereas, depletion of ACLY and ACSS2 abolished pyruvate-induced EPHB2 transcription activation by decreasing the levels of nuclear acetyl-CoA content and histone acetylation (Fig. 3E-H). Similarly to pyruvate treatment, glucose significantly augmented, while 2-Deoxy-D-glucose (2-DG), a competitive inhibitor of glycolysis to reducing pyruvate production, treatment decreased nuclear acetyl-CoA content (Fig. 3I). To further assess the levels of histone acetylation in EPHB2 promoter, the CRISPR affinity purification in situ of regulatory elements (CAPTURE) approach [28] were used to enrich histone H3 and H4 targeting EPHB2 in Tu212 (Fig. 3J). We found that histone H3 and H4 acetylation were increased on the EPHB2 promoter under pyruvate and glucose treatment (Fig. 3J). The p300, H3K27ac and RNA polymerase II (pol II) enriched on the EPHB2 promoter under pyruvate and glucose treatment, which substantially activated EPHB2 transcription (Fig. 3K-L). Whereas, 2-Deoxy-D-glucose (2-DG) suppressed glycolysis-induced EPHB2 transcriptional activation (Fig. 3K-L). These results demonstrated that highly activated carbohydrate metabolism with increase of acetyl-CoA production transcriptionally upregulated EPHB2 expression via promoting histone acetylation in HNSCC.

Fig. 3

Enhanced carbohydrate metabolism upregulated EPHB2 in HNSCC. A. The mRNA (upper) and protein (lower) expression levels of EPHB2 in the normal epithelial cell line from nasopharynx (NP69) and larynx (HuLa-PC) and HNSCC cell lines. a-Tubulin was used as a loading control. B. The protein expression levels of EPHB2 in lymphatic metastatic HNSCC and non-lymphatic metastatic HNSCC. a-Tubulin was used as a loading control. C. qRT-PCR (upper) and western blot analysis (lower) of EPHB2 under different concentration of pyruvate with or without TM-1. a-Tubulin was used as a loading control. D. The source of nuclear acetyl-CoA. E. Nuclear acetyl-CoA levels in the indicated groups. F. The levels of histone H3 and H4 acetylation. G. The luciferase activities of EPHB2 reporter in the indicated groups. (H) qRT-PCR (upper) and western blot analysis (lower) of EPHB2 in the indicated groups. (I) Nuclear acetyl-CoA levels under pyruvate, glucose and 2-DG treatment. (J) Schematic (left) of the CRISPR CAPTURE approach following the western blot (right) to assess the levels of histone H3 and H4 acetylation in EPHB2 promoter in LECs under pyruvate, glucose and 2-DG treatment. K. ChIP assays examining the enrichment of P300, H327ac and Pol II on the EPHB2 promoter under pyruvate, glucose and 2-DG treatment. L. qRT-PCR (upper) and western blot analysis (lower) of EPHB2 under pyruvate, glucose and 2-DG treatment. NC, non-targeting negative control

EPHB2 promoted lymphatic metastasis

To investigate whether EPHB2 was involved in lymphatic metastasis, we performed knockdown and overexpression modifications in the HNSCC cell lines with a high-level of EPHB2 expression (Fig. 3A). Western blot and ELISA assays showed that overexpressing EPHB2 significantly increased, while silencing EPHB2 decreased the abundance of EPHB2 both in whole cell lysis and supernatant (Fig. 4A-B), indicated that EPHB2 contributed to lymphatic metastasis as a secretory protein. Notably, the conditioned medium (CM) from EPHB2-overexpressing Tu212 and Tu686 cells significantly increased, while which from silencing EPHB2 cells suppressed lymphangiogenesis in vitro, as indicated by the migration and tube formation of LECs (Fig. 4C-D). We further injected with control, EPHB2-overexpressing or EPHB2-silengcing Tu212 cells at the footpad to investigate the role of EPHB2 in lymphatic metastasis in vivo. After 45 days of inoculation, we sacrificed the mice, extracted, and analyzed the primary footpad tumors and popliteal lymph nodes. Results showed that EPHB2 promoted, while depletion of EPHB2 inhibited, LN metastasis as indicated by the volumes of LNs and the number of metastatic Tu212 cells in LNs (Fig. 4E-F). Additionally, the relative mRNA ratio of human CK to mouse ACTB indicated a higher proportion of colonized tumor cells in LNs in the EPHB2-ovexpressing group, which was significantly decreased in the EPHB2-silencing group (Fig. 4G-H). These results suggested that EPHB2 was contributed to lymphatic metastasis in laryngeal cancer.

Fig. 4

EPHB2 promoted lymphatic metastasis. A. Western blot analysis of EPHB2 in whole cell lysate (WCL) and supernatant (Sup) in control, EPHB2-silenced, and EPHB2-overexpressing Tu212 and Tu686 cells. B. ELISA analysis of EPHB2 levels in cell culture medium of control, EPHB2-silenced, and EPHB2-overexpressing Tu212 and Tu686 cells. C. LECs were incubated with the indicated medium and then subjected to transwell migration assays. Scale bars: 100 μm. D. Tube-formation assays of LECs incubated with the indicated medium. Scale bars: 200 μm. E. Image of the inguinal lymph nodes (LNs) from each group. F. Quantification of LNs volumes in each group. G. qRT-PCR analysis of hCK relative to mACTB in the LNs from each group. The ratio indicated the proportion of metastatic cells. H. The analysis of metastasis LNs ration in each group. CM, conditional medium

EPHB2 inhibited HIPPO signaling activation to triggering lymphangiogenesis via EFNB1-induced YAP/TAZ cytoplasmic retention

Interestingly, VGX-100, a highly specific human monoclonal antibody for VEGF-C, or SAR131675, a selective VEGFR3 inhibitor, treatments only slightly impaired EPHB2-induced lymphangiogenesis (Fig. 5A). Similarly, silencing VEGFR3 had slight effect on EPHB2-induced lymphangiogenesis (Fig. 5B and Supplementary Fig. 3A). These results suggested that EPHB2 may promote lymphangiogenesis in a VEGF-C independent manner. Recently, several studies have identified Hippo-YAP/TAZ signaling components as novel players in lymphatic vascular development by regulating LECs specification, differentiation, and sprouting during early lymphatic development and maintaining lymphatic integrity during adulthood, which modulate lymphatic plasticity throughout life via regulating PROX1 activity [35]. Next, the 8×GTIIC luciferase reporter, which harbors eight TEAD-binding sites, was utilized to indicate the transactivation of YAP/TAZ in LECs, and we found that EPHB2 treatment dramatically decreased the YAP/TAZ transcriptional activity on TEAD-binding elements in LECs (Fig. 5C). EPHB2 was the member of Eph receptors whose extracellular signal transduction relied on the binding of their cognate membrane-tethered ligands, known as ephrin B ligands [36]. As shown in Fig. 5D, only silencing EFNB1 did not impair YAP/TAZ transcriptional activity on TEAD-binding elements, indicated that EPHB2 affected on YAP/TAZ transcriptional function in LECs via binding to EFNB1. Furthermore, western blot analysis showed that EPHB2 enhanced YAP and TAZ phosphorylation, which lose the function of nuclear translocation, and reduced PROX1 expression (Fig. 5E). PROX1, a master transcription factor of inducing a lymphatic phenotype, was reported to transcriptional upregulation of the lymphatic markers [37], such as LYVE-1, podoplanin (PDPN) and secondary lymphoid chemokine (SLC), and the cell-matrix receptor integrin alpha 9 (ITGA9) [38]. Indeed, we found that EPHB2 increased the expression of LYVE-1, PDPN, SLC and ITGA9 (Fig. 5F and Supplementary Fig. 3B). However, depletion EFNB1 of LECs inhibited EPHB2-induced YAP phosphorylation and PROX1 upregulation, resulted in decrease of LYVE-1, PDPN, SLC and ITGA9 expression (Fig. 5E-F and Supplementary Fig. 3B). In addition, the nucleus-localized YAP-5SA mutant, which was resistance to YAP phosphorylation, robustly suppressed EPHB2-induced upregulation of PROX1 and its downstream target genes (Fig. 5G-H), which also impaired EPHB2-medicated LECs migration and tube formation (Fig. 5I). Furthermore, immunofluorescence (IF) staining revealed that EPHB2 enhanced YAP cytoplasmic localization, while YAP showed more nuclear localization in EFNB1-silencing LECs (Fig. 5J). Co-IP assays indicated that EPHB2 promoted the interaction between EFNB1 and YAP/TAZ, while neither YAP nor TAZ interacted with EFNB1 without EPHB2 treatment (Fig. 5K-L), suggesting that EPHB2 inhibited YAP/TAZ nuclear translocation via promoting the interaction between EFNB1 and YAP/TAZ. To further investigate whether EPHB2-mediated YAP cytoplasmic localization relied on EFNB1 phosphorylation, EFNB1-Y317A/Y331A mutants were transfected into the LECs, in which both Y317 and Y331 were EPHB2-induced tyrosine phosphorylation sites [39]. Co-IP assays showed that dephosphorylated EFNB1 had no effect on the combination with EPHB2, while abolished the interaction between EFNB1 and YAP or TAZ (Fig. 5M). And the EFNB1-Y317A/Y331A mutants inhibited EPHB2-induced YAP and TAZ phosphorylation and PROX1 upregulation (Fig. 5N). Take together, these results demonstrated that EPHB2 binding with EFNB1 promoted YAP/TAZ cytoplasm retention, resulted in alleviating YAP/TAZ activation-induced PROX1 transcriptional repression in LECs, finally triggered lymphangiogenesis.

Fig. 5

EPHB2 inhibited HIPPO signaling activation to triggering lymphangiogenesis via EFNB1-induced YAP/TAZ cytoplasmic localization A. LECs were incubated with the indicated medium and then subjected to transwell migration assays (upper) and tube-formation assays (lower) under treated with VGX-100 (30 µg/mL) or SAR131675 (40 nM). B. Control and VEGFR3-silenced LECs were incubated with the indicated medium and then subjected to transwell migration assays (upper) and tube-formation assays (lower). C. The 8×GTIIC luciferase reporter assays of Hippo pathway in LECs with or without EPHB2 addition. D. The 8×GTIIC luciferase reporter assays of Hippo pathway in LECs with silencing EFNB1, EFNB2, and EFNB3, respectively. E. Western blot of the expression of critical genes of Hippo signaling pathway and PROX1 in control and EFNB1-silenced LECs incubated with control or EPHB2-overexpressed Tu212 and Tu686 cells. a-Tubulin was used as a loading control. F. The mRNA levels of LYVE-1, PDPN, SLC and ITGA9 in control and EFNB1-silenced LECs incubated with control or EPHB2-overexpressed Tu212 cells. G. Western blot of the expression of phosphorylated-YAP and PROX1 in indicated groups. H. The mRNA levels of LYVE-1, PDPN, SLC and ITGA9 in LECs transfected with YAP-Wt or YAP-5SA. I. YAP-Wt or YAP-5SA-transfected LECs were incubated with or without rEPHB2 treatment and then subjected to transwell migration assays (upper) and tube-formation assays (lower). J. The immunofluorescence (IF) staining of YAP was performed in indicated groups. K. Immunoprecipitation (IP) assays of EFNB1, EPHB2, YAP, TAZ and TEAD in LECs with or without rEPHB2 treatment. L. Co-IP assays of EFNB1, EPHB2, YAP or TAZ in LECs transfected with Flag-YAP or HA-TAZ. M. Co-IP assays of EFNB1, EPHB2, YAP and TAZ in LECs transfected with EFNB1-Wt or EFNB1-Y317A/Y331A. N. Western blot of the expression of critical genes of Hippo signaling pathway and PROX1 in LECs transfected with EFNB1-Wt or EFNB1-Y317A/Y331A. NC, non-targeting negative control

Combined with EFNB1-Fc promoted VEGFR3 inhibitor efficiency in laryngeal cancer with lymphatic metastasis

Finally, we assessed the therapeutic effect of blocking EPHB2 by EFNB1-Fc on Tu212 and Tu686 cells with lymphatic metastasis. As shown in Fig. 6A-B, EFNB1-Fc reduced the migration and tube formation of LECs in CM derived from EPHB2-overexpressing Tu212 and Tu686 cells. Further, we injected 2 PDOs from laryngeal cancer tissues with lymphatic metastasis to NSG mice. The patients’ clinical information was showed in Supplementary Fig. 4A, and EPHB2 expression of the two PDOs was measured by western blot and ELISA assays, which levels was comparable with Tu212 and Tu686 cells (Supplementary Fig. 4B). The administration of EFNB1-Fc slightly abolished tumor volumes, while dramatically decreased intra-tumoral lymphangiogenesis in vivo as indicated by IF and immunohistochemical (IHC) staining of LEC markers LYVE-1 (Fig. 6C-D). The therapeutic effect of combined EFNB1-Fc and SAR131675 on lymphatic metastasis in laryngeal cancer was further validated using in vivo mouse models. Mice were randomly divided into four groups, and equally injected Tu212 cells at the footpad. After 10 days of tumor inoculation, we then treated them with either PBS, EFNB1-Fc, SAR131675 or EFNB1-Fc + SAR131675. After a month of treatment, mice were euthanized and the metastasis in LNs was analyzed. Remarkably, the coadministration of EFNB1-Fc and SAR131675 dramatically shrank the volumes of LNs and the number of metastatic cancer cells in LNs compared with EFNB1-Fc or SAR131675 monotherapy (Fig. 6E-G).

Fig. 6

Combined with EFNB1-Fc promoted VEGFR3 inhibitor efficiency in laryngeal cancer with lymphatic metastasis. A. LECs were incubated with the indicated CM with or without EFNB1-Fc treatment and then subjected to transwell migration assays. Scale bars: 100 μm. B. Tube-formation assays of LECs incubated with the indicated CM with or without EFNB1-Fc treatment. Scale bars: 200 μm. C. The tumor volume of 2 patients-derived organoids (PDO) treated with or without EFNB1-Fc. D. The IHC staining of lymphatic endothelial cell markers LYVE-1 was performed in tumors from 2 PDO. E. Image (left) and quantification (right) of LNs volumes in each group. F. qRT-PCR analysis of hCK relative to mACTB in the LNs from each group. The ratio indicated the proportion of metastatic cells. G. The analysis of metastasis LNs ration in each group. H. Study model: pyruvate-mediated EPHB2 upregulation promotes lymphatic metastasis in HNSCC

Therefore, our results demonstrated that, in laryngeal cancer, highly activated pyruvate metabolism increase acetyl-CoA production, and then transcriptionally upregulated EPHB2 levels, which inhibited YAP/TAZ access to the nucleus and resulted in augment of PROX1 expression, eventually induced lymphangiogenesis (Fig. 6H). In addition, combined treatment with EFNB1-Fc and SAR131675 exerted synergistic effects on blocking lymphatic metastasis, suggesting that targeting EPHB2 might be a potential strategy to patients who do not respond to VEGFR3 inhibitor.

Discussion

In this study, we firstly established a new classification of HNSCC patients based on metabolism gene expression profiles related to lymphatic metastases. Four metabolic subtypes were identified, and their clinical characteristics, metabolic signatures were also explored. Especially Cluster D, which has a higher level of carbohydrate metabolism, was highly related to lymphatic metastasis and significantly associated with a worse RFS outcome. We found that EPHB2 was the key gene of the Cluster D subtype, in which highly activated carbohydrate metabolism with increase of acetyl-CoA production transcriptionally upregulated EPHB2 expression via promoting histone acetylation. Additionally, tumor-derived EPHB2 promoted lymphatic metastasis both in vitro and in vivo. Interestingly, we found that LECs cultured in EPHB2-overexpressed cells conditional medium were relevant to resistance of multiple VEGFR3 inhibitors, indicating that EPHB2 promoted lymphangiogenesis in a VEGF-C independent manner. Furthermore, EPHB2 activated EFNB1 to promote YAP/TAZ cytoplasmic retention, resulted in inhibiting HIPPO activation-induced PROX1 transcriptional repression, finally triggered lymphangiogenesis.

Recently, genetic signature profiling based on high-throughput sequencing date has become the informative approach to elucidate tumor heterogeneity [40, 41]. Previously published study recruited a cohort of 9125 TCGA samples across 33 cancer types and characterised tumor subtypes based on the expression of seven metabolic pathways [42]. Some researchers have identified and validated multiple highly distinct metabolic subtypes in pancreatic ductal adenocarcinoma and lower-grade glioma patients, developing a metabolic signature with better performance of prognosis prediction[43, 44]. In this study, based on the metabolite profiling, we classified the HNSCC patients into four subtypes, and found that Cluster D with mainly enriched in carbohydrate metabolism has the worst RFS and was highly related to lymphatic metastasis. Given that the current deficiency of rapid and convenient approaches to detect the metabolic characteristics of tumors in clinic, we further employed bioinformatics methods to screen out that EPHB2 as the key gene of Cluster D. Moreover, we observed the correlation between EPHB2 and lymphatic metastasis in multiple HNSCC tissues, which discovered that EPHB2 may be regarded as a potential biomarker to indicate lymphatic metastasis potential and a novel therapeutic target against lymphatic metastasis in HNSCC.

Meanwhile, metabolic reprogramming and epigenetic modifications are the hallmarks of cancer cells [45, 46]. Previous studies showed that epigenetic modifications were closely related to cancer cell metabolism regulations, such as histone methylation, acetylation, DNA methylation and RNA N6-methyladenosine (m6A) methylation. Glycolysis is the main source of energy for tumor cell proliferation, and its main product, acetyl-CoA, also plays an important role in tumor invasion and migration. Previous research has shown that an increase rate of the acetyl-CoA synthesis promotes post-translational histone acetylation resulted in tumor malignant progression [47, 48]. For example, the production of ACLY-dependent acetyl-CoA was proved to play an vital role in early stage pancreatic cancer development, and targeting the acetyl-CoA-dependent pathway inhibited tumor proliferation [49]. Another study revealed that maintaining ACLY expression to upregulate acetyl-CoA levels promoted the proliferation, metastasis, and even drug resistance in nasopharyngeal carcinoma cells [50]. In this study, we also found that pyruvate treatment increase the acetyl-CoA levels in nucleus, which promotes H3 and H4 acetylation to transcriptional upregulate EPHB2. While depletion of ACLY and ACSS2 abolished pyruvate-induced EPHB2 transcription activation by impairing nuclear acetyl-CoA production, demonstrating that the subtype with highly activated carbohydrate metabolism promotes HNSCC lymphatic metastasis via transcriptional upregulation of EPHB2.

In previously published studies related to lymphangiogenesis, VEGF-C and VEGF-D played significant roles by activating its receptor VEGFR3 on LECs to activate a protein kinase C/ERK signaling cascade, which ultimately triggers the phosphorylation of AKT and contributes to the proliferation and migration of LECs [51]. Additionally, there are many enzymes, bioactive lipids, chemokines, adhesion molecules, and noncoding RNAs that also participate in lymphangiogenesis by functioning in either a VEGF-C/D-dependent or -independent manner. Recently, several studies have identified Hippo-YAP/TAZ signaling components as novel players in lymphatic vascular development. YAP/TAZ depletion or hyperactivation in LECs during embryonic development results in structurally aberrant and poorly functional lymphatics and lethality, highlighting its importance in early development of the lymphatic system [35, 52]. In this study, we found that EPHB2 promoted lymphangiogenesis in a VEGF-C independent manner, which binding to EFNB1 dramatically relieved the YAP/TAZ-mediated transcriptional repression on PROX1 in LECs and then promoted the expression of lymphatic markers and the cell-matrix receptor, eventually triggered lymphangiogenesis in HNSCC. Moreover, EPHB2-mediated YAP/TAZ cytoplasmic localization relied on EFNB1 phosphorylation. However, the detail mechanism of EPHB2/EFNB1 sequester YAP/TAZ in the cytoplasm remains unclear and whether this process required other proteins and involved phase separation still needs to be further explored.

Although the clinical application of single-target VEGFR3 inhibitor in HNSCC treatment was lack, the application of multi-target tyrosine kinase inhibitors containing targeting VEGFR3 has only achieved certain efficacy in a small group of patients with tumor [53,54,55]. For example, Apatinib has shown modest antitumor activities with manageable side effects in the recurrent or metastatic HNSCC patients [54]. But still some patients couldn’t benefit from the VEGFR inhibitor, and need other treatment strategy. In this study, we proved the combination treatment of EFNB1-Fc and SAR131675 could dramatically shrank the volumes of LNs and the number of metastatic cancer cells in mice, suggesting that the combination therapy may have good clinical application prospects especially for those HNSCC patients who do not respond to VEGFR3 inhibitor. Given that the wild expression of ENFB1 and VEGFR3 in normal tissues such as lung, thyroid, breast, prostate, bone marrows, we still need further study to investigate the side effects of EFNB1-Fc and SAR131675, especially when they used together.

Conclusion

In summary, our study firstly reclassified HNSCC from the metabolic perspective and proposed four subtypes with distinct prognosis and metabolic phenotype. Among these four subtypes, Cluster D was highly related to lymphatic metastasis and significantly associated with a worse RFS outcome, which has a higher level of carbohydrate metabolism. Then we demonstrated that this subtype with highly activated pyruvate metabolism increased acetyl-CoA production and then transcriptionally upregulated EPHB2 levels, which promoted YAP/TAZ cytoplasmic localization and resulted in augment of PROX1 expression, eventually induced lymphangiogenesis. Our findings extended the molecular subtyping of HNSCC patients and showed that combination of EFNB1-Fc and VEGFR3 inhibitor exerted synergistic effects on blocking lymphatic metastasis, suggesting that targeting EPHB2 might be a potential strategy to patients who do not respond to VEGFR3 inhibitor.