Su-Su Zheng1†
Bo-Heng Zhang
Cheng Huang
Tong-Chun Xue- 1Department of Hepatic Oncology, Xiamen Clinical Research Center for Cancer Therapy, Zhongshan Hospital, Fudan University (Xiamen Branch), Xiamen, China
- 2The Liver Cancer Institute, Zhongshan Hospital and Shanghai Medical School, Fudan University, Key Laboratory for Carcinogenesis and Cancer Invasion, The Chinese Ministry of Education, Shanghai, China
- 3Center for Evidence-based Medicine, Shanghai Medical School, Fudan University, Shanghai, China
- 4Department of Liver Surgery, Liver Cancer Institute, Zhongshan Hospital, and Key Laboratory of Carcinogenesis and Cancer Invasion (Ministry of Education), Fudan University, Shanghai, China
Myeloid cells are physiologically related to innate immunity and inflammation. Tumor-associated myeloid cells gained increasing interest because of their critical roles in tumor progression and anticancer immune responses in human malignancies. However, the associations between tumor-associated myeloid cell-related genes and hepatocellular carcinoma have yet to be revealed. Here, through the integrating analysis of bulk and single-cell RNA (scRNA) sequencing of public HCC samples, we developed a gene signature to investigate the role of HCC-specific myeloid signature genes in HCC patients. We firstly defined 317 myeloid cell marker genes through analyzing scRNA data of HCC from the GEO dataset. After selecting the differentially expressed genes, eleven genes were also proved prognostic. Then we built a gene signature from the TCGA cohort and verified further with the ICGC dataset by applying the LASSO Cox method. An eight genes signature (FABP5, C15orf48, PABPC1, TUBA1B, AKR1C3, NQO1, AKR1B10, SPP1) was achieved finally. Patients in the high risk group correlated with higher tumor stages and poor survival than those in the low-risk group. The risk score was proved to be an independent risk factor for prognosis. The high risk group had higher infiltrations of dendritic cells, macrophages and Tregs. And the APC co-inhibition, T cell co-inhibition pathways were also activated. Besides, the risk score positively correlated with multidrug resistance proteins. In conclusion, our myeloid cell marker genes related signature can predict patients’ survival and may also indicate the levels of immune infiltration and drug resistance.
Introduction
Hepatocellular carcinoma accounts for the majority of liver cancer and is one of the leading cause for cancer-related deaths globally (Sung et al., 2021). Further, most patients lost the chance of surgery because of vascular invasion or even metastasis. For these advanced HCC, systemic treatments including immune treatment are highly recommended (Sharafi et al., 2022). Among various solid tumors, immune treatment is especially widely used in HCC. Immune checkpoint inhibitors treatment was highly recommended to combine with TKI or bevacizumab. The combination of ICIs and tyrosine kinase inhibitor/bevacizumab therapy has gained synergistic effects (Finn et al., 2020; Roskoski, 2022). However, until now, there are no well-established indicators of immunotherapy response for HCC because of limited knowledge about the immune-related tumor microenvironment.
Myeloid cells are physiologically related to innate immunity and inflammation for cancer. Tumor associated myeloid cells, as an important part of the TME, can be divided into monocytes/macrophages, myeloid-derived suppressor cells, dendritic cells, and granulocytes (Bassler et al., 2019). In recent decades, the exact roles of TAMCs have gained increasing interest since these cells are important indicators for treatment efficacy and disease prognosis in patients. For example, the density of TAMCs are commonly negatively correlated with patients’ survival, while some special subtypes of TAMCs (M1 subtype macrophage, N2 subtype neutrophil) are proved to be associated with better prognosis (DeNardo and Ruffell, 2019; Engblom et al., 2016; Fridlender et al., 2009; Locati et al., 2020; Nakamura and Smyth, 2020). In HCC, our previous work showed that the high post-treatment neutrophil-to-lymphocyte ratio indicated a higher risk of metastasis for patients undergoing transarterial chemoembolization (Xue et al., 2015). Especially, the TME infiltrating neutrophils, as one kind of innate immune cells have been shown to be involved in treatment and prognosis among HCC patients (Chen et al., 2021). Due to the flexibility of TAMCs, it is necessary to develop gene signatures to evaluate their roles in immune infiltration and predict patients’ prognosis. Through high-throughput sequencing technologies, we are now able to perform analysis to comprehensively catalog myeloid cells related genes in cancers. Liu et al. developed a myeloid cells related gene signature containing 5 genes, which was valuable in evaluating the prognosis and immunity for head and neck squamous cell carcinoma (Liu et al., 2021). However, their myeloid signature genes were generated from published literature of various tumors. To the best of our knowledge, the HCC specific myeloid cells marker genes related signature has not yet been reported.
Through single-cell RNA-sequencing (scRNA-seq) method, we may further uncover the molecular characteristics of HCC associated myeloid cells in the TME (Chen et al., 2019). In the present study, we firstly identified the tumor associated myeloid cells related marker genes with scRNA-seq analysis of HCC samples from the Gene Expression Omnibus (GEO) dataset. And then, a myeloid cell marker genes related signature was built and validated as an independent risk factor. Besides, the risk score was proved to be associated with immune infiltration and immune-suppressive microenvironment. Moreover, the relationship between risk scores and drug sensitivity was further evaluated.
Materials and methods
Data source
ScRNA-seq data of 7 tumor samples from 2 HCC patients was obtained from the GEO dataset, namely, GSE112271. RNA sequencing data of (fragment per kilobase million, FPKM) was obtained from the TCGA database (https://portal.gdc.cancer.gov/repository), which included 374 HCC patients. The International Cancer Genome Consortium (ICGC) database (https://icgc.org/) contained 231 HCC samples. And the TCGA database was used for model construction, while the ICGC database for model validation. All datasets were available from public websites, and ethics approval was confirmed to be obtained from original studies.
Achieving HCC-specific myeloid cell marker genes
The scRNA-seq data from GSE112271 has been fully described by Bojan Losic et al. (Losic et al., 2020). Here, we combined data from all 7 samples derived from 2 HCC patients. The data was then analyzed with the R software by using the Seurat package (Butler et al., 2018). Quality control was performed as the following: gene numbers more than 500, mitochondrial gene percentage no more than 20, and total UMI counts more than 1,000 (Guan et al., 2022). The “RunHarmony” function in R package “harmony” was then applied to remove the batch effects (Korsunsky et al., 2019). Principal component analysis (PCA) was conducted with the top 2000 variable genes. And cell clusters were defined with the method of shared nearest neighbor (SNN), and the resolution was set as 0.5. The uniform manifold approximation and projection (UMAP) analysis was introduced for visualization. The differentially expressed genes (DEGs) were identified with the “FindAllMarkers” function. The cutoff threshold was set as an adjusted p-value <0.01 and |log2 (fold change)| >1. Cell clusters were annotated by using data from published literature as well as the Human Primary Cell Atlas (Losic et al., 2020; Mabbott et al., 2013; Sun et al., 2021).
Construction and validation of a myeloid marker genes related gene signature
The cutoff threshold was set as an adjusted p-value <0.01 and |log2 (fold change)| >1.5 for identifying DEGs between tumor tissues and non-tumor tissues with the “limma” R package in TCGA dataset. The prognostic value of myeloid marker genes was further analyzed by Univariate Cox. Intersection genes were acquired from DEGs and prognostic genes and were then used for model construction with the method of least absolute shrinkage and selection operator (LASSO) Cox regression (Tibshirani, 1997).Then, an eight myeloid marker genes related gene signature was constructed. For each patient, a risk score was calculated by using his gene expression and the regression coefficient. According to the median risk score, all patients were divided into high and low-risk groups. The “survminer” and “timeROC” R package was adopted for survival analysis and predictive accuracy. The ICGC dataset was used for verification with the same formula. Besides, clinical factors such as age, gender, and tumor stage were manually extracted from both the TCGA and ICGC datasets, and were applied for univariate and multivariable Cox regression analysis.
Patients and immunohistochemical staining
Seventy tumor tissues were obtained from 70 HCC patients in our liver cancer institute. SPP1, CD11b, CD4, and CD8 were selected for protein validation by immunohistochemistry. The expression intensity of SPP1 was evaluated according to a previous research (Wu et al., 2010). The number of positive cells were counted as previously reported (Deng et al., 2017). The study was approved by the Ethics Committee of Zhongshan Hospital, Fudan University.
Functional enrichment analysis
Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed by gene set enrichment analysis (GSEA) with GSEA software 4.0.1 according to the high and low-risk group. And the top 50 pathways between different groups were exhibited.
Immune infiltration, TME and drug sensitivity analysis
Immune infiltration analysis was conducted with the “GSVA” R package. Single-sample gene set enrichment analysis (ssGSEA) was adopted for calculating the infiltration scores of immune cells and immune-related pathways (Newman et al., 2019). The association between risk score and pan-cancer immune infiltrates was also compared (Tamborero et al., 2018). Stromal and immune cell infiltration scores were also analyzed with the “ESTIMATE” R package (Yoshihara et al., 2013). Correlations between multidrug resistance genes and risk scores were also analyzed.
Statistical analysis
DEGs between different groups were compared by Wilcoxon test. The ssGSEA scores between different groups were compared by Mann-Whitney test, and the adjusted p-values were generated by using Benjamini–Hochberg method. Survival analysis was conducted with the Kaplan-Meier method. Univariate and multivariate Cox analyses were utilized for identifying prognostic factors. p-value < 0.05 was considered to be statistical significance if not specified. All data analyses and figures were accomplished with the R software (version 4.1.3).
Results
Identification of HCC-specific myeloid cell marker genes
Through scRNA-seq data analysis of GSE112271, we obtained 30054 eligible cells from seven samples of two HCC patients and were further divided into 15 clusters (Figures 1A, B). Next, we annotated the cell clusters with data from published literature and the Human Primary Cell Atlas. Similar to the results of Bojan Losic et al., 7 cell types were discovered, including Malignant cell (ALB, FGG), myeloid-derived cell (HLA-DQB1, CD68), cancer-associated fibroblasts (ACTA2, TAGL), endothelial cell (VWF, KDR), T cell (CD3D, CD2), B cell (CD79A, IGJ) and NK cell (CD69, GNLY). Cells in cluster 2, 13, and 14 were defined as myeloid-derived cells (Figure 1C). By using the “FindAllMarkers” function, marker genes belonging to different cell type were defined, the top 10 genes were shown in Figure 1D. Finally, we identified 317 HCC-specific myeloid cell marker genes (Supplementary Table S1). The flow chart of the study design and analysis was indicated in Supplementary Figure S1.
FIGURE 1. Myeloid cell marker genes defined by Single-cell RNA sequencing (A) UMAP plot colored by 7 HCC samples. (B) UMAP plot colored by 15 cell clusters. (C) Cell types defined by known cell marker genes. (D) Heatmap exhibiting the top 10 marker genes for each cell type.
Construction and validation of a myeloid cell marker genes related gene signature
Based on the TCGA cohort, we firstly identified 32 DEGs between the tumor and non-tumor tissues. At the same time, by univariate Cox analysis, 11 genes were finally selected to be both differently expressed and prognostic (Figure 2A). The expression of these genes was shown as a heatmap (Figure 2B). The hazard ratio of each gene was shown in Figure 2C. The correlations among different genes were exhibited in Figure 2D. Next, through the LASSO Cox regression analysis, we constructed an 8 genes signature, the calculation function of risk score was set as follows: Riskscore = (0.06*FABP5 expression) + (0.064*C15orf48 expression) + (0.063*PABPC1 expression) + (0.103*TUBA1B expression) + (0.019*AKR1C3 expression) + (0.015*NQO1 expression) + (0.005*AKR1B10 expression) + (0.07*SPP1 expression). We next selected the top 3 genes with highest weights (TUBA1B/SPP1/C15orf48) and validated the correlation between their expression and the infiltration of myeloid cell in HCC with the online tool Timer (https://cistrome.shinyapps.io/timer/). And the results showed that all three genes were positively correlated with the infiltration of macrophage, neutrophil, and dendritic cell (Supplementary Figure S2). Moreover, SPP1 protein level was further investigated by immunohistochemistry, and the association between SPP1 and tumor-infiltrating immune cells was explored. As shown in Supplementary Figure S3, SPP1 was positively correlated with the infiltration of myeloid cells (indicated as CD11b positive cells), and negatively correlated with CD4/CD8 cells.
FIGURE 2. Identification of the candidate Myeloid cell marker genes. (A) Venn plot showing the intersection genes. (B) Heatmap showing the expression of the 11 intersection genes. (C) Forest plots showing the prognostic value of the 11 intersection genes. (D) The correlations between 11 selected genes.
According to the median risk score, patients from the TCGA cohort were divided into the two groups. Figures 3A, B showed the distribution between risk scores and patients’ prognosis. By Kaplan-Meier analysis, the results showed that patients in the high-risk group had worse overall survival (OS) (Figure 3C). Further, the time-dependent ROC results indicated that the 1, 2, and 3-year AUC values were 0.745, 0.661, and 0.626, respectively (Figure 3D).
FIGURE 3. Construction and validation of the gene signature. (A, B) The distribution of risk score and survival status for the TCGA dataset. (C) Survival analysis between the high and low risk groups for the TCGA dataset. (D) ROC curves showing the predictive value of the risk score for the TCGA dataset. (E, F) The distribution of risk score and survival status for the ICGC dataset. (G) Survival analysis between the high and low risk groups for the ICGC dataset. (H) ROC curves showing the predictive value of the risk score for the ICGC dataset.
Furthermore, we assessed the predictive value of the gene signature in the ICGC HCC dataset. ICGC patients were also divided into different groups by applying the same median risk score from the TCGA dataset. Similar to the result of the TCGA dataset, the high risk group indicated poorer OS (Figures 3E–G), and the AUC value was 0.669 for 1 year, 0.677 for 2 years, and 0.650 for 3 years (Figure 3H).
Independent prognostic value of the risk score
By univariate Cox analysis, our results showed that the risk score was significantly correlated with patients’ survival (For TCGA cohort: HR = 1.68, 95% CI = 1.369–2.062, p < 0.001; For ICGC cohort: HR = 2.203, 95% CI = 1.519–3.195, p < 0.001; Figures 4A, B). By multivariable Cox regression analysis, the risk score was proved to be an independent risk factor for both the training and validation cohort (For TCGA cohort: HR = 3.203, 95% CI = 2.016–5.088, p < 0.001; For ICGC cohort: HR = 2.395, 95% CI = 1.164–4.928, p = 0.018; Figures 4C, D). We also evaluated the associations between the risk score and common clinical parameters. Our results showed that the risk score was significantly associated with tumor grade and tumor stage (Figures 5C, D, G).
FIGURE 4. The independent prognostic value of the risk score. (A) Univariate cox regression analysis for the TCGA dataset. (B) Univariate cox regression analysis for the ICGC dataset. (C) Multivariate cox regression analysis for the TCGA dataset. (D) Multivariate cox regression analysis for the ICGC dataset.
FIGURE 5. Box plots showing the correlations between risk score and clinicopathological factors. (A–D) Age, gender, grade and stage from TCGA dataset. (E–G) Age, gender and stage from ICGC dataset.
Functional enrichment analysis between the low and high risk groups
Results of the GSEA analysis showed the 50 most significantly enriched KEGG pathways. The endocytosis, chemokine signaling, leukocyte transendothelial migration, and VEGF signaling pathways were enriched in the high risk group. While peroxisome, peroxisome proliferator-activated receptors, fatty acid metabolism, and nitrogen metabolism pathways were enriched in the low risk group (Figure 6).
FIGURE 6. GSEA analysis exhibiting the top 5 enriched KEGG pathways in the high and low risk groups.
Immune infiltration and TME analysis
By ssGSEA, we further analyzed the correlations between risk score and immune infiltration. We calculated the scores of different immune cells and immune-related pathways. Our results indicated that higher levels of infiltration of aDCs, DCs, iDCs, pDCs, macrophages and Treg were observed in the high risk group both in the TCGA and ICGC cohorts, while the APC co-inhibition, T cell co-inhibition pathways were also activated in these groups. Interestingly, the type II IFN response was inhibited in the high risk group both in the TCGA and ICGC cohorts (Figures 7A–D).
FIGURE 7. Immune infiltration between the high and low risk groups and drug resistance analysis. (A, C) Immune cells and immune-related pathways between the high-risk and low-risk groups in the TCGA dataset. (B, D) Immune cells and immune-related pathways between the high-risk and low-risk groups in the ICGC dataset. (E) Correlations between risk score and different immune subtypes. (F, G) Correlations between risk score and immune score, stromal score. (p values were represented as: ns not significant; *p < 0.05; **p < 0.01; ***p < 0.001.)
The correlations between risk score and pan-cancer immune infiltration subtypes was further evaluated. As reported by David Tamborero, immune infiltrates in solid tumors may be divided into six types (C1-C6 types), which indicated the role from tumor promotion to tumor inhibition (Tamborero et al., 2018). Our results showed that high risk score correlated with tumor promotion subtype, while low risk score implied tumor inhibition subtype (Figure 7E). By applying the ESTIMATE analysis, our results showed that the risk score was weak positive correlated with the immune score and stromal score (Figures 7F, G). In addition, multidrug resistance (MDR) is one of the main causes for treatment failure in HCC, and multidrug resistance proteins (MRPs) mediated multidrug resistance in various cancers including HCC (Ding et al., 2022). Here, we also analyzed the correlation between the risk score and MRPs (MRP1-MRP9). The results showed that the risk score was positively correlated with MRP1, MRP4, MRP5, and MRP7, which indicated that patients in the low risk group may be more likely to derive benefit from chemotherapy or targeted therapy (Figure 8).
FIGURE 8. Scatter plots showing the relationship between risk score and MRP1, MRP4, MRP5, MRP7. (A–D).
Discussion
ScRNA-seq technologies have deeply changed the paradigm to explore the TME. TAMCs, as the main immune cell types, are enriched and highly heterogeneous in tumor tissues. Here, we adopted the scRNA-seq method to explore the roles of HCC related TAMCs. We firstly defined 317 HCC-specific myeloid cell marker genes through scRNA-seq analysis of the GEO dataset. Then a novel prognostic gene signature was developed and validated by the TCGA and ICGC dataset. Similar to our results, another myeloid cell marker gene related gene signature was proved to be prognostic and predictive of immunotherapy response for patients with HNSCC (Liu et al., 2021). Further, the risk score was proved to be correlated with tumor grade and stage. Patients in the high-risk group had shorter OS. Besides, the risk score was associated with a tumor promotion environment. In addition, the risk score was also proved to be positively associated with drug resistance.
Our gene signature contained 8 myeloid cell related genes (FABP5, C15orf48, PABPC1, TUBA1B, AKR1C3, NQO1, AKR1B10, SPP1). All these genes were shown to be overexpressed in HCC tissues and also to be prognostic. 1) FABP5 takes part in fatty acids delivery as an intracellular carrier (Kaczocha et al., 2012). Accumulating evidence has suggested that FABP5 is commonly upregulated in most human malignancies. FABP5 was overexpressed in HCC and involved in the proliferation, migration, and invasion of HCC cells through a FABP5/CREB/miR-889-5p/KLF9 axis (Tang et al., 2022). More recently, FABP5 in monocytes/macrophages was shown to promote lipid accumulation and induction of the inhibitory tumor microenvironment of HCC(Liu J. et al., 2022). 2) C15orf48 is also known as NMES1, the function of this gene is poorly understood. Recently, C15orf48 was found to be overexpressed in response to activation of macrophages and involved in regulating inflammatory cytokines expression (Liu et al., 2009). C15orf48 was also proved to reduce tissue inflammation and immunity and proved to be protective during infection and inflammation (Lee et al., 2021). 3) PABPC1 was proved to take part in miRNA-mediated gene silencing including HCC(Tritschler et al., 2010; Zhang et al., 2015). Interestingly, PABPC1 was shown to regulate immunoglobulin secretion in immune cells (Peng et al., 2017). 4) TUBA1B, encodes the protein of tubulin alpha-1B chain, is the major constituent of microtubules. Study has shown that changes in the relative content of tubulin may regulate neutrophils activation (Rothwell et al., 1993). 5&6) AKR1C3 and AKR1B10 are enzymes that catalyzes redox transformation and play important role in tumor progression (Penning, 2015). AKR1C3 promoted HCC cells proliferation and metastasis through the AKR1C3/NF-κB/STAT3 axis, and was upregulated in HCC tissues. High expression of AKR1C3 correlated with poor survival (Zhou et al., 2021). AKR1C3 also inhibited the ubiquitination of PARP1 and thus resulting in HCC cell proliferation and resistance to Cisplatin (Pan et al., 2022). Similarly, AKR1B10 also promoted HCC progression and drug resistance (Zhang et al., 2022). 7) NQO1 is a flavoprotein, which is important in the cellular response to numerous stresses (Lee et al., 2015). NQO1 was shown to be overexpressed in HCC and correlated with poor survival (Lin et al., 2017). Mechanistically, NQO1 functioned as an agonist at pathways of PI3K/Akt and MAPK/ERK, which promoted HCC cell proliferation and tumor growth (Dimri et al., 2020). 8) SPP1 is a multifunctional gene, which takes part in a variety of cellular processes including cell epithelial transformation, a cytokine participating in the regulation of IL-12/IFN-γ, and the promotion of tumorigenesis (Ashkar et al., 2000; Han et al., 2019; Kashani-Sabet et al., 2017). Interestingly, SPP1 was also proved to be overexpressed in HCC and correlated with patients’ survival as well as an indicator for immunotherapy. More importantly, SPP1 may promote macrophages transition via SPP1-CD44 signaling (Liu L. et al., 2022;Liu Z. et al., 2022;Wang et al., 2019). In glioblastoma, SPP1 was found to be responsible for neutrophil and macrophage infiltration (Atai et al., 2011). Similarly, in this study, we found that SPP1 was positively correlated with myeloid cell infiltration but negatively correlated with CD4/CD8 cell infiltration. In summary, our findings of the signature genes may provide new therapeutic and prognostic targets for HCC.
The association between risk score and immune infiltrations were further explored. We firstly explored its role with the previously reported pan-cancer immune infiltration subtypes (Tamborero et al., 2018). Our results indicated that higher risk score positively correlated with C1 tumor promoting subtype. Next, we adopted the ssGSEA for evaluating the correlations between risk score and immune cell infiltration. Our study indicated that higher levels of DCs, Treg, and macrophages infiltrations were observed in the high-risk group. Besides, the APC co-inhibition, T cell co-inhibition pathways were also activated in the high-risk group. However, the type II IFN response was inhibited. It is widely recognized that a high abundance of Treg contributes to immunosuppression and tumor progression, as well as the tumor associated macrophages (Deng et al., 2021; Zhang et al., 2019). Consistently, type II IFN response has played important roles in tumor surveillance, and IFN-γ, as the only executor of type II IFN response and the most important macrophage stimuli, can induce direct antimicrobial and antitumor effects (Petermann et al., 2019; Schroder et al., 2004). Moreover, through ESTIMATE analysis, our results showed that the risk score did not significantly correlated with both the immune score and stromal score, which needed to be verified in future studies. We next analyzed the pathways changed between the different groups by GSEA analysis. And the data showed that the endocytosis, chemokine signaling, leukocyte transendothelial migration and VEGF signaling pathways activated in the high-risk group. While fatty acid metabolism, peroxisome, peroxisome proliferator-activated receptors and nitrogen metabolism pathways were enriched in the other group.
Chemical or targeted drugs are important method for HCC treatment. For example, sorafenib and lenvatinib were all recommended as first-line therapies for unresectable hepatocellular carcinoma (Cheng et al., 2009; Kudo et al., 2018). However, MDR seriously affected the treatment efficacy. The MRPs are the well-known transporters, which takes part in multidrug resistance through extruding drugs from cancer cells (Sodani et al., 2012). Here, our results showed that the risk score was positively correlated with MRP1, MRP4, MRP5, and MRP7. Our data indicated that the low-risk group may be more sensitive to chemical or targeted therapy.
Nevertheless, our study has some limitations. Firstly, we only validated the expression of SPP1 in HCC tissues. As all the analyses retrieved is from public sources, further experimental studies are still needed to verify these findings. Secondly, due to the limitation of scRNA-seq technique, a depth coverage of rare cell types was restricted. We could not clearly distinguish neutrophils/macrophage subtypes with known marker genes. Besides, all the included genes were derived from myeloid cell marker genes, while the TME was rather heterogeneous. So, the prognostic predicting value may be inhibited. Last but not least, the underlying connections for these genes between TME and tumor cells deserved to be investigated further.
In conclusion, we constructed a gene signature from HCC specific myeloid cell marker genes and this gene signature was proved to be valuable in immune infiltration analysis, drug resistance prediction and prognostic prediction. The study exhibits valuable knowledge regarding the roles of myeloid cell marker genes in HCC. Our findings may be useful for individualized treatment decisions making.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.
Ethics statement
The studies involving human participants were reviewed and approved by ethics committee from original studies the immunohistochemical staining study was approved by the Ethics Committee of Zhongshan Hospital, Fudan University. The patients/participants provided their written informed consent to participate in this study.
Author contributions
B-HZ, CH, and T-CX contributed to the study design and critical revision of the manuscript. S-SZ and Y-FW carried out the study and drafted the manuscript. S-SZ and Y-FW analyzed the data. All authors read and approved the final manuscript.
Funding
This research was supported by grants from the Medical and Health Project of Xiamen City (Key Project, No. 3502Z20191105); Natural Science Foundation of Fujian Province (grant number 2022J05332); Foundation of Xiamen Branch, Zhongshan Hospital, Fudan University of China (2020ZSXMYS01); Science and Technology Bureau of Xiamen City (3502Z20224ZD1081); Foundation of Zhongshan Hospital, Fudan University of China (2020ZSLC60).
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2023.1118377/full#supplementary-material
References
Ashkar, S., Weber, G. F., Panoutsakopoulou, V., Sanchirico, M. E., Jansson, M., Zawaideh, S., et al. (2000). Eta-1 (osteopontin): An early component of type-1 (cell-mediated) immunity. Science 287 (5454), 860–864. doi:10.1126/science.287.5454.860
Atai, N. A., Bansal, M., Lo, C., Bosman, J., Tigchelaar, W., Bosch, K. S., et al. (2011). Osteopontin is up-regulated and associated with neutrophil and macrophage infiltration in glioblastoma. Immunology 132 (1), 39–48. doi:10.1111/j.1365-2567.2010.03335.x
Bassler, K., Schulte-Schrepping, J., Warnat-Herresthal, S., Aschenbrenner, A. C., and Schultze, J. L. (2019). The myeloid cell compartment-cell by cell. Annu. Rev. Immunol. 37, 269–293. doi:10.1146/annurev-immunol-042718-041728
Butler, A., Hoffman, P., Smibert, P., Papalexi, E., and Satija, R. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36 (5), 411–420. doi:10.1038/nbt.4096
Chen, H., Ye, F., and Guo, G. (2019). Revolutionizing immunology with single-cell rna sequencing. Cell. Mol. Immunol. 16 (3), 242–249. doi:10.1038/s41423-019-0214-4
Chen, H., Zhou, X. H., Li, J. R., Zheng, T. H., Yao, F. B., Gao, B., et al. (2021). Neutrophils: Driving inflammation during the development of hepatocellular carcinoma. Cancer Lett. 522, 22–31. doi:10.1016/j.canlet.2021.09.011
Cheng, A. L., Kang, Y. K., Chen, Z., Tsao, C. J., Qin, S., Kim, J. S., et al. (2009). Efficacy and safety of sorafenib in patients in the asia-pacific region with advanced hepatocellular carcinoma: A phase iii randomised, double-blind, placebo-controlled trial. Lancet Oncol. 10 (1), 25–34. doi:10.1016/S1470-2045(08)70285-7
DeNardo, D. G., and Ruffell, B. (2019). Macrophages as regulators of tumour immunity and immunotherapy. Nat. Rev. Immunol. 19 (6), 369–382. doi:10.1038/s41577-019-0127-6
Deng, L., He, K., Pan, Y., Wang, H., Luo, Y., and Xia, Q. (2021). The role of tumor-associated macrophages in primary hepatocellular carcinoma and its related targeting therapy. Int. J. Med. Sci. 18 (10), 2109–2116. doi:10.7150/ijms.56003
Deng, Y., Cheng, J., Fu, B., Liu, W., Chen, G., Zhang, Q., et al. (2017). Hepatic carcinoma-associated fibroblasts enhance immune suppression by facilitating the generation of myeloid-derived suppressor cells. Oncogene 36 (8), 1090–1101. doi:10.1038/onc.2016.273
Dimri, M., Humphries, A., Laknaur, A., Elattar, S., Lee, T. J., Sharma, A., et al. (2020). Nad(p)h quinone dehydrogenase 1 ablation inhibits activation of the phosphoinositide 3-kinase/akt serine/threonine kinase and mitogen-activated protein kinase/extracellular signal-regulated kinase pathways and blocks metabolic adaptation in hepatocellular carcinoma. Hepatology 71 (2), 549–568. doi:10.1002/hep.30818
Ding, P., Gao, Y., Wang, J., Xiang, H., Zhang, C., Wang, L., et al. (2022). Progress and challenges of multidrug resistance proteins in diseases. Am. J. Cancer Res. 12 (10), 4483–4501.
Engblom, C., Pfirschke, C., and Pittet, M. J. (2016). The role of myeloid cells in cancer therapies. Nat. Rev. Cancer. 16 (7), 447–462. doi:10.1038/nrc.2016.54
Finn, R. S., Qin, S., Ikeda, M., Galle, P. R., Ducreux, M., Kim, T. Y., et al. (2020). Atezolizumab plus bevacizumab in unresectable hepatocellular carcinoma. N. Engl. J. Med. 382 (20), 1894–1905. doi:10.1056/NEJMoa1915745
Fridlender, Z. G., Sun, J., Kim, S., Kapoor, V., Cheng, G., Ling, L., et al. (2009). Polarization of tumor-associated neutrophil phenotype by tgf-beta: "n1" versus "n2" tan. Cancer Cell 16 (3), 183–194. doi:10.1016/j.ccr.2009.06.017
Guan, J., Wang, G., Wang, J., Zhang, Z., Fu, Y., Cheng, L., et al. (2022). Chemical reprogramming of human somatic cells to pluripotent stem cells. Nature 605 (7909), 325–331. doi:10.1038/s41586-022-04593-5
Han, X., Wang, W., He, J., Jiang, L., and Li, X. (2019). Osteopontin as a biomarker for osteosarcoma therapy and prognosis. Oncol. Lett. 17 (3), 2592–2598. doi:10.3892/ol.2019.9905
Kaczocha, M., Vivieca, S., Sun, J., Glaser, S. T., and Deutsch, D. G. (2012). Fatty acid-binding proteins transport n-acylethanolamines to nuclear receptors and are targets of endocannabinoid transport inhibitors. J. Biol. Chem. 287 (5), 3415–3424. doi:10.1074/jbc.M111.304907
Kashani-Sabet, M., Nosrati, M., Miller, J. R., Sagebiel, R. W., Leong, S., Lesniak, A., et al. (2017). Prospective validation of molecular prognostic markers in cutaneous melanoma: A correlative analysis of e1690. Clin. Cancer Res. 23 (22), 6888–6892. doi:10.1158/1078-0432.CCR-17-1317
Korsunsky, I., Millard, N., Fan, J., Slowikowski, K., Zhang, F., Wei, K., et al. (2019). Fast, sensitive and accurate integration of single-cell data with harmony. Nat. Methods. 16 (12), 1289–1296. doi:10.1038/s41592-019-0619-0
Kudo, M., Finn, R. S., Qin, S., Han, K. H., Ikeda, K., Piscaglia, F., et al. (2018). Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: A randomised phase 3 non-inferiority trial. Lancet 391 (10126), 1163–1173. doi:10.1016/S0140-6736(18)30207-1
Lee, C., Kerouanton, B., Chothani, S., Zhang, S., Chen, Y., Mantri, C. K., et al. (2021). Coding and non-coding roles of mocci (c15orf48) coordinate to regulate host inflammation and immunity. Nat. Commun. 12 (1), 2130. doi:10.1038/s41467-021-22397-5
Lee, H., Oh, E. T., Choi, B. H., Park, M. T., Lee, J. K., Lee, J. S., et al. (2015). Nqo1-induced activation of ampk contributes to cancer cell death by oxygen-glucose deprivation. Sci. Rep. 5, 7769. doi:10.1038/srep07769
Lin, L., Sun, J., Tan, Y., Li, Z., Kong, F., Shen, Y., et al. (2017). Prognostic implication of nqo1 overexpression in hepatocellular carcinoma. Hum. Pathol. 69, 31–37. doi:10.1016/j.humpath.2017.09.002
Liu, G., Friggeri, A., Yang, Y., Park, Y. J., Tsuruta, Y., and Abraham, E. (2009). Mir-147, a microrna that is induced upon toll-like receptor stimulation, regulates murine macrophage inflammatory responses. Proc. Natl. Acad. Sci. U. S. A. 106 (37), 15819–15824. doi:10.1073/pnas.0901216106
Liu, J., Sun, B., Guo, K., Yang, Z., Zhao, Y., Gao, M., et al. (2022). Lipid-related fabp5 activation of tumor-associated monocytes fosters immune privilege via pd-l1 expression on treg cells in hepatocellular carcinoma. Cancer Gene Ther. 29, 1951–1960. doi:10.1038/s41417-022-00510-0
Liu, L., Zhang, R., Deng, J., Dai, X., Zhu, X., Fu, Q., et al. (2022). Construction of tme and identification of crosstalk between malignant cells and macrophages by spp1 in hepatocellular carcinoma. Cancer Immunol. Immunother. 71 (1), 121–136. doi:10.1007/s00262-021-02967-8
Liu, Z., Zhang, D., Liu, C., Li, G., Chen, H., Ling, H., et al. (2021). Comprehensive analysis of myeloid signature genes in head and neck squamous cell carcinoma to predict the prognosis and immune infiltration. Front. Immunol. 12, 659184. doi:10.3389/fimmu.2021.659184
Liu, Z., Zhang, S., Ouyang, J., Wu, D., Chen, L., Zhou, W., et al. (2022). Single-cell rna-seq analysis reveals dysregulated cell-cell interactions in a tumor microenvironment related to hcc development. Dis. Markers. 2022, 4971621. doi:10.1155/2022/4971621
Locati, M., Curtale, G., and Mantovani, A. (2020). Diversity, mechanisms, and significance of macrophage plasticity. Annu. Rev. Pathol. 15, 123–147. doi:10.1146/annurev-pathmechdis-012418-012718
Losic, B., Craig, A. J., Villacorta-Martin, C., Martins-Filho, S. N., Akers, N., Chen, X., et al. (2020). Intratumoral heterogeneity and clonal evolution in liver cancer. Nat. Commun. 11 (1), 291. doi:10.1038/s41467-019-14050-z
Mabbott, N. A., Baillie, J. K., Brown, H., Freeman, T. C., and Hume, D. A. (2013). An expression atlas of human primary cells: Inference of gene function from coexpression networks. BMC Genomics 14, 632. doi:10.1186/1471-2164-14-632
Nakamura, K., and Smyth, M. J. (2020). Myeloid immunosuppression and immune checkpoints in the tumor microenvironment. Cell. Mol. Immunol. 17 (1), 1–12. doi:10.1038/s41423-019-0306-1
Newman, A. M., Steen, C. B., Liu, C. L., Gentles, A. J., Chaudhuri, A. A., Scherer, F., et al. (2019). Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat. Biotechnol. 37 (7), 773–782. doi:10.1038/s41587-019-0114-2
Pan, D., Yang, W., Zeng, Y., Qin, H., Xu, Y., Gui, Y., et al. (2022). Akr1c3 regulated by nrf2/mafg complex promotes proliferation via stabilizing parp1 in hepatocellular carcinoma. Oncogene 41 (31), 3846–3858. doi:10.1038/s41388-022-02379-7
Peng, Y., Yuan, J., Zhang, Z., and Chang, X. (2017). Cytoplasmic poly(a)-binding protein 1 (pabpc1) interacts with the rna-binding protein hnrnpll and thereby regulates immunoglobulin secretion in plasma cells. J. Biol. Chem. 292 (29), 12285–12295. doi:10.1074/jbc.M117.794834
Penning, T. M. (2015). The aldo-keto reductases (akrs): Overview. Chem. Biol. Interact. 234, 236–246. doi:10.1016/j.cbi.2014.09.024
Petermann, F., Pekowska, A., Johnson, C. A., Jankovic, D., Shih, H. Y., Jiang, K., et al. (2019). The magnitude of ifn-gamma responses is fine-tuned by dna architecture and the non-coding transcript of ifng-as1. Mol. Cell. 75 (6), 1229–1242. doi:10.1016/j.molcel.2019.06.025
Roskoski, R. J. (2022). Properties of fda-approved small molecule protein kinase inhibitors: A 2022 update. Pharmacol. Res. 175, 106037. doi:10.1016/j.phrs.2021.106037
Rothwell, S. W., Nath, J., and Wright, D. G. (1993). Rapid and reversible tubulin tyrosination in human neutrophils stimulated by the chemotactic peptide, fmet-leu-phe. J. Cell. Physiol. 154 (3), 582–592. doi:10.1002/jcp.1041540317
Schroder, K., Hertzog, P. J., Ravasi, T., and Hume, D. A. (2004). Interferon-gamma: An overview of signals, mechanisms and functions. J. Leukoc. Biol. 75 (2), 163–189. doi:10.1189/jlb.0603252
Sharafi, F., Hasani, S. A., Alesaeidi, S., Kahrizi, M. S., Adili, A., Ghoreishizadeh, S., et al. (2022). A comprehensive review about the utilization of immune checkpoint inhibitors and combination therapy in hepatocellular carcinoma: An updated review. Cancer Cell Int. 22 (1), 269. doi:10.1186/s12935-022-02682-z
Sodani, K., Patel, A., Kathawala, R. J., and Chen, Z. S. (2012). Multidrug resistance associated proteins in multidrug resistance. Chin. J. Cancer 31 (2), 58–72. doi:10.5732/cjc.011.10329
Sun, Y., Wu, L., Zhong, Y., Zhou, K., Hou, Y., Wang, Z., et al. (2021). Single-cell landscape of the ecosystem in early-relapse hepatocellular carcinoma. Cell 184 (2), 404–421. e16. doi:10.1016/j.cell.2020.11.041
Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71 (3), 209–249. doi:10.3322/caac.21660
Tamborero, D., Rubio-Perez, C., Muinos, F., Sabarinathan, R., Piulats, J. M., Muntasell, A., et al. (2018). A pan-cancer landscape of interactions between solid tumors and infiltrating immune cell populations. Clin. Cancer Res. 24 (15), 3717–3728. doi:10.1158/1078-0432.CCR-17-3509
Tang, Y., Li, K., Hu, B., Cai, Z., Li, J., Tao, H., et al. (2022). Fatty acid binding protein 5 promotes the proliferation, migration, and invasion of hepatocellular carcinoma cells by degradation of kruppel-like factor 9 mediated by mir-889-5p via camp-response element binding protein. Cancer Biol. Ther. 23 (1), 424–438. doi:10.1080/15384047.2022.2094670
Tibshirani, R. (1997). The lasso method for variable selection in the cox model. Stat. Med. 16 (4), 385–395. doi:10.1002/(sici)1097-0258(19970228)16:4<385::aid-sim380>3.0.co;2-3
Tritschler, F., Huntzinger, E., and Izaurralde, E. (2010). Role of gw182 proteins and pabpc1 in the mirna pathway: A sense of deja vu. Nat. Rev. Mol. Cell Biol. 11 (5), 379–384. doi:10.1038/nrm2885
Wang, J., Hao, F., Fei, X., and Chen, Y. (2019). Spp1 functions as an enhancer of cell growth in hepatocellular carcinoma targeted by mir-181c. Am. J. Transl. Res. 11 (11), 6924–6937.
Wu, J. C., Sun, B. S., Ren, N., Ye, Q. H., and Qin, L. X. (2010). Genomic aberrations in hepatocellular carcinoma related to osteopontin expression detected by array-cgh. J. Cancer Res. Clin. Oncol. 136 (4), 595–601. doi:10.1007/s00432-009-0695-0
Xue, T. C., Jia, Q. A., Ge, N. L., Chen, Y., Zhang, B. H., and Ye, S. L. (2015). Imbalance in systemic inflammation and immune response following transarterial chemoembolization potentially increases metastatic risk in huge hepatocellular carcinoma. Tumour Biol. 36 (11), 8797–8803. doi:10.1007/s13277-015-3632-7
Yoshihara, K., Shahmoradgoli, M., Martinez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612
Zhang, H., Sheng, C., Yin, Y., Wen, S., Yang, G., Cheng, Z., et al. (2015). Pabpc1 interacts with ago2 and is responsible for the microrna mediated gene silencing in high grade hepatocellular carcinoma. Cancer Lett. 367 (1), 49–57. doi:10.1016/j.canlet.2015.07.010
Zhang, Q., He, Y., Luo, N., Patel, S. J., Han, Y., Gao, R., et al. (2019). Landscape and dynamics of single immune cells in hepatocellular carcinoma. Cell 179 (4), 829–845. doi:10.1016/j.cell.2019.10.003
Zhang, T., Guan, G., Zhang, J., Zheng, H., Li, D., Wang, W., et al. (2022). E2f1-mediated auf1 upregulation promotes hcc development and enhances drug resistance via stabilization of akr1b10. Cancer Sci. 113 (4), 1154–1167. doi:10.1111/cas.15272
Keywords: tumor associated myeloid cells, single-cell RNA sequencing, hepatocellular carcinoma, prognosis, immune infiltration, drug resistance
Citation: Zheng S-S, Wu Y-F, Zhang B-H, Huang C and Xue T-C (2023) A novel myeloid cell marker genes related signature can indicate immune infiltration and predict prognosis of hepatocellular carcinoma: Integrated analysis of bulk and single-cell RNA sequencing. Front. Mol. Biosci. 10:1118377. doi: 10.3389/fmolb.2023.1118377
Received: 07 December 2022; Accepted: 27 February 2023;
Published: 07 March 2023.
Edited by:
Yunpeng Hua, The First Affiliated Hospital of Sun Yat-sen University, ChinaReviewed by:
Yong Jiang, East Tennessee State University, United StatesXiaoyun Bu, Sun Yat-sen University Cancer Center (SYSUCC), China
Copyright © 2023 Zheng, Wu, Zhang, Huang and Xue. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Cheng Huang, huang.cheng@zs-hospital.sh.cn; Tong-Chun Xue, xue.tongchun@zs-hospital.sh.cn
†These authors have contributed equally to this work and share first authorship