Tumor Immune Microenvironment Characterization in Hepatocellular Carcinoma Identifies Four Prognostic and Immunotherapeutically Relevant Subclasses

Purpose: The tumor microenvironment (TME) plays a critical role in the pathogenesis of hepatocellular carcinoma (HCC). However, underlying compositions and functions that drive the establishment and maintenance of the TME classifications are less-well understood.

Methods: A total of 766 HCC patients from three public cohorts were clustered into four immune-related subclasses based on 13 TME signatures (11 immune-related cells and 2 immune-related pathways) calculated by MCP-counter. After analyzing the landscapes of functional annotation, methylation, somatic mutation, and clinical characteristics, we built a TME-based Support Vector Machine of 365 patients (discovery phase) and 401 patients (validation phase). We applied this SVM model on another two independent cohorts of patients who received sorafenib/pembrolizumab treatment.

Results: About 33% of patients displayed an immune desert pattern. The other subclasses were different in abundance of tumor infiltrating cells. The Immunogenic subclass (17%) associated with the best prognosis presented a massive T cell infiltration and an activation of immune checkpoint pathway. The 13 TME signatures showed a good potential to predict the TME classification (average AUC = 88%). Molecular characteristics of immunohistochemistry from Zhejiang cohort supported our SVM classification. The optimum response to pembrolizumab (78%) and sorafenib (81%) was observed in patients belonging to the Immunogenic subclass.

Conclusions: The HCC patients from distinct immune subclass showed significant differences in clinical prognosis and response to personalized treatment. Based on tumor transcriptome data, our workflow can help to predict the clinical outcomes and to find appropriate treatment strategies for HCC patients.

Introduction

Hepatocellular carcinoma (HCC), the predominant type of primary liver cancer, is the fourth leading cause of cancer mortality worldwide with about 782,000 deaths annually (1). HCC is strongly influenced by the tumor microenvironment (TME), as well as reported to benefit from the immune-checkpoint blockade treatment. As an inflammation-related tumor, the specific TME of HCC can influence the immune tolerance and evasion by mixed mechanisms. Numerous studies have been reported that the TME plays a critical role in tumor initiation, progression, and outcome (2, 3). The TME is an intricate system, which coexists and interacts with cancer cells, immune cell subsets, extracellular matrix, various cytokine, and other unknown components to maintain the tumorigenesis of HCC (4). The complexity of the TME relies on the immune infiltration, as tumors could be classified into the tumor immunity continuum (5). Hegde et al. suggested that human tumors could be categorized as inflamed, immune desert, or immune-excluded phenotypes correspond to different mechanisms of immune response and escape (6). Thorsson et al. proposed a classification of six immune subclasses (wound healing, lymphocyte depleted, TGF-β dominant, inflammatory, immunologically quiet, IFN-γ dominant), based on extensive immunogenomic analysis of 33 cancer types compiled by TCGA (7).

Meanwhile, estimating the cellular composition of the TME requires accurate and robust methods. Fluorescence-activated cell sorting (FACS) operates only a small number of cell type-specific markers and requires large amounts of fresh tumor tissues, which limit the applications on tumor biopsies (8). Single-cell sequencing has a high precision, but currently it is too expensive for large-scale clinical application (9). In order to overcome the above shortcomings, we turn to the methods for high-throughput technologies applied in clinical settings, which TME are inferred using computational algorithms. High throughput technologies, such as RNA-Seq and microarray, provide large-scaled transcriptome data and offer opportunity for estimation of the abundance of tumor infiltrated immune cells. Several methods like Microenvironment Cell Population-counter (MCP-counter), CIBERSORT and TIMER have been developed to robustly and precisely quantify immune cells using transcriptome data obtained from bulk tissue specimens (10–12).

After exploring the distinct compositions and functions of the TME by MCP-counter, a total of 766 HCC patients from three public cohorts were clustered into four subclasses (namely Immune desert, Immunogenic, Innate immune and Mesenchymal) based on 13 TME signatures. Furthermore, a Support Vector Machine (SVM) was constructed to predict the HCC classification (average AUC = 88%). Finally, by applying our SVM model on another two independent cohorts of patients who received sorafenib/pembrolizumab treatment, we found that patients classified into Immunogenic subclass showed the highest response rate to sorafenib (81%) and pembrolizumab (78%). Thus, we suggested that HCC patients may benefit from identifying the immune subclass which infer clinical outcomes and guide personalized treatment strategies (Supplementary Figure 1).

Materials and Methods

Ethics Statement and Consent for Publication

The studies involving human participants were reviewed and approved by Clinical Research Ethics Committee of the First Affiliated Hospital College of Medicine, Zhejiang University (2014-334). Operation informed consents and Informed consent form for scientific research were obtained from all participants for the publication of any potentially identifiable images or data included in this article.

Clinical Cohorts and Preprocessing

