The datasets of HCC were acquired from GEO, including 4 scRNA-seq cohorts (GSE156625 (24), GSE149614 (25), GSE151530 (26), GSE189903 (27)). Specifically, GSE156625 (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA658541&o=acc_s%3Aa) was the raw sequencing data and the other 3 datasets are processed read count matrices. To accurately analyze the impact of tumor subclonal complexity on the remodeling of the tumor microenvironment, we established the following criteria for sample selection: 1) all samples are from the 10X Genomics platform; 2) availability of raw sequencing data or read count expression matrices; 3) samples are untreated; 4) samples are primary tumors without invasion, recurrence and metastasis; 5) for multi-region sequencing data, we selected samples from the tumor core region whenever possible. Finally, only 30 samples were retained for further analysis (Supplementary Table S1).

Raw reads from GSE156625 were processed to generate gene expression matrices by using the standard internal pipeline based on the Cell Ranger toolkit (3.1.1). Then the expression matrices of each individual sample were converted to Seurat object by the “CreateSeuratObject” function from Seurat package (4.3.0) (Hao et al., 2021). We firstly filtered out cells that had either lower than 200 or higher than 8,000 expressed genes. Furthermore, the cells with the percent of mitochondrial genes over 40% of total expressed genes were discarded. Additionally, we initially removed doublet cells by DoubletFinder (2.0.3) (McGinnis et al., 2019) to avoid the potential influence on analytical results. A total of 67,125 cells were retained for further analysis.

To better reflect the biological features from the data, we used the “NormalizeData” function to normalize the expression matrices for the Seurat object. And the top 2,000 highly variable genes (HVGs) identified by “FindVariableFeatures” function were used to scale data through “ScaleData” function. The principal component analysis (PCA) was carried out by “RunPCA” function to generate 50 PCs, followed by “RunHarmony” function implanted in Harmony package (0.1.1) (Korsunsky et al., 2019) with the top 30 PCs to eliminate batch effect. The clustering analysis was implemented depended on the integrated joint embedding produced by Harmony with the Louvain algorithm after computing a shared nearest-neighbor (SNN) graph. The UMAP technique was used to create a 2D map on which the discovered clusters were displayed. Then we used “FindAllMarkers” function to identify over-expressed genes in specific cluster when compared with the other clusters (adjusted P-value <0.05, only.pos = T and logfc.threshold = 0.25). The well-known cell markers (Ma et al., 2021; Liu et al., 2023) were used to annotate the clusters: Endothelial (ENG, VWF, PECAM1), Fibroblasts (ACTA2, COL1A2, PDGFRB), T cells (CD2, CD3D, CD3E, CD3G), B cells (CD79A, CD79B, MS4A1, BANK1), Myeloid (CD14, CD68), Plasma (MZB1, DERL3, SDC1), Hepatocytes (APOA2, ALB, APOA1, AMBP), proliferative T (TOP2A, MKI67, TUBB).

We further clustered T cells, endothelial cells and myeloid cells following a similar pipeline as described above, including normalization, highly variable genes identification, dimensionality reduction, batch effect correction with Harmony and clustering. Within each cell subtype, an iterative procedure that a cell cluster displaying at least two canonical marker expression of 7 major cell type was determined as doublet cells to be removed, if any, and reclustered the remaining cells, was adopted to ensure the reliability of the study. Specifically, we used the following well-known markers (Cheng et al., 2021; Zheng et al., 2021) for subtype identification. NK (GNLY, KLRF1, NKG7), CD4+_CCR7 (CD4, CCR7), CD4+_IFNG+ (CD4, IFNG), CD4+_FOXP3_Treg (CD4, FOXP3, IL2RA), CD8+_GZMK (CD8A, GZMK), CD8+_TNFSF9 (CD8A, TNFSF9), macrophage (C1QC, APOE, CD68), monocyte (FCN1, LYZ), cDC1_CLEC9A (CLEC9A, FLT3), cDC2_CD1C (CD1C, FCER1A), cDC3_LAMP3 (LAMP3, CCR7), mast (KIT). Similarly, endothelial cell were also confirmed by canonical markers (Goveia et al., 2020; Schupp et al., 2021), such as activated_PCV (ACKR1, SELP, VCAM1, POSTN), Capillaries_EDNRB (EDNRB, CA4, HPGD, IL1RL1), Capillaries_FCN3 (FCN3, BTNL9, NOSTRIN, EDN), arteries_GJA5 (GJA5, EFNB2, SOX17, DKK2), tip_CXCR4 (CXCR4, ADM, ANGPT2, APLN), RBP4+_EC (RBP4, IGFBP7), CXCL10+_EC (CXCL10, GBP1, SOD2). To accurately annotate proliferative T cells, we calculated the averaged Pearson correlation coefficient and averaged p-value for each proliferative T cell against other T cell subtypes and NK subtypes based on single-cell expression profiles. Proliferative T cells were excluded if the mean correlation coefficients were all less than 0.4 or the mean p-values all exceeded 0.05. Otherwise, they were assigned to the subtype with the highest average correlation.

