Single-cell transcriptomic data from HCC patients were accessed from the GEO database (GSE242889) and included 5 samples from patients who underwent surgical resection. Three patients had MVI and two did not. Inclusion criteria were: the absence of pre-surgical metastasis, the presence of a single tumor, 3 - 5 cm in diameter and the absence of prior treatment. Spatial transcriptomics (ST) data, including 5 samples from 2 HCC patients, one of whom had combined portal vein tumor thrombus (PVTT), were downloaded from Wu et al. (20, 21).

HCC tissue samples with (n = 2) and without MVI (n = 2) were collected from Eastern Hepatobiliary Surgery Hospital, fixed in 10% neutral-buffered formalin for 24 - 48 hours (h), embedded in paraffin and 4 µm sections made for IHC analysis using a rabbit polyclonal anti-TM4SF1 antibody (dilution 1: 200). Sections were deparaffinized, rehydrated and antigen retrieval performed in citrate buffer (potential of hydrogen 6.0) at 95°C for 20 minutes. Endogenous peroxidase activity was blocked with 3% hydrogen peroxide and sections incubated overnight at 4°C with primary antibody. Biotinylated secondary antibody and streptavidin-peroxidase conjugate was added, developed with 3,3’-diaminobenzidine (DAB) substrate and counterstained with hematoxylin. Images of five random high-power fields (HPFs) within the vascular regions were captured and staining intensity analyzed by ImageJ software with IHC Profiler macro as high positive, positive, low positive or negative to generate histoscores. Sample staining intensity was represented as the sum of high positive and positive staining. This study was approved by the Ethic Committee of the the Third Affiliated Hospital of Naval Medical University (EHBHKY2021-Y-007).

scRNA-seq data was imported via the Seurat function Read10X which removed cells with fewer than 400 unique genes or more than 7,000 genes (22). Cells were combined into a Seurat object and those with counts greater than 10,000 or more than 5% mitochondrial RNA were removed. Data processing was performed using Seurat’s standard pipeline. Data were normalized using LogNormalize with a scale factor of 10,000. Variable features were identified by the findVariableFeatures function with vst, selecting 2000 features. ScaleData was applied to all genes and principal component analysis performed using RunPCA. Cells were clustered with the FindNeighbors function, using dimensions 1 - 20, and FindClusters at a resolution of 1.0. Cell annotation was performed via singleR and tumor cells identified by inferCNV (23).

High-dimensional weighted correlation network analysis (hdWGCNA) was performed on high-dimensional transcriptomic data, including scRNA-seq and spatial transcriptomics data, using the hdWGCNA package (24), to identify interconnected genes with contextual information from biological knowledge sources. Network modules and genes with significant MVI association were indicated.

Gene Ontology (GO) enrichment analyses were performed on MVI modules using the clusterProfiler R package (25) and the top 5 GO terms identified.

Candidate genes were analyzed via Gene Expression Profiling Interactive Analysis (GEPIA), an open platform that integrates data from The Cancer Genome Atlas (TCGA) and the Genotype Tissue Expression (GTEx) databases (26). Genes with high expression in tumor tissues and impact on HCC prognosis were identified.

To explore the dynamic developmental trajectories of cell populations, we performed trajectory analysis utilizing CytoTRACE and the Monocle2 R package. By quantifying the similarity in gene expression profiles among individual cells, CytoTRACE allowed us to infer the hierarchical order of cellular differentiation. Furthermore, Monocle2 was applied to delve deeper into the trajectory patterns of specific subpopulations. These complementary analyses provided a comprehensive view of the progression of cell states and the lineage relationships (27, 28).

A total of 9,061 cells remained after quality filtration, dimensionality reduction and cluster analysis and were divided into 20 clusters. 5,360 cells were derived from MVI tissues and 3,701 from non-MVI tissues. Each cluster was assigned a single cell type by singleR: B cells, T cells, dendritic cells (DCs), endothelial cells, hepatocytes, macrophages, monocytes or smooth muscle cells ( Figure 1A ). Hepatocytes were confirmed as the cells accounting for malignancy by InferCNV analysis ( Supplementary Figure S1 ). 1,478 endothelial cells were identified and reclustered, at a resolution of 0.2, into 7 distinct clusters. The top 5 marker genes in each cluster were identified by the FindAllMarkers algorithm ( Figure 1B ) and the first two genes used to assign cellular phenotype ( Figure 1C ). Marker genes differentially expressed in individual clusters were subjected to GO analysis and TECs cluster functions inferred ( Figures 1D-K and Supplementary Tables S1-S8 ).

Functional cluster analysis produced the following results: the TECs0 cluster showed high expression of the genes, TPT1/VWF, associated with endothelial cell differentiation and angiogenesis regulation; TECs1 had high expression of CLEC1B/CLEC4G and TECs4 of MT1G/RGCC both of which were involved in stress responses to metal ions and detoxification; TECs2 had high expression of DUSP1/ZFP36, associated with epithelial cell migration and leukocyte differentiation; TECs3 had high expression of FBLN5/TMEM100, linked to the EMT; TECs5 had high expression of HLA-DRB5/RBP4, related to cytoplasmic translation and oxidative phosphorylation; TECs6 had high expression of NID2/MK167, involved in cell cycle regulation and TECs7 had high expression of LYZ/TYROBP, associated with immune response regulation.