Three public transcriptome data sets were enrolled in our study, including the TCGA-LIHC cohort of the Cancer Genome Atlas (n = 365), the CHCC-HBV cohort of Gao et al. (n = 159) and the GSE14520 cohort of Roessler et al. (n = 242) (13, 14). Among above data sets, any case with null value of survival information had been excluded. As to TCGA-LIHC and CHCC-HBV cohort, the fragments per kilobase per million (FPKM) data and clinic information were downloaded from the UCSC Xena (xenabrowser.net) and The National Omics Data Encycolpedia (www.biosino.org/node), respectively. The “CEL” files of GSE14520 were downloaded and normalized by the “frma” function using frma (R package). Besides, GSE78220 (28 melanoma patients received pembrolizumab) and GSE109211 (67 HCC patients received sorafenib) were also included in our study (15, 16). All FPKM values were transformed into transcripts per kilobase million (TPM) values. Raw data files of GSE109211 were normalized by lumi (R package). cBioPortal (www.cbioportal.org) was used to download the beta value of DNA methylation status of checkpoint genes from TCGA-LIHC cohort. Gene mutation data (MAF files) of TCGA-LIHC cohort was achieved from TCGA database.

In our previous study, liver tumor tissues from 32 patients in First Affiliated Hospital, School of Medicine, Zhejiang University were collected from November 2013 to July 2014 (GSE138485/PRJNA576155) (17). Only a few of the 32 patients have the available Formalin Fixed Paraffin Embedded (FFPE) specimens. Detailed information of patients is described in Supplementary Tables S7, S8.

Quantification of TME Infiltration

The abundances of immune and stromal cells in TME were quantified by MCP-counter based on cell-type specific transcriptome signatures (10). According to Sylvie’s study (18), a total of 13 TME signatures, which contained 11 stromal and immune cell populations (Lymphoid, B_derived, T_adaptive, Cytotoxic, Monocyte_derived, Myeloid, NK_or_T, Fibroblast, HSCactivated, HSCquiescent, and Myofibroblast) and two functional signatures representing the immune checkpoints (named Checkpoint) and the immunosuppression pathways (named Immunosuppression), were included (Supplementary Table S1). In addition, we used the CIBERSORT to validate the immune characterization (11).

Evaluation of the Immune Score and Stromal Score