We inferred copy number variations (CNVs) for each cell through the InferCNV (1.10.1) package of R. The non-malignant cells were applied as the reference to estimate the CNVs values of the hepatocytes (observation). The genes were sorted by their genomic locations on each chromosome. Specifically, The InferCNV analysis was implemented with parameters “denoise = T, HMM = T, analysis_mode = “subclusters”, cutoff = 0.1, cluster_by_groups = T”.

To identify the malignant cells from the hepatocytes, we designed the CNVs score quantifying the degree of CNVs fluctuation for each cell as follows: C N V s c o r e = ∑ i = 1 X i − X mean 2   /   n where X i represented the CNVs value of gene i in the cell, X mean was the mean CNVs value of gene i in the cell, and n was the number of the genes in the cell. If the CNVs score of a hepatocyte exceeded the threshold determined by the intersection of the distributions of CNVs score in hepatocytes and non-malignant cells, then the hepatocyte was considered malignant (Supplementary Figure S1C).

The non-negative Matrix Factorization (NMF) algorithm was applied to identify the underlying expression programs from the malignant cells based on the NMF R package. We employed NMF (rank = 2:6, nrun = 10) to the relative expression matrices of the top 2000 HVGs in each sample with all negative values converted to 0. We chose the optimal rank value at which the cophenetic coefficient started producing the maximum descent (Supplementary Figure S1E). In total, we identified 98 programs from 30 samples (Supplementary Table S2). The 30 genes of each expression program with the top NMF scores were input into “AddModuleScore” function to generate the program score for malignant cells in each sample. We calculated the correlations between the 98 program scores in each sample individually. After excluding clustering groups containing less than half of samples, we finally extracted 3 meta-programs from the 98 programs according to hierarchical clustering of averaged correlations of pairs of programs across all samples (Figure 1F; Supplementary Table S3). Each meta-program, retaining 30 genes with the highest average NMF score, was used to perform pathway enrichment analysis.

Shannon entropy was used to quantify the CNV complexity for each sample based on discrete CNVs profiles inferred from InferCNV. Firstly, we calculated the frequency of copy number gains and losses for each gene within the tumor as follows: L V F = n l o s s / n G V F = n g a i n / n where n represented the total number of malignant cells, n l o s s was the number of cells in which the copy number of the gene was lost, and n g a i n was the number of cells in which the copy number of the gene was gained. Specifically, Discrete CNVs value of the gene greater than (less than) 3 was considered to be gained (lost). The frequency of copy number gain and loss were divided into 10 equal-length frequency intervals separately (i.e., Loss:0–0.1, 0.1–0.2, …, 0.9–1; Gain:0–0.1, 0.1–0.2, …, 0.9–1). After mapping all the genes to equal-length frequency intervals according to the L V F and G V F of the genes, we then counted the number of genes belonging to each equal-length frequency interval and calculated the frequency of the number of genes in each interval. Finally, the CNV complexity score could be calculated for each sample as follow: C N V _ c o m p l e x i t y _ s c o r e = − ∑ i P X i log 2 ⁡ P X i where i was the ith interval,   X i was the number of genes falling into the ith interval, and P X i was the frequency of the number of genes in the ith interval.