Gene modules linked to MVI and significant module-specific genes were identified by hdWGCNA. Transcriptomic samples from MVI and non-MVI tissues were integrated and an optimal soft threshold of 4 determined from the scale-free topology model ( Figure 2A ). Gene expression patterns were used to generate a hierarchical clustering dendrogram with different colors indicating module membership. Nine co-expression modules (SM) related to scRNA-seq were identified ( Figure 2B ).

All modules apart from SM1 and SM4 were associated with MVI ( Figure 2C ) and clusters, TECs0, TECs2, TECs3 and TECs6, belonged to MVI-related modules ( Figure 2D ). Genes expressed in the SMs were assigned kME values (calculated using hdWGCNA software) to indicate the degree of connectivity within the module ( Figure 2E ). Those genes showing greater degree of connectivity were considered more central to module function and the top 10 genes for each SM are shown in Figure 2F .

Genes that were differentially expressed (DEGs) between MVI and non-MVI TECs were defined according to the criteria of logFC > 0.5, adjusted p value < 0.05 and KME > 0.7 ( Figure 3A , Supplementary Table S9 and S10 ). Kaplan-Meier (K-M) survival analysis indicated the association of high expression of TM4SF1, PRL23 and PRS12 with poor prognosis in HCC patients ( Figures 3B-D ). TM4SF1 was highly expressed in the TECs3 cluster whereas PRL23 and PRS12 were expressed at similar levels throughout the subgroups ( Figures 3E-G ). scRNA-seq analysis showed a greater proportion of TECs3 cells in MVI samples than in non-MVI ( Figures 3H-J ), implying that TECs undergoing the EMT are involved in MVI pathogenesis in HCC. hdWGCNA showed TECs3 to be associated with SM2 and SM6 ( Figure 2D ). SM2 was shown to be associated with cell-matrix adhesion and epithelial cell migration and SM6 with angiogenesis and EMT regulation by GO analysis ( Figures 3K, L , Supplementary Tables S11 and S12 ). Thus, the TECs3 subgroup of cells and their high expression of TM4SF1 may facilitate HCC cell invasion and MVI.

To understand the developmental trajectory of endothelial cells in hepatocellular carcinoma (HCC) and its relationship with microvascular invasion (MVI), we utilized CytoTRACE and Monocle2 for analysis.

CytoTRACE, which utilizes a probabilistic model, estimated the progression of endothelial cells along a trajectory and concentrated on capturing the dynamic alterations and lineage relationships within the endothelial cell population. The results of the CytoTRACE analysis are presented in Figures 4A, E . In MVI samples ( Figure 4A ), the CytoTRACE scores and phenotypes of the endothelial cells showed a particular distribution pattern, which provided initial insights into the possible developmental states of these cells. In non - MVI samples ( Figure 4E ), a different pattern was observed, suggesting differences in the developmental characteristics of endothelial cells between MVI and non - MVI conditions.

Monocle2, with its reverse graph embedding approach, arranged the cells in a sequential order along a trajectory, permitting the visualization of continuous developmental processes. Monocle2 was applied to evaluate the differentiation of endothelial cell subtypes. This enabled us to comprehend the differentiation kinetics of different endothelial cell subpopulations in HCC and infer the dynamic changes during the differentiation process. In MVI samples, as depicted in Figures 4B-D , the proposed temporal distribution of endothelial cell subpopulations was distinctly differentiated. TECs3 was situated at the branch inception for MVI samples ( Figures 4B, C ), implying its possible role in the early stages of endothelial cell development related to MVI. In contrast, for non - MVI samples ( Figures 4F, G ), TECs3 was located at the branch terminus.

Pseudotime dynamics analysis revealed that TM4SF1 was predominantly expressed in TECs3. In MVI samples, the expression of TM4SF1 declined during ontogeny ( Figure 4D ), while in non - MVI samples, it increased ( Figure 4H ). The changes in TM4SF1 expression during the developmental process of TECs3 might play a crucial role in modulating the function of endothelial cells and their interaction with tumor cells, thereby influencing the occurrence and progression of MVI in HCC.

MVI during HCC may lead to macrovascular invasion, when the tumor extends into a major vessel, and portal vein tumor thrombus (PVTT). The spatial enrichment of TM4SF1 expression and of TECs3 as a proportion of the total endothelial cell content were analyzed in control (HCCN), tumor peripheral (HCCL), central tumor (HCCT) and PVTT tumor (HCCP) tissues. TM4SF1 expression was higher in HCCT and HCCP than in HCCL. Control (HCCN) tissues showed the lowest expression of all samples ( Figures 5A, B ). TM4SF1 expression was significantly higher in HCCT with PVTT than in HCCT1 without PVTT ( Figures 5C, D ). A similar pattern of spatial TECs3 enrichment was shown to that for spatial expression of TM4SF1. TECs3s were more abundant in HCCL, HCCT and HCCP tissues than in control tissues ( Figures 5E, F ). Moreover, a greater proportion of TECs3s was found in tumor tissues with PVTT than in those without PVTT ( Figures 5G, H ). Thus, TECs3 abundance and high TM4SF1 expression were associated with the progression of MVI to more advanced stages of vascular invasion.

Quantitative analysis showed MVI samples to show a greater overall expression of TM4SF1 than non-MVI samples, producing a higher overall histoscore and indicating increased TM4SF1 expression in the presence of MVI ( Figure 6 ).