Myeloid signature reveals immune contexture and predicts the prognosis of hepatocellular carcinoma

Hepatocellular carcinoma (HCC) is one of the most prevalent malignancies worldwide, with an increasing incidence and a dismal prognosis (1–4). In contrast to most other malignancies, 90% of HCCs arise from chronic inflammation (3–7). A variety of immune cells contribute to the development and progression of HCC (5, 6, 8, 9), constituting a highly complex, dynamic, and interactive immune contexture. This ever-evolving immune cell landscape results in dramatic heterogeneity among immune infiltrates in HCC tumors and great variability between patients (5, 8, 10). Although our knowledge of the tumor microenvironment (TME) has substantially increased over the past several decades, challenges remain in interpreting the immune microenvironment and translating this information into clinical benefits (11). Thus, it would be ideal to construct a simple and reliable signature of local immune responses for individual HCC patients that may depict the immune contexture, reveal prognostic information, and predict therapeutic response.

Myeloid cells are a group of heterogeneous innate immune cells, including monocytes/macrophages (Mo/Mφs), polymorphonuclear leukocytes (PMNs), myeloid-derived suppressor cells (MDSCs), and dendritic cells. These cells are important components of the tumor tissue and are key regulators of the immune milieu. In recent years, myeloid cells have attracted particular attention as a result of a growing body of evidence indicating that they strikingly influence therapeutic responses (12). For example, programed death-ligand 1 (PD-L1) expression on myeloid cells is essential for anti–PD-L1–mediated tumor regression (13, 14), and successful anti-CTLA4 and anti–programmed cell death protein 1 (anti–PD-1) treatments require that the myeloid infiltrate expands and sustains the therapeutic response, even though these cells are not the direct targets of the treatments (15, 16). However, in contrast to the identification of T cell–related biomarkers, it has been difficult to establish a signature that can determine the protumor or antitumor local myeloid response, largely because of the high heterogeneity/plasticity of the myeloid compartment. In many cases, a myeloid cell subtype identified by general markers consists of subsets with distinct or even opposing impacts on tumor progression. Moreover, these pro- and antitumor myeloid subtypes often simultaneously coexist in a tumor with correlated cell numbers, potentially competing for control of the balance of the local myeloid response. Owing in part to this correlation, the clinical implications of particular myeloid cells are not always in accordance with their expected roles in cancer. For example, we and others have found that the presence of immune cells (e.g., Mφs) expressing the suppressive immune checkpoint protein PD-L1 in cancer nests is associated with a favorable, rather than a negative, prognosis in HCC (17) and ovarian serous carcinoma (18). Therefore, a valid myeloid signature that can reveal the systemic effect of the local myeloid response, specifically to attack, ignore, or support the tumor, should be a prerequisite for understanding individual myeloid cell behaviors and may facilitate precision therapy targeting the TME.

In the present study, we screened myeloid markers that cover major myeloid subtypes in human tumors and then constructed and validated a simple prognostic myeloid signature, called the myeloid response score (MRS), based on a data set from a multicenter cohort of 1177 HCC patients. Estimation of the MRS relies on only 2 intratumoral immunohistochemistry (IHC) features, namely, CD11b and CD169. HCCs with different MRSs showed distinct intratumoral myeloid contextures that were associated with different T lymphocyte activation statuses. Moreover, MRS-based nomograms exhibited notable discriminatory ability, accuracy, and clinical usefulness in predicting postsurgery HCC prognosis, and evidence suggests that the MRS may be associated with sorafenib efficacy for recurrent HCC and with the immune class.

Construction of the MRS. To build a technically simple and reliable myeloid signature that reflects the myeloid contexture and has a potential clinical impact, we first screened potential myeloid markers in a primary cohort of 488 HCC patients (see schematic overview of the study design in Figure 1). Nine common myeloid markers associated with HCC progression were assessed, covering the major tumor-infiltrating myeloid subtypes (Supplemental Figure 1 and Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/JCI135048DS1). Among these, CD11b and CD33 are common myeloid markers widely expressed on myeloid subsets and are often used to mark MDSCs. CD15 is a granulocyte-specific marker, and S100 is expressed on dendritic cells. CD68, CD169, CD163, CD204, and CD206 are common Mo/Mφ markers with different prognostic implications. In general, a lesser extent of myeloid cell infiltration, with higher variances across patients, was present in the tumor than in the peritumoral tissue (Figure 2A). Correlation analyses showed that despite their different prognostic implications (Supplemental Figures 2 and 3), the expressions of various myeloid markers were closely associated with each other in the same region, in both intratumoral and peritumoral tissues (Figure 2B).

Figure 1

Study design. A schematic overview of the study (more cohort details are provided in Supplemental Table 2 and Supplemental Methods). FUZH, Fudan University Zhongshan Hospital; IF, immunofluorescence; SYSUCC, Sun Yat-sen University Cancer Center; THZP, Taizhou Hospital of Zhejiang Province.

Figure 2