To evaluate whether the CNV complexity score could accurately distinguish between polyclonal and monoclonal samples, we used computer simulations to generate single-cell CNV profiles for 100 polyclonal samples and 100 monoclonal samples as the ground truth data, respectively. Briefly, the size of cell population N ∼ U [200, 2000] and the total number of clones C ∼ U [5, 20] were generated for each sample. For each clone, the size of clone was generated through binomial distribution B (N, p). For monoclonal samples, we assumed that p = 0.9 for one dominant clone and p = 0.1/(c-1) for other rare clones. We considered int (N/2) subclones in the tumor cells, where int takes the nearest integer in the polyclonal samples. We used p = 2/N to generate the size of subclone. For detailed steps of the entire simulation data generation process, please refer to Wang et al. (2023). We then applied the three methods (CNV complexity score, Guo et al.’s and Ma et al.’s) to the simulated single-cell CNV profiles which had been labeled as either monoclonal or polyclonal to calculate scores. A higher score for this sample means that the algorithm is more likely to predict the sample as a polyclonal sample. Finally, we used R package pROC (Robin et al., 2011) to calculate precision, recall, and the area under the receiver operating characteristic curve (AUC).

Gene set enrichment analysis (GSEA) and pathway enrichment analysis were performed using R package clusterProfiler (Yu et al., 2012). The ssGSEA algorithm was used to calculate the activity score of the gene sets for each cell based on GSVA package (Hanzelmann et al., 2013). We obtained 13 cancer functional states from CancerSEA (Yuan et al., 2019), as well as a collection of 50 cancer Hallmark gene sets downloaded from MsigDB (Liberzon et al., 2015). The signature gene lists of macrophages (e.g. Angiogenesis, Phagocytosis, MHC I molecules, MHC II molecules; Supplementary Table S5) and T cells (e.g. Cytotoxic, Exhaustion, Terminally exhaustion, Proliferation, TCR signaling pathway, Type II interferon response; Supplementary Table S6) were collected from previous studies (Rooney et al., 2015; Long et al., 2022; Li et al., 2020; Jiang et al., 2021; Zhang Y. et al., 2022).

The R package CytoTRACE v.0.3.3 (47) was applied to predict the differentiation state of cells from the single-cell RNA-seq (scRNA-seq) profiles. The CytoTRACE score range from 0 to 1, with higher scores being associated with greater stemness (less differentiation) and vice versa, which help us determine the starting point of differentiation. We employed Monocle2 (Qiu et al., 2017) to discover the differentiation trajectory for T cells and endothelial cells. The top 2,000 DEGs between the cell subgroups were used to create the DDRTree. Utilizing the root state determined by the CytoTRACE, each cell was given a pseudotime value using the “order_cells” function. The cells were ordered and visualized with the “plot_cell_trajectory” function.

The SCENIC (Aibar et al., 2017) was used to construct TF regulatory network for malignant cells in different subclonal complexity. Briefly, GRNBoost2 estimated the co-expression network and RcisTarget was used to identify the regulons. The regulon activity was then measured by AUCell for every cell. The differently activated TFs regulons between different subclonal complexity were identified by wilcox.test(). The TF-gene regulatory network was visualized by Cytoscape (Shannon et al., 2003).