The immune score and stromal score of each patient were calculated by ESTIMATE algorithm based on transcriptomic data (19). The R code of ESTIMATE was downloaded from the public source website (https://sourceforge.net/projects/estimateproject).

Unsupervised Clustering Based on 13 TME Signatures

Based on above 13 signatures, consensus clustering method was used to classify HCC patients into distinct immune subclasses by ConsensusClusterPlus (R package) (20). Detailed settings were as followed: repetitions = 500 times; pItem = 0.8; pFeature = 0.8. The number of the clusters was determined by consensus cumulative distribution function (CDF) curve and the delta area (relative change in CDF area). Because the CDF curve and delta area plot showed that the delta area increased slightly for k = 5 compared to k = 4, we finally selected k = 4 (four immune subclasses) as the best solution (Supplementary Figure 2).

Functional Characterization of Immune Subclass

For pathway analysis of transcriptomic data among immune subclasses, we performed Gene Set Variation Analysis (GSVA) analysis with the “GSVA” R package (21). The gene sets of “c2.cp.kegg.v7.1.symbols”, “c2.cp.bp.v7.1.symbols”, “c2.cp.biocarta.v7.1.symbols”, and “c2.cp.pid.v7.1.symbols” were downloaded from Molecular Signatures Database (MSigDB). Adjusted ANOVA model q value <0.05 was considered as statistically significance.

Classifier Model Construction and Validation

Based on the 13 TME signatures, we developed a classifier model using Support Vector Machine, as implemented in python package “scikit-learn” (version 0.21.3). The TCGA-LIHC cohort was used as discovery phase, as well as the CHCC-HBV and GSE14520 cohort were validation phases. Detailed information for model construction is described in Supplementary Methods. The efficiency of the classifier model was evaluated by receiver operating characteristic (ROC) curve and the area under the curve (AUC).

Statistical Analysis

Continuous variables conforming to normal distribution were compared with Student t test, otherwise the Wilcoxon rank sum test was used. One-way ANOVA models and Kruskal-Wallis tests were used for multigroup comparison. The association between immune subclasses and the clinical parameters were evaluated by chi-squared test or Fisher-exact test. Overall survival (OS) curves were calculated according to the Kaplan-Meier method (R package survival) and differences between curves were assessed using the log-rank test. Statistical analyses were performed on R 3.6.2 software and SPSS V26.0 for Windows.

Results

Identified Four Immune Subclasses Based on TME of HCC

Three public HCC data sets with clinical information (TCGA, CHCC-HBV, GSE14520) were enrolled in this study. According to Sylvie’s study (18), a total of 13 TME signatures represented the major infiltrated cell composition and several components of tumor-stroma interaction were included (Supplementary Table S1). First, we performed spearman correlation analysis on the 11 cell-type signatures to find the interdependent relationship. Figure 1A and Supplementary Figure 3A showed that these 11 cell-type signatures were clustered into three distinct clusters (named ACTIVATED_FIBROBLASTS, INNATE_IMMUNITY, and ADAPTIVE_IMMUNITY). In addition, two functional pathways were added, namely, a signature of immune checkpoint related to immune therapy and a signature of genes involved in immunosuppression. Then consensus clustering was performed on the three data sets based on above 13 signatures, and four distinct immune patterns were finally identified, named C1 to C4 (Figure 1B, Supplementary Figure 3B). Subclass C1 showed an immune desert pattern distinguished by low abundances of all the TME signatures. Subclass C2 displayed an immunogenic pattern distinguished by high abundances of both INNATE_IMMUNITY and ADAPTIVE_IMMUNITY, activation of immune checkpoint pathway and low abundances of ACTIVATED_FIBROBLASTS. Subclass C3 showed an innate immune pattern distinguished by moderate to high abundances of INNATE_IMMUNITY and immunosuppression pathway, rather low abundances of ADAPTIVE_IMMUNITY. Subclass C4 displayed a mesenchymal pattern which was characterized by high abundances of ACTIVATED_FIBROBLASTS and immunosuppression pathway. Principal component analysis (PCA) also showed a significant spatial separation among these four subclasses in TCGA cohort (Figure 1C).

FIGURE 1

Figure 1 Classification based on tumor microenvironment stratifies HCCs into four subclasses. (A) Correlation heatmap of 11 TME cell signatures in two data sets. Color scale: Spearman correlation coefficient from 0 (blue) to 1 (orange). (B) Consensus clustering analysis of two data sets revealed four HCC subclasses based on 13 TME signatures. Color scale: Z score from -2 (blue) to +2 (orange). (C) Principal-component analysis based on 13 TME signatures separated different subclasses in TCGA cohort. Kaplan-Meier curves of overall survival for TCGA cohort (D) and CHCC-HBV cohort (E) based on immune subclasses (log-rank test).

Next, we investigated the association between clinical outcomes and immune subclass. In all cohorts, significant differences of prognosis were existed among four immune subclasses, indicating that they could be clinically relevant subclasses (TCGA cohort: log-rank test P = 0.0005; CHCC-HBV cohort: log-rank test P = 0.0052; GSE14520 cohort: log-rank test P = 0.0073) (Figures 1D, E, Supplementary Figure 3C). In TCGA cohort, C2 Immunogenic subclass showed the longest median survival time (MST) (MST = 82.8 months), followed by C1 Immune desert subclass (MST = 71.0 months), thirdly C4 Mesenchymal subclass (MST = 52.0 months), lastly C3 Innate immune subclass (MST = 21.3 months). Similar results were found in other two data sets. In summary, subclass C2 showed a survival advantage with respect to the other subclasses.

Immune Functional Characteristics of the Immune Subclasses

To refine the immune characterization, we performed both ESTIMATE and CIBERSORT to calculate the Immune/Stromal scores and the proportion of 22 tumor infiltrating immune cells. C1 Immune desert subclass showed the lower Immune score and Stromal score, while C2 Immunogenic subclass and C4 Mesenchymal subclass showed the highest Immune score and Stromal score, respectively (Supplementary Figure 4). CIBERSORT analysis revealed significant difference in 16 out of 22 tumor-infiltrating immune cells, especially an enrichment of CD4⁺ and CD8⁺ T cells in C2, as well as an enrichment of M2 macrophages in C3 (Supplementary Table S2). The results of ESTIMATE and CIBERSORT further confirmed the characteristics of the four subclasses defined by MCP-counter.

To explore the functional differences among these four subclasses, we performed GSVA enrichment analysis on TCGA cohort (Figure 2A, Supplementary Table S3). C1 Immune desert subclass showed a highly attenuation of stromal and immune pathways. C2 Immunogenic subclass showed an enrichment of immune response, such as major histocompatibility complex (MHC) class I and class II biosynthesis, B cell mediated immunity and chemotaxis, T cell cytotoxicity, CD8⁺ T cell activation and differentiation, T cell survival, and immune checkpoint pathway (CTLA4 and PD-1). C3 Innate immune subclass was enriched in macrophage activation, M2 macrophage polarization, TLR3 and LPS signaling. C3 subclass also showed an enrichment of T cell activation and cytotoxicity pathway, but not of T cell survival (different with C2 subclass). These results may account for the lack of adaptive immunity in C3 subclass. C4 subclass was remarkably enriched in activated HSC and stromal pathways such as ECM assembly, EMT, angiogenesis, TGF beta, integrin signaling pathway. The expressions of genes belonging to several immune pathways confirmed the differences among the four subclasses (Figure 2B). Markers of macrophage chemotaxis were increased in both C2 and C3, while C3 showed higher expression of macrophage activation. Moreover, the increases in markers of T cell survival, T cell chemotaxis and activation were observed in C2. Markers of all the pathways were low expressed in C1. Specially, compared to the other subclass, C2 subclass showed the both higher expressions and methylation levels of checkpoint-related genes obtained from bulk tumor tissues, which indicated that the overexpression of checkpoints in C2 subclass was potentially triggered by hypomethylation of these genes (Figures 2C, D) (22).

FIGURE 2

Figure 2 Functional characteristics of four immune subclasses. (A) Heatmap of GSVA scores for indicated functional signatures. Color scale: GSVA score from −1 (blue) to +1(red). (B) Boxplot plot of the expression levels for selected immune-related pathways. (One-way ANOVA test). Boxplot of the expression (C) and methylation (D) levels of immune checkpoint-related genes between C2 and non-C2 subclass. (Student’s t test). All P values labels: ns P > 0.05, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001). Error bars are presented as the standard deviation (SD).

