A comprehensive dataset and clinical samples for 76 pairs of tumor and adjacent non-tumor liver tissues was collected from January 2020 to September 2021 at the Affiliated Hospital of Qingdao University (Supplementary Table 1). The collection of specimens and associated research procedures were approved by the Research Ethics Committee of the Affiliated Hospital of Qingdao University (Approval Identifier: QYFYWZLL26589). Informed consent was obtained from participants or their legal representatives before the start of any research activity. For the animal experiments, 5-week old male Balb/c nude mice were procured from SPF Biotechnology (Beijing, China), ensuring compliance with ethical guidelines. The mice were maintained in an individually ventilated caging (IVC) system to promote their well-being. Subsequently, 5 $\times$ 10⁶ enhanced green fluorescent protein labeled lentivirus (EGFP-Lv)-infected Huh7 cells were subcutaneously administered to each mouse over a period of two weeks. Oversight was provided by the Research Ethics Committee of Qingdao University (Approval Identifier: 20220309BALB/c-nu3020221009059). Mice were sacrificed using a 1% solution of pentobarbital administered intraperitoneally (50–60 mg/kg). All protocols concerning animal experimentation were rigorously followed at the designated research facility of the Affiliated Hospital of Qingdao University, and in strict accordance with ethical standards.

In addition, 32 liver cancer patients who were administered sorafenib were subjected to computed tomography (CT) imaging and their tissue specimens quantitatively assessed for ENG expression (Supplementary Table 2). Patients with high and low ENG expression in liver cancer were selected, and CT was used to assess changes in tumor diameter before and after treatment with sorafenib. All results were evaluated by three specialized radiologists.

Proteins from cellular extracts and tissue samples were isolated using Radio-Immunoprecipitation Assay (d) solution (R0010, Solarbio, Beijing, China). The protein samples were separated by 8% SDS-PAGE and immediately transferred onto Polyvinylidene fluoride (PVDF) membranes. Prior to the application of antibody, the membranes were first blocked using 5% non-fat milk. For the detection of specific proteins, the membranes were probed with a series of primary antibodies as follows: anti-ENG (1:2000 dilution, ab169545, Abcam, Cambridge, UK), anti- Collagen type I$\alpha$1 (COL1A1) (1:5000 dilution, 67288-1-Ig, Proteintech, Rosemont, IL, USA), anti-CD31 (1:1000 dilution, ab9498, Abcam, Cambridge, UK), anti-Beta Actin (1:20,000 dilution, 66009-1-Ig, Proteintech, Rosemont, USA), and anti-GAPDH (1:50,000 dilution, 60004-1-Ig, Proteintech, Rosemont, USA). Following incubation with primary antibodies, the membranes were subsequently exposed to horseradish peroxidase-conjugated secondary antibody (SA00001-2/SA00001-1, Proteintech, Rosemont, USA). The visualization of protein bands was achieved using an advanced enhanced chemiluminescence (ECL) method, and images were captured using the Tanon-5200 chemiluminescence imaging system. Total RNA was isolated from 76 pairs of tumor and adjacent non-tumor liver tissue samples using TRIzol reagent (R1100, Solarbio, Beijing, China). The RNA was assessed by quantitative real-time PCR (qRT-PCR) with a SYBR Green PCR kit (RK21204, ABclonal, Wuhan, China) and using an ABI Prism 7500 Sequence Detection system (Applied Biosystems, Foster City, USA). Detailed primer sequences for qRT-PCR are listed in Supplementary Table 3. The relative expression levels of genes were calculated using the 2^{-Δ⁢Δ⁢Ct} method. $\mathrm{\Delta }$Ct was calculated as Ct^(ENG) - Ct^(GAPDH), and $\mathrm{\Delta }$$\mathrm{\Delta }$Ct was calculated as $\mathrm{\Delta }$Ct^(tumor) - $\mathrm{\Delta }$Ct^(normal). Graphs were generated using GraphPad Prism 8.0 software (GraphPad Software, Inc., La Jolla, CA, USA).