Construction of the MRS. (A) Expression of 18 myeloid features (9 markers in 2 regions) in peritumor liver and intratumor tissues of HCC (n = 488). ***P < 0.001 (paired t test). Error bars, median and IQR. (B) Correlation between the 18 myeloid immune features. Correlation coefficients shown in black indicate a significant correlation (P < 0.05; Spearman rank correlation); coefficients in gray indicate an insignificant correlation (P ≥ 0.05; Spearman rank correlation). (C) Construction of the MRS using a lasso Cox model. Lasso coefficient profiles of the 18 HCC-associated myeloid features are shown. (D) Tenfold cross-validation for fine-tuning parameter selection in the lasso model. Red dots and solid gray lines represent the partial likelihood deviance ± SE. The dotted vertical lines are drawn at the optimal values by minimum criteria (black) and 1-SE criteria (red).

To construct a myeloid signature with minimal IHC features, we used an L1-penalized (lasso) Cox regression model to fit the time-to-recurrence (TTR) data from the training subset, which was randomly selected from the primary cohort (Supplemental Tables 2 and 3). In the optimized model, 2 of the 18 myeloid features — the expression of 9 markers in 2 regions, peritumoral liver tissue (P) and intratumoral tissue (T) — were selected, namely, CD11b_T and CD169_T (Figure 2, C and D). A risk score was built on the basis of the penalized coefficients (Figure 2, C and D) and was then rescaled to range from 0 to 100 (Supplemental Figure 4). This risk score, referred to hereafter as the MRS, was calculated as follows: MRS = 0.161 × CD11b_T – 0.106 × CD169_T + 35 (0 ≤ MRS ≤ 100). To assess methodological consistency, we found that the automated enumeration method used in this study was highly consistent with the counting results obtained by a pathologist (Supplemental Figure 5), suggesting that this MRS algorithm is readily applicable to data from both counting methods. Moreover, using a separate set of 10 HCCs, we evaluated the MRS obtained from at least 6 different intratumoral tissue blocks from each tumor. According to the results, the variance in MRS values was significantly smaller (coefficients of variation, 0%–11.4%) than the variance in individual markers (CD11b and CD169) in the random sampling of tumor blocks (Supplemental Figure 6).

To examine the prognostic power of the MRS, we performed survival analysis using both the training and test subsets. Time-dependent receiver operating characteristic (ROC) curve analysis showed that the MRS, as a continuous score, had a high discriminatory ability for predicting the TTR and overall survival (OS) of patients after HCC resection (Figure 3, A and B). Additionally, the MRS allowed us to classify the HCC patients into subgroups with distinct clinical outcomes (MRS^lo, 0–37.9; MRS^int, 38.0–60.6; and MRS^hi, 60.7–100; cutoffs were determined with the X-Tile program [ref. 19], as shown in Supplemental Figure 7). In the test subset, the 2-year recurrence rates of patients with low (MRS^lo; 39.3% of patients) and high (MRS^hi; 23.8% of patients) scores were 24.5% and 76.2%, respectively. The 3-year OS rates of patients in the MRS^lo and MRS^hi subgroups were 91.7% and 33.8%, respectively. Kaplan-Meier curves and Cox univariate regression analyses further demonstrated that the MRS correlated strongly with both TTR and OS (Figure 3, A and B). These data suggest that the IHC-based MRS may serve as a novel and reliable marker for predicting the prognosis of HCC patients after curative resection.

Figure 3

Prognostic performance of the MRS in the primary cohort. (A) The discriminatory power of the MRS for TTR and OS in the training subset of the primary cohort. (B) The discriminatory power of the MRS for TTR and OS in the test subset of the primary cohort. P values were calculated using the log-rank test for trend. AUC, area under the receiver operating characteristic curve; RF, recurrence free.

We also used a similar formula to calculate the peritumoral MRS (pMRS). The vast majority of the pMRS values were below the MRS^lo/MRS^int cutoff (>84% in the training subset and >77% in the test subset), suggesting that the peritumoral myeloid cell composition was similar to the MRS^lo environment (Supplemental Figure 8A). Interestingly, the pMRS and the intratumoral MRS derived from the same tissue were not significantly correlated (Supplemental Figure 8, B and C).

Validating the prognostic performance of the MRS. The prognostic power of the MRS was further examined in an independent internal cohort and in 2 external cohorts (Supplemental Figures 4 and 9). Time-dependent ROC curve analysis indicated that the MRS, as a continuous parameter, maintained its discriminatory power in these validation cohorts (Figure 4, A–C, left, and Supplemental Figure 10). After classification of the patients into subgroups using the same criteria derived from the primary training subset, these MRS subgroups presented distinct TTR and OS in the validation cohorts (Figure 4, A–C, right). To determine the relationship between the MRS and different clinicopathological characteristics, we performed χ² tests and found a significant association between the MRS and the presence of vascular invasion, Barcelona Clinic Liver Cancer (BCLC) stage, and tumor node metastasis (TNM) stage (Supplemental Tables 4 and 5). Univariate Cox regression models (Supplemental Tables 6–9), ROC curve analyses (Supplemental Figure 11), and Kaplan-Meier curves (Supplemental Figures 12 and 13) demonstrated the consistent ability of the MRS to distinguish patients with distinct clinical outcomes in subgroups with different levels of risk factors. The MRS also maintained its predictive power in different subgroups stratified by BCLC or TNM stage (Supplemental Figure 14). Moreover, multivariate Cox proportional hazards regression analysis further verified that the MRS could serve as an independent predictor for tumor recurrence and OS of HCC patients in the primary, internal validation, and external validation cohorts (Supplemental Table 10).

Figure 4

Validation of the MRS in independent internal and external cohorts. (A) The prognostic performance of the MRS in the internal validation cohort. (B) The prognostic performance of the MRS in the FUZH cohort. (C) The prognostic performance of the MRS in the THZP cohort. P values were calculated using the log-rank test for trend.

The MRS reflects the differential composition of myeloid cells. As the prognostic value of the MRS indicated that this 2-feature-based signature might reflect the distinct myeloid immune contexture that affects HCC progression, we analyzed the distribution of myeloid infiltrates in the different MRS subgroups. The density of CD169⁺ cells was significantly higher in MRS^lo tumors (Figure 5A) and was associated with a good prognosis (Supplemental Figures 2 and 3). These CD169⁺ cells in HCC constitute a subset of Mφs with relatively high expression of HLA-DR and costimulatory signal CD86 (20, 21). In contrast, the density of CD11b⁺ cells, along with those of CD15⁺, CD163⁺, CD33⁺, CD68⁺, and CD204⁺ infiltrates, was increased in MRS^hi HCC (Figure 5A). Among these molecules, CD11b, CD33, and CD15 have been described as markers for tumor-associated PMNs/PMN-MDSCs (22), whereas CD163, CD204, and CD68 are thought to be markers for tumor-associated Mφs (23). Interestingly, compared with MRS^lo and MRS^hi HCCs, MRS^int HCCs showed less infiltration of most myeloid subsets (Figure 5A).

Figure 5

MRS reflects different myeloid immune contextures between HCC tumors. (A) Densities of myeloid subsets in the intratumoral (T) area of tumors with low (n = 191), intermediate (n = 182), and high (n = 115) MRSs. Error bars, median and IQR. *P < 0.05, **P < 0.01, ***P < 0.001 (2-way ANOVA followed by Tukey’s test). (B) Representative multiplexed immunofluorescence images show different expression patterns of 5 common myeloid markers in subgroups of HCCs. Scale bar: 100 μm. (C) Different compositions of myeloid contextures in MRS^lo (n = 6), MRS^int (n = 5), and MRS^hi (n = 5) tissues. Data are shown in an UpSet plot. Each row corresponds to 1 myeloid marker, and each column (a subgroup including 3 bars) corresponds to a subset of cells with the indicated pattern of marker coexpression. Circles are either light gray, indicating that this subset was not, or was marginally, stained for that marker, or black, showing that this subset expressed that marker. Error bars, mean and SEM. *P < 0.05, **P < 0.01 (2-way ANOVA followed by Tukey’s test).

To further characterize the phenotypes of myeloid cells in HCC, we performed multiplexed immunofluorescence analysis. The peritumoral myeloid contextures were similar between HCCs with different MRSs (Supplemental Figure 15). In the tumor, however, HCCs with low, intermediate, and high MRSs displayed patterns of immunocompetent, immunodeficient, and immunosuppressive TMEs, respectively (Figure 5, B and C). Specifically, MRS^lo HCC exhibited a myeloid cell contexture (Figure 5, B and C) that was comparable to that of peritumoral tissue (Supplemental Figure 15, B and C), suggesting that the myeloid subset composition in these tissues was closest to that of normal tissue. Furthermore, CD169⁺CD163⁺CD14⁺CD11b^lo/– Mφs were predominant in MRS^lo HCC (35.9% ± 6.1%, mean ± SEM), but were almost absent in MRS^hi HCC (0.03% ± 0.03%). Instead, CD11b⁺CD15⁺ PMNs/PMN-MDSCs (24.9% ± 4.0%) and CD169^–CD11b⁺CD163⁺ Mo/Mφs (32.5% ± 13.9%; including CD169^–CD11b⁺CD163⁺CD14^– [15.9% ± 3.9%] and CD169^–CD11b⁺CD163⁺CD14⁺ cells [16.6% ± 2.5%]) represented the majority of myeloid cells in MRS^hi tumors, in contrast to their rarity in peritumoral liver tissue (Supplemental Figure 15, B and C) and in both MRS^lo and MRS^int HCCs (Figure 5, B and C). Taken together, these results suggest that low and high MRSs are associated with 2 opposing compositions of myeloid cell subtypes.

Since myeloid cells in MRS^lo tumors were phenotypically similar to their counterparts in the peritumoral tissue, the altered phenotype of myeloid cells in MRS^hi tumors might suggest the impact of tumor environment. Thus, we examined the influence of tumor cell culture supernatant (TSN) on the phenotype of human PMNs and monocyte-derived Mφs. Indeed, CD11b expression was significantly upregulated on PMNs after exposure to TSN from 2 hepatoma cell lines (Supplemental Figure 16), consistent with the increased CD11b⁺CD15⁺ PMNs/PMN-MDSCs in MRS^hi tumors (Figure 5, B and C). Notably, TSN-treated PMNs markedly elevated the intracellular level of reactive oxygen species, which is essential for PMN-MDSCs to induce antigen-specific T cell tolerance (8, 9). In contrast, exposure to TSN induced profound downregulation of CD169 on Mφs. This was accompanied by decreased expression of type II MHC molecule HLA-DR and costimulatory signal CD86, suggesting an impaired antigen-presenting capability. These in vitro data support our in situ finding that the MRS^hi tumors accumulate CD11b⁺CD15⁺ PMNs/PMN-MDSCs and Mφs with low CD169 expression.

The MRS is linked to the T cell activation status in HCC. Given that CD8⁺ T cells are the ultimate tumor-eliminating lymphoid effectors, we expected an inverse correlation between the MRS and intratumoral CD8⁺ T cell density. To our surprise, the number of infiltrating CD8⁺ T cells in MRS^hi tumors was comparable to that in MRS^lo tumors (Figure 6A), and a negative correlation between the MRS and the CD8⁺ T cell number was observed only for the MRS^lo group (Figure 6B). Moreover, a trend of a positive correlation between the MRS and CD8⁺ T cells was observed for the MRS^hi group. The coexistence of abundant CD8⁺ T cells and tumor-supporting myeloid subsets in this group with a poor prognosis further suggests the presence of a highly immunosuppressive TME in MRS^hi HCCs.

Figure 6

A high MRS indicates a more suppressive immune microenvironment. (A) Densities of intratumoral CD8⁺ cells in subgroups of HCCs with low (n = 191), intermediate (n = 182), and high (n = 115) MRSs. **P < 0.01, ***P < 0.001 (ANOVA followed by Tukey’s test). (B) Correlation between MRS values and CD8⁺ T cell infiltration levels in HCC. Spearman’s rank correlation coefficients and P values are shown. (C–E) Enrichment of genes related to tumor-associated Mφs (TAMs) (C), MDSCs (D), and CD8⁺ T cell exhaustion (E) in HCCs with high MRS values. The gene sets are shown in Supplemental Table 12. The P value of the normalized enrichment score (NES) was calculated using permutation test, adjusted for multiple hypothesis testing. (F) Flow cytometric analysis of PD-1 expression on CD8⁺ T cells freshly isolated from HCC tissues. *P < 0.05, **P < 0.01 (Mann-Whitney test).

To characterize the cellular activities underlying the change in the MRS between tumors, we performed gene expression microarray analysis using 21 HCC samples (MRS^lo, n = 10; MRS^hi, n = 11; Supplemental Figure 17A and Supplemental Table 11). According to Gene Ontology (GO) analysis, 24 of the top 25 enriched GO terms among upregulated genes (FDR < 0.05, logFC > 1) in the MRS^hi group were immune related (Supplemental Figure 17B). The gene-concept network categorized these immune-related terms into 3 interrelated key processes: enhanced recruitment of myeloid leukocytes into tumors, activation of suppressive myeloid cells, and induction of CD8⁺ T cell exhaustion (Supplemental Figure 18A). Consistently, gene set enrichment analysis with gene sets featuring tumor-associated Mφs (TAMs) (24), MDSCs (25), and exhausted CD8⁺ T cells (26) showed that these genes related to immunosuppression were markedly overrepresented in MRS^hi HCC tumors (Figure 6, C–E, and Supplemental Table 12).

Among the leading-edge core-enriched genes related to exhausted CD8⁺ T cells were those encoding important immune checkpoint proteins and exhaustion markers, including PDCD1, the gene encoding PD-1 protein. Along this line, in a subset of 410 patients within the primary cohort, we found that in MRS^lo, MRS^int, and MRS^hi subgroups the PD-L1⁺ rates were 26.35% (39/148), 18.90% (31/164), and 41.84% (41/98), respectively. PD-L1–expressing HCC tissues were enriched in the MRS^hi subgroup (χ² = 5.096, P = 0.024), suggesting a correlation between the MRS and PD-L1 expression in HCC. This correlation was more robust between the MRS and the presence of PD-L1⁺ Mφs (Supplemental Figure 18, B and C). Consistently, flow cytometry confirmed that CD8⁺ T cells freshly isolated from MRS^hi HCCs had a significantly higher proportion of the PD-1^hi exhausted subset than did those from MRS^lo tumors (Figure 6F). The associated elevation of PD-L1 expression on Mφs and PD-1 expression on infiltrating CD8⁺ T cells further suggest that MRS^hi tumors represent a type of highly immunosuppressive HCC microenvironment.

TCGA Liver Hepatocellular Carcinoma data support the immunological and prognostic relevance of the MRS. Recently, several computational methods have been proposed for quantitating tumor-infiltrating immune subsets from the transcriptomic profiling data of bulk tumors; however, bioinformatics quantitation of myeloid subsets is not as accurate as that of lymphoid cells in HCC tissues (Supplemental Figure 19), partially because of a relatively minor contribution of myeloid cell–derived RNA to the bulk tissue and the high plasticity of their transcriptomes. Thus, instead of quantifying specific myeloid subpopulations, we evaluated microenvironmental differences associated with different MRSs in reference to differential gene expression. Based on transcriptional similarities identified by the Monocle 2 algorithm (27), the 371 HCCs from The Cancer Genome Atlas Liver Hepatocellular Carcinoma (TCGA-LIHC) project were ordered in “pseudotime” (27), displaying a “climb trajectory” of MRS estimates across these tumors (Figure 7A). The HCC samples were subjected to unsupervised clustering into 5 subgroups (Figure 7A), revealing stepwise increases and decreases in the expression levels of MRS-associated upregulated and downregulated genes, respectively (Figure 7B). In accordance with our findings above, we observed a decrease in T cell cytotoxic activity along with higher enrichment scores for T cell exhaustion genes associated with pseudotime ordering (Figure 7C) and an inverse correlation between MRS estimates and HCC prognosis (Figure 7D). These transcriptomic profiling data provide additional support that a higher MRS correlates with enhanced CD8⁺ T cell exhaustion and worse patient survival.

Figure 7

An increase in the MRS correlates with CD8⁺ T cell dysfunction and poor survival among HCC patients in the TCGA-LIHC data set. (A) The trajectory of MRS increases along pseudotime in a 2-dimensional state-space defined by Monocle 2. The 371 HCCs from the TCGA-LIHC project are placed in the “trajectory” inferred from expression of MRS-associated differentially expressed genes (FDR < 0.01 and |logFC| > 1; Supplemental Table 11). State 1: n = 90; State 2: n = 49; State 3: n = 72; State 4: n = 52; State 5: n = 102. (B) Enrichment scores for genes positively (left) and negatively (right) correlated with MRS were determined by gene set variation analysis (GSVA). (C) Cytotoxic activity and exhaustion levels of T cells inferred from bulk RNA-Seq data. The cytotoxic activity of T cells was calculated as the geometrical mean of PRF1 and GZMA divided by the geometrical mean of CD8A and CD8B. Exhaustion levels are defined by GSVA enrichment scores for CD8⁺ T cell exhaustion–specific genes in HCC. (D) OS of HCC patients in subgroups of different myeloid response states.

Building predictive MRS-based nomograms. The above data show that the MRS is a simple myeloid signature that can reveal significant immune discrepancies between tumors with validated prognostic implication. Therefore, this signature may contribute to an improved prognostic model with a combination of relevant biomarkers. To exemplify the use of the MRS in assessing the postsurgery risk for individual patients, we built 2 Cox proportional hazards regression models to predict postsurgery TTR and OS (see Supplemental Methods for details). The fitting models were simplified and presented with 2 nomograms (Figure 8, A and B). The MRS, as a continuous variable, and several clinicopathological covariates (i.e., tumor size, vascular invasion, cirrhosis state, tumor number, and Child-Pugh score) were included in the final nomograms. These nomograms showed notable discriminatory abilities in predicting recurrence and OS among postsurgical HCC patients (Figure 8C). In addition, based on ROC curve analysis, the AUCs of the nomogram predicting 2-year recurrence and 3-year OS were 0.79 (95% CI, 0.73–0.84) and 0.87 (95% CI, 0.83–0.90) in the primary cohort, 0.75 (95% CI, 0.68–0.80) and 0.84 (95% CI, 0.79–0.88) in the internal validation cohort, and 0.75 (95% CI, 0.69–0.80) and 0.79 (95% CI, 0.74–0.84) in the external validation cohort, respectively (Figure 8D). Although time-dependent ROC curve analysis showed that the MRS alone had a higher prognostic accuracy than both the BCLC and TNM staging systems for the entire cohort of HCC patients, the prognostic performance of the MRS-based nomograms was significantly higher than that of the MRS itself (Supplemental Figure 20). The calibration analysis confirmed the reliability of the nomograms for all cohorts (Figure 8E). Additionally, decision curve analysis was used to assess whether nomogram-assisted decisions on applying a medical intervention may improve patient outcomes (28, 29). The results showed the net benefit of nomogram-assisted decisions at a wide range of threshold probabilities, suggesting potential clinical usefulness of the nomograms in all 3 independent cohorts (Figure 8F). In the entire cohort, these 2 nomograms exhibited superior discriminatory ability and clinical usefulness compared with the BCLC and TNM staging systems for a wide range of medical decisions after HCC resection (Supplemental Figure 21).

Figure 8

Nomograms to predict recurrence and survival probabilities. (A and B) Nomograms to predict recurrence (A) and OS (B) after curative resection for HCC. (C) The indicated cohorts were divided into 3 subgroups according to the total points given by the TTR nomogram shown in A or by the OS nomogram shown in B. Kaplan-Meier analysis was performed for TTR or OS in subgroups of HCC patients. P values were calculated using the log-rank test for trend. (D–F) The time-dependent ROC curve (D), calibration curve (E), and decision curve analysis of the clinical usefulness (F) for the TTR (top) and OS (bottom) nomograms in the 3 independent cohorts. The 2 external cohorts were combined to meet the data size requirement of the analysis. The AUC is indicated as the mean and 95% CI.

MRS and sorafenib treatment for recurrent HCC. In addition to the impact of the myeloid contexture on cancer progression, emerging evidence suggests that myeloid cells may affect the efficacy of antitumor therapies (13–16, 22). Therefore, we sought to assess whether MRS might be associated with the therapeutic response in HCC. Sorafenib is the recommended first-line systemic therapy for advanced HCC worldwide (1, 30, 31), and the efficacy of this treatment is closely associated with activities of immune components in the microenvironment (32–34). In a cohort of 51 patients who had recurrent HCC after curative resection, a daily dose of 800 mg sorafenib was prescribed after the diagnosis of tumor recurrence (see Supplemental Methods and Supplemental Table 13 for cohort details). Recurrence was diagnosed 0.8 to 121.4 months after resection (IQR, 4.9 to 30.1 months), and the duration of sorafenib treatment was 1.7 to 60.8 months (IQR, 6.2 to 17.9 months). After stratifying patients into subgroups according to MRS values using the same criteria derived from the primary training subset, we found that objective responses occurred almost exclusively in the MRS^lo subgroup (Figure 9A). Time-dependent ROC curve analysis showed that the AUCs of the MRS predicting 1-year progression-free survival (PFS) and 2-year PFS were 0.76 (95% CI, 0.62–0.89) and 0.74 (95% CI, 0.58–0.88; Figure 9B), respectively. The median PFS in the MRS^lo, MRS^int, and MRS^hi subgroups was 7.0, 4.4, and 2.8 months after treatment initiation, respectively, based on Kaplan-Meier analysis (Figure 9C). Moreover, log-rank test results showed that PFS was closely associated with the MRS (P = 0.010); specifically, MRS^hi patients had a hazard ratio of 2.45 (95% CI, 1.06–5.65; P = 0.006) relative to MRS^lo patients. These results suggest that the MRS may also serve as a signature to predict the efficacy of postrecurrence sorafenib treatment.

Figure 9

The association of MRS with response to sorafenib and immune class. (A) The MRS and the response to sorafenib treatment for recurrent HCC. Individuals’ OS after sorafenib treatment initiation (diagnosed recurrence), time to response, and duration of response to sorafenib (n = 51). Patients with confirmed progressive disease received further treatment, including transcatheter arterial chemoembolization, radiofrequency ablation, or percutaneous ethanol injection. (B) ROC curve analysis of the MRS, as a continuous score, predicting the progression-free survival (PFS) of patients receiving postrecurrence sorafenib treatment. (C) Kaplan-Meier curves for PFS in different MRS subgroups are shown. The P value was calculated using the log-rank test for trend. (D) The 21 HCC samples with gene expression microarray data were classified into the immune class or the rest subgroups (see more details in Supplemental Table 14). Percentages of samples in each subgroup are shown. (E) The 371 HCC samples from the TCGA-LIHC data set were classified into the immune class or the rest subgroups (see more details in Supplemental Table 15). Percentages of samples in each subgroup are shown. (F) The virtual MRS pseudotime values in the immune class (n = 141) and the rest (n = 230) subgroups are shown.

MRS and immune class. Recently, a molecular classification of HCC based on inflammatory response markers was proposed, according to which about 25% of HCCs could be classified into an “immune class” that might be susceptible to immunotherapy (35). We thus performed the immune molecular classification using the gene expression profiles of 21 HCC samples in our cohort to probe the relationship between the MRS and the immune class. The result showed that 8 of the 11 MRS^hi samples (72.7%) and 2 of 10 MRS^lo samples (20%) were categorized into the immune class (P = 0.016; Figure 9D and Supplemental Table 14). Moreover, in the TCGA-LIHC data set, we identified that HCC samples of the immune class were highly enriched in the subgroup with high MRS estimates (Figure 9E and Supplemental Table 15). Conversely, HCC tissues of the immune class had markedly higher values of MRS pseudotime compared with other HCC samples (Figure 9F). These results suggest the close correlation between the simple myeloid IHC signature MRS and the gene expression–based, immune-specific class of HCCs.

Myeloid cells are crucial regulators of tumor progression and affect virtually all types of cancer therapy (12). However, distinguishing antitumor versus protumor myeloid responses in patients remains a major challenge because of the high heterogeneity and functional plasticity of myeloid cells (12, 36, 37). In this study, we constructed a simple and reliable IHC-based myeloid signature, the MRS, that was able to describe the immune microenvironment of HCC and thereby predict disease prognosis. HCCs with different MRSs differ not only in the intratumoral myeloid contexture but also in adaptive immune reaction status patterns. Moreover, we provide examples of the potential application of MRS, such as building predictive nomograms and linking MRS with the efficacy of postrecurrence sorafenib treatment. The findings indicate the valuable clinical potential of myeloid-related biomarkers in oncology.

Myeloid cells are the first-line defenders of the innate immune system and can orchestrate divergent immune responses according to different pathological conditions. By directly interacting with tumor cells and indirectly modulating the tumor stroma, myeloid cells play various, and sometimes seemingly opposing, roles in cancer progression and therapeutic responses (12, 38, 39). In general, assessing the exact myeloid contexture of a tumor can provide valuable information for predicting cancer prognosis and identifying more precise therapeutic strategies. However, the use of single-cell and “omics” technologies to accurately decipher the myeloid contexture remains challenging, in part because of the substantial quantitative biases and phenotypic alterations that occur during myeloid cell isolation. Thus, in situ immunophenotyping is still the optimal approach for a valid assessment of myeloid infiltrates in a tumor and is a readily adaptable technique for clinical practice. In this study, we constructed the MRS, a myeloid signature that is able to reveal the balance of local myeloid responses in HCC. Irrespective of their distinct prognostic implications, individual myeloid features correlated closely with each other in the same region, in both intratumoral and peritumoral tissues. In contrast, the MRS not only indicated the relative abundance of antitumor and protumor myeloid subsets but also exhibited clear immunological relevance and robust prognostic power. Moreover, the MRS calculated from different intratumoral tissue blocks of individual tumors showed a significantly smaller variance than individual myeloid features, supporting that the MRS will have high reproducibility and reliability in clinical practice.

Recent advancements in HCC therapy (40–42) have stressed the need for valid prognostic and predictive algorithms that define subclasses of patients and ultimately lead to better curative strategies (43). Most prognostic signatures developed to date were generated using transcriptomic (44, 45) or epigenetic data (46, 47) derived from the tumor. In contrast to the development of the Immunoscore (48, 49) in colon cancer, few immune-specific biomarkers have been identified in HCC (35, 50), and most of them are primarily based on features of T lymphocytes. To build a simple and clinically reliable myeloid signature, we selected 9 myeloid markers covering known significant myeloid cell types in HCC. We constructed the MRS using a penalized Cox regression model that incorporated 18 myeloid features (the expression of the 9 markers in both peritumoral and intratumoral regions). This myeloid signature showed remarkable prognostic value, as validated in a large, multicenter cohort of HCC patients. In the present data set, MRS alone exhibited a higher discriminatory ability and greater clinical usefulness than did current staging systems in predicting the TTR and OS of HCC patients. It should be noted that the MRS should be kept in the form of a continuous score when it is clinically used, since any stratification could cause loss of relevant information for predicting prognosis. Moreover, its predictive power for the prognosis of individual patients could be further improved with MRS-based nomograms. In contrast to many nomograms in oncology that require complicated variables, these MRS-based nomograms require only the 2-IHC-feature-based MRS and common clinicopathological information, but exhibited notable discriminatory ability, accuracy, and clinical usefulness. Thus, the MRS not only adds to the growing collection of prognostic algorithms for HCC but also highlights the largely unexplored potential of assessing the balance of the intratumoral myeloid response to predict disease and therapeutic outcomes.

Our data suggest that HCCs with low, intermediate, and high MRS estimates correspond to tumors with immunocompetent, immunodeficient, and immunosuppressive TME, respectively. Despite their distinct clinical prognosis, MRS^hi and MRS^lo tumors have comparable levels of CD8⁺ T cell infiltration, though CD8⁺ T cells in MRS^hi tumors displayed notable exhaustion patterns. These findings suggest that the immune responses in “hot tumors” can be very different and that myeloid cells may act as crucial regulators. Indeed, the increase in the MRS between tumors correlated with a shift in the predominant myeloid populations from CD11b^–CD169⁺ Mφs in MRS^lo tumors to CD11b⁺CD169^– Mo/Mφs and CD15⁺CD11b⁺ PMNs/PMN-MDSCs in MRS^hi tumors. These findings demonstrate the plasticity and heterogeneity of Mφs in the TME. It is now clear that Mφ activation and polarization consist of a range of nonterminal phenotypes rather than 2 binary states of classical (M1) or alternative (M2) activation (51–53), particularly in human tumors. In some animal models and in human HCC, CD169⁺ cells constituted a subset of potent antitumor Mφs (20, 21, 54, 55). Interestingly, CD169 was recently identified as a specific marker for suppressive tumor-associated Mφs in breast and endometrial cancers (56), further highlighting the phenotypic and functional plasticity of Mφs in different organs/cancers. On the other hand, CD11b can negatively regulate proinflammatory pathways (57) and is often expressed by immunosuppressive myeloid cells in the blood and tumors of cancer patients (58). Furthermore, targeting CD11b has recently been shown to attenuate the immunosuppressive myeloid response and thus to sensitize tumors to checkpoint blockade (59). Hence, the MRS, derived from the CD169 and CD11b expression in HCC, may provide a measure for classifying the “hot tumors” with very different myeloid contextures, paving the way for the development of more precise therapeutic strategies for each of these subtypes.

In support of the notion that myeloid cells are broadly involved in treatment responses (12, 39, 60), we found that the MRS also correlated with the efficacy of sorafenib treatment for recurrent HCC. Notably, objective responses were mostly restricted to the MRS^lo subgroup, whereas the MRS^hi subgroup had a significantly higher risk of disease progression. In addition, we found that the MRS was closely associated with the immune class, a molecular classification of HCCs that might be susceptible to immunotherapy (35). These findings support the potential of the MRS as a predictive biomarker for response to therapy, which is consistent with the idea that the myeloid response balance can impact the efficacy of treatments (51, 53, 61–64). Thus, future models incorporating the MRS and related parameters as candidate variables may help to create a more precise criterion to tailor treatment options.

Our study has limitations that include the retrospective design and the insufficient sample size of HBV^– HCC patients in our cohort (Supplemental Figures 12 and 13). Additionally, the biology underlying MRS discrepancies remains unclear. In future studies, the combination of a variety of “omics” profiling methods and mechanistic investigations may help to reveal the mechanism regulating the myeloid contexture in HCC and may provide an opportunity for overcoming resistance to targeted and immune-based therapies by rerouting the local myeloid response. Moreover, although high-dimensional profiling is still challenging when applied to a large number of samples, the MRS may provide a useful method for classifying the intratumoral immune reaction, thus enabling patient stratification.

In conclusion, we defined a prognostic signature that requires only 2 intratumoral myeloid IHC features. The MRS can serve as a measure of the myeloid response in HCC, reflecting the significant shift in the immune contexture. The prognostic power of the MRS and MRS-based nomograms for predicting the TTR and OS of HCC after resection were independently validated in both internal and external cohorts. The clinical potential of the MRS was also demonstrated by linking of this signature to postrecurrence sorafenib treatment and to the immune class of HCC. Together, the data highlight the largely unexplored potential of myeloid signatures in oncology and suggest that MRS can provide the basis for the stratification of tumor immune subtypes that may aid in a wide range of clinical decisions.

Patient populations. Four independent cohorts of HCC patients who underwent primary curative resection at 3 centers between 2001 and 2010 were enrolled. The resection procedure and postoperative surveillance for recurrence were performed as described in our previous articles (60, 65, 66). Details regarding the patient information and treatments are described in Supplemental Methods and Supplemental Table 1.

IHC and evaluation of myeloid features. Paraffin-embedded samples were cut into 4 μm sections, which were processed for IHC as previously described (65, 67), with minor modifications. Details regarding the antibodies used in this study are provided in Supplemental Tables 16 and 17. Using tissue microarray (TMA) of HCC samples from the primary cohort, which contained 488 patients, the densities of myeloid cell subsets in both peritumoral and intratumoral regions of HCC were determined (Figure 2, A and B, and Supplemental Figure 1). IHC-stained tissues were scanned at ×20 magnification and evaluated by a computerized image analysis system built with Vectra 2.0 Automated Quantitative Pathology Imaging System and InForm Cell Analysis 2.2 (PerkinElmer). Images were processed using Nuance Image Analysis Software (PerkinElmer) for unmixing of diaminobenzidine and hematoxylin signals and to establish a spectral library. To segment different tissue regions, a pathologist manually drew the training regions of intratumor, peritumor, or blank areas to train the machine-learning tissue segmentation, which was later applied to all tissue images. Following tissue segmentation, images were subject to color deconvolution using the established spectral library, and cell numbers were estimated with the counting object module of InForm Cell Analysis 2.2 (PerkinElmer). The results of the automated measurements were verified and compared with the counting results obtained by a pathologist (Supplemental Figure 5). Quantification of CD11b_T and CD169_T across cohorts is shown in Supplemental Figure 9. A separate set of 10 HCC patients who underwent curative resection at the Sun Yat-sen University Cancer Center was included to evaluate the reproducibility of the TMA core-based IHC features (Supplemental Figure 6).

Penalized Cox regression model construction. We used an L1-penalized (lasso) Cox regression model to select the most prognostic features from the 18 myeloid features using the 1–standard error (1-SE) criterion and then constructed a 2-myeloid-feature-based scoring algorithm using the penalized coefficients to predict the time to recurrence (TTR) in the training set (Figure 2, C and D). Stratification of patients into subgroups based on this risk score was achieved using the optimal cutoffs generated by X-Tile plots (Supplemental Figure 7).

Data and R code availability. The gene expression microarray data of 21 samples are available from the NCBI’s Gene Expression Omnibus database (GEO GSE134921). R code files used in this study are provided in the supplemental materials.

Statistics. Statistical tests were performed with R software (version 3.5.0), GraphPad Prism 8 (GraphPad Software), or SPSS Statistics 25 (IBM). TTR and OS were defined as the time from resection to the first documented instance of disease recurrence by an independent radiological assessment and to death by any cause, respectively. PFS was defined as the time from the initiation of treatment to the first documented instance of disease progression by an independent radiological assessment or death by any cause, whichever occurred first. Univariate and multivariate analyses were performed using the Cox proportional hazards model. Two-tailed Student’s t test was used to compare the means of 2 groups, and ANOVA was applied for multiple comparisons. Associations between parametric, nonparametric, and stratified variables were evaluated using Pearson, Spearman, and χ² tests (or Fisher’s exact test, when appropriate), respectively. P values less than 0.05 were considered significant.

Study approval. For experiments using human samples, all samples were anonymously coded in accordance with the local ethical guidelines (as stipulated by the Declaration of Helsinki). This study was approved by the Institutional Review Boards of the Sun Yat-sen University Cancer Center (Guangzhou, China), the Taizhou Hospital of Zhejiang Province (Taizhou, China), and Fudan University Zhongshan Hospital (Shanghai, China). Written informed consent was obtained from all patients.