Mutational Landscape of the Immune Subclasses

Somatic alterations have been proven to be correlated with TME (23). We analyzed somatic mutation data from the whole tumor of several genes with high frequency of mutation and in specific pathways, such as P53-pathway, Wnt-pathway, Chromatin modifiers pathway, and hepatic differentiation from TCGA-LIHC cohort (Figure 3A, Supplementary Table S4) (24). C1 subclass showed the highest mutation frequency of CTNNB1 (P < 0.0001), while the highest mutation frequency of ARID2 (P = 0.049) was observed in C2 subclass (Figure 3B, Supplementary Table S4). The highest mutation frequency of TP53 was observed in C3 subclass (P = 0.015). Then, we compared the tumor mutation burden and predicted neoantigens among these four subclasses (Figures 3C, D). The lowest tumor mutation burden and numbers of predicted neoantigens were detected in C4 subclass. Moreover, in CHCC-HBV cohort, Gao et al. found that signature for aristolochic acids (AA signature) was correlated to tumor mutation burden and response to immunotheropy. In CHCC-HBV cohort, C4 subclass also showed a lower proportion of AA signature than other three subclasses, which indicated a lower benefit from checkpoint blockade therapy (Figure 3E).

FIGURE 3

Figure 3 Differences in the mutational landscape among distinct immune subclasses. (A) Oncoplot of tumor somatic mutation of genes in P53 pathway, Wnt/beta-catenin pathway, Chromatin modifiers pathway and hepatic differentiation based on TCGA cohort. (B) Comparisons of the frequently mutated genes among four immune subclasses based on TCGA cohort. (Fisher’s exact test). Comparison of tumor mutation burden (C) and predicted neoantigens (D) among four immune subclasses. (Wilcoxon rank sum test). (E) The proportion of patients with AA signature among distinct subclasses in CHCC-HBV cohort. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Associations Between Immune Subclass and Clinical Characteristic in TCGA and CHCC-HBV Cohort

Next, we discovered the associations between clinical characteristics and immune subclass in TCGA and CHCC-HBV cohort (Figure 4, Supplementary Table S5, S6). The patients from C2 subclass had a higher proportion of pathologic stage I/II (TCGA: P = 0.002, CHCC-HBV: P = 0.011). Patients from C3 subclass showed a higher proportion of HBV infection (P = 0.048). Most of the clinicopathological characteristics did not show significant differences, which suggested that the main factors distinguishing distinct immune subclass were the TME signatures, rather than the above-mentioned clinical features.

FIGURE 4

Figure 4 Clinical characteristics of four immune subclasses. Correlation of the immune subclass with clinical characteristics and previously reported HCC classification in TCGA cohort (A) and CHCC-HBV cohort (B).

Furthermore, we compared our classification with several previous reported classifications based on transcriptomic, including Lee’s classification (High/Low survival), Boyault’s classification (G1 to G6), Chiang’s classification (five classes), Hoshida’s classification (S1 to S3), and Lachenmayer’s classification (CTNNB1 class/Wnt-TGF-beta class) (25–29). In TCGA and CHCC-HBV cohort, C1 subclass was co-clustered with the better-prognosis subclasses (Lee’s High survival, Boyault’s G5&6, Chiang’s CTNNB1 class, Hoshida’s S3, Lachenmayer’s CTNNB1-class). C3 subclass was largely co-clustered with poor-prognosis subclasses (Lee’s Low survival, Boyault’s G3, Chiang’s Proliferation class, Hoshida’s S1, Lachenmayer’s Wnt-TGF-beta class). C2 subclass was linked to both better-prognosis subclass (Lee’s High survival, Hoshida’s S3) and poor-prognosis subclasses (Boyault’s G1&G2, Chiang’s Proliferation, Wnt-TGF beta class). C4 subclass was co-clustered with Lee’s Low survival, Boyault’s G1&G2, Chiang’s Proliferation and Interferon class, Hoshida’s S3, Lachenmayer’s Wnt-TGF-beta class.

Construction and Validation of a Classifier Based on TME