The cell lines Huh7 and HepG2, representative of liver cancer, as well as human umbilical vein endothelial cells (HUVECs), were purchased from the China Center for Type Culture Collection (CCTCC, Wuhan, China). Cells were grown in Dulbecco’s Modified Eagle Medium (DMEM) containing 10% Fetal Bovine Serum (FBS) and 100$\times$ concentration of Penicillin-Streptomycin (P1400, Solarbio, Beijing, China). The Huh7 cell line was genetically modified to stably express EGFP using an EGFP-Lv. All cell lines were validated by Short Tandem Repeat (STR) profiling and tested negative for mycoplasma. Cells were all cultured in a humidified incubator at 37 °C and 5% CO${}_2$.

Lentiviral shRNA and vectors for the knockdown of ENG (sh-ENG) and overexpression of COL1A1 (OE-COL1A1) were designed and constructed by Genechem (Shanghai, China). The interference sequences are as follows: sh-ENG-1: GATCTGGACCACTGGAGAATATT CAAGAGATATTCTCCAGTGGTCC AGATCTTTTTT; sh-ENG-2: GCTCAGGATGACATGGACATCTT CAAGAGAGATGTCCATGTCATCC TGAGCTTTTTT; sh-ENG-3: CGCCATGACCCTGGTACTAAATT CAAGAGATTTAGTACCAGGGTCA TGGCGTTTTTT; and negative control (NC) shRNA, TTCTCCGAACGTGTCACGTAA TTCAAGAGATTACGTGACACGTTCG GAGAATTTTTT. The ENG knockdown vector is GV248, with the element sequence: hU6-MCS-Ubiquitin-EGFP-IRES-puromycin, and the restriction sites are EcoRⅠ and AgeⅠ. The COL1A1 overexpression vector is GV492, with the element sequence: Ubi-MCS-3FLAG-CBh-gcGFP-IRES-puromycin, and the restriction sites are BamHⅠ and AgeⅠ. The plasmid carries a 3$\times$ FLAG fusion protein. The expression of the target genes was detected by Western blotting.

The density of neovascularization in tumor tissue was evaluated in vivo using an ultrasound system (VINNO 6 LAB, VINNO Technology Co., Ltd., Suzhou, China) designed for small animals and coupled with a laser confocal microscope (LSM 900 & Axio Imager M2, Carl Zeiss AG, Oberkochen, Germany). Following sedation with 0.3% pentobarbital sodium (50–60 mg/kg), mice were positioned on the imaging platform of the ultrasound system. Subsequently, a high-resolution 40 MHz probe designed specifically for small animal research was prepared. The imaging process was begun by activating the software and applying a conductive gel to the tumor region. The ultrasound probe was then finely adjusted for optimal positioning to capture detailed images of neovascularization density.

Frozen tissue sections were preserved in 75% ethanol and permeabilized with 0.1% Triton X-100. To reduce non-specific binding, sections were first blocked with 1% goat serum. The following primary antibodies were applied: anti-EGFP (1:50 dilution, ab184601, Abcam, Cambridge, UK), anti-CD31 (1 µg/mL dilution, ab9498, Abcam, Cambridge, UK), and anti-Alpha-fetoprotein (AFP) (1:50 dilution, 14550-1-AP, Proteintech, Rosemont, USA). After washing with PBS, secondary antibodies conjugated to Fluorescein Isothiocyanate (FITC) or Tetramethylrhodamine Isothiocyanate (TRITC) (1:200 dilution) were used for visualization. Nuclear staining was performed using 5 µL of 4${}^{\prime }$,6-diamidino-2-phenylindole (DAPI) for 10 minutes. Dual immunofluorescence was performed using laser scanning confocal microscopy (A1R HD25, NIKON, Tokyo, Japan).

The IHC protocol followed the instructions given by the manufacturer (ZSGB, Beijing, China). Tissue slices were first subjected to dehydration at 56 °C for 30 minutes and subsequently rehydrated. To minimize non-specific interactions, blocking was initially performed using serum, followed by the addition of 100 µL of diluted (1:100) primary antibody solution. Samples were incubated at 37 °C for one hour, then treated with a biotinylated secondary antibody (goat anti-rabbit or anti-mouse IgG diluted 1:1000) followed by horseradish peroxidase-conjugated streptavidin. 3,3${}^{\prime }$-diaminobenzidine (DAB) solution was applied for 5 minutes to develop color, and images of stained slides were captured under a microscope for analysis. Antibodies: ENG (ab169545, Abcam, Cambridge, UK), CD34 (HPA036723, Sigma-Aldrich, St. Louis, MO, USA), CD3 (ab16669, Abcam, Cambridge, UK), CD8 (66868-1-Ig, Proteintech, Rosemont, USA), CD45 (20103-1-AP, Proteintech, Rosemont, USA).

