Combined molecular analyses identify a new type of liver cancer with important biological and clinical features

Six HCC samples were obtained from six patients who underwent hepatectomy as the initial treatment, along with one normal liver sample provided by a hepatic hemangioma patient through surgical resection, at the Cancer Hospital of Guangxi Medical University. These samples were utilized for proteomic analysis, scRNA-seq, and ST-seq. The patients were enrolled at the Cancer Hospital of Guangxi Medical University from June to September 2021. Detailed information on the diagnostic criteria for HCC, along with patient inclusion and exclusion criteria, has been reported previously (19). In summary, all enrolled patients with HCC were newly diagnosed, pathologically confirmed, and free from other cancer types. Additionally, tumor and adjacent tissues were collected from 40 HCC patients who were diagnosed and treated with radical surgery between January 2021 and January 2024 at the Cancer Hospital of Guangxi Medical University for immunohistochemical (IHC) experiments. The detailed clinical information is present in Supplementary Table S1.

We screened the HCCDB database (http://lifeome.net/database/hccdb/home.html↗) (20) to identify the candidate datasets. The inclusion criteria were as follows: 1) the dataset included both gene expression profiles and the prognosis of patients with HCC, 2) the number of patients with a survival of more than 30 days should be more than 100, and 3) the gene expression profile of the dataset should contain more than 10,000 genes. In the HCCDB database, four datasets met the above criteria. We selected and downloaded the three largest datasets by sample size (GSE14520_GPL3921, TCGA-LIHC, and LIRI-JP) for analysis. The dataset GSE14520_GPL3921 (21) containing 225 HCC and 220 tumor-adjacent liver tissue samples was utilized to develop our subtyping systems. The TCGA-LIHC dataset, which contains RNA-seq data and clinical information for 356 HCC patients from The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/tcga↗), and the LIRI-JP data set containing RNA-seq data and clinical information for 212 HCC patients from the JP Project from the International Cancer Genome Consortium (https://dcc.icgc.org/↗), were used to validate the subtyping systems.

Finally, we gathered proteomic analysis data of tumor and tumor adjacent tissues from 159 cases of HBV-associated HCC reported in Gao's study (22). This served as a validation of the proteomic analysis to confirm our findings. We also collected the ST-seq cohort GSE238264 (23), which was diagnosed with HCC and treated with a combination of tyrosine kinase inhibitors (TKIs) and programmed cell death protein 1 (PD-1) inhibitors, serving as a validation cohort sourced from the GEO database. The workflow of the present study is illustrated in Figure 1.

The gene expression profiles of GSE14520_GPL3921 were first utilized to calculate TP via the ESTIMATE package (24). GSE14520_GPL3921 was also used to calculate the TME score via the xCell tool (https://xcell.ucsf.edu/↗) (25) with the xCell gene signature. The DEGs in HCC compared to tumor-adjacent liver tissue were identified via the R package limma (v3.54.2) (26). Genes with fold changes > 1.5 and P (adjusted by false discovery rate) values < 0.05 were considered significant.

The TP and TME scores were separately analyzed via the Shapiro-Wilk test. Spearman or Pearson correlation analyses were performed to calculate the correlation between DEGs and the TP and the TME scores. A DEG that was positively correlated with TP and negatively correlated with the TME score was considered a TP-related gene, whereas a DEG that was negatively correlated with TP and positively correlated the TME score was considered a TME-related gene. In addition, the TME-related genes do not include marker genes of TME cells in the xCell signature.

The online Metascape tool (https://metascape.org/↗) was used for functional enrichment analysis of DEGs in HCC and tumor-adjacent liver tissues. We performed GSEA (27, 28) on the GSE14520_GPL3921 dataset via the GSEA Java software (http://www.gsea-msigdb.org/gsea/index.jsp↗). Hallmark and canonical pathway gene sets derived from the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database were downloaded from the Molecular Signatures Database (MSigDB) (28, 29) and used as reference gene sets. The threshold was set to a nominal P (NOM P) value < 0.05 and FDR q value < 0.25. Functional enrichment analyses of DEGs between single-cell clusters were performed via the R package clusterProfiler (v.4.6.2) (29), which is based on the Gene Ontology (GO) or MsigDB. The corrected enrichment terms with P<0.05 were considered statistically significant. Functional scoring was performed via AddModuleScore in the R package Seurat for gene sets in the MsigDB. Functional enrichment analyses of DEGs among single-cell clusters were conducted via the R package clusterProfiler (version 4.6.2) (30), which utilizes the GO framework or data from MsigDB. Enrichment terms with a corrected p-value of less than 0.05 were deemed statistically significant. Additionally, functional scoring was performed via the AddModuleScore function in the R package Seurat for gene sets derived from MsigDB.

The PPI networks of TP-related and TME-related genes were obtained from the STRING database (v11.5) (31) to preliminarily reveal the crosstalk between tumor cells and TME. The interactions with high confidence (>0.7) were included in the present study and visualized via Cytoscape software (v 3.8.0) (32).

First, to develop the TP-related polygenic risk score (PRS), overall survival (OS)-associated TP-related genes were identified via univariate Cox regression analysis. Second, the expression profiles of the OS-associated TP-related genes were used to carry out least absolute shrinkage and selection operator (LASSO) Cox regression model analysis with leave-one-out cross validation via the glmnet package (33). The genes with nonzero coefficients were considered the optimal features and subjected to multivariate Cox regression and stepwise regression analysis. The TP-related PRS was subsequently developed via the following formula: TP-related PRS = Σ (Expression_i ∗ Coeffient_i) where "Coeffient" and "Expression" represent the risk coefficient and expression of each gene in the multivariate Cox regression and stepwise regression analysis, respectively. The TME-related PRS was also developed according to the same method as above.

The optimal cutoffs of the TP-related and TME-related PRS were identified vis the surv_cutpoint function from the Survminer package (https://CRAN.R-project.org/package=survminer↗) to separately divide patients into high and low TP- and TME-related PRS groups. Each individual received a TP- and a TME-related PRS levels, and we developed the TP-TME subtype according to the TP- and TME-related PRS levels. Patients with high TP- and TME-related PRS were considered the high-risk subtype, those possessing low TP- and TME-related PRS were considered the low-risk subtype, and the remaining patients with high TP-related and low TME-related PRS or a low TP-related and high TME-related PRS were considered the intermediate-risk subtype.

The protein samples were extracted, digested, and labeled with Tandem Mass Tag (TMT) according to the experimental specifications. A 10 μL aliquot of the supernatant was injected into a nanoflow HPLC system (Thermo Scientific) linked to an Orbitrap Fusion Lumos mass spectrometer (Thermo Scientific). The extracts were then applied to an Acclaim PepMap100 C18 column and separated on an EASY-Spray C18 column. In the Orbitrap mode, the mass spectrometer performed a comprehensive mass spectrometry (MS) scan across the 300-1500 m/z range in positive ion mode (with a source voltage fixed at 2.1 kV) and achieved a resolution of 120,000. After the complete MS scan, the 20 most abundant ions with different charge states were selected for high-energy collisional dissociation fragmentation analysis. For this experiment, the UniProt HUMAN database, which was downloaded on April 20, 2019, served as the database. MS/MS data were analyzed via Proteome Discoverer 1.4.

Raw FASTQ data were processed via Cell Ranger (version 8.0.0; 10× Genomics, USA), generating gene count matrices on the basis of the human genome reference set GRCh38 with all default parameter settings. The output filtered gene expression matrix for each sample was analyzed via the R package Seurat (v.4.3.0) (34). We calculated the doublet fraction for each cell via the R package DoubletFinder (v.2.0.4) (35) with aim of removing potential doublets with a target value of 7.6% per 1000 droplets. Additionally, for each sample, cells with fewer than 300 unique molecular identifiers (UMIs), or expressing more than 7000 or fewer than 300 genes were excluded. To eliminate dead or dying cells, we further removed cells with more than 10% UMIs originating from the mitochondrial genome. Next, the "FindVariableFeatures" function in Seurat and the vst method were employed to screen for the 2000 variable genes exhibiting the largest normalized variance, which were subsequently processed for principal component analysis (PCA). The "RunHarmony" function from the Harmony package (v.1.2.0) (36) was then utilized for sample batch correction. The "RunUMAP" function was applied to perform UMAP downscaling of the first 20 principal components according to the aforementioned steps. Finally, the "FindNeighbors" and "FindClusters" functions were employed to identify cell clusters.

We collected typical marker genes for the identification of the major cell types. Epithelial cells were identified by the expression of ALB, APOA2, and EPCAM; fibroblasts were characterized by COL1A1, COL1A2, and DCN; endothelial cells were identified via PECAM1, VWF, and PLVAP; T/NK cells were marked via CD3D, CD3E, and CD2; B cells were identified via CD79A and MZB1; monocytes/macrophages were characterized via LYZ, CD86, and C1QC; mast cells were identified via TPSB2, CPA3, and KIT; neutrophils were marked via G0S2, CXCL8, and CSF3R; and dividing cells were identified via MKI67 and STMN1.

Following tissue collection, samples from patients were fixed in 10% formaldehyde for 12 hours. The tissues were subsequently dehydrated, clarified, paraffin embedded, and sectioned at a thickness of 4 µm for both H&E and IHC staining. The staining procedures were performed according to the manufacturer's guidelines. Sections stained with H&E were examined under a light microscope (OLYMPUS BX43) via a ×10 eyepiece and a ×40 objective lens, with images captured using ImageView 4.15 (Pooher) software. For immunohistochemical analysis, the sections were scanned with a Pannoramic digital section scanner (3DHISTECH). Two pathologists, blinded to the clinical characteristics and findings of the patients, independently evaluated all the sections. Scoring was conducted on a 4-point scale on the basis of the intensity of cellular staining: no positive staining (negative) received a score of 0, yellowish (weakly positive) scored a score of 1, tan (positive) scored a score of 2, and tan (strongly positive) scored a score of 3. Additionally, a 4-point scale based on the percentage of positive cells was employed, with ≤25% scoring 1, 26%-50% scoring 2, 51%-75% scoring 3, and >75% scoring 4. The final score was derived by multiplying the two individual scores. The following primary antibodies were used to bind specific IHC proteins: XPO1 (Proteintech, 27917-1-AP) and RCN2 (Proteintech, 10193-2-AP), and the secondary antibody used was horse anti-mouse/rabbit IgG (Vector, ZF1028). Raw data for 40 HCC patients were uploaded (). 2

The R package InferCNV (v1.21.0) (https://github.com/broadinstitute/inferCNV↗) was used to infer CNV changes in the scRNA-seq data. The raw gene expression counts extracted from the Seurat object were imported into the Infercnv object via the "InfercnvObject()" function. T/NK cells and B cells were selected as control datasets for reference. The CNV value for each cell was estimated via the "run()" function in InferCNV with a cutoff value of 0.1.

Activity analyses of TFs were performed to pinpoint key regulatory TFs in the selected cell clusters, and SCENIC analysis was conducted via the pySCENIC (37) package. The necessary databases for executing SCENIC, which include the TF database (cisTarget.hg38.mc9nr.feather) and the subject annotation database (hgnc.v9.1.0), were acquired from the pySCENIC website (https://github.com/aertslab/pySCENIC↗). The normalized expression matrix generated by Seurat served as the input matrix for the pySCENIC. TF activity was determined by the area under the recovery curve (AUC), and detailed findings from the transcription factor activity analyses are shown in Supplementary Table S2.

Cell-cell communication was inferred via the R software version of CellChat (v1.6.1) (38) and an existing database of ligand-receptor interactions. The apparent overexpression of ligands and receptors in specific cell clusters was initially identified via CellChat. The probability of communication occurring between two interacting clusters was quantified on the basis of the average expression level of the ligand in one cluster and the average expression level of the receptor in the other. The significance of communication was assessed via a permutation test. Interaction pairs with a P value < 0.05 were selected to assess intercellular communication. The detailed results of the cell-cell communication analysis are listed in Supplementary Table S3.

The expression matrices from the ST-seq data were processed via Seurat. "SCTransform()" normalised the values at each point, and "RunPCA()" retained the first 20 principal components (PCs) to reduce dimensionality. The top 30 DEGs for each cell cluster in the scRNA-seq cohort were used as input. Scores were assigned to individual points in the ST-seq dataset via the "gsva()" function of the R package GSVA (v1.46) (39) with default parameters. The spatial feature expression plots were generated via the "SpatialFeaturePlot()" function in Seurat.

Gene mutation data from the TCGA-LIHC dataset were extracted from mutation annotation format (MAF) files via the GDCquery_Maf function in the "TCGAbiolinks" package (40). The gene mutation frequencies of each risk subtype were visualized as a waterfall plot via the oncoplot function in the "TCGAbiolinks" package. The tumor mutational burden (TMB) of each sample was obtained from a previous study (41). The stemness score (42) was calculated for each individual in the TCGA-LIHC dataset via the TCGA analyze_Stemness function in the "TCGAbiolinks" package.

For the TP-TME risk subtypes, we compared the expression levels of two immune checkpoints (PDL1 and CTLA4) and five antigens (CD133, EPCAM, GCP3, MSLN, and MUC1) (43) to predict the potential response to these treatments. In addition, we performed drug repositioning analysis for the high-risk subtype via the PATHOME-Drug (http://statgen.snu.ac.kr/software/pathome/↗) web tool.

Unless otherwise stated, all analyses were performed in R (version 4.2.3). We identified DEGs via unpaired t-tests provided via the limma package. The Shapiro-Wilk test was used for the normality test. Time-dependent receiver operating characteristic (tROC) curve analysis was carried out using the tROC package (44). Kaplan–Meier survival curves for OS and progression free survival (PFS) were compared in different subtypes using the log-rank method in the "survival" package (https://CRAN.R-project.org/package=survival↗) and the "survminer" package (https://CRAN.R-project.org/package=survminer↗). Intergroup differences in continuous variables were assessed for significance via Wilcoxon, Kruskal–Wallis, or unpaired t tests. All tests were two-sided, and unless otherwise stated, we set P value < 0.05 to indicate statistical significance. Visualization was done vis the "ggplot2" R package.

Using the GSE14520_GPL3921 dataset, we identified a total of 2,263 differentially expressed genes (DEGs) in HCC compared with tumor-adjacent liver tissue (Figure 2A). A total of 451 TP-related genes and 121 TME-related genes were identified through Spearman correlation analysis. Bidirectional hierarchical clustering revealed that the expression patterns of these genes could largely distinguish between HCC and tumor-adjacent liver tissue (Figure 2B). Unsurprisingly, the TP-related genes were enriched in mainly cancer-related Gene Ontology (GO) terms and pathways, such as the cell cycle and mismatch repair (Figure 2C). In contrast, TME-related genes were associated primarily with immune system processes (Figure 2D). The PPI networks of the TP-related and TME-related genes contain 342 nodes and 1177 edges (Supplementary Figure S1). In the PPI networks, red nodes indicate genes whose expression is upregulated, and blue nodes indicate genes whose expression is downregulated in tumors. Circular nodes represent TP-related genes, whereas rhombic nodes represent TME-related genes.

Fifty TP-related genes were identified as OS-associated genes, twenty-two of which presented nonzero coefficients (Figure 3A), and 11 genes (ALG6, ATP5MF, CNIH4, ESM1, HEY1, LANCL1, P2RX4, PEX11B, POP7, RCN2, and XPO1) were used to generate the TP-related PRS (Supplementary Table S4). The TP-related PRS was significantly associated with OS [P < 0.0001, hazard ratio (HR) = 2.718 (95% CI for HR = 2.147-3.442)], and the area under the curve (AUC) of the tROC analysis was stable at approximately 0.8 (Figure 3B). The HCC patients with high TP-related PRS had shorter OS than did those with low TP-related PRS (P < 0.0001) (Figure 3C). Among the TME-related genes, twelve genes were identified as OS-associated genes, ten of which had nonzero coefficients (Figure 3D), and seven genes (ALDH1B1, CTSC, GUCY1A1, MRC1, SPRY2, TARP, and TRIM22) were used to generate the TME-related PRS (Supplementary Table S5). The TME-related PRS was also significantly associated with OS [P < 0.0001, HR = 2.718 (95%CI for HR = 1.978-3.735)], and the AUC of tROC analysis was 0.7-0.8 (Figure 3E). HCC patients with a high TME-related PRS had shorter OS than did those with low TP-related PRS (P < 0.0001) (Figure 3F). Our TP-TME risk subtype was generated on the basis of two PRSs, and 34, 52, and 123 patients with HCC were divided into high-, intermediate-, and low-risk subtypes, respectively. Patients in the high-risk subtype had the poorest survival, those in the low-risk subtype had the best survival, and those in the intermediate-risk cases had a better prognosis than did those in the high-risk subtype and worse prognosis than did those in the low-risk subtype did (Figure 3G). A similar trend was also observed for PFS (Figure 3H). Furthermore, the TP-TME risk subtyping system demonstrated superior prognostic predictive power compared to routine clinicopathological features and remained independent of these features (Figure 3I). As in the GSE14520_GPL3921 dataset, the TP-related PRS and the TME-related PRS in the TCGA-LIHC and LIRI-JP datasets were calculated according to the abovementioned formulas. We found similar results in the TCGA-LIHC cohort (Supplementary Figures S2A–E) and the LIRI-JP cohort (Supplementary Figures S2F–I). Overall, TP-TME risk subtype is a robust prognostic prediction system.

Compared with those in the TP-TME intermediate- and high-risk subtypes, the liver function-related HALLMARK (Figure 4A) and metabolism-related KEGG (Figure 4B) gene sets were significantly enriched in the TP-TME low-risk subtype. These findings suggest that the TP-TME low-risk subtype of HCC is well- differentiated. The TP-TME intermediate-risk subtype was characterized by enrichment of transcription factor E2F and MYC targets (Figure 4C) and cell cycle pathways (Figure 4D), whereas the TP-TME high-risk subtype was characterized by enrichment of hypoxia, Wnt/beta-catenin signaling (Figure 4E), the Notch signaling pathway and the TGF-beta signaling pathway (Figure 4F). These results indicate that there is significant biological heterogeneity among these three subtypes.

To validate the expression of genes characterizing the TP-TME risk subtype at the protein level, we collected tumor and adjacent-tumor tissues from 6 HCC patients for proteomic analysis. Analysis of the protein expression levels of genes associated with the HCC-TP-TME risk subtype showed, that XPO1, RCN2, PEX11B, P2RX4, LANCL1, ATP5MF, ALG6, TRIM22, GUCY1A3, CTSC, and ALDH1B1 were significantly elevated in tumor tissues compared with adjacent tumor tissues (Figure 4G). Furthermore, we validated these findings via published data from Gao et al. (22), which corroborated our results in an independent cohort (Figure 4H). Collectively, these results indicate that the characteristic genes of the TP-TME risk subtypes identified in our study have important clinical implications for the prognostic assessment of HCC and warrant further investigation into the biological functions of these genes.

By analyzing the scRNA-seq data for HCC, we obtained 62,163 high-quality cells. Nine distinct cell types were identified on the basis of known markers: epithelial cells, fibroblasts, endothelial cells, T/NK cells, B cells, monocytes/macrophages, mast cells, neutrophils, and cycling cells (Figures 5A–C). Furthermore, we analyzed the proportions of various cell types across different patients and found that, although all cell types were present in each patient, the predominant infiltrating cell types varied, which may reflect heterogeneity among HCC patients (Figure 5D).

Next, we re-clustered the epithelial cells on the basis of their differentially expressed genes, identifying 4 clusters: XPO1+Epithelial, CYP2E1+Epithelial, S100A6+Epithelial, and STMN1+Epithelial (Figure 5E). The bubble diagram illustrates the highly expressed genes in each cluster (Figure 5F). Functional enrichment analysis revealed that XPO1+Epithelial was predominantly associated with the acute inflammatory response, cell growth, positive regulation of angiogenesis, and epithelial cell proliferation, indicating its potential role in tumor progression. In contrast, CYP2E1+Epithelial were primarily linked to material and energy metabolism, whereas S100A6+Epithelial were associated with the regulation of ubiquitin-protein ligase activity, among other functions. STMN1+Epithelial exhibited characteristics related to cytokinesis (Figure 5F). We subsequently analyzed the chromosome copy number variation in these epithelial cells via inferCNV. The results indicated that all four clusters presented significant chromosome amplifications and deletions (Supplementary Figure S3A). Additionally, we observed minimal capture of epithelial cells in normal liver tissue samples (Figure 5D), leading us to conclude that all four clusters consisted of malignant epithelial cells.

We found that 7 genes characterizing TME risk subtypes were expressed at varying levels across multiple cell types within the TME (Supplementary Figure S3B). However, 11 genes associated with the TP risk subtype signatures exhibited significantly increased expression in the XPO1+Epithelial. Specifically, the expression levels of XPO1, RCN2, P2RX4, PEX11B, LANCL1, HEY1, and ESM1 were significantly increase in this cluster than in the other clusters (Figures 5G, H). Consequently, we focused on this cluster characterized by high expression of genes in the XPO1+Epithelial group. Next, we collected tumor tissues and adjacent-tumor tissues from 40 HCC patients for immunohistochemical staining (Figure 5I) to validate our findings. We observed that the expression levels of XPO1 and RCN2 were elevated in tumor tissues compared with adjacent-tumor tissues (Figure 5J). Further analysis indicated that the expression of XPO1 and RCN2 was greater in tumor tissues with high Ki67 expression (Figures 5I, K, L). These results suggest that the elevated expression of XPO1 and RCN2 is closely associated with the proliferation of HCC.

Further analysis revealed that the XPO1+Epithelial cluster scored significantly increase than other clusters did in terms of functions associated with the TGF-beta signaling pathway, WNT/beta-catenin signaling, Notch signaling, Hedgehog signaling, and hypoxia signaling (Figure 6A). This finding was consistent with the upregulation of these functions observed in high-risk patients within the TP-TME risk model (Figures 4A–F). Additionally, compared with the other clusters, the XPO1+Epithelial cluster presented significantly increase scores for proliferation, migration, epithelial-mesenchymal transition (EMT), and angiogenesis (Figure 6A). Our analysis of the key TFs driving the distinct functions of these clusters revealed that the ten most active transcription factors in XPO1+Epithelial were predominantly associated with tumor proliferation, migration, and invasion (Figure 6B). For example, MEF2A may play a dual role in promoting tumor proliferation and metastasis by inducing the activation of EMT and WNT/beta -catenin signaling (45), whereas TCF7L2 serves as a core TF of the WNT signaling pathway, and is involved in regulating tumor cell proliferation and migration (46). Furthermore, we observed that patients exhibiting high expression of XPO1+Epithelial signatures had significantly shorter OS and PFS in the TCGA cohort (Figure 6C).

In summary, we found that XPO1+Epithelial constitute a cluster of cells characterized by TP high-risk subtypes, which exhibit elevated expression in cancer tissues. This cluster shows higher expression levels in patients with tumors that have high proliferation rates and is positively correlated with poor patient prognosis. The underlying mechanism may involve the contribution of this cluster to tumor progression through pathways such as EMT and WNT/beta-catenin signaling.

To further explore the crosstalk between XPO1+Epithelial and the TME in HCC, we performed intercellular communication analysis via 'CellChat'. The results indicated that among the four epithelial cell clusters, XPO1+Epithelial exhibited the highest degree of communication with TME components (Figure 7A). XPO1+Epithelial was regulated primarily by fibroblasts (Figure 7B). Additionally, signals sent by XPO1+Epithelial predominantly regulate monocyte/macrophages, endothelial cells, and T/NK cells (Figure 7C). Our further analysis of key ligand-receptor pairs that interact with XPO1+Epithelial (Figures 7D, E) revealed that fibroblasts regulate XPO1+Epithelial mainly through ligand-receptor pairs such as CD99-CD99 and FN1-(ITGAV+ITGB1) (Figure 7D). Conversely, XPO1+Epithelial promotes monocyte/macrophage recruitment by regulating monocytes through MIF-(CD74+CXCR4) and MIF-(CD74+CD44), as well as C3-(ITGAX+ITGB2). Furthermore, the iso-ligand receptors MIF-(CD74+CXCR4) and MIF-(CD74+CD44) regulate T/NK cells (Figure 7E), potentially facilitating tumor cell evasion of immune surveillance (47, 48). Additionally, XPO1+Epithelial regulates endothelial cells via VEGFA-VEGFR1R2 and VEGFA-VEGFR1, which may be associated with promoting tumor vascular production (49, 50).

To confirm the results observed in the cell-cell communication analyses, we validated our findings via the TCGA-LIHC cohort at the bulk RNA-seq level. These data corroborated the positive correlation between the expression levels of specific receptors in the target cells (Figure 7F; Supplementary Figure S4A). Further evidence was provided by spatial transcriptome analysis (Figure 7G; Supplementary Figure S4B), which consistently demonstrated that fibroblasts co-localized with XPO1+Epithelial in specific physical locations. Additionally, the scores for both cell types were positively correlated, and fibroblasts regulated the co-localization of the primary ligand receptor for XPO1+Epithelial cells. Similarly, we observed physical positional co-localization between XPO1+Epithelial and monocytes/macrophages, T/NK cells, endothelial cells, and their corresponding ligands and receptors, suggesting communication exchanges among these cell types. Thus, we validated the interaction network between XPO1+Epithelial cells and the TME via multi-omics.

In summary, our findings indicate that XPO1+Epithelial cells are key components in the remodeling of the TME, and are regulated primarily by ligand signaling from fibroblasts. This interaction may modulate endothelial cells, monocyte macrophages, and T-cells through seeded ligand receptors, potentially influencing immune cell recruitment, immunosuppression, and pro-angiogenesis.

In light of the analyses conducted at both the scRNA-seq and ST-seq levels, we determined that XPO1+Epithelial in HCC significant interacts with endothelial cells and T cells. This observation led us to speculate that XPO1+Epithelial may serve as a potential target for anti-angiogenic therapies and immunotherapy. To further investigate this hypothesis, we validated our findings via spatial transcriptome samples from the HCC treatment cohort. Compared with the nonresponsive group, the group that responded to the combination of TKIs and PD-1 treatment presented significantly increase signature gene scores (Figures 8A, B) for TP-TME risk subtypes, as well as elevated scores for XPO1+Epithelial (Figure 8C). Additionally, T-cell and endothelial cell infiltration was notably more pronounced in the combination treatment response group (Figures 8D, E).

Through PATHOME-Drug analysis, we constructed drug-target networks to identify potentially effective drugs for the high-risk subtypes (). The identified drugs included recommended agents, such as sorafenib, regorafenib, cabozantinib, and pembrolizumab (). Collectively, these results suggest that TP-TME risk subtypes may be used to predict the efficacy of targeted and immunological therapies, warranting further investigation in follow-up cohort studies. 5 1

We conducted a preliminary exploration of the potential mechanisms associated with the genes that characterize the TP-TME risk subtypes and their implications for immunotherapy efficacy. Analysis of the top 30 mutated genes (Figures 8F–H) across the high-, intermediate-, and low-risk subtypes revealed that the genes with the highest mutation percentages in the different risk groups included TP53, TTN, and CTNNB1. Notably, the mutation percentages are greater in the intermediate- and high-risk groups, and these mutations are closely linked to the onset and progression of various tumors. There are currently few known drugs that target these mutated genes. Additionally, we observed that the tumor mutational burden (TMB) was low in HCC patients and did not differ significantly among the three subtypes (Supplementary Figure S5B), indicating that TMB may not serve as a valid biomarker for selecting HCC patients for treatment with immune checkpoint inhibitors (ICIs). Furthermore, our analyses revealed no significant differences in the stemness scores among the three risk subtypes (Supplementary Figure S5C). Although no significant differences were observed in the expression of CD274 (also known as PDL1) among the three risk subtypes, the expression levels of CTLA4 and PDCD1 were significantly increase in the intermediate- and high-risk subtypes than in compared to the low-risk subtype (Figure 8I). Consequently, the response rates of the intermediate- and high-risk subtypes to ICIs treatment may be greater than that those of the low-risk subtype. Furthermore, the three risk subtypes of the TP-TME presented distinct expression levels of the 5 cancer antigens targeted in chimeric antigen receptor-modified T cell (CAR-T) therapy (Figure 8I). Specifically, GPC3 expression was elevated in the intermediate-risk subtype relative to the low-risk subtype, whereas MSLN expression was higher in both the intermediate- and high-risk subtypes than in the low-risk subtype. The highest expression levels of MUC1, EPCAM, and PROM1 were observed in the high-risk subtype. Therefore, the three TP-TME risk subtypes may exhibit varying therapeutic responses to the corresponding CAR-T therapies.