Characteristics of the clinical traits and biological behaviors among four immune subclasses supported our classification. To apply this classification on clinical use, we developed a Support Vector Machine model to classify HCC patients into above four immune subclasses. The input was the 13 TME signatures of each case from the above data sets, and the output was the immune subclass of each case calculated by this model. The ROC curve represents the accuracy between the subclass clustered in Figure 1B and Supplementary Figure 3B and the subclass predicted by this SVM model. As Figures 5A–C showed, these 13 TME signatures revealed a great classification performance in both discovery phase (TCGA cohort: AUC = 0.98) and validation phases (CHCC-HBV cohort: AUC = 0.91; GSE14520 cohort: AUC = 0.85). Furthermore, we applied this classifier model on 32 HBV-related HCC patients (Zhejiang cohort) to divide these patients into four subclasses. We selected a random sample from each subclass to perform immunohistochemical staining for verifying the accuracy of our classifier model (Figure 5D, Supplementary Figure 5). Several markers of immune and stromal cells were selected, specifically, CD4 and CD8 for T-lymphocytes, CD20 for B-lymphocytes, CD68 for macrophages, αSMA for fibroblastic cells and Vimentin for mesenchymal cells. These markers varied markedly among four subclasses. The expression of Vimentin and αSMA was low in the patient classified into subclass C1, moderate in subclass C2 and C3, high in subclass C4. The patient classified into subclass C2 was characterized by the massive infiltration CD4⁺, CD8⁺ and CD20⁺ lymphocytes. Innate immune cells (macrophages) were also observed in subclass C2. Subclass C3 displayed a high infiltration of macrophages. Subclass C4 contained a low density of macrophages and CD4⁺ T cells. Thus, based on immunohistochemistry, the phenotypic features of HCC tumors were consistent with the classification of our SVM model. The results of immunohistochemistry not only partially proved the accuracy of our model, but also supported the rationality of the classification of four immune-related subclasses.

FIGURE 5

Figure 5 Construction of support vector machine model and performance validation. ROC curves for classifiers designed to predict the immune subclass for TCGA (A), CHCC-HBV (B), and GSE14520 (C). (D) Representative immunohistochemical pictures of HCC samples belonging to each subclass (100X) (Zhejiang cohort).

Different Sensitivity to Personalized Treatment Among Four Immune Subclasses

In HCC patients with Child-Pugh Class A or B, the multi-kinase inhibitor sorafenib has become the first-line systemic therapy (30). However, there were still no effective clinical characteristics to predict the response to sorafenib so far. Several recent studies suggested that sorafenib may exert the anti-tumor effect by regulating the TME of HCC (31, 32). By our SVM classifier, 67 patients from GSE109211 were divided into four subclasses to explore the associations between the immune subclass and the response to sorafenib. We found that 81% cases of C2 subclass showed a significant response to sorafenib, indicating that patients from Immunogenic subclass were more likely to benefit from sorafenib treatment (Figure 6A). As Pinyol and colleagues divided these 67 patients into Good/Poor Prognosis subgroups, C2 subclass was also coclustered with the Good Prognosis subgroup (95%) (Figure 6B).

FIGURE 6

Figure 6 The role of immune subclass in personalized treatment and schematic summary of each immune subclass. (A) The number of patients with response to sorafenib. (B) The number of patients belonged to good or poor prognosis subgroup. (C) The number of patients with response to anti-PD-1 therapy. (D) Kaplan-Meier curves of overall survival for GSE78220 cohort (log-rank test, P = 0.046). (E) The mechanisms of immune escape for each immune subclass. (F) Potential therapeutic strategies for each immune subclass.

In recent years, some immunotherapies like PD-1 blockade have achieved success in HCC (33). Different immune cell infiltrations and expressions of checkpoint-related genes suggested that four immune subclasses could have the distinct response to immunotherapy (34). We tried to apply our model on another pembrolizumab-treated cohort (GSE78220). The highest response rate (77.8%) to pembrolizumab was observed in patients belonging to C2 subclass with the best outcome (Figures 6C, D). The results indicated that HCC patients belonging to immunogenic subclass may benefit from anti-PD-1 therapy inferentially.

Discussion