HUVECs in the logarithmic growth phase were suspended in conditioned medium for accurate counting. A 200 µL cell suspension was added to each Transwell chamber well, with 200 µL of 10% FBS growth medium added to the culture plates beneath. After 12 or 48 h of incubation, the adherent cells were fixed, stained, and imaged to assess their morphology and distribution. Cells were counted in six randomly chosen fields per group using microscopy. This experiment was repeated in triplicate.

Corning Matrigel was carefully dispensed into each well of a 96-well plate (60 µL per well), then placed in an incubator for 30 min to ensure thorough solidification of the Matrigel. Following this, HUVECs were prepared in a specialized conditioned medium. Each well was inoculated with 10,000 HUVECs, with three replicates for each experimental set. After incubation for 6 h, the progression and morphology of cells was evaluated by microscopic examination. Cells were systematically photographed and quantitatively counted.

The RNA sequencing (RNA-seq) dataset for single-cell analysis (GSE149614), along with 9 additional RNA-seq datasets focused on liver cancer, were retrieved directly from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE149614). Concurrently, 33 tumor RNA-seq datasets accompanied by comprehensive clinicopathological characteristics were obtained from The Cancer Genome Atlas (TCGA). The processing and analysis of these single-cell RNA-seq datasets was performed using the Seurat R package (v5.1.0), which facilitated detailed cellular profiling. To understand the dynamic progression and fate determination of tumor cells, the Monocle2 R package (v2.32.0) was employed to efficiently analyze pseudo-time trajectories. The criteria for identifying differentially expressed genes (DEGs) were set with the stringent thresholds of $|$logFold-Change$|$ $>$1 and a p-value 0.05. These were applied to specifically analyze ENG-high and weakly expressed liver hepatocellular carcinoma (LIHC) RNA-Seq data from TCGA (https://portal.gdc.cancer.gov/). Gene Set Enrichment Analysis (GSEA) was performed using the MSigDB to identify DEGs that were statistically significant. The tumor microenvironment (TME) was quantitatively estimated using the Estimate R package to obtain a better understanding of the cellular context. Furthermore, Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) was employed for the detailed assessment of immune infiltration in LIHC. Differences between several variables were analyzed using two-tailed Student’s t-tests, and survival rates were assessed by log-rank tests. Spearman’s rank correlation coefficient was utilized to investigate associations between ENG expression and various clinical parameters, such as StromalScore, ImmuneScore, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data Score (ESTIMATScore), TumorPurity, immune infiltration, the effectiveness of sorafenib, and alterations in tumor size. Log-rank tests were used to evaluate the statistical significance of differences in overall and disease-free survival rates between groups. The outcomes of these analyses were reported as the mean $\pm$ standard deviation. A threshold value of p 0.05 was set for statistical significance, below which findings were deemed to have meaningful implications for prognosis.

A correlation between specific types of liver cancer cells and the trajectory of disease progression is well documented [20]. However, the intricate details of how tumor cells evolve and progress in liver cancer have not been comprehensively elucidated. To address this gap in knowledge, we retrieved the single-cell dataset GSE149614 from the GEO database. Detailed cluster analysis of this dataset revealed the presence of 8 distinct tumor cell subtypes (Supplementary Fig. 1A,B), each identified by unique biomarkers (Supplementary Fig. 1C,D). These were present in both the primary liver cancer tissue and adjacent non-tumor tissue. A pseudo-time analysis was performed to unravel the evolutionary trajectories of each tumor cell subtype. This revealed that during migration from the central tumor mass to the adjacent normal tissue, the tumor cells underwent significant retro-differentiation at three major branch points (Fig. 1A). Retro-differentiation involved substantial changes in the expression profiles of specific markers in the tumor cell subpopulations (Fig. 1B). During this transformative process, a notable increase was observed in the expression of markers typically associated with vascular ECs, including ENG, Alpha-2-Macroglobulin (A2M) and Platelet Endothelial Cell Adhesion Molecule 1 (PECAM1) (Fig. 1C–E). To further substantiate these observations, immunofluorescence staining was conducted on biopsy samples from liver cancer patients. This confirmed the elevated expression of CD31 (n = 6), a marker of vascular ECs, predominantly in the tumor cells located at the edges of the tumor mass (Fig. 1F). The findings suggest that as liver cancer cells migrate from the central tumor region to surrounding normal tissues, they transform into vascular endothelial-like phenotypes, indicating the presence of a complex mechanism of cancer progression and metastasis.

Fig. 1.

Retro-differentiation of tumor cells to vascular endothelial-like cells during the process of infiltrating para-carcinoma. (A) During infiltration from within the tumor to adjacent normal tissue, the tumor cellsretro-differentiate at three distinct branch points. (B) Heat map showing that specific markers of the subpopulation are reprogrammed during this retro-differentiation process. Markers in different subgroups exhibit different expression profiles as the cells undergo de-differentiation. From left to right represents the dynamic changes in gene expression levels from the starting point to the endpoint. (C–E) Markers of vascular endothelial cells (ECs) (Endoglin (ENG), Alpha-2-Macroglobulin (A2M), Platelet Endothelial Cell Adhesion Molecule 1 (PECAM1)) increase during the retro-differentiation process. (C) Each dot represents a single cell, with the colors indicating different cell clusters. From left to right indicates the cell differentiation trajectory. (D) Colored dots indicate different tissue clusters. The vertical axis represents the expression level, with each point in the figure representing a single cell, and showing the overall expression level of the gene across these cells. (E) t-Distributed Stochastic Neighbor Embedding (t-SNE) plot showing gene expression data for the marker gene. (F) Pathological tissue at the cancer-peritumoral junction was collected and subjected to immunofluorescent staining, revealing a dense distribution of tissue cells at this junction. Scale bars, 25 μm. The vascular EC marker, CD31, was highly expressed in tumor cells at the tumor junction. The entire trajectory is shown, from malignant cells in tumor tissue to differentiated ECs in adjacent normal tissue. The red line represents the boundary between tumor tissue and adjacent tissue.

To investigate tumor cell transformation into ECs during invasion, Huh7 liver cancer cells were genetically modified to express EGFP via stable infection with an EGFP-shRNA lentivirus (Fig. 2A). A total of 5 $\times$ 10⁶ EGFP-tagged cells were injected subcutaneously into nude mice housed in an IVC system for two weeks. Following tumor formation, ultrasonographic evaluation using the S-Sharp Prospect small animal imaging system revealed significant microvascular structures within the tumors (Fig. 2B). To validate these findings, tumors were pre-injected with 1% Evans blue dye to enhance blood vessel visualization, and then examined with an in vivo laser confocal microscope. This confirmed the presence of enriched vascularization (Fig. 2C). The above procedures also confirmed the successful establishment of a subcutaneous tumor model in nude mice, as shown in Fig. 2D. To determine the origins of the new vascular structures, dual immunofluorescence staining was targeted to EGFP and CD31. This showed that tumor cells mainly expressed EGFP and AFP, while some ECs co-expressed EGFP, AFP and CD31, thereby confirming the endothelial nature of some of the tumor-derived cells (Fig. 2E,F). Furthermore, a comparative analysis showed that Huh7 cells expressed significantly lower levels of EC biomarkers compared to Huh7 cell subcutaneous tumor tissues, as depicted in Fig. 2G. This differential expression suggests that the ECs in the subcutaneous tumor tissues may originate from transformed tumor cells, supporting the hypothesis that tumor cells can differentiate into ECs under certain conditions.

Fig. 2.

Transformation of liver cancer cells into neovascular ECs in a nude mouse model of cancer. (A) The Huh7 cell line was stably infected with exogenous EGFP gene. Scale bars, 25 µm. (B) The S-Sharp Prospect small animal ultrasound imaging system was used to examine the distribution of blood vessels (box). (C) After injecting 1% Evans blue into the tail vein, the vascular morphology of tumor-bearing tissue was observed with a probe for small animals using in vivo laser confocal microscopy. Scale bars, 20 µm. (D) Subcutaneous tumor-bearing tissues were successfully generated in nude mice. (E,F) A representative image from confocal microscopy. Case-1 and case-2 in (E) show that vascular ECs in tumor-bearing tissues were derived from Huh7 tumor cells, especially ECs located in the vascular cavity and which expressed both the EC marker CD31 and the tumor marker enhanced green fluorescent protein (EGFP). Similarly, the two liver cancer cases from patients shown in (F) also show that vascular ECs express both the liver cancer marker Alpha-fetoprotein (AFP) and the EC marker CD31. 4${}^{\prime }$,6-diamidino-2-phenylindole (DAPI) is the nuclear marker, and the merged panel represents the incorporated image. Arrows indicate the vascular endothelial cells within the tumors. Scale bars, 25 µm. (G) Western blot analysis showing differences in expression of the EC markers ENG and CD31 between Huh7 cells (left, case-1 and case-2) and tumor-bearing tissue (right, case-1 and case-2). Results are presented as mean $\pm$ SEM, n = 3. **p 0.01, ***p 0.001.

We next incorporated transcriptomic data from 33 distinct cancer types, supplemented with 9 datasets from the GEO. Analysis of these datasets revealed specific expression patterns for ENG, FLT1, A2M, and PECAM1. Of note, these genes were found to be significantly overexpressed in 12, 15, 15, and 17 of the 33 surveyed tumor types (Fig. 3A and Supplementary Fig. 2A–C). In-depth analysis was carried out of the specific GEO datasets targeted at liver cancer, including GSE102097, GSE22058, GSE25097, GSE12128, GSE36376, GSE46444, GSE76297, GSE57957, and GSE76427. This revealed notable expression of ENG, FLT1, and A2M in 7, 4, and 9 datasets respectively, as shown in Fig. 3B and Supplementary Fig. 2D,E. The ENG mRNA level in liver cancer samples was evaluated using qRT-PCR. This method confirmed that the ENG mRNA level was substantially upregulated in non-tumor tissue adjacent to the primary tumor, relative to the tumor tissue (Fig. 3C). To further validate these findings, Western blot analysis was performed on 20 pairs of human liver cancer samples, thereby allowing correlation of ENG protein and mRNA levels. A consistent pattern was found between ENG protein and mRNA expression, thus confirming the transcription data (Fig. 3D). This study found that vascular EC biomarkers are highly expressed in the tissue adjacent to liver tumors, suggesting their involvement in the surrounding pathological landscape of liver cancer.

Fig. 3.

The vascular EC biomarker ENG is highly expressed in adjacent tissues. (A) Differential expression of ENG in different cancer types from the The Cancer Genome Atlas (TCGA) database. Green represents tumor tissue, while red represents non-tumor tissue. (B) Differences in the expression of ENG between tumor (green) and non-tumor (red) tissue in 9 different GSE databases. (C) The mRNA level of ENG in 76 pairs of liver cancer tissue samples (tumor and adjacent normal tissue) was analyzed by quantitative real-time PCR (qRT-PCR). (D) Western blot analysis showing ENG protein expression in tumor tissue (T) and paired adjacent tissue (A) from 20 liver cancer patients. Results are presented as mean $\pm$ SEM, n = 3. Significant differences compared to the respective adjacent tissue group are indicated as: *p 0.05, **p 0.01, ***p 0.001, ****p 0.0001. TCGA, The Cancer Genome Atlas; 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; HNSC, 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, 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; HCC, hepatocellular carcinoma.

We next analyzed correlations between the expression of vascular EC markers and various clinicopathological parameters, including microvascular invasion, tumor stage and grade, age, gender, Body Mass Index (BMI), AFP levels, and Hepatitis B Virus (HBV) infection status. This analysis found no significant correlations between ENG expression and gender, age, BMI, or HBV status (Fig. 4). However, strong positive correlations were found between ENG expression and the severity of microvascular invasion and higher tumor grades. Additionally, significant positive correlations were found between A2M gene expression and both patient age (p = 0.028) and AFP level (Supplementary Fig. 3, p = 0.036). Moreover, PECAM1 expression was significantly associated with tumor grade (Supplementary Fig. 4, p 0.05). ENG was evaluated as a biomarker for microvascular infiltration in liver cancer through multivariate analysis. This identified ENG expression, microvascular invasion, AFP level, viral infection, BMI, age, gender, tumor grade, and stage as independent predictors of various survival metrics. Detailed analyses found that ENG expression was significantly associated with the outcome of liver cancer patients, including with overall survival (OS; HR = 0.795, p = 0.007), disease-specific survival (DSS; HR = 0.735, p = 0.005), disease-free interval (DFI; HR = 0.7618, p 0.001), and progression-free interval (PFI; HR = 0.7922, p = 0.001), as shown in Fig. 5A–D. Lower ENG expression in liver tumor tissue was associated with poorer OS (p = 0.00059), reduced DSS (p = 0.0019), and shorter PFI (p = 0.0014) (Fig. 5E–H). Furthermore, elevated ENG expression in peritumoral areas was consistently associated with negative outcomes, such as worse OS (p = 0.0017) and DSS (p = 0.046) (Supplementary Fig. 5). These findings indicate that ENG is an independent prognostic marker and is significantly associated with liver cancer progression and patient outcomes.

Fig. 4.

The expression of ENG in hepatocellular carcinoma is associated with clinical characteristics. (A–H) Correlation of ENG expression in liver cancer with clinicopathological characteristics such as AFP, age, Body Mass Index (BMI), gender, Hepatitis B Virus (HBV), microvascular infiltration, pathological grade, and pathological stage of tumor specimens in the TCGA database. p 0.05 indicates statistical significance.

Fig. 5.

Elevated expression of ENG in adjacent tissue promotes the progression of liver cancer. (A–D) Forrest plot of stepwise Cox univariate and multivariate proportional hazard regression models used to analyze the association between ENG expression and the outcome of patients with liver cancer. (A) overall survival (OS). (B) disease-specific survival (DSS). (C) disease-free interval (DFI). (D) progression-free interval (PFI). AFP, $\alpha$-fetoprotein; BMI, body mass index. (E–H) Analysis of the association between ENG expression in adjacent normal tissue and (E) OS, (F) DSS, (G) DFI, and (H) PFI among patients with liver cancer in the TCGA cohort. p 0.05 indicates statistical significance.

This investigation initially used the Gene Expression Profiling Interactive Analysis (GEPIA) database to analyze the relationship between ENG and the well-established vascular marker CD34. Analysis confirmed a strong positive correlation between the expression of ENG and CD34, suggesting a potential interactive role in tumor vascularization (Fig. 6A). To further substantiate this finding, the microvessel density (MVD) was assessed in liver tumor samples from a cohort of 65 patients. IHC staining for CD34 was used to quantify the vascular structures. Samples with high ENG expression exhibited a significantly higher MVD compared to samples with low ENG levels (Fig. 6B–D). Lentiviral vectors were used to introduce ENG shRNA into the liver cancer cell lines Huh7 and HepG2, with an empty vector serving as a control. Following 72 h incubation, protein extraction and subsequent Western blot analysis revealed a marked decrease in ENG expression in cells treated with ENG shRNA compared to controls (n = 3) (Fig. 6E). Transwell migration and HUVEC tubule formation assays were conducted to assess the impact of ENG downregulation. A significant reduction in tubular structure formation (n = 3) (Fig. 6F) was observed following the suppression of ENG, as well as a significant decrease in migration (n = 3) (Fig. 6G). Together, these results suggest that ENG plays a role in the angiogenesis of liver cancer.

Fig. 6.

ENG expression can promote microvessel formation in liver cancer. (A) Analysis of the GEPIA database showed that ENG and CD34 expression were positively correlated. (B) Immunohistochemistry (IHC) revealed significantly elevated expression of CD34 in tumor tissues with high expression of ENG. Scale bars, 100 µm. (C,D) The MVD in the group with high ENG expression was significantly higher than in the group with low ENG expression. (E) Western blot assay showed that the protein level of ENG decreased significantly after downregulating ENG. (F) The tubule formation assay with HUVECs revealed significantly reduced ability for tubule formation after interference of ENG expression. Scale bars, 100 µm. (G) Transwell migration assay indicated that the migration ability of HUVECs cultured with conditioned medium was significantly reduced following downregulation of ENG expression. Scale bars, 100 µm. Results are presented as mean $\pm$ SEM, n = 3. TPM, Transcripts Per Million; ns, not significant; ***p 0.001, ****p 0.0001. NC, negative control.

Initially, we conducted an enrichment analysis that focused on pathways involved in different cancer types. This analysis consistently indicated that angiogenesis is a prevalent pathway in almost all the cancer types examined (Fig. 7A,B). We next examined the direct relationship between ENG and angiogenesis. A strong positive correlation was found between ENG expression and angiogenic activity in tumors (Fig. 7C). Given the crucial role of COL1A1 in the development and progression of various cancer types [21, 22, 23], the study was expanded to examine interactions between ENG and COL1A1. Analysis of extensive datasets from TCGA, LIHC, Genotype-Tissue Expression (GTEx), and Cancer Cell Line Encyclopedia (CCLE) (https://portals.broadinstitute.org/ccle) revealed a strong positive correlation between the levels of ENG and COL1A1 expression (Fig. 7D–G). To empirically verify these correlations, Western blot assays were performed on cell lines that express low levels of ENG. This revealed that downregulation of ENG led to a marked decrease in COL1A1 expression in both Huh7 and HepG2 cells, indicating a potential regulatory relationship (n = 3) (Fig. 7H, Supplementary Fig. 10). Based on these insights, we hypothesized that ENG may regulate HUVEC migration and tubule formation by influencing COL1A1 expression. To test this hypothesis, the Huh7 and HepG2 cell lines, which typically show low ENG expression, were genetically modified to overexpress COL1A1. Following introduction of the COL1A1 overexpression plasmid and a 72 h incubation period, the results showed that overexpression of COL1A1 reversed the significant reduction in COL1A1 expression caused by ENG downregulation (Fig. 7I). Subsequently, Transwell migration and tubule formation assays were performed on the cells. The results showed significant restoration of migration ability and tubule formation in HUVECs, with a reversal of the deficits observed in the low ENG expression group (n = 3) (Fig. 7J–M).

Fig. 7.

ENG promotes microvessel formation in liver cancer by upregulating Collagen type I$\alpha$1 (COL1A1) expression. (A) Pathway enrichment for ENG in different tumor types. (B) GSEA results showing the enriched pathways in activated angiogenesis. (C) Pearson’s correlation analysis showed that ENG expression was positively correlated with angiogenesis. (D–G) The correlation between ENG and COL1A1 expression in TCGA, Liver Hepatocellular Carcinoma (LIHA), GTEx, and Cancer Cell Line Encyclopedia (CCLE) databases. Note that each dot corresponds to one cancer type (D), one patient (E), one tissue type (F), and one cell line (G). (H) Western blot assay shows significantly decreased COL1A1 protein expression compared to the control group following interference with ENG expression. (I) Western blot analysis showed that overexpression of COL1A1 reversed the significant reduction in COL1A1 expression caused by ENG downregulation. (J) Transwell cell migration assay results showing that HUVEC migration ability was significantly reduced following the infection of Huh7 and HepG2 cells with shENG lentivirus. In contrast, the migration ability of HUVECs was significantly enhanced after transfection with the COL1A1 overexpression plasmid. Scale bars, 100 µm. (K) Three fields were randomly selected under the microscope for cell counting and statistical analysis. (L) The tubule formation ability of HUVECs was significantly inhibited after Huh7 and HepG2 cells were infected with the shENG lentivirus. In contrast, the number of tubules formed by HUVECs increased significantly after transfection with the COL1A1 overexpression plasmid. Scale bars, 100 µm. (M) Three fields were randomly selected under the microscope for counting and statistical analysis in the tubule formation assay with HUVECs. Results are presented as mean $\pm$ SEM, n = 3. *p 0.05, **p 0.01, ***p 0.001. NES, Normalized Enrichment Score; PANCOR, Pan-Cancer (multiple cancer types); UVM, Uveal Melanoma; MYC, Myelocytomatosis Oncogene; E2F_TARGETS, E2F Transcription Factor Targets; G2M_CHECKPOINT, G2M Cell Cycle Checkpoint; OXIDATIVE_PHOSPHORYLATION, Oxidative Phosphorylation Pathway; KRAS, Kirsten Rat Sarcoma Viral Oncogene Homolog; UV_RESPONSE_UP, Upregulated Response to Ultraviolet (UV) Radiation; TGF_BETA_SIGNALING, Transforming Growth Factor Beta Signaling Pathway.

Bioinformatics analysis indicated that ENG potentially regulates immune cell infiltration and the extracellular matrix (ECM) in liver cancer, suggesting it is significantly involved (Supplementary Fig. 6). Subsequent investigations using the ESTIMATE algorithm assessed the immune, stromal and estimated score, as well as overall tumor purity (Supplementary Table 4). These results showed that elevated levels of ENG expression were directly associated with higher immune, stromal and estimate scores, but inversely related to tumor purity (Supplementary Fig. 7A–D). The CIBERSORT algorithm was subsequently utilized to further characterize the immune cell composition in patients with varying levels of ENG expression (Supplementary Table 5). This profiling identified substantial infiltration of various immune cell types in patients with elevated ENG levels, including central memory CD4⁺ T cells, memory CD8⁺ T cells, natural killer cells, and natural killer T cells (Fig. 8A,B). To substantiate these findings, Spearman’s correlation analysis was employed to examine correlations between ENG expression and the prevalence of various immune cell types. This revealed significant positive correlations with Type 1 T helper cells (r = 0.55), central memory CD4 T cells (r = 0.57), activated CD8 T cells (r = 0.2), Regulatory T cell (r = 0.51), T follicular helper cell (r = 0.44), Central memory CD8 T cell (r = 0.29), Gamma delta T cell (r = 0.21), Effector memeory CD8 T cell (r = 0.5), Effector memeory CD4 T cell (r = 0.31), CD56dim natural killer cell (r = 0.35), Natural killer T cell (r = 0.36), Natural killer cell (r = 0.65), Plasmacytoid dendritic cell (r = 0.48), Immature dendritic cell (r = 0.33), Activated dendritic cell (r = 0.23), and Macrophage (r = 0.55) (Fig. 8C–R). To further confirm the association of ENG with immune cell dynamics, IHC was performed on paraffin-embedded tissue sections of liver cancer. Patients with high ENG expression exhibited significantly increased levels of immune markers such as CD3, CD8 and CD45 compared to those with low ENG expression, thus substantiating the bioinformatic predictions (Supplementary Fig. 8). These findings validate the role of ENG in tumor angiogenesis and immune cell infiltration, as well as highlighting its potential as a therapeutic target.

Fig. 8.

Correlation of ENG expression with specific tumor-infiltrating immune cell types. (A,B) The heatmap represents tumor-infiltrating immune cells, memory CD4 cells, memory CD8 cells, natural killer cells, and natural killer T cells in patients with high expression of ENG. (C–R) Spearman’s correlation analysis reveals the correlations between ENG expression and the level of infiltrating immune cells using quantitative data. MDSC, Myeloid-Derived Suppressor Cells.

Detailed analyses involving data extracted from the TCGA, LIHA, GTEx, and CCLE databases were performed to explore interactions between ENG and several crucial target genes for sorafenib involved in angiogenesis, including VEGFR, PDGFR, KIT Proto-Oncogene Receptor Tyrosine Kinase (KIT), and Vascular Endothelial Growth Factor A (VEGFA). Consistent and significant positive correlations were found between ENG expression and these key angiogenic indicators (Supplementary Fig. 9). A cohort of 32 liver cancer patients undergoing sorafenib treatment was used to further investigate the clinical significance of these findings, consisting of 19 responders and 13 non-responders. Analysis revealed that non-responders typically expressed lower levels of ENG (Fig. 9D). Patients with higher levels of ENG expression generally experienced a reduction in tumor diameter after receiving sorafenib treatment (Fig. 9E). Furthermore, the level of ENG expression was negatively correlated with changes in tumor diameter (Fig. 9F). Two cases showing changes in tumor diameter following treatment with sorafenib are presented in Fig. 9A, while four cases with different levels of ENG expression are presented in Fig. 9B,C. These findings indicate the efficacy of sorafenib for suppression of tumor growth may be influenced by the expression level of ENG.

Fig. 9.

The expression level of ENG correlated with the efficacy of sorafenib in patients with liver cancer. (A) Tumor diameter was recorded by a radiologist based on the CT imaging with a black line. (B,C) Representative images of IHC staining of ENG expression in tumors from liver cancer patients. Scale bars, 100 µm, n = 32. (D) The expression level of ENG in liver cancer patients who were responders or non-responders to sorafenib, n = 32. (E) Changes in tumor diameter (mm) in liver cancer patients treated with sorafenib. Increased tumor diameter is shown in red, and decreased tumor diameter in blue. (F) Spearman’s rank correlation was performed to analyze the correlation between changes in tumor diameter and the expression level of ENG. *p 0.05, ****p 0.0001.