Molecular characterization of HDF in CRC
The colorectal cancer (CRC) single cell transcriptome dataset was the data base for exploring the cell types and interactions between cells in the TME, and it could be found that the main cell types in the TME were Epithelial cells, Stromal cells, Myeloid cells, B lymphocyte cells, T lymphocyte cells and Mast cells. CRC tumor cells were mainly derived from the malignant transformation of epithelial cells (Fig. 1A). Subsequently, assessment of the interactions among these cells using the CellChat tool revealed strong interactions between Epithelial cells and Stromal cells. Moreover, the interactions between these microenvironment cells were widespread (Fig. 1B). HDAC was a kind of proteases that played an important role in structural modification of chromosomes and regulation of gene expression. The mRNA expression of genes involved in the HDAC pathway was found to differ significantly between cell types from CRC single cell transcriptome data. High mRNA expression of many related genes was showed in the cell subtypes of Epithelial cells (Fig. 1C), which was consistent with the previously reported overexpression of HDAC related genes in cancer4. Based on HDAC related genes, cell subtypes were further identified using NMF for the major cell types in the TME, including four cell subtypes for stromal cells, five cell subtypes for macrophage cells, four cell subtypes for B cells, and five cell subtypes for T cells. It could be found that inter-tumor heterogeneity among different patients was reflected both in the proportional differences of major cell types and cell subtypes (Fig. 1D).
Characterization of HDF-mediated stromal cell subtypes in CRC
Pseudotime trajectories analysis is an important indicator of the distance of cell differentiation using the single cell transcriptome, and the cells of TME transfers by corresponding state transitions during tumor progression. In order to investigate the role of HDAC in the transition of stromal cell state, we used the Monocle2 R package for pseudotime trajectories analysis of stromal cells, and then analyzed the relationship between HDAC-related genes and pseudotime trajectories, it could be found that the high expression of some genes significantly appeared at the head or tail of the pseudotime trajectories. The main highly expressed genes that appeared in the tail of the pseudotime trajectories were BRMS1, CHD3, HDAC2, MORF4L1, PRKD1, HOPX, MAPK8, DR1, MORF4L2 and PHB. On the contrary, the main highly expressed genes that appeared in the head of the pseudotime trajectories were HDAC10, PRMT6, C6orf89, RCOR2, SIRT3, MTA3, TRERF1, LRRK2, HDAC7 and PRDM5 (Fig. 2A). A comparison of HDAC-mediated interactions between subtypes of stromal cell and epithelial cells showed no significant differences between each of the four fibroblast subtypes and epithelial cells (Fig. 2B). The proportional distribution of the four fibroblast subtypes in tumor tissues compared to normal tissues was also not significantly different (Fig. 2C). However, Fib_C2 and Fib_C4 subtypes of fibroblast had a slightly up-regulated trend in the proportional distribution. Based on the KEGG database, functional enrichment analysis of genes specifically highly expressed in the four fibroblast subtypes revealed that the Fib_C2 subtype was enriched for TNF signaling pathway, IL-17 signaling pathway and Osteoclast differentiation. Fib_C4 subtype was enriched for Autophagy animal. Fib_C1 subtype was enriched for Protein digestion and absorption, ECM-receptor interaction, and Proteoglycans in cancer. Fib_C3 subtype was enriched for Focal adhesion, Hypertrophic cardiomyopathy and Dilated cardiomyopathy (Fig. 2D). Further correlation analysis of HDAC-mediated stromal cell subtypes revealed that the marker proteins significantly associated with Fib_C2 subtype and Fib_C4 subtype were mainly ADAMDEC1 and PLVAP. Moreover, the marker proteins significantly associated with Fib_C1 subtype and Fib_C1 subtype were mainly CREM (Fig. 2E). The genes were highly expressed in Fib_C2 subtype were mainly IRF1, CCL2, SOCS3, JUNB and FOSB (Fig. 2F)
Characterization of HDF-mediated myeloid cell subtypes in CRC
To investigate the function of HDAC in myeloid cell state transition, the CellChat tool was used to compare the interaction between five HDAC-mediated macrophage subtypes and epithelial cells, and it could be found that the Mac_C3 subtype of macrophages with SFPQ as a highly expressed gene had a strong interaction with epithelial cells (Fig. 3A), and macrophage Mac_C3 subtype had a strong correlation with macrophage M1 type (anti-cancer and pro-inflammatory type) (Fig. 3B), which suggested that HDAC had the potential ability to regulate macrophage functional polarization in the microenvironment of CRC tumors. Mac_C1 subtype with MORF4L1 as a highly expressed gene, macrophage Mac_C2 subtype with VEGFA as a highly expressed gene, and macrophage Mac_C4 subtype with PHB as a highly expressed gene also had some interactions with epithelial cells (Fig. 3A), and macrophage Mac_C1 subtype and Mac_C4 subtype had a strong correlation with macrophage M2 type (pro-oncogenic and anti-inflammatory type) (Fig. 3B), which further suggested the important role of HDAC in the regulation of macrophage functional polarization in the microenvironment of CRC tumors. The functional polarization of macrophages had a very important role in CRC tumor progression, and it could be found that the distribution of M1-type and M2-type features of macrophages in TME was widespread and did not have strong cluster distribution characteristics (Fig. 3C), which suggested that the functional polarization of macrophages in TME was dynamic The macrophage subtypes with GZMB as the highly expressed gene, the macrophage subtype with IL32 as the highly expressed gene and the macrophage subtype with S100A8 as the highly expressed gene showed strong cluster distribution characteristics in the tSNE descending plot (Fig. 3C), which may be related to their corresponding fixed biological functions. Based on the KEGG database, functional enrichment analysis of genes specifically highly expressed in five adult macrophage subtypes revealed that the pathways enriched in the Mac_C3 subtype of macrophages were mainly SPLICEOSOME, while the pathways enriched in the Mac_C1 and Mac_C4 subtypes of macrophages were mainly PYRUVATE METABOLISM, LIMONENE AND PINENE DEGRADATION, GLUTATHIONE METABOLISM, PROTEASOME, PROPANOATE METABOLISM, and FRUCTOSE AND MANNOSE METABOLISM (Fig. 3D).
Characterization of HDF-mediated B-cell subtypes in CRC
In order to investigate the role of HDAC in the transition of B cells state, we used the Monocle2 R package for pseudotime trajectories analysis of B cells, and then analyzed the relationship between HDAC-related genes and pseudotime trajectories, it could be found that the high expression of some genes significantly appeared at the tail of the pseudotime trajectories. The main highly expressed genes that appear in the tail of the pseudotime trajectories were SDR16C5, PER1, CAMK2D, LRRK2, NCAPG2, SUV39H1, BCL6, BRMS1L, UCN, CHD3, MTA2 (Fig. 4A). Comparison of the interactions between the five HDAC-mediated B-cell subtypes and epithelial cells using the CellChat tool revealed that the interactions were all relatively weak (Fig. 4B). Subsequent comparison of the correlation between the five HDAC-mediated B-cell subtypes and protein marker-labeled B-cell subtypes revealed that the five B-cell subtypes showed a more consistent correlation, all showing strong correlation with IgG and IgA proteins, and weak correlation with CD19 and CD20 proteins (Fig. 4C). Intercellular interactions were further explored using the NichNet tool, which considereed not only ligand-receptor interactions but also target genes activated after ligand-receptor interactions, and it could be found that the more active ligands in B cells were MANF, LGALS3, ARF1, RPS19, TIMP1, CD99 and GSTP1, while the more active target genes in epithelial cells were CCND1, BID, CCL2, EDN1, TIMP1, NQO1and JUNB. The activation of target genes had some ligand genes, suggesting that the B cell-epithelial cell interaction process had a more complex multiple network relationship (Fig. 4D). Assessment of the differentiation potential of B cells using the CytoTRACE tool revealed that the differentiation potential of B cells was characterized by a more significant cluster distribution, but there was a poor agreement between the cluster characteristics of B cells differentiation potential and HDAC-mediated clusters (Fig. 4E). Based on the KEGG database, functional enrichment analysis of genes specifically highly expressed in five adult B-cell subtypes revealed that the main enriched pathways in B-cell subtypes Bcell_C1 and Bcell_C4 were Ribosome, Coronavirus disease-COVID-19, Allograft rejection, Leishmaniasis, Asthma, Intestinal immune network for IgA production, Viral myocarditis, while B-cell subtypes Bcell_C3 and Bcell_C5 were mainly enriched for Thyroid hormone synthesis, Protein export, Protein processing in endoplasmic reticulum (Fig. 4F).
Characterization of HDF-mediated T-cell subtypes in CRC
Five subtypes of T cells were annotated: KLRB1-positive T cells, CD8-positive T cells, regulatory T cells, cytotoxic T cells, and KRT18-positive T cells (Fig. 5A). To investigate the role of HDAC in T cell state transition, we compared the correlation between HDAC-mediated T-cell subtypes and marker protein-annotated T-cell subtypes, and found that HDAC-mediated T-cell subtype Tcell_C3 was mainly distributed in CD8-positive T cells and cytotoxic T cells, while other T-cell subtypes did not have significant distribution differences (Fig. 5B), suggesting that HDAC might play a regulatory role in the killing process of T cells against tumor cells. Further comparison of the interactions between T-cell subtypes and epithelial cells mediated by the combination of 25 HDAC-related genes and marker proteins using the CellChat tool revealed both CD8-positive T cells and cytotoxic T cells had strong interactions with epithelial cells (Fig. 5C), and there was no significant difference in the interactions between the five HDAC-mediated T-cell subtypes and epithelial cells in the cluster of CD8-positive T cells, with the same trend in cytotoxic T cells. These results further emphasized the important role of HDAC in the process of killing function of T cells. Subsequently, five HDAC-mediated T-cell subtypes were compared for immune response-related scores and gene expression differences, and it could be found that the HDAC-mediated T-cell subtype Tcell_C3 had a strong signal on the Cytotoxic score, Exhaustion score, and Teffect score. High expression of genes in Co -inhibitors included HAVCR2, TGFB1, TIGIT, CD247, CD160, CD244, CD96, LAG3, Co-stimuiations included CD226, TMIGD2, ENTPD1, TNFRSF9, and genes highly expressed in T-function Genes were CD8B, GZMB, PRF1, TBX21. (Fig. 5D). Finally, we compared the correlation between HDAC-mediated B-cell subtypes and marker protein-annotated B-cell subtypes, it could be found that NRT18-positive B-cell subtypes mainly contain HDAC-mediated B-cell subtype Bcell_C3 (Fig. 5E).
The prognosis assessment of HDF-mediated microenvironmental cell subtype features in CRC
To further assess the clinical prognostic significance of genes characterizing the subtypes of major cell types in the HDAC-mediated CRC tumor microenvironment, six GEO datasets (GSE72970, GSE39582, GSE39084, GSE17538, GSE17536, and GSE103479) with overall survival (OS) information were collected, and using univariate Cox analysis to assess their relationship on OS, it could be found that most of the genetic characteristics had some prognostic significance (Fig. 6ABCD). Eight GEO datasets (GSE72970, GSE39582, GSE38832, GSE33113, GSE17538, GSE17536, GSE14333, GSE103479) with progression-free survival (PFS) information were collected, and similar results to OS assessment could be found (Fig. 6EFGH).
The predictive power assessment of immunotherapy benefit for HDF-mediated micro-environmental cell subtype characterization in CRC
To further assessed the clinical immunotherapeutic significance of HDAC-mediated subtype-specific genes of major cell types in TME, the GSE78220 and IMvigor210 immunotherapy datasets were collected and evaluated accordingly. In the GSE78220 dataset, the risk score allowed to divide the sample into two groups, High_RiskScore and Low_RiskScore (Fig. 7A). There was a significant difference in the survival curves between the high and low risk groups (Fig. 7B), and the proportion of immunotherapy responsed between the high and low risk groups trend (Fig. 7C), suggesting that patients in the low-risk score group were more likely to benefit from immunotherapy. A consistent trend (Fig. 7DEF) can be found in the IMvigor210 dataset.