To explore the prognostic power of TFs, we downloaded two gene expression datasets from cBioPortal (Gao et al., 2013) and GEO (GSE76427) for HCC patients, which contained clinical information. HCC patients from TCGA were treated as a training dataset, and HCC patients from GSE76427 were used as an independent validation dataset. A univariate Cox proportional hazards regression analysis was performed to evaluate the association between the critical TFs and patients’ OS. Only TFs with p-value <0.05 were selected to further conduct variable selection according to stepwise Cox proportional hazards regression analysis. Finally, we created a risk-score formula based on TFs expression weighted by Cox regression coefficients. R i s k   s c o r e = ∑ i N c o e f i × e x p r i where N was the number of selected TFs, expr i was the expression value of the ith gene, and coef i was the Cox regression coefficient of the ith gene in the univariate Cox regression analysis. Patients were divided into high-risk and low-risk group based on the median risk score. KM curves were performed to compare the OS and DFS between two groups. To evaluate whether risk score could be an independent prognostic factor, Multivariate Cox proportional hazards regression model was conducted with risk score, age, sex, height, weight, race as covariates.

We applied CellPhoneDB (Efremova et al., 2020) to infer cell-cell communication between different cell types. The specific ligand-receptor interactions between cell types were identified based on permutation test. The interaction score represented the total mean of the individual ligand-receptor partner average expression value in the corresponding interacting pairs of cell types. To predict ligands driving the transcriptomic changes of target cells, we further performed NicheNet analysis (Browaeys et al., 2020) between RBP4+_EC and T cells with the “nichenet_seuratobj_aggregate” function. Specifically, the differentially expressed genes (DEGs) among T cells between different subclonal complexity with high complexity group as control were used as gene sets of interest and all expressed genes in T cells were used as background of genes.

We used CytoSig algorithm (Jiang et al., 2021) to quantify the signaling activity of cytokines for T cells. The CytoSig was a linear model predicting cytokines signaling activity score for each cell and covered 51 cytokines provided by the curated training data. We implemented CytoSig for single-cell transcriptional profiles of T cells in Python 3.

Continuous variables were compared by the Wilcoxon rank-sum test and categorical variables were compared by the chi-squared test. Pearson correlation was calculated by cor () and cor. test (). The log-rank test was used to compare Kaplan–Meier curves. When p < 0.05, the differences were considered to be significant. All statistical analyses were conducted using R software version 4.1.2 (http://www.rproject.org).

We downloaded 30 untreated primary HCC samples from four cohorts of single-cell RNA-seq data (Sharma et al., 2020; Lu et al., 2022; Ma et al., 2021; Ma et al., 2022) (Supplementary Table S1). After quality control, we obtained a total of 66,609 high-quality cells with a median of 2,220 cells per sample. Cell types were annotated by canonical cell markers, identifying 31,832 hepatocytes and 34,777 non-malignant cells (Figure 1A; Supplementary Figure S1A). We inferred copy number variations (CNVs) based on scRNA-seq data using InferCNV((Project. iotTC)) (Supplementary Figure S1B). Malignant and non-malignant hepatocytes were distinguished by their accumulated CNV scores (see Section Methods, Supplementary Figure S1C), ultimately identifying 30,955 malignant cells.

We performed clustering and UMAP visualization separately for malignant and non-malignant cells. Consistent with previous findings, non-malignant cells exhibited distinct separation by cell type but showed a mixture of different samples (Supplementary Figure S1D; Figure 1B), whereas malignant cells exhibited high inter-tumor heterogeneity (Figure 1C) (Puram et al., 2017). The transcriptomic similarity among inter-tumor cells was lower than that among intra-tumor cells (Figure 1D). However, the similarity of intra-tumor malignant cells varied across patients; for instance, HCC04T showed relatively high similarity, while H72 exhibited lower similarity among malignant cells (Figure 1E).

To identify underlying transcriptional programs, we applied non-negative matrix factorization (NMF) (Gaujoux and Seoighe, 2010) to each sample. The identified gene programs tend to be mutually exclusive among malignant cells within sample, and were negatively correlated (Supplementary Figure S1E–G). A total of 98 gene programs were identified across the 30 tumor samples (Supplementary Table S2). We then performed hierarchical clustering on these 98 expression programs, extracting three meta-programs (Figure 1F; Supplementary Table S3) that encompassed more than half of the samples. These three meta-programs were associated with distinct signaling pathways: Meta-program 1 was enriched in cell cycle-related pathways, Meta-program 2 was enriched in metabolic pathways, and Meta-program 3 was enriched in stress response-related pathways (Figure 1G).

To explore intra-tumor subclonal complexity, we proposed a metric to quantify CNV complexity based on the discrete CNVs profiles inferred from InferCNV (Figure 2A, method). We quantified CNV complexity for each HCC sample and classified samples into high (H) and low (L) groups based on the median CNV complexity score. We found that CNV complexity scores were significantly correlated with two other heterogeneity scores proposed by (Ma et al., 2019) and (Guo et al., 2022) respectively (Figures 2B, C). Moreover, the CNV complexity score was more robust than the other two heterogeneity scores, both of which were negatively correlated with the number of the cells (Figures 2D, E). Simulation data further demonstrated the robustness of the CNV complexity score (Supplementary Figure S2A).

Finally, we evaluated the three methods on simulated datasets generated in silico (Supplementary Figure S2B, C). The CNV complexity score showed excellent performance in reflecting the degree of subclonal complexity, as demonstrated by a receiver operating characteristic (ROC) curve with the highest AUC value, which was superior to the other methods (Figure 2F). Overall, our constructed CNV complexity score demonstrates excellent stability and sensitivity, and more importantly, it is not influenced by the number of cells, providing a more accurate quantification of subclonal complexity in samples.

To explore the phenotypic plasticity underlying different subclonal complexity, we used cancer functional gene sets (Supplementary Table S4) along with ssGSEA and identified 37 gene sets with significant differences between the H and L groups, with 34 functional gene sets being significantly upregulated in the L group (p < 0.05, Wilcoxon rank-sum test, Figure 3A; Supplementary Figure S3). For instance, the L group was significantly enriched for hypoxia, invasiveness, metastasis, emphasizing varying signaling activity between the H group and L group. Malignant cells are susceptible to epithelial-mesenchymal transition under hypoxic conditions, gaining invasive and metastatic abilities, accompanied by a transition from low to high stemness. Consistent with hypoxia, invasion and metastasis activities mentioned above, it was discovered that the stemness of malignant cells in L group was also significantly higher than that in the H group (p = 0.0049, Wilcoxon rank-sum test, Figure 3B) (Gulati et al., 2020), which highlights that stemness of malignant cells is closely related to hypoxia, invasiveness and metastasis.

We further investigated the cell lineage origin of malignant cells using scHCL(Han et al., 2020), and found significantly different distributions of cell lineages between the H and L group (p < 0.05, chi-squared test, Figure 3C), with a significantly higher proportion of cell lineage overexpressing FGB in the H group (p < 0.05, Wilcoxon rank-sum test, Figure 3D). FGB was significantly upregulated in malignant cells from the H group (p < 2.22e-16, Wilcoxon rank-sum test, Supplementary Figure S4A). Altered expression of this gene leads to several disorders, including afibrinogenemia, dysfibrinogenemia, and thrombotic tendency. In summary, malignant cells from samples with varying CNV complexity might undergo diverse transcriptional reprogramming, resulting in different cancer functional states and cell lineage origins.

To investigate the regulatory mechanisms contributing to the different cell lineages and functional phenotypes of tumor cells between the H and L group, we employed SCENIC to identify transcription factors (TFs) regulatory modules (Figure 3E). We then constructed TFs regulatory networks corresponding to the H and L groups using differentially expressed and differential CNV genes targeted by identified TFs (Supplementary Figure S4B; Supplementary Figures S5A, B). Functional enrichment analysis of the TF regulatory network revealed that the upregulated genes in the L group’s TFs regulatory network were enriched in pathways such as hypoxia, apoptosis, TNFA signaling via NFkB and regulation of Wnt signaling. In contrast, the upregulated genes in the H group were involved in pathways regulating epithelial cell proliferation and metabolism-associated pathways (Supplementary Figure S4C).

Among these TFs, we identified four transcription factors MAFG, YY1, EZH2, AR, associated with overall survival (OS) which were upregulated in the H group. Samples were stratified into high-risk and low-risk groups based on the median risk score quantified by these four TFs. Kaplan-Meier curve analysis demonstrated poorer survival in the high-risk group compared to the low-risk group in terms of OS and disease-free survival (DFS) (OS: p < 0.0001, DFS: p < 0.0001, log-rank test, Figure 3F). Further multivariate Cox proportional hazards regression analysis confirmed the independent prognostic value of the risk score (Figure 3G; Supplementary Figure S4D). Finally, we validated the prognostic power of the 4-TF signatures in an external validation dataset (GSE76427). The KM curves of the two groups were significantly different (p = 0.0016, log-rank test, Supplementary Figure S4E) and the risk score showed the independent prognostic value (Supplementary Figure S4F).

To investigate the characteristics of cell-cell communications within the tumor microenvironment under different complexity of subclonal structure, we performed cell-cell interaction analysis using CellPhoneDB (Efremova et al., 2020) between malignant cells and non-malignant cells in the H and L groups, respectively. We found that fibroblast and endothelial cells strongly interacted with tumor cells in the L group rather than the H group (Figure 4A). Additionally, myeloid cells, including macrophage, monocytes and dendritic cells (Supplementary Figure S6A, B), showed strong interactions with tumor cells in both the H and L groups, with increased communications observed in the L group (Figure 4B).

We further explored the ligand-receptor interactions between tumor cells and myeloid cells in the L group and H group. We found that many chemokine-mediated interactions were specifically present in the L group, such as the CCL4L2-VSIR, FCGR2A-CXCL9, CCL3L3-CCR1, CCL3L1-CCR11, IDE-CCL23, and CXCL11-DPP4, CCL2-CCR10 pathways which suggests that tumor cells in the L group might recruit more myeloid cells. In the H group, there were also a few specific interactions, such as the interactions between the dendritic cell-expressed ligands CCR6 and CXCR3 with the tumor cells-expressed receptor CCL20 (Figure 4C). Additionally, we confirmed that coinhibitory interactions mediated by pathways like CSF1-SIRPA, LAG3-FGL1, CCL4L2-VSIR, and LGALS9-HAVCR2 were specifically activated in the L group (Figure 4D). Similar results were also observed in the costimulatory-mediated interactions between malignant cells and myeloid cells (Supplementary Figure S6C). Interestingly, the L group was enriched for both coinhibitor- and costimulator-mediated ligand-receptor interaction pairs. In summary, these results revealed that varying subclonal complexity could alter the cell-cell interactions landscape.

Given the differential interaction strength of myeloid cells and tumor cells between the H and L groups, we investigated the functional characteristics of myeloid cells. Myeloid cells from the L group showed upregulation of MHC I and MHC II molecules and higher angiogenesis and phagocytosis signature scores (p < 0.05, Wilcoxon rank-sum test, Supplementary Table S5; Figures 4E, F). Based on gene set enrichment analysis (GSEA) of the differentially expressed genes identified in myeloid cells between the L and H groups (Supplementary Figure S6D), more hallmark signaling pathways were enriched in myeloid cells from the L group, including hypoxia, EMT, TNFA signaling via NFkB, DNA repair, apoptosis, and angiogenesis (Figure 4G; Supplementary Figure S6E, F). These results were consistent with what we observed in tumor cells from samples with different complexity of subclonal structure.

Although there was no significant difference in tumor-T cells interaction strength between the H and L groups, T cells exhibited significantly different activation statuses based on ssGSEA analysis of functional gene sets (Supplementary Table S6). Both CD4⁺ and CD8⁺ T cells from the L group had higher scores for TCR signaling pathway, exhaustion and proliferation but lower type II interferon response score (p < 0.05, Wilcoxon rank-sum test, Figure 5A). Moreover, cytotoxicity-related genes (e.g., GZMA, GZMB, GZMH, NKG7, PRF1) and immune checkpoint genes (e.g., LAG3, PDCD1, LCK) were significantly upregulated in T cells from the L group (p < 0.05, Wilcoxon rank-sum test, Supplementary Figure S7A, B). The scores of these functional gene sets were significantly correlated with CNV complexity score, implying the potential of subclonal complexity in predicting the activation status of T cells (Figure 5B).

We performed subclustering on T/NK cells, resulting in seven cell subpopulations: CD4_CCR7, CD4_IFNG+, CD4_FOXP3+_Treg, CD8_GZMK+, CD8+_TNSF9+, proliferative T cells, and NK cells (Figure 5C; Supplementary Figure S7C). The proliferative T cells were further annotated according to their similarity with other T subtypes and NK subtype (see methods, Figure 5D). Based on trajectory analysis of CD4⁺ and CD8⁺ T cells, we found that both CD4⁺ and CD8⁺ T cells from the L group tended to locate in the intermediate of the trajectory whereas those from the H group were primarily located at the terminal of trajectory (Figure 5E, F). Additionally, T cells from L group showed increased CytoTRACE score, indicating a heightened activation status of T cells in the L group (CD8+: p = 0.043, CD4+: p = 0.075, Wilcoxon rank-sum test, Figures 5E, F; Supplementary Figure S7D, E).

To investigate the regulatory mechanisms underlying the different activation states of T cells, we estimated the signaling activities of 51 cytokines on T cells and their correlation with T cell proliferation using the CytoSig model (Jiang et al., 2021). The activity of the cytokine HMGB1 showed a significant positive correlation with the proliferative score of T cells (p = 0.014, R = 0.45, Figure 6A). Moreover, HMGB1 signaling activity was significantly higher in the L group compared to the H group (p = 0.011, Wilcoxon rank-sum test, Figure 6B). HMGB1 was primarily expressed in endothelial and fibroblast cells (Figure 6C). However, the expression levels of HMGB1 in endothelial cells, rather than fibroblast cells, were upregulated in the L group (p = 0.0023, Wilcoxon rank-sum test) and significantly correlated with T cell proliferation (Figures 6D, E; Supplementary Figure S8A). Additionally, the expression levels of HMGB1 in endothelial cells was positively correlated with the activity scores of T cell exhaustion and TCR signaling pathway, and negatively correlated with the activity score of type II interferon response score (Supplementary Figure S8B).

We further dissected the subtypes of endothelial cells, identifying seven distinct endothelial cell subtypes (Figure 6F; Supplementary Figure S8C). The proportion of RBP4 positive endothelial cells (RBP4+EC) was significantly higher in the L group compared to the H group (p < 0.05, Wilcoxon rank-sum test, Figure 6G). Moreover, HMGB1 expression which was present only in RBP4+ ECs was significantly increased in the L group compared to the H group (p = 0.009, Wilcoxon rank-sum test, Figure 6H; Supplementary Figure S8D). Trajectory analysis demonstrated that RBP4+ECs persisted throughout the differentiation trajectory and showed highly expressed HMGB1 in the L group along with the pseudotime (Figures 6I–K). Moreover, CytoTRACE analysis suggested that RBP4+ECs in the L group showed a less differentiated status (p = 0.0079, Wilcoxon rank-sum test, Supplementary Figure S8E–G).

To further investigate the regulatory mechanisms of RBP4+ECs on T cell activation, we identified regulators from RBP4+ECs for T cell activation through Nichenet analysis (Browaeys et al., 2020) (Figure 7A). We found that RBP4+ECs regulated the expression of 29 target genes in T cells by upregulating HMGB1 in the L group, including IFNG, which enhances tumor immune function and promotes T cell exhaustion (Benci et al., 2019), and CCL3, which facilitates dendritic cells (DC) recruitment for antigen presentation and T cell activation (Castellino et al., 2006). Interaction analysis between RBP4+ECs and T cells showed more and stronger chemokine-mediated interactions between RBP4+ECs and CD8+T in the L group including CXCR6-CXCL16, NR3C1-CCL2, CXCL12-CXCR3, CXCL12-CXCR4, and PGRMC2-CCL4L2 (Figure 7B). Similar results were observed in the interactions between RPB4+ECs and CD4⁺ T cells (Figure 7C).