Although there is a strong heterogeneity in the tumor immune microenvironment of each HCC patient, a clinical benefit could be made from classifying a patient into a specific immune subclass. After analyzing the landscapes of transcriptome, methylation, somatic mutation, and clinical characteristics, we found that these four subclasses may correspond to different mechanisms of immune escape (Figure 6E). Immune desert subclass (C1) is characterized by immune ignorance and a lack of priming T cell, corresponding to immune-desert phenotype. The activation of the β-catenin caused by CTNNB1 mutation might account for the low immune infiltration represented in C1 (35). Immunogenic subclass (C2) is characterized by a massive immune cell (CD4⁺ T cell, CD8⁺ T cell, B cell, and macrophage) infiltration in tumor corresponding to immune-inflamed phenotype. Negative regulators of the immune response (PD-L1, CTLA4, etc.) might be involved in counteraction of anti-tumor immune response (36). A mount of patients belonging to C2 subclass showed a low pathological stage (TNM stage). In line with our study, several studies demonstrated that the tumors in low pathologic stage usually infiltrated with numerous immune cells. The patients belonging to C2 subclass showed the highest mutation frequency of ARID2 which was related to the efficacy of checkpoint blockade immunotherapy in clear cell renal cell carcinoma (37). Innate immune subclass (C3) is characterized by the activation of M2 macrophages (related to innate immunity). M2 macrophage, which exerts the anti-inflammatory and immunosuppressive effects, might promote the immune escape represented in C3 subclass through inhibiting the infiltration of adaptive immune cells (38). Mesenchymal subclass (C4) shows a large number of activated fibroblasts (including HSC and myofibroblast) which influence EMT and the sensitivity of drug treatment through synthesizing growth factors, chemokines and adhesion molecules (39).

Kaplan-Meier analysis based on 766 participants showed that significant differences in overall survival were discovered to exist among our four immune subclasses. Immunogenic subclass (C2) represented the best clinical outcome, while innate immune subclass (C3) have the worst. This suggested that the different TME continuously and chronically affects the progression of HCC, which is ultimately reflected in the different clinical outcome. Based on the RNA-seq data from bulk tumor tissues, our convenient classification dividing the patients into four subclasses may infer the prognosis. What is more, the conversions of the immune subclasses by external interventions may benefit the long-term clinical results and outcomes of HCC patients.

Additionally, we established an SVM model based on the 13 TME signatures (11 immune-related cells and 2 immune-related pathways) and confirmed its predictive value (CHCC-HBV cohort: AUC = 0.91; GSE14520 cohort: AUC = 0.85). The input (13 TME signatures) of the SVM model were calculated by MCP-counter based on transcriptomic data, while giving a specific output (which immune subclass). Thanks to the wide applications of RNA-seq, our workflow only required frozen/fresh tissue samples (< 100 mg), as well the model constructed on python was convenient and efficient. The immunohistochemistry from each subclass proved not only the rationality of the TME classification but also the accuracy of the SVM model (Figure 5D). Accordingly, this suggested that for any HCC patient undergoing liver biopsy or liver resection, our SVM model can be used to infer the prognosis and guide the follow-up treatment.

HCC patients may benefit from identifying immune subclass which may guide personalized treatment strategies (Figure 6F). In our study, we found that the patients belonging to C2 subclass might be more suitable for sorafenib and anti-PD-1 therapy. According to the reported study, the patients belonging to C3 subclass could be treated with colony-stimulating factor-1 inhibitor which improved the efficacy of immunotherapy through inhibiting the intertumoral accumulation of M2 macrophages (40, 41). For the abundant fibrous stroma observed in C4 subclass, anti-fibrosis drugs (like NOX4 inhibitor) suppressed the activation of cancer-associated fibroblasts and promoted the infiltration of CD8⁺ T cells, ultimately improving the efficacy of immunotherapy (42). As for C1 subclass, the application of cytotoxic and modulating agents which can convert cold tumors to inflamed tumors was a potential strategy (43).

Our workflow is limited by the HCC patients obtained specimens for the first time, as well as the influences of confounding variables such as HBV/HCV infection, alcoholic fatty liver, non-alcoholic fatty liver, and cirrhosis were not considered. We will improve them in the future work.

In conclusion, our study dementated a new landscape for the composition of HCC tumor microenvironment. We identified four immune subclasses with distinct mechanisms of immune escape. The patients from distinct subclasses showed a significant difference in clinical prognosis and response to personalized treatment. Based on transcriptome data, our workflow might help to predict the clinical outcome and to find appropriate treatment strategies for HCC patients.

Data Availability Statement

The data sets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the Clinical Research Ethics Committee of the First Affiliated Hospital College of Medicine, Zhejiang University (2014-334). The patients/participants provided written informed consent to participate in this study. Written informed consent was obtained from the individuals for the publication of any potentially identifiable images or data included in this article.

Author Contributions

Study concept and design: SZ, LZ, and XG. Analysis and interpretation of data: XG. Technical and material support: SY, YW, and CP. Drafting of the manuscript: SZ, LZ, and HH. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by Zhejiang Provincial Public Welfare Technology Research Program (LGF18C100001). This study was also supported by Innovative Research Groups of National Natural Science Foundation of China (no. 81721091), National S&T Major Project of China (no. 2017ZX100203205) and Research Unit Project of Chinese Academy of Medical Sciences (2019-I2M-5-030). This study was also supported by Zhejiang International Science and Technology Cooperation Project (no. 2016C04003).

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.610513/full#supplementary-material

References

1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2018) 68:394–424. doi: 10.3322/caac.21492

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Hou J, Zhang H, Sun B, Karin M. The immunobiology of hepatocellular carcinoma in humans and mice: Basic concepts and therapeutic implications. J Hepatol (2020) 72:167–82. doi: 10.1016/j.jhep.2019.08.014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Prieto J, Melero I, Sangro B. Immunological landscape and immunotherapy of hepatocellular carcinoma. Nat Rev Gastroenterol Hepatol (2015) 12:681–700. doi: 10.1038/nrgastro.2015.173

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Lu C, Rong D, Zhang B, Zheng W, Wang X, Chen Z, et al. Current perspectives on the immunosuppressive tumor microenvironment in hepatocellular carcinoma: challenges and opportunities. Mol Cancer (2019) 18:130. doi: 10.1186/s12943-019-1047-6

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Hegde PS, Chen DS. Top 10 Challenges in Cancer Immunotherapy. Immunity (2020) 52:17–35. doi: 10.1016/j.immuni.2019.12.011

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Hegde PS, Karanikas V, Evers S. The Where, the When, and the How of Immune Monitoring for Cancer Immunotherapies in the Era of Checkpoint Inhibition. Clin Cancer Res (2016) 22:1865–74. doi: 10.1158/1078-0432.CCR-15-1507

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Thorsson V, Gibbs DL, Brown SD, Wolf D, Bortone DS, Ou Yang TH, et al. The Immune Landscape of Cancer. Immunity (2018) 48:812–30.e14. doi: 10.1016/j.immuni

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Sturm G, Finotello F, Petitprez F, Zhang JD, Baumbach J, Fridman WH, et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics (2019) 35:i436–45. doi: 10.1093/bioinformatics/btz363

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Zhang Q, He Y, Luo N, Patel SJ, Han Y, Gao R, et al. Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma. Cell (2019) 179:829–45.e20. doi: 10.1016/j.cell.2019.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol (2016) 17:218. doi: 10.1186/s13059-016-1070-5

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res (2017) 77:e108–10. doi: 10.1158/0008-5472.CAN-17-0307

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Gao Q, Zhu H, Dong L, Shi W, Chen R, Song Z, et al. Integrated Proteogenomic Characterization of HBV-Related Hepatocellular Carcinoma. Cell (2019) 179:561–77.e22. doi: 10.1016/j.cell.2019.08.052

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Roessler S, Jia HL, Budhu A, Forgues M, Ye QH, Lee JS, et al. A unique metastasis gene signature enables prediction of tumor relapse in early-stage hepatocellular carcinoma patients. Cancer Res (2010) 70:10202–12. doi: 10.1158/0008-5472.CAN-10-2607

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al. Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell (2016) 165:35–44. doi: 10.1016/j.cell.2016.02.065

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Pinyol R, Montal R, Bassaganyas L, Sia D, Takayama T, Chau GY, et al. Molecular predictors of prevention of recurrence in HCC with sorafenib as adjuvant treatment and prognostic factors in the phase 3 STORM trial. Gut (2019) 68:1065–75. doi: 10.1136/gutjnl-2018-316408

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Huang H, Ren Z, Gao X, Hu X, Zhou Y, Jiang J, et al. Integrated analysis of microbiome and host transcriptome reveals correlations between gut microbiota and clinical outcomes in HBV-related hepatocellular carcinoma. Genome Med (2020) 12:102. doi: 10.1186/s13073-020-00796-5

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Job S, Rapoud D, Dos Santos A, Gonzalez P, Desterke C, Pascal G, et al. Identification of Four Immune Subtypes Characterized by Distinct Composition and Functions of Tumor Microenvironment in Intrahepatic Cholangiocarcinoma. Hepatology (2020) 72:965–81. doi: 10.1002/hep.31092

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun (2013) 4:2612. doi: 10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Hanzelmann S, Castelo R, Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinf (2013) 14:7. doi: 10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

22. Emran AA, Chatterjee A, Rodger EJ, Tiffen JC, Gallagher SJ, Eccles MR, et al. Targeting DNA Methylation and EZH2 Activity to Overcome Melanoma Resistance to Immunotherapy. Trends Immunol (2019) 40:328–44. doi: 10.1016/j.it.2019.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Rooney MS, Shukla SA, Wu CJ, Getz G, Hacohen N. Molecular and genetic properties of tumors associated with local immune cytolytic activity. Cell (2015) 160:48–61. doi: 10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Llovet JM, Montal R, Sia D, Finn RS. Molecular therapies and precision medicine for hepatocellular carcinoma. Nat Rev Clin Oncol (2018) 15:599–616. doi: 10.1038/s41571-018-0073-4

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Lee JS, Chu IS, Heo J, Calvisi DF, Sun Z, Roskams T, et al. Classification and prediction of survival in hepatocellular carcinoma by gene expression profiling. Hepatology (2004) 40:667–76. doi: 10.1002/hep.20375

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Boyault S, Rickman DS, de Reyniès A, Balabaud C, Rebouissou S, Jeannot E, et al. Transcriptome classification of HCC is related to gene alterations and to new therapeutic targets. Hepatology (2007) 45:42–52. doi: 10.1002/hep.21467

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Chiang DY, Villanueva A, Hoshida Y, Peix J, Newell P, Minguez B, et al. Focal gains of VEGFA and molecular classification of hepatocellular carcinoma. Cancer Res (2008) 68:6779–88. doi: 10.1158/0008-5472.CAN-08-0742

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Hoshida Y, Nijman SM, Kobayashi M, Chan JA, Brunet JP, Chiang DY, et al. Integrative transcriptome analysis reveals common molecular subclasses of human hepatocellular carcinoma. Cancer Res (2009) 69:7385–92. doi: 10.1158/0008-5472.CAN-09-1089

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Lachenmayer A, Alsinet C, Savic R, Cabellos L, Toffanin S, Hoshida Y, et al. Wnt-pathway activation in two molecular classes of hepatocellular carcinoma and experimental modulation by sorafenib. Clin Cancer Res (2012) 18:4997–5007. doi: 10.1158/1078-0432.CCR-11-2322

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Benson AB3, D’Angelica MI, Abbott DE, Abrams TA, Alberts SR, Saenz DA, et al. NCCN Guidelines Insights: Hepatobiliary Cancers, Version 1.2017. J Natl Compr Canc Netw (2017) 15:563–73. doi: 10.6004/jnccn.2017.0059

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Sprinzl MF, Reisinger F, Puschnik A, Ringelhan M, Ackermann K, Hartmann D, et al. Sorafenib perpetuates cellular anticancer effector functions by modulating the crosstalk between macrophages and natural killer cells. Hepatology (2013) 57:2358–68. doi: 10.1002/hep.26328

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Chen ML, Yan BS, Lu WC, Chen MH, Yu SL, Yang PC, et al. Sorafenib relieves cell-intrinsic and cell-extrinsic inhibitions of effector T cells in tumor microenvironment to augment antitumor immunity. Int J Cancer (2014) 134:319–31. doi: 10.1002/ijc.28362

PubMed Abstract | CrossRef Full Text | Google Scholar

33. El-Khoueiry AB, Sangro B, Yau T, Crocenzi TS, Kudo M, Hsu C, et al. Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet (2017) 389:2492–502. doi: 10.1016/S0140-6736(17)31046-2

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Gibney GT, Weiner LM, Atkins MB. Predictive biomarkers for checkpoint inhibitor-based immunotherapy. Lancet Oncol (2016) 17:e542–51. doi: 10.1016/S1470-2045(16)30406-5

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Berraondo P, Ochoa MC, Olivera I, Melero I. Immune Desertic Landscapes in Hepatocellular Carcinoma Shaped by β-Catenin Activation. Cancer Discovery (2019) 9:1003–5. doi: 10.1158/2159-8290.CD-19-0696

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Pardoll DM. The blockade of immune checkpoints in cancer immunotherapy. Nat Rev Cancer (2012) 12:252–64. doi: 10.1038/nrc3239

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Miao D, Margolis CA, Gao W, Voss MH, Li W, Martini DJ, et al. Genomic correlates of response to immune checkpoint therapies in clear cell renal cell carcinoma. Science (2018) 359:801–6. doi: 10.1126/science.aan5951

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Noy R, Pollard JW. Tumor-associated macrophages: from mechanisms to therapy. Immunity (2014) 41:49–61. doi: 10.1016/j.immuni.2014.06.010

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Trédan O, Galmarini CM, Patel K, Tannock IF. Drug resistance and the solid tumor microenvironment. J Natl Cancer Inst (2007) 99:1441–54. doi: 10.1093/jnci/djm135

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Mok S, Koya RC, Tsui C, Xu J, Robert L, Wu L, et al. Inhibition of CSF-1 receptor improves the antitumor efficacy of adoptive cell transfer immunotherapy. Cancer Res (2014) 74:153–61. doi: 10.1158/0008-5472.CAN-13-1816

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Paulus P, Stanley ER, Schäfer R, Abraham D, Aharinejad S. Colony-stimulating factor-1 antibody reverses chemoresistance in human MCF-7 breast cancer xenografts. Cancer Res (2006) 66:4349–56. doi: 10.1158/0008-5472.CAN-05-3523

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ford K, Hanley CJ, Mellone M, Szyndralewiez C, Heitz F, Wiesel P, et al. NOX4 Inhibition Potentiates Immunotherapy by Overcoming Cancer-Associated Fibroblast-Mediated CD8 T-cell Exclusion from Tumors. Cancer Res (2020) 80:1846–60. doi: 10.1158/0008-5472.CAN-19-3158

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Galon J, Bruni D. Approaches to treat immune hot, altered and cold tumours with combination immunotherapies. Nat Rev Drug Discovery (2019) 18:197–218. doi: 10.1038/s41573-018-0007-y

CrossRef Full Text | Google